update unittest module-PSF

tests/PSFInterpTest/test_Convolve.py

+'''
+test galsim.interpolatedImage & galsim.convolve
+'''
+
+import os
+import unittest
+
+import numpy as np
+import matplotlib.pyplot as plt
+import galsim
+import scipy.io
+from scipy import ndimage
+
+
+def psfEncircle(img, fraction=0.8, psfSampleSizeInMicrons=2.5, focalLengthInMeters=28, cenPix=None):
+    #imgMaxPix_x, imgMaxPix_y = findMaxPix(img)
+    y,x = ndimage.center_of_mass(img)  #y-rows, x-cols
+    imgMaxPix_x = x #int(x)
+    imgMaxPix_y = y #int(y)
+    if cenPix != None:
+        imgMaxPix_x = cenPix[0]
+        imgMaxPix_y = cenPix[1]
+
+    im1 = img.copy()
+    im1size = im1.shape
+
+    dis = np.zeros_like(img)
+    for irow in range(im1size[0]):
+        for icol in range(im1size[1]):
+            dx = icol - imgMaxPix_x
+            dy = irow - imgMaxPix_y
+            dis[irow, icol] = np.hypot(dx, dy)
+
+    nn = im1size[1]*im1size[0]
+    disX = dis.reshape(nn)
+    disXsortId = np.argsort(disX)
+
+    imgX = img.reshape(nn)
+    imgY = imgX[disXsortId]
+    psfFrac = np.cumsum(imgY)/np.sum(imgY)
+    ind = np.where(psfFrac > fraction)[0][0]
+
+    REE80 = np.rad2deg(dis[np.where(img == imgY[ind])]*psfSampleSizeInMicrons*1e-6/focalLengthInMeters)*3600
+    return REE80
+
+
+
+def check_galsimConvolve(path=None, plotImage=True):
+    #load psf data
+    data=scipy.io.loadmat(path)
+    imPSF = data['psf']
+    pixSize = np.rad2deg(5.*1e-6/28)*3600
+    imPSF = imPSF/np.sum(imPSF)
+
+    #psf -> galsimInterpolatedImage
+    img = galsim.ImageF(imPSF, scale=pixSize)
+    imgt= galsim.InterpolatedImage(img)
+
+    #imPSFt = imgt.drawImage(nx=256, ny=256, scale=pixSize, method='no_pixel')
+    imPSFt = imgt.drawImage(nx=256, ny=256, scale=pixSize)
+
+    ree80 = psfEncircle(imPSF, fraction=0.8, psfSampleSizeInMicrons=5.)
+    ree80_pix = ree80/(np.rad2deg((5.*1e-6/28))*3600)
+
+    sliceX = slice(128-int(np.round(ree80_pix[0])), 128+int(np.round(ree80_pix[0]))+1, 1)
+
+    #set a point sorce
+    src    = galsim.DeltaFunction(flux=1.0)
+
+    result = galsim.Convolve(src, imgt)
+
+    #drawImage with same pixSize
+    #tmp    = result.drawImage(nx=256, ny=256, scale=pixSize, method='no_pixel')
+    tmp    = result.drawImage(nx=256, ny=256, scale=pixSize)
+
+
+    if plotImage != True:
+        res = (imPSFt.array - imPSF)/imPSF
+        d0 = np.mean(res[sliceX, sliceX].flatten())
+
+        res = (tmp.array - imPSFt.array)/imPSFt.array
+        d1 = np.mean(res[sliceX, sliceX].flatten())
+        return d0, d1
+
+
+    #plot images
+    fig = plt.figure(figsize=(22, 5))
+    ax=plt.subplot(1,3,1)
+    plt.imshow(imPSF[128-10:128+10, 128-10:128+10])
+    plt.annotate("ORG", [0.1, 0.9],xycoords="axes fraction", fontsize=16, color="w")
+    plt.colorbar()
+
+    ax=plt.subplot(1,3,2)
+    plt.imshow(imPSFt.array[128-10:128+10, 128-10:128+10])
+    plt.annotate("InterpolatedImage", [0.1, 0.9],xycoords="axes fraction", fontsize=16, color="w")
+    plt.colorbar()
+
+    ax=plt.subplot(1,3,3)
+    plt.imshow(tmp.array[128-10:128+10, 128-10:128+10])
+    plt.annotate("ConvolvedImage", [0.1, 0.9],xycoords="axes fraction", fontsize=16, color="w")
+    plt.colorbar()
+    plt.savefig(OUTPUTPATH+'/fig_check_galsimConvolve_1.pdf')
+
+    fig = plt.figure(figsize=(13, 10))
+    ax=plt.subplot(2,2,1)
+    res = (imPSFt.array - imPSF)/imPSF
+    plt.imshow(res[128-10:128+10, 128-10:128+10])
+    plt.annotate("$\Delta_1$", [0.1, 0.9],xycoords="axes fraction", fontsize=16, color="w")
+    plt.colorbar()
+
+    ax=plt.subplot(2,2,2)
+    plt.hist(res[sliceX, sliceX].flatten(), alpha=0.75, bins=4)
+    #plt.annotate("$\Delta_1^{\\rm REE80}$", [0.1, 0.9],xycoords="axes fraction", fontsize=16, color="k")
+    plt.xlabel("$\Delta_1^{\\rm REE80}$", fontsize=16)
+    plt.ylabel("PDF", fontsize=16)
+
+    ax=plt.subplot(2,2,3)
+    res = (tmp.array - imPSFt.array)/imPSFt.array
+    plt.imshow(res[128-10:128+10, 128-10:128+10])
+    plt.annotate("$\Delta_2$", [0.1, 0.9],xycoords="axes fraction", fontsize=16, color="w")
+    plt.colorbar()
+
+    ax=plt.subplot(2,2,4)
+    plt.hist(res[sliceX, sliceX].flatten(), alpha=0.75, bins=4)
+    #plt.annotate("$\Delta_2^{\\rm REE80}$", [0.1, 0.9],xycoords="axes fraction", fontsize=16, color="k")
+    plt.xlabel("$\Delta_2^{\\rm REE80}$", fontsize=16)
+    plt.ylabel("PDF", fontsize=16)
+
+    plt.savefig(OUTPUTPATH+'/fig_check_galsimConvolve_2.pdf')
+
+
+def check_galsimConvolveALL(dataPath):
+    d0 = np.zeros(900)
+    d1 = np.zeros(900)
+    for ipsf in range(1,901,1):
+        print("ipsf={:}".format(ipsf), end="\r")
+        psfPath = dataPath+"/ccd1-w1/psf_{:}_centroidWgt_BC.mat".format(ipsf)
+        t0,t1=check_galsimConvolve(path = psfPath, plotImage=False)
+        d0[ipsf-1] = t0
+        d1[ipsf-1] = t1
+      
+    fig = plt.figure(figsize=(12,6))
+    ax = plt.subplot(1,2,1)
+    #plt.scatter(np.linspace(1,900,900), d0)
+    plt.hist(d0, bins=8, alpha=0.75)
+    plt.xlabel("mean($\Delta_1^{\\rm REE80}$)", fontsize=16)
+    plt.ylabel("PDF", fontsize=16)
+
+    ax = plt.subplot(1,2,2)
+    #plt.scatter(np.linspace(1,900,900), d1)
+    plt.hist(d1, bins=8, alpha=0.75)
+    plt.xlabel("mean($\Delta_2^{\\rm REE80}$)", fontsize=16)
+    plt.ylabel("PDF", fontsize=16)
+
+    plt.savefig(OUTPUTPATH+'/fig_check_galsimConvolveALL.pdf')
+ 
+
+class testConvolve(unittest.TestCase):
+    def __init__(self, methodName='runTest'):
+        super(testConvolve,self).__init__(methodName)
+        self.dataPath = os.path.join(os.getenv('UNIT_TEST_DATA_ROOT'), 'csst_fz_gc1')
+
+        global OUTPUTPATH
+        OUTPUTPATH = os.path.join(self.dataPath, 'outputs')
+
+    def test_galsimConvolve(self):
+        ipsf = 1
+        psfPath = self.dataPath+"/ccd1-w1/psf_{:}_centroidWgt_BC.mat".format(ipsf)
+        check_galsimConvolve(path = psfPath)
+
+    def test_galsimConvolveALL(self):
+        check_galsimConvolveALL(dataPath=self.dataPath)
+
+
+if __name__ == "__main__":
+    unittest.main()

