FN2_Segment_and_Intensity.py

"""
******************************************************************************************
		Written by Jim Rowe in collaboration with Alexander Jones
								Started: 2022-08-01		
							 		@BotanicalJim
							james.rowe at sheffield.ac.uk
									Version v2.01



******************************************************************************************
"""


# *******************************import libraries*****************************************
from ij 			import IJ, ImageStack, ImagePlus,CompositeImage
from fiji.util.gui 	import GenericDialogPlus
from ij.process 	import ImageProcessor, StackStatistics, ImageConverter, FloatProcessor, ImageConverter
from ij.measure 	import ResultsTable
from array 			import array, zeros
from java.lang 		import Thread
from ij.plugin 		import Slicer
import json
from os.path 		import exists

import math
from java.awt.event	import ActionListener 
from java.awt 		import GridLayout

# *******************************Classes************************************************

class buttonClick(ActionListener):
	"""Class which unique function is to handle the button clics"""
	def __init__(self):
		self.bypass = 0 
		self.quickload = 0
		self.save = 0 
		self.load = 0
	def actionPerformed(self, event): # self to state that the method will be associated to the class instances
		self.bypass=0
		# Check from where comes the event
		source = event.getSource()
		

		# Do an action depending on the button clicked
		if source.label == "Quickload nuclei settings":
			self.thresholdMethod = 'Otsu'
			self.gaussianSigma =1.2 
			self.largeDoGSigma = 6.0
			self.DoG =True
			self.manualSegment= False
			self.dilation=0
			self.sizeExclude= False 
			self.watershed = True
			self.pixelByPixel = False
			self.bypass = 1
			self.quickload=1
			# Do an action depending on the button clicked
		if source.label == "Quickload diffuse settings":
			self.thresholdMethod = 'Otsu'
			self.gaussianSigma =2
			self.largeDoGSigma = 6.0
			self.DoG =False
			self.manualSegment= False
			self.dilation=0
			self.sizeExclude= False 
			self.watershed = False
			self.pixelByPixel = True
			self.bypass = 1 
			self.quickload=1
			self.load=0
		if source.label == "Load settings file":
			loadFilepath=self.fileSelectDialog()
			if exists(loadFilepath):
				print('exists')
				sf=file(loadFilepath, 'rb')
				self.options=json.load(sf)

				print 'options loaded'
			else: 
				print 'failed to load'
			self.bypass = 1			
			self.load = 1
			self.quickload = 0
		if source.label == "Save settings file":
			self.saveFilepath=self.fileSelectDialog()
			self.save = 1
			self.bypass = 1
			self.quickload = 0
			print self.saveFilepath
			
	def fileSelectDialog(self):
		"""file location"""
		gd = GenericDialogPlus("Select file location")
		gd.addFileField("Select settings file .json file", "")
		gd.showDialog()
		
		if gd.wasCanceled():
			IJ.exit()
		settingFilePath =gd.getNextString()
	
		return 	settingFilePath
# *******************************functions************************************************
def globalBackSub(labelGFX, quantGFX, otherGFX):
	"""Requires a labelGFX image, the quantGFX to be quantified (also the output) and, one sacrificial otherGFX images"""
	results=ResultsTable()
	clij2.statisticsOfBackgroundAndLabelledPixels(quantGFX, labelGFX, results)
	IJ.log(str(results.getValue("MEAN_INTENSITY",0)))
	clij2.addImageAndScalar(quantGFX, otherGFX, -results.getValue("MEAN_INTENSITY",0))
	clij2.copy(otherGFX,quantGFX)
	return quantGFX
	
	
def createSubtractionLabels(labelGFX, gfx2, outputGFX):
	"""Requires a labelGFX image, two sacrificial otherGFX images, one of which will be the output"""
	clij2.dilateLabels(labelGFX, gfx2 , 4)
	clij2.subtractImages(gfx2, labelGFX, outputGFX)
	return outputGFX
	
	
def localLabelBackSub(dilatedLabelGFX, labelGFX, quantGFX, otherGFX1, otherGFX2):
	"""Requires a dilatedlabelGFX image, labelGFX, a quantGFX image, and two sacrificial otherGFX images"""
	results4=ResultsTable()
	clij2.statisticsOfBackgroundAndLabelledPixels(quantGFX, dilatedLabelGFX, results4)
	intensities=results4.getColumn(12)
	fp= FloatProcessor(len(intensities), 1, intensities, None)
	intensitiesImp= ImagePlus("IntensitiesImp", fp)
	intGFX=clij2.push(intensitiesImp)
	clij2.replaceIntensities(labelGFX, intGFX, otherGFX1)
	intGFX.close()
	clij2.subtractImages(quantGFX, otherGFX1, otherGFX2)
	clij2.copy(otherGFX2, quantGFX)
	return quantGFX
			
							
def extractChannel(imp, nChannel, nFrame):
	"""extract a channel from the image, at a given frame returning a new imagePlus labelled with the channel name"""
	stack = imp.getImageStack()
	ch=ImageStack(imp.width, imp.height)
	for i in range(imp.getNSlices()):
		index = imp.getStackIndex(nChannel, i, nFrame)
		ch.addSlice(str(i), stack.getProcessor(index))
	imp3 = ImagePlus("Channel " + str(nChannel), ch).duplicate()
	stats =StackStatistics(imp3) 
	IJ.setMinAndMax(imp3, stats.min, stats.max)
	return imp3

