#Code to compute ROI-based profiles over time from the live-view image in GMS. 
#Run the code with the live-view image front-most in GMS. 
#To stop calculation, close the live view image. 
#Lines of code between #XXXXXXXX... lines are specific to computing the ROI-based profile
	#All other lines of code are general, and can be re-used to produce other kinds of profiles from a live-view image or its FFT

#Requires Scipy. 
#To install packages like scipy, see instructions in GMS Help:Python:Installation and Configuration:Additional Packages

#Code written by Ben Miller. Last Updated 2020-06

import numpy as np
import time
import sys
if (DM.IsScriptOnMainThread() == False):
	print( ' MatplotLib and scipy scripts require to be run on the main thread.',
		'\n Uncheck the "Execute on Background Thread"',
		'checkbox at the bottom of the Script Window' )
	exit()
import scipy
from scipy import ndimage
from scipy import signal
from scipy import fftpack
from scipy.ndimage.interpolation import geometric_transform
from numpy import linalg as LA

integration_width = 100 #width of profile (perpendicular to profile direction) in pixels
initial_result_image_width = 100 #how many profiles can be displayed in the intial result window (window is automatically expanded as needed)
time_interval_mult=1 #multiple for the median time interval between frames (set higher to ensure time between each frame is exactly identical)
N_init=9 #number of intial frames to use to estimate median time interval

