Skip to content
Snippets Groups Projects
Commit 28060ac0 authored by Francois POLACK's avatar Francois POLACK
Browse files

Feat: 2D-PSD added to the available tools

parent 4ce72c21
Branches main
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
.gitignore
+ 1
− 0
View file @ 28060ac0
......@@ -8,6 +8,7 @@ doxygen/*
docs/
dabam/
dabam.*
newhtml/
#CB files
*.cbp
......
heightmap.py
+ 116
− 7
View file @ 28060ac0
......@@ -260,6 +260,31 @@ class HeightMap(object):
format(self.rawZ.shape[1], self.rawZ.shape[0], 1e3*size[0], 1e3*size[1], self.num_invalid))
return True
def read_h5file(self, filename=''):
if not Path(filename).is_file():
tkinter.Tk().withdraw()
filename=filedialog.askopenfilename(title="Open matlab data file", filetypes=( ("hdf5 file",("*.h5",".hdf5") ), ))
if Path(filename).suffix==".hdf5":
file = h5.File(filename,'r')
self.filepath=filename
self.x=self.rawX=np.copy(file['/Simu/x'])[::-1]
self.y=self.rawY=np.copy(file['/Simu/y'])[::-1]
self.ze=self.z=self.rawZ=np.ma.masked_invalid(file['/Simu/height'])
else:
return False
# print("x,y after closing file\n", self.x, '\n', self.y)
# print("heights\n", self.ze)
self.origin=[self.x[0], self.y[0]]
self.pixel=[abs(self.x[1]-self.x[0]), abs(self.y[1]-self.y[0])]
self.ROIorigin=(0,0)
self.ROIsize=[self.z.shape[1],self.z.shape[0]]
self.num_invalid=self.z[self.z.mask].size
size=(self.pixel[0]*self.rawZ.shape[1],self.pixel[1]*self.rawZ.shape[0])
print("Map size {:d} x {:d} pixels^2 dimension {:.1f} x {:.1f} mm^2 {:d} invalid data points".
format(self.rawZ.shape[1], self.rawZ.shape[0], 1e3*size[0], 1e3*size[1], self.num_invalid))
def set_scale(self, offset=(0,0)):
"""!
------- THIS function is no longer used -------------
......@@ -505,20 +530,21 @@ class HeightMap(object):
(self.sf_sag, self.srms_sag, self.hrms_sag, self.PSD_sag)=self.CRMSslope(axis=0)
def CRMSslope(self, axis=-1):
def CRMSslope(self, axis=-1, Z=None):
if(type(Z)==type(None)): Z=self.ze
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))
windata=spline_interpolate(Z, 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),
(sf, h2psd) = scipy.signal.periodogram(Z, 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)
print(" Size of the PSD ", h2psd.shape, " size of input ", Z.shape)
# uncomment to reverse cumulation direction
# s1csd = np.cumsum(s1psd[::-1])[::-1] * (sf[1] - sf[0])
# ajout 20/03/24
......@@ -530,6 +556,40 @@ class HeightMap(object):
print("Total cumulated rms=", srms[srms.shape[0]-1]*1e6, "µrad")
return (sf,srms, hrms,h1psd)
def PSD_2D(self, Z=None, winsize=(20,5)):
if(type(Z)==type(None)): Z=self.ze
fmax = (0.5 / self.pixel[0], 0.5/ self.pixel[1])
# we convolute both direction using alwaysspline_interpolate
windata=spline_interpolate(Z, rowwise=True, window=('tukey', 0.2))
Z=spline_interpolate(windata, rowwise=False, window=('tukey', 0.2))
# apply 2D Fourier transform
psize=(int((Z.shape[0]+2)/2), int((Z.shape[1]+2)/2))
spect=np.square(np.abs(scipy.fft.fft2(Z, norm="ortho")))[:psize[0], :psize[1]]
spect[1:-1, 1:-1]*=4.
spect[0, 1:-1]*=2.
spect[1:-1, 0]*=2.
fsag=np.linspace(0, fmax[0], psize[0])
ftan=np.linspace(0, fmax[1], psize[1])
# here we have the 2D periodogram
# We now average the X frequencies on a sliding window of size winsize
cs=np.zeros(ftan.shape[0]+1)
cs[1:]=ftan.cumsum()
fftan=(cs[winsize[0]:]-cs[:-winsize[0]])/winsize[0]
cs=np.zeros(fsag.shape[0]+1)
cs[1:]=fsag.cumsum()
ffsag=(cs[winsize[1]:]-cs[:-winsize[1]])/winsize[1]
# print('fftan',fftan)
# print('ffsag',ffsag)
cs=np.zeros((spect.shape[0], spect.shape[1]+1))
cs[:, 1:]=spect.cumsum(axis=1)
fspect=(cs[:,winsize[0]:]-cs[:,:-winsize[0]])/winsize[0]
cs=np.zeros((fspect.shape[0]+1, fspect.shape[1]))
cs[1:,:]=fspect.cumsum(axis=0)
fspect=(cs[winsize[1]:,:]-cs[:-winsize[1],:])/winsize[1]
return (ftan, fsag, spect, fftan, ffsag, fspect)
def print_statistics(self, height = True, raw = False, percent=0):
"""! Print statistics excluding the specified percentage of outlier pixels.
@param height (boolean) If **True**, prints height statistics, otherwise (**False**) prints slope statistics
......@@ -858,7 +918,7 @@ class HeightMap(object):
def save(self, filepath=''):
if not Path(filepath).parent.is_dir():
tkinter.Tk().withdraw()
filepath=filedialog.asksavefilename(title="open a hmap file", filetypes=(("superflat hmap","*.hmap"), ("all","*.*") ) )
filepath=filedialog.asksavefilename(title="save a hmap file", filetypes=(("superflat hmap","*.hmap"), ("all","*.*") ) )
print("writing to", filepath)
with open(filepath, 'wb') as f:
pickle.dump(self, f)
......@@ -920,7 +980,7 @@ class HeightMap(object):
print ('ROI origin =(', self.ROIorigin[0], ', ', self.ROIorigin[1], ')')
else:
origin_x=0
origin_y=rawZ.shape[0]-1
origin_y=self.rawZ.shape[0]-1
# Create raw_data
xraw=np.linspace(-origin_x*xpix,(self.rawZ.shape[1]-origin_x)*xpix,num=self.rawZ.shape[1], endpoint=False)
yraw=np.linspace(-origin_y*ypix,(self.rawZ.shape[0]-origin_y)*ypix,num=self.rawZ.shape[0], endpoint=False)
......@@ -1042,3 +1102,52 @@ def load(filepath=''):
print("loading from", filepath)
with open(filepath, 'rb') as f:
return pickle.load(f)
# Changed NAMING error degree nx vs ny
def Legendre_fit(nx, ny, datin, maxdeg=-1):
import sys
np.set_printoptions(precision=3, threshold=sys.maxsize, linewidth=150)
# compute mask of valid Legendres
Lmask=np.ones((nx+1,ny+1)); # changed
if maxdeg>0:
nt=nx+ny-maxdeg+1
Lmask[nx-nt+1:,ny-nt+1:]=np.triu(np.ones((nt, nt)))[:, ::-1] # changed
Lmask[:3,:3]=np.tril(np.ones((3,3)))[:, ::-1]
print('Lmask\n', Lmask)
print("data in" , datin.shape)
print("RMS", np.sqrt(np.square(datin).mean()))
xin=np.linspace(-1.,1., datin.shape[1])
yin=np.linspace(-1,1,datin.shape[0])
x, y = np.meshgrid(xin, yin)
print("x & y shape",x.shape)
V=np.polynomial.legendre.legvander2d(x, y, [nx,ny]) # changed
print("fitmat", V.shape)
A=V.reshape(V.shape[0]*V.shape[1], V.shape[2])
print(A.shape)
# print(A)
# Compute Legendre normalization matrix
cx=1./np.sqrt(2*np.linspace(0,nx,nx+1)+1)
cy=1./np.sqrt(2*np.linspace(0,ny,ny+1)+1)
print(cx, "\n", cy)
Mxy=cx.reshape(nx+1,1)@cy.reshape(1,ny+1) # changed 2
# print("Mxy\n", Mxy)
#compute fit
sol= np.linalg.lstsq(A, datin.ravel(), rcond=None)
# normalize solution
S=sol[0].reshape(nx+1,ny+1)*Mxy # changed
fitres=datin-V@sol[0].ravel()
print("RMS total fit",np.sqrt(np.square(S).sum()),
" RMS residuals", np.sqrt(np.square(fitres).mean()))
# imprime les coefficients du fit coef Y en ligne , coeff de X en colonne dan a disposintion eigen << list
print("fit\n",(S*Lmask)) #.transpose())
fitmap=V@(sol[0]*Lmask.ravel())
print("L sum", np.sqrt(np.square(S*Lmask).sum()),
"fitmap rms", np.sqrt(np.square(fitmap).mean()))
res=datin-fitmap
print("RMS residual", np.sqrt(np.square(res).mean()))
print("residuals", res.shape)
return (fitmap, res)
processing.py
+ 18
− 18
View file @ 28060ac0
......@@ -307,7 +307,7 @@ def fit_conic_cylinder(xin,yin,z, guess=None,variability=(0,1,1,1,0,1)):
"""! fit the Z function tabulated with grid values xin and yin, by a conic cylinder defined by a tuplet ( 1/p,1/q, theta, x0,y0, alpha)
@param xin the value of X sampling of the array Z
@param yin the value of Y sampling of the array Z
@para z the arrays of values to fit with a conic cylinder
@param z the arrays of values to fit with a conic cylinder
@param guess the tuplet of 6 initial values ( 1/p,1/q, theta, x0,y0, alpha) where
p is a the signed distance of entrance focus to apex point, q is the exist distance of exit focus to apex,
theta is the grazing angle of focal segments to the tangent at apex
......
test_legendre.py 0 → 100644
+ 253
− 0
View file @ 28060ac0
# -*- coding: utf-8 -*-
## @file test_Legendre.py
# This file is is a supplement to the Superflat Evaluation Toolkit
#
# This file is a first draft of a toolset aimed at analyzing a surface in termes of
# - Legendre polynomials for the low frequency part, and
# - 2D-PSD for the high frequency parts
#
# It is still in a early development stage and is given as it is, in case it might be useful
__author__ = "François Polack"
__copyright__ = "2024, Superflat-PCP project"
__email__ = "superflat-pcp@synchrotron-soleil.fr"
__version__ = "0.1.0"
import sys,warnings
# sys.path.append('D:\projets\ProjetsPython\Superflat')
import tkinter
from tkinter import filedialog
import heightmap
from pathlib import Path
import matplotlib.pyplot as plt
import matplotlib.figure as fig
import numpy as np
# from scipy import ndimage as img
#from edit_results import *
from crmsfig import create_crms_figure
def plot_map(X, Y, Z):
warnings.simplefilter('ignore', category=UserWarning)
errmap=plt.figure(figsize=fig.figaspect(Z), constrained_layout=True)
errmap.gca().invert_yaxis()
bounds=np.array([0.1, 99.9])
vbounds=np.percentile(Z,bounds)
plt.pcolormesh(X, Y, Z, cmap=plt.get_cmap('inferno'),
vmin=vbounds[0], vmax=vbounds[1], shading='auto')
warnings.resetwarnings()
cbar = plt.colorbar(aspect=25)
zlabel = "height (nm)"
cbar.set_label(zlabel) #, rotation=270)
plt.xlabel('x (mm)')
plt.ylabel('y (mm)')
return errmap
if __name__ == "__main__":
tkinter.Tk().withdraw()
filepath=filedialog.askopenfilename(title="open a hmap file", filetypes=(("superflat hmap","*.hmap"), ("all","*.*") ) )
if filepath!='' :
hmp=heightmap.load(filepath)
if hmp!= None:
print("data dims:", hmp.ze.shape)
print('pixel',hmp.pixel)
print(" size=", np.multiply(hmp.pixel, np.subtract(hmp.ze.shape,1.)))
# attention ordre degré y, degré x (ne pas se fier aux noms de ariable
(legfit,legres)=heightmap.Legendre_fit(10,10, hmp.ze, maxdeg=10) #8,6
plot_map(hmp.x*1e3, hmp.y*1e3, hmp.ze*1e9)
plt.title(Path(filepath).stem+" height error")
plot_map(hmp.x*1e3, hmp.y*1e3, legfit*1e9)
plt.title(Path(filepath).stem+" Legendre fit")
plot_map(hmp.x*1e3, hmp.y*1e3, legres*1e9)
plt.title(Path(filepath).stem+" residuals")
# ************ PROFILS TANGENTIELS **************
# *************************
(ftan, fsag,spect, fftan, ffsag, fspect)=hmp.PSD_2D(legres,(1,5))
# print('fspectrum size', fspect.shape,
# '\ntan size', fftan.shape[0], 'sag size ', ffsag.shape[0])
mp=plt.figure()
plt.loglog()
plt.grid(visible=True, which='both', axis='both')
plt.pcolormesh(fftan[1:],ffsag, fspect[:,1:], norm='log', cmap=plt.get_cmap('seismic'), shading='auto')
cbar = plt.colorbar(aspect=25)
plt.contour(fftan, ffsag, np.log10(fspect), 5, colors='k')
fvalues=(100 ,450, 800, 1000, 2000,4000)
sagcuts=np.searchsorted(ffsag, fvalues)
graph=create_crms_figure(figsize=fig.figaspect(0.5)) # new plot window
graph.alphcoeff=1
plt.loglog()
for cut in sagcuts :
plt.plot(fftan*1e-3, fspect[cut,:]*1e21, label='f_y {:.3f}'.format(ffsag[cut]*1e-3), lw=2, ls='-')
plt.grid(visible=True, which='both', axis='both')
plt.gca().set_xlim(left=0.05, right=6)
plt.xlabel(r'spatial frequency $(mm^{-1})$')
plt.ylabel('PSD $(nm^{2}/mm^{-1})$')
plt.title('Tangential PSDs {}'.format(''))
plt.legend()
# ** tangential PSD **
(sf,srms, hrms,h1psd)=hmp.CRMSslope(axis=-1, Z=legres)
# graph=plt.figure(figsize=fig.figaspect(0.5)) # new plot window
graph=create_crms_figure(figsize=fig.figaspect(0.5)) # new plot window
graph.alphcoeff=1
plt.plot(sf*1e-3, h1psd*1e21, label='Sagd avg',color='blue', lw=2, ls='-')
plt.plot(fftan[1:]*1e-3, fspect[0,1:]*1e21, label='local sag',color='red', lw=2, ls='-')
plt.grid(visible=True, which='both', axis='both')
plt.loglog()
plt.xlabel(r'spatial frequency $(mm^{-1})$')
plt.ylabel('PSD $(nm^{2}/mm^{-1})$')
plt.title('Tangential PSD {}'.format(''))
plt.legend()
# ************ PROFILS SAGITTAUX **************
# *************************
(ftan, fsag,spect, fftan, ffsag, fspect)=hmp.PSD_2D(legres,(20,1))
# print('ffsag',ffsag)
mp=plt.figure()
plt.loglog()
plt.grid(visible=True, which='both', axis='both')
plt.pcolormesh(fftan, ffsag[1:], fspect[1:,:], norm="log", cmap=plt.get_cmap('seismic'), shading='auto')
cbar = plt.colorbar(aspect=25)
plt.contour(fftan, ffsag, np.log10(fspect), 5, colors='k')
fvalues=(100, 150, 280 ,540, 800, 1500,4000)
tancuts=np.searchsorted(fftan, fvalues)
graph=create_crms_figure(figsize=fig.figaspect(0.5)) # new plot window
graph.alphcoeff=1
plt.loglog()
for cut in tancuts :
plt.plot(ffsag*1e-3, fspect[:, cut]*1e24, label='f_y {:.3f}'.format(fftan[cut]*1e-3), lw=2, ls='-')
plt.grid(visible=True, which='both', axis='both')
plt.gca().set_xlim(left=0.05, right=6)
plt.xlabel(r'spatial frequency $(mm^{-1})$')
plt.ylabel('PSD $(nm^{2}/mm^{-1})$')
plt.title('Sagittal PSDs {}'.format(''))
plt.legend()
# ** Sagittal CRMS height **
(sf,srms, hrms,h1psd)=hmp.CRMSslope(axis=0, Z=legres)
graph=create_crms_figure(figsize=fig.figaspect(0.5)) # new plot window
graph.alphcoeff=1
plt.plot(sf*1e-3, h1psd*1e21, label='tan avg',color='blue', lw=2, ls='-')
plt.plot(ffsag[1:]*1e-3, fspect[1:,0]*1e21, label='local avgg',color='red', lw=2, ls='-')
plt.grid(visible=True, which='both', axis='both')
plt.loglog()
plt.xlabel(r'spatial frequency $(mm^{-1})$')
plt.ylabel('PSD $(nm^{2}/mm^{-1})$')
plt.title('Sagittal PSD {}'.format(''))
plt.legend()
# ********** Carte 2D *********
# ***************
(ftan, fsag,spect, fftan, ffsag, fspect)=hmp.PSD_2D(legres,(10,3))
padspect=np.pad(fspect,((1,0),(1,0)), mode='edge')
ft=np.pad(fftan,(1,0),mode='constant', constant_values=ftan[1])
fs=np.pad(ffsag,(1,0),mode='constant', constant_values=fsag[1])
# print('fs', fs, "\n", ffsag[0], fsag[0])
mp=plt.figure()
plt.loglog()
plt.grid(visible=True, which='both', axis='both')
plt.pcolormesh(ft*1e-3,fs*1e-3, padspect*1e24, norm='log', cmap=plt.get_cmap('rainbow'), shading='auto')
cbar = plt.colorbar(aspect=25)
cbar.set_label(r'PSD $(nm^{2}mm^{2})$')
plt.xlabel(r'Tangential frequency $(mm^{-1})$')
plt.ylabel(r'Sagittal frequency $(mm^{-1})$')
plt.contour(ft*1e-3, fs*1e-3, np.log10(padspect*1e21),(-1,1,3,5), colors='k')
# *************************************************************************************
# mp=plt.figure()
# plt.contour(fftan, ffsag, np.log10(fspect), 20, color='k') # , cmap=plt.get_cmap('plasma'), shading='auto')
# cbar = plt.colorbar(aspect=25)
# plt.loglog()
# (sf,srms, hrms,h1psd)=hmp.CRMSslope(axis=-1, Z=legres)
# # ** tangential CSrms ***
# # graph=plt.figure(figsize=fig.figaspect(0.5)) # new plot window
# graph=create_crms_figure(figsize=fig.figaspect(0.5)) # new plot window
# plt.plot(hmp.sf*1e-3, hmp.srms*1e6, label='det 2°deg',color='blue', lw=2, ls='-')
# plt.plot(sf*1e-3, srms*1e6, label='det Leg',color='red', lw=2, ls='-')
# plt.grid(visible=True, which='major', axis='both')
# plt.loglog()
# plt.xlabel(r'spatial frequency $(mm^{-1})$')
# plt.ylabel('rms slope (µrad)')
# plt.title('Cumulated RMS tangential slopes {}'.format(''))
# plt.legend()
# # ** Tangential CRMS height **
# # graph=plt.figure(figsize=fig.figaspect(0.5)) # new plot window
# from crmsfig import create_crms_figure
# graph=create_crms_figure(figsize=fig.figaspect(0.5)) # new plot window
# plt.plot(hmp.sf*1e-3, hmp.hrms*1e9, label='det 2°deg',color='blue', lw=2, ls='-')
# plt.plot(sf*1e-3, hrms*1e9, label='det Leg',color='red', lw=2, ls='-')
# plt.grid(visible=True, which='both', axis='both')
# plt.loglog()
# plt.xlabel(r'spatial frequency $(mm^{-1})$')
# plt.ylabel('PSD (nm)')
# plt.title(' Tangential height CRMS {}'.format(''))
# plt.legend()
# (sf,srms, hrms,h1psd)=hmp.CRMSslope(axis=0, Z=legres)
# # ** Sagittal CRMS height **
# graph=create_crms_figure(figsize=fig.figaspect(0.5)) # new plot window
# plt.plot(sf*1e-3, hmp.hrms_sag*1e9, label='det 2°deg',color='blue', lw=2, ls='-')
# plt.plot(sf*1e-3, hrms*1e9, label='det Leg',color='red', lw=2, ls='-')
# plt.grid(visible=True, which='both', axis='both')
# plt.loglog()
# plt.xlabel(r'spatial frequency $(mm^{-1})$')
# plt.ylabel('PSD (nm)')
# plt.title(' Sagittal height CRMS {}'.format(''))
# plt.legend()
# # ** sagittal PSD **
# # graph=plt.figure(figsize=fig.figaspect(0.5)) # new plot window
# graph=create_crms_figure(figsize=fig.figaspect(0.5)) # new plot window
# graph.alphcoeff=1
# plt.plot(sf*1e-3, hmp.PSD_sag*1e21, label='det 2°deg',color='blue', lw=2, ls='-')
# plt.plot(sf*1e-3, h1psd*1e21, label='det Leg',color='red', lw=2, ls='-')
# plt.grid(visible=True, which='both', axis='both')
# plt.loglog()
# plt.xlabel(r'spatial frequency $(mm^{-1})$')
# plt.ylabel('PSD $(nm^{2}/mm^{-1})$')
# plt.title('Sagittal PSD {}'.format(''))
# plt.legend()
plt.show()
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment

The place to be for your source code