def extractFrame(imp, nFrame):
	"""extract a frame from the image, returning a new 16 bit imagePlus labelled with the channel name"""
	stack = imp.getImageStack()
	fr=ImageStack(imp.width, imp.height)
	for i in range(1, imp.getNSlices() + 1):
		for nChannel in range(1, imp.getNChannels()+1):
			index = imp.getStackIndex(nChannel, i, nFrame)
			fr.addSlice(str(i), stack.getProcessor(index))
	imp3 = ImagePlus("Frame " + str(nFrame), fr).duplicate()
	imp3.setDimensions(imp.getNChannels(), imp.getNSlices(), 1)
	comp = CompositeImage(imp3, CompositeImage.COMPOSITE)  
	comp.show()
	return comp

def errorDialog(message):
	"""Outputs a given error for end users"""
	gd = GenericDialogPlus("Error")
	gd.addMessage(message)
	gd.showDialog()
	return
	
def concatStacks(masterStack, impToAdd):
	"""takes an IMP and adds it to a stack, returning the concatenated stack"""
	impToAddStack=impToAdd.getImageStack()
	for i in xrange(1, impToAdd.getNSlices()+1):
		try:	
			masterStack.addSlice(impToAddStack.getProcessor(i))	
		except: print "FAILED To addto stack for: "+ impToAdd.getTitle() +" " + str(i)	
	return masterStack