class CListen(DM.Py_ScriptObject):

	#Initialization Function
	def __init__(self, img):
		self.name = "ImageO"
		self.i = 0
		self.sleeptime = 0
		self.N_init=N_init
		self.profiletime=0
		self.time_outliers=[]
		self.duration=[]
		self.j = True #this is only for IS player testing
		#get the original image and assign it to self.imgref
		self.imgref = img
		#get the shape and calibration of the original image
		(input_sizex, input_sizey) = img.GetNumArray().shape
		origin, self.scale_orig, self.unit_stringo = self.imgref.GetDimensionCalibration(0, 0)
		#XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
		#Get a line ROI from the image
		self.line_ROI=self.get_line_ROI_coords(img)
		# Process the first image
		processedimagedata = self.processimage(img.GetNumArray(), self.line_ROI,integration_width)
		#Set the size and calibration of the result image
		(r_img_size,) = processedimagedata.shape
		(self.scale,self.unit_string) = (self.scale_orig, self.unit_stringo)
		#XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
		#create empty set for result images
		self.result_images = {}
		#create 1st result image and set calibration
		self.result_images[self.name] = DM.CreateImage(np.zeros((r_img_size,initial_result_image_width)))
		self.result_images[self.name].SetDimensionCalibration(1,0,self.scale,self.unit_string,0)
		self.result_images[self.name].ShowImage()		
		#get numpy array from result image
		self.result_array = self.result_images[self.name].GetNumArray()
		DM.Py_ScriptObject.__init__(self)
		self.end=time.perf_counter()
		self.stop = 0
	
	#Function to end Image Listener
	def __del__(self):
		DM.Py_ScriptObject.__del__(self)
		
	#XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX	
	#Function to Get Line ROI Coordinates from an image
	def get_line_ROI_coords(self,img):
		roi_disp = img.GetImageDisplay(0)
		num_rois = roi_disp.CountROIs()
		#If an ROI is selected, get this one, otherwise get 1st ROI
		for r in range(num_rois):
			if roi_disp.IsROISelected(roi_disp.GetROI(r)):
				roi = roi_disp.GetROI(r)
				break
			else:
				roi = roi_disp.GetROI(0)
		#If line ROI is found, get coordinates, otherwise prompt user and stop script
		try:
			#left, top, right, bottom = roi.GetLine()
			line_ROI=roi.GetLine()
		except:
			DM.OkDialog( 'Error: No ROI Found. \n\nBefore Running this Script, Select an Image with an ROI' )
			print('\nError: No ROI Found. Script Stopped. \nBefore Running this Script, Select an Image with an ROI')
			sys.exit()
		return line_ROI
		
	# Function to Process the Image Data in Each Image and Produce a Vector (integration width =1)
	#(This function is no longer used)
	def processimage_old(self, numpy_data, lr):
		#Produce a profile line from an image and 4 coordinates
		x_y_length=[lr[2]-lr[0],lr[3]-lr[1]]
		num=np.trunc(LA.norm(x_y_length)).astype('int')
		#Produce line coordinates
		x, y = np.linspace(lr[0], lr[2], num), np.linspace(lr[1], lr[3], num)
		coords=np.vstack((y,x))
		# Extract the values along the line, using cubic interpolation
		processed_output=scipy.ndimage.map_coordinates(numpy_data, coords)
		return(processed_output)
		
	# Function to Process the Image Data in Each Image and Produce a Vector	
	def processimage(self, numpy_data, lr, integration_width):
		#Get Coordinates function
		def get_coordinates(start, end, n):
			length_f=LA.norm([end[0] - start[0],end[1] - start[1]])
			num_samples = int(np.trunc(length_f))
			a = np.linspace(0, num_samples - 1, num_samples)  # along line
			t = np.linspace(-(n-1)*0.5, (n-1)*0.5, n)  # orthogonal to line
			x, y = np.meshgrid(a, t)
			dy = (end[0] - start[0]) / num_samples
			dx = (end[1] - start[1]) / num_samples
			yy = start[0] + dy * x + dx * y
			xx = start[1] + dx * x - dy * y
			return yy, xx
		yy, xx = get_coordinates(lr[:2], lr[2:], int(np.trunc(integration_width)))
		interpdata = scipy.ndimage.map_coordinates(numpy_data, (xx, yy))
		if len(interpdata.shape) > 1:
			return np.sum(interpdata, 0)/interpdata.shape[0]
		else:
			return interpdata
	#XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX	

	
	#This function is run each time the image changes
	def HandleDataChangedEvent(self, flags, image):
		if not self.stop:
			if self.j: #this is only for IS player testing
				#start timing
				start=self.end
				self.j = not self.j #this is only for IS player testing
				(result_sizey, result_sizex) = self.result_array.shape
				#if the result image is nearly full, make it 2x larger
				if self.i > result_sizex-2:
					#create a new numpy array 2x larger
					self.result_array_temp = np.append(self.result_array, np.zeros_like(self.result_array), axis=1)
					#close the old results image in DM
					DM.DeleteImage(self.result_images[self.name])
					print("ImageDeleted")
					#create a new results image and calibrate it
					self.name="Image{0}".format(self.i)
					self.result_images[self.name] = DM.CreateImage(np.copy((self.result_array_temp)))
					self.result_images[self.name].SetDimensionCalibration(1,0,self.scale,self.unit_string,0)
					print("ImageCreated")
					#display new result image in DM
					self.result_images[self.name].ShowImage()	
					#get numpy array from new result image
					self.result_array = (self.result_images[self.name].GetNumArray())
				#XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
				#compute ROI profile from the image, and place this profile into results image
				self.result_array[:,self.i] = self.processimage(self.imgref.GetNumArray(), self.line_ROI,integration_width)
				#XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
				self.result_images[self.name].UpdateImage()
				#end timing and output time to process this frame
				self.end=time.perf_counter()
				self.duration.append(self.end-start-self.sleeptime)
				print("Processed Image "+str(self.i)+"  Processzing Time= "+str(self.duration[-1]))
				if self.i==self.N_init: self.profiletime = time_interval_mult*np.median(np.asarray(self.duration))
				self.sleeptime = max(0,self.profiletime-self.duration[-1])
				if self.sleeptime == 0: self.time_outliers.append((self.i,self.profiletime-self.duration[-1]))
				DM.Sleep(self.sleeptime)
				self.i = self.i+1
			else: 
				self.j = not self.j #this is only for IS player testing
	#Function to end script if source image window is closed
	def HandleWindowClosedEvent(self, event_flags, window):
		print("Window Closed")
		self.stop = 1
		self.__del__()
		print("Script Ended")
		DM.Sleep(2)
		print(self.time_outliers[self.N_init:])
		#del self.imgref
				

#Main Code Starts Here
#Get front image in GMS		
img1 = DM.GetFrontImage()
#Get the image window, so we can check if it gets closed
imageDoc = DM.GetFrontImageDocument()
imDocWin = imageDoc.GetWindow()
#initiate the image listener
listener = CListen(img1)
#check if the source window closes
listener.WindowHandleWindowClosedEvent(imDocWin, 'pythonplugin')
#check if the source image changes
listener.ImageHandleDataChangedEvent(img1, 'pythonplugin')