import sys import numpy as np import scipy.io from itertools import islice # from PSF.PSFInterp.PSFUtil import * # from PSFUtil import findNeighbors, findNeighbors_hoclist ###加载PSF信息### def LoadPSF(iccd, iwave, ipsf, psfPath, InputMaxPixelPos=True, PSFCentroidWgt=False, PixSizeInMicrons=2.5, PSFBinning=False, PSFConvGauss=False): """ load psf informations from psf matrix. Parameters: iccd (int): ccd number [1,30]. iwave(int): wave-index [1,4]. ipsf (int): psf number [1, 100]. psfPath (int): path to psf matrix InputMaxPixelPos(bool-optional): only True for 30*30 psf-matrix PSFCentroidWgt(bool-optional): True for using psfMat with libCentroid.lib Returns: psfInfo (dirctionary) """ if iccd not in np.linspace(1, 30, 30, dtype='int'): print('Error - iccd should be in [1, 30].') sys.exit() if iwave not in np.linspace(1, 4, 4, dtype='int'): print('Error - iwave should be in [1, 4].') sys.exit() if ipsf not in np.linspace(1, 900, 900, dtype='int'): print('Error - ipsf should be in [1, 900].') sys.exit() psfInfo = {} fpath = psfPath +'/' +'ccd{:}'.format(iccd) +'/' + 'wave_{:}'.format(iwave) #获取ipsf矩阵 ''' if not PSFCentroidWgt: ##读取PSF原数据 fpathMat = fpath +'/' +'5_psf_array' +'/' +'psf_{:}.mat'.format(ipsf) data = scipy.io.loadmat(fpathMat) psfInfo['psfMat'] = data['psf'] if PSFCentroidWgt: ##读取PSFCentroidWgt ffpath = psfPath +'_proc/' +'ccd{:}'.format(iccd) +'/' + 'wave_{:}'.format(iwave) ffpathMat = ffpath +'/' +'5_psf_array' +'/' +'psf_{:}_centroidWgt.mat'.format(ipsf) data = scipy.io.loadmat(ffpathMat) psfInfo['psfMat'] = data['psf'] ''' if not PSFCentroidWgt: ffpath = fpath ffpathMat = ffpath +'/5_psf_array/psf_{:}.mat'.format(ipsf) if PSFCentroidWgt: ffpath = psfPath +'_proc/ccd{:}/wave_{:}'.format(iccd, iwave) ffpathMat = ffpath +'/5_psf_array/psf_{:}_centroidWgt.mat'.format(ipsf) if PSFBinning and (not PSFConvGauss): ffpath = psfPath +'_bin256x256/ccd{:}/wave_{:}'.format(iccd, iwave) ffpathMat = ffpath +'/5_psf_array/psf_{:}_centroidWgt.mat'.format(ipsf) if (not PSFBinning) and PSFConvGauss: ffpath = psfPath +'_convGauss/ccd{:}/wave_{:}'.format(iccd, iwave) ffpathMat = ffpath +'/5_psf_array/psf_{:}_centroidWgt.mat'.format(ipsf) if PSFBinning and PSFConvGauss: ffpath = psfPath +'_convGauss/ccd{:}/wave_{:}'.format(iccd, iwave) ffpathMat = ffpath +'/5_psf_array/psf_{:}_centroidWgt_BC.mat'.format(ipsf) data = scipy.io.loadmat(ffpathMat) psfInfo['psfMat'] = data['psf'] psfInfo['pixsize'] = PixSizeInMicrons #microns if PSFBinning: psfInfo['pixsize'] = PixSizeInMicrons*2.0 #获取ipsf波长 fpathWave = fpath +'/' +'1_wavelength.txt' f = open(fpathWave, 'r') wavelength = np.float(f.readline()) f.close() psfInfo['wavelength'] = wavelength #获取ipsf位置 fpathCoordinate = fpath +'/' +'4_PSF_coordinate.txt' f = open(fpathCoordinate, 'r') header = f.readline() for line in islice(f, ipsf-1, ipsf): line = line.strip() columns = line.split() f.close() icol = 0 psfInfo['field_x'] = float(columns[icol]) #deg, 视场采样位置 icol+= 1 psfInfo['field_y'] = float(columns[icol]) #deg icol+= 1 psfInfo['centroid_x'] = float(columns[icol]) #mm, psf质心相对主光线的偏移量 icol+= 1 psfInfo['centroid_y'] = float(columns[icol]) #mm icol+= 1 if InputMaxPixelPos == True: psfInfo['max_x'] = float(columns[icol]) #mm, max pixel postion icol+= 1 psfInfo['max_y'] = float(columns[icol]) #mm icol+= 1 psfInfo['image_x'] = float(columns[icol]) #mm, 主光线位置 icol+= 1 psfInfo['image_y'] = float(columns[icol]) #mm if PSFCentroidWgt: psfInfo['centroid_x'] = data['cx'][0][0] #mm, psfCentroidWgt质心相对主光线的偏移量 psfInfo['centroid_y'] = data['cy'][0][0] #mm # if PSFBinning: # psfInfo['centroid_x'] = data['cx'][0][0] +0.0025/2 #binning引起的主光线漂移 # psfInfo['centroid_y'] = data['cy'][0][0] +0.0025/2 return psfInfo ###插值PSF图像-IDW方法### def psfMaker_IDW(px, py, PSFMat, cen_col, cen_row, IDWindex=2, OnlyNeighbors=True, hoc=None, hoclist=None, PSFCentroidWgt=False): """ psf interpolation by IDW Parameters: px, py (float, float): position of the target PSFMat (numpy.array): image cen_col, cen_row (numpy.array, numpy.array): potions of the psf centers IDWindex (int-optional): the power index of IDW OnlyNeighbors (bool-optional): only neighbors are used for psf interpolation Returns: psfMaker (numpy.array) """ minimum_psf_weight = 1e-8 ref_col = px ref_row = py ngy, ngx = PSFMat[0, :, :].shape npsf = PSFMat[:, :, :].shape[0] psfWeight = np.zeros([npsf]) from .PSFUtil import findNeighbors, findNeighbors_hoclist if OnlyNeighbors == True: if hoc is None: neigh = findNeighbors(px, py, cen_col, cen_row, dr=5., dn=4, OnlyDistance=False) if hoc is not None: neigh = findNeighbors_hoclist(cen_col, cen_row, tx=px,ty=py, dn=4, hoc=hoc, hoclist=hoclist) neighFlag = np.zeros(npsf) neighFlag[neigh] = 1 for ipsf in range(npsf): if OnlyNeighbors == True: if neighFlag[ipsf] != 1: continue dist = np.sqrt((ref_col - cen_col[ipsf])**2 + (ref_row - cen_row[ipsf])**2) if IDWindex == 1: psfWeight[ipsf] = dist if IDWindex == 2: psfWeight[ipsf] = dist**2 if IDWindex == 3: psfWeight[ipsf] = dist**3 if IDWindex == 4: psfWeight[ipsf] = dist**4 psfWeight[ipsf] = max(psfWeight[ipsf], minimum_psf_weight) psfWeight[ipsf] = 1./psfWeight[ipsf] psfWeight /= np.sum(psfWeight) psfMaker = np.zeros([ngy, ngx], dtype=np.float32) for ipsf in range(npsf): if OnlyNeighbors == True: if neighFlag[ipsf] != 1: continue iPSFMat = PSFMat[ipsf, :, :].copy() ipsfWeight = psfWeight[ipsf] psfMaker += iPSFMat * ipsfWeight psfMaker /= np.nansum(psfMaker) return psfMaker ###TEST### if __name__ == '__main__': iccd = 1 iwave= 1 ipsf = 1 psfPath = '/data/simudata/CSSOSDataProductsSims/data/csstPSFdata/CSSOS_psf_20210108/CSST_psf_ciomp_2p5um_cycle3_ccr90' psfInfo = LoadPSF(iccd, iwave, ipsf, psfPath, InputMaxPixelPos=True, PSFCentroidWgt=True) print('psfInfo.keys:', psfInfo.keys()) print(psfInfo)