fix: Invalid data in ROI no longer prevent computing PSD and CPSD

     After applying the Tukey window, Invalid data points are interpolated
     by cubic splines. If the number of bad data points in a line exceeds 20%
     this line is discarded from further computation.
     The percentage of invalid points in the ROI is printed to the console

activefigure.py

@@ -279,7 +279,7 @@ def new_active_figure(hmap, raw=False, percent=5):
     plt.connect('button_release_event', on_release)
     plt.connect('close_event', on_close)
     activefig.selectaxes=plt.gca().axes
-    print(activefig.selectaxes)
+
 
     activefig.canvas.manager.toolmanager.add_tool('Validate marked ROI', ROItool)
     activefig.canvas.manager.toolbar.add_tool('Validate marked ROI', 'Superflat', 0)

heightmap.py

@@ -93,6 +93,7 @@ class HeightMap(object):
         self.sf = 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.
         @param filename  [optional] full path to the file. If not given or en empty string, a fileopen dialog will pop-up.
@@ -286,14 +287,18 @@ class HeightMap(object):
         The returned gradient hence has the same shape as the input array.
         """
         if self.num_invalid !=0:
-            print("\nWARNING --  ROI contains invalid values  -- WARNING\n            Results might be inaccurate")
+            print("\nWARNING --  ROI contains invalid data  -- WARNING\n          some slope values might be inaccurate")
 
         self.dzdx = np.ma.masked_invalid(np.gradient(self.z, self.pixel[0], axis=-1 ) )  # derivates with respect of faster axis and mask for safety
 #end compute_x_slope
 
 
-    def compute_CPSD(self):
+    def compute_CPSD(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
+            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)
 
             The slope PSD is deduced by multiplying by the squared  frequency  and cumulated.
@@ -301,17 +306,32 @@ class HeightMap(object):
             The square root of the cumulated slope PSD is assigned to HeightMap.srms and represents the cumulated rms slope over frequencies.\n
             The frequency sampling is stored in HeightMap.sf.
         """
-        if self.num_invalid !=0:
-            print("\nWARNING --  ROI contains invalid values  -- WARNING\n            Results might be inaccurate")
 
         xf = 1.0 / self.pixel[0]  # sampling rate in x
-
+        print ("xf=", xf)
+        print("FT input type", self.ze.dtype)
         #compute from heights
+        if method=='scipy':
+            if self.num_invalid !=0:
+                # print("\nWARNING --  ROI contains invalid values  -- WARNING\n            Results might be inaccurate")
+                print("\n --  ROI contains invalid values  -- Using spline interpolation\n")
+                windata=spline_interpolate(self.ze, window=('tukey', 0.2))
+                (self.sf, h2psd) = scipy.signal.periodogram(windata, xf, window='boxcar', 
+                                               return_onesided=True)
+            else:
                 (self.sf, h2psd) = scipy.signal.periodogram(self.ze, xf, window=('tukey', 0.2), 
-                                               return_onesided=True)  # one sided psd with tukey (0.2) taper
-        h1psd = np.mean(h2psd, axis=0)       # average over l
-        s1psd = h1psd * (2 * np.pi * self.sf)**2  # height to slope in fourier space
+                                               return_onesided=True)  # one sided psd with tukey (0.2) taper along last (fast varying) axis (=-1, default)
+        elif method=='NFFT':
+            (self.sf, h2psd) = holey_periodogram(self.ze, xf, window=('tukey', 0.2),
+                                               return_onesided=True)
+        else :
+            print("invalid method. Cannot compute the PSD")
+            return
+        print("FToutput type", h2psd.dtype)
         
+        h1psd = np.mean(h2psd, axis=0)       # average over slow varying axis
+        s1psd = h1psd * (2 * np.pi * self.sf)**2  # height to slope in fourier space
+        print("    TAILLE de la PSD ", h2psd.shape, "    taille input ", self.ze.shape)
 
         # compute from slopes
         # (sf, s2psd) = scipy.signal.periodogram(self.dzdx, wf, window=('tuckey',0.2),
@@ -323,6 +343,7 @@ class HeightMap(object):
         #  s1csd = np.cumsum(s1psd[::-1])[::-1] * (self.sf[1] - self.sf[0])
         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
 
 
@@ -501,3 +522,31 @@ class HeightMap(object):
         plt.title('Cumulated RMS tangential slopes {}'.format(comment))
 # end plot_slope_rms
 