def previewDialog(imp, options):
	"""Generates the settings dialog and preview window, which live updates dependent on chosen settings"""
	gd = GenericDialogPlus("FRETENATOR2: Segment and Quantify")
	#unpack default settings
	segmentChannel, donorChannel, acceptorChannel, acceptorChannel2, thresholdMethod, maxIntensity, gaussianSigma, largeDoGSigma, DoG,  manualSegment, manualThreshold, makeNearProj, dilation, sizeExclude, minSize, maxSize, watershed, backsubVal, pixelByPixel, saveSettings =options
	#create a list of the channels in the provided imagePlus
	types = []
	for i in xrange(1, imp.getNChannels()+1):
		types.append(str(i))
	gd.setInsets(0,0,0)
	gd.addMessage("""Channel choices:""")
	#user can pick which channel to base the segmentation on
	gd.addChoice("Segmentation", types, str(segmentChannel))

	gd.setInsets(0,0,0)
	gd.addMessage("""Segmentation settings:""")
	methods=["Otsu","Default", "Huang", "Intermodes", "IsoData", "IJ_IsoData", "Li", "MaxEntropy", "Mean", "MinError", "Minimum", "Moments", "Percentile", "RenyiEntropy", "Shanbhag", "Triangle", "Yen"]
	gd.setInsets(0,15,0)
	gd.addCheckbox("Difference of Gaussian instead of Gaussian?", DoG)
	gd.addSlider("Gaussian / DoG 1 sigma", 0.5, 10, gaussianSigma, 0.1)
	gd.addSlider("DoG 2 sigma", 0.5, 20, largeDoGSigma,0.1)
	gd.setModal(False)
	gd.addChoice("Autosegmentation method", methods, thresholdMethod)
	gd.addCheckbox("Manually set threshold? ", manualSegment)
	gd.addSlider("Manual threshold", 10, 65534, manualThreshold, 1)
	
	dilationOptions=["0", "1", "2","3", "4", "5", "6"]
	
	gd.addChoice("Dilation?", dilationOptions, str(dilation))
	gd.addCheckbox("Size exclusion of ROI? ", sizeExclude)
	gd.addSlider("Minimum ROI size", 0, 9999, minSize, 1)
	gd.addSlider("Maximum ROI size", 1, 10000, maxSize, 1)
	gd.addCheckbox("Watershed object splitting? ", watershed)
	
	gd.setInsets(0,0,0)
	gd.addMessage("""Analysis settings:""")
	backsubOpts=["Off", "Local label based", "Global mean"]
	intensities=["254", "4094", "65534"]
	gd.addChoice("Saturation removal threshold", intensities, str(maxIntensity))
	gd.addChoice("Background subtraction", backsubOpts, backsubOpts[backsubVal])
	gd.addCheckbox("""Use pixel by pixel analysis?""", pixelByPixel)
	gd.addCheckbox("Create nearest point projection? ", makeNearProj)
	gd.addCheckbox("Save settings as default? ", False)
	buttonListener = buttonClick()
	gd.addButton("Quickload nuclei settings", buttonListener)
	gd.addButton("Quickload diffuse settings", buttonListener)
	
	gd.addButton("Load settings file", buttonListener)
	gd.addButton("Save settings file", buttonListener)
	gd.setInsets(0,0,0)
	gd.addMessage("""For tutorials, click 'help'.
	
	Rowe (2022) Met Mol Bio | doi.org/10.1007/978-1-0716-2297-1_17
	Rowe (2023) Nat Plants | doi.org/10.1038/s41477-023-01447-4 """)
	gd.addHelp("https://github.com/JimageJ/FRETENATOR2/blob/main/README.md")
	gd.setLocation(0,0)
	#gd.setLayout(GridLayout(0,8))

	gd.showDialog()

		
	cal = imp.getCalibration()
	pixelAspect=(cal.pixelDepth/cal.pixelWidth)
	
	originalTitle=imp1.getTitle()
	
	choices=gd.getChoices()
	sliders=gd.getSliders()
	checkboxes=gd.getCheckboxes()		
	segmentChannel=int(choices.get(0).getSelectedItem())
	thresholdMethod=choices.get(1).getSelectedItem()
	dilation=int(choices.get(2).getSelectedItem())
	maxIntensity=int(choices.get(3).getSelectedItem())
	backSub= choices.get(4).getSelectedItem()
	
	gaussianSigma=sliders.get(0).getValue()/10.0
	largeDoGSigma = gd.sliders.get(1).getValue()/10.0
	manualThreshold=gd.sliders.get(2).getValue()
	minSize = gd.sliders.get(3).getValue()
	maxSize = gd.sliders.get(4).getValue()
	
	
	

	DoG = gd.checkboxes.get(0).getState()
	manualSegment = gd.checkboxes.get(1).getState()
	sizeExclude=gd.checkboxes.get(2).getState()
	watershed = gd.checkboxes.get(3).getState()
	#backSub = gd.checkboxes.get(4).getState()
	pixelByPixel = gd.checkboxes.get(4).getState()
	#print dir(gd.sliders.get(5))
	#print maxSize
	
	segmentChannelOld=segmentChannel
	thresholdMethodOld=thresholdMethod
	maxIntensityOld=maxIntensity
	gaussianSigmaOld=gaussianSigma
	largeDoGSigmaOld= largeDoGSigma
	DoGOld=DoG
	manualSegmentOld= manualSegment
	manualThresholdOld=manualThreshold
	dilationOld=dilation
	sizeExcludeOld=sizeExclude
	minSizeOld=minSize
	maxSizeOld=maxSize
	watershedOld=watershed
	clij2.clear()
	
	segmentImp=extractChannel(imp1, segmentChannel, 0)

	try:
		gfx1=clij2.push(segmentImp)
		gfx2=clij2.create(gfx1)
		gfx3=clij2.create(gfx1)
		gfx4=clij2.create(gfx1)
		gfx5=clij2.create(gfx1)
		gfx7=clij2.create([imp.getWidth(), imp.getHeight()])
	except:	
		try:
		
			Thread.sleep(500)
			print("Succeeded to sending to graphics card on the second time...")
			gfx1=clij2.push(segmentImp)
			gfx2=clij2.create(gfx1)
			gfx3=clij2.create(gfx1)
			gfx4=clij2.create(gfx1)
			gfx5=clij2.create(gfx1)
			gfx7=clij2.create([imp.getWidth(), imp.getHeight()])
		except:
			errorDialog("""Could not send image to graphics card, it may be too large!
		
			Easy solutions: Try	processing as 8-bit, cropping or scaling the image, or
			select a different CLIJ2 GPU.

			This issue is often intermittent, so trying again may also work! 

			See the "Big Images on x graphics cards' notes at:
			https://clij2.github.io/clij2-docs/troubleshooting for more solutions
			
			"""	+ str(clij2.reportMemory()) )


	gfx1,gfx2,gfx3,gfx4,gfx5 = segment(gfx1,gfx2,gfx3,gfx4,gfx5, gaussianSigma, thresholdMethod, maxIntensity, largeDoGSigma, pixelAspect, originalTitle, DoG,  manualSegment, manualThreshold, dilation, sizeExclude, minSize, maxSize, watershed)
	clij2.maximumZProjection(gfx5, gfx7)

	labelPrevImp= clij2.pull(gfx7)
	IJ.setMinAndMax(labelPrevImp, 0,clij2.getMaximumOfAllPixels(gfx7))
	labelPrevImp.setTitle("Preview segmentation")
	labelPrevImp.show()
	
	IJ.run("glasbey_inverted")
	buttonListener.bypass = 0
	while ((not gd.wasCanceled()) and not (gd.wasOKed())):
		

	#IF SETTINGS ARE LOADED, REMAKE THE DIALOG TO APPLY THEM
		if buttonListener.bypass ==1:

			if buttonListener.quickload ==1:
				thresholdMethod= buttonListener.thresholdMethod
				gaussianSigma = buttonListener.gaussianSigma
				largeDoGSigma = buttonListener.largeDoGSigma
				DoG = buttonListener.DoG
				manualSegment = buttonListener.manualSegment
				dilation = buttonListener.dilation
				sizeExclude= buttonListener.sizeExclude
				watershed = buttonListener.watershed
				pixelByPixel = buttonListener.pixelByPixel
				buttonListener.quickload ==0
			if buttonListener.load ==1:
				segmentChannel, donorChannel, acceptorChannel, acceptorChannel2, thresholdMethod, maxIntensity, gaussianSigma, largeDoGSigma, DoG,  manualSegment, manualThreshold, makeNearProj, dilation, sizeExclude, minSize, maxSize, watershed, backSub, pixelByPixel, saveSettings =buttonListener.options
			gd.dispose()
			if buttonListener.save==1:
				options= segmentChannel, donorChannel, acceptorChannel, acceptorChannel2, thresholdMethod, maxIntensity, gaussianSigma, largeDoGSigma, DoG,  manualSegment, manualThreshold, makeNearProj, dilation, sizeExclude, minSize, maxSize, watershed, backsubVal, pixelByPixel, saveSettings
				sf=file(buttonListener.saveFilepath, 'wb')
				json.dump(options, sf)
				sf.close()
			gd=GenericDialogPlus("FRETENATOR2: 2FRET2FURIOUSLY")
	
			gd.addMessage("Channel choices:")
			#user can pick which channel to base the segmentation on

			gd.setInsets(0,0,0)
			gd.addChoice("Segmentation", types, str(segmentChannel))


			gd.setInsets(0,0,0)
			gd.addMessage("""Segmentation settings:""")
			methods=["Otsu","Default", "Huang", "Intermodes", "IsoData", "IJ_IsoData", "Li", "MaxEntropy", "Mean", "MinError", "Minimum", "Moments", "Percentile", "RenyiEntropy", "Shanbhag", "Triangle", "Yen"]
			gd.setInsets(0,15,0)
			gd.addCheckbox("Difference of Gaussian instead of Gaussian?", DoG)
			gd.addSlider("Gaussian / DoG 1 sigma", 0.5, 10, gaussianSigma, 0.1)
			gd.addSlider("DoG 2 sigma", 0.5, 20, largeDoGSigma,0.1)
			gd.setModal(False)
			gd.addChoice("Autosegmentation method", methods, thresholdMethod)
			gd.addCheckbox("Manually set threshold? ", manualSegment)
			gd.addSlider("Manual threshold", 10, 65534, manualThreshold, 1)
			
			dilationOptions=["0", "1", "2","3", "4", "5", "6"]
			
			gd.addChoice("Dilation?", dilationOptions, str(dilation))
			gd.addCheckbox("Size exclusion of ROI? ", sizeExclude)
			gd.addSlider("Minimum ROI size", 0, 9999, minSize, 1)
			gd.addSlider("Maximum ROI size", 1, 10000, maxSize, 1)
			gd.addCheckbox("Watershed object splitting? ", watershed)
			
			gd.setInsets(0,0,0)
			gd.addMessage("""Analysis settings:""")
			backsubOpts=["Off", "Local label based", "Global mean"]
			intensities=["254", "4094", "65534"]
			gd.addChoice("Saturation removal threshold", intensities, str(maxIntensity))
			gd.addChoice("Background subtraction", backsubOpts, backsubOpts[backsubVal])
			gd.addCheckbox("""Use pixel by pixel analysis?""", pixelByPixel)
			gd.addCheckbox("Create nearest point projection? ", makeNearProj)
			gd.addCheckbox("Save settings as default? ", False)
			buttonListener = buttonClick()
			gd.addButton("Quickload nuclei settings", buttonListener)
			gd.addButton("Quickload diffuse settings", buttonListener)
			
			gd.addButton("Load settings file", buttonListener)
			gd.addButton("Save settings file", buttonListener)
			gd.setInsets(0,0,0)
			gd.addMessage("""For tutorials, click 'help'.
			
			https://doi.org/10.1007/978-1-0716-2297-1_17
			https://doi.org/10.1038/s41477-023-01447-4 """)
			gd.addHelp("https://github.com/JimageJ/FRETENATOR2/blob/main/README.md")
			gd.setLocation(0,0)

			gd.showDialog()
			buttonListener.bypass = 0

			choices=gd.getChoices()
			sliders=gd.getSliders()
			checkboxes=gd.getCheckboxes()
		segmentChannel=int(choices.get(0).getSelectedItem())
		thresholdMethod=choices.get(1).getSelectedItem()
		dilation=int(choices.get(2).getSelectedItem())
		maxIntensity=int(choices.get(3).getSelectedItem())
		
		gaussianSigma=sliders.get(0).getValue()/10.0
		largeDoGSigma = gd.sliders.get(1).getValue()/10.0
		manualThreshold=gd.sliders.get(2).getValue()
		minSize = gd.sliders.get(3).getValue()
		maxSize = gd.sliders.get(4).getValue()
		
		
		
		
		DoG = gd.checkboxes.get(0).getState()
		manualSegment = gd.checkboxes.get(1).getState()
		sizeExclude=gd.checkboxes.get(2).getState()
		watershed = gd.checkboxes.get(3).getState()
		#backSub = gd.checkboxes.get(4).getState()
		pixelByPixel = gd.checkboxes.get(4).getState()
		
		if (segmentChannelOld !=segmentChannel or
			thresholdMethodOld !=thresholdMethod or
			maxIntensityOld !=maxIntensity or
			gaussianSigmaOld !=gaussianSigma or
			largeDoGSigmaOld != largeDoGSigma or
			DoGOld !=DoG or
			manualSegmentOld != manualSegment or
			manualThresholdOld !=manualThreshold or
			dilation != dilationOld or
			sizeExcludeOld!=sizeExclude or
			minSizeOld!=minSize or
			maxSizeOld!=maxSize or
			watershedOld!=watershed
			):
				if minSizeOld!=minSize:
					if minSize>=maxSize:
						maxSize=minSize+1
						gd.sliders.get(5).setValue(maxSize)
				if maxSizeOld!=maxSize:
					if minSize>=maxSize:
						minSize=maxSize-1
						gd.sliders.get(4).setValue(minSize)
				if segmentChannelOld!=segmentChannel:
						clij2.clear()
						print('eh')
						segmentImp=extractChannel(imp1, segmentChannel, 0)
						gfx1=clij2.push(segmentImp)
						gfx2=clij2.create(gfx1)
						gfx3=clij2.create(gfx1)
						gfx4=clij2.create(gfx1)
						gfx5=clij2.create(gfx1)
						gfx7=clij2.create([imp.getWidth(), imp.getHeight()])
				gfx1,gfx2,gfx3,gfx4,gfx5 = segment(gfx1,gfx2,gfx3,gfx4,gfx5, gaussianSigma, thresholdMethod, maxIntensity, largeDoGSigma, pixelAspect, originalTitle, DoG, manualSegment, manualThreshold, dilation,sizeExclude, minSize, maxSize, watershed)
				clij2.maximumZProjection(gfx5, gfx7)
				labelPrevImp.close()
				labelPrevImp= clij2.pull(gfx7)
				IJ.setMinAndMax(labelPrevImp, 0,clij2.getMaximumOfAllPixels(gfx7))
				labelPrevImp.setTitle("Preview segmentation")
				labelPrevImp.show()
				
				IJ.run("glasbey_inverted")
		
		segmentChannelOld=segmentChannel
		thresholdMethodOld=thresholdMethod
		maxIntensityOld=maxIntensity
		gaussianSigmaOld=gaussianSigma
		largeDoGSigmaOld = largeDoGSigma
		DoGOld=DoG
		manualSegmentOld= manualSegment
		manualThresholdOld=manualThreshold
		dilationOld=dilation
		sizeExcludeOld=sizeExclude
		minSizeOld=minSize
		maxSizeOld=maxSize
		watershedOld=watershed
		buttonListener.bypass=0 
		buttonListener.quickload = 0
		buttonListener.save = 0 
		buttonListener.load = 0
		buttonListener.options = segmentChannel, donorChannel, acceptorChannel, acceptorChannel2, thresholdMethod, maxIntensity, gaussianSigma, largeDoGSigma, DoG,  manualSegment, manualThreshold, makeNearProj, dilation, sizeExclude, minSize, maxSize, watershed, backSub, pixelByPixel, saveSettings
		Thread.sleep(200)
	if gd.wasCanceled():
		clij2.clear()
		IJ.exit()
	
	labelPrevImp.close()
	makeNearProj = gd.checkboxes.get(5).getState()
	backSub= choices.get(4).getSelectedItem()
	saveSettings = gd.checkboxes.get(6).getState()
	backsubVal=0
	if backSub=="Local label based":
		backsubVal=1
	if backSub=="Global mean":
		backsubVal=2
	if gd.wasCanceled():
		clij2.clear()
		IJ.exit()
	return segmentChannel, donorChannel, acceptorChannel, acceptorChannel2, thresholdMethod, maxIntensity, gaussianSigma, largeDoGSigma, DoG,  manualSegment, manualThreshold, makeNearProj, dilation, sizeExclude, minSize, maxSize, watershed, backsubVal, pixelByPixel, saveSettings
	
