Commit 883b7a77 authored by Fang Yuedong's avatar Fang Yuedong
Browse files

bug fixed

parent e82a5dcc
#include <math.h>
#define NRANSI
#include "nrutil.h"
void svdcmp(double **a, int m, int n, double *w, double **v)
{
int flag,i,its,j,jj,k,l,nm;
double anorm,c,f,g,h,s,scale,x,y,z,*rv1;
double pythag(double a, double b);
rv1=dvector(1,n);
g=scale=anorm=0.0;
for (i=1;i<=n;i++) {
l=i+1;
rv1[i]=scale*g;
g=s=scale=0.0;
if (i <= m) {
for (k=i;k<=m;k++) scale += fabs(a[k][i]);
if (scale) {
for (k=i;k<=m;k++) {
a[k][i] /= scale;
s += a[k][i]*a[k][i];
}
f=a[i][i];
g = -SIGN(sqrt(s),f);
h=f*g-s;
a[i][i]=f-g;
for (j=l;j<=n;j++) {
for (s=0.0,k=i;k<=m;k++) s += a[k][i]*a[k][j];
f=s/h;
for (k=i;k<=m;k++) a[k][j] += f*a[k][i];
}
for (k=i;k<=m;k++) a[k][i] *= scale;
}
}
w[i]=scale *g;
g=s=scale=0.0;
if (i <= m && i != n) {
for (k=l;k<=n;k++) scale += fabs(a[i][k]);
if (scale) {
for (k=l;k<=n;k++) {
a[i][k] /= scale;
s += a[i][k]*a[i][k];
}
f=a[i][l];
g = -SIGN(sqrt(s),f);
h=f*g-s;
a[i][l]=f-g;
for (k=l;k<=n;k++) rv1[k]=a[i][k]/h;
for (j=l;j<=m;j++) {
for (s=0.0,k=l;k<=n;k++) s += a[j][k]*a[i][k];
for (k=l;k<=n;k++) a[j][k] += s*rv1[k];
}
for (k=l;k<=n;k++) a[i][k] *= scale;
}
}
anorm=FMAX(anorm,(fabs(w[i])+fabs(rv1[i])));
}
for (i=n;i>=1;i--) {
if (i < n) {
if (g) {
for (j=l;j<=n;j++)
v[j][i]=(a[i][j]/a[i][l])/g;
for (j=l;j<=n;j++) {
for (s=0.0,k=l;k<=n;k++) s += a[i][k]*v[k][j];
for (k=l;k<=n;k++) v[k][j] += s*v[k][i];
}
}
for (j=l;j<=n;j++) v[i][j]=v[j][i]=0.0;
}
v[i][i]=1.0;
g=rv1[i];
l=i;
}
for (i=IMIN(m,n);i>=1;i--) {
l=i+1;
g=w[i];
for (j=l;j<=n;j++) a[i][j]=0.0;
if (g) {
g=1.0/g;
for (j=l;j<=n;j++) {
for (s=0.0,k=l;k<=m;k++) s += a[k][i]*a[k][j];
f=(s/a[i][i])*g;
for (k=i;k<=m;k++) a[k][j] += f*a[k][i];
}
for (j=i;j<=m;j++) a[j][i] *= g;
} else for (j=i;j<=m;j++) a[j][i]=0.0;
++a[i][i];
}
for (k=n;k>=1;k--) {
for (its=1;its<=30;its++) {
flag=1;
for (l=k;l>=1;l--) {
nm=l-1;
if ((double)(fabs(rv1[l])+anorm) == anorm) {
flag=0;
break;
}
if ((double)(fabs(w[nm])+anorm) == anorm) break;
}
if (flag) {
c=0.0;
s=1.0;
for (i=l;i<=k;i++) {
f=s*rv1[i];
rv1[i]=c*rv1[i];
if ((double)(fabs(f)+anorm) == anorm) break;
g=w[i];
h=pythag(f,g);
w[i]=h;
h=1.0/h;
c=g*h;
s = -f*h;
for (j=1;j<=m;j++) {
y=a[j][nm];
z=a[j][i];
a[j][nm]=y*c+z*s;
a[j][i]=z*c-y*s;
}
}
}
z=w[k];
if (l == k) {
if (z < 0.0) {
w[k] = -z;
for (j=1;j<=n;j++) v[j][k] = -v[j][k];
}
break;
}
if (its == 30) nrerror("no convergence in 30 svdcmp iterations");
x=w[l];
nm=k-1;
y=w[nm];
g=rv1[nm];
h=rv1[k];
f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
g=pythag(f,1.0);
f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x;
c=s=1.0;
for (j=l;j<=nm;j++) {
i=j+1;
g=rv1[i];
y=w[i];
h=s*g;
g=c*g;
z=pythag(f,h);
rv1[j]=z;
c=f/z;
s=h/z;
f=x*c+g*s;
g = g*c-x*s;
h=y*s;
y *= c;
for (jj=1;jj<=n;jj++) {
x=v[jj][j];
z=v[jj][i];
v[jj][j]=x*c+z*s;
v[jj][i]=z*c-x*s;
}
z=pythag(f,h);
w[j]=z;
if (z) {
z=1.0/z;
c=f*z;
s=h*z;
}
f=c*g+s*y;
x=c*y-s*g;
for (jj=1;jj<=m;jj++) {
y=a[jj][j];
z=a[jj][i];
a[jj][j]=y*c+z*s;
a[jj][i]=z*c-y*s;
}
}
rv1[l]=0.0;
rv1[k]=f;
w[k]=x;
}
}
free_dvector(rv1,1,n);
}
#undef NRANSI
/* (C) Copr. 1986-92 Numerical Recipes Software )1!. */
#OPTS += -D
CC = gcc
OPTIMIZE = -fPIC -shared -g -O3 #-Wall -wd981 #-wd1419 -wd810
#GSLI = -I/home/alex/opt/gsl/include
#GSLL = -L/home/alex/opt/gsl/lib -lgsl -lgslcblas
#FFTWI = -I/home/alex/opt/fftw/include
#FFTWL = -L/home/alex/opt/fftw/lib -lfftw3 -lfftw3f
#HDF5I = -I/home/alex/opt/hdf5/include
#HDF5L = -L/home/alex/opt/hdf5/lib -lhdf5_hl -lhdf5
#FITSI = -I/home/alex/opt/cfitsio/include
#FITSL = -L/home/alex/opt/cfitsio/lib -lcfitsio
#EXTRACFLAGS =
#EXTRACLIB =
CLINK=$(CC)
CFLAGS=$(OPTIMIZE) #$(EXTRACFLAGS) $(OPTS)
CLIB=-lm #$(EXTRACLIB)
OBJS = test.o
EXEC = libtest.so
all: $(EXEC)
$(EXEC): $(OBJS)
$(CLINK) $(CFLAGS) -o $@ $(OBJS) $(CLIB)
$(OBJS): Makefile
.PHONY : clean
clean:
rm -f *.o $(EXEC)
#include <stdio.h>
int sum(int a,int b){
return a + b;
}
import numpy as np
import ctypes
libPCA = ctypes.CDLL('../libPCA.so') # CDLL加载库
print('load libPCA')
npsf = 2
npix = 3
libPCA.psfPCA.argtypes = [ctypes.POINTER(ctypes.c_float), ctypes.c_int, ctypes.c_int, ctypes.POINTER(ctypes.c_double), ctypes.POINTER(ctypes.c_double)]
Nstar = npsf
Mp = npix*npix
NM = Nstar*Mp
NN = Nstar*Nstar
arr = (ctypes.c_float*NM)()
basef = (ctypes.c_double*NM)()
coeff = (ctypes.c_double*NN)()
#psf1 = np.random.random([npix, npix])
#psf2 = np.random.random([npix, npix])
psfT = np.random.random(Nstar*Mp)
arr[:] = psfT
libPCA.psfPCA(arr, Nstar, Mp, basef, coeff)
print('haha')
print(arr[:])
"""
CSST image simulation module (in python3): Point Spread Function (PSF)
author:: Wei Chengliang <chengliangwei@pmo.ac.cn>
"""
import sys
from itertools import islice
import numpy as np
import matplotlib.pyplot as plt
import scipy.io
from scipy.io import loadmat
#import xlrd
from scipy import ndimage
from scipy.interpolate import RectBivariateSpline
#from astropy.modeling.models import Ellipse2D
#from astropy.coordinates import Angle
#import matplotlib.patches as mpatches
import ctypes
import galsim
def setupPSFimg(iccd, iwave, psfPath="/data/simudata/CSSOSDataProductsSims/data/csstPSFdata/CSSOS_psf_ciomp"):
"""
psf model setup for csst-sim
Parameters:
iccd, iwave (int, int): psf model on iccd & iwave
psfPath (string, optional): path to psf matrix
Returns:
psf_model (psf_class): psf model
Methods:
psf_model.PSFinplace(self, px, py, interpScheme=1): psf interpolation
psf_model.PSFspin(self, psf, sigSpin, sigGauss, dx, dy): psf rotation (from Yudong)
"""
psf_model = PSFimg(iccd, iwave, psfPath)
return psf_model
##################################################
# A. psf matrix loading & checking #
##################################################
def psfPixelLayout(nrows, ncols, cenPosRow, cenPosCol, pixSizeInMicrons=5.0):
"""
convert psf pixels to physical position
Parameters:
nrows, ncols (int, int): psf sampling with [nrows, ncols].
cenPosRow, cenPosCol (float, float): A physical position of the chief ray for a given psf.
pixSizeInMicrons (float-optional): The pixel size in microns from the psf sampling.
Returns:
psfPixelPos (numpy.array-float): [posx, posy] in mm for [irow, icol]
Notes:
1. show positions on ccd, but not position on image only [+/- dy]
"""
psfPixelPos = np.zeros([2, nrows, ncols])
if nrows % 2 != 0:
sys.exit()
if ncols % 2 != 0:
sys.exit()
cenPix_row = nrows/2 + 1 #中心主光线对应pixle [由长光定义]
cenPix_col = ncols/2 + 1
for irow in range(nrows):
for icol in range(ncols):
delta_row = ((irow + 1) - cenPix_row)*pixSizeInMicrons*1e-3
delta_col = ((icol + 1) - cenPix_col)*pixSizeInMicrons*1e-3
psfPixelPos[0, irow, icol] = cenPosCol + delta_col
psfPixelPos[1, irow, icol] = cenPosRow - delta_row #note-1
return psfPixelPos
def imSigRange(img, fraction=0.80):
"""
extract the image within x-percent (DISCARD)
Parameters:
img (numpy.array-float): image
fraction (float-optional): a percentage
Returns:
im1 (numpy.array-float): image
"""
im1 = img.copy()
im1size = im1.shape
im2 = np.sort(im1.reshape(im1size[0]*im1size[1]))
im2 = im2[::-1]
im3 = np.cumsum(im2)/np.sum(im2)
loc = np.where(im3 > fraction)
#print(im3[loc[0][0]], im2[loc[0][0]])
im1[np.where(im1 <= im2[loc[0][0]])]=0
return im1
def imPlotInRange(img):
"""
plot image within a selected range
Parameters:
img (numpy.array-float): image
Returns:
"""
im1 = img.copy()
im1size = im1.shape
X,Y = np.meshgrid(range(im1size[1]),range(im1size[0]))
Z = im1
resolution = 25
f = lambda x,y: Z[int(y),int(x) ]
g = np.vectorize(f)
x = np.linspace(0,Z.shape[1], Z.shape[1]*resolution)
y = np.linspace(0,Z.shape[0], Z.shape[0]*resolution)
X2, Y2= np.meshgrid(x[:-1],y[:-1])
Z2 = g(X2,Y2)
#plt.pcolormesh(X,Y, Z)
#plt.imshow(img, origin='lower')
plt.contour(X2-0.5,Y2-0.5,Z2, [0.], colors='w', linestyles='--', linewidths=[1])
return
def findMaxPix(img):
"""
get the pixel position of the maximum-value
Parameters:
img (numpy.array-float): image
Returns:
imgMaxPix_x, imgMaxPix_y (int, int): pixel position in columns & rows
"""
maxIndx = np.argmax(img)
maxIndx = np.unravel_index(maxIndx, np.array(img).shape)
imgMaxPix_x = maxIndx[1]
imgMaxPix_y = maxIndx[0]
return imgMaxPix_x, imgMaxPix_y
def psfTailor(img, apSizeInArcsec=0.5, psfSampleSizeInMicrons=5, focalLengthInMeters=28):
"""
psf tailor within a given aperture size
Parameters:
img (numpy.array-float): image
apSizeInArcsec (float-optional): aperture size in arcseconds.
psfSampleSizeInMicrons (float-optional): psf pixel size in microns.
focalLengthInMeters (float-optional): csst focal length im meters.
Returns:
imgT (numpy.array-float): image
"""
imgMaxPix_x, imgMaxPix_y = findMaxPix(img)
apSizeInMicrons = np.deg2rad(apSizeInArcsec/3600.)*focalLengthInMeters*1e6
apSizeInPix = apSizeInMicrons/psfSampleSizeInMicrons
apSizeInPix = np.int(np.ceil(apSizeInPix))
imgT = np.zeros_like(img)
imgT[imgMaxPix_y-apSizeInPix:imgMaxPix_y+apSizeInPix+1,
imgMaxPix_x-apSizeInPix:imgMaxPix_x+apSizeInPix+1] = \
img[imgMaxPix_y-apSizeInPix:imgMaxPix_y+apSizeInPix+1,
imgMaxPix_x-apSizeInPix:imgMaxPix_x+apSizeInPix+1]
return imgT
def psfEncircle(img, fraction=0.8, psfSampleSizeInMicrons=5, focalLengthInMeters=28):
"""
psf tailor within a given percentage.
Parameters:
img (numpy.array-float): image
fraction (float-optional): a percentage for psf tailor.
psfSampleSizeInMicrons (float-optional): psf pixel size in microns.
focalLengthInMeters (float-optional): csst focal length im meters.
Returns:
img*wgt (numpy.array-float): image
REE80 (float): radius of REE80 in arcseconds.
"""
imgMaxPix_x, imgMaxPix_y = findMaxPix(img)
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]
wgt = np.ones_like(dis)
wgt[np.where(dis > dis[np.where(img == imgY[ind])])] = 0
REE80 = np.rad2deg(dis[np.where(img == imgY[ind])]*psfSampleSizeInMicrons*1e-6/focalLengthInMeters)*3600
return img*wgt, REE80
def psfSecondMoments(psfMat, cenX, cenY, pixSize=1):
"""
estimate the psf ellipticity by the second moment of surface brightness
Parameters:
psfMat (numpy.array-float): image
cenX, cenY (float, float): pixel position of the psf center
pixSize (float-optional): pixel size
Returns:
sz (float): psf size
e1, e2 (float, float): psf ellipticity
"""
I = psfMat
ncol = I.shape[1]
nrow = I.shape[0]
w = 0.0
w11 = 0.0
w12 = 0.0
w22 = 0.0
for icol in range(ncol):
for jrow in range(nrow):
x = icol*pixSize - cenX
y = jrow*pixSize - cenY
w += I[jrow, icol]
w11 += x*x*I[jrow, icol]
w12 += x*y*I[jrow, icol]
w22 += y*y*I[jrow, icol]
w11 /= w
w12 /= w
w22 /= w
sz = w11 + w22
e1 = (w11 - w22)/sz
e2 = 2.0*w12/sz
return sz, e1, e2
def LoadPSF(iccd, iwave, ipsf, psfPath, psfSampleSize=5, CalcPSFsize=True, CalcPSFcenter=True, SigRange=False, TailorScheme=1, InputMaxPixelPos=False):
'''加载psf信息'''
"""
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
psfSampleSize (float-optional): psf size in microns.
CalcPSFsize (bool-optional): whether calculate psf size & ellipticity. Default: True
CalcPSFcenter (bool-optional): whether calculate psf center. Default: True
SigRange (bool-optional): whether use psf tailor. Default: False
TailorScheme (int-optional): which method for psf tailor. Default: 1
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矩阵
fpathMat = fpath +'/' +'5_psf_array' +'/' +'psf_{:}.mat'.format(ipsf)
data = scipy.io.loadmat(fpathMat)
psfInfo['psfMat'] = data['psf']
#获取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
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
#nrows = 180 #psf采样大小, in pixels
#ncols = 180
nrows, ncols = psfInfo['psfMat'].shape
psfPos = psfPixelLayout(nrows, ncols, psfInfo['image_y'], psfInfo['image_x'], pixSizeInMicrons=5.0)
imgMaxPix_x, imgMaxPix_y = findMaxPix(psfInfo['psfMat'])
psfInfo['imgMaxPosx_ccd'] = psfPos[0, imgMaxPix_y, imgMaxPix_x] #cx, psf最大值位置, in mm
psfInfo['imgMaxPosy_ccd'] = psfPos[1, imgMaxPix_y, imgMaxPix_x] #cy
#计算psf size & ellipticity
if CalcPSFsize is True:
psfMat = psfInfo['psfMat'].copy()
cenX, cenY, sz, e1, e2, REE80 = psfSizeCalculator(psfMat, psfSampleSize=psfSampleSize, CalcPSFcenter=CalcPSFcenter, SigRange=SigRange, TailorScheme=TailorScheme)
psfInfo['psfCenX_img'] = cenX #in local pixels, psf质心位置, in pixels
psfInfo['psfCenY_img'] = cenY #in local pixels
psfInfo['psfSize'] = sz
psfInfo['psf_e1'] = e1
psfInfo['psf_e2'] = e2
psfInfo['REE80'] = REE80
return psfInfo
def psfSizeCalculator(psfMat, psfSampleSize=5, CalcPSFcenter=True, SigRange=False, TailorScheme=1):
"""
calculate psf size & ellipticity
Parameters:
psfMat (numpy.array): image
psfSampleSize (float-optional): psf size in microns.
CalcPSFcenter (bool-optional): whether calculate psf center. Default: True
SigRange (bool-optional): whether use psf tailor. Default: False
TailorScheme (int-optional): which method for psf tailor. Default: 1
Returns:
cenX, cenY (float, float): the pixel position of the mass center
sz (float): psf size
e1, e2 (float, float): psf ellipticity
REE80 (float): radius of REE80 in arcseconds
"""
psfSampleSize = psfSampleSize*1e-3 #mm
REE80 = -1.0 ##encircling 80% energy
if SigRange is True:
if TailorScheme == 1:
psfMat = imSigRange(psfMat, fraction=0.80)
psfInfo['psfMat'] = psfMat #set on/off
if TailorScheme == 2:
img = psfTailor(psfMat, apSizeInArcsec=0.5)
imgX, REE80 = psfEncircle(psfMat)
psfMat = img
REE80 = REE80[0]
if CalcPSFcenter is True:
img = psfMat/np.sum(psfMat)
y,x = ndimage.center_of_mass(img) #y-rows, x-cols
cenX = x
cenY = y
if CalcPSFcenter is False:
cenPix_X = psfMat.shape[1]/2 #90
cenPix_Y = psfMat.shape[0]/2 #90
cenX = cenPix_X + psfInfo['centroid_x']/psfSampleSize
cenY = cenPix_Y - psfInfo['centroid_y']/psfSampleSize
pixSize = 1
sz, e1, e2 = psfSecondMoments(psfMat, cenX, cenY, pixSize=pixSize)
return cenX, cenY, sz, e1, e2, REE80
def psfStack(*psfMat):
"""
stacked image from the different psfs
Parameters:
*psfMat (numpy.array): the different psfs for stacking
Returns:
img (numpy.array): image
"""
nn = len(psfMat)
img = np.zeros_like(psfMat[0])
for ii in range(nn):
img += psfMat[ii]/np.sum(psfMat[ii])
img /= np.sum(img)
return img
##################################################
# B. psf interpolation #
##################################################
def img2fits(img, fitsName=None):
"""
saving image to fits file
Parameters:
img (numpy.array): image
fitsName (string): path+filename of fits
Returns:
"""
from astropy.io import fits
grey = fits.PrimaryHDU(img)
greyHDU = fits.HDUList([grey])
if fitsName != None:
greyHDU.writeto(fitsName)
def psfMatLoad(iccd, iwave, psfPath, psfSampleSize=5, CalcPSFsize=False, CalcPSFcenter=True):
"""
load psf for interpolation
Parameters:
iccd, iwave, psfPath: # of ccd/wave and path for psfs
CalcPSFsize (bool-optional): whether calculate psf size & ellipticity. Default: False
CalcPSFcenter (bool-optional): whether calculate psf center. Default: True
Returns:
PSFMat (numpy.array): images
cen_col, cen_row (numpy.array, numpy.array): position of psf center in the view field
"""
psfSet = []
for ipsf in range(1, 101):
psfInfo = LoadPSF(iccd, iwave, ipsf, psfPath, CalcPSFsize=CalcPSFsize, CalcPSFcenter=CalcPSFcenter, InputMaxPixelPos=False)
psfSet.append(psfInfo)
npsf = len(psfSet)
ngy, ngx = psfSet[0]['psfMat'].shape
PSFMat = np.zeros([npsf, ngy, ngx])
cen_col= np.zeros(npsf)
cen_row= np.zeros(npsf)
FieldPos = False
for ipsf in range(npsf):
PSFMat[ipsf, :, :] = psfSet[ipsf]['psfMat']
if FieldPos == True:
cen_col[ipsf] = psfSet[ipsf]['field_x'] #cx
cen_row[ipsf] = psfSet[ipsf]['field_y'] #cy
if FieldPos == False:
cen_col[ipsf] = psfSet[ipsf]['imgMaxPosx_ccd'] #cx
cen_row[ipsf] = psfSet[ipsf]['imgMaxPosy_ccd'] #cy
return PSFMat, cen_col, cen_row
def findNeighbors(tx, ty, px, py, dr=0.1, dn=1, OnlyDistance=True):
"""
find nearest meighbors by 2D-KDTree
Parameters:
tx, ty (float, float): a given position
px, py (numpy.array, numpy.array): position data for tree
dr (float-optional): distance
dn (int-optional): nearest-N
OnlyDistance (bool-optional): only use distance to find neighbors. Default: True
Returns:
dataq (numpy.array): index
"""
import scipy.spatial as spatial
datax = px
datay = py
tree = spatial.KDTree(list(zip(datax.ravel(), datay.ravel())))
dataq=[]
rr = dr
if OnlyDistance == True:
dataq = tree.query_ball_point([tx, ty], rr)
if OnlyDistance == False:
while len(dataq) < dn:
dataq = tree.query_ball_point([tx, ty], rr)
rr += dr
dd = np.hypot(datax[dataq]-tx, datay[dataq]-ty)
ddSortindx = np.argsort(dd)
dataq = np.array(dataq)[ddSortindx[0:dn]]
return dataq
def psfCentering(img, apSizeInArcsec=4., psfSampleSizeInMicrons=5, focalLengthInMeters=28, CenteringMode=1):
"""
centering psf within an aperture
Parameters:
img (numpy.array): image
apSizeInArcsec (float-optional): aperture size in arcseconds.
psfSampleSizeInMicrons (float-optional): psf pixel size in microns.
focalLengthInMeters (float-optional): csst focal length im meters.
CenteringMode (int-optional): how to center psf images
Returns:
imgT (numpy.array)
"""
if CenteringMode == 1:
imgMaxPix_x, imgMaxPix_y = findMaxPix(img)
if CenteringMode == 2:
y,x = ndimage.center_of_mass(img) #y-rows, x-cols
imgMaxPix_x = int(x)
imgMaxPix_y = int(y)
apSizeInMicrons = np.deg2rad(apSizeInArcsec/3600.)*focalLengthInMeters*1e6
apSizeInPix = apSizeInMicrons/psfSampleSizeInMicrons
apSizeInPix = np.int(np.ceil(apSizeInPix))
imgT = np.zeros_like(img)
ngy, ngx = img.shape
cy = int(ngy/2)
cx = int(ngx/2)
imgT[cy-apSizeInPix:cy+apSizeInPix+1,
cx-apSizeInPix:cx+apSizeInPix+1] = \
img[imgMaxPix_y-apSizeInPix:imgMaxPix_y+apSizeInPix+1,
imgMaxPix_x-apSizeInPix:imgMaxPix_x+apSizeInPix+1]
return imgT
def psfMaker_IDW(px, py, PSFMat, cen_col, cen_row, IDWindex=2, OnlyNeighbors=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])
if OnlyNeighbors == True:
neigh = findNeighbors(px, py, cen_col, cen_row, dr=5., dn=9, OnlyDistance=False)
neighFlag = np.zeros(npsf)
neighFlag[neigh] = 1
# print("neigh:", neigh)
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='float64')
for ipsf in range(npsf):
if OnlyNeighbors == True:
if neighFlag[ipsf] != 1:
continue
iPSFMat = PSFMat[ipsf, :, :].copy()
iPSFMat = psfCentering(iPSFMat, CenteringMode=1)
ipsfWeight = psfWeight[ipsf]
psfMaker += iPSFMat * ipsfWeight
psfMaker /= np.nansum(psfMaker)
return psfMaker
def psfMaker_PCA(px, py, PSFMat, cen_col, cen_row, OnlyNeighbors=False, libPCApath='libPCA/libPCA.so'):
"""
psf interpolation by PCA
Parameters:
Returns:
"""
ref_col = px
ref_row = py
ngy, ngx = PSFMat[0, :, :].shape
npsf = PSFMat[:, :, :].shape[0]
neigh = findNeighbors(px, py, cen_col, cen_row, dr=0.3, dn=5, OnlyDistance=False)
npsfX = len(neigh)
psfMatX = np.zeros([npsfX, ngy, ngx])
cen_colX= np.zeros(npsfX)
cen_rowX= np.zeros(npsfX)
for ipsf in range(npsfX):
psfMatX[ipsf, :, :] = PSFMat[neigh[ipsf], :, :]
cen_colX[ipsf] = cen_col[neigh[ipsf]]
cen_rowX[ipsf] = cen_row[neigh[ipsf]]
psfMaker = np.zeros((ngy, ngx), dtype='float64')
if OnlyNeighbors == True:
PCAbasef, PCAcoeff = psfPCA_generator(psfMatX, npsfX, ngx, libPCApath)
nPCA = npsfX
for iPCA in range(nPCA):
coeffX = fitPoly(ref_col, ref_row, cen_colX, cen_rowX, PCAcoeff[:, iPCA], order=2)
psfMaker += coeffX*PCAbasef[iPCA, :, :]
return psfMaker
def psfPCA_generator(psfMat, npsf, npix, libPCApath):
"""
generate PCs from psfMat
Parameters:
Returns:
"""
libPCA = ctypes.CDLL(libPCApath) # CDLL加载库
libPCA.psfPCA.argtypes = [ctypes.POINTER(ctypes.c_float), ctypes.c_int, ctypes.c_int, ctypes.POINTER(ctypes.c_double), ctypes.POINTER(ctypes.c_double)]
Nstar = npsf
Mp = npix*npix
NM = Nstar*Mp
NN = Nstar*Nstar
arr = (ctypes.c_float*NM)()
basef = (ctypes.c_double*NM)()
coeff = (ctypes.c_double*NN)()
psfT = np.zeros(NM)
for ipsf in range(npsf):
lp = 0 + ipsf*Mp
up = Mp+ ipsf*Mp
ipsfMat = psfMat[ipsf, :, :]
psfT[lp:up] = ipsfMat.reshape(Mp)
arr[:] = psfT
libPCA.psfPCA(arr, Nstar, Mp, basef, coeff)
PCAbasef = np.zeros([npsf, npix, npix])
PCAcoeff = np.zeros([npsf, npsf])
for ipsf in range(npsf):
lp = 0 + ipsf*Mp
up = Mp+ ipsf*Mp
PCAbasef[ipsf, :, :] = np.array(basef[lp:up]).reshape(npix, npix)
lp = 0 + ipsf*npsf
up = npsf+ ipsf*npsf
PCAcoeff[ipsf, :] = np.array(coeff[lp:up])
return PCAbasef, PCAcoeff
def fitPoly(px, py, datax, datay, dataz, order = 2):
if order == 1:
# best-fit linear plane
A = np.c_[datax, datay, np.ones(datax.shape[0])]
C,_,_,_ = scipy.linalg.lstsq(A, dataz) # coefficients
pz = C[0]*px + C[1]*py + C[2]
elif order == 2:
# best-fit quadratic curve
A = np.c_[np.ones(datax.shape[0]), np.c_[datax, datay], np.prod(np.c_[datax, datay], axis=1), np.c_[datax, datay]**2]
C,_,_,_ = scipy.linalg.lstsq(A, dataz)
pz = np.dot(np.c_[1, px, py, px*py, px**2, py**2], C)
"""
elif order == 3:
# best-fit cubic curve
A = np.c_[np.ones(datax.shape[0]), np.c_[datax, datay], np.prod(np.c_[datax, datay], axis=1), np.c_[datax, datay]**2, np.c_[datax, datay]**3]
C,_,_,_ = scipy.linalg.lstsq(A, dataz)
pz = np.dot(np.c_[1, px, py, px*py, px**2, py**2, px**3, py**3], C)
"""
return pz
"""
############################
### not used temporarily ###
############################
def psfSplineMake(px, py, PSFMat, cen_col, cen_row, OnlyNeighbors=False):
minimum_psf_weight = 1e-8
ref_col = px
ref_row = py
cdelt1p = 1
cdelt2p = 1
ngy, ngx = PSFMat[0, :, :].shape
psfx = np.linspace(0, ngx-1, ngx)
psfy = np.linspace(0, ngy-1, ngy)
npsf = PSFMat[:, :, :].shape[0]
psfWeight = np.zeros([npsf])
for ipsf in range(npsf):
psfWeight[ipsf] = np.sqrt((ref_col - cen_col[ipsf])**2 + (ref_row - cen_row[ipsf])**2)
psfWeight[ipsf] = max(psfWeight[ipsf], minimum_psf_weight)
psfWeight[ipsf] = 1./psfWeight[ipsf]
psfWeight /= np.sum(psfWeight)
psf = np.zeros((ngy, ngx), dtype='float64')
for ipsf in range(npsf):
iPSFMat = PSFMat[ipsf, :, :]
ipsfWeight = psfWeight[ipsf]
psf += iPSFMat * ipsfWeight
psf /= (np.nansum(psf) * cdelt1p * cdelt2p)
psfSpline = RectBivariateSpline(psfy, psfx, psf)
return psf, psfSpline
def psfToImage(psfSpline, cutoff_radius=180):
ng = 180
img = np.zeros([ng, ng], dtype='float64')
for i in range(ng):
for j in range(ng):
star_row = 5
star_column = 5
if np.sqrt((j-star_column)**2 + (i-star_row)**2) <= cutoff_radius:
star_flux = 8
column_cen = j #j - star_column
row_cen = i #i - star_row
img[i,j] += star_flux * psfSpline.integral(row_cen-0.5, row_cen+0.5, column_cen-0.5, column_cen+0.5)
return img
"""
##################################################
# C. csstPSF class #
##################################################
class PSFimg(object):
def __init__(self, iccd, iwave, psfPath):
self.iccd = iccd
self.iwave= iwave
self.psfPath = psfPath
#loading psfSet >>>
"""
psfSet = []
for ipsf in range(1, 901):
psfInfo = LoadPSF(iccd, iwave, ipsf, psfPath, CalcPSFsize=True, CalcPSFcenter=True, SigRange=False)
psfSet.append(psfInfo)
self.psfSet = psfSet
"""
a, b, c = psfMatLoad(iccd, iwave, psfPath)
self.psfMat = a
self.cenPosx= b
self.cenPosy= c
def PSFinplace(self, px, py, interpScheme=1):
if interpScheme == 1:
idwIndx = 2
psf = psfMaker_IDW(px, py, self.psfMat, self.cenPosx, self.cenPosy, IDWindex=idwIndx, OnlyNeighbors=True)
if interpScheme ==2:
libPCA = "/Users/chengliangwei/Desktop/csstPSF/libPCA/libPCA.so"
psf = psfMaker_PCA(px, py, self.psfMat, self.cenPosx, self.cenPosy, OnlyNeighbors=True, libPCApath=libPCA)
img = galsim.ImageF(psf, scale=0.074/2)
xpsf = galsim.InterpolatedImage(img)
return xpsf
"""
def gcPlot(self, psf,pscale=0.074,figout="GC.png"):
size = np.size(psf,axis=0)
cxy = 0.5*(size-1)
width = 0.5*size
# log scale
radius = np.arange(np.log10(0.2),np.log10(width),0.01)
radius = 10.0**radius
nr = len(radius)
gc = []
for i in range(nr):
iflux, iferr, xflag = sep.sum_circle(psf,cxy,cxy,radius[i],subpix=0)
gc += [iflux.tolist()]
# Estimate the radius for a given flux ratio
fratio = 0.8
mid = [i for i in range(nr) if gc[i]<=fratio and gc[i+1]>fratio][0]
r0, r1 = radius[mid], radius[mid+1]
gc0, gc1 = gc[mid], gc[mid+1]
r5 = (fratio-gc0)/(gc1-gc0)*(r1-r0) + r0
hlf = r5*pscale
# plot
pfit = interp1d(radius, gc, kind='cubic')
fig = pl.figure(figsize=(5.5,4.0))
ax = fig.add_axes([0.16,0.15,0.80,0.81])
ax.plot(radius*pscale, pfit(radius), "k", linewidth=2.0)
ax.plot(radius*pscale, gc, "*r", markersize=5.0,mec="r",alpha=0.3)
ax.plot([hlf,hlf],[0,fratio],"k",linewidth=2.5)
ax.plot([0,hlf],[fratio,fratio],"k",linewidth=2.5)
ax.text(radius[10]*pscale,0.6,"$r_{%.1f}$=%.2f\""%(fratio,hlf))
ax.set_xlabel("Radius (arcsec)",fontsize=15)
ax.set_ylabel("Growth of Curve",fontsize=15)
ax.set_xscale("log")
ax.set_xlim(radius[0]*pscale,radius[-1]*pscale)
ax.set_ylim(0.0,1.0)
for tick in ax.xaxis.get_major_ticks(): tick.label.set_fontsize(15)
for tick in ax.yaxis.get_major_ticks(): tick.label.set_fontsize(15)
pl.savefig(figout)
pl.clf()
pl.close()
return
"""
def PSFspin(self, psf, sigSpin, sigGauss, dx, dy):
"""
The PSF profile at a given image position relative to the axis center
Parameters:
theta : spin angles in a given exposure in unit of [arcsecond]
dx, dy: relative position to the axis center in unit of [pixels]
Return:
Spinned PSF: g1, g2 and axis ratio 'a/b'
"""
a2Rad = np.pi/(60.0*60.0*180.0)
ff = sigGauss * 0.107 * (1000.0/10.0) # in unit of [pixels]
rc = np.sqrt(dx*dx + dy*dy)
cpix = rc*(sigSpin*a2Rad)
beta = (np.arctan2(dy,dx) + np.pi/2)
ell = cpix**2/(2.0*ff**2+cpix**2)
#ell *= 10.0
qr = np.sqrt((1.0+ell)/(1.0-ell))
#psfShape = galsim.Shear(e=ell, beta=beta)
#g1, g2 = psfShape.g1, psfShape.g2
#qr = np.sqrt((1.0+ell)/(1.0-ell))
#return ell, beta, qr
PSFshear = galsim.Shear(e=ell, beta=beta*galsim.radians)
return psf.shear(PSFshear), PSFshear(base)
##################################################
# D. TEST #
##################################################
def psfMaker_IDW_test(tpsf, px, py, PSFMat, cen_col, cen_row, IDWindex=2, OnlyNeighbors=False):
minimum_psf_weight = 1e-8
ref_col = px
ref_row = py
ngy, ngx = PSFMat[0, :, :].shape
npsf = PSFMat[:, :, :].shape[0]
psfWeight = np.zeros([npsf])
if OnlyNeighbors == True:
neigh = findNeighbors(px, py, cen_col, cen_row, dr=0.1, dn=9, OnlyDistance=False)
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='float64')
for ipsf in range(npsf):
"""
if OnlyNeighbors == True:
iy, ix = np.unravel_index(ipsf, (10,10))
ty, tx = np.unravel_index(tpsf, (10,10))
if np.abs(iy - ty) > 1 or np.abs(ix - tx) > 1:
continue
"""
if OnlyNeighbors == True:
if neighFlag[ipsf] != 1:
continue
if ipsf == tpsf:
continue
iPSFMat = PSFMat[ipsf, :, :].copy()
iPSFMat = psfCentering(iPSFMat, CenteringMode=1)
ipsfWeight = psfWeight[ipsf]
psfMaker += iPSFMat * ipsfWeight
psfMaker /= np.nansum(psfMaker)
return psfMaker
def psfMaker_PCA_test(tpsf, px, py, PSFMat, cen_col, cen_row, OnlyNeighbors=False, libPCApath='libPCA/libPCA.so'):
"""
psf interpolation by PCA
Parameters:
Returns:
"""
ref_col = px
ref_row = py
ngy, ngx = PSFMat[0, :, :].shape
npsf = PSFMat[:, :, :].shape[0]
neigh = findNeighbors(px, py, cen_col, cen_row, dr=0.3, dn=9, OnlyDistance=False)
npsfX = len(neigh)
#去掉tpsf,neigh中排在第一个是最近的psf
print("CHECK:::", cen_col[neigh[0]], cen_row[neigh[0]], cen_col[tpsf], cen_row[tpsf], cen_col[neigh[0]]-cen_col[tpsf], cen_row[neigh[0]]-cen_row[tpsf])
psfMatX = np.zeros([npsfX-1, ngy, ngx])
cen_colX= np.zeros(npsfX-1)
cen_rowX= np.zeros(npsfX-1)
for ipsf in range(npsfX):
if ipsf == 0:
continue
psfMatX[ipsf-1, :, :] = PSFMat[neigh[ipsf], :, :]
cen_colX[ipsf-1] = cen_col[neigh[ipsf]]
cen_rowX[ipsf-1] = cen_row[neigh[ipsf]]
psfMaker = np.zeros((ngy, ngx), dtype='float64')
if OnlyNeighbors == True:
PCAbasef, PCAcoeff = psfPCA_generator(psfMatX, npsfX-1, ngx, libPCApath)
nPCA = npsfX-1
for iPCA in range(nPCA):
coeffX = fitPoly(ref_col, ref_row, cen_colX, cen_rowX, PCAcoeff[:, iPCA], order=2)
psfMaker += coeffX*PCAbasef[iPCA, :, :]
return psfMaker
def test_loadPSF():
iccd = 1 #[1, 30]
iwave= 1 #[1, 4]
ipsf = 1 #[1, 100]
psfPath = '/Users/chengliangwei/csstPSFdata/CSSOS_psf_ciomp'
psfSet = []
for ipsf in range(1, 901):
psfInfo = LoadPSF(iccd, iwave, ipsf, psfPath, CalcPSFsize=True, CalcPSFcenter=True, SigRange=False)
psfSet.append(psfInfo)
print('psfSet has been loaded.')
print('Usage: psfSet[i][keys]')
print('psfSet.keys:', psfSet[0].keys())
return psfSet
def test_psfPCA():
#load psf
print('load psf...')
psfSet = test_loadPSF()
#set input for psfPCA calc.
print('PCA calc...')
npsf = 5
npix = 180
psfMat = np.zeros([npsf, npix, npix])
libPCApath = './libPCA/libPCA.so'
for ipsf in range(5):
psfMat[ipsf, :, :] = psfSet[ipsf]['psfMat']
PCAbasef, PCAcoeff = psfPCA_generator(psfMat, npsf, npix, libPCApath)
#plot check
print('plot...')
fig = plt.figure(figsize=(20, 10))
cc = 90
dcc= 15
ax = plt.subplot(2, 4, 1)
plt.imshow(PCAbasef[0, :, :][cc-dcc:cc+dcc, cc-dcc:cc+dcc], origin='lower')
ax = plt.subplot(2, 4, 2)
plt.imshow(PCAbasef[1, :, :][cc-dcc:cc+dcc, cc-dcc:cc+dcc], origin='lower')
ax = plt.subplot(2, 4, 3)
plt.imshow(PCAbasef[2, :, :][cc-dcc:cc+dcc, cc-dcc:cc+dcc], origin='lower')
ax = plt.subplot(2, 4, 4)
plt.imshow(PCAbasef[3, :, :][cc-dcc:cc+dcc, cc-dcc:cc+dcc], origin='lower')
ax = plt.subplot(2, 4, 5)
plt.imshow(PCAbasef[4, :, :][cc-dcc:cc+dcc, cc-dcc:cc+dcc], origin='lower')
ax = plt.subplot(2, 4, 6)
ipsf = 1
im = PCAcoeff[ipsf,0]*PCAbasef[0, :, :]
im+= PCAcoeff[ipsf,1]*PCAbasef[1, :, :]
im+= PCAcoeff[ipsf,2]*PCAbasef[2, :, :]
im+= PCAcoeff[ipsf,3]*PCAbasef[3, :, :]
im+= PCAcoeff[ipsf,4]*PCAbasef[4, :, :]
plt.imshow(im[cc-dcc:cc+dcc, cc-dcc:cc+dcc], origin='lower')
plt.colorbar()
ax = plt.subplot(2, 4, 8)
plt.imshow(psfMat[ipsf,:,:][cc-dcc:cc+dcc, cc-dcc:cc+dcc], origin='lower')
plt.colorbar()
ax = plt.subplot(2, 4, 7)
plt.imshow((psfMat[ipsf,:,:]-im)[cc-dcc:cc+dcc, cc-dcc:cc+dcc], origin='lower')
plt.colorbar()
plt.show()
def test_fitPoly():
datax = np.random.random(10)
datay = np.random.random(10)*2
dataz = datay**2 + datay*datax+10*datay**2
px = 0.28
py = 10.34
pz = fitPoly(px, py, datax, datay, dataz, order=2)
print('check: pz-out:', pz, 'pz-in', py**2+px*py+10*py**2)
##################################################
# __main__ #
##################################################
if __name__ == '__main__':
print('PSF modules for csst-imgsim')
#test_loadPSF()
#test_psfPCA()
test_fitPoly()
"""
CSST image simulation module (in python3): Point Spread Function (PSF)
author:: Wei Chengliang <chengliangwei@pmo.ac.cn>
"""
import sys
from itertools import islice
import numpy as np
import matplotlib.pyplot as plt
import scipy.io
from scipy.io import loadmat
#import xlrd
from scipy import ndimage
from scipy.interpolate import RectBivariateSpline
#from astropy.modeling.models import Ellipse2D
#from astropy.coordinates import Angle
#import matplotlib.patches as mpatches
import ctypes
import galsim
def setupPSFimg(iccd, iwave, psfPath="/data/simudata/CSSOSDataProductsSims/data/csstPSFdata/CSSOS_psf_ciomp"):
"""
psf model setup for csst-sim
Parameters:
iccd, iwave (int, int): psf model on iccd & iwave
psfPath (string, optional): path to psf matrix
Returns:
psf_model (psf_class): psf model
Methods:
psf_model.PSFinplace(self, px, py, interpScheme=1): psf interpolation
psf_model.PSFspin(self, psf, sigSpin, sigGauss, dx, dy): psf rotation (from Yudong)
"""
psf_model = PSFimg(iccd, iwave, psfPath)
return psf_model
##################################################
# A. psf matrix loading & checking #
##################################################
def psfPixelLayout(nrows, ncols, cenPosRow, cenPosCol, pixSizeInMicrons=5.0):
"""
convert psf pixels to physical position
Parameters:
nrows, ncols (int, int): psf sampling with [nrows, ncols].
cenPosRow, cenPosCol (float, float): A physical position of the chief ray for a given psf.
pixSizeInMicrons (float-optional): The pixel size in microns from the psf sampling.
Returns:
psfPixelPos (numpy.array-float): [posx, posy] in mm for [irow, icol]
Notes:
1. show positions on ccd, but not position on image only [+/- dy]
"""
psfPixelPos = np.zeros([2, nrows, ncols])
if nrows % 2 != 0:
sys.exit()
if ncols % 2 != 0:
sys.exit()
cenPix_row = nrows/2 + 1 #中心主光线对应pixle [由长光定义]
cenPix_col = ncols/2 + 1
for irow in range(nrows):
for icol in range(ncols):
delta_row = ((irow + 1) - cenPix_row)*pixSizeInMicrons*1e-3
delta_col = ((icol + 1) - cenPix_col)*pixSizeInMicrons*1e-3
psfPixelPos[0, irow, icol] = cenPosCol + delta_col
psfPixelPos[1, irow, icol] = cenPosRow - delta_row #note-1
return psfPixelPos
def imSigRange(img, fraction=0.80):
"""
extract the image within x-percent (DISCARD)
Parameters:
img (numpy.array-float): image
fraction (float-optional): a percentage
Returns:
im1 (numpy.array-float): image
"""
im1 = img.copy()
im1size = im1.shape
im2 = np.sort(im1.reshape(im1size[0]*im1size[1]))
im2 = im2[::-1]
im3 = np.cumsum(im2)/np.sum(im2)
loc = np.where(im3 > fraction)
#print(im3[loc[0][0]], im2[loc[0][0]])
im1[np.where(im1 <= im2[loc[0][0]])]=0
return im1
def imPlotInRange(img):
"""
plot image within a selected range
Parameters:
img (numpy.array-float): image
Returns:
"""
im1 = img.copy()
im1size = im1.shape
X,Y = np.meshgrid(range(im1size[1]),range(im1size[0]))
Z = im1
resolution = 25
f = lambda x,y: Z[int(y),int(x) ]
g = np.vectorize(f)
x = np.linspace(0,Z.shape[1], Z.shape[1]*resolution)
y = np.linspace(0,Z.shape[0], Z.shape[0]*resolution)
X2, Y2= np.meshgrid(x[:-1],y[:-1])
Z2 = g(X2,Y2)
#plt.pcolormesh(X,Y, Z)
#plt.imshow(img, origin='lower')
plt.contour(X2-0.5,Y2-0.5,Z2, [0.], colors='w', linestyles='--', linewidths=[1])
return
def findMaxPix(img):
"""
get the pixel position of the maximum-value
Parameters:
img (numpy.array-float): image
Returns:
imgMaxPix_x, imgMaxPix_y (int, int): pixel position in columns & rows
"""
maxIndx = np.argmax(img)
maxIndx = np.unravel_index(maxIndx, np.array(img).shape)
imgMaxPix_x = maxIndx[1]
imgMaxPix_y = maxIndx[0]
return imgMaxPix_x, imgMaxPix_y
def psfTailor(img, apSizeInArcsec=0.5, psfSampleSizeInMicrons=5, focalLengthInMeters=28):
"""
psf tailor within a given aperture size
Parameters:
img (numpy.array-float): image
apSizeInArcsec (float-optional): aperture size in arcseconds.
psfSampleSizeInMicrons (float-optional): psf pixel size in microns.
focalLengthInMeters (float-optional): csst focal length im meters.
Returns:
imgT (numpy.array-float): image
"""
imgMaxPix_x, imgMaxPix_y = findMaxPix(img)
apSizeInMicrons = np.deg2rad(apSizeInArcsec/3600.)*focalLengthInMeters*1e6
apSizeInPix = apSizeInMicrons/psfSampleSizeInMicrons
apSizeInPix = np.int(np.ceil(apSizeInPix))
imgT = np.zeros_like(img)
imgT[imgMaxPix_y-apSizeInPix:imgMaxPix_y+apSizeInPix+1,
imgMaxPix_x-apSizeInPix:imgMaxPix_x+apSizeInPix+1] = \
img[imgMaxPix_y-apSizeInPix:imgMaxPix_y+apSizeInPix+1,
imgMaxPix_x-apSizeInPix:imgMaxPix_x+apSizeInPix+1]
return imgT
def psfEncircle(img, fraction=0.8, psfSampleSizeInMicrons=5, focalLengthInMeters=28):
"""
psf tailor within a given percentage.
Parameters:
img (numpy.array-float): image
fraction (float-optional): a percentage for psf tailor.
psfSampleSizeInMicrons (float-optional): psf pixel size in microns.
focalLengthInMeters (float-optional): csst focal length im meters.
Returns:
img*wgt (numpy.array-float): image
REE80 (float): radius of REE80 in arcseconds.
"""
imgMaxPix_x, imgMaxPix_y = findMaxPix(img)
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]
wgt = np.ones_like(dis)
wgt[np.where(dis > dis[np.where(img == imgY[ind])])] = 0
REE80 = np.rad2deg(dis[np.where(img == imgY[ind])]*psfSampleSizeInMicrons*1e-6/focalLengthInMeters)*3600
return img*wgt, REE80
def psfSecondMoments(psfMat, cenX, cenY, pixSize=1):
"""
estimate the psf ellipticity by the second moment of surface brightness
Parameters:
psfMat (numpy.array-float): image
cenX, cenY (float, float): pixel position of the psf center
pixSize (float-optional): pixel size
Returns:
sz (float): psf size
e1, e2 (float, float): psf ellipticity
"""
I = psfMat
ncol = I.shape[1]
nrow = I.shape[0]
w = 0.0
w11 = 0.0
w12 = 0.0
w22 = 0.0
for icol in range(ncol):
for jrow in range(nrow):
x = icol*pixSize - cenX
y = jrow*pixSize - cenY
w += I[jrow, icol]
w11 += x*x*I[jrow, icol]
w12 += x*y*I[jrow, icol]
w22 += y*y*I[jrow, icol]
w11 /= w
w12 /= w
w22 /= w
sz = w11 + w22
e1 = (w11 - w22)/sz
e2 = 2.0*w12/sz
return sz, e1, e2
def LoadPSF(iccd, iwave, ipsf, psfPath, psfSampleSize=5, CalcPSFsize=True, CalcPSFcenter=True, SigRange=False, TailorScheme=1, InputMaxPixelPos=False):
'''加载psf信息'''
"""
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
psfSampleSize (float-optional): psf size in microns.
CalcPSFsize (bool-optional): whether calculate psf size & ellipticity. Default: True
CalcPSFcenter (bool-optional): whether calculate psf center. Default: True
SigRange (bool-optional): whether use psf tailor. Default: False
TailorScheme (int-optional): which method for psf tailor. Default: 1
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矩阵
fpathMat = fpath +'/' +'5_psf_array' +'/' +'psf_{:}.mat'.format(ipsf)
data = scipy.io.loadmat(fpathMat)
psfInfo['psfMat'] = data['psf']
#获取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
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
#nrows = 180 #psf采样大小, in pixels
#ncols = 180
nrows, ncols = psfInfo['psfMat'].shape
psfPos = psfPixelLayout(nrows, ncols, psfInfo['image_y'], psfInfo['image_x'], pixSizeInMicrons=5.0)
imgMaxPix_x, imgMaxPix_y = findMaxPix(psfInfo['psfMat'])
psfInfo['imgMaxPosx_ccd'] = psfPos[0, imgMaxPix_y, imgMaxPix_x] #cx, psf最大值位置, in mm
psfInfo['imgMaxPosy_ccd'] = psfPos[1, imgMaxPix_y, imgMaxPix_x] #cy
#计算psf size & ellipticity
if CalcPSFsize is True:
psfMat = psfInfo['psfMat'].copy()
cenX, cenY, sz, e1, e2, REE80 = psfSizeCalculator(psfMat, psfSampleSize=psfSampleSize, CalcPSFcenter=CalcPSFcenter, SigRange=SigRange, TailorScheme=TailorScheme)
psfInfo['psfCenX_img'] = cenX #in local pixels, psf质心位置, in pixels
psfInfo['psfCenY_img'] = cenY #in local pixels
psfInfo['psfSize'] = sz
psfInfo['psf_e1'] = e1
psfInfo['psf_e2'] = e2
psfInfo['REE80'] = REE80
return psfInfo
def psfSizeCalculator(psfMat, psfSampleSize=5, CalcPSFcenter=True, SigRange=False, TailorScheme=1):
"""
calculate psf size & ellipticity
Parameters:
psfMat (numpy.array): image
psfSampleSize (float-optional): psf size in microns.
CalcPSFcenter (bool-optional): whether calculate psf center. Default: True
SigRange (bool-optional): whether use psf tailor. Default: False
TailorScheme (int-optional): which method for psf tailor. Default: 1
Returns:
cenX, cenY (float, float): the pixel position of the mass center
sz (float): psf size
e1, e2 (float, float): psf ellipticity
REE80 (float): radius of REE80 in arcseconds
"""
psfSampleSize = psfSampleSize*1e-3 #mm
REE80 = -1.0 ##encircling 80% energy
if SigRange is True:
if TailorScheme == 1:
psfMat = imSigRange(psfMat, fraction=0.80)
psfInfo['psfMat'] = psfMat #set on/off
if TailorScheme == 2:
img = psfTailor(psfMat, apSizeInArcsec=0.5)
imgX, REE80 = psfEncircle(psfMat)
psfMat = img
REE80 = REE80[0]
if CalcPSFcenter is True:
img = psfMat/np.sum(psfMat)
y,x = ndimage.center_of_mass(img) #y-rows, x-cols
cenX = x
cenY = y
if CalcPSFcenter is False:
cenPix_X = psfMat.shape[1]/2 #90
cenPix_Y = psfMat.shape[0]/2 #90
cenX = cenPix_X + psfInfo['centroid_x']/psfSampleSize
cenY = cenPix_Y - psfInfo['centroid_y']/psfSampleSize
pixSize = 1
sz, e1, e2 = psfSecondMoments(psfMat, cenX, cenY, pixSize=pixSize)
return cenX, cenY, sz, e1, e2, REE80
def psfStack(*psfMat):
"""
stacked image from the different psfs
Parameters:
*psfMat (numpy.array): the different psfs for stacking
Returns:
img (numpy.array): image
"""
nn = len(psfMat)
img = np.zeros_like(psfMat[0])
for ii in range(nn):
img += psfMat[ii]/np.sum(psfMat[ii])
img /= np.sum(img)
return img
##################################################
# B. psf interpolation #
##################################################
def img2fits(img, fitsName=None):
"""
saving image to fits file
Parameters:
img (numpy.array): image
fitsName (string): path+filename of fits
Returns:
"""
from astropy.io import fits
grey = fits.PrimaryHDU(img)
greyHDU = fits.HDUList([grey])
if fitsName != None:
greyHDU.writeto(fitsName)
def psfMatLoad(iccd, iwave, psfPath, psfSampleSize=5, CalcPSFsize=False, CalcPSFcenter=True):
"""
load psf for interpolation
Parameters:
iccd, iwave, psfPath: # of ccd/wave and path for psfs
CalcPSFsize (bool-optional): whether calculate psf size & ellipticity. Default: False
CalcPSFcenter (bool-optional): whether calculate psf center. Default: True
Returns:
PSFMat (numpy.array): images
cen_col, cen_row (numpy.array, numpy.array): position of psf center in the view field
"""
psfSet = []
for ipsf in range(1, 901):
psfInfo = LoadPSF(iccd, iwave, ipsf, psfPath, CalcPSFsize=CalcPSFsize, CalcPSFcenter=CalcPSFcenter, InputMaxPixelPos=True)
psfSet.append(psfInfo)
npsf = len(psfSet)
ngy, ngx = psfSet[0]['psfMat'].shape
PSFMat = np.zeros([npsf, ngy, ngx])
cen_col= np.zeros(npsf)
cen_row= np.zeros(npsf)
FieldPos = False
for ipsf in range(npsf):
PSFMat[ipsf, :, :] = psfSet[ipsf]['psfMat']
if FieldPos == True:
cen_col[ipsf] = psfSet[ipsf]['field_x'] #cx
cen_row[ipsf] = psfSet[ipsf]['field_y'] #cy
if FieldPos == False:
cen_col[ipsf] = psfSet[ipsf]['imgMaxPosx_ccd'] #cx
cen_row[ipsf] = psfSet[ipsf]['imgMaxPosy_ccd'] #cy
return PSFMat, cen_col, cen_row
def findNeighbors(tx, ty, px, py, dr=0.1, dn=1, OnlyDistance=True):
"""
find nearest meighbors by 2D-KDTree
Parameters:
tx, ty (float, float): a given position
px, py (numpy.array, numpy.array): position data for tree
dr (float-optional): distance
dn (int-optional): nearest-N
OnlyDistance (bool-optional): only use distance to find neighbors. Default: True
Returns:
dataq (numpy.array): index
"""
import scipy.spatial as spatial
datax = px
datay = py
tree = spatial.KDTree(list(zip(datax.ravel(), datay.ravel())))
dataq=[]
rr = dr
if OnlyDistance == True:
dataq = tree.query_ball_point([tx, ty], rr)
if OnlyDistance == False:
while len(dataq) < dn:
dataq = tree.query_ball_point([tx, ty], rr)
rr += dr
dd = np.hypot(datax[dataq]-tx, datay[dataq]-ty)
ddSortindx = np.argsort(dd)
dataq = np.array(dataq)[ddSortindx[0:dn]]
return dataq
def psfCentering(img, apSizeInArcsec=0.5, psfSampleSizeInMicrons=5, focalLengthInMeters=28, CenteringMode=1):
"""
centering psf within an aperture
Parameters:
img (numpy.array): image
apSizeInArcsec (float-optional): aperture size in arcseconds.
psfSampleSizeInMicrons (float-optional): psf pixel size in microns.
focalLengthInMeters (float-optional): csst focal length im meters.
CenteringMode (int-optional): how to center psf images
Returns:
imgT (numpy.array)
"""
if CenteringMode == 1:
imgMaxPix_x, imgMaxPix_y = findMaxPix(img)
if CenteringMode == 2:
y,x = ndimage.center_of_mass(img) #y-rows, x-cols
imgMaxPix_x = int(x)
imgMaxPix_y = int(y)
apSizeInMicrons = np.deg2rad(apSizeInArcsec/3600.)*focalLengthInMeters*1e6
apSizeInPix = apSizeInMicrons/psfSampleSizeInMicrons
apSizeInPix = np.int(np.ceil(apSizeInPix))
imgT = np.zeros_like(img)
ngy, ngx = img.shape
cy = int(ngy/2)
cx = int(ngx/2)
imgT[cy-apSizeInPix:cy+apSizeInPix+1,
cx-apSizeInPix:cx+apSizeInPix+1] = \
img[imgMaxPix_y-apSizeInPix:imgMaxPix_y+apSizeInPix+1,
imgMaxPix_x-apSizeInPix:imgMaxPix_x+apSizeInPix+1]
return imgT
def psfMaker_IDW(px, py, PSFMat, cen_col, cen_row, IDWindex=2, OnlyNeighbors=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])
if OnlyNeighbors == True:
neigh = findNeighbors(px, py, cen_col, cen_row, dr=5., dn=9, OnlyDistance=False)
neighFlag = np.zeros(npsf)
neighFlag[neigh] = 1
print("neigh:", neigh)
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='float64')
for ipsf in range(npsf):
if OnlyNeighbors == True:
if neighFlag[ipsf] != 1:
continue
iPSFMat = PSFMat[ipsf, :, :].copy()
iPSFMat = psfCentering(iPSFMat, CenteringMode=1)
ipsfWeight = psfWeight[ipsf]
psfMaker += iPSFMat * ipsfWeight
psfMaker /= np.nansum(psfMaker)
return psfMaker
def psfMaker_PCA(px, py, PSFMat, cen_col, cen_row, OnlyNeighbors=False, libPCApath='libPCA/libPCA.so'):
"""
psf interpolation by PCA
Parameters:
Returns:
"""
ref_col = px
ref_row = py
ngy, ngx = PSFMat[0, :, :].shape
npsf = PSFMat[:, :, :].shape[0]
neigh = findNeighbors(px, py, cen_col, cen_row, dr=0.3, dn=5, OnlyDistance=False)
npsfX = len(neigh)
psfMatX = np.zeros([npsfX, ngy, ngx])
cen_colX= np.zeros(npsfX)
cen_rowX= np.zeros(npsfX)
for ipsf in range(npsfX):
psfMatX[ipsf, :, :] = PSFMat[neigh[ipsf], :, :]
cen_colX[ipsf] = cen_col[neigh[ipsf]]
cen_rowX[ipsf] = cen_row[neigh[ipsf]]
psfMaker = np.zeros((ngy, ngx), dtype='float64')
if OnlyNeighbors == True:
PCAbasef, PCAcoeff = psfPCA_generator(psfMatX, npsfX, ngx, libPCApath)
nPCA = npsfX
for iPCA in range(nPCA):
coeffX = fitPoly(ref_col, ref_row, cen_colX, cen_rowX, PCAcoeff[:, iPCA], order=2)
psfMaker += coeffX*PCAbasef[iPCA, :, :]
return psfMaker
def psfPCA_generator(psfMat, npsf, npix, libPCApath):
"""
generate PCs from psfMat
Parameters:
Returns:
"""
libPCA = ctypes.CDLL(libPCApath) # CDLL加载库
libPCA.psfPCA.argtypes = [ctypes.POINTER(ctypes.c_float), ctypes.c_int, ctypes.c_int, ctypes.POINTER(ctypes.c_double), ctypes.POINTER(ctypes.c_double)]
Nstar = npsf
Mp = npix*npix
NM = Nstar*Mp
NN = Nstar*Nstar
arr = (ctypes.c_float*NM)()
basef = (ctypes.c_double*NM)()
coeff = (ctypes.c_double*NN)()
psfT = np.zeros(NM)
for ipsf in range(npsf):
lp = 0 + ipsf*Mp
up = Mp+ ipsf*Mp
ipsfMat = psfMat[ipsf, :, :]
psfT[lp:up] = ipsfMat.reshape(Mp)
arr[:] = psfT
libPCA.psfPCA(arr, Nstar, Mp, basef, coeff)
PCAbasef = np.zeros([npsf, npix, npix])
PCAcoeff = np.zeros([npsf, npsf])
for ipsf in range(npsf):
lp = 0 + ipsf*Mp
up = Mp+ ipsf*Mp
PCAbasef[ipsf, :, :] = np.array(basef[lp:up]).reshape(npix, npix)
lp = 0 + ipsf*npsf
up = npsf+ ipsf*npsf
PCAcoeff[ipsf, :] = np.array(coeff[lp:up])
return PCAbasef, PCAcoeff
def fitPoly(px, py, datax, datay, dataz, order = 2):
if order == 1:
# best-fit linear plane
A = np.c_[datax, datay, np.ones(datax.shape[0])]
C,_,_,_ = scipy.linalg.lstsq(A, dataz) # coefficients
pz = C[0]*px + C[1]*py + C[2]
elif order == 2:
# best-fit quadratic curve
A = np.c_[np.ones(datax.shape[0]), np.c_[datax, datay], np.prod(np.c_[datax, datay], axis=1), np.c_[datax, datay]**2]
C,_,_,_ = scipy.linalg.lstsq(A, dataz)
pz = np.dot(np.c_[1, px, py, px*py, px**2, py**2], C)
"""
elif order == 3:
# best-fit cubic curve
A = np.c_[np.ones(datax.shape[0]), np.c_[datax, datay], np.prod(np.c_[datax, datay], axis=1), np.c_[datax, datay]**2, np.c_[datax, datay]**3]
C,_,_,_ = scipy.linalg.lstsq(A, dataz)
pz = np.dot(np.c_[1, px, py, px*py, px**2, py**2, px**3, py**3], C)
"""
return pz
"""
############################
### not used temporarily ###
############################
def psfSplineMake(px, py, PSFMat, cen_col, cen_row, OnlyNeighbors=False):
minimum_psf_weight = 1e-8
ref_col = px
ref_row = py
cdelt1p = 1
cdelt2p = 1
ngy, ngx = PSFMat[0, :, :].shape
psfx = np.linspace(0, ngx-1, ngx)
psfy = np.linspace(0, ngy-1, ngy)
npsf = PSFMat[:, :, :].shape[0]
psfWeight = np.zeros([npsf])
for ipsf in range(npsf):
psfWeight[ipsf] = np.sqrt((ref_col - cen_col[ipsf])**2 + (ref_row - cen_row[ipsf])**2)
psfWeight[ipsf] = max(psfWeight[ipsf], minimum_psf_weight)
psfWeight[ipsf] = 1./psfWeight[ipsf]
psfWeight /= np.sum(psfWeight)
psf = np.zeros((ngy, ngx), dtype='float64')
for ipsf in range(npsf):
iPSFMat = PSFMat[ipsf, :, :]
ipsfWeight = psfWeight[ipsf]
psf += iPSFMat * ipsfWeight
psf /= (np.nansum(psf) * cdelt1p * cdelt2p)
psfSpline = RectBivariateSpline(psfy, psfx, psf)
return psf, psfSpline
def psfToImage(psfSpline, cutoff_radius=180):
ng = 180
img = np.zeros([ng, ng], dtype='float64')
for i in range(ng):
for j in range(ng):
star_row = 5
star_column = 5
if np.sqrt((j-star_column)**2 + (i-star_row)**2) <= cutoff_radius:
star_flux = 8
column_cen = j #j - star_column
row_cen = i #i - star_row
img[i,j] += star_flux * psfSpline.integral(row_cen-0.5, row_cen+0.5, column_cen-0.5, column_cen+0.5)
return img
"""
##################################################
# C. csstPSF class #
##################################################
class PSFimg(object):
def __init__(self, iccd, iwave, psfPath):
self.iccd = iccd
self.iwave= iwave
self.psfPath = psfPath
#loading psfSet >>>
"""
psfSet = []
for ipsf in range(1, 901):
psfInfo = LoadPSF(iccd, iwave, ipsf, psfPath, CalcPSFsize=True, CalcPSFcenter=True, SigRange=False)
psfSet.append(psfInfo)
self.psfSet = psfSet
"""
a, b, c = psfMatLoad(iccd, iwave, psfPath)
self.psfMat = a
self.cenPosx= b
self.cenPosy= c
def PSFinplace(self, px, py, interpScheme=1):
if interpScheme == 1:
idwIndx = 2
psf = psfMaker_IDW(px, py, self.psfMat, self.cenPosx, self.cenPosy, IDWindex=idwIndx, OnlyNeighbors=True)
if interpScheme ==2:
libPCA = "/Users/chengliangwei/Desktop/csstPSF/libPCA/libPCA.so"
psf = psfMaker_PCA(px, py, self.psfMat, self.cenPosx, self.cenPosy, OnlyNeighbors=True, libPCApath=libPCA)
img = galsim.ImageF(psf, scale=0.074/2)
xpsf = galsim.InterpolatedImage(img)
return xpsf
"""
def gcPlot(self, psf,pscale=0.074,figout="GC.png"):
size = np.size(psf,axis=0)
cxy = 0.5*(size-1)
width = 0.5*size
# log scale
radius = np.arange(np.log10(0.2),np.log10(width),0.01)
radius = 10.0**radius
nr = len(radius)
gc = []
for i in range(nr):
iflux, iferr, xflag = sep.sum_circle(psf,cxy,cxy,radius[i],subpix=0)
gc += [iflux.tolist()]
# Estimate the radius for a given flux ratio
fratio = 0.8
mid = [i for i in range(nr) if gc[i]<=fratio and gc[i+1]>fratio][0]
r0, r1 = radius[mid], radius[mid+1]
gc0, gc1 = gc[mid], gc[mid+1]
r5 = (fratio-gc0)/(gc1-gc0)*(r1-r0) + r0
hlf = r5*pscale
# plot
pfit = interp1d(radius, gc, kind='cubic')
fig = pl.figure(figsize=(5.5,4.0))
ax = fig.add_axes([0.16,0.15,0.80,0.81])
ax.plot(radius*pscale, pfit(radius), "k", linewidth=2.0)
ax.plot(radius*pscale, gc, "*r", markersize=5.0,mec="r",alpha=0.3)
ax.plot([hlf,hlf],[0,fratio],"k",linewidth=2.5)
ax.plot([0,hlf],[fratio,fratio],"k",linewidth=2.5)
ax.text(radius[10]*pscale,0.6,"$r_{%.1f}$=%.2f\""%(fratio,hlf))
ax.set_xlabel("Radius (arcsec)",fontsize=15)
ax.set_ylabel("Growth of Curve",fontsize=15)
ax.set_xscale("log")
ax.set_xlim(radius[0]*pscale,radius[-1]*pscale)
ax.set_ylim(0.0,1.0)
for tick in ax.xaxis.get_major_ticks(): tick.label.set_fontsize(15)
for tick in ax.yaxis.get_major_ticks(): tick.label.set_fontsize(15)
pl.savefig(figout)
pl.clf()
pl.close()
return
"""
def PSFspin(self, psf, sigSpin, sigGauss, dx, dy):
"""
The PSF profile at a given image position relative to the axis center
Parameters:
theta : spin angles in a given exposure in unit of [arcsecond]
dx, dy: relative position to the axis center in unit of [pixels]
Return:
Spinned PSF: g1, g2 and axis ratio 'a/b'
"""
a2Rad = np.pi/(60.0*60.0*180.0)
ff = sigGauss * 0.107 * (1000.0/10.0) # in unit of [pixels]
rc = np.sqrt(dx*dx + dy*dy)
cpix = rc*(sigSpin*a2Rad)
beta = (np.arctan2(dy,dx) + np.pi/2)
ell = cpix**2/(2.0*ff**2+cpix**2)
#ell *= 10.0
qr = np.sqrt((1.0+ell)/(1.0-ell))
#psfShape = galsim.Shear(e=ell, beta=beta)
#g1, g2 = psfShape.g1, psfShape.g2
#qr = np.sqrt((1.0+ell)/(1.0-ell))
#return ell, beta, qr
PSFshear = galsim.Shear(e=ell, beta=beta*galsim.radians)
return psf.shear(PSFshear), PSFshear(base)
##################################################
# D. TEST #
##################################################
def psfMaker_IDW_test(tpsf, px, py, PSFMat, cen_col, cen_row, IDWindex=2, OnlyNeighbors=False):
minimum_psf_weight = 1e-8
ref_col = px
ref_row = py
ngy, ngx = PSFMat[0, :, :].shape
npsf = PSFMat[:, :, :].shape[0]
psfWeight = np.zeros([npsf])
if OnlyNeighbors == True:
neigh = findNeighbors(px, py, cen_col, cen_row, dr=0.1, dn=9, OnlyDistance=False)
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='float64')
for ipsf in range(npsf):
"""
if OnlyNeighbors == True:
iy, ix = np.unravel_index(ipsf, (10,10))
ty, tx = np.unravel_index(tpsf, (10,10))
if np.abs(iy - ty) > 1 or np.abs(ix - tx) > 1:
continue
"""
if OnlyNeighbors == True:
if neighFlag[ipsf] != 1:
continue
if ipsf == tpsf:
continue
iPSFMat = PSFMat[ipsf, :, :].copy()
iPSFMat = psfCentering(iPSFMat, CenteringMode=1)
ipsfWeight = psfWeight[ipsf]
psfMaker += iPSFMat * ipsfWeight
psfMaker /= np.nansum(psfMaker)
return psfMaker
def psfMaker_PCA_test(tpsf, px, py, PSFMat, cen_col, cen_row, OnlyNeighbors=False, libPCApath='libPCA/libPCA.so'):
"""
psf interpolation by PCA
Parameters:
Returns:
"""
ref_col = px
ref_row = py
ngy, ngx = PSFMat[0, :, :].shape
npsf = PSFMat[:, :, :].shape[0]
neigh = findNeighbors(px, py, cen_col, cen_row, dr=0.3, dn=9, OnlyDistance=False)
npsfX = len(neigh)
#去掉tpsf,neigh中排在第一个是最近的psf
print("CHECK:::", cen_col[neigh[0]], cen_row[neigh[0]], cen_col[tpsf], cen_row[tpsf], cen_col[neigh[0]]-cen_col[tpsf], cen_row[neigh[0]]-cen_row[tpsf])
psfMatX = np.zeros([npsfX-1, ngy, ngx])
cen_colX= np.zeros(npsfX-1)
cen_rowX= np.zeros(npsfX-1)
for ipsf in range(npsfX):
if ipsf == 0:
continue
psfMatX[ipsf-1, :, :] = PSFMat[neigh[ipsf], :, :]
cen_colX[ipsf-1] = cen_col[neigh[ipsf]]
cen_rowX[ipsf-1] = cen_row[neigh[ipsf]]
psfMaker = np.zeros((ngy, ngx), dtype='float64')
if OnlyNeighbors == True:
PCAbasef, PCAcoeff = psfPCA_generator(psfMatX, npsfX-1, ngx, libPCApath)
nPCA = npsfX-1
for iPCA in range(nPCA):
coeffX = fitPoly(ref_col, ref_row, cen_colX, cen_rowX, PCAcoeff[:, iPCA], order=2)
psfMaker += coeffX*PCAbasef[iPCA, :, :]
return psfMaker
def test_loadPSF():
iccd = 1 #[1, 30]
iwave= 1 #[1, 4]
ipsf = 1 #[1, 100]
psfPath = '/Users/chengliangwei/csstPSFdata/CSSOS_psf_ciomp'
psfSet = []
for ipsf in range(1, 901):
psfInfo = LoadPSF(iccd, iwave, ipsf, psfPath, CalcPSFsize=True, CalcPSFcenter=True, SigRange=False)
psfSet.append(psfInfo)
print('psfSet has been loaded.')
print('Usage: psfSet[i][keys]')
print('psfSet.keys:', psfSet[0].keys())
return psfSet
def test_psfPCA():
#load psf
print('load psf...')
psfSet = test_loadPSF()
#set input for psfPCA calc.
print('PCA calc...')
npsf = 5
npix = 180
psfMat = np.zeros([npsf, npix, npix])
libPCApath = './libPCA/libPCA.so'
for ipsf in range(5):
psfMat[ipsf, :, :] = psfSet[ipsf]['psfMat']
PCAbasef, PCAcoeff = psfPCA_generator(psfMat, npsf, npix, libPCApath)
#plot check
print('plot...')
fig = plt.figure(figsize=(20, 10))
cc = 90
dcc= 15
ax = plt.subplot(2, 4, 1)
plt.imshow(PCAbasef[0, :, :][cc-dcc:cc+dcc, cc-dcc:cc+dcc], origin='lower')
ax = plt.subplot(2, 4, 2)
plt.imshow(PCAbasef[1, :, :][cc-dcc:cc+dcc, cc-dcc:cc+dcc], origin='lower')
ax = plt.subplot(2, 4, 3)
plt.imshow(PCAbasef[2, :, :][cc-dcc:cc+dcc, cc-dcc:cc+dcc], origin='lower')
ax = plt.subplot(2, 4, 4)
plt.imshow(PCAbasef[3, :, :][cc-dcc:cc+dcc, cc-dcc:cc+dcc], origin='lower')
ax = plt.subplot(2, 4, 5)
plt.imshow(PCAbasef[4, :, :][cc-dcc:cc+dcc, cc-dcc:cc+dcc], origin='lower')
ax = plt.subplot(2, 4, 6)
ipsf = 1
im = PCAcoeff[ipsf,0]*PCAbasef[0, :, :]
im+= PCAcoeff[ipsf,1]*PCAbasef[1, :, :]
im+= PCAcoeff[ipsf,2]*PCAbasef[2, :, :]
im+= PCAcoeff[ipsf,3]*PCAbasef[3, :, :]
im+= PCAcoeff[ipsf,4]*PCAbasef[4, :, :]
plt.imshow(im[cc-dcc:cc+dcc, cc-dcc:cc+dcc], origin='lower')
plt.colorbar()
ax = plt.subplot(2, 4, 8)
plt.imshow(psfMat[ipsf,:,:][cc-dcc:cc+dcc, cc-dcc:cc+dcc], origin='lower')
plt.colorbar()
ax = plt.subplot(2, 4, 7)
plt.imshow((psfMat[ipsf,:,:]-im)[cc-dcc:cc+dcc, cc-dcc:cc+dcc], origin='lower')
plt.colorbar()
plt.show()
def test_fitPoly():
datax = np.random.random(10)
datay = np.random.random(10)*2
dataz = datay**2 + datay*datax+10*datay**2
px = 0.28
py = 10.34
pz = fitPoly(px, py, datax, datay, dataz, order=2)
print('check: pz-out:', pz, 'pz-in', py**2+px*py+10*py**2)
##################################################
# __main__ #
##################################################
if __name__ == '__main__':
print('PSF modules for csst-imgsim')
#test_loadPSF()
#test_psfPCA()
test_fitPoly()
import numpy as np
import matplotlib.pyplot as plt
import scipy.io
import galsim
vc_A = 2.99792458e+18 # speed of light: A/s
vc_m = 2.99792458e+8 # speed of light: m/s
h_Plank = 6.626196e-27 # Plank constant: erg s
def photonEnergy(lambd):
nu = vc_A / lambd
eph = h_Plank * nu
return eph
pixSize = 0.037
mag_star = 15.
Nx = 256
Ny = 256
###加载PSF信息###
def LoadPSF(iccd, iwave, ipsf, psfPath, psfSampleSize=5, PSFCentroidWgt=False):
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']
return psfInfo
def psfCentering(img, apSizeInArcsec=0.5, psfSampleSizeInMicrons=5, focalLengthInMeters=28, CenteringMode=1):
imgMaxPix_x, imgMaxPix_y = findMaxPix(img)
apSizeInMicrons = np.deg2rad(apSizeInArcsec/3600.)*focalLengthInMeters*1e6
apSizeInPix = apSizeInMicrons/psfSampleSizeInMicrons
apSizeInPix = np.int(np.ceil(apSizeInPix))
imgT = np.zeros_like(img)
ngy, ngx = img.shape
cy = int(ngy/2)
cx = int(ngx/2)
imgT[cy-apSizeInPix:cy+apSizeInPix+1,
cx-apSizeInPix:cx+apSizeInPix+1] = \
img[imgMaxPix_y-apSizeInPix:imgMaxPix_y+apSizeInPix+1,
imgMaxPix_x-apSizeInPix:imgMaxPix_x+apSizeInPix+1]
return imgT
def findMaxPix(img):
maxIndx = np.argmax(img)
maxIndx = np.unravel_index(maxIndx, np.array(img).shape)
imgMaxPix_x = maxIndx[1]
imgMaxPix_y = maxIndx[0]
return imgMaxPix_x, imgMaxPix_y
def magToFlux(mag):
flux = 10**(-0.4*(mag+48.6))
return flux
def getElectronFluxFilt(mag, exptime=150.):
pE = photonEnergy(lambd=6199.8)
flux = magToFlux(mag)
factor = 1.0e4 * flux/pE * vc_A * (1.0/5370.0 - 1.0/7030.0)
return factor * 0.5040 * np.pi * (0.5 * 2.0)**2 * exptime
def radial_profile(img, cx, cy, nbins=100, Rmin=16, Rmax=128):
nx = img.shape[0]
ny = img.shape[1]
x = np.arange(0, img.shape[0], 1)
y = np.arange(0, img.shape[1], 1)
xx, yy = np.meshgrid(x, y)
dist = np.sqrt((xx - cx)**2 + (yy - cy)**2)
img_flat = img.flatten()
dist_flat = dist.flatten()
a, bins = np.histogram(dist_flat, range=(Rmin, Rmax), bins=nbins, weights=img_flat)
b, bins = np.histogram(dist_flat, range=(Rmin, Rmax), bins=nbins)
b[b==0] = 1.
mid = (bins[0:-1]+bins[1:])/2.
return a / b, mid
if __name__ == '__main__':
iccd = 1
iwave = 1
ipsf = 1
psfPath= '/data/simudata/CSSOSDataProductsSims/data/csstPSFdata/CSSOS_psf_ciomp_30X30'
psfInfo= LoadPSF(iccd, iwave, ipsf, psfPath, psfSampleSize=5, PSFCentroidWgt=False)
ipsfMat = psfInfo['psfMat']
print(ipsfMat[0:100,0:100])
cgrid = 128
dgrid = 15
cx = int(Nx/2)
cy = int(Ny/2)
with np.printoptions(precision=5, suppress=True):
fig = plt.figure(figsize=(18,12))
ax = plt.subplot(2,2,1)
plt.imshow(ipsfMat[cgrid-dgrid:cgrid+dgrid, cgrid-dgrid:cgrid+dgrid], origin='lower')
plt.annotate('originalPSF', [0.1, 0.85], xycoords='axes fraction', color='w')
ax = plt.subplot(2,2,2)
img = galsim.ImageF(ncol=Nx, nrow=Ny, scale=pixSize)
psf_img = galsim.ImageF(ipsfMat, scale=pixSize)
psf = galsim.InterpolatedImage(psf_img)
psf_list = [psf] * 4
mag_star = 15
nphotons_tot = getElectronFluxFilt(mag=mag_star)
print(nphotons_tot)
for i in range(4):
nphotons = nphotons_tot / 4.
# star = galsim.Gaussian(sigma=1.e-8, flux=1.)
star = galsim.DeltaFunction()
star = star.withFlux(nphotons)
star = galsim.Convolve(psf_list[i], star)
img = star.drawImage(image=img, method='phot', center=galsim.PositionD(cx, cy), add_to_image=True)
print(img.array.shape)
plt.imshow(img.array[cx-dgrid:cx+dgrid, cy-dgrid:cy+dgrid], origin='lower')
plt.annotate('photon shooting 4 times', [0.1, 0.85], xycoords='axes fraction', color='w')
plt.colorbar()
ax = plt.subplot(2,2,3)
val, bins = radial_profile(img=img.array, cx=cx, cy=cy, nbins=100, Rmax=cx)
# plt.hist(val, bins=bins)
plt.plot(bins, val, label='mag=%d'%(mag_star))
img = galsim.ImageF(ncol=Nx, nrow=Ny, scale=pixSize)
psf_img = galsim.ImageF(ipsfMat, scale=pixSize)
psf = galsim.InterpolatedImage(psf_img)
psf_list = [psf] * 4
mag_star = 13
nphotons_tot = getElectronFluxFilt(mag=mag_star)
print(nphotons_tot)
for i in range(4):
nphotons = nphotons_tot / 4.
# star = galsim.Gaussian(sigma=1.e-8, flux=1.)
star = galsim.DeltaFunction()
star = star.withFlux(nphotons)
star = galsim.Convolve(psf_list[i], star)
img = star.drawImage(image=img, method='phot', center=galsim.PositionD(cx, cy), add_to_image=True)
print(img.array.shape)
val, bins = radial_profile(img=img.array, cx=cx, cy=cy, nbins=100, Rmax=cx)
# plt.hist(val, bins=bins)
plt.plot(bins, val, label='mag=%d'%(mag_star))
plt.legend(loc='upper right', fancybox=True)
plt.xlabel("R [pix]", size='x-large')
plt.ylabel("photon count", size='x-large')
plt.ylim([0, 500])
# print(img.array[0:100,0:100])
#psf Convolve galsim.DeltaFunction
#photon shooting ?
#plot image?
# ax = plt.subplot(2,2,2)
# plt.imshow(ipsfMat[cgrid-dgrid:cgrid+dgrid, cgrid-dgrid:cgrid+dgrid], origin='lower')
# plt.annotate('originalPSF', [0.1, 0.85], xycoords='axes fraction', color='w')
# ax = plt.subplot(2,2,4)
# img = galsim.ImageF(ncol=Nx, nrow=Ny, scale=pixSize)
# psf_img = galsim.ImageF(ipsfMat, scale=pixSize)
# psf = galsim.InterpolatedImage(psf_img)
# psf_list = [psf] * 4
# nphotons_tot = getElectronFluxFilt(mag=mag_star)
# print(nphotons_tot)
# obj_list = []
# for i in range(4):
# nphotons = nphotons_tot / 4.
# # star = galsim.Gaussian(sigma=1.e-8, flux=1.)
# star = galsim.DeltaFunction()
# star = star.withFlux(nphotons)
# star = galsim.Convolve(psf_list[i], star)
# obj_list.append(star)
# star = galsim.Sum(obj_list)
# img = star.drawImage(image=img, method='phot', center=galsim.PositionD(cx, cy), add_to_image=True)
# plt.annotate('photon shooting once', [0.1, 0.85], xycoords='axes fraction', color='w')
# print(img.array.shape)
# plt.imshow(img.array[cx-dgrid:cx+dgrid, cy-dgrid:cy+dgrid], origin='lower')
# plt.colorbar()
# print(img.array[0:100,0:100])
# ax = plt.subplot(2,3,2)
# imy = psfCentering(ipsfMat, apSizeInArcsec=2.0)
# plt.imshow(imy[cgrid-dgrid:cgrid+dgrid, cgrid-dgrid:cgrid+dgrid], origin='lower')
# plt.annotate('PSFCentroidOld', [0.1, 0.85], xycoords='axes fraction', color='w')
# ax = plt.subplot(2,3,5)
# #psf Convolve galsim.DeltaFunction
# #photon shooting ?
# #plot image?
# ax = plt.subplot(2,3,3)
# psfInfo= LoadPSF(iccd, iwave, ipsf, psfPath, psfSampleSize=5, PSFCentroidWgt=True)
# ipsfMat = psfInfo['psfMat']
# # print(ipsfMat[0:100,0:100])
# plt.imshow(ipsfMat[cgrid-dgrid:cgrid+dgrid, cgrid-dgrid:cgrid+dgrid], origin='lower')
# plt.annotate('PSFCentroidNew', [0.1, 0.85], xycoords='axes fraction', color='w')
# ax = plt.subplot(2,3,6)
# #psf Convolve galsim.DeltaFunction
# #photon shooting ?
# #plot image?
plt.savefig('testPlot.pdf')
......@@ -16,8 +16,4 @@ NP=40
date
echo $NP
# mpirun -np $NP --oversubscribe -H comput101 python /public/home/fangyuedong/CSST/ObservationSim/runExposure.py
# python /public/home/fangyuedong/sim_code_release/CSSTSim/ObservationSim/preprocess.py
# mpirun -np $NP python /public/home/fangyuedong/sim_code_release/CSSTSim/ObservationSim/runExposure.py
mpirun -np $NP python /public/home/fangyuedong/sim_code_release/CSSTSim/ObservationSim/runExposure.py config_sim.yaml -c /public/home/fangyuedong/sim_code_release/CSSTSim/config
mpirun -np $NP python /public/home/fangyuedong/test_release/CSST/ObservationSim/runExposure.py config_sim.yaml -c /public/home/fangyuedong/test_release/CSST/config
......@@ -10,7 +10,7 @@
# Can add some of the command-line arguments here as well;
# OK to pass either way or both, as long as they are consistent
# work_dir: "/public/home/fangyuedong/sim_code_release/CSST/test/"
work_dir: "/public/home/fangyuedong/sim_code_release/CSSTSim/workplace/"
work_dir: "/public/home/fangyuedong/test_release/CSST/workplace/"
data_dir: "/data/simudata/CSSOSDataProductsSims/data/"
run_name: "TEST"
......@@ -39,9 +39,9 @@ obs_setting:
# Options for survey types:
# "Photometric": simulate photometric chips only
# "Spectroscoplic": simulate slitless spectroscopic chips only
# "Spectroscopic": simulate slitless spectroscopic chips only
# "All": simulate full focal plane
survey_type: "Photometric"
survey_type: "All"
# Exposure time [seconds]
exp_time: 150.
......@@ -61,15 +61,15 @@ obs_setting:
# Number of calibration pointings
# Note: only valid when a pointing list is specified
np_cal: 1
np_cal: 0
# (Optional) only run specific pointing(s).
# Note: only valid when a pointing list is specified
run_pointings: [0, 1, 2, 3, 4]
run_pointings: [5, 6]
# (Optional) only run specific chip(s)
# Note: for all pointings
run_chips: [6, 25]
# run_chips: [6, 25]
###############################################
# Input path setting
......@@ -116,7 +116,11 @@ psf_setting:
# path to PSF data
# NOTE: only valid for "Interp" PSF
psf_dir: "csstPSFdata/CSSOS_psf_20210108/CSST_psf_ciomp_2p5um_cycle3_ccr90_proc"
psf_dir: "csstPSFdata/CSSOS_psf_20210108/CSST_psf_ciomp_2p5um_cycle3_ccr90"
# path to field-distortion model
# Note: only valid when ins_effects: field_dist is "ON"
fd_path: "FieldDistModelv2.0.pickle"
# sigma_spin: 0.0 # psf spin?
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment