Center Diffraction Disk

Edit Delete

Details

Name	Center Diffraction Disk
Permalink	https://serialemscripts.nexperion.net/script/86
Category	Calibrations
Language	Python
Description	A Python script to center the Diffraction Disk on four segments of the Panther detector
Author	Peter Kirchweger
Literature reference	https://doi.org/10.1101/2025.11.26.690630

Version History

Version 1.0 (December 08, 2025; 96 download(s); SerialEM version ≥ 4.2)

Version 1.0

Download

Upload date	December 08, 2025
Download count	96
SerialEM version ≥	4.2

#!Python
### NEW!!!!!
#=============================
#ScriptName CenterDisk
#  Author: Peter Kirchweger & Shahar Seifer
# (C)2025 Elbaum lab, Weizmann Institute of Science
#=============================

from autoscript_tem_microscope_client import TemMicroscopeClient
from autoscript_tem_microscope_client.enumerations import *
from autoscript_tem_microscope_client.structures import *
import numpy as np
import math
import time
from matplotlib import pyplot as plt
import pandas as pd
from sklearn import linear_model

import serialem as sem

factor_sensitivity=10
attenuation=0.3

def solve_theta(C: float):
# Solve theta-sin(theta)=C by Newton's method
Cnorm = C
if Cnorm > 2 * math.pi:
   Cnorm = 2 * math.pi
   if Cnorm < -2 * math.pi:
      Cnorm = -2 * math.pi
      theta = Cnorm
      for n in range(0, 10):
      if abs(-1 + math.cos(theta)) > 0.001:
         theta = theta - (Cnorm - theta + math.sin(theta)) / (-1 + math.cos(theta));
         return theta
         
         microscope = TemMicroscopeClient()
         microscope.connect("192.168.0.1", 7521)
         time.sleep(1)
         
         bf_s = microscope.detectors.get_scanning_detector(DetectorType.BF_S)
         df_s = microscope.detectors.get_scanning_detector(DetectorType.DF_S)
         
         # Enable the inner segments on the BF-S detector
         bf_s.set_enabled_segments(DetectorSegmentType.INNER_RING)
         
         # Enable all segments on the DF-S detector
         df_s.set_enabled_segments(DetectorSegmentType.ALL)
         
         # Join the two lists together
         segments = bf_s.get_enabled_segments() + df_s.get_enabled_segments()
         RR = 40.0
         df_x = pd.DataFrame(columns=['con_x', 'con_y', 'x'], index=[0, 1, 2])
         df_y = pd.DataFrame(columns=['con_x', 'con_y', 'y'], index=[0, 1, 2])
         arbstep = 0.05
         
         # Convergence parameters
         convergence_threshold = 0.10  # 10% difference threshold
         large_difference_threshold = 1.00  # 100% difference threshold
         large_step_multiplier = 3.0  # Make steps 3x larger when far from center
         
         # Track best position throughout entire run
         best_difference = float('inf')
         best_con_x = 0
         best_con_y = 0
         best_iteration = -1
         
         for iteration in range(0, 18):
         
         # Acquire one image per segment, a total of 12 images in this case
         images = microscope.acquisition.acquire_stem_segment_images(segments, 256, 3E-6)
         
         # Plot the images in max 4 columns wide
         n = len(images)
         columns = min(4, n)
         rows = math.ceil(n / columns)
         fig, plots = plt.subplots(rows, columns)
         fig.tight_layout()
         # Create a new dictionary for the images and the mean
         image_mean_dict = {}
         mean_intensity=np.zeros(12)
         
         for i in range(0, 12):
         gray = np.mean(images[i].data, axis=0)
         mean_intensity[i] = np.mean(gray)
         title = segments[i].detector_name + "-" + segments[i].name
         image_mean_dict[title] = round(mean_intensity[i], 2)
         column = i % 4
         row = math.floor(i / 4)
         plot = plots[row, column]
         plot.axis('off')
         plot.imshow(images[i].data, cmap='gray')
         plot.set_title(title)
         plot.text(100, 0, 'Mean: ' + str(round(mean_intensity[i], 2)), fontsize=8, ha='center')
         sem.Echo()
         print("################################################")
         print("Starting round " + str(iteration))    
         print(image_mean_dict)                 
         
         # Read dif.alignement from SerialEM
         (con_x, con_y) = sem.ReportDiffractionShift()
         con_x=con_x*factor_sensitivity
         con_y=con_y*factor_sensitivity
         
         # Get channel values directly without bias subtraction
         ch1 = image_mean_dict.get('DF-S-Outer1')
         ch2 = image_mean_dict.get('DF-S-Outer2')
         ch3 = image_mean_dict.get('DF-S-Outer3')
         ch4 = image_mean_dict.get('DF-S-Outer4')
         
         print(round(ch1,3), round(ch2,3), round(ch3,3), round(ch4,3))
         
         # Check convergence and track best position
         ch_mean = (ch1 + ch2 + ch3 + ch4) / 4.0
         ch_max = max(ch1, ch2, ch3, ch4)
         ch_min = min(ch1, ch2, ch3, ch4)
         
         if ch_mean > 0:
            relative_difference = (ch_max - ch_min) / ch_mean
            print(f'Relative difference: {relative_difference:.3f} ({relative_difference*100:.1f}%)')
            
            # Track best position at ANY iteration
            if relative_difference < best_difference:
               best_difference = relative_difference
               best_con_x = con_x
               best_con_y = con_y
               best_iteration = iteration
               print(f'NEW BEST POSITION! con_x={con_x:.3f}, con_y={con_y:.3f}, diff={relative_difference:.3f}, iteration={iteration}')
               
               # Check if converged
               if relative_difference < convergence_threshold:
                  print(f'CONVERGED! Channel difference < {convergence_threshold*100}%')
                  print(f'Final channel values: ch1={ch1:.2f}, ch2={ch2:.2f}, ch3={ch3:.2f}, ch4={ch4:.2f}')
                  break
                  
                  # Adaptive step size based on difference
                  if relative_difference > large_difference_threshold:
                     current_attenuation = attenuation * large_step_multiplier
                     print(f'Large difference detected - using {large_step_multiplier}x step size')
                   else:
                     current_attenuation = attenuation
                   else:
                     current_attenuation = attenuation
                     
                     sumright = ch1 + ch2
                     sumleft = ch3 + ch4
                     sumup = ch1 + ch4
                     sumdown = ch2 + ch3
                     theta1 = solve_theta(math.pi * (1 + abs(sumdown - sumup) / (abs(sumdown) + abs(sumup) + 100)))
                     theta2 = solve_theta(math.pi * (1 + abs(sumright - sumleft) / (abs(sumright) + abs(sumleft) + 100)))
                     y = RR * (math.sqrt(1.0 - pow(math.sin(theta2 / 2), 2)))
                     x = RR * (math.sqrt(1.0 - pow(math.sin(theta1 / 2), 2)))
                     # mind the signs
                     if (sumup < sumdown):
                        x = -x
                        if (sumright < sumleft):
                           y = -y
                           print('x= {0:.6}   y= {1:.6}'.format(x,y))
                           df_x.loc[iteration] = [con_x, con_y, x]
                           df_y.loc[iteration] = [con_x, con_y, y]
                           
                           regx_ready = True
                           regy_ready = True
                           if iteration < 8:
                              regx_ready = False
                              regy_ready = False
                              
                              if regx_ready:
                                 regrx = linear_model.LinearRegression()
                                 regrx.fit(df_x[['con_x', 'con_y']], df_x['x'])
                                 print('regrx coef= {0:.6} , {1:.6}'.format(regrx.coef_[0],regrx.coef_[1]))
                                 if (abs(regrx.coef_[0]) < 0.05 and abs(regrx.coef_[1]) < 0.05):
                                    regx_ready = False
                                    if regy_ready:
                                       regry = linear_model.LinearRegression()
                                       regry.fit(df_y[['con_x', 'con_y']], df_y['y'])
                                       print('regry_coef= {0:.6} , {1:.6}'.format(regry.coef_[0],regry.coef_[1]))
                                       if (abs(regry.coef_[0]) < 0.05 and abs(regry.coef_[1]) < 0.05):
                                          regy_ready = False
                                          
                                          if not (regx_ready) and iteration % 2 == 0:
                                             con_x = con_x +  arbstep
                                             elif not (regy_ready) and iteration % 2 == 1:
                                             con_y = con_y +  arbstep
                                             if not iteration==3:
                                                arbstep = -arbstep
                                                elif regx_ready and iteration % 2 == 0:
                                                if abs(regrx.coef_[0]) > 0.2 * abs(regrx.coef_[1]):
                                                   con_x = con_x - current_attenuation* x / regrx.coef_[0]
                                                   if abs(regrx.coef_[1]) > 0.2 * abs(regrx.coef_[0]):
                                                      con_y = con_y - current_attenuation* x / regrx.coef_[1]
                                                      elif regy_ready and iteration % 2 == 1:
                                                      if abs(regry.coef_[0]) > 0.2 * abs(regry.coef_[1]):
                                                         con_x = con_x - current_attenuation*y / regry.coef_[0]
                                                         if abs(regry.coef_[1]) > 0.2 * abs(regry.coef_[0]):
                                                            con_y = con_y - current_attenuation*y / regry.coef_[1]
                                                            
                                                            #limit shift
                                                            con_x=max(min(con_x, 2), -2)
                                                            con_y=max(min(con_y, 2), -2)
                                                            
                                                            # Modify dif.alignement by SerialEM
                                                            message = "sem.SetDiffractionShift({0:.6f},{1:.6f})".format(con_x/factor_sensitivity, con_y/factor_sensitivity)
                                                            exec(message)
                                                            print(con_x/factor_sensitivity, con_y/factor_sensitivity)
                                                            
                                                            # After all iterations, return to best position found
                                                            print('################################################')
                                                            print('################################################')
                                                            print(f'Alignment complete. Best position was at iteration {best_iteration}')
                                                            print(f'Best difference: {best_difference:.3f} ({best_difference*100:.1f}%)')
                                                            print(f'Moving to best position: con_x={best_con_x:.3f}, con_y={best_con_y:.3f}')
                                                            message = "sem.SetDiffractionShift({0:.6f},{1:.6f})".format(best_con_x/factor_sensitivity, best_con_y/factor_sensitivity)
                                                            exec(message)
                                                            print('################################################')
                                                            print('################################################')
                                                            print(f'Final position set to: {best_con_x/factor_sensitivity:.6f}, {best_con_y/factor_sensitivity:.6f}')
                                                            print('################################################')
                                                            print('################################################')

Related Scripts

The following scripts provide a similar or related functionality, but are not required to run this script:

License

This code can be shared and adapted under the Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International (CC BY-NC-SA 4.0) License (https://creativecommons.org/licenses/by-nc-sa/4.0/)