tests/PSFInterpTest/PSFInterpModule_coverage.py → tests/PSFInterpTest/test_PSFInterpModule_coverage.py

@@ -35,7 +35,7 @@ def psfSecondMoments(psfMat, cenX, cenY, pixSize=1):
     pxs = 2.5 #microns
     apr = np.deg2rad(apr/3600.)*fl*1e6
     apr = apr/pxs
-    apr = np.int(np.ceil(apr))
+    apr = int(np.ceil(apr))
 
     I = psfMat
     ncol = I.shape[1]
@@ -90,7 +90,7 @@ def test_psfEll(iccd, iwave, psfMat):
         #print('ell======', ipsf, np.sqrt(e1**2 + e2**2))
     #######
     arr = [imx, imy, psf_e1, psf_e2, psf_sz]
-    np.save('data/psfEll{:}_{:}_{:}'.format(int(np.sqrt(npsf)),iccd, iwave), arr)
+    np.save(OUTPUTPATH+'/psfEll{:}_{:}_{:}'.format(int(np.sqrt(npsf)),iccd, iwave), arr)
 
 
 def test_psfEllPlot(OVERPLOT=False):
@@ -99,7 +99,7 @@ def test_psfEllPlot(OVERPLOT=False):
         prefix = 'psfEll30'
         iccd = 1
         iwave= 1
