Live radial profile
Python Script
Code to compute radial-max profiles of the FFT of a live-view image in DigitalMicrograph, producing a 2D display of these profiles over time.
Preview
'''
Code to compute radial-max profiles of the FFT of live-view image in GMS,
producing a 2D display of these profiles over time.
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.
Profiles are computed as often as possible.
Lines of code between #XXXXXXXX... lines are specific to computing a radial-max profile
All other lines of code are general, and can be re-used to produce other kinds of profiles from a live-view image
Code written by Ben Miller. Last Updated Apr 2022
'''
import time
import numpy as np
import traceback
if not DM.IsScriptOnMainThread(): print('Scipy scripts cannot be run on Background Thread.'); exit()
import scipy
from scipy import ndimage
from scipy import signal
from scipy import fftpack
from scipy.ndimage.interpolation import geometric_transform
#User editable variables are set here
#XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
median = 0 #set to 1 to apply median filter to the FFT prior to profile creation (slows calculation, especially for large input images)
initial_result_image_width = 200 #how many profiles can be displayed in the intial result window (window is automatically expanded as needed)
profile_result_length_ratio = 4 #Set to some integer 2^N, N=>0. Smaller N will make calculation slower. Default: 4
profile_angular_sampling_resolution = 256 #set how many samples are taken around the circumference of the radial profile
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
#XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
#Function to Determine Result Profile Length
def profile_length(self,sx,sf):
return int(sx/sf)
#Function to Set Result Image Calibration
def calibrate_result(self,scale,sx,sf,ustr):
diff_scale = sf/scale/2/sx
unit_str = ustr+"-1"
return(diff_scale,unit_str)
#Funtion to convert cartesian-coordinate image to polar-coordinate image
def topolar(self, img, r_size, theta_size, order=1):
sx, sy = img.shape
max_radius = int(sx/2)
#define transform
def transform(coords):
theta = 2.0*np.pi*coords[1] / (theta_size - 1.)
radius = max_radius * coords[0] / r_size
i = int(sx/2) - radius*np.sin(theta)
j = radius*np.cos(theta) + int(sx/2)
return i,j
#perform transform
polar = geometric_transform(img, transform, output_shape=(r_size,theta_size), order=order,mode='constant',cval=1.0,prefilter=False)
return polar
#Function to calculate radial profile of FFT from image
def FFT_radial_profile(self, image_o, profile_ang_res,length_ratio,do_median):
#compute FFT
fft_im = np.absolute(scipy.fftpack.fftshift(np.fft.fft2(image_o)))
#Median-Filter FFT to remove single-pixel outliers
if do_median: fft_im_median = scipy.ndimage.median_filter(fft_im, size=3)
else: fft_im_median = fft_im
#determine profile size
sx, sy = fft_im.shape
profile_size = int(sx/length_ratio)
#convert FFT image to polar coordinates
polar_im = self.topolar(fft_im_median, profile_size, profile_ang_res, order=1)
#compute radial mean and maximum profiles
radial_max=np.amax(polar_im,1)
radial_mean=np.mean(polar_im,1)
#median-filter the radial mean profile to smooth this further
radial_mean_median = scipy.signal.medfilt(radial_mean)
#radial profile is radial-max minus radial-mean
radial_profile = np.atleast_2d(radial_max-radial_mean_median)
return radial_profile
#XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
#Initialization Function
def __init__(self, img):
try:
#Create an index that is incremented each time data is processed. (This index is not utilized in this script.)
self.i = 0
self.name = "ImageO"
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 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 and calibration 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'
#XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
#Set the size and calibration of the result image
(self.scale,self.unit_string) = self.calibrate_result(x_scale,input_sizex,profile_result_length_ratio,x_unit)
r_img_size=self.profile_length(input_sizex,profile_result_length_ratio)
#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)
#Copy tags from original to processed image
Tag_Copy(self.imgref, self.result_images[self.name],subPath='OriginalImageTags')
self.result_images[self.name].ShowImage()
#get numpy array from result image
self.result_array = self.result_images[self.name].GetNumArray()
#Set the image name which will be displayed in the image window's title bar
self.result_images[self.name].SetName("Radial Max Profiles of "+img.GetName())
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()
(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])
#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)
#Copy tags from original to processed image
Tag_Copy(self.imgref, self.result_images[self.name],subPath='OriginalImageTags')
#Set the image name which will be displayed in the image window's title bar
self.result_images[self.name].SetName("Radial Max Profiles of "+self.imgref.GetName())
#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
#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)]
#compute radial FFT profile from the image, and place this profile into results image
self.result_array[:,self.i] = self.FFT_radial_profile(self.data,profile_angular_sampling_resolution,profile_result_length_ratio, median)
#XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
#Update the result image display
self.result_images[self.name].UpdateImage()
#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.