<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Patrick provided a pair of images from AuxTel.
Let's look at how those images work with our cwfs code
load the modules
Step1:  Define the image objects. Input arguments
Step2:  Define the instrument. Input arguments
Step3:  Define the algorithm being used. Input arguments
Step4:  Run it
Step5:  Print the Zernikes Zn (n>=4)
Step6:  plot the Zernikes Zn (n>=4)
Step7:  We check that the optical parameters provided are consistent with the image diameter. Otherwise the numerical solutions themselves do not make much sense.
Step8:  Patrick asked the question
Step9:  Now we do the forward raytrace using our wavefront solutions
The code is simply borrowed from existing cwfs code.
We first set up the pupil grid. Oversample means how many ray to trace from each grid point on the pupil.
Step10:  We now trace the rays to the image plane. Lutxp and Lutyp are image coordinates for each (oversampled) ray. showProjection() makes the intensity image. Then, to down sample the image back to original resolution, we want to use the function downResolution() which is defined for the image class.
Step11:  Now do the same thing for extra focal image
<ASSISTANT_TASK:>
Python Code:
from lsst.cwfs.instrument import Instrument
from lsst.cwfs.algorithm import Algorithm
from lsst.cwfs.image import Image, readFile, aperture2image, showProjection
import lsst.cwfs.plots as plots
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
Explanation: Patrick provided a pair of images from AuxTel.
Let's look at how those images work with our cwfs code
load the modules
End of explanation
fieldXY = [0,0]
I1 = Image(readFile('../tests/testImages/AuxTel/I1_intra_20190912_HD21161_z05.fits'), fieldXY, Image.INTRA)
I2 = Image(readFile('../tests/testImages/AuxTel/I2_extra_20190912_HD21161_z05.fits'), fieldXY, Image.EXTRA)
plots.plotImage(I1.image,'intra')
plots.plotImage(I2.image,'extra')
Explanation: Define the image objects. Input arguments: file name, field coordinates in deg, image type
The colorbar() below may produce a warning message if your matplotlib version is older than 1.5.0
( https://github.com/matplotlib/matplotlib/issues/5209 )
End of explanation
inst=Instrument('AuxTel',I1.sizeinPix)
Explanation: Define the instrument. Input arguments: instrument name, size of image stamps
End of explanation
algo=Algorithm('exp',inst,0)
Explanation: Define the algorithm being used. Input arguments: baseline algorithm, instrument, debug level
End of explanation
algo.runIt(inst,I1,I2,'paraxial')
Explanation: Run it
End of explanation
print(algo.zer4UpNm)
Explanation: Print the Zernikes Zn (n>=4)
End of explanation
plots.plotZer(algo.zer4UpNm,'nm')
Explanation: plot the Zernikes Zn (n>=4)
End of explanation
print("Expected image diameter in pixels = %.0f"%(inst.offset/inst.fno/inst.pixelSize))
plots.plotImage(I1.image0,'original intra', mask=algo.pMask)
plots.plotImage(I2.image0,'original extra', mask=algo.pMask)
Explanation: We check that the optical parameters provided are consistent with the image diameter. Otherwise the numerical solutions themselves do not make much sense.
End of explanation
nanMask = np.ones(I1.image.shape)
nanMask[I1.pMask==0] = np.nan
fig, ax = plt.subplots(1,2, figsize=[10,4])
img = ax[0].imshow(algo.Wconverge*nanMask, origin='lower')
ax[0].set_title('Final WF = estimated + residual')
fig.colorbar(img, ax=ax[0])
img = ax[1].imshow(algo.West*nanMask, origin='lower')
ax[1].set_title('residual wavefront')
fig.colorbar(img, ax=ax[1])
fig, ax = plt.subplots(1,2, figsize=[10,4])
img = ax[0].imshow(I1.image, origin='lower')
ax[0].set_title('Intra residual image')
fig.colorbar(img, ax=ax[0])
img = ax[1].imshow(I2.image, origin='lower')
ax[1].set_title('Extra residual image')
fig.colorbar(img, ax=ax[1])
Explanation: Patrick asked the question: can we show the results of the fit in intensity space, and also the residual?
Great question. The short answer is no.
The long answer: the current approach implemented is the so-called inversion approach, i.e., to inversely solve the Transport of Intensity Equation with boundary conditions. It is not a forward fit. If you think of the unperturbed image as I0, and the real image as I, we iteratively map I back toward I0 using the estimated wavefront. Upon convergence, our "residual images" should have intensity distributions that are nearly uniform. We always have an estimated wavefront, and a residual wavefront. The residual wavefront is obtained from the two residual images.
However, using tools availabe in the cwfs package, we can easily make forward prediction of the images using the wavefront solution. This is basically to take the slope of the wavefront at any pupil position, and raytrace to the image plane. We will demostrate these below.
End of explanation
oversample = 10
projSamples = I1.image0.shape[0]*oversample
luty, lutx = np.mgrid[
        -(projSamples / 2 - 0.5):(projSamples / 2 + 0.5),
        -(projSamples / 2 - 0.5):(projSamples / 2 + 0.5)]
lutx = lutx / (projSamples / 2 / inst.sensorFactor)
luty = luty / (projSamples / 2 / inst.sensorFactor)
Explanation: Now we do the forward raytrace using our wavefront solutions
The code is simply borrowed from existing cwfs code.
We first set up the pupil grid. Oversample means how many ray to trace from each grid point on the pupil.
End of explanation
lutxp, lutyp, J = aperture2image(I1, inst, algo, algo.converge[:,-1], lutx, luty, projSamples, 'paraxial')
show_lutxyp = showProjection(lutxp, lutyp, inst.sensorFactor, projSamples, 1)
I1fit = Image(show_lutxyp, fieldXY, Image.INTRA)
I1fit.downResolution(oversample, I1.image0.shape[0], I1.image0.shape[1])
Explanation: We now trace the rays to the image plane. Lutxp and Lutyp are image coordinates for each (oversampled) ray. showProjection() makes the intensity image. Then, to down sample the image back to original resolution, we want to use the function downResolution() which is defined for the image class.
End of explanation
luty, lutx = np.mgrid[
        -(projSamples / 2 - 0.5):(projSamples / 2 + 0.5),
        -(projSamples / 2 - 0.5):(projSamples / 2 + 0.5)]
lutx = lutx / (projSamples / 2 / inst.sensorFactor)
luty = luty / (projSamples / 2 / inst.sensorFactor)
lutxp, lutyp, J = aperture2image(I2, inst, algo, algo.converge[:,-1], lutx, luty, projSamples, 'paraxial')
show_lutxyp = showProjection(lutxp, lutyp, inst.sensorFactor, projSamples, 1)
I2fit = Image(show_lutxyp, fieldXY, Image.EXTRA)
I2fit.downResolution(oversample, I2.image0.shape[0], I2.image0.shape[1])
#The atmosphere used here is just a random Gaussian smearing. We do not care much about the size at this point
from scipy.ndimage import gaussian_filter
atmSigma = .6/3600/180*3.14159*21.6/1.44e-5
I1fit.image[np.isnan(I1fit.image)]=0
a = gaussian_filter(I1fit.image, sigma=atmSigma)
fig, ax = plt.subplots(1,3, figsize=[15,4])
img = ax[0].imshow(I1fit.image, origin='lower')
ax[0].set_title('Forward prediction (no atm) Intra')
fig.colorbar(img, ax=ax[0])
img = ax[1].imshow(a, origin='lower')
ax[1].set_title('Forward prediction (w atm) Intra')
fig.colorbar(img, ax=ax[1])
img = ax[2].imshow(I1.image0, origin='lower')
ax[2].set_title('Real Image, Intra')
fig.colorbar(img, ax=ax[2])
I2fit.image[np.isnan(I2fit.image)]=0
b = gaussian_filter(I2fit.image, sigma=atmSigma)
fig, ax = plt.subplots(1,3, figsize=[15,4])
img = ax[0].imshow(I2fit.image, origin='lower')
ax[0].set_title('Forward prediction (no atm) Extra')
fig.colorbar(img, ax=ax[0])
img = ax[1].imshow(b, origin='lower')
ax[1].set_title('Forward prediction (w atm) Extra')
fig.colorbar(img, ax=ax[1])
img = ax[2].imshow(I2.image0, origin='lower')
ax[2].set_title('Real Image, Extra')
fig.colorbar(img, ax=ax[2])
Explanation: Now do the same thing for extra focal image
End of explanation
<END_TASK>