Live FFT color map

'''
Code to process the live-view image in GMS, producing a map of crystalline regions
This requires a module FFTArrayAnalysis, which requires scipy, skimage, matplotlib, and tqdm 
The module is part of a package which can be installed with pip:
    pip install BenMillerScripts

WARNING: Due to a bug in GMS 3.5.0 and 3.5.1, this script will not run in those versions

Run the code with the live-view image, containing a rectangular ROI, front-most in GMS
    The code also works with an IS video played back with the IS player
    To stop calculation, delete the ROI. 
New frames are computed as often as possible.
Lines of code between #XXXXXXXX... lines can be changed to modify the data processing
    All other lines of code are general, and can be re-used to produce other kinds of processed maps from a live-view image
#Code written by Ben Miller. Last Updated May 2022
'''
if not DM.IsScriptOnMainThread(): print('Scipy scripts cannot be run on Background Thread.'); exit()
from benmillerscripts import FFTArrayAnalysis as FAA
import time 
import numpy as np
import traceback

#User editable variables are set here
#XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
#Set FFT Size and Spacing
FAA.FFTsize = 128 #(Default 128)
FAA.spacing = 32 #(Default 32)
map_var = "theta" #Parameter to be mapped ("theta" or "radius") (default "theta")
FAA.maskP = 10 #Set Percentage of Radial Profile (from center) to be Ignored (Default 10)
FAA.maskC_width = 2 #Set to an integer value to also mask the center vertical and horizontal lines in the FFT, set to 0 to not mask (Default 2)
FAA.binning = 1 #Set how much to bin the raw image data prior to computing FFTs higher values save time, but may result in loss of information (Default 1)
Scale = 2 #Upscaling factor for the FFT Map. Either an integer >0 or "match" to output same size as input (Default 2)
#a scale value of "match" is required for the overlay option
Overlay=False #Set whether to overlay the original image intensity on top of the color map (must also set Scale to "match") (Default False)
show_cube = False #This should only be used for troubleshooting. This returns a 4D cube of FFTs, which can be processed like a 4DSTEM dataset in GMS (Default False)
Display_Initial_Diffractogram = True #If True, this displays a diffractogram of the first frame. This can be useful for determining whether the binning, maskC_width, and maskP are suitable. (Default True)
display_scale = True #If True, this displays the color scale for the map, indicating which direction corresponds to each color.
print_timing = True      # (Default True) Select whether to output the time it takes to compute each frame
#XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