-        data = np.load('data/'+prefix+'_1_1.npy')
+        data = np.load(OUTPUTPATH+'/'+prefix+'_1_1.npy')
         imx= data[0]
         imy= data[1]
         psf_e1 = data[2]
@@ -136,7 +136,7 @@ def test_psfEllPlot(OVERPLOT=False):
 
         if OVERPLOT == True:
             prefix = 'psfEll20'
-            data = np.load('data/'+prefix+'_1_1.npy')
+            data = np.load(OUTPUTPATH+'/'+prefix+'_1_1.npy')
             imx= data[0]
             imy= data[1]
             psf_e1 = data[2]
@@ -156,7 +156,7 @@ def test_psfEllPlot(OVERPLOT=False):
         plt.gca().set_aspect(1)
         if OVERPLOT == True:
             prefix = 'psfEllOP'
-        plt.savefig('figs/'+prefix+'_iccd{:}.pdf'.format(iccd))
+        plt.savefig(OUTPUTPATH+'/'+prefix+'_iccd{:}.pdf'.format(iccd))
 
 
 def test_psfIDW(iccd, iwave, psfMatA, chip, psfMatB):
@@ -176,7 +176,7 @@ def test_psfIDW(iccd, iwave, psfMatA, chip, psfMatB):
         print('ipsf:', ipsf, end='\r', flush=True)
         tpos_img = pos_img(psfMatA.cen_col[iwave-1, ipsf-1], psfMatA.cen_row[iwave-1, ipsf-1])
         psfIDW = psfMatB.get_PSF(chip, tpos_img, bandpass, galsimGSObject=False, findNeighMode='treeFind')
-        np.save('figs/psfIDW/psfIDW_{:}_{:}_{:}'.format(iccd, iwave, ipsf), psfIDW)
+        np.save(OUTPUTPATH+'/psfIDW_{:}_{:}_{:}'.format(iccd, iwave, ipsf), psfIDW)
 
         cenX = 256
         cenY = 256
@@ -185,14 +185,14 @@ def test_psfIDW(iccd, iwave, psfMatA, chip, psfMatB):
         psf_e2[ipsf-1] = e2
         psf_sz[ipsf-1] = sz
     arr = [psf_e1, psf_e2, psf_sz]
-    np.save('data/psfEll20IDW_{:}_{:}'.format(iccd, iwave), arr)
+    np.save(OUTPUTPATH+'/psfEll20IDW_{:}_{:}'.format(iccd, iwave), arr)
 
 
 def test_psfResidualPlot(iccd, iwave, ipsf, psfMatA):
     psfMat_iwave = psfMatA.psfMat[iwave-1, :,:,:]
 
     psfMatORG = psfMat_iwave[ipsf-1, :, :]