def segment(gfx1,gfx2,gfx3,gfx4,gfx5, gaussianSigma, thresholdMethod, maxIntensity, largeDoGSigma, pixelAspect, originalTitle, DoG,  manualSegment, manualThreshold, dilation, sizeExclude, minSize, maxSize, watershed):
	"""Segmentation based on user settings"""
	# DoG filter for background normalised segmentation. NB. Kernel is Z-normalised to pixel aspect ratio
	if DoG == True :	
		clij2.differenceOfGaussian3D(gfx1, gfx2, gaussianSigma, gaussianSigma, 1+(gaussianSigma-1)/pixelAspect, largeDoGSigma, largeDoGSigma,largeDoGSigma/pixelAspect)
	else:
		clij2.gaussianBlur3D(gfx1, gfx2, gaussianSigma,gaussianSigma, 1+(gaussianSigma-1)/pixelAspect)

	if manualSegment == True :
		clij2.threshold(gfx2, gfx3, manualThreshold)
	else:
		#auto threshold and watershed to seed the object splitting
		clij2.automaticThreshold(gfx2, gfx3, thresholdMethod)

	if watershed:
		clij2.watershed(gfx3,gfx2)
	else:
		clij2.copy(gfx3,gfx2)
	
	# add watershed to original threshold, and then use this to generate a binary image of any ROI lost in watershed process
	clij2.addImages(gfx2, gfx3, gfx5)
	clij2.floodFillDiamond(gfx5, gfx4, 1, 2)
	clij2.replaceIntensity(gfx4, gfx5, 2, 0)
	
	#label watershed image
	clij2.connectedComponentsLabelingDiamond(gfx2, gfx4)

	#dilate all the labeled watershed ROI out (only onto zero labeled pixels), then multiply this by original binary map, to get labeled ROI
	clij2.dilateLabels(gfx4, gfx2, 6)

	clij2.multiplyImages(gfx2,gfx3, gfx4)
	
	#label the missed ROI then add on the largest value from the other labelled image (so they can be combined)
	watershedLabelMax =clij2.getMaximumOfAllPixels(gfx4)
	clij2.connectedComponentsLabelingDiamond(gfx5, gfx2)
	clij2.addImageAndScalar(gfx2, gfx5, (1 + watershedLabelMax))
	
	#delete the background and combine the two images
	clij2.replaceIntensity(gfx5, gfx2,(1 + watershedLabelMax), 0)
	clij2.maximumImages(gfx4, gfx2, gfx5)
	
	#remove labeled objects that are too big or too small
	if sizeExclude:
		clij2.excludeLabelsOutsideSizeRange(gfx5, gfx4, minSize, maxSize)
		clij2.copy(gfx4, gfx5)
		
	#dilate images
	if dilation:
		clij2.dilateLabels(gfx5, gfx4, dilation)
		clij2.copy(gfx4, gfx5)
	#gfx3 = threshold channel, gfx5 = label image, gfx1=original image, gfx2 & gfx4  are junk
	clij2.closeIndexGapsInLabelMap(gfx5,gfx4)
	return gfx1,gfx2,gfx3,gfx5, gfx4	


