Fix: General update and bug fix

	Specialized ZygoDatx added. Datx mechanism to find Surface date
	improved for robustness
	Graphic ROI selection display coordinate when dragging the rectangle

ZygoDat.py

@@ -43,7 +43,7 @@ class ZygoSoftwareVersion(BigEndianStructure):
 
 class ZygoSoftwareInfo(BigEndianStructure):
     _pack_ = 2
-    _fields_ = [("Type", c_short),
+    _fields_ = [("Type", c_short),      # The softwares are: unknown or MX (0), MetroPro (1), MetroBASIC (2), and d2bug (3).
                 ("Date", c_char * 30),
                 ("Version", ZygoSoftwareVersion) ]
 
@@ -282,7 +282,7 @@ class DatFile(object) :
         self.filepath=""
         """ (`string`) The full path to access the file"""
         self.header=Header()
-        """ (`Header`) File header conataining all metadata relative to a measurement"""
+        """ (`Header`) File header containing all metadata relative to a measurement"""
 
 
 
@@ -316,7 +316,7 @@ class DatFile(object) :
             return False
         return True
 
-    def getPhase(self, multiplier=1.):
+    def getHeight(self, multiplier=1.):
         """ Reads the phase data from thefile and returns it in a `numpy.array` \n
             Possible invalid data are replaced by NaN in the returned array\n
             Parameters:
@@ -362,9 +362,10 @@ class DatFile(object) :
         the pixel size in  \\(m^{-1}\\) (as a float)
 
         """
+
         return self.header.BEattributes.LateralResolution
 
-    def get_origin(self):
+    def getOrigin(self):
         """
         Return
         ------

ZygoDatx.py

+# -*- coding: utf-8 -*-
+
+## @file ZygoDatx.py
+# This file provides tools for analysing surface height data with statistical and PSD tools.
+#
+# The file defines the HeightMap class, which is the main storage element of surface height measurements.
+#   @author François Polack, Uwe Flechsig
+#   @copyright 2022, Superflat-PCP project  <mailto:superflat-pcp@synchrotron-soleil.fr>
+#   @version 0.1.0
+
+
+import tkinter
+from tkinter import filedialog
+from pathlib import Path
+import h5py as h5
+import numpy as np
+
+class Node:
+
+
+    def __init__(self, name, value=None):
+        self.name=name
+        self.value=value
+        self.parent=None
+        self.children=[]
+
+    def append_node(self,node):
+        if not isinstance(node, Node):
+            raise TypeError("Argument of append_node mus be an instance of class Node")
+        self.children.append(node)
+        node.parent=self
+        return node
+
+    def __getitem__(self, key):
+        for node in self.children:
+            if node.name==key:
+                return node
+        raise KeyError
+
+    def __iter__(self):
+        return self.children.__iter__()
+
+    def __len__(self):
+        return len(self.children)
+
+    def __setitem__(self, key, value):
+        for nnd in self.children:
+            if nnd.name==key:
+                self.children.remove(nnd)
+        if not isinstance(key,type("")):
+            raise TypeError("Key must be of 'str' type")
+        node=Node(key,value)
+        self.children.append(node)
+        node.parent=self
+
+    def dump(self, level=0):
+        indent="\t" * level
+        for node in self.children:
+            print(indent, node.name,':' , node.value)
+            if len(node)!=0:
+                node.dump(level+1)
+
+class DatxFile(object):
+    def __init__(self):
+        ## (str) The full path to access the file
+        self.filepath=""
+        ## (Node) The tree of metadata constructed from the 'MetaData' group of the datx file
+        self.metatree=Node('Root')
+        ## (int) the layout version of the file (0 if not defined)
+        self.fileversion=0
+
+        self.attributes=None
+
+
+    def openFile(self, filename=''):
+
+        if Path(filename).is_file():
+            self.filepath=filename
+        else:
+            tkinter.Tk().withdraw()
+            self.filepath=filedialog.askopenfilename(title="open a Zygo.dat file", filetypes=(("Zygo datx","*.datx")) )
+        print("Opening:", self.filepath)
+
+
+        h5file=h5.File(self.filepath, 'r')
+
+        if h5file['Attributes'].attrs.__contains__('File Layout Version'):
+            self.fileversion=h5file['Attributes'].attrs['File Layout Version'][0]
+        else:
+            self.fileversion=0
+
+        print("File version=", self.fileversion)
+
+        nodedict=dict(Root=self.metatree)
+
+        metadata=h5file["/MetaData"] # open the MetaData DS to get the link to the Surface DS
+
+        for line in metadata:
+            src=nodedict[line['Source'].decode('ascii')]
+            dest=line['Destination'].decode('ascii')
+            link= line['Link'].decode('ascii')
+            src[link]=dest
+            if dest[0]=='{':  # this is a group node we need to add it to the node list
+                nodedict[dest]=src[link]
+
+        self.attributes=h5file[self.metatree['Measurement']['Attributes'].value ].attrs
+        self.info()
+        return True
+
+    def getHeight(self, multiplier=1.):
+        h5file=h5.File(self.filepath, 'r')
+        dsSurf=h5file[self.metatree['Measurement']['Surface']['Path'].value]
+        scale_factor=multiplier*dsSurf.attrs['Z Converter'][0][2][1]*dsSurf.attrs['Z Converter'][0][2][2]*dsSurf.attrs['Z Converter'][0][2][3]
+
+
+       # print("shape="+str(dsSurf.shape)+ " type="+ str(dsSurf.dtype)+ " compression="+ dsSurf.compression+ " chunks="+ str(dsSurf.chunks))
+
+        try:
+            Z=dsSurf[()]
+        except OSError as e:
+            if dsSurf.compression == 'szip' :
+                msg="The Surface dataset of this file was compressed with szip and cannot be read by this version of h5py.\n"
+                msg+="You may consider to update the hdf5.dll and hdf5_hl.dll library files located in the Lib/site-packages/h5py directory of your Conda environment"
+            else:
+                msg="The Surface dataset cannot be read by this version of h5py.\n"
+                msg+="shape="+str(dsSurf.shape)+ " type="+ str(dsSurf.dtype)+ " compression="+ dsSurf.compression+ " chunks="+ str(dsSurf.chunks)
+            raise OSError(msg)
+        Z[np.greater_equal(Z, dsSurf.attrs['No Data'][0]) ]=np.nan  #  invalid data set to NaN
+        return scale_factor*Z
+
+
+    def getPixel(self):
+        h5file=h5.File(self.filepath, 'r')
+        dsSurf=h5file[self.metatree['Measurement']['Surface']['Path'].value]
+        return [dsSurf.attrs['X Converter'][0][2][1], dsSurf.attrs['Y Converter'][0][2][1] ]
+
+
+    def getOrigin(self):
+        h5file=h5.File(self.filepath, 'r')
+        dsSurf=h5file[self.metatree['Measurement']['Surface']['Path'].value]
+        return [dsSurf.attrs['Coordinates'][0][0],dsSurf.attrs['Coordinates'][0][1]]
+
+
+    def info(self):
+
+        softlist=('Mx','MetroPro', 'MetroBASIC', 'd2bug')
+        print("Content of file", self.filepath)
+        self.metatree.dump(1)
+        if self.attributes.__contains__('Data Context.Data Attributes.Software Info Name'):
+            print("Acquisition software", self.attributes['Data Context.Data Attributes.Software Info Name'][0], "version",
+                   self.attributes['Data Context.Data Attributes.Software Info Version'][0] )
+
+

ZygoXYZ.py

@@ -255,7 +255,7 @@ class XYZfile(object):
 
 
 
-    def getPhase(self, multiplier=1.):
+    def getHeight(self, multiplier=1.):
         """ Reads the phase data from thefile and returns it in a `numpy.array` \n
             Possible invalid data are replaced by NaN in the returned array\n
             Parameters:
@@ -312,7 +312,7 @@ class XYZfile(object):
         """
         return self.PixelSize
 
-    def get_origin(self):
+    def getOrigin(self):
         """
         Return
         ------
@@ -320,4 +320,3 @@ class XYZfile(object):
             the oordinates (in pixels) of the upper left corner of the phase array in the interferometer camera frame
         """
         return self.PhaseRect[0:2]
-   
\ No newline at end of file

activefigure.py

@@ -18,7 +18,7 @@ import warnings
 import matplotlib.pyplot as plt
 import matplotlib.figure as fig
 import matplotlib.patches as patches
-from matplotlib.backend_tools import ToolBase   #, ToolToggleBase
+from matplotlib.backend_tools import ToolBase   #, ToolCursorPosition   #, ToolToggleBase
 from processing import *
 
 
@@ -76,7 +76,9 @@ class ActiveFigure(fig.Figure):
         self.r1=None
         """ Region defining white rectangle """
         self.r2=None
-        """ Black dashe rectangle superimposed to `ActiveFigure.r1` """
+        """ Black dashed rectangle superimposed to `ActiveFigure.r1` """
+        self.pxorg=None
+
         self.ROI=[0,0,0,0]
         """ ROI rectangle set-up by the mouse up event """
         self.selectaxes=None
@@ -94,6 +96,14 @@ class ActiveFigure(fig.Figure):
                 org=(event.xdata,event.ydata)
                 self.r1.set_xy(org)
                 self.r2.set_xy(org)
+
+                self.pxorg=(int(round(1e-3*org[0]/self.hmap.pixel[0],0)-self.hmap.origin[0]),
+                            int(round(1e-3*org[1]/self.hmap.pixel[1],0)-self.hmap.origin[1]) )
+
+                smsg="origin:({:d},{:d}) size:({:d}, {:d})  [{:.1f}, {:.1f} : {:.1f}, {:.1f}]mm".format(
+                self.pxorg[0],self.pxorg[1], 0, 0,  org[0],org[1],0, 0)
+                event.canvas.manager.toolbar.set_message(smsg)
+
                 self.mvcid = event.canvas.mpl_connect('motion_notify_event', on_mouse_move)
 
     def on_mouse_move(self,event):
@@ -107,6 +117,15 @@ class ActiveFigure(fig.Figure):
             self.r2.set_width(width)
             self.r2.set_height(height)
 
+
+            pxwidth=int(round(1e-3*width/self.hmap.pixel[0],0))
+            pxheight=int(round(1e-3*height/self.hmap.pixel[1],0))
+
+
+            smsg="origin:({:d},{:d}) size:({:d}, {:d})  [{:.1f}, {:.1f} : {:.1f}, {:.1f}]mm".format(
+            self.pxorg[0],self.pxorg[1], pxwidth, pxheight,  org[0],org[1],width, height)
+            event.canvas.manager.toolbar.set_message(smsg)
+
             event.canvas.draw()
 
     def on_release(self, event):
@@ -155,10 +174,14 @@ def new_active_figure(hmap, raw=False, percent=5):
     if raw:
         z=hmap.rawZ
     else:
-        (fitparams, z) = fit2D_order2(x,y,hmap.rawZ)
-
+        #  (fitparams, z) = fit2D_order2(x,y,hmap.rawZ)
+        ((Z0, pitch, roll, x2, xy, y2), z) = fit2D_order2(x,y,hmap.rawZ)
+        """
+        print("Raw fit coeffs: Z0 = {:.3e}m, pitch = {:.4e}rad,  roll = {:.4e}rad, x2=  {:.3e}m-1, xy=  {:.3e}m-1, y2=  {:.3e}m-1".
+            format(Z0, pitch, roll, x2, xy, y2))
+        """
     scale= 'mm' if raw  else  'nm'
-    title= 'Height' if raw else "Height error"
+    title= 'Height' if raw else 'Height error'
 
     if(scale=="nm"):
         zscale=1.e9
@@ -177,7 +200,12 @@ def new_active_figure(hmap, raw=False, percent=5):
 
     bounds=np.array([percent/2, 100-percent/2])
 
-    vbounds=np.percentile(uu,bounds)
+    #  vbounds=np.percentile(uu,bounds)
+    #percentile do no work with masked arrays we need to force it
+
+    vmask=~np.ma.getmaskarray(uu)
+    vbounds=np.percentile(uu[vmask],bounds)
+
     #print( "map range: min/max = ", np.min(uu), np.max(uu), " percentile :" ,vbounds )
 
     (w, h)=fig.figaspect(z)
@@ -223,6 +251,9 @@ def new_active_figure(hmap, raw=False, percent=5):
 
     activefig.canvas.manager.toolmanager.add_tool('Validate marked ROI', ROItool)
     activefig.canvas.manager.toolbar.add_tool('Validate marked ROI', 'Superflat', 0)
+   # activefig.canvas.manager.toolmanager.add_tool('position', ToolCursorPosition)
+   # activefig.canvas.manager.toolbar.add_tool('position', 'Superflat', -1)
+   # activefig.canvas.manager.toolbar.set_message('my message')
 
     activefig.hmap=hmap
 

heightmap.py

@@ -37,6 +37,7 @@ import matplotlib.pyplot as plt
 import matplotlib.figure as fig
 import ZygoDat
 import ZygoXYZ
+import ZygoDatx
 from processing import *
 from activefigure import *
 
@@ -115,56 +116,48 @@ class HeightMap(object):
             filename=filedialog.askopenfilename(title="Open a Zygo data file", filetypes=(("Zygo dat","*.dat"), ("Zygo datx","*.datx"),
                 ("Zygo XYZ", "*.xyz"), ("Zygo ascii","*.ascii") ))
 
-        if Path(filename).suffix==".datx":
-            h5file=h5.File(filename, 'r')
-            # directly open the Surface dataset
-            # dsSurf=h5file["/Measurement/Surface"] Unfortunately we can only do that on the most recent datx files,
-            # the older ones do not have the Measurement group with links
-            metadata=h5file["/MetaData"] # open the MetaData DS to get the link to the Surface DS
-
-            for i in range(0,len(metadata)):
-                if metadata[i]['Destination'].decode('ascii') == 'Surface':
-                   break
-            if i==len(metadata):
-                print("The Surface group location was not found in file", filename)
-                return False
-            print(metadata[i]['Destination'],": ", metadata[i-1]['Destination']) # .decode('ascii'))
-
-
-            dsSurf=h5file[metadata[i-1]['Destination'].decode('ascii')]
-            #    scale_factor=dsSurf.attrs['Interferometric Scale Factor'][0]*dsSurf.attrs['Obliquity Factor'][0]*dsSurf.attrs['Wavelength'][0]
-            scale_factor= dsSurf.attrs['Z Converter'][0][2][1]*dsSurf.attrs['Z Converter'][0][2][2]*dsSurf.attrs['Z Converter'][0][2][3]
 
-            # self.z=np.ma.masked_greater_equal(dsSurf, dsSurf.attrs['No Data'][0]).set_fill_value(np.nan)   # this masks the invalid data but not really changing their value to nan
-            print("\nfile", filename)
-            # print(dsSurf)
-            print("shape=", dsSurf.shape, "type", dsSurf.dtype, "compression", dsSurf.compression, "chunks", dsSurf.chunks)
-            print("no data=",dsSurf.attrs['No Data'][0])
 
-            self.z=dsSurf[()]
-            self.z[np.greater_equal(self.z, dsSurf.attrs['No Data'][0]) ]=np.nan  #  invalid data set to NaN
-            self.rawZ=self.z=scale_factor*np.ma.masked_invalid(self.z)  # and then masked
-            # set to scale
-            self.rawZ=self.z=scale_factor*self.z
-
-            self.pixel=[dsSurf.attrs['X Converter'][0][2][1], dsSurf.attrs['Y Converter'][0][2][1] ]
-            fileorigin=[dsSurf.attrs['Coordinates'][0][0],dsSurf.attrs['Coordinates'][0][1] ]
-        else:
         if Path(filename).suffix==".dat":
             datfile=ZygoDat.DatFile()
         elif Path(filename).suffix==".xyz" or Path(filename).suffix==".ascii":
             datfile=ZygoXYZ.XYZfile()
+        elif Path(filename).suffix==".datx":
+            datfile=ZygoDatx.DatxFile()
         else:
             print("Unknown file type:", Path(filename).suffix)
             return False
         if datfile.openFile(filename):      #if false invalid file message is printed
-                self.rawZ=self.z=np.ma.masked_invalid( datfile.getPhase())
+            self.rawZ=self.z=np.ma.masked_invalid( datfile.getHeight())
             if self.z.size==0:
                 return False
+
+            fileorigin=datfile.getOrigin()
+
+            pixel=datfile.getPixel() #datx return a tuple , not dat and xyz
+            if type(pixel)==list :
+                self.pixel=pixel
+            else:
                 self.pixel[1]=self.pixel[0]=datfile.getPixel()
-                fileorigin=datfile.get_origin()
+
+            if self.pixel[0]==0 or self.pixel[1]==0 : # this may happen in  processed data
+                print("The pixel size is not documented in this file.")
+                if self.pixel[0]==0:
+                    instr=input("Enter X size in mm or empty input to cancel loading >")
+                    if instr=="":
+                        return False
+                    else:
+                        self.pixel[0]=float(instr)*1e-3
+                if self.pixel[1]==0:
+                    instr=input("Enter Y size in mm or empty input to cancel loading >")
+                    if instr=="":
+                        return False
+                    else:
+                        self.pixel[1]=float(instr)*1e-3
+
         else:
             return False
+
         self.ze=None
         print("pixel size=",1e3*self.pixel[0], 'x', 1e3*self.pixel[1], " mm^2")
         print("file origin of phase image: ", self.origin)
@@ -224,6 +217,10 @@ class HeightMap(object):
             It prints the fitting parameters on standard output and return a tuple containing the parameter list and  the fit residuals.
                    """
         ((Z0, pitch, roll, x2, xy, y2), residuals) = fit2D_order2(self.x,self.y,self.z)
+        """
+        print("Raw fit coeffs: Z0 = {:.3e}m, pitch = {:.4e}rad,  roll = {:.4e}rad, x2=  {:.3e}m-1, xy=  {:.3e}m-1, y2=  {:.3e}m-1".
+            format(Z0, pitch, roll, x2, xy, y2))
+        """
         d=x2-y2
         angle=math.atan(xy/d)
         d/=math.cos(angle)
@@ -239,7 +236,7 @@ class HeightMap(object):
             ry = 0.5/ Cy
         except ZeroDivisionError:
             ry=math.inf
-        print("debug: fit params:" , Z0, roll, pitch, x2, xy, y2)
+        # print("debug: fit params:" , Z0, roll, pitch, x2, xy, y2)
         print("angle = {:.4e}rad, rx = {:.3e}m, pitch = {:.4e}rad, ry = {:.3e}m, roll = {:.4e}rad, Z0 = {:.3e}m".
               format(angle, rx, pitch, ry, roll, Z0))
         #verification  antidiagonal should be almost 0
@@ -363,7 +360,7 @@ class HeightMap(object):
             bounds=np.array([percent/2, 100-percent/2])
             print("statistics on", 100-percent, "% of ROI surface")
             #percentile do no work with masked we need to force it
-            vmask=~tu.mask
+            vmask=~np.ma.getmaskarray(tu)
             vbounds=np.percentile(tu[vmask],bounds)
             uu=np.ma.masked_outside(tu, vbounds[0],vbounds[1], copy=True)
             print("{:s}: rms= {:.3f} {:s}, pv= {:.3f} {:s}".
@@ -427,7 +424,7 @@ class HeightMap(object):
 
         bounds=np.array([percent/2, 100-percent/2])
         #percentile do no work with masked arrays we need to force it
-        vmask=~uu.mask
+        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 )
@@ -499,7 +496,7 @@ class HeightMap(object):
         plt.plot(self.sf*1e-3, self.srms*1e6)
         plt.grid(visible=True, which='major', axis='both')
         plt.xlabel(r'spatial frequency $(mm^{-1})$')
-        plt.ylabel('rms (µrad)')
+        plt.ylabel('rms slope (µrad)')
         plt.semilogx()
         plt.title('Cumulated RMS tangential slopes {}'.format(comment))
 # end plot_slope_rms

processing.py

@@ -51,8 +51,10 @@ def fit2D_poly(xin,yin,z, powers,rot=0):
     #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
+    #vmask=~z.mask  Not valid if z has no invalid data
+    vmask=~np.ma.getmaskarray(z)
     valid=z[vmask].size
+    print("Mask shape", z.mask.shape ,"Vmask shape",vmask.shape, "unmasked shape",z[vmask].shape)
 
     print ("valid data " , valid, "  num coeffs ", len(powers))
 
@@ -121,7 +123,7 @@ def fit_toroid1(xin,yin,z):
     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
+    vmask=~np.ma.getmaskarray(z)     # z.mask
     valid=z[vmask].size
     #coeffs=np.ones(5)       # 1, x, y, x^2, y^2
     powers=((0,0),(1,0),(0,1),(2,0),(0,2))