-    psfMatIDW = np.load('figs/psfIDW/psfIDW_{:}_{:}_{:}.npy'.format(iccd, iwave, ipsf))
+    psfMatIDW = np.load(OUTPUTPATH+'/psfIDW_{:}_{:}_{:}.npy'.format(iccd, iwave, ipsf))
 
     npix = psfMatORG.shape[0]
     pixCutEdge= int(npix/2-15)
@@ -238,7 +238,7 @@ def test_psfResidualPlot(iccd, iwave, ipsf, psfMatA):
     cbar.ax.set_yticklabels(['$10^{-5}$','$10^{-4}$','$10^{-3}$','$10^{-2}$', '$10^{-1}$'])
 
     print(np.max((psfMatORG-psfMatIDW)))
-    plt.savefig('figs/psfResidual_iccd{:}.pdf'.format(iccd))
+    plt.savefig(OUTPUTPATH+'/psfResidual_iccd{:}.pdf'.format(iccd))
 
 
 def test_psfEllIDWPlot(OVERPLOT=False):
@@ -247,7 +247,7 @@ def test_psfEllIDWPlot(OVERPLOT=False):
         prefix = 'psfEll20'
         iccd = 1
         iwave= 1
-        data = np.load('data/'+prefix+'_1_1.npy')
+        data = np.load(OUTPUTPATH+'/'+prefix+'_1_1.npy')
         imx= data[0]
         imy= data[1]
         psf_e1 = data[2]
@@ -284,7 +284,7 @@ def test_psfEllIDWPlot(OVERPLOT=False):
 
         if OVERPLOT == True:
             prefix = 'psfEll20IDW'
-            data = np.load('data/'+prefix+'_1_1.npy')
+            data = np.load(OUTPUTPATH+'/'+prefix+'_1_1.npy')
             #imx= data[0]
             #imy= data[1]
             psf_e1 = data[0]
@@ -304,7 +304,7 @@ def test_psfEllIDWPlot(OVERPLOT=False):
         plt.gca().set_aspect(1)
         if OVERPLOT == True:
             prefix = 'psfEllOPIDW'
-        plt.savefig('figs/'+prefix+'_iccd{:}.pdf'.format(iccd))
+        plt.savefig(OUTPUTPATH+'/'+prefix+'_iccd{:}.pdf'.format(iccd))
 
 
 def test_psfdEllabsPlot(iccd):
@@ -313,7 +313,7 @@ def test_psfdEllabsPlot(iccd):
         prefix = 'psfEll20'
         #iccd = 1
         #iwave= 1
-        data = np.load('data/'+prefix+'_{:}_1.npy'.format(iccd))
+        data = np.load(OUTPUTPATH+'/'+prefix+'_{:}_1.npy'.format(iccd))
         imx= data[0]
         imy= data[1]
         psf_e1 = data[2]
@@ -329,7 +329,7 @@ def test_psfdEllabsPlot(iccd):
 
         ##############################
         prefix = 'psfEll20IDW'
-        data = np.load('data/'+prefix+'_{:}_1.npy'.format(iccd))
+        data = np.load(OUTPUTPATH+'/'+prefix+'_{:}_1.npy'.format(iccd))
         #imx= data[0]
         #imy= data[1]
         psf_e1 = data[0]
@@ -360,28 +360,33 @@ def test_psfdEllabsPlot(iccd):
         #plt.ylim([-0.0018, 0.0018])
         plt.xlabel('$\epsilon_{\\rm ORG}$')
         plt.ylabel('$\Delta$')
-        plt.savefig('figs/psfEllOPIDWPDF_{:}.pdf'.format(iccd))
+        plt.savefig(OUTPUTPATH+'/psfEllOPIDWPDF_{:}.pdf'.format(iccd))
 
         fig=plt.figure(figsize=(6, 6))
         plt.hist((szY-szX)/szX, bins=20, color='r', alpha=0.5)
         plt.xlabel('$(R_{\\rm IDW}-R_{\\rm ORG})/R_{\\rm ORG}$')
         plt.ylabel('PDF')
