Advanced Interferogram Processing#

Here we will go over some of the more advanced interferometer data processing methods available in prysm. Many unrelated techniques will be covered here. A “master interferogram” is going to be the starting point for each process. This also demonstrates how you can create checkpoints in your own processing routines, if you wish.

As always, we begin with a few imports.

[1]:
import numpy as np

from prysm.interferogram import Interferogram
from prysm.sample_data import sample_files
from prysm.geometry import circle

Now we’re going to make the master dataset,

[2]:
path = sample_files('dat')
master = Interferogram.from_zygo_dat(path)
master.recenter()
master.mask(circle(20, master.r))
master.crop()
master.remove_piston()
master.remove_tiptilt()
master.plot2d()
[2]:
(<Figure size 640x480 with 2 Axes>, <AxesSubplot:>)
../_images/how-tos_Advanced-Interferogram-Processing_3_1.png

Two things should be noted here:

  • The area outside the clear aperture is filled with NaN values

  • There is a region with data dropout within the clear aperture

For reference, the PVr and RMS in units of nanometer are:

[3]:
master.pvr(), master.rms
[3]:
(112.99312233839987, 17.786794973956194)

PVr is a method because the user may wish to control the normalization radius used in the Zernike fit that is part of the definition of PVr. Before continuing, let’s look at all the things we can do with our interferogram:

[4]:
[s for s in dir(master) if not s.startswith('_')]
[4]:
['Sa',
 'bandlimited_rms',
 'copy',
 'crop',
 'data',
 'dropout_percentage',
 'dx',
 'exact_polar',
 'exact_x',
 'exact_xy',
 'exact_y',
 'fill',
 'filter',
 'from_zygo_dat',
 'intensity',
 'interferogram',
 'interpf_2d',
 'interpf_x',
 'interpf_y',
 'latcal',
 'mask',
 'meta',
 'pad',
 'plot2d',
 'psd',
 'pv',
 'pvr',
 'r',
 'recenter',
 'remove_piston',
 'remove_power',
 'remove_tiptilt',
 'render_from_psd',
 'rms',
 'save_zygo_ascii',
 'save_zygo_dat',
 'shape',
 'size',
 'slices',
 'spike_clip',
 'std',
 'strehl',
 'strip_latcal',
 'support',
 'support_x',
 'support_y',
 't',
 'total_integrated_scatter',
 'wavelength',
 'x',
 'y']

Some of these things (x,y,r,t) represent the coordinate grid. Some others (Sa, pv, PVr, rms, strehl, std, dropout_percentage, total_integrated_scatter) are statistical descriptions of the data. The low-order removal methods were already discussed. We have one alternative visualization method:

[5]:
master.interferogram(tilt_waves=(1,1))
[5]:
(<Figure size 640x480 with 2 Axes>, <AxesSubplot:>)
../_images/how-tos_Advanced-Interferogram-Processing_9_1.png

Some like to view these synthetic interferograms. The method allows the visibility, number of passes, and any extra tilt fringes to be controlled.

The first thing you may want to do is evaluate the bandlimited RMS value of the data. We can do so by first filling our NaNs with zero and then using the method. Here we’ll look in the 1 to 10 mm spatial period bandpass. Equivalent arguments are provided for frequencies, instead of periods.

[6]:
scratch = master.copy()
scratch.fill()
scratch.bandlimited_rms(1, 10)
[6]:
8.821682771742305

This value is in nanometers, and is roughly half the total RMS of our part. We can filter the data to the asme spatial period range and see that we get a similar answer:

[7]:
# filter only takes frequencies
scratch.filter((1/10, 1), typ='bandpass')
mask = np.isfinite(master.data)
scratch.mask(mask)
scratch.plot2d()
scratch.rms
/home/docs/checkouts/readthedocs.org/user_builds/prysm/envs/latest/lib/python3.8/site-packages/prysm/mathops.py:94: RuntimeWarning: invalid value encountered in divide
  out = special.j1(r) / r
[7]:
8.832251904041032
../_images/how-tos_Advanced-Interferogram-Processing_14_2.png

