feat: added computation and plot of cumulated height rms

.gitignore

@@ -6,6 +6,8 @@ unused/
 #doxygen 
 doxygen/*
 docs/
+dabam/
+dabam.*
 
 #CB files 
 *.cbp

edit_results.py

@@ -51,6 +51,32 @@ def plot_slope_rms(maplist, ymax=1., comment='' ) :
     plt.title('Cumulated RMS tangential slopes {}'.format(comment))
     plt.legend()
 
+
+def plot_height_rms(maplist, ymax=0.5, comment='' ) :
+    """! Display a graph of the cumulative rms height error versus spatial frequency
+        @param comment an optional comment which will be appended to the graph title
+
+        The cumulative rms Height error is the square root of the Height CPSD
+    """
+    graph=plt.figure(figsize=fig.figaspect(0.5))  # new plot window
+    for hmap in maplist :
+        if type(hmap.hrms) == type(None):
+            return
+        (label,color,linewidth, linestyle)=hmap.get_curveparms()  
+        plt.plot(hmap.sf*1e-3, (hmap.hrms[-1]-hmap.hrms)*1e9, label=label,color=color, lw=linewidth, ls=linestyle)
+    plt.grid(visible=True, which='major', axis='both')
+    # plt.axis([.01,10., 0, ymax])
+    plt.semilogx()
+    # plt.loglog()    
+    graph.axes[0].set_ylim([0,ymax])
+    plt.xlabel(r'spatial frequency $(mm^{-1})$')
+    plt.ylabel('rms height (nm)')
+
+    plt.title('Cumulated RMS tangential heights {}'.format(comment))
+    plt.legend()
+
+
+
 def slope_map(maplist,direction='x', filter=None, vbounds=None, color='rainbow'):
     ''' Display all slop maps in black and white with identical scale
         @param direction the derivation direction 'x'  or  'y'

heightmap.py

@@ -122,6 +122,9 @@ class HeightMap(object):
         #added 15/11/2023
         self.tanCSrms_line=None
         self.sagCSrms_line=None
+        #added 20/03/2024
+        self.hrms=None
+        self.hrms_sag=None
 
     def read_Zygofile(self, filename='', origin='zero' ):
         """! Reads a Zygo datafile in .datx, .dat or xyz format, and loads it into the HeightMap object.
@@ -494,9 +497,9 @@ class HeightMap(object):
 
     def compute_CPSD(self):
         print("tagential RMS slope error")
-        (self.sf, self.srms)=self.CRMSslope(axis=-1)
+        (self.sf, self.srms, self.hrms)=self.CRMSslope(axis=-1)
         print("sagittal RMS slope error")
-        (self.sf_sag, self.srms_sag)=self.CRMSslope(axis=0)
+        (self.sf_sag, self.srms_sag, self.hrms_sag)=self.CRMSslope(axis=0)
         
 
     def CRMSslope(self, axis=-1):
@@ -515,10 +518,14 @@ class HeightMap(object):
         print("    Size of the PSD ", h2psd.shape, "   size of input ", self.ze.shape)
         # uncomment to reverse cumulation direction
         #  s1csd = np.cumsum(s1psd[::-1])[::-1] * (sf[1] - sf[0])
+        # ajout 20/03/24
+        hcpsd=np.cumsum(h1psd) * (sf[1] - sf[0])
+        hrms=np.sqrt(hcpsd)
+        #
         s1csd = np.cumsum(s1psd) * (sf[1] - sf[0])
         srms = np.sqrt(s1csd)
         print("Total cumulated rms=", srms[srms.shape[0]-1]*1e6, "µrad")
-        return (sf,srms)
+        return (sf,srms, hrms)
 
     def print_statistics(self, height = True, raw = False, percent=0):
         """! Print statistics excluding the specified percentage of outlier pixels.
@@ -576,7 +583,13 @@ class HeightMap(object):
 
 # end print_statistics
     def clip(self, origin, size):
+        '''
+            clip raw data to avoid full rows and columns of invalid data
+        '''
         self.rawZ=self.rawZ[origin[1]:origin[1]+size[1], origin[0]:origin[0]+size[0]]
+        # we need also to clip self.rawX and self.rawY
+        self.rawX=self.rawX[origin[0]:origin[0]+size[0]]
+        self.rawY=self.rawY[origin[1]:origin[1]+size[1]]
 
     def set_ROI(self, firstX, firstY, nwidth, nheight, flip=None):
         """! define a rectangular ROI and set the `HeightMap.z`array accordingly
@@ -783,6 +796,31 @@ class HeightMap(object):
             plt.semilogx()  
 # end plot_slope_rms
 
+    def plot_height_rms(self, tangential=True, scale='log', comment='') :
+        """! Display a graph of the cumulative rms height error versus spatial frequency
+            @param comment an optional comment which will be appended to the graph title
+
+            The cumulative rms height error is the square root of the height CPSD
+        """
+        from crmsfig import create_crms_figure
+        csfig=create_crms_figure(figsize=fig.figaspect(0.5))  # new plot window
+
+        if tangential:
+            plt.plot(self.sf*1e-3, self.hrms*1e9)
+            plt.title('Cumulated RMS tangential heights {}'.format(comment))
+        else:
+            plt.plot(self.sf_sag*1e-3, self.hrms_sag*1e9)
+            plt.title('Cumulated RMS sagittal heights {}'.format(comment))
+        plt.grid(visible=True, which='major', axis='both')
+        plt.xlabel(r'spatial frequency $(mm^{-1})$')
+        plt.ylabel('rms height (nm)')
+        if scale=='log':
+            plt.loglog()
+            csfig.canvas.mpl_connect('close_event', self.on_close_fig)
+        else:
+            plt.semilogx()  
+# end plot_height_rms
+
     def on_close_fig(self,event):
         from crmsfig import Crms_fig
         if isinstance(event.canvas.figure, Crms_fig):

processing.py

@@ -42,15 +42,15 @@ def fit2D_poly(xin,yin,z, powers,rot=0):
 
     # grid coords
     x, y = np.meshgrid(xin, yin)
-    if(rot !=0 ):
     # print("fit2D_poly called with rotation angle:", rot)
+    if(rot !=0 ):
         R=np.array([[np.cos(rot), np.sin(rot)],[-np.sin(rot),np.cos(rot)]])
         V=np.array([x.ravel(),y.ravel()])
         #print(R.shape, V.shape)
         U=R@V
         x=np.reshape(U[0,:], x.shape)
         y=np.reshape(U[1,:], y.shape)
-    #print( "shape X", x.shape, "shape Y", y.shape, "shape Z" , z.shape)
+    print( "shape X", x.shape, "shape Y", y.shape, "shape Z" , z.shape)
     x=np.ma.array(x,mask=z.mask)
     y=np.ma.array(y,mask=z.mask)
     #vmask=~z.mask  Not valid if z has no invalid data