class CListen(DM.Py_ScriptObject):
    #Function to find an ROI placed on an image by the user, returning the ROI ID. 
    #If no ROI found, create a new one covering the entire image.
    def find_ROI(self,image):
        imageDisplay = image.GetImageDisplay(0)
        numROIs = imageDisplay.CountROIs()
        id = None
        for n in range(numROIs):
            roi = imageDisplay.GetROI(n) 
            if roi.IsRectangle():
                roi.SetVolatile(False)
                roi.SetResizable(False)  
                id = roi.GetID()
                break
        if id is None:
            #If No ROI is found, create one that covers the whole image. 
            print("\nRectangular ROI not found... using whole image")
            data_shape = image.GetNumArray().shape
            roi=DM.NewROI()
            roi.SetRectangle(0, 0, data_shape[0], data_shape[1])
            imageDisplay.AddROI(roi)
            roi.SetVolatile(False)
            roi.SetResizable(False) 
            id = roi.GetID()
        return id

    #Initialization Function
    def __init__(self,img):
        try:
            #Create an index that is incremented each time data is processed.
            self.i = 0
            #get the original image and assign it to self.imgref
            self.imgref = img
            #Get the data from the region within an ROI
            self.roi = DM.GetROIFromID(self.find_ROI(self.imgref))
            val, val2, val3, val4 = self.roi.GetRectangle()
            self.data = self.imgref.GetNumArray()[int(val):int(val3),int(val2):int(val4)]
            #get the shape of the original image
            (input_sizex, input_sizey) = self.data.shape
            x_origin = self.imgref.GetDimensionOrigin(0); y_origin = self.imgref.GetDimensionOrigin(1)
            x_scale = self.imgref.GetDimensionScale(0);   y_scale = self.imgref.GetDimensionScale(1)
            #scale unit of microns causes problems for python in DM
            try: 
                x_unit = self.imgref.GetDimensionUnitString(0) ; y_unit = self.imgref.GetDimensionUnitString(1)
            except: 
                x_unit = self.imgref.GetDimensionUnitInfo(0)[0]; y_unit = self.imgref.GetDimensionUnitInfo(1)[0]
                if x_unit == 'micrometer': x_unit = 'um'        
                if y_unit == 'micrometer': y_unit = 'um'
            #XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
            #First plot the color scale
            if display_scale: FAA.DrawColorScale(map_var = map_var)
            # use the FAA module to calculate the shape of the map which will be produced
            result_size_x,result_size_y = FAA.FFT_Result_Shape(self.data, scale=Scale)
            #Use the FAA module to create a new color image to contain the results of processing (it is currently difficult to work with color images in GMS using Python)
            self.result_image = FAA.CreateDM_RGB(np.zeros((result_size_x,result_size_y,3)))
            #Process the first frame
            (processedframe,direction_image,intensity_image,spacing_image,diffractogram_max) = FAA.processframe(
                                    self.imgref,map_var,show_cube, scale=Scale, overlay=Overlay, im_data=self.data)
            #In addition to the live color map, optionally produce and display a static maximum diffractogram image from the first frame
            if Display_Initial_Diffractogram: 
                diffim = DM.CreateImage(np.copy(diffractogram_max))
                diffim.SetDimensionCalibration(0,-1/(2*(x_scale)),1/(x_scale)/diffractogram_max.shape[0],x_unit+"-1",0)
                diffim.SetDimensionCalibration(1,-1/(2*(y_scale)),1/(y_scale)/diffractogram_max.shape[1],y_unit+"-1",0)
                diffim.SetName("Initial Maximum Diffractogram of "+img.GetName())
                diffim.ShowImage()
                del diffim
            #Initialize the color map intensity scale maximum
            self.RGB_scale_max  = np.percentile(intensity_image, 95)
            #Set the calibration based on the original data
            self.result_image.SetDimensionCalibration(0,x_origin,x_scale*input_sizex/result_size_x,x_unit,0)
            self.result_image.SetDimensionCalibration(1,y_origin,y_scale*input_sizey/result_size_y,y_unit,0)
            #Display the result image in GMS
            self.result_image.ShowImage()
            #Set the image name which will be displayed in the image window's title bar
            self.result_image.SetName("FFT spot " + map_var +" map of "+img.GetName())
            #Set the display zoom of the map, so it is not too small
            zoom = 4
            DM.ExecuteScriptString("GetFrontImage().setzoom("+str(zoom)+")\nImageDocumentOptimizeWindow(GetFrontImageDocument())")
            #Get the ID of the result image
            self.result_ID = self.result_image.GetID()
            #Use the FAA module to update the color map (it is currently difficult to work with color images in GMS using Python)
            FAA.UpdateRGB(self.result_ID, processedframe)
            #XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
            #Copy tags from original to processed image
            Tag_Copy(self.imgref, self.result_image,subPath='OriginalImageTags')
            DM.Py_ScriptObject.__init__(self)
            self.stop = 0
        except: print(traceback.format_exc())
    
    #This function is run each time the image changes
    def HandleDataChangedEvent(self, flags, image):
        try:
            if not self.stop:
                #start timing
                start=time.perf_counter() 
                #Get an (updated) ROI position
                val, val2, val3, val4 = self.roi.GetRectangle()
                #Get the data from the ROI area as a numpy array
                self.data = self.imgref.GetNumArray()[int(val):int(val3),int(val2):int(val4)]                
                #XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
                #compute color map from the image, and place this profile into results image
                (processedframe,direction_image,intensity_image,spacing_image,diffractogram_max) = FAA.processframe(
                                    self.imgref,map_var,show_cube, scale=Scale, overlay=Overlay,im_data = self.data, RGB_scale_max=self.RGB_scale_max )
                #Update the color map intensity scale maximum for use with the next frame
                self.RGB_scale_max  = max(self.RGB_scale_max,np.percentile(intensity_image, 95))
                #Update the color map with the latest data
                FAA.UpdateRGB(self.result_ID, processedframe)
                #XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
                #end timing and output time to process this frame
                end=time.perf_counter()
                if print_timing: print("Processed Image "+str(self.i)+"  Processing Time= "+str(end-start))                                         
                #Increment an index each time data is processed. 
                self.i = self.i+1
        except: print(traceback.format_exc())

    #Function to Delete Image Listener
    def __del__(self):
        print("Listener Deleted")
        DM.Py_ScriptObject.__del__(self)
        
    #Function to end processing by deleting or unregistering listener
    def RemoveListeners(self):
        try: 
            if not self.stop:
                self.stop = 1
                DM.DoEvents()
                global listener
                #DM 3.5.2 and higher have new function for unregistering listeners. 
                #DM 3.4.3 and lower should delete the listener instead
                #DM 3.5.0 and 3.5.1 have a fatal flaw regarding listeners, so this script is not compatible with those versions
                if (get_DM_version()[1][0] == "4" and get_DM_version()[0] == "3"):
                    del listener
                else:
                    listener.UnregisterAllListeners()
                print("Live Processing Script Ended")
        except: print(traceback.format_exc())
        
    #Remove listeners if source image window is closed
    def HandleWindowClosedEvent(self, event_flags, window):
        print("Window Closed")
        self.RemoveListeners()
                
    #Remove listeners if the ROI is deleted
    def HandleROIRemovedEvent(self, img_disp_event_flags, img_disp, roi_change_flag, roi_disp_change_flags, roi):
        print("ROI Removed")
        self.RemoveListeners()
        