The value we get by this computation is a bit lower than the value we got with the bandlimited RMS function (about 15% lower). The reason for this is because spectral methods have finite out-of-band rejection. While prysm has significantly higher out of band rejection than the software sold with interferometers (> 60 dB higher), it is still finite, especially when the critical frequencies are near the lower or upper sampling limits of the data. We can view the PSD before and after filtering to see things more clearly:

[8]:
scratch2 = master.copy()
scratch2.fill()
psd_no_filter = scratch2.psd()

fig, ax = psd_no_filter.slices().plot('azavg')
scratch.fill()
psd_filter = scratch.psd()
psd_filter.slices().plot('azavg', fig=fig, ax=ax)
ax.set(xlabel='Spatial frequency, cy/mm', ylabel='PSD, nm^2/(cy/mm)^2', yscale='log', xscale='log', ylim=(1e-4,1e5), xlim=(1e-3,10))
ax.legend(['unfiltered', 'filtered'])
ax.grid(True)
../_images/how-tos_Advanced-Interferogram-Processing_16_0.png

In this case, we can see about three orders of magnitude rejection in both out-of-band regions. This would be considerably larger if the data had more samples (pixels), but the sample file is low resolution:

[9]:
print(master.shape)
(339, 339)

If we use only low or highpass filters far from the low and high frequency cutoffs, we can achieve stronger rejection:

[10]:
scratch = master.copy()
scratch.fill()
scratch.filter(0.1, typ='lp')
fig, ax = psd_no_filter.slices().plot('azavg')
scratch.psd().slices().plot('azavg', fig=fig,ax=ax)

ax.set(yscale='log', xscale='log', ylim=(1e-8,1e5), xlim=(1e-3,10))
ax.legend(['unfiltered', 'filtered'])
ax.grid(True)
ax.axvline(0.1)
[10]:
<matplotlib.lines.Line2D at 0x7f892ab19b20>
../_images/how-tos_Advanced-Interferogram-Processing_20_1.png

The small gain in power in the bandpass is a computational artifact (spectral leakage) and once again related to the low resolution of this interferogram. We can see a rejection from about 10^2 to 10^-7 by the time we reach 2x the cutoff frequency, or -80dB.

The last processing feature built into the Interferogram class is for spike clipping. This works the same way it does in MetroPro and Mx:

[11]:
scratch = master.copy()
scratch.spike_clip(3)  # 3 sigma is the default, too.
[11]:
<prysm.interferogram.Interferogram at 0x7f8928ec2400>

A thoughtful API for polynomial fitting as part of the interferogram interface has not been designed yet. If you strongly desire one, please do a design and submit a pull request on github. This does not mean polynomial fitting is not possible. Here we show fitting some low order Zernike polynomials,

[12]:
from prysm.polynomials import (
    fringe_to_nm,
    zernike_nm_sequence,
    lstsq,
    sum_of_2d_modes
)
from prysm.polynomials.zernike import barplot_magnitudes, zernikes_to_magnitude_angle

from prysm.util import rms
[13]:
r, t, data = master.r, master.t, master.data
normalization_radius = master.support/2
r = r / normalization_radius
fringe_indices = range(1,37)
nms = [fringe_to_nm(j) for j in fringe_indices]
modes = list(zernike_nm_sequence(nms, r, t))
fit = lstsq(modes, data)

pak = [[*nm, c] for nm, c in zip(nms, fit)]
magnitudes = zernikes_to_magnitude_angle(pak)
barplot_pak = {k: v[0] for k, v in magnitudes.items()}
barplot_magnitudes(barplot_pak)
[13]:
(<Figure size 640x480 with 1 Axes>, <AxesSubplot:>)
../_images/how-tos_Advanced-Interferogram-Processing_26_1.png

We can view the projection of various Zernike bandpasses:

[14]:
from matplotlib import pyplot as plt
low_order_projection = sum_of_2d_modes(modes[:10], fit[:10])
low_order_projection[~mask] = np.nan
plt.imshow(low_order_projection)
[14]:
<matplotlib.image.AxesImage at 0x7f8928dc2a00>
../_images/how-tos_Advanced-Interferogram-Processing_28_1.png
[15]:
mid_order_projection = sum_of_2d_modes(modes[10:22], fit[10:22])
mid_order_projection[~mask] = np.nan
plt.imshow(mid_order_projection)
[15]:
<matplotlib.image.AxesImage at 0x7f8928dac5e0>
../_images/how-tos_Advanced-Interferogram-Processing_29_1.png
[16]:
high_order_projection = sum_of_2d_modes(modes[22:], fit[22:])
high_order_projection[~mask] = np.nan
plt.imshow(high_order_projection)
[16]:
<matplotlib.image.AxesImage at 0x7f8928d17190>
../_images/how-tos_Advanced-Interferogram-Processing_30_1.png

