notebook.community

Edit and run

Tiago provided a pair of images from AuxTel.

Let's look at how those images work with our cwfs code

load the modules





In [1]:



    


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

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 )

It would probably be better if we do background subtraction before feeding the images to cwfs. Seems like when background is low (in our case, ~300 vs ~40000 for signal) we are still fine.
Another thing to note is that we don't wan the image stamps to be much larger than the donuts themselves, otherwise things becomes very slow. Below we will re-cut the image stamps.




In [2]:



    


fieldXY = [0,0]

I1 = Image(readFile('../tests/testImages/AuxTel2001/1579925613-16Pup_intra-0-1.fits'), fieldXY, Image.INTRA)
I2 = Image(readFile('../tests/testImages/AuxTel2001/1579925662-16Pup_extra-0-1.fits'), fieldXY, Image.EXTRA)

I1p = Image(readFile('../tests/testImages/AuxTel2001/1579925833-16Pup_intra-0-1.fits'), fieldXY, Image.INTRA)
I2p = Image(readFile('../tests/testImages/AuxTel2001/1579925882-16Pup_extra-0-1.fits'), fieldXY, Image.EXTRA)



    


















    






original intra image size = (880, 1164)
we cut it to (880, 880)
original extra image size = (880, 1164)
we cut it to (880, 880)
original intra image size = (880, 1164)
we cut it to (880, 880)
original extra image size = (880, 1164)
we cut it to (880, 880)

















In [3]:



    


plots.plotImage(I1.image,'intra')
plots.plotImage(I2.image,'extra')
plots.plotImage(I1p.image,'intra')
plots.plotImage(I2p.image,'extra')

Define the instrument. Input arguments: instrument name, size of image stamps





In [4]:



    


inst=Instrument('AuxTel',I1.sizeinPix)



    











In [5]:



    


print("Expected image diameter in pixels = %.0f"%(inst.offset/inst.fno/inst.pixelSize))



    


















    






Expected image diameter in pixels = 79

















In [6]:



    


I1.image = I1.image[300-80:300+80,400-80:400+80]
I1.sizeinPix = I1.image.shape[0]
I2.image = I2.image[300-80:300+80,400-80:400+80]
I2.sizeinPix = I2.image.shape[0]

I1p.image = I1p.image[350-80:350+80,400-80:400+80]
I1p.sizeinPix = I1p.image.shape[0]
I2p.image = I2p.image[350-80:350+80,400-80:400+80]
I2p.sizeinPix = I2p.image.shape[0]

inst=Instrument('AuxTel',I1.sizeinPix)



    











In [7]:



    


plots.plotImage(I1.image,'intra')
plots.plotImage(I2.image,'extra')
plots.plotImage(I1p.image,'intra')
plots.plotImage(I2p.image,'extra')

Define the algorithm being used. Input arguments: baseline algorithm, instrument, debug level





In [8]:



    


algo=Algorithm('exp',inst,0)
algop=Algorithm('exp',inst,0)

Run it





In [9]:



    


algo.runIt(inst,I1,I2,'paraxial')
algop.runIt(inst,I1p,I2p,'paraxial')



    


















    






WARNING: negative scale parameter,             image is within caustic, zcCol (in um)=


















In [10]:



    


print(algo.zer4UpNm)
print(algop.zer4UpNm)



    


















    






[-61.08548695  -6.26481943 137.74661812  40.99011616  37.98610521
 -38.14227693  -2.27884537  -7.45553831  -5.57344792  -3.53182548
   0.7031806   -2.72815556   0.28420168   1.30709954   4.03710218
  -0.29270525   3.39658136 -15.54449884  -7.20569604]
[-117.43286629   41.73833098   42.61086748   29.83620146   -6.23734934
  -27.33510745  -20.39833327  -16.31023311   -2.31753517    5.43923026
  -16.19910404   -4.70514527    7.64033705    0.51187898    6.12856685
    4.02476456   -2.04109311    1.77882406   -6.07392003]

plot the Zernikes Zn (n>=4)





In [11]:



    


plt.plot(range(4,23), algo.zer4UpNm,'b.-',label = '1')
plt.plot(range(4,23), algop.zer4UpNm,'r.-',label = '2')
plt.legend()
plt.grid()
#the Zernike do look a bit different, 
#but, 
#the images above (especially the intra focal images do look kind of different?)



    


















    
























In [12]:



    


plots.plotImage(I1.image0,'original intra',mask=I1.cMask)



    


















    
























In [13]:



    


plots.plotImage(I2.image0,'original extra', mask=I2.cMask)

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.





In [14]:



    


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])



    


















    
Out[14]:








<matplotlib.colorbar.Colorbar at 0x1c1dfdd1d0>










    
























In [15]:



    


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])



    


















    
Out[15]:








<matplotlib.colorbar.Colorbar at 0x1c1e053e50>

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.





In [16]:



    


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)

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.





In [17]:



    


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])

Now do the same thing for extra focal image





In [18]:



    


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])



    











In [19]:



    


#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])



    


















    
Out[19]:








<matplotlib.colorbar.Colorbar at 0x1c1e2dc450>










    
























In [20]:



    


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])



    


















    
Out[20]:








<matplotlib.colorbar.Colorbar at 0x1c1da9bd10>










    
























In [ ]:

Content source: bxin/cwfs

Similar notebooks:

notebook.community | gallery | about