-        plt.savefig('figs/psfEllOPIDWPDF_dsz_{:}.pdf'.format(iccd))
+        plt.savefig(OUTPUTPATH+'/psfEllOPIDWPDF_dsz_{:}.pdf'.format(iccd))
 
 
 class PSFInterpModule_coverage(unittest.TestCase):
+    def __init__(self, methodName='runTest'):
+        super(PSFInterpModule_coverage,self).__init__(methodName)
+        self.dataPath = os.path.join(os.getenv('UNIT_TEST_DATA_ROOT'), 'csst_fz_gc1')
+        global OUTPUTPATH
+        OUTPUTPATH = os.path.join(self.dataPath, 'outputs')
+
     def test_psfEll_(self):
         iccd = 1
         iwave= 1
 
-        config_file = "/public/home/weichengliang/CSST_git/newVersion/CSST/config/config_C3.yaml"
+        config_file = os.path.join(self.dataPath, 'config_test.yaml') 
         chip = defineCCD(iccd, config_file)
         print(chip.chipID)
         print(chip.cen_pix_x, chip.cen_pix_y)
 
-        ipath = '/data/simudata/CSSOSDataProductsSims/data/csstPSFdata/psfCubeTest'
-        psfMatA = PSFInterp(chip, npsf=400, PSF_data_file=ipath, PSF_data_prefix="S20x20_")
-        psfMatB = PSFInterp(chip, npsf=900, PSF_data_file=ipath, PSF_data_prefix="S30x30_")
+        psfMatA = PSFInterp(chip, npsf=400, PSF_data_file=self.dataPath, PSF_data_prefix="S20x20_")
+        psfMatB = PSFInterp(chip, npsf=900, PSF_data_file=self.dataPath, PSF_data_prefix="S30x30_")
 
         test_psfEll(iccd, iwave, psfMatA)
         test_psfEll(iccd, iwave, psfMatB)

tests/PSFInterpTest/loadPSFSet.py → tests/PSFInterpTest/test_loadPSFSet.py

@@ -26,15 +26,14 @@ def defineCCD(iccd, config_file):
     #chip = Chip(chipID=iccd, ccdEffCurve_dir=path_dict["ccd_dir"], CRdata_dir=path_dict["CRdata_dir"], normalize_dir=path_dict["normalize_dir"], sls_dir=path_dict['sls_dir'], config=config)
     return chip
 
-def loadPSFSet(iccd):
-    config_file = "/public/home/weichengliang/CSST_git/newVersion/CSST/config/config_C3.yaml"
+def loadPSFSet(iccd, dataPath):
+    config_file = os.path.join(dataPath, 'config_test.yaml')
     chip = defineCCD(iccd, config_file)
     print(chip.chipID)
     print(chip.cen_pix_x, chip.cen_pix_y)
 
-    ipath = '/data/simudata/CSSOSDataProductsSims/data/csstPSFdata/psfCube'
-    psfMat= PSFInterp(chip, npsf=900, PSF_data_file=ipath, PSF_data_prefix="")
-    psfSet= psfMat._loadPSF(iccd, ipath, PSF_data_prefix="")
+    psfMat= PSFInterp(chip, npsf=900, PSF_data_file=dataPath, PSF_data_prefix="S30x30_")
+    psfSet= psfMat._loadPSF(iccd, dataPath, PSF_data_prefix="S30x30_")
 
     twave = 0 #[0...3]
     tpsf  = 0 #[0...899]
@@ -50,9 +49,13 @@ def loadPSFSet(iccd):
 
 
 class PSFInterpModule_coverage(unittest.TestCase):
+    def __init__(self, methodName='runTest'):
+        super(PSFInterpModule_coverage,self).__init__(methodName)
+        self.dataPath = os.path.join(os.getenv('UNIT_TEST_DATA_ROOT'), 'csst_fz_gc1')
+
     def test_psfEll_(self):
         iccd = 1  #[1...30]
-        psfSet = loadPSFSet(iccd)
+        psfSet = loadPSFSet(iccd, dataPath=self.dataPath)
 
 
 if __name__ == '__main__':