As well as the total fit Zernike component:

[17]:
total_projection = sum_of_2d_modes(modes, fit)
total_projection[~mask] = np.nan
plt.imshow(total_projection)
[17]:
<matplotlib.image.AxesImage at 0x7f8928cf6d30>
../_images/how-tos_Advanced-Interferogram-Processing_32_1.png

And the fit error:

[18]:
fit_err_map = master.data - total_projection
plt.imshow(fit_err_map, clim=(-50,50), cmap='RdBu')
rms(fit_err_map) # nm
[18]:
5.832225080017871
../_images/how-tos_Advanced-Interferogram-Processing_34_1.png

We can do the same with other polynomial bases,

[19]:
from prysm.polynomials import Q2d_sequence
[20]:
modesQ = list(Q2d_sequence(nms, r, t))
fitQ = lstsq(modesQ, data)
[21]:
total_projection = sum_of_2d_modes(modesQ, fitQ)
total_projection[~mask] = np.nan
plt.imshow(total_projection)
[21]:
<matplotlib.image.AxesImage at 0x7f8928bd6520>
../_images/how-tos_Advanced-Interferogram-Processing_38_1.png
[22]:
fit_err_map = master.data - total_projection
plt.imshow(fit_err_map, clim=(-50,50), cmap='RdBu')
rms(fit_err_map) # nm
[22]:
11.600414095971875
../_images/how-tos_Advanced-Interferogram-Processing_39_1.png

We can see that the common polynomial framework of prysm made it trivial to swap out one polynomial basis for another.

As a final note, the metadata from the dat file is available in a python-friendly format:

