feat: added sagittalRMS slope error computation

refact:	computation of height and slope CPSDs

heightmap.py

@@ -98,10 +98,13 @@ class HeightMap(object):
         self.fitparams = None
 
         ## slope rms CPSD
-        self.srms = None
+        self.srms = None # cumulated tangential slope error
+        self.srms_sag=None #cumulated sagittal slope error
 
-        ## slope PSD frequency array
+        ## tangential slope PSD frequency array
         self.sf = None
+        ## sagittal slope PSD frequency array
+        self.sf_sag = None
         
         ## detrending method used
         self.detrend= None
@@ -428,10 +431,10 @@ class HeightMap(object):
 #end compute_y_slope
 
 
-    def compute_CPSD(self, method='scipy'):
+    def compute_CPSD_old(self,  method='scipy'):
         """! @brief Compute PSDs and CPSDs of the data and store the cumulated rms slope in member variable HeightMap.srms
         
-            @variable method can be 'scipy' meaning PSD computed with scipy.signal.periodogram, or 'NFFT' meaning that PSD will be
+            The method can be 'scipy' meaning PSD computed with scipy.signal.periodogram, or 'NFFT' meaning that PSD will be
             computed by the function processing.holey_periodogram which uses NFFT and can handle NaN containing arrays
             
             The HeightMap.compute_CPSD function  compute first `the PSD of the height errors (HeightMap.ze)
@@ -479,8 +482,35 @@ class HeightMap(object):
         s1csd = np.cumsum(s1psd) * (self.sf[1] - self.sf[0])
         self.srms = np.sqrt(s1csd)
         print("Total cumulated rms=", self.srms[self.srms.shape[0]-1]*1e6, "µrad")
-#€nd compute_CPSD
+#€nd compute_CPSD_old
 
+    def compute_CPSD(self):
+        print("tagential RMS slope error")
+        (self.sf, self.srms)=self.CRMSslope(axis=-1)
+        print("sagittal RMS slope error")
+        (self.sf_sag, self.srms_sag)=self.CRMSslope(axis=0)
+        
+
+    def CRMSslope(self, axis=-1):
+        xf = 1.0 / self.pixel[1+axis]  # sampling rate in x if axis=-1
+        #print ("xf=", xf)
+        if self.num_invalid !=0:
+            print("\n --  ROI contains invalid values  -- Using spline interpolation on rows/cols with missing values\n")
+            windata=spline_interpolate(self.ze, rowwise=(axis==-1), window=('tukey', 0.2))
+            (sf, h2psd) = scipy.signal.periodogram(windata, xf, window='boxcar', 
+                           return_onesided=True, axis=axis)
+        else:
+            (sf, h2psd) = scipy.signal.periodogram(self.ze, xf, window=('tukey', 0.2), 
+                           return_onesided=True, axis=axis)  # one sided psd with tukey (0.2) taper along axis (=-1, default)
+        h1psd = np.mean(h2psd, axis=(1+axis) )      # average over opposite axis
+        s1psd = h1psd * (2 * np.pi * sf)**2  # height to slope in fourier space
+        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])
+        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)
 
     def print_statistics(self, height = True, raw = False, percent=0):
         """! Print statistics excluding the specified percentage of outlier pixels.
@@ -713,19 +743,27 @@ class HeightMap(object):
 #end plot2dslope
 
 
-    def plot_slope_rms(self, comment='') :
+    def plot_slope_rms(self, tangential=True, scale='lin', comment='') :
         """! Display a graph of the cumulative rms slope error versus spatial frequency
             @param comment an optional comment which will be appended to the graph title
 
             The cumulative rms slope error is the square root of the slope CPSD
         """
         plt.figure(figsize=fig.figaspect(0.5))  # new plot window
+        if tangential:
             plt.plot(self.sf*1e-3, self.srms*1e6)
+            plt.title('Cumulated RMS tangential slopes {}'.format(comment))
+        else:
+            plt.plot(self.sf_sag*1e-3, self.srms_sag*1e6)
+            plt.title('Cumulated RMS sagittal slopes {}'.format(comment))
         plt.grid(visible=True, which='major', axis='both')
         plt.xlabel(r'spatial frequency $(mm^{-1})$')
         plt.ylabel('rms slope (µrad)')
+        if scale=='log':
+            plt.loglog()
+        else:
             plt.semilogx()
-        plt.title('Cumulated RMS tangential slopes {}'.format(comment))
+        
 # end plot_slope_rms
 
     def save(self, filepath=''):
@@ -800,7 +838,7 @@ class HeightMap(object):
         x=nx.NXfield(xraw, name='x', unit='meter')
         y=nx.NXfield(yraw, name='y', unit='meter')
         rawdata=nx.NXdata(nx.NXfield(self.rawZ, name='height', unit='meter'),[y,x],name='raw_data')
-        #create detrended data group wit height and slope fields
+        #create detrended data group with height and slope fields
         height=nx.NXfield(self.ze, name='height', unit='meter')
         xview=np.linspace((self.ROIorigin[0]-origin_x)*xpix,(self.ROIorigin[0]+self.ROIsize[0]-origin_x)*xpix,num=self.z.shape[1], endpoint=False)
         yview=np.linspace((self.ROIorigin[1]-origin_y)*ypix,(self.ROIorigin[1]+self.ROIsize[1]-origin_y)*ypix,num=self.z.shape[0], endpoint=False)
@@ -835,6 +873,7 @@ class HeightMap(object):
         process.data=nx.NXlink(rawdata)
         
         detrended_slopes=nx.NXdata(slope_x,[nx.NXlink(y),nx.NXlink(x)], slope_y=slope_y, name='detrended_slopes')
+        detrended_slopes.auxiliary_signals=nx.NXattr('slope_y')
         #create the process group with detrending parameters and ROI definition
         process=nx.NXprocess()
         process.program=nx.NXattr('numerical_derivation')
@@ -848,6 +887,7 @@ class HeightMap(object):
         
         # add the data group require by the NXmetrology definition 
         entry.data=nx.NXdata(nx.NXlink(height), [nx.NXlink(y), nx.NXlink(x)], slope_x=nx.NXlink(slope_x), slope_y=nx.NXlink(slope_y))
+        entry.data.auxiliary_signals=nx.NXattr(['slope_x', 'slope_y'])
         entry.data.set_default()
         #add a drawing in the sample group if defined
         if self.sample_drawing!=None:

processing.py

@@ -148,8 +148,8 @@ def fit_toroid1(xin,yin,z):
 
 
 # Replacement of periodogram with NFFT skipping NaN values
-# (self.sf, h2psd) = scipy.signal.periodogram(self.ze, xf, window=('tukey', 0.2),
-#                                      return_onesided=True)  # one sided psd with tukey (0.2) taper along last (fast varying) axis (=-1, default)
+
+
 def holey_periodogram(data, fs=1.0, window='boxcar', FFTlength=None, detrend='constant', return_onesided=True, scaling='density', axis=-1):
     
     windata=data*scipy.signal.get_window(window, data.shape[1])