def Calculations(imp1, nFrame, table, gfx1, gfx2, gfx3, gfx4, gfx5, originalTitle,  pixelByPixel):
	"""Perform FRET calculations (Sorry! Complex function!)"""
	
	#Extract appropriate channels
	for i in range(imp1.getNChannels()):
		Imp=extractChannel(imp1, donorChannel, nFrame)
	
		#push donor and acceptor channels to gpu threshold
		gfx4=clij2.push(Imp)
		gfx6=clij2.create(gfx4)
		

		#thresholds to create a mask to remove saturated pixels from donor and acceptor images	
		clij2.threshold(gfx4,gfx6, maxIntensity)
		clij2.binarySubtract(gfx3, gfx2, gfx6)
		clij2.threshold(gfx5,gfx2, maxIntensity)
		clij2.binarySubtract(gfx6, gfx2, gfx3)
		clij2.threshold(gfx3,gfx6, 0.5)
	

	
	
	#Mask to remove the saturated pixels
	#donor is gfx2, acceptor FRET is gfx4, segment channel (acceptor normal) is gfx6, gfx3 is triple threshold
	clij2.mask(gfx4,gfx6, gfx2)
	clij2.mask(gfx5,gfx6, gfx4)
			
	#NB have to push the acceptor image now...
	gfx6=clij2.push(acceptorImp2)

	#extract the intensity of each nucleus for each channel
	results=ResultsTable()
	clij2.statisticsOfBackgroundAndLabelledPixels(gfx2, gfx1, results)
	donorChIntensity=results.getColumn(13)
	results2=ResultsTable()
	clij2.statisticsOfBackgroundAndLabelledPixels(gfx4, gfx1, results2)
	acceptorChIntensity=results2.getColumn(13)
	results3=ResultsTable()
	clij2.statisticsOfBackgroundAndLabelledPixels(gfx6, gfx1, results3)
	
	#calculate the fret ratios, removing any ROI with intensity of zero
	FRET =[]
	
	for i in xrange(len(acceptorChIntensity)):
		if (acceptorChIntensity[i]>0) and (donorChIntensity[i]>0):
			#don't write in the zeros to the results
			FRET.append((1000*acceptorChIntensity[i]/donorChIntensity[i]))
			table.incrementCounter()
			#frame, label and ER
			table.addValue("Frame (Time)", nFrame)
			table.addValue("Label", i)
			table.addValue("Emission ratio", acceptorChIntensity[i]/donorChIntensity[i])

			#mean emission
			table.addValue("Mean donor emission", results.getValue("MEAN_INTENSITY",i))
			table.addValue("Mean acceptor emission (FRET)", results2.getValue("MEAN_INTENSITY",i))
			table.addValue("Mean acceptor emission", results3.getValue("MEAN_INTENSITY",i))
			
			#sum emission
			table.addValue("Sum donor emission", donorChIntensity[i])
			table.addValue("Sum acceptor emission (FRET)", acceptorChIntensity[i])
			table.addValue("Sum acceptor emission", results3.getValue("SUM_INTENSITY",i))
			#shape and location descriptors
			table.addValue("Volume", cal.pixelWidth * cal.pixelHeight * cal.pixelDepth * results.getValue("PIXEL_COUNT",i))
			table.addValue("Pixel count", results.getValue("PIXEL_COUNT",i))
			table.addValue("x", cal.pixelWidth*results.getValue("CENTROID_X",i))
			table.addValue("y", cal.pixelHeight*results.getValue("CENTROID_Y",i))
			table.addValue("z", cal.pixelDepth*results.getValue("CENTROID_Z",i))
			
			#File name for traceability 
			table.addValue("File name", originalTitle)
		else:
			#must write in the zeros as this array is used to generate the map of emission ratios
			FRET.append(0)
			
	table.show("Results of " + originalTitle)
	#export label image
	labelImp = clij2.pull(gfx1)
	if pixelByPixel==0:
		#write all the emission ratios to an array, push to an GFX image, use this to map emission ratios
		FRET[0]=0
		FRETarray= array( "f", FRET)
		fp= FloatProcessor(len(FRET), 1, FRETarray, None)
		FRETImp= ImagePlus("FRETImp", fp)
		gfx4=clij2.push(FRETImp)
		clij2.replaceIntensities(gfx1, gfx4, gfx5)
		maxProj=clij2.create(gfx5.getWidth(), gfx5.getHeight(), 1)
		clij2.maximumZProjection(gfx5, maxProj)
		
		
		#pull the images
		FRETimp2=clij2.pull(gfx5)
		FRETProjImp=clij2.pull(maxProj)
		
	else:
		
		#donor is gfx2, acceptor FRET is gfx4, segment channel (acceptor normal) is gfx6, threshold image:gfx3
		#blur the acceptor and donor channels, then remask to prevent NaN/infinite values after division
		clij2.gaussianBlur3D(gfx4, gfx5, 1.1, 1.1, 1.1)
		clij2.mask(gfx5, gfx1, gfx4)
		clij2.gaussianBlur3D(gfx2, gfx5, 1.1, 1.1, 1.1)
		clij2.mask(gfx5, gfx1, gfx2)
		
		#create Z sum projected donor and acceptor images for Z-proj ratio calc -> may replace with a different technique later
		donorSum=clij2.create(gfx4.getWidth(), gfx4.getHeight(), 1)
		acceptorFSum=clij2.create(gfx4.getWidth(), gfx4.getHeight(), 1)
		clij2.sumZProjection(gfx2, donorSum)
		clij2.sumZProjection(gfx4, acceptorFSum)
		
		#Divide Z proj Acceptor by Z proj Donor and pull image
		maxProj=clij2.create(gfx4.getWidth(), gfx4.getHeight(), 1)
		clij2.divideImages(acceptorFSum, donorSum, maxProj)
		clij2.multiplyImageAndScalar(maxProj, donorSum,1000)
		FRETProjImp=clij2.pull(donorSum)
		
		#pull acceptor and donor stacks to convert to 32 bit
		acceptorImp=clij2.pull(gfx4)
		donorImp=clij2.pull(gfx2)
		ImageConverter(acceptorImp).convertToGray32()
		ImageConverter(donorImp).convertToGray32()
		
		#clean up GPU memory
		clij2.clear()
		
		#push 32bit images and perform ratio calc
		gfx4=clij2.push(acceptorImp)
		gfx2=clij2.push(donorImp)
		gfx1=clij2.create(gfx2.getDimensions(), clij2.Float)
		clij2.divideImages(gfx4, gfx2, gfx1)
		clij2.multiplyImageAndScalar(gfx1, gfx2,1000)
		#pull ratio stack
		FRETimp2=clij2.pull(gfx1)

	#clean up
	clij2.clear()
	donorImp.close()
	acceptorImp.close()
	acceptorImp2.close()
	
	return table, FRETimp2, FRETProjImp, labelImp
	
	