[23]:
master.meta
[23]:
{'magic_number': 2283471727,
 'header_format': 1,
 'header_size': 834,
 'swtype': 1,
 'swdate': 'Mon Jan 04 08:42:09 2010',
 'swmaj': 8,
 'swmin': 3,
 'swpatch': 5,
 'ac_x': 0,
 'ac_y': 0,
 'ac_width': 640,
 'ac_height': 480,
 'ac_n_buckets': 1,
 'ac_range': 255,
 'ac_n_bytes': 614400,
 'cn_x': 139,
 'cn_y': 38,
 'cn_width': 425,
 'cn_height': 432,
 'cn_n_bytes': 734400,
 'timestamp': 1429905292,
 'comment': '',
 'source': 0,
 'scale_factor': 0.5,
 'wavelength': 6.327999813038332e-07,
 'numerical_aperture': 0.0,
 'obliquity_factor': 1.0,
 'magnification': 0.0,
 'lateral_resolution': 0.00011792420264100656,
 'acq_type': 0,
 'intensity_average_count': 0,
 'sfac_limit': 3,
 'ramp_cal': 0,
 'ramp_gain': 1753,
 'part_thickness': 0.0,
 'sw_llc': 1,
 'target_range': 0.10000000149011612,
 'rad_crv_measure_seq': 0,
 'min_mod': 17,
 'min_mod_count': 50,
 'phase_res': 1,
 'min_area': 20,
 'discontinuity_action': 1,
 'discontinuity_filter': 60.0,
 'connect_order': 0,
 'sign': 0,
 'camera_width': 640,
 'camera_height': 480,
 'sys_type': 23,
 'sys_board': 0,
 'sys_serial': 62398,
 'sys_inst_id': 0,
 'obj_name': 'Sm Aperture',
 'part_name': '',
 'codev_type': 0,
 'phase_avg_count': 10,
 'sub_sys_err': 0,
 'part_sn': '',
 'refractive_index': 1.0,
 'remove_tilt': 0,
 'remove_fringes': 0,
 'max_area': 9999999,
 'setup_type': 0,
 'wrapped': 0,
 'pre_connect_filter': 0.0,
 'wavelength_in_1': 6.327999813038332e-07,
 'wavelength_in_2': 0.0,
 'wavelength_in_3': 0.0,
 'wavelength_select': '1',
 'fda_res': 0,
 'scan_description': '',
 'n_fiducials': 0,
 'fiducial_1': 0.0,
 'fiducial_2': 0.0,
 'fiducial_3': 0.0,
 'fiducial_4': 0.0,
 'fiducial_5': 0.0,
 'fiducial_6': 0.0,
 'fiducial_7': 0.0,
 'fiducial_8': 0.0,
 'fiducial_9': 0.0,
 'fiducial_10': 0.0,
 'fiducial_11': 0.0,
 'fiducial_12': 0.0,
 'fiducial_13': 0.0,
 'fiducial_14': 0.0,
 'pixel_width': 7.3999999585794285e-06,
 'pixel_height': 7.3999999585794285e-06,
 'exit_pupil_diameter': 0.0,
 'light_level_percent': 55.49248504638672,
 'coords_state': 0,
 'coords_x': 0.0,
 'coords_y': 0.0,
 'coords_z': 0.0,
 'coords_a': 0.0,
 'coords_b': 0.0,
 'coords_c': 0.0,
 'cohrence_mode': 0,
 'surface_filter': 0,
 'sys_err_filename': '',
 'zoom_descr': '   1X ',
 'alpha_part': 0.0,
 'beta_part': 0.0,
 'dist_part': 0.0,
 'cam_split_loc_x': 0,
 'cam_split_loc_y': 0,
 'cam_split_trans_x': 0,
 'cam_split_trans_y': 0,
 'material_a': '',
 'material_b': '',
 'dmi_center_x': 0.0,
 'dmi_center_y': 0.0,
 'sph_distortion_correction': 0,
 'sph_dist_part_na': 0.0,
 'sph_dist_part_radius': 0.0,
 'sph_dist_cal_na': 0.0,
 'sph_dist_cal_radius': 0.0,
 'surface_type': 0,
 'ac_surface_type': 0,
 'z_pos': 0.0,
 'power_mul': 0.0,
 'focus_mul': 0.0,
 'roc_focus_cal_factor': 0.0,
 'roc_power_cal_factor': 0.0,
 'ftp_pos_left': 0.0,
 'ftp_pos_right': 0.0,
 'ftp_pos_pitch': 0.0,
 'ftp_pos_roll': 0.0,
 'min_mod_percent': 7.0,
 'max_intens': 0,
 'ring_of_fire': 0,
 'rc_orientation': b'\x00',
 'rc_distance': 0.0,
 'rc_angle': 0.0,
 'rc_diameter': 0.0,
 'rem_fringes_mode': 0,
 'ftpsi_phase_res': 0,
 'frames_acquired': 0,
 'cavity_type': 0,
 'cam_frame_rate': 0.0,
 'tune_range': 0.0,
 'cal_pix_x': 0,
 'cal_pix_y': 0,
 'test_cal_pts_1': 0.0,
 'test_cal_pts_2': 0.0,
 'test_cal_pts_3': 0.0,
 'test_cal_pts_4': 0.0,
 'ref_cal_pts_1': 0.0,
 'ref_cal_pts_2': 0.0,
 'ref_cal_pts_3': 0.0,
 'ref_cal_pts_4': 0.0,
 'test_cal_pix_opd': 0.0,
 'test_ref_pix_opd': 0.0,
 'flash_phase_cd_mask': 9.139576869988608e-40,
 'flash_phase_alias_mask': 0.0,
 'flash_phase_filter': 0.0,
 'scan_direction': 0,
 'ftpsi_res_factor': 0}

As well, the actual intensity camera data is available:

[24]:
plt.imshow(master.intensity)
[24]:
<matplotlib.image.AxesImage at 0x7f8929037f40>
../_images/how-tos_Advanced-Interferogram-Processing_43_1.png

Wrapping up, in this how-to we explored the various advanced processing routines for interferometer data present in prysm. We did not cover computing a PSF, MTF, or other downstream optical data products from the data. The .data and .dx attributes can be used to import the numerical data into the propagation routines of prysm. The facilities here can be combined to replace the software that comes with an interferometer to perform both basic and advanced processing alike.