Loading ASONG_software_Manual.pdf 0 → 100644 +7.29 MiB File added.No diff preview for this file type. View file asong_slope.pro 0 → 100644 +149 −0 Original line number Diff line number Diff line ;+ Author: Laura Schreiber - August 2022 ; ; ROUTINE NAME: ;asong_slope ; ;PURPOSE: ;Given a frame and the binary masks defining the four quadrants of the ASONG WFS, compute the wavefront derivatives in X and Y. ; ;CALLING SEQUENCE: ;asong_slope, frame, qmask, pup, centres, sx, sy ; ;INPUT: ;image: measured frame with subtracted background and dovided by the flat (if required) ;qmask: 3D array with four planes, each representing a binary mask ; associated to one quadrant of the sensor ;pup: footprint of sx and sy ;centres: 4 X 2-elements vector of coordinates of pupils centroids ; ;KEYWORDS: ;gain: when this keyword is set, the gain is not computed and the value set in the keyword is used instead. ; It can come from calibration procedures. DEFAULT VALUE = 0 ;pyr: Set this keyword in case a Pyramid WFS is used in stead of ASONG. DEFAULT VALUE = 0 ;dxf, dyf: GTF x and y period in mm. If only one of the two keywords is set, the period in X and Y is considered to be the same. ; DEFAULT VALUE = 9.05 mm ;foc: distance from the aperture and the fourier plane ;display: set this keyword to visualise the computed slopes in a window ;pup_scale: System magnification. Theoretical value is given by the ratio between the focalising optic focal length and ; and the MLA focal length. This parameter should be measured in the lab by reimaging in the pupil an object of known lenght ;det_px: detector pixel size in microns. DEFAULT VALUE = 2.4 ; ;OUTPUT: ;sx, sy: array of signals in X and Y ;- ; NOTES: this procedure is adapted for the MLA with 4 lenses. In case of more lenses, minor modifications needs to be done in the code ; ;dir = 'C:\Users\schreibl\IDLWorkspace\Default\ASONG_software\TEST\' -- directory of your test folder ;name = 'test_' ;restore, dir + 'ASONG_mask_test.sav' ;restore,dir + 'test_image.sav' ;asong_slope, ima_avg, qmask, pup, centres, sx, sy ,/disp,/save, dir = dir, name = name PRO asong_slope, image, qmask, pup, centres, $ ;inputs GAIN = GAIN , $ dxf = dxf , $ dyf = dyf , $ foc = foc , $ display = display , $ pup_scale = pup_scale , $ det_px = det_px , $ save_file = save_file , $ dir = dir , $ wscale = wscale , $ name = name , $ ;keywords sx, sy ;outputs ;on_error, 2 if keyword_set(save_file) and not keyword_set(dir) then begin print, 'files will be saved in your working directory' pwd endif if not keyword_set(wscale) then wscale = 1 if not keyword_set(name) then name = '' if not keyword_set(PYR) then begin if not keyword_set(foc) then foc = 500. ; focal lengh in mm if not keyword_set(dxf) then dxf = 9.05 ; period in mm if not keyword_set(pup_scale) then pup_scale = 12.10 ; pupil object to pupil image scale factor if not keyword_set(det_px) then det_px = 2.4 ; detector pixel size in microns dyf = dxf gx = dxf / foc / 2. / !dpi / sqrt(2.) * (det_px * pup_scale) gy = dyf / foc / 2. / !dpi / sqrt(2.) * (det_px * pup_scale) endif else begin if not keyword_set(GAIN) then begin gx = 1.0 & gy = 1.0 endif else begin gx = gain & gy = gain endelse endelse dimpup = (size(pup,/dim))[0] c=centres ; Extract four quadrants from the frame nquad = 4 quad = fltarr(dimpup, dimpup, nquad) lx = round(c[*,0]) - dimpup/2 ux = round(c[*,0]) + dimpup/2 ly = round(c[*,1]) - dimpup/2 uy = round(c[*,1]) + dimpup/2 for q = 0, nquad - 1 do begin quad[*,*,q] = (image * qmask[*,*,q])[lx[q]:ux[q]-1,ly[q]:uy[q]-1]*pup endfor ; Compute signal and derivative sx = fltarr(dimpup, dimpup) sy = fltarr(dimpup, dimpup) norm = total(quad, 3) quadmask = (pup and (norm ne 0)) and 1B w = where(quadmask) if keyword_set(PYR) then begin sx[w] = ((quad[*,*,3] + quad[*,*,1]) - (quad[*,*,2] + quad[*,*,0]))[w] sy[w] = ((quad[*,*,2] + quad[*,*,3]) - (quad[*,*,0] + quad[*,*,1]))[w] sx[w] = gx * sx[w] / norm[w] sy[w] = gy * sy[w] / norm[w] endif else begin sx[w] = ((quad[*,*,3] + quad[*,*,1]) - (quad[*,*,2] + quad[*,*,0]))[w];/visx sy[w] = ((quad[*,*,0] + quad[*,*,1]) - (quad[*,*,2] + quad[*,*,3]))[w];/visy sx[w] = sx[w] / norm[w] sy[w] = sy[w] / norm[w] if total(where(abs(sx[w]) gt 1.)) gt 0 then begin print,'WARNING! Detected' + strcompress(string(n_elements(where(abs(sx[w]) gt 1.))))+ ' infinite values in the slope computation' print,'Automathic sobstitution of infinite values...' sx[where(abs(sx*quadmask) gt 1.)] = signum((sx[w])[where(abs(sx[w]) gt 1.)]) endif if total(where(abs(sy[w]) gt 1.)) gt 0 then begin print,'WARNING! Detected' + strcompress(string(n_elements(where(abs(sy[w]) gt 1.))))+ ' infinite values in the slope computation' print,'Automathic sobstitution of infinite values...' sy[where(abs(sy*quadmask) gt 1.)] = signum((sy[w])[where(abs(sy[w]) gt 1.)]) endif sx[w] = asin(sx[w]) * gx sy[w] = asin(sy[w]) * gy endelse if keyword_set(display) then begin window, 0 , title = 'Slope X and Y', xs = wscale * dimpup*2, ys = wscale * dimpup tvscl,congrid([sx,sy],wscale *dimpup*2,wscale *dimpup) endif if keyword_set(save_file) then begin if keyword_set(dir) then begin save, sx, sy, pup, filename = dir + 'ASONG_slopes' +name + '.sav' WRITE_BMP, dir + 'ASONG_slopes' + name + '_preview.bmp', BYTSCL([sx,sy]) endif else begin save, sx, sy, pup, filename = 'ASONG_slopes' +name + '.sav' WRITE_BMP,'ASONG_slopes' + name + '_preview.bmp', BYTSCL([sx,sy]) endelse endif return end asong_wfe_reconstruction.pro 0 → 100644 +253 −0 Original line number Diff line number Diff line ;+ Author: Laura Schreiber - August 2022 ; ; ROUTINE NAME: ;asong_wfe_reconstruction ; ;PURPOSE: ;Given the Wavefront X and Y slopes and the pupil mask, it reconstructs the phase through iterative method proposed by ;F. Roddier & C. Roddier (“Wavefront reconstruction using iterative Fourier transforms,”Applied Optics, 30, 1325, (1991)) ;and through (low order) modal reconstruction using Zernike modes. ; ;CALLING SEQUENCE: ;asong_wfe_reconstruction, sx, sy, pup, wfe_it, wfe_zer, zer_coeff ; ;INPUT: ;sx : Wavefront slope X ;sy : Wavefront slope Y ;pup : footprint of sx and sy ; ;KEYWORDS: ;save_file = save_file : Variable type: LOGICAL. Use this keyword to save the mask and the output variables. ;save_rmat = save_rmat : Variable type: LOGICAL. Use this keyword to save the reconstruction matrix and the Zernike cube. ;dir = dir : Variable type: STRING. Set this keyword to define the directory in which the variables must be stored. ;name = name : Variable type: STRING. Set this keyword to set the name of the file. DEFAULT VALUE = 'ASONG_mask' ;NOTILT = NOTILT : Variable type: LOGICAL. Set this keyword to subract the avarage value to slope maps. DEFAULT VALUE = 1 ;DISPLAY = DISPLAY : Variable type: LOGICAL. Set this keyword to have displayed the wavefronts once computed. DEFAULT VALUE = 1 ;RMAT = RMAT : Variable type: 2D array. This keyword can be used to pass to the procedure the reconstruction matrix if already computed. ; If the reconstruction matrix RMAT is not available, this keywork can be used to have it as an output and use it for a second iteration. ;n_zer = n_zer : Variable type: SCALAR. Number of computed Zernike modes. DEFAULT VALUE = 20. ; WARNING!!!! Since the pupil resolution is very high (around 1000 X 1000 pixels) ; the nominal value of Zernike modes that can be used for the reconstruction is very high. The computation time and the memory required for an high ; number of Zernike modes can be prohibitive. If possible, use this only with low nomber of Zernike (less than few hundreds) ;sigma_cut = sigma_cut : Data more than this number of standard deviations from the median is ignored. Suggested values: 2.0 and up. DEFAULT VALUE = 3 ;calibration = calibration :If set, the slopes are multiplied by a coefficient computed during calibration operations. DEFAULT VALUE = 1 ;lambda = lambda : Variable type: SCALAR. DEFAULT VALUE: 0.633 microns ;waves = waves : Variable type: LOGICAL. If set, returns the values in waves ;NOZER = NOZER : Variable type: LOGICAL. If set, only iterative reconstruction will be performed and all concerning modal reconstruction ; will be ignored ;ONLY_ZER = ONLY_ZER : Variable type: LOGICAL. If set, only modal reconstruction will be performed and all concerning iterative reconstruction ; will be ignored ; ;OUTPUT: ;wfe_it : Calibrated reconstructed WAVEFRONT throught the method of R&R ;wfe_zer : Low order reconstructed WAVEFRONT throught modal reconstruction ;zer_coeff : Zernike coefficients coming from modal reconstruction ;OPTIONAL OUTPUTS: ;RMAT : See the keyword RMAT. Precomputed Reconstruction matrix ;zermat : 3D Array containing the Zernike modes normalised on the pupil. To be visualised they have to be multiplied by the pupil array ;- ; NOTES and DISCLAIMER: This procedures uses many libraries, personal and pubblic. Please, check for permissions and disclaimar before distribution ; ;dir = 'C:\Users\schreibl\IDLWorkspace\Default\ASONG_software\TEST\' -- directory of your test folder ;name = 'test_' ;restore, dir + 'ASONG_slopes_test.sav' ;restore,dir + 'test_image.sav' ;asong_wfe_reconstruction, sx, sy, pup, wfe_it, wfe_zer, zer_coeff, zermat, rmat = rmat,/sav, dir=dir, name = name ;MODIFICATION HISTORY ;29/09/2022 Added keyword ONLY_ZER and relative header PRO asong_wfe_reconstruction, sx, sy, pup , $ NOTILT = NOTILT , $ DISPLAY = DISPLAY , $ rmat = rmat , $ NOZER = NOZER , $ n_zer = n_zer , $ sigma_cut = sigma_cut , $ calibration = calibration , $ lambda = lambda , $ waves = waves , $ save_file = save_file , $ dir = dir , $ name = name , $ silent = silent , $ save_rmat = save_rmat , $ ONLY_ZER = ONLY_ZER , $ _EXTRA = extra , $ wfe_it, wfe_zer, zer_coeff, zermat ;on_error,2 if not keyword_set(n_zer) then n_zer = 20 if not keyword_set(calibration) then calibration = 1 if not keyword_set(sigma_cut) then sigma_cut = 3 if not keyword_set(lambda) then lambda = 0.633 ;wavelengh in microns if not keyword_set(name) then name = '' if n_elements(display) eq 0 then display = 1 if n_elements(notilt) eq 0 then notilt = 1 if keyword_set(save_file) and not keyword_set(dir) then begin print, 'files will be saved in your working directory' pwd endif if keyword_set(only_zer) and keyword_set(nozer) then begin print,'Keywords ONLY_ZER and NOZER set at the same time.' print,'At least one reconstruction method should be chosen.' print,'Return' return endif dimpup = (size(pup,/dim))[0] sx = sx * calibration & sy = sy * calibration ; removing tilt from slopes by subtracting the average value if keyword_set(notilt) then begin ww = where(pup) sx_ = fltarr(dimpup,dimpup) & sy_ = fltarr(dimpup,dimpup) sx_[ww]=sx[ww]-mean(sx[ww]) & sy_[ww]=sy[ww]-mean(sy[ww]) endif else begin sx_ = sx & sy_ = sy endelse ; Modal reconstruction in case required if not keyword_set(NOZER) then begin if keyword_set(rmat) then begin if n_elements(zermat) eq 0 then begin print,'Zernike cube not defined.' print,'Re-computing RMAT and Zernike cube...' endif if n_zer ne (size(rmat,/dim))[1] or n_elements(where(pup)) ne (size(rmat,/dim))[0]/2 then begin print,'Incompatible reconstruction matrix encountered.' print,'Re-computing RMAT...' rmat = 0 endif endif if not keyword_set(rmat) then begin ;or n_zer ne (size(rmat,/dim))[1] rmat = recmat(n_zer, dimpup, pup, zermat) endif else begin n_zer = (size(rmat,/dim))[1] endelse calc_wfe_ls, sx_, sy_, pup, rmat, zermat, wfe_zer , zer_coeff, notilt = notilt if not keyword_set(silent) then print, 'number of zernike modes computed: ', n_zer endif ; R&R wavefront reconstruction if not keyword_set(ONLY_ZER) then wfe_it = slope_to_wavefront_ft(sx_, sy_, pup) else wfe_it = 0 ; RMS and PtV computation if keyword_set(waves) then begin if not keyword_set(ONLY_ZER) then wfe_it = wfe_it / lambda if not keyword_set(NOZER) then wfe_zer = wfe_zer / lambda endif if not keyword_set(ONLY_ZER) then begin wfe_it_rms = stddev(wfe_it[where(pup)]) & wfe_it_PtV = max(wfe_it[where(pup)]) - min(wfe_it[where(pup)]) if not keyword_set(silent) then begin if not keyword_set(waves) then print, 'WFE RMS ='+ strcompress(string(wfe_it_rms)) + ' microns' else $ print, 'WFE RMS ='+ strcompress(string(wfe_it_rms)) + ' waves (@0.633 microns)' if not keyword_set(waves) then print, 'WFE PtV ='+ strcompress(string(wfe_it_PtV)) + ' microns' else $ print, 'WFE PtV ='+ strcompress(string(wfe_it_PtV)) + ' waves (@0.633 microns)' endif endif if not keyword_set(NOZER) then begin wfe_zer_rms = stddev(wfe_zer[where(pup)]) & wfe_zer_PtV = max(wfe_zer[where(pup)]) - min(wfe_zer[where(pup)]) endif if not keyword_set(ONLY_ZER) then begin wfe_diff_rms = stddev(wfe_it[where(pup)]-wfe_zer[where(pup)]) wfe_diff_PtV = max(wfe_it[where(pup)]-wfe_zer[where(pup)]) - min(wfe_it[where(pup)]-wfe_zer[where(pup)]) endif ;resistant_mean, wfe_it[where(pup)], sigma_cut, wfe_it_mean, wfe_it_rms_clip, _EXTRA = extra ;wfe_it_PtV_clip = max((wfe_it[where(pup)])(where(abs(wfe_it[where(pup)]-wfe_it_mean) lt 3 * wfe_it_rms_clip ))) - min((wfe_it[where(pup)])(where(abs(wfe_it[where(pup)]-wfe_it_mean) lt 3 * wfe_it_rms_clip ))) ;resistant_mean, wfe_zer[where(pup)], sigma_cut, wfe_zer_mean, wfe_zer_rms_clip, _EXTRA = extra ;wfe_zer_PtV_clip = max((wfe_zer[where(pup)])(where(abs(wfe_zer[where(pup)]-wfe_zer_mean) lt 3 * wfe_zer_rms_clip ))) - min((wfe_zer[where(pup)])(where(abs(wfe_zer[where(pup)]-wfe_zer_mean) lt 3 * wfe_zer_rms_clip ))) ; display options if keyword_set(display) then begin if not keyword_set(NOZER) and not keyword_set(ONLY_ZER) then begin im1 = image(wfe_it*pup,image_dim = [dimpup,dimpup],image_loc = [0,0], dim = [dimpup*3,dimpup*1.2],xrange = [0,dimpup*3], yrange = [-0.2*dimpup,dimpup]) im2 = image(wfe_zer*pup,/over,image_dim = [dimpup,dimpup],image_loc = [dimpup,0]) im3 = image(wfe_it*pup - wfe_zer*pup, /over, image_dim = [dimpup,dimpup],image_loc = [2*dimpup,0]) text_tit_it = text(0.1, 0.15, 'Reconstructed Wavefront Roddier & Roddier',/current) text_tit_zer = text(0.4, 0.15, ' Modal Wavefront (first ' +strcompress(string(n_zer))+' modes)',/current) text_tit_diff = text(0.7, 0.15, ' Difference',/current) if not keyword_set(waves) then text_rms_it = text(0.1,0.1,'WFE RMS ='+ strcompress(string(wfe_it_rms)) + ' $\mu m$',/current) else $ text_rms_it = text(0.1,0.1,'WFE RMS ='+ strcompress(string(wfe_it_rms)) + ' waves (@0.633 $\mu m$)',/current) if not keyword_set(waves) then text_ptv_it = text(0.1,0.05,'WFE PtV ='+ strcompress(string(wfe_it_ptv)) + ' $\mu m$',/current) else $ text_ptv_it = text(0.1,0.05,'WFE PtV ='+ strcompress(string(wfe_it_ptv)) + ' waves (@0.633 $\mu m$)',/current) if not keyword_set(waves) then text_rms_zer = text(0.4,0.1,'WFE RMS ='+ strcompress(string(wfe_zer_rms)) + ' $\mu m$',/current) else $ text_rms_zer = text(0.4,0.1,'WFE RMS ='+ strcompress(string(wfe_zer_rms)) + ' waves (@0.633 $\mu m$)',/current) if not keyword_set(waves) then text_ptv_zer = text(0.4,0.05,'WFE PtV ='+ strcompress(string(wfe_zer_ptv)) + ' $\mu m$',/current) else $ text_ptv_zer = text(0.4,0.05,'WFE PtV ='+ strcompress(string(wfe_zer_ptv)) + ' waves (@0.633 $\mu m$)',/current) if not keyword_set(waves) then text_rms_diff = text(0.7,0.1,'WFE RMS ='+ strcompress(string(wfe_diff_rms)) + ' $\mu m$',/current) else $ text_rms_diff = text(0.7,0.1,'WFE RMS ='+ strcompress(string(wfe_diff_rms)) + ' waves (@0.633 $\mu m$)',/current) if not keyword_set(waves) then text_ptv_diff = text(0.7,0.05,'WFE PtV ='+ strcompress(string(wfe_diff_ptv)) + ' $\mu m$',/current) else $ text_ptv_diff = text(0.7,0.05,'WFE PtV ='+ strcompress(string(wfe_diff_ptv)) + ' waves (@0.633 $\mu m$)',/current) if not keyword_set(waves) then p = plot(findgen(n_zer), zer_coeff,dim = [2000,1000],xtitle = 'zernike mode', ytitle = 'WFE rms [um]', title = 'Zernike Coefficients', thick = 3) else $ p = plot(findgen(n_zer), zer_coeff,dim = [2000,1000],xtitle = 'zernike mode', ytitle = 'WFE rms [waves (@0.633 $\mu m$)]', title = 'Zernike Coefficients', thick = 3) if keyword_set(save_file) then im3.save, dir + 'ASONG_wfe'+ name + '_preview.jpg' endif if keyword_set(NOZER) then begin im1 = image(wfe_it*pup,image_dim = [dimpup,dimpup],image_loc = [0,0], dim = [dimpup,dimpup*1.2],xrange = [0,dimpup], yrange = [-0.2*dimpup,dimpup]) if not keyword_set(waves) then text_rms_it = text(0.1,0.1,'WFE RMS ='+ strcompress(string(wfe_it_rms)) + ' $\mu m$',/current) else $ text_rms_it = text(0.1,0.1,'WFE RMS ='+ strcompress(string(wfe_it_rms)) + ' waves (@0.633 $\mu m$)',/current) if not keyword_set(waves) then text_ptv_it = text(0.1,0.05,'WFE PtV ='+ strcompress(string(wfe_it_ptv)) + ' $\mu m$',/current) else $ text_ptv_it = text(0.1,0.05,'WFE PtV ='+ strcompress(string(wfe_it_ptv)) + ' waves (@0.633 $\mu m$)',/current) if keyword_set(save_file) then im1.save, dir + 'ASONG_wfe'+ name + '_preview.jpg' endif if keyword_set(ONLY_ZER) then begin im1 = image(wfe_zer*pup,image_dim = [dimpup,dimpup],image_loc = [0,0], dim = [dimpup,dimpup*1.2],xrange = [0,dimpup], yrange = [-0.2*dimpup,dimpup]) if not keyword_set(waves) then text_rms_zer = text(0.1,0.1,'WFE RMS ='+ strcompress(string(wfe_zer_rms)) + ' $\mu m$',/current) else $ text_rms_zer = text(0.1,0.1,'WFE RMS ='+ strcompress(string(wfe_zer_rms)) + ' waves (@0.633 $\mu m$)',/current) if not keyword_set(waves) then text_ptv_it = text(0.1,0.05,'WFE PtV ='+ strcompress(string(wfe_zer_ptv)) + ' $\mu m$',/current) else $ text_ptv_it = text(0.1,0.05,'WFE PtV ='+ strcompress(string(wfe_zer_ptv)) + ' waves (@0.633 $\mu m$)',/current) if keyword_set(save_file) then im1.save, dir + 'ASONG_wfe'+ name + '_preview.jpg' endif endif if not keyword_set(silent) then begin if not keyword_set(NOZER) then begin for i = 0,n_zer - 1 do begin print, 'zernike coefficient ' + strcompress(string(i)) + ' = ' + strcompress(string(zer_coeff[i])) endfor endif endif ;saving data if keyword_set(save_file) then begin if keyword_set(dir) then begin if not keyword_set(NOZER) or not keyword_set(ONLY_ZER) then begin if not keyword_set(ONLY_ZER) then save, wfe_it, wfe_zer, zer_coeff, pup, wfe_zer_rms, wfe_it_rms, wfe_zer_PtV, wfe_it_PtV, filename = dir + 'ASONG_wfe'+ name + '.sav' $ else save, wfe_zer, zer_coeff, pup, wfe_zer_rms, wfe_zer_PtV, filename = dir + 'ASONG_wfe'+ name + '.sav' if keyword_set(save_rmat) then save,rmat, zermat, pup, filename = dir + 'rmat' + name + '.sav' id = findgen(n_zer) openw, lun, dir + 'ASONG_Zernike_coeff'+ name + '.txt', /GET_LUN out=[transpose(id),transpose(zer_coeff)] if not keyword_set(waves) then printf,lun, '# id Value in microns' else printf,lun, '# id Value in waves (@0.633 microns)'' printf, lun, out, format='(I5,3x, F9.4)' & free_lun, lun endif if keyword_set(NOZER) then $ save, wfe_it, pup, wfe_it_rms, wfe_it_PtV, filename = dir + 'ASONG_wfe'+ name + '.sav' endif else begin if not keyword_set(NOZER) or not keyword_set(ONLY_ZER) then begin if not keyword_set(ONLY_ZER) then save, wfe_it, wfe_zer, zer_coeff, pup, wfe_zer_rms, wfe_it_rms, wfe_zer_PtV, wfe_it_PtV, filename = 'ASONG_wfe' + name + '.sav' $ else save, wfe_zer, zer_coeff, pup, wfe_zer_rms, wfe_zer_PtV, filename = 'ASONG_wfe' + name + '.sav' if keyword_set(save_rmat) then save, rmat, zermat, pup, filename = 'rmat' + name + '.sav' id = findgen(n_zer) openw, lun, 'ASONG_Zernike_coeff'+ name + '.txt', /GET_LUN out=[transpose(id),transpose(zer_coeff)] if not keyword_set(waves) then printf,lun, '# id Value in microns' else printf,lun, '# id Value in waves (@0.633 microns)'' printf, lun, out, format='(I5,3x, F9.4)' & free_lun, lun endif if keyword_set(NOZER) then save, wfe_it, pup, wfe_it_rms, wfe_it_PtV, filename = 'ASONG_wfe' + name + '.sav' endelse endif return end No newline at end of file Loading
ASONG_software_Manual.pdf 0 → 100644 +7.29 MiB File added.No diff preview for this file type. View file
asong_slope.pro 0 → 100644 +149 −0 Original line number Diff line number Diff line ;+ Author: Laura Schreiber - August 2022 ; ; ROUTINE NAME: ;asong_slope ; ;PURPOSE: ;Given a frame and the binary masks defining the four quadrants of the ASONG WFS, compute the wavefront derivatives in X and Y. ; ;CALLING SEQUENCE: ;asong_slope, frame, qmask, pup, centres, sx, sy ; ;INPUT: ;image: measured frame with subtracted background and dovided by the flat (if required) ;qmask: 3D array with four planes, each representing a binary mask ; associated to one quadrant of the sensor ;pup: footprint of sx and sy ;centres: 4 X 2-elements vector of coordinates of pupils centroids ; ;KEYWORDS: ;gain: when this keyword is set, the gain is not computed and the value set in the keyword is used instead. ; It can come from calibration procedures. DEFAULT VALUE = 0 ;pyr: Set this keyword in case a Pyramid WFS is used in stead of ASONG. DEFAULT VALUE = 0 ;dxf, dyf: GTF x and y period in mm. If only one of the two keywords is set, the period in X and Y is considered to be the same. ; DEFAULT VALUE = 9.05 mm ;foc: distance from the aperture and the fourier plane ;display: set this keyword to visualise the computed slopes in a window ;pup_scale: System magnification. Theoretical value is given by the ratio between the focalising optic focal length and ; and the MLA focal length. This parameter should be measured in the lab by reimaging in the pupil an object of known lenght ;det_px: detector pixel size in microns. DEFAULT VALUE = 2.4 ; ;OUTPUT: ;sx, sy: array of signals in X and Y ;- ; NOTES: this procedure is adapted for the MLA with 4 lenses. In case of more lenses, minor modifications needs to be done in the code ; ;dir = 'C:\Users\schreibl\IDLWorkspace\Default\ASONG_software\TEST\' -- directory of your test folder ;name = 'test_' ;restore, dir + 'ASONG_mask_test.sav' ;restore,dir + 'test_image.sav' ;asong_slope, ima_avg, qmask, pup, centres, sx, sy ,/disp,/save, dir = dir, name = name PRO asong_slope, image, qmask, pup, centres, $ ;inputs GAIN = GAIN , $ dxf = dxf , $ dyf = dyf , $ foc = foc , $ display = display , $ pup_scale = pup_scale , $ det_px = det_px , $ save_file = save_file , $ dir = dir , $ wscale = wscale , $ name = name , $ ;keywords sx, sy ;outputs ;on_error, 2 if keyword_set(save_file) and not keyword_set(dir) then begin print, 'files will be saved in your working directory' pwd endif if not keyword_set(wscale) then wscale = 1 if not keyword_set(name) then name = '' if not keyword_set(PYR) then begin if not keyword_set(foc) then foc = 500. ; focal lengh in mm if not keyword_set(dxf) then dxf = 9.05 ; period in mm if not keyword_set(pup_scale) then pup_scale = 12.10 ; pupil object to pupil image scale factor if not keyword_set(det_px) then det_px = 2.4 ; detector pixel size in microns dyf = dxf gx = dxf / foc / 2. / !dpi / sqrt(2.) * (det_px * pup_scale) gy = dyf / foc / 2. / !dpi / sqrt(2.) * (det_px * pup_scale) endif else begin if not keyword_set(GAIN) then begin gx = 1.0 & gy = 1.0 endif else begin gx = gain & gy = gain endelse endelse dimpup = (size(pup,/dim))[0] c=centres ; Extract four quadrants from the frame nquad = 4 quad = fltarr(dimpup, dimpup, nquad) lx = round(c[*,0]) - dimpup/2 ux = round(c[*,0]) + dimpup/2 ly = round(c[*,1]) - dimpup/2 uy = round(c[*,1]) + dimpup/2 for q = 0, nquad - 1 do begin quad[*,*,q] = (image * qmask[*,*,q])[lx[q]:ux[q]-1,ly[q]:uy[q]-1]*pup endfor ; Compute signal and derivative sx = fltarr(dimpup, dimpup) sy = fltarr(dimpup, dimpup) norm = total(quad, 3) quadmask = (pup and (norm ne 0)) and 1B w = where(quadmask) if keyword_set(PYR) then begin sx[w] = ((quad[*,*,3] + quad[*,*,1]) - (quad[*,*,2] + quad[*,*,0]))[w] sy[w] = ((quad[*,*,2] + quad[*,*,3]) - (quad[*,*,0] + quad[*,*,1]))[w] sx[w] = gx * sx[w] / norm[w] sy[w] = gy * sy[w] / norm[w] endif else begin sx[w] = ((quad[*,*,3] + quad[*,*,1]) - (quad[*,*,2] + quad[*,*,0]))[w];/visx sy[w] = ((quad[*,*,0] + quad[*,*,1]) - (quad[*,*,2] + quad[*,*,3]))[w];/visy sx[w] = sx[w] / norm[w] sy[w] = sy[w] / norm[w] if total(where(abs(sx[w]) gt 1.)) gt 0 then begin print,'WARNING! Detected' + strcompress(string(n_elements(where(abs(sx[w]) gt 1.))))+ ' infinite values in the slope computation' print,'Automathic sobstitution of infinite values...' sx[where(abs(sx*quadmask) gt 1.)] = signum((sx[w])[where(abs(sx[w]) gt 1.)]) endif if total(where(abs(sy[w]) gt 1.)) gt 0 then begin print,'WARNING! Detected' + strcompress(string(n_elements(where(abs(sy[w]) gt 1.))))+ ' infinite values in the slope computation' print,'Automathic sobstitution of infinite values...' sy[where(abs(sy*quadmask) gt 1.)] = signum((sy[w])[where(abs(sy[w]) gt 1.)]) endif sx[w] = asin(sx[w]) * gx sy[w] = asin(sy[w]) * gy endelse if keyword_set(display) then begin window, 0 , title = 'Slope X and Y', xs = wscale * dimpup*2, ys = wscale * dimpup tvscl,congrid([sx,sy],wscale *dimpup*2,wscale *dimpup) endif if keyword_set(save_file) then begin if keyword_set(dir) then begin save, sx, sy, pup, filename = dir + 'ASONG_slopes' +name + '.sav' WRITE_BMP, dir + 'ASONG_slopes' + name + '_preview.bmp', BYTSCL([sx,sy]) endif else begin save, sx, sy, pup, filename = 'ASONG_slopes' +name + '.sav' WRITE_BMP,'ASONG_slopes' + name + '_preview.bmp', BYTSCL([sx,sy]) endelse endif return end
asong_wfe_reconstruction.pro 0 → 100644 +253 −0 Original line number Diff line number Diff line ;+ Author: Laura Schreiber - August 2022 ; ; ROUTINE NAME: ;asong_wfe_reconstruction ; ;PURPOSE: ;Given the Wavefront X and Y slopes and the pupil mask, it reconstructs the phase through iterative method proposed by ;F. Roddier & C. Roddier (“Wavefront reconstruction using iterative Fourier transforms,”Applied Optics, 30, 1325, (1991)) ;and through (low order) modal reconstruction using Zernike modes. ; ;CALLING SEQUENCE: ;asong_wfe_reconstruction, sx, sy, pup, wfe_it, wfe_zer, zer_coeff ; ;INPUT: ;sx : Wavefront slope X ;sy : Wavefront slope Y ;pup : footprint of sx and sy ; ;KEYWORDS: ;save_file = save_file : Variable type: LOGICAL. Use this keyword to save the mask and the output variables. ;save_rmat = save_rmat : Variable type: LOGICAL. Use this keyword to save the reconstruction matrix and the Zernike cube. ;dir = dir : Variable type: STRING. Set this keyword to define the directory in which the variables must be stored. ;name = name : Variable type: STRING. Set this keyword to set the name of the file. DEFAULT VALUE = 'ASONG_mask' ;NOTILT = NOTILT : Variable type: LOGICAL. Set this keyword to subract the avarage value to slope maps. DEFAULT VALUE = 1 ;DISPLAY = DISPLAY : Variable type: LOGICAL. Set this keyword to have displayed the wavefronts once computed. DEFAULT VALUE = 1 ;RMAT = RMAT : Variable type: 2D array. This keyword can be used to pass to the procedure the reconstruction matrix if already computed. ; If the reconstruction matrix RMAT is not available, this keywork can be used to have it as an output and use it for a second iteration. ;n_zer = n_zer : Variable type: SCALAR. Number of computed Zernike modes. DEFAULT VALUE = 20. ; WARNING!!!! Since the pupil resolution is very high (around 1000 X 1000 pixels) ; the nominal value of Zernike modes that can be used for the reconstruction is very high. The computation time and the memory required for an high ; number of Zernike modes can be prohibitive. If possible, use this only with low nomber of Zernike (less than few hundreds) ;sigma_cut = sigma_cut : Data more than this number of standard deviations from the median is ignored. Suggested values: 2.0 and up. DEFAULT VALUE = 3 ;calibration = calibration :If set, the slopes are multiplied by a coefficient computed during calibration operations. DEFAULT VALUE = 1 ;lambda = lambda : Variable type: SCALAR. DEFAULT VALUE: 0.633 microns ;waves = waves : Variable type: LOGICAL. If set, returns the values in waves ;NOZER = NOZER : Variable type: LOGICAL. If set, only iterative reconstruction will be performed and all concerning modal reconstruction ; will be ignored ;ONLY_ZER = ONLY_ZER : Variable type: LOGICAL. If set, only modal reconstruction will be performed and all concerning iterative reconstruction ; will be ignored ; ;OUTPUT: ;wfe_it : Calibrated reconstructed WAVEFRONT throught the method of R&R ;wfe_zer : Low order reconstructed WAVEFRONT throught modal reconstruction ;zer_coeff : Zernike coefficients coming from modal reconstruction ;OPTIONAL OUTPUTS: ;RMAT : See the keyword RMAT. Precomputed Reconstruction matrix ;zermat : 3D Array containing the Zernike modes normalised on the pupil. To be visualised they have to be multiplied by the pupil array ;- ; NOTES and DISCLAIMER: This procedures uses many libraries, personal and pubblic. Please, check for permissions and disclaimar before distribution ; ;dir = 'C:\Users\schreibl\IDLWorkspace\Default\ASONG_software\TEST\' -- directory of your test folder ;name = 'test_' ;restore, dir + 'ASONG_slopes_test.sav' ;restore,dir + 'test_image.sav' ;asong_wfe_reconstruction, sx, sy, pup, wfe_it, wfe_zer, zer_coeff, zermat, rmat = rmat,/sav, dir=dir, name = name ;MODIFICATION HISTORY ;29/09/2022 Added keyword ONLY_ZER and relative header PRO asong_wfe_reconstruction, sx, sy, pup , $ NOTILT = NOTILT , $ DISPLAY = DISPLAY , $ rmat = rmat , $ NOZER = NOZER , $ n_zer = n_zer , $ sigma_cut = sigma_cut , $ calibration = calibration , $ lambda = lambda , $ waves = waves , $ save_file = save_file , $ dir = dir , $ name = name , $ silent = silent , $ save_rmat = save_rmat , $ ONLY_ZER = ONLY_ZER , $ _EXTRA = extra , $ wfe_it, wfe_zer, zer_coeff, zermat ;on_error,2 if not keyword_set(n_zer) then n_zer = 20 if not keyword_set(calibration) then calibration = 1 if not keyword_set(sigma_cut) then sigma_cut = 3 if not keyword_set(lambda) then lambda = 0.633 ;wavelengh in microns if not keyword_set(name) then name = '' if n_elements(display) eq 0 then display = 1 if n_elements(notilt) eq 0 then notilt = 1 if keyword_set(save_file) and not keyword_set(dir) then begin print, 'files will be saved in your working directory' pwd endif if keyword_set(only_zer) and keyword_set(nozer) then begin print,'Keywords ONLY_ZER and NOZER set at the same time.' print,'At least one reconstruction method should be chosen.' print,'Return' return endif dimpup = (size(pup,/dim))[0] sx = sx * calibration & sy = sy * calibration ; removing tilt from slopes by subtracting the average value if keyword_set(notilt) then begin ww = where(pup) sx_ = fltarr(dimpup,dimpup) & sy_ = fltarr(dimpup,dimpup) sx_[ww]=sx[ww]-mean(sx[ww]) & sy_[ww]=sy[ww]-mean(sy[ww]) endif else begin sx_ = sx & sy_ = sy endelse ; Modal reconstruction in case required if not keyword_set(NOZER) then begin if keyword_set(rmat) then begin if n_elements(zermat) eq 0 then begin print,'Zernike cube not defined.' print,'Re-computing RMAT and Zernike cube...' endif if n_zer ne (size(rmat,/dim))[1] or n_elements(where(pup)) ne (size(rmat,/dim))[0]/2 then begin print,'Incompatible reconstruction matrix encountered.' print,'Re-computing RMAT...' rmat = 0 endif endif if not keyword_set(rmat) then begin ;or n_zer ne (size(rmat,/dim))[1] rmat = recmat(n_zer, dimpup, pup, zermat) endif else begin n_zer = (size(rmat,/dim))[1] endelse calc_wfe_ls, sx_, sy_, pup, rmat, zermat, wfe_zer , zer_coeff, notilt = notilt if not keyword_set(silent) then print, 'number of zernike modes computed: ', n_zer endif ; R&R wavefront reconstruction if not keyword_set(ONLY_ZER) then wfe_it = slope_to_wavefront_ft(sx_, sy_, pup) else wfe_it = 0 ; RMS and PtV computation if keyword_set(waves) then begin if not keyword_set(ONLY_ZER) then wfe_it = wfe_it / lambda if not keyword_set(NOZER) then wfe_zer = wfe_zer / lambda endif if not keyword_set(ONLY_ZER) then begin wfe_it_rms = stddev(wfe_it[where(pup)]) & wfe_it_PtV = max(wfe_it[where(pup)]) - min(wfe_it[where(pup)]) if not keyword_set(silent) then begin if not keyword_set(waves) then print, 'WFE RMS ='+ strcompress(string(wfe_it_rms)) + ' microns' else $ print, 'WFE RMS ='+ strcompress(string(wfe_it_rms)) + ' waves (@0.633 microns)' if not keyword_set(waves) then print, 'WFE PtV ='+ strcompress(string(wfe_it_PtV)) + ' microns' else $ print, 'WFE PtV ='+ strcompress(string(wfe_it_PtV)) + ' waves (@0.633 microns)' endif endif if not keyword_set(NOZER) then begin wfe_zer_rms = stddev(wfe_zer[where(pup)]) & wfe_zer_PtV = max(wfe_zer[where(pup)]) - min(wfe_zer[where(pup)]) endif if not keyword_set(ONLY_ZER) then begin wfe_diff_rms = stddev(wfe_it[where(pup)]-wfe_zer[where(pup)]) wfe_diff_PtV = max(wfe_it[where(pup)]-wfe_zer[where(pup)]) - min(wfe_it[where(pup)]-wfe_zer[where(pup)]) endif ;resistant_mean, wfe_it[where(pup)], sigma_cut, wfe_it_mean, wfe_it_rms_clip, _EXTRA = extra ;wfe_it_PtV_clip = max((wfe_it[where(pup)])(where(abs(wfe_it[where(pup)]-wfe_it_mean) lt 3 * wfe_it_rms_clip ))) - min((wfe_it[where(pup)])(where(abs(wfe_it[where(pup)]-wfe_it_mean) lt 3 * wfe_it_rms_clip ))) ;resistant_mean, wfe_zer[where(pup)], sigma_cut, wfe_zer_mean, wfe_zer_rms_clip, _EXTRA = extra ;wfe_zer_PtV_clip = max((wfe_zer[where(pup)])(where(abs(wfe_zer[where(pup)]-wfe_zer_mean) lt 3 * wfe_zer_rms_clip ))) - min((wfe_zer[where(pup)])(where(abs(wfe_zer[where(pup)]-wfe_zer_mean) lt 3 * wfe_zer_rms_clip ))) ; display options if keyword_set(display) then begin if not keyword_set(NOZER) and not keyword_set(ONLY_ZER) then begin im1 = image(wfe_it*pup,image_dim = [dimpup,dimpup],image_loc = [0,0], dim = [dimpup*3,dimpup*1.2],xrange = [0,dimpup*3], yrange = [-0.2*dimpup,dimpup]) im2 = image(wfe_zer*pup,/over,image_dim = [dimpup,dimpup],image_loc = [dimpup,0]) im3 = image(wfe_it*pup - wfe_zer*pup, /over, image_dim = [dimpup,dimpup],image_loc = [2*dimpup,0]) text_tit_it = text(0.1, 0.15, 'Reconstructed Wavefront Roddier & Roddier',/current) text_tit_zer = text(0.4, 0.15, ' Modal Wavefront (first ' +strcompress(string(n_zer))+' modes)',/current) text_tit_diff = text(0.7, 0.15, ' Difference',/current) if not keyword_set(waves) then text_rms_it = text(0.1,0.1,'WFE RMS ='+ strcompress(string(wfe_it_rms)) + ' $\mu m$',/current) else $ text_rms_it = text(0.1,0.1,'WFE RMS ='+ strcompress(string(wfe_it_rms)) + ' waves (@0.633 $\mu m$)',/current) if not keyword_set(waves) then text_ptv_it = text(0.1,0.05,'WFE PtV ='+ strcompress(string(wfe_it_ptv)) + ' $\mu m$',/current) else $ text_ptv_it = text(0.1,0.05,'WFE PtV ='+ strcompress(string(wfe_it_ptv)) + ' waves (@0.633 $\mu m$)',/current) if not keyword_set(waves) then text_rms_zer = text(0.4,0.1,'WFE RMS ='+ strcompress(string(wfe_zer_rms)) + ' $\mu m$',/current) else $ text_rms_zer = text(0.4,0.1,'WFE RMS ='+ strcompress(string(wfe_zer_rms)) + ' waves (@0.633 $\mu m$)',/current) if not keyword_set(waves) then text_ptv_zer = text(0.4,0.05,'WFE PtV ='+ strcompress(string(wfe_zer_ptv)) + ' $\mu m$',/current) else $ text_ptv_zer = text(0.4,0.05,'WFE PtV ='+ strcompress(string(wfe_zer_ptv)) + ' waves (@0.633 $\mu m$)',/current) if not keyword_set(waves) then text_rms_diff = text(0.7,0.1,'WFE RMS ='+ strcompress(string(wfe_diff_rms)) + ' $\mu m$',/current) else $ text_rms_diff = text(0.7,0.1,'WFE RMS ='+ strcompress(string(wfe_diff_rms)) + ' waves (@0.633 $\mu m$)',/current) if not keyword_set(waves) then text_ptv_diff = text(0.7,0.05,'WFE PtV ='+ strcompress(string(wfe_diff_ptv)) + ' $\mu m$',/current) else $ text_ptv_diff = text(0.7,0.05,'WFE PtV ='+ strcompress(string(wfe_diff_ptv)) + ' waves (@0.633 $\mu m$)',/current) if not keyword_set(waves) then p = plot(findgen(n_zer), zer_coeff,dim = [2000,1000],xtitle = 'zernike mode', ytitle = 'WFE rms [um]', title = 'Zernike Coefficients', thick = 3) else $ p = plot(findgen(n_zer), zer_coeff,dim = [2000,1000],xtitle = 'zernike mode', ytitle = 'WFE rms [waves (@0.633 $\mu m$)]', title = 'Zernike Coefficients', thick = 3) if keyword_set(save_file) then im3.save, dir + 'ASONG_wfe'+ name + '_preview.jpg' endif if keyword_set(NOZER) then begin im1 = image(wfe_it*pup,image_dim = [dimpup,dimpup],image_loc = [0,0], dim = [dimpup,dimpup*1.2],xrange = [0,dimpup], yrange = [-0.2*dimpup,dimpup]) if not keyword_set(waves) then text_rms_it = text(0.1,0.1,'WFE RMS ='+ strcompress(string(wfe_it_rms)) + ' $\mu m$',/current) else $ text_rms_it = text(0.1,0.1,'WFE RMS ='+ strcompress(string(wfe_it_rms)) + ' waves (@0.633 $\mu m$)',/current) if not keyword_set(waves) then text_ptv_it = text(0.1,0.05,'WFE PtV ='+ strcompress(string(wfe_it_ptv)) + ' $\mu m$',/current) else $ text_ptv_it = text(0.1,0.05,'WFE PtV ='+ strcompress(string(wfe_it_ptv)) + ' waves (@0.633 $\mu m$)',/current) if keyword_set(save_file) then im1.save, dir + 'ASONG_wfe'+ name + '_preview.jpg' endif if keyword_set(ONLY_ZER) then begin im1 = image(wfe_zer*pup,image_dim = [dimpup,dimpup],image_loc = [0,0], dim = [dimpup,dimpup*1.2],xrange = [0,dimpup], yrange = [-0.2*dimpup,dimpup]) if not keyword_set(waves) then text_rms_zer = text(0.1,0.1,'WFE RMS ='+ strcompress(string(wfe_zer_rms)) + ' $\mu m$',/current) else $ text_rms_zer = text(0.1,0.1,'WFE RMS ='+ strcompress(string(wfe_zer_rms)) + ' waves (@0.633 $\mu m$)',/current) if not keyword_set(waves) then text_ptv_it = text(0.1,0.05,'WFE PtV ='+ strcompress(string(wfe_zer_ptv)) + ' $\mu m$',/current) else $ text_ptv_it = text(0.1,0.05,'WFE PtV ='+ strcompress(string(wfe_zer_ptv)) + ' waves (@0.633 $\mu m$)',/current) if keyword_set(save_file) then im1.save, dir + 'ASONG_wfe'+ name + '_preview.jpg' endif endif if not keyword_set(silent) then begin if not keyword_set(NOZER) then begin for i = 0,n_zer - 1 do begin print, 'zernike coefficient ' + strcompress(string(i)) + ' = ' + strcompress(string(zer_coeff[i])) endfor endif endif ;saving data if keyword_set(save_file) then begin if keyword_set(dir) then begin if not keyword_set(NOZER) or not keyword_set(ONLY_ZER) then begin if not keyword_set(ONLY_ZER) then save, wfe_it, wfe_zer, zer_coeff, pup, wfe_zer_rms, wfe_it_rms, wfe_zer_PtV, wfe_it_PtV, filename = dir + 'ASONG_wfe'+ name + '.sav' $ else save, wfe_zer, zer_coeff, pup, wfe_zer_rms, wfe_zer_PtV, filename = dir + 'ASONG_wfe'+ name + '.sav' if keyword_set(save_rmat) then save,rmat, zermat, pup, filename = dir + 'rmat' + name + '.sav' id = findgen(n_zer) openw, lun, dir + 'ASONG_Zernike_coeff'+ name + '.txt', /GET_LUN out=[transpose(id),transpose(zer_coeff)] if not keyword_set(waves) then printf,lun, '# id Value in microns' else printf,lun, '# id Value in waves (@0.633 microns)'' printf, lun, out, format='(I5,3x, F9.4)' & free_lun, lun endif if keyword_set(NOZER) then $ save, wfe_it, pup, wfe_it_rms, wfe_it_PtV, filename = dir + 'ASONG_wfe'+ name + '.sav' endif else begin if not keyword_set(NOZER) or not keyword_set(ONLY_ZER) then begin if not keyword_set(ONLY_ZER) then save, wfe_it, wfe_zer, zer_coeff, pup, wfe_zer_rms, wfe_it_rms, wfe_zer_PtV, wfe_it_PtV, filename = 'ASONG_wfe' + name + '.sav' $ else save, wfe_zer, zer_coeff, pup, wfe_zer_rms, wfe_zer_PtV, filename = 'ASONG_wfe' + name + '.sav' if keyword_set(save_rmat) then save, rmat, zermat, pup, filename = 'rmat' + name + '.sav' id = findgen(n_zer) openw, lun, 'ASONG_Zernike_coeff'+ name + '.txt', /GET_LUN out=[transpose(id),transpose(zer_coeff)] if not keyword_set(waves) then printf,lun, '# id Value in microns' else printf,lun, '# id Value in waves (@0.633 microns)'' printf, lun, out, format='(I5,3x, F9.4)' & free_lun, lun endif if keyword_set(NOZER) then save, wfe_it, pup, wfe_it_rms, wfe_it_PtV, filename = 'ASONG_wfe' + name + '.sav' endelse endif return end No newline at end of file