def nearestZProject(imp1):
	relicedImp=Slicer().reslice(imp1)
	relicedStack=relicedImp.getStack()
	width=imp1.getWidth()
	height=imp1.getHeight()
	depth=imp1.getNSlices()
	
	topPixels=zeros('f', width * height)  
	
	stack2=ImageStack( width, height)
	for i in range(1,relicedImp.getNSlices()):
		pixels= relicedStack.getPixels(i)
		for x in xrange(width):
			for pixel in xrange(x, x+width*(depth-1),width):
				#after finding the first pixel above the threshold value, add the value to the list
				if pixels[pixel] != 0:
				
					topPixels[i*width+x]=pixels[pixel]
					#break from looping the y when 1st threshold pixel is found is met -> increases speed drastically! Otherwise need an if statement every loop...
					break
	
	ip2=FloatProcessor(width, height, topPixels, None)
	imp2=ImagePlus("Nearest point proj",ip2)
	imp3= imp2.resize(imp2.getWidth()*2, imp2.getHeight()*2, 'none')
	return imp3



def outline(imp3, originalTitle):
	
	
	#clij implementation -thicker lines
	"""
	src=clij2.push(imp3)
	dst=clij2.create(src)
	dst2=clij2.create(src)
	clij2.detectLabelEdges(src,dst)
	clij2.binaryNot(dst,dst2)
	clij2.multiplyImages(src, dst2, dst)
	imp4=clij2.pull(dst)
	imp4.show()
	clij2.clear()
	"""
	
	imp2=imp3.duplicate()
	stack1=imp3.getStack()
	width=imp3.getWidth()
	height=imp3.getHeight()
	stack2=ImageStack(width, height)
	pixlist=[]
	
	for i in range(imp3.getNSlices()):
		pixlist=[]
		pixels1=stack1.getPixels(i+1)
		#if pixel is different to the pixel to the left or above, set it to 0
		pixels2=map(lambda j: pixels1[j] if pixels1[j]-pixels1[j-1]==0 and pixels1[j]-pixels1[j-width]==0 else 0, xrange(len(pixels1)))
		processor=FloatProcessor(width, height, pixels2, None)
		stack2.addSlice(processor)
	imp2=ImagePlus("Nearest point emission ratios of "+ originalTitle, stack2)
	imp2.show()
	return imp2