+
+    def testplot(self, data, scale=1, percent=5):
+        uu=data*scale
+        
+        bounds=np.array([percent/2, 100-percent/2])
+        #percentile do no work with masked arrays we need to force it
+        vmask=~np.ma.getmaskarray(uu)     # uu.mask
+        vbounds=np.percentile(uu[vmask],bounds)
+
+        #print( "map range: min/max = ", np.min(uu), np.max(uu), " percentile :" ,vbounds )
+
+        (w, h)=fig.figaspect(uu)
+
+        warnings.simplefilter('ignore', category=UserWarning)
+        plotfig=plt.figure(figsize=(w, h), constrained_layout=True)
+        warnings.resetwarnings()
+
+        plotfig.gca().invert_yaxis()  # images are scanned in line order, from top left to bottom right
+
+        plt.pcolormesh(self.x*1e-3, self.y*1e3, uu, cmap=plt.get_cmap('rainbow'),
+                  vmin=vbounds[0], vmax=vbounds[1], shading='auto')
+                 # vmin=np.min(uu), vmax=np.max(uu), shading='auto')
+
+        cbar= plt.colorbar(aspect=50)
+        cbar.set_label('windowed heights (nm)') #, rotation=270)
+        plt.xlabel('x (mm)')
+        plt.ylabel('y (mm)')
+        plt.title("testplot")

processing.py

@@ -14,6 +14,8 @@ __version__ = "0.1.0"
 
 
 import numpy as np
+import scipy.signal
+import nfft
 
 
 def fit2D_poly(xin,yin,z, powers,rot=0):
@@ -144,3 +146,85 @@ def fit_toroid1(xin,yin,z):
 
     return (coeffs, residues) # return  the coeffs of the solution and the residual matrix
 
+
+# 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])
+    xbase=np.linspace(-0.5, 0.5, num= data.shape[1], endpoint=False)
+    # print(xbase)
+    if FFTlength:
+        N=2*(FFTlength//2)
+    else :
+        N= 2*(data.shape[1]//2)
+    N2=N//2
+    F = np.ma.empty(dtype=np.cdouble, shape=(data.shape[0], N) )
+    for row in range(data.shape[0]):
+        segment=np.ma.array(scipy.signal.detrend(windata[row, :], type=detrend), mask=windata[row, :].mask)
+        num_invalid=segment[segment.mask].size
+        if num_invalid :
+            print (row, ":" , num_invalid, "invalid data points")
+            if num_invalid > 0.2 * data.shape[1] :
+                print( "   too many datapoints, line will be skipped from computation")
+                F[row,:]=np.ma.masked_all(shape=(N,))
+                continue
+        x=xbase[~ segment.mask]
+        y=segment[~ segment.mask]
+        # F[row,:]=np.ma.zeros(shape=(N,))
+        F[row,:]=nfft.nfft_adjoint(x,y,N, sigma=20, tol=1E-8, m=None, kernel='bspline')
+        # if num_invalid :
+            # print (x.shape,y.shape)
+
+    if scaling == 'density' :
+        if return_onesided :
+            psd=np.square(np.abs(F[:,N2:]))*2/(N*fs)
+            psd[:,0]=np.square(np.abs(F[:,N2]))/(N*fs)
+            return ( np.linspace(0, 0.5*fs, num= N2, endpoint=False), psd)
+        else :
+            return ( np.linspace((-0.5+1/N)*fs, 0.5*fs, num= N2, endpoint=True),   np.square(np.abs(F))/(N*fs))
+    elif scaling == 'spectrum' :
+        #Not shure this is what scipy.signal calls spectrum
+        if return_onesided :
+            return ( np.linspace(0, 0.5*fs, num= N2, endpoint=False),   np.square(np.abs(F[:,N2:])/N))
+        else :
+            return ( np.linspace((-0.5+1/N)*fs, 0.5*fs, num= N2, endpoint=True),   np.square(np.abs(F)/N))
+
+def spline_interpolate(data, window=('tukey', 0.2) ):
+    windata=data*scipy.signal.get_window(window, data.shape[1])
+    xbase=np.arange(0., data.shape[1],dtype=float)
+    outrow=0
+    invalids=0
+    outdata=np.ma.empty(dtype=float, shape=data.shape)
+    for inrow in range(data.shape[0]):
+        segment=np.ma.array(windata[inrow, :], mask=windata[inrow, :].mask)
+        if segment.mask[0]:
+            segment[0]=0
+            segment.mask[0]=False
+        if segment.mask[-1]:
+            segment[-1]=0
+            segment.mask[-1]=False
+                
+        num_invalid=segment[segment.mask].size
+        if num_invalid :
+            print (inrow, ":" , num_invalid, "invalid data points")
+            if num_invalid > 0.2 * data.shape[1] :
+                print( "   too many datapoints, line will be skipped from computation")
+                invalids+=data.shape[1]
+                continue
+                
+            invalids+=num_invalid
+            
+            x=xbase[~ segment.mask]
+            y=segment[~ segment.mask]
+            # defining 1st derivative null at ends (periodic is not working)
+            spline3=scipy.interpolate.make_interp_spline(x, y, bc_type=([(1, 0.0)], [(1, 0.0)]))
+            outdata[outrow,:]=spline3(xbase, extrapolate='periodic')
+        else:
+            outdata[outrow,:]=windata[inrow,:]
+        outrow+=1
+    np.ma.resize(outdata, (outrow,))
+    print ("percent of invalid data=", 100*invalids/(data.shape[0]*data.shape[1]) )
+    print("output shape", outdata.shape)
+    return outdata
\ No newline at end of file