#Function to get the currently used version of DigitalMicrograph 
def get_DM_version():
    #No Python script command exists to get the DM version, 
    #so we first run a DM script to put the values in the global tags
    dm_script = ('number minor, major, bugVersion\n'
        'GetApplicationVersion(major, minor, bugVersion)\n'
        'GetPersistentTagGroup().TagGroupSetTagAsLong("Python_Temp:DM_Version_Major",major)\n'
        'GetPersistentTagGroup().TagGroupSetTagAsLong("Python_Temp:DM_Version_Minor",minor)\n'
        'GetPersistentTagGroup().TagGroupSetTagAsLong("Python_Temp:DM_Version_bugVersion",bugVersion)')
    DM.ExecuteScriptString(dm_script)
    #Now get the information stored in the global tags by the DM script
    version = [0,0,0]
    _,version[0] = DM.GetPersistentTagGroup().GetTagAsString("Python_Temp:DM_Version_Major")
    _,version[1] = DM.GetPersistentTagGroup().GetTagAsString("Python_Temp:DM_Version_Minor")
    _,version[2] = DM.GetPersistentTagGroup().GetTagAsString("Python_Temp:DM_Version_bugVersion")
    return version  
          
def Tag_Copy(image_source, image_dest, subPath = None ):
    '''
    Copy all tags between source and destination.
    If no destination subPath is provided, the destination tags will be replaced.
    '''
    #Copy Tags
    tg_source = image_source.GetTagGroup()
    tg_dest = image_dest.GetTagGroup()
    if ( subPath != None ):
        tg_dest.SetTagAsTagGroup(subPath,tg_source.Clone())
    else:
        tg_dest.DeleteAllTags()
        tg_dest.CopyTagsFrom(tg_source.Clone())

#Main Code Starts Here

#Check that we are not running 3.5.0 or 3.5.1 which have a known bug affecting this script.
if (((get_DM_version()[1] == '51') or (get_DM_version()[1] == '50')) and get_DM_version()[0] == "3"):
    DM.OkDialog("Due to a bug in DigitalMicrograph 3.5.0 and 3.5.1, this script would cause DM to crash in those versions. \n\nScript Aborted.")
    exit()

#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()
#Get the image display, for the ROI-removed listener
imageDisplay = img1.GetImageDisplay(0)

#Listeners are started here
#initiate the image listener
listener = CListen(img1)
#check if the source window closes
WindowClosedListenerID = listener.WindowHandleWindowClosedEvent(imDocWin, 'pythonplugin')
#check if the ROI has been deleted
ROIRemovedListenerID = listener.ImageDisplayHandleROIRemovedEvent(imageDisplay,'pythonplugin')
#check if the source image changes
DataChangedListenerID = listener.ImageHandleDataChangedEvent(img1, 'pythonplugin')
#IDs are not used in this script, but could be used to unregister individual listeners in DM 3.5.2 and higher.