# *****************************body of code starts****************************************


#give install instructions for CLIJ if not installed

try: 
	from net.haesleinhuepf.clij2 import CLIJ2

except:
	errorDialog("""This plugin requires clij2 to function. 
	
	To install please follow these instructions: 
	
	1. Click Help>Update> Manage update sites
	2. Make sure the "clij2" update site is selected.
	3. Click Close> Apply changes.
	4. Close and reopen ImageJ""")


clij2 = CLIJ2.getInstance()
clij2.clear()


#get the current image
imp1= IJ.getImage()


if exists('FRETENATOR2SegSettings.json'):
	print('exists')
	sf=file('FRETENATOR2SegSettings.json', 'rb')
	options=json.load(sf)
	sf.close()
	print(options)
else:
	options=(3, 1, 2, 3, 'Otsu', 65534, 0.8, 4.0, True, False, 3000, True, 0, False, 10, 10000, True, 0, False, True)
	
options= previewDialog(imp1, options)

#get the pixel aspect for use in zscaling kernels for filters
cal = imp1.getCalibration()
pixelAspect=(cal.pixelDepth/cal.pixelWidth)
originalTitle=imp1.getTitle()



IJ.log(originalTitle +" settings:")
IJ.log("segmentChannel, donorChannel, acceptorChannel, acceptorChannel2, thresholdMethod, maxIntensity, gaussianSigma, largeDoGSigma, DoG, manualSegment, manualThreshold, makeNearProj, dilation, sizeExclude, minSize, maxSize, watershed, backSub:")
IJ.log(str(options))
segmentChannel, donorChannel, acceptorChannel, acceptorChannel2, thresholdMethod, maxIntensity, gaussianSigma, largeDoGSigma, DoG,  manualSegment, manualThreshold, makeNearProj, dilation, sizeExclude, minSize, maxSize, watershed, backsubVal, pixelByPixel, saveSettings  =options
if saveSettings==1:
	sf=file('FRETENATOR2SegSettings.json', 'wb')
	json.dump(options, sf)
	sf.close()
if pixelByPixel==1:
	makeNearProj =0
totalFrames=imp1.getNFrames() +1

#table is the final results table
table = ResultsTable()

clij2 = CLIJ2.getInstance()
clij2.clear()



conThresholdStack=ImageStack(imp1.width, imp1.height)
conlabelImpStack=ImageStack(imp1.width, imp1.height)
for nFrame in xrange(1, totalFrames):
	clij2.clear()
	segmentImp=extractChannel(imp1, segmentChannel, nFrame)
	gfx1=clij2.push(segmentImp)
	gfx2=clij2.create(gfx1)
	gfx3=clij2.create(gfx1)
	gfx4=clij2.create(gfx1)
	gfx5=clij2.create(gfx1)
	
	gfx1,gfx2,gfx3,gfx4,gfx5 = segment(gfx1,gfx2,gfx3,gfx4,gfx5, gaussianSigma, thresholdMethod,maxIntensity, largeDoGSigma, pixelAspect, originalTitle, DoG,manualSegment, manualThreshold, dilation,sizeExclude, minSize, maxSize, watershed)
	#gfx1 = original image, gfx3=threshold image, gfx4= label map


	labelImp = clij2.pull(gfx4)
	IJ.setMinAndMax(labelImp, 0,255)
	labelImp.setCalibration(cal)
	labelImp.setTitle("Label map of "+originalTitle)
	

	labelImp.show()
	IJ.run("glasbey_inverted")	
	ltable=ResultsTable()		
	table=ResultsTable()
	clij2.statisticsOfBackgroundAndLabelledPixels(gfx4,gfx4, ltable)
	for i in range(ltable.size()-1):
		table.setValue("Label",i, ltable.getValue("MEAN_INTENSITY",i))
		table.setValue("Volume",i, cal.pixelWidth * cal.pixelHeight * cal.pixelDepth * ltable.getValue("PIXEL_COUNT",i))
		table.setValue("Pixel count",i, ltable.getValue("PIXEL_COUNT",i))
		table.setValue("x",i, cal.pixelWidth*ltable.getValue("CENTROID_X",i))
		table.setValue("y",i, cal.pixelHeight*ltable.getValue("CENTROID_Y",i))
		table.setValue("z",i, cal.pixelDepth*ltable.getValue("CENTROID_Z",i))
	
	

	for nChannel in xrange(1, imp1.getNChannels()):
		chtable=ResultsTable()
		gfx2.close()
		gfx2=clij2.push(extractChannel(imp1, nChannel, nFrame))
		if backsubVal==1:
			gfx7=clij2.create(gfx1)
			gfx3=createSubtractionLabels(gfx4, gfx5 , gfx3)
			gfx2=localLabelBackSub(gfx3, gfx4, gfx2, gfx1, gfx5) #	"""Requires a dilatedlabelGFX image, labelGFX, a quantGFX image, and two sacrificial otherGFX images"""
			
		if backsubVal==2:
			gfx2=globalBackSub(gfx4, gfx2, gfx1) # 	"""Requires a labelGFX image, the quantGFX to be quantified (also the output) and, one sacrificial otherGFX images"""
		clij2.statisticsOfBackgroundAndLabelledPixels(gfx2,gfx4, chtable)

		for i in range(chtable.size()-1):
			table.setValue("Mean Intensity Ch" +str(nChannel) ,i, chtable.getValue("MEAN_INTENSITY",i))
			table.setValue("Min Intensity Ch" +str(nChannel) ,i, chtable.getValue("MINIMUM_INTENSITY",i))
			table.setValue("Max Intensity Ch" +str(nChannel) ,i, chtable.getValue("MAXIMUM_INTENSITY",i))
			table.setValue("Sum Intensity Ch" +str(nChannel) ,i, chtable.getValue("SUM_INTENSITY",i))
#chtable.show('test2')
table.show('Results')

print backsubVal
clij2.clear()