Commit 18bdecdb authored by Fang Yuedong's avatar Fang Yuedong
Browse files

encapsulate instrumental data

parent 72f88213
import os, sys
import numpy as np
import scipy.io
import mpi4py.MPI as MPI
import ObservationSim.PSF.PSFInterp.PSFConfig as PSFConfig
import ObservationSim.PSF.PSFInterp.PSFUtil as PSFUtil
def mkdir(path):
isExists = os.path.exists(path)
if not isExists:
os.mkdir(path)
############################################
comm = MPI.COMM_WORLD
ThisTask = comm.Get_rank()
NTasks = comm.Get_size()
npsf = 900
psfPath = '/data/simudata/CSSOSDataProductsSims/data/csstPSFdata/CSSOS_psf_20210108/CSST_psf_ciomp_2p5um_cycle3_ccr90'
npsfPerTasks = int(npsf/NTasks)
iStart= 0 + npsfPerTasks*ThisTask
iEnd = npsfPerTasks + npsfPerTasks*ThisTask
if ThisTask == NTasks:
iEnd = npsf
for iccd in range(1, 31):
iccdPath = psfPath + '_proc/ccd{:}'.format(iccd)
if ThisTask == 0:
mkdir(iccdPath)
comm.barrier()
for iwave in range(1, 5):
iwavePath = iccdPath + '/wave_{:}'.format(iwave)
if ThisTask == 0:
mkdir(iwavePath)
comm.barrier()
psfMatPath = iwavePath + '/5_psf_array'
if ThisTask == 0:
mkdir(psfMatPath)
comm.barrier()
for ii in range(iStart, iEnd):
ipsf = ii+1
if ThisTask ==0:
print('iccd-iwave-ipsf: {:4}{:4}{:4}'.format(iccd, iwave, ipsf), end='\r',flush=True)
#if iccd != 1 or iwave !=1 or ipsf != 1:
# continue
ipsfOutput = psfMatPath + '/psf_{:}_centroidWgt.mat'.format(ipsf)
psfInfo = PSFConfig.LoadPSF(iccd, iwave, ipsf, psfPath, InputMaxPixelPos=True)
ipsfMat = psfInfo['psfMat']
npix_y, npix_x = ipsfMat.shape
#if npsf == 100:
# ncut = 160
#if npsf == 900:
# ncut = 200
ncut = int(npix_y*0.90)
ncut = ncut + ncut%2 #even pixs
img, cx, cy = PSFUtil.centroidWgt(ipsfMat, nt=ncut)
dcx = cx - npix_x/2 #pixel coords -> global coords
dcy = cy - npix_y/2 #pixel coords -> global coords
dcx*= psfInfo['pixsize']*1e-3 #5e-3 #pixels -> mm
dcy*= psfInfo['pixsize']*1e-3 #5e-3 #pixels -> mm
nn = npix_y
dn = int((nn - ncut)/2)
imgt = np.zeros([nn, nn], dtype=np.float32)
imgt[dn:-dn, dn:-dn] = img
scipy.io.savemat(ipsfOutput, {'cx':dcx, 'cy':dcy, 'psf':imgt})
if iccd != 1 or iwave !=1 or ipsf != 1:
if ThisTask == 0:
print('CHECK::', dcx, dcy, psfInfo['centroid_x'],psfInfo['centroid_y'])
import sys
import ctypes
import numpy as np
from scipy import ndimage
import scipy.spatial as spatial
###定义PSF像素的全局坐标###
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:in CCD全局坐标
return psfPixelPos
###查找最大pixel位置###
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
###查找neighbors位置###
def findNeighbors(tx, ty, px, py, dr=0.1, dn=1, OnlyDistance=True):
"""
find nearest neighbors 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
"""
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
###查找neighbors位置-hoclist###
def hocBuild(partx, party, nhocx, nhocy, dhocx, dhocy):
if np.max(partx) > nhocx*dhocx:
print('ERROR')
sys.exit()
if np.max(party) > nhocy*dhocy:
print('ERROR')
sys.exit()
npart = partx.size
hoclist= np.zeros(npart, dtype=np.int32)-1
hoc = np.zeros([nhocy, nhocx], dtype=np.int32)-1
for ipart in range(npart):
ix = int(partx[ipart]/dhocx)
iy = int(party[ipart]/dhocy)
hoclist[ipart] = hoc[iy, ix]
hoc[iy, ix] = ipart
return hoc, hoclist
def hocFind(px, py, dhocx, dhocy, hoc, hoclist):
ix = int(px/dhocx)
iy = int(py/dhocy)
neigh=[]
it = hoc[iy, ix]
while it != -1:
neigh.append(it)
it = hoclist[it]
return neigh
def findNeighbors_hoclist(px, py, tx=None,ty=None, dn=4, hoc=None, hoclist=None):
nhocy = nhocx = 20
pxMin = np.min(px)
pxMax = np.max(px)
pyMin = np.min(py)
pyMax = np.max(py)
dhocx = (pxMax - pxMin)/(nhocx-1)
dhocy = (pyMax - pyMin)/(nhocy-1)
partx = px - pxMin +dhocx/2
party = py - pyMin +dhocy/2
if hoc is None:
hoc, hoclist = hocBuild(partx, party, nhocx, nhocy, dhocx, dhocy)
return hoc, hoclist
if hoc is not None:
tx = tx - pxMin +dhocx/2
ty = ty - pyMin +dhocy/2
itx = int(tx/dhocx)
ity = int(ty/dhocy)
ps = [-1, 0, 1]
neigh=[]
for ii in range(3):
for jj in range(3):
ix = itx + ps[ii]
iy = ity + ps[jj]
if ix < 0:
continue
if iy < 0:
continue
if ix > nhocx-1:
continue
if iy > nhocy-1:
continue
#neightt = myUtil.hocFind(ppx, ppy, dhocx, dhocy, hoc, hoclist)
it = hoc[iy, ix]
while it != -1:
neigh.append(it)
it = hoclist[it]
#neigh.append(neightt)
#ll = [i for k in neigh for i in k]
if dn != -1:
ptx = np.array(partx[neigh])
pty = np.array(party[neigh])
dd = np.hypot(ptx-tx, pty-ty)
idx = np.argsort(dd)
neigh= np.array(neigh)[idx[0:dn]]
return neigh
###PSF中心对齐###
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
###插值对齐-fft###
def psfCentering_FFT(image):
"""
centering psf within an aperture by FFT
"""
ny, nx = image.shape
py, px = ndimage.center_of_mass(image)
dx = (px - nx/2)
dy = (py - ny/2)
k=np.zeros((nx,ny,2),dtype=float)
fg=np.fft.fft2(image)
ge =np.zeros_like(fg)
inx = int(nx/2)
jny = int(ny/2)
#prepare for the phase multiply matrix
#left bottom
for i in range(inx+1):
for j in range(jny+1):
k[i][j][0]=i;
k[i][j][1]=j;
#top right
for i in range(inx-1):
for j in range(jny-1):
k[i+inx+1][j+jny+1][0]=(-(nx/2-1)+i)
k[i+inx+1][j+jny+1][1]=(-(ny/2-1)+j)
#bottom right
for i in range(inx+1):
for j in range(jny-1):
k[i][j+jny+1][0]=i
k[i][j+jny+1][1]=(-(ny/2-1)+j)
#top left
for i in range(inx-1):
for j in range(int(jny+1)):
k[i+inx+1][j][0]=(-(nx/2-1)+i)
k[i+inx+1][j][1]=j
for i in range(nx):
for j in range(ny):
ge[i][j]=fg[i][j]*np.exp(2.*np.pi*(dx*k[i][j][0]/nx+dy*k[i][j][1]/ny)*1j)
get=np.fft.ifft2(ge).real
return(get)
###图像叠加###
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
###计算PSF椭率-接口###
def psfSizeCalculator(psfMat, psfSampleSize=2.5, CalcPSFcenter=True, SigRange=True, TailorScheme=2, cenPix=None):
"""
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, cenPix=cenPix)
#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
###计算PSF椭率###
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
"""
apr = 0.5 #arcsec, 0.5角秒内测量
fl = 28. #meters
pxs = 5.0 #microns
apr = np.deg2rad(apr/3600.)*fl*1e6
apr = apr/pxs
apr = np.int(np.ceil(apr))
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
rr = np.sqrt(x*x + y*y)
wgt= 0.0
if rr <= apr:
wgt = 1.0
w += I[jrow, icol]*wgt
w11 += x*x*I[jrow, icol]*wgt
w12 += x*y*I[jrow, icol]*wgt
w22 += y*y*I[jrow, icol]*wgt
w11 /= w
w12 /= w
w22 /= w
sz = w11 + w22
e1 = (w11 - w22)/sz
e2 = 2.0*w12/sz
return sz, e1, e2
###计算REE80###
def psfEncircle(img, fraction=0.8, psfSampleSizeInMicrons=2.5, focalLengthInMeters=28, cenPix=None):
"""
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)
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]
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 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 psfTailor(img, apSizeInArcsec=0.5, psfSampleSizeInMicrons=5, focalLengthInMeters=28, cenPix=None):
"""
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)
y,x = ndimage.center_of_mass(img) #y-rows, x-cols
imgMaxPix_x = int(x)
imgMaxPix_y = int(y)
if cenPix != None:
imgMaxPix_x = int(cenPix[0])
imgMaxPix_y = int(cenPix[1])
apSizeInMicrons = np.deg2rad(apSizeInArcsec/3600.)*focalLengthInMeters*1e6
apSizeInPix = apSizeInMicrons/psfSampleSizeInMicrons
apSizeInPix = np.int(np.ceil(apSizeInPix))
print('apSizeInPix=',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
###centroid with a window###
def centroidWgt(img, nt=160):
#libCentroid = ctypes.CDLL('/public/home/weichengliang/lnData/CSST_new_framwork/csstPSF/libCentroid/libCentroid.so') # CDLL加载库
libCentroid = ctypes.CDLL('/public/home/weichengliang/lnData/CSST_new_framwork/csstPSF_20210108/libCentroid/libCentroid.so') # CDLL加载库
libCentroid.centroidWgt.argtypes = [ctypes.POINTER(ctypes.c_float), ctypes.c_int, ctypes.c_int, ctypes.POINTER(ctypes.c_double)]
libCentroid.imSplint.argtypes = [ctypes.POINTER(ctypes.c_float), ctypes.c_int, ctypes.c_int, ctypes.POINTER(ctypes.c_double), ctypes.c_int, ctypes.c_int, ctypes.POINTER(ctypes.c_float)]
imx = img/np.sum(img)
ny, nx = imx.shape
#imx centroid
nn = nx*ny
arr = (ctypes.c_float*nn)()
arr[:] = imx.reshape(nn)
para = (ctypes.c_double*10)()
libCentroid.centroidWgt(arr, ny, nx, para)
imx_cy = para[3] #irow
imx_cx = para[4] #icol
#imx -> imy
nxt=nyt=nt
nn=nxt*nyt
yat = (ctypes.c_float*nn)()
libCentroid.imSplint(arr, ny, nx, para, nxt, nyt, yat)
imy = np.array(yat[:]).reshape([nxt, nyt])
return imy, imx_cx, imx_cy
'''
def psfCentering_wgt(ipsfMat, psf_image_x, psf_image_y, psfSampleSizeInMicrons=5.0):
img, cx, cy = centroidWgt(ipsfMat, nt=160)
nrows, ncols = ipsfMat.shape
cyt = (cy + nrows/2)*psfSampleSizeInMicrons*1e-3 +psf_image_y
cxt = (cx + ncols/2)*psfSampleSizeInMicrons*1e-3 +psf_image_x
return img, cxt, cyt
'''
###binning image###
def binningPSF(img, ngg):
imgX = img.reshape(ngg, img.shape[0]//ngg, ngg, img.shape[1]//ngg).mean(-1).mean(1)
return imgX
from .PSFInterp import PSFInterp
\ No newline at end of file
#OPTS += -D
CC = gcc
OPTIMIZE = -fPIC -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= -shared -lm #$(EXTRACLIB)
OBJS = centroidWgt.o nrutil.o
EXEC = libCentroid.so
all: $(EXEC)
$(EXEC): $(OBJS)
$(CLINK) $(CFLAGS) -o $@ $(OBJS) $(CLIB)
$(OBJS): nrutil.h Makefile
.PHONY : clean
clean:
rm -f *.o $(EXEC)
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "nrutil.h"
//y is the input image, nx and ny are the pixels along x and y axis;
//yt and residu are the required matrix for the function
//para contains the needed parameters: centx=para[3], centy=para[4]
void star_gaus(float **y,int nx,int ny,double *para);
void Interp_bicubic(int nx,int ny,float **a0,int nbound, int nxt,int nyt, float **at,double *xcen);
void splie2(float **ya, int m, int n, float **y2a);
void spline(float y[], int n, float yp1, float ypn, float y2[]);
void splint(float ya[], float y2a[], int n, float x, float *y);
void centroidWgt(float *y, int nx, int ny, double *para)
{
int i, j, k;
float **yy;
double **basef;
double **coff;
double ccol, crow;
yy = matrix(0, nx-1, 0, ny-1);
for(i=0; i<nx; i++)
for(j=0; j<ny; j++)
{
k = j + i*nx;
yy[i][j] = y[k];
}
star_gaus(yy, nx, ny, para);
//ccol = para[4];
//crow = para[3];
}
void imSplint(float *y, int nx, int ny, double *para, int nxt, int nyt, float *yat)
{
int i, j, k;
float **yy;
double **basef;
double **coff;
double ccol, crow;
yy = matrix(0, nx-1, 0, ny-1);
for(i=0; i<nx; i++)
for(j=0; j<ny; j++)
{
k = j + i*nx;
yy[i][j] = y[k];
}
crow = para[3];
ccol = para[4];
int nbound;
float **at;
double xcen[2];
at = matrix(0, nxt-1, 0, nyt-1);
nbound = (int)((nx - nxt)/2);
xcen[0]= crow;
xcen[1]= ccol;
Interp_bicubic(nx, ny, yy, nbound, nxt, nyt, at, xcen);
for(i=0; i<nxt; i++)
for(j=0; j<nyt; j++)
{
k = j + i*nxt;
yat[k] = at[i][j];
}
}
//void star_gaus(float **y,int nx,int ny,float **yt,float **residu,double *para)
void star_gaus(float **y,int nx,int ny,double *para)
{
void gasfit_2D(double **x, double *y,int np,double *para,double bg0, int fbg);
double **xs,*ys,ymax;
int i,j,np,npt,im,jm,**xi,k;
double kc,bgc,sigmac,xc,yc,sigma2v,det;
double bg0=0;
int fbg=0;
np=nx*ny;
xs=dmatrix(0,np-1,0,1);
ys=dvector(0,np-1);
xi=imatrix(0,np-1,0,1);
ymax=y[0][0];
for(i=0;i<nx;i++)for(j=0;j<ny;j++){
if(ymax<y[i][j]){ymax=y[i][j];im=i;jm=j;}
}
int cutPix;
cutPix = 23;
npt=0;
for(i=-cutPix;i<=cutPix;i++){
for(j=-cutPix;j<=cutPix;j++){
xs[npt][0]=xi[npt][0]=i+im;
xs[npt][1]=xi[npt][1]=j+jm;
ys[npt]=y[im+i][jm+j];
npt++;
}
}
gasfit_2D(xs, ys,npt,para,bg0,fbg);
kc=para[0];sigmac=para[1];bgc=para[2];xc=para[3];yc=para[4];
// printf("%e %e %e %e %e\n",kc,sigmac,bgc,xc,yc);
/*
sigma2v=-1./(2.*sigmac*sigmac);
for(i=0;i<nx;i++){
for(j=0;j<ny;j++){
det=DSQR(i-xc)+DSQR(j-yc);
yt[i][j]=kc*exp(det*sigma2v)+bgc;
residu[i][j]=y[i][j]-yt[i][j];
}
}
*/
free_dmatrix(xs,0,np-1,0,1);
free_imatrix(xi,0,np-1,0,1);
free_dvector(ys,0,np-1);
}
void gasfit_2D(double **x, double *y,int np,double *para,double bg0,int fbg)
{
void search_2D(double **x,double *y,int np,double kc,double kd,double sigmac,
double sigmad,double bgc,double bgd,double xc, double xd,
double yc, double yd,double *para,double bg0,int fbg);
int i,j,k,imax=0,isigma;
double ymax,ymin,kc,kd,sigmac,sigmad,bgc,bgd,ysigma,xc,yc,xd,yd;
double det,dett;
ymin=ymax=y[imax];
for(i=1;i<np;i++){
if(ymax<y[i]){imax=i;ymax=y[i];}
if(ymin>y[i])ymin=y[i];
}
kc=ymax;kd=ymax/12.;
det=ysigma=kc*exp(-0.5);
for(i=0;i<np;i++){
dett=fabs(ysigma-y[i]);
if(dett<det && y[i]!=kc){det=dett;isigma=i;}
//if(dett<det){det=dett;isigma=i;}
}
xc=x[imax][0];yc=x[imax][1];
sigmac=sqrt(DSQR(xc-x[isigma][0])+DSQR(yc-x[isigma][1]))*1.;
xd=yd=sigmac*0.25;
sigmad=0.25*sigmac;
bgc=0.;bgd=fabs(ymin);
for(i=0;i<4;i++){
search_2D(x,y,np,kc,kd,sigmac,sigmad,bgc,bgd,xc,xd,yc,yd,para,bg0,fbg);
kd*=0.33;sigmad*=0.33;bgd*=0.33;xd*=0.33;yd*=0.33;
kc=para[0];sigmac=para[1];bgc=para[2];xc=para[3];yc=para[4];
}
if(fbg==0)para[2]=bg0;
}
void search_2D(double **x,double *y,int np,double kc,double kd,double sigmac,
double sigmad,double bgc,double bgd,double xc, double xd,
double yc, double yd,double *para,double bg0,int fbg)
{
double k,sigma,bg,k2,k20,sigma2v,det,xt,yt;
int i,j,l,m,p,q;
sigma2v=-1./(2.*sigmac*sigmac);
bg=bgc;
k20=0;
for(m=0;m<np;m++){
det=DSQR(x[m][0]-xc)+DSQR(x[m][1]-yc);
det=kc*exp(det*sigma2v)+bgc-y[m];
k20+=det*det;
}
for(i=-4;i<=4;i++){
k=kc+i*kd;
if(k>0){
for(j=-4;j<=4;j++){
sigma=sigmac+j*sigmad;
if(sigma>0){
sigma2v=-1./(2.*sigma*sigma);
for(p=-4;p<=4;p++){
xt=xc+p*xd;
for(q=-4;q<=4;q++){
yt=yc+q*yd;
k2=0;
if(fbg==0){
bg=bg0;
for(m=0;m<np;m++){
det=DSQR(x[m][0]-xt)+DSQR(x[m][1]-yt);
det=k*exp(det*sigma2v)+bg-y[m];
k2+=det*det;
}
}
else{
for(l=-4;l<=4;l++){
bg=bgc+l*bgd;
for(m=0;m<np;m++){
det=DSQR(x[m][0]-xt)+DSQR(x[m][1]-yt);
det=k*exp(det*sigma2v)+bg-y[m];
k2+=det*det;
}
}
}
//printf("k20=%e k2=%e\n",k20,k2);
if(k2<=k20){k20=k2;para[5]=k2;
para[0]=k;para[1]=sigma;para[2]=bg;para[3]=xt;para[4]=yt;}
}
}
}
}
}
}
}
/////////////////////////////
//#include <math.h>
//#include <stdio.h>
//#include <stdlib.h>
//#include "nrutil.h"
//nx,ny is the pixels of x and y direction
//a0 is the input image in nx*ny
//nbound is the boundary limit
//nxt,nyt is out put image dimentions
//at is the output image and cen represent the required center
void Interp_bicubic(int nx,int ny,float **a0,int nbound,
int nxt,int nyt, float **at,double *xcen)
{
void splie2(float **ya, int m, int n, float **y2a);
void spline(float y[], int n, float yp1, float ypn, float y2[]);
void splint(float ya[], float y2a[], int n, float x, float *y);
int i,j,m,n,ic,jc;
float *ytmp,*yytmp,**y2a,x1,x2,shift1,shift2;
y2a=matrix(0,nx,0,ny);
ytmp=vector(0,ny);
yytmp=vector(0,ny);
splie2(a0,nx,ny,y2a);
ic=nx*0.5;
jc=ny*0.5;
shift1=xcen[0]-ic;//-0.5;
shift2=xcen[1]-jc;//-0.5;
if(fabs(shift1)>nbound || fabs(shift2)>nbound){
printf("bicubic shifts too much %e %e\n",shift1,shift2);
getchar();
}
for (n=0;n<nyt;n++){
x2=n+nbound+shift2;
for(i=0;i<nx;i++){
splint(a0[i],y2a[i],ny,x2,&yytmp[i]);
spline(yytmp,ny,1.0e30,1.0e30,ytmp);
}
for (m=0;m<nxt;m++){
x1=m+nbound+shift1;
splint(yytmp,ytmp,ny,x1,&at[m][n]);
}
}
free_vector(yytmp,0,ny);
free_vector(ytmp,0,ny);
free_matrix(y2a,0,nx,0,ny);
}
void splie2(float **ya, int m, int n, float **y2a)
{
void spline(float y[], int n, float yp1, float ypn, float y2[]);
int j;
for (j=0;j<m;j++)spline(ya[j],n,1.0e30,1.0e30,y2a[j]);
}
void spline(float y[], int n, float yp1, float ypn, float y2[])
{
int i,k;
float p,qn,sig,un,*u;
u=vector(0,n-1);
if (yp1 > 0.99e30)
y2[0]=u[0]=0.0;
else {
y2[0] = -0.5;
u[0]=3.0*(y[1]-y[0]-yp1);
}
sig=0.5;
for (i=1;i<=n-2;i++) {
p=sig*y2[i-1]+2.0;
y2[i]=(sig-1.0)/p;
u[i]=(y[i+1]-2.*y[i]+y[i-1]);
u[i]=(3.0*u[i]-sig*u[i-1])/p;
}
if (ypn > 0.99e30)
qn=un=0.0;
else {
qn=0.5;
un=3.0*(ypn-y[n-1]+y[n-2]);
}
y2[n-1]=(un-qn*u[n-2])/(qn*y2[n-2]+1.0);
for (k=n-2;k>=0;k--)y2[k]=y2[k]*y2[k+1]+u[k];
free_vector(u,0,n-1);
}
void splint(float ya[], float y2a[], int n, float x, float *y)
{
void nrerror(char error_text[]);
int klo,khi,k;
float h,b,a;
klo=x;
if( klo<0)klo=0;
if(klo==n-1)klo=n-2;
khi=klo+1;
if (klo<0 || khi>=n) printf("error:klo, khi, n: %d %d %d\n", klo, khi, n);
if (klo<0 || khi>=n) nrerror("Bad xa input to routine splint");
a=(khi-x);
b=(x-klo);
*y=a*ya[klo]+b*ya[khi]+((a*a*a-a)*y2a[klo]+(b*b*b-b)*y2a[khi])/6.0;
}
#if defined(__STDC__) || defined(ANSI) || defined(NRANSI) /* ANSI */
#include <stdio.h>
#include <stddef.h>
#include <stdlib.h>
#define NR_END 1
#define FREE_ARG char*
void nrerror(char error_text[])
/* Numerical Recipes standard error handler */
{
fprintf(stderr,"Numerical Recipes run-time error...\n");
fprintf(stderr,"%s\n",error_text);
fprintf(stderr,"...now exiting to system...\n");
exit(1);
}
float *vector(long nl, long nh)
/* allocate a float vector with subscript range v[nl..nh] */
{
float *v;
v=(float *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(float)));
if (!v) nrerror("allocation failure in vector()");
return v-nl+NR_END;
}
int *ivector(long nl, long nh)
/* allocate an int vector with subscript range v[nl..nh] */
{
int *v;
v=(int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int)));
if (!v) nrerror("allocation failure in ivector()");
return v-nl+NR_END;
}
unsigned char *cvector(long nl, long nh)
/* allocate an unsigned char vector with subscript range v[nl..nh] */
{
unsigned char *v;
v=(unsigned char *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(unsigned char)));
if (!v) nrerror("allocation failure in cvector()");
return v-nl+NR_END;
}
unsigned long *lvector(long nl, long nh)
/* allocate an unsigned long vector with subscript range v[nl..nh] */
{
unsigned long *v;
v=(unsigned long *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(long)));
if (!v) nrerror("allocation failure in lvector()");
return v-nl+NR_END;
}
double *dvector(long nl, long nh)
/* allocate a double vector with subscript range v[nl..nh] */
{
double *v;
v=(double *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(double)));
if (!v) nrerror("allocation failure in dvector()");
return v-nl+NR_END;
}
float **matrix(long nrl, long nrh, long ncl, long nch)
/* allocate a float matrix with subscript range m[nrl..nrh][ncl..nch] */
{
long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
float **m;
/* allocate pointers to rows */
m=(float **) malloc((size_t)((nrow+NR_END)*sizeof(float*)));
if (!m) nrerror("allocation failure 1 in matrix()");
m += NR_END;
m -= nrl;
/* allocate rows and set pointers to them */
m[nrl]=(float *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float)));
if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
m[nrl] += NR_END;
m[nrl] -= ncl;
for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
/* return pointer to array of pointers to rows */
return m;
}
double **dmatrix(long nrl, long nrh, long ncl, long nch)
/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
{
long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
double **m;
/* allocate pointers to rows */
m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*)));
if (!m) nrerror("allocation failure 1 in matrix()");
m += NR_END;
m -= nrl;
/* allocate rows and set pointers to them */
m[nrl]=(double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double)));
if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
m[nrl] += NR_END;
m[nrl] -= ncl;
for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
/* return pointer to array of pointers to rows */
return m;
}
int **imatrix(long nrl, long nrh, long ncl, long nch)
/* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */
{
long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
int **m;
/* allocate pointers to rows */
m=(int **) malloc((size_t)((nrow+NR_END)*sizeof(int*)));
if (!m) nrerror("allocation failure 1 in matrix()");
m += NR_END;
m -= nrl;
/* allocate rows and set pointers to them */
m[nrl]=(int *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(int)));
if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
m[nrl] += NR_END;
m[nrl] -= ncl;
for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
/* return pointer to array of pointers to rows */
return m;
}
float **submatrix(float **a, long oldrl, long oldrh, long oldcl, long oldch,
long newrl, long newcl)
/* point a submatrix [newrl..][newcl..] to a[oldrl..oldrh][oldcl..oldch] */
{
long i,j,nrow=oldrh-oldrl+1,ncol=oldcl-newcl;
float **m;
/* allocate array of pointers to rows */
m=(float **) malloc((size_t) ((nrow+NR_END)*sizeof(float*)));
if (!m) nrerror("allocation failure in submatrix()");
m += NR_END;
m -= newrl;
/* set pointers to rows */
for(i=oldrl,j=newrl;i<=oldrh;i++,j++) m[j]=a[i]+ncol;
/* return pointer to array of pointers to rows */
return m;
}
float **convert_matrix(float *a, long nrl, long nrh, long ncl, long nch)
/* allocate a float matrix m[nrl..nrh][ncl..nch] that points to the matrix
declared in the standard C manner as a[nrow][ncol], where nrow=nrh-nrl+1
and ncol=nch-ncl+1. The routine should be called with the address
&a[0][0] as the first argument. */
{
long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1;
float **m;
/* allocate pointers to rows */
m=(float **) malloc((size_t) ((nrow+NR_END)*sizeof(float*)));
if (!m) nrerror("allocation failure in convert_matrix()");
m += NR_END;
m -= nrl;
/* set pointers to rows */
m[nrl]=a-ncl;
for(i=1,j=nrl+1;i<nrow;i++,j++) m[j]=m[j-1]+ncol;
/* return pointer to array of pointers to rows */
return m;
}
float ***f3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
/* allocate a float 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] */
{
long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
float ***t;
/* allocate pointers to pointers to rows */
t=(float ***) malloc((size_t)((nrow+NR_END)*sizeof(float**)));
if (!t) nrerror("allocation failure 1 in f3tensor()");
t += NR_END;
t -= nrl;
/* allocate pointers to rows and set pointers to them */
t[nrl]=(float **) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float*)));
if (!t[nrl]) nrerror("allocation failure 2 in f3tensor()");
t[nrl] += NR_END;
t[nrl] -= ncl;
/* allocate rows and set pointers to them */
t[nrl][ncl]=(float *) malloc((size_t)((nrow*ncol*ndep+NR_END)*sizeof(float)));
if (!t[nrl][ncl]) nrerror("allocation failure 3 in f3tensor()");
t[nrl][ncl] += NR_END;
t[nrl][ncl] -= ndl;
for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep;
for(i=nrl+1;i<=nrh;i++) {
t[i]=t[i-1]+ncol;
t[i][ncl]=t[i-1][ncl]+ncol*ndep;
for(j=ncl+1;j<=nch;j++) t[i][j]=t[i][j-1]+ndep;
}
/* return pointer to array of pointers to rows */
return t;
}
double ***d3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
/* allocate a double 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] */
{
long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
double ***t;
/* allocate pointers to pointers to rows */
t=(double ***) malloc((size_t)((nrow+NR_END)*sizeof(double**)));
if (!t) nrerror("allocation failure 1 in f3tensor()");
t += NR_END;
t -= nrl;
/* allocate pointers to rows and set pointers to them */
t[nrl]=(double **) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double*)));
if (!t[nrl]) nrerror("allocation failure 2 in f3tensor()");
t[nrl] += NR_END;
t[nrl] -= ncl;
/* allocate rows and set pointers to them */
t[nrl][ncl]=(double *) malloc((size_t)((nrow*ncol*ndep+NR_END)*sizeof(double)));
if (!t[nrl][ncl]) nrerror("allocation failure 3 in f3tensor()");
t[nrl][ncl] += NR_END;
t[nrl][ncl] -= ndl;
for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep;
for(i=nrl+1;i<=nrh;i++) {
t[i]=t[i-1]+ncol;
t[i][ncl]=t[i-1][ncl]+ncol*ndep;
for(j=ncl+1;j<=nch;j++) t[i][j]=t[i][j-1]+ndep;
}
/* return pointer to array of pointers to rows */
return t;
}
unsigned char **cmatrix(long nrl, long nrh, long ncl, long nch)
/* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */
{
long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
unsigned char **m;
/* allocate pointers to rows */
m=(unsigned char **) malloc((size_t)((nrow+NR_END)*sizeof(unsigned char*)));
if (!m) nrerror("allocation failure 1 in matrix()");
m += NR_END;
m -= nrl;
/* allocate rows and set pointers to them */
m[nrl]=(unsigned char *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(unsigned char)));
if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
m[nrl] += NR_END;
m[nrl] -= ncl;
for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
/* return pointer to array of pointers to rows */
return m;
}
void free_vector(float *v, long nl, long nh)
/* free a float vector allocated with vector() */
{
free((FREE_ARG) (v+nl-NR_END));
}
void free_ivector(int *v, long nl, long nh)
/* free an int vector allocated with ivector() */
{
free((FREE_ARG) (v+nl-NR_END));
}
void free_cvector(unsigned char *v, long nl, long nh)
/* free an unsigned char vector allocated with cvector() */
{
free((FREE_ARG) (v+nl-NR_END));
}
void free_lvector(unsigned long *v, long nl, long nh)
/* free an unsigned long vector allocated with lvector() */
{
free((FREE_ARG) (v+nl-NR_END));
}
void free_dvector(double *v, long nl, long nh)
/* free a double vector allocated with dvector() */
{
free((FREE_ARG) (v+nl-NR_END));
}
void free_matrix(float **m, long nrl, long nrh, long ncl, long nch)
/* free a float matrix allocated by matrix() */
{
free((FREE_ARG) (m[nrl]+ncl-NR_END));
free((FREE_ARG) (m+nrl-NR_END));
}
void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch)
/* free a double matrix allocated by dmatrix() */
{
free((FREE_ARG) (m[nrl]+ncl-NR_END));
free((FREE_ARG) (m+nrl-NR_END));
}
void free_imatrix(int **m, long nrl, long nrh, long ncl, long nch)
/* free an int matrix allocated by imatrix() */
{
free((FREE_ARG) (m[nrl]+ncl-NR_END));
free((FREE_ARG) (m+nrl-NR_END));
}
void free_submatrix(float **b, long nrl, long nrh, long ncl, long nch)
/* free a submatrix allocated by submatrix() */
{
free((FREE_ARG) (b+nrl-NR_END));
}
void free_convert_matrix(float **b, long nrl, long nrh, long ncl, long nch)
/* free a matrix allocated by convert_matrix() */
{
free((FREE_ARG) (b+nrl-NR_END));
}
void free_f3tensor(float ***t, long nrl, long nrh, long ncl, long nch,
long ndl, long ndh)
/* free a float f3tensor allocated by f3tensor() */
{
free((FREE_ARG) (t[nrl][ncl]+ndl-NR_END));
free((FREE_ARG) (t[nrl]+ncl-NR_END));
free((FREE_ARG) (t+nrl-NR_END));
}
void free_d3tensor(double ***t, long nrl, long nrh, long ncl, long nch,
long ndl, long ndh)
/* free a double f3tensor allocated by f3tensor() */
{
free((FREE_ARG) (t[nrl][ncl]+ndl-NR_END));
free((FREE_ARG) (t[nrl]+ncl-NR_END));
free((FREE_ARG) (t+nrl-NR_END));
}
void free_cmatrix(unsigned char **m, long nrl, long nrh, long ncl, long nch)
/* free a character matrix allocated by matrix() */
{
free((FREE_ARG) (m[nrl]+ncl-NR_END));
free((FREE_ARG) (m+nrl-NR_END));
}
#else /* ANSI */
/* traditional - K&R */
#include <stdio.h>
#define NR_END 1
#define FREE_ARG char*
void nrerror(error_text)
char error_text[];
/* Numerical Recipes standard error handler */
{
void exit();
fprintf(stderr,"Numerical Recipes run-time error...\n");
fprintf(stderr,"%s\n",error_text);
fprintf(stderr,"...now exiting to system...\n");
exit(1);
}
float *vector(nl,nh)
long nh,nl;
/* allocate a float vector with subscript range v[nl..nh] */
{
float *v;
v=(float *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(float)));
if (!v) nrerror("allocation failure in vector()");
return v-nl+NR_END;
}
int *ivector(nl,nh)
long nh,nl;
/* allocate an int vector with subscript range v[nl..nh] */
{
int *v;
v=(int *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(int)));
if (!v) nrerror("allocation failure in ivector()");
return v-nl+NR_END;
}
unsigned char *cvector(nl,nh)
long nh,nl;
/* allocate an unsigned char vector with subscript range v[nl..nh] */
{
unsigned char *v;
v=(unsigned char *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(unsigned char)));
if (!v) nrerror("allocation failure in cvector()");
return v-nl+NR_END;
}
unsigned long *lvector(nl,nh)
long nh,nl;
/* allocate an unsigned long vector with subscript range v[nl..nh] */
{
unsigned long *v;
v=(unsigned long *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(long)));
if (!v) nrerror("allocation failure in lvector()");
return v-nl+NR_END;
}
double *dvector(nl,nh)
long nh,nl;
/* allocate a double vector with subscript range v[nl..nh] */
{
double *v;
v=(double *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(double)));
if (!v) nrerror("allocation failure in dvector()");
return v-nl+NR_END;
}
float **matrix(nrl,nrh,ncl,nch)
long nch,ncl,nrh,nrl;
/* allocate a float matrix with subscript range m[nrl..nrh][ncl..nch] */
{
long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
float **m;
/* allocate pointers to rows */
m=(float **) malloc((unsigned int)((nrow+NR_END)*sizeof(float*)));
if (!m) nrerror("allocation failure 1 in matrix()");
m += NR_END;
m -= nrl;
/* allocate rows and set pointers to them */
m[nrl]=(float *) malloc((unsigned int)((nrow*ncol+NR_END)*sizeof(float)));
if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
m[nrl] += NR_END;
m[nrl] -= ncl;
for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
/* return pointer to array of pointers to rows */
return m;
}
double **dmatrix(nrl,nrh,ncl,nch)
long nch,ncl,nrh,nrl;
/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
{
long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
double **m;
/* allocate pointers to rows */
m=(double **) malloc((unsigned int)((nrow+NR_END)*sizeof(double*)));
if (!m) nrerror("allocation failure 1 in matrix()");
m += NR_END;
m -= nrl;
/* allocate rows and set pointers to them */
m[nrl]=(double *) malloc((unsigned int)((nrow*ncol+NR_END)*sizeof(double)));
if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
m[nrl] += NR_END;
m[nrl] -= ncl;
for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
/* return pointer to array of pointers to rows */
return m;
}
int **imatrix(nrl,nrh,ncl,nch)
long nch,ncl,nrh,nrl;
/* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */
{
long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
int **m;
/* allocate pointers to rows */
m=(int **) malloc((unsigned int)((nrow+NR_END)*sizeof(int*)));
if (!m) nrerror("allocation failure 1 in matrix()");
m += NR_END;
m -= nrl;
/* allocate rows and set pointers to them */
m[nrl]=(int *) malloc((unsigned int)((nrow*ncol+NR_END)*sizeof(int)));
if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
m[nrl] += NR_END;
m[nrl] -= ncl;
for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
/* return pointer to array of pointers to rows */
return m;
}
float **submatrix(a,oldrl,oldrh,oldcl,oldch,newrl,newcl)
float **a;
long newcl,newrl,oldch,oldcl,oldrh,oldrl;
/* point a submatrix [newrl..][newcl..] to a[oldrl..oldrh][oldcl..oldch] */
{
long i,j,nrow=oldrh-oldrl+1,ncol=oldcl-newcl;
float **m;
/* allocate array of pointers to rows */
m=(float **) malloc((unsigned int) ((nrow+NR_END)*sizeof(float*)));
if (!m) nrerror("allocation failure in submatrix()");
m += NR_END;
m -= newrl;
/* set pointers to rows */
for(i=oldrl,j=newrl;i<=oldrh;i++,j++) m[j]=a[i]+ncol;
/* return pointer to array of pointers to rows */
return m;
}
float **convert_matrix(a,nrl,nrh,ncl,nch)
float *a;
long nch,ncl,nrh,nrl;
/* allocate a float matrix m[nrl..nrh][ncl..nch] that points to the matrix
declared in the standard C manner as a[nrow][ncol], where nrow=nrh-nrl+1
and ncol=nch-ncl+1. The routine should be called with the address
&a[0][0] as the first argument. */
{
long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1;
float **m;
/* allocate pointers to rows */
m=(float **) malloc((unsigned int) ((nrow+NR_END)*sizeof(float*)));
if (!m) nrerror("allocation failure in convert_matrix()");
m += NR_END;
m -= nrl;
/* set pointers to rows */
m[nrl]=a-ncl;
for(i=1,j=nrl+1;i<nrow;i++,j++) m[j]=m[j-1]+ncol;
/* return pointer to array of pointers to rows */
return m;
}
double ***d3tensor(nrl,nrh,ncl,nch,ndl,ndh)
long nch,ncl,ndh,ndl,nrh,nrl;
/* allocate a float 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] */
{
long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
double ***t;
/* allocate pointers to pointers to rows */
t=(double ***) malloc((unsigned int)((nrow+NR_END)*sizeof(double**)));
if (!t) nrerror("allocation failure 1 in f3tensor()");
t += NR_END;
t -= nrl;
/* allocate pointers to rows and set pointers to them */
t[nrl]=(double **) malloc((unsigned int)((nrow*ncol+NR_END)*sizeof(double*)));
if (!t[nrl]) nrerror("allocation failure 2 in f3tensor()");
t[nrl] += NR_END;
t[nrl] -= ncl;
/* allocate rows and set pointers to them */
t[nrl][ncl]=(double *) malloc((unsigned int)((nrow*ncol*ndep+NR_END)*sizeof(double)));
if (!t[nrl][ncl]) nrerror("allocation failure 3 in f3tensor()");
t[nrl][ncl] += NR_END;
t[nrl][ncl] -= ndl;
for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep;
for(i=nrl+1;i<=nrh;i++) {
t[i]=t[i-1]+ncol;
t[i][ncl]=t[i-1][ncl]+ncol*ndep;
for(j=ncl+1;j<=nch;j++) t[i][j]=t[i][j-1]+ndep;
}
/* return pointer to array of pointers to rows */
return t;
}
float ***f3tensor(nrl,nrh,ncl,nch,ndl,ndh)
long nch,ncl,ndh,ndl,nrh,nrl;
/* allocate a float 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] */
{
long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
float ***t;
/* allocate pointers to pointers to rows */
t=(float ***) malloc((unsigned int)((nrow+NR_END)*sizeof(float**)));
if (!t) nrerror("allocation failure 1 in f3tensor()");
t += NR_END;
t -= nrl;
/* allocate pointers to rows and set pointers to them */
t[nrl]=(float **) malloc((unsigned int)((nrow*ncol+NR_END)*sizeof(float*)));
if (!t[nrl]) nrerror("allocation failure 2 in f3tensor()");
t[nrl] += NR_END;
t[nrl] -= ncl;
/* allocate rows and set pointers to them */
t[nrl][ncl]=(float *) malloc((unsigned int)((nrow*ncol*ndep+NR_END)*sizeof(float)));
if (!t[nrl][ncl]) nrerror("allocation failure 3 in f3tensor()");
t[nrl][ncl] += NR_END;
t[nrl][ncl] -= ndl;
for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep;
for(i=nrl+1;i<=nrh;i++) {
t[i]=t[i-1]+ncol;
t[i][ncl]=t[i-1][ncl]+ncol*ndep;
for(j=ncl+1;j<=nch;j++) t[i][j]=t[i][j-1]+ndep;
}
/* return pointer to array of pointers to rows */
return t;
}
unsigned char **cmatrix(nrl,nrh,ncl,nch)
long nch,ncl,nrh,nrl;
/* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */
{
long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
unsigned char **m;
/* allocate pointers to rows */
m=(unsigned char **) malloc((unsigned int)((nrow+NR_END)*sizeof(unsigned char*)));
if (!m) nrerror("allocation failure 1 in matrix()");
m += NR_END;
m -= nrl;
/* allocate rows and set pointers to them */
m[nrl]=(unsigned char *) malloc((unsigned int)((nrow*ncol+NR_END)*sizeof(unsigned char)));
if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
m[nrl] += NR_END;
m[nrl] -= ncl;
for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
/* return pointer to array of pointers to rows */
return m;
}
void free_vector(v,nl,nh)
float *v;
long nh,nl;
/* free a float vector allocated with vector() */
{
free((FREE_ARG) (v+nl-NR_END));
}
void free_ivector(v,nl,nh)
int *v;
long nh,nl;
/* free an int vector allocated with ivector() */
{
free((FREE_ARG) (v+nl-NR_END));
}
void free_cvector(v,nl,nh)
long nh,nl;
unsigned char *v;
/* free an unsigned char vector allocated with cvector() */
{
free((FREE_ARG) (v+nl-NR_END));
}
void free_lvector(v,nl,nh)
long nh,nl;
unsigned long *v;
/* free an unsigned long vector allocated with lvector() */
{
free((FREE_ARG) (v+nl-NR_END));
}
void free_dvector(v,nl,nh)
double *v;
long nh,nl;
/* free a double vector allocated with dvector() */
{
free((FREE_ARG) (v+nl-NR_END));
}
void free_matrix(m,nrl,nrh,ncl,nch)
float **m;
long nch,ncl,nrh,nrl;
/* free a float matrix allocated by matrix() */
{
free((FREE_ARG) (m[nrl]+ncl-NR_END));
free((FREE_ARG) (m+nrl-NR_END));
}
void free_dmatrix(m,nrl,nrh,ncl,nch)
double **m;
long nch,ncl,nrh,nrl;
/* free a double matrix allocated by dmatrix() */
{
free((FREE_ARG) (m[nrl]+ncl-NR_END));
free((FREE_ARG) (m+nrl-NR_END));
}
void free_imatrix(m,nrl,nrh,ncl,nch)
int **m;
long nch,ncl,nrh,nrl;
/* free an int matrix allocated by imatrix() */
{
free((FREE_ARG) (m[nrl]+ncl-NR_END));
free((FREE_ARG) (m+nrl-NR_END));
}
void free_submatrix(b,nrl,nrh,ncl,nch)
float **b;
long nch,ncl,nrh,nrl;
/* free a submatrix allocated by submatrix() */
{
free((FREE_ARG) (b+nrl-NR_END));
}
void free_convert_matrix(b,nrl,nrh,ncl,nch)
float **b;
long nch,ncl,nrh,nrl;
/* free a matrix allocated by convert_matrix() */
{
free((FREE_ARG) (b+nrl-NR_END));
}
void free_f3tensor(t,nrl,nrh,ncl,nch,ndl,ndh)
float ***t;
long nch,ncl,ndh,ndl,nrh,nrl;
/* free a float f3tensor allocated by f3tensor() */
{
free((FREE_ARG) (t[nrl][ncl]+ndl-NR_END));
free((FREE_ARG) (t[nrl]+ncl-NR_END));
free((FREE_ARG) (t+nrl-NR_END));
}
void free_d3tensor(t,nrl,nrh,ncl,nch,ndl,ndh)
double ***t;
long nch,ncl,ndh,ndl,nrh,nrl;
/* free a float f3tensor allocated by f3tensor() */
{
free((FREE_ARG) (t[nrl][ncl]+ndl-NR_END));
free((FREE_ARG) (t[nrl]+ncl-NR_END));
free((FREE_ARG) (t+nrl-NR_END));
}
void free_cmatrix(m,nrl,nrh,ncl,nch)
unsigned **m;
long nch,ncl,nrh,nrl;
/* free a double matrix allocated by dmatrix() */
{
free((FREE_ARG) (m[nrl]+ncl-NR_END));
free((FREE_ARG) (m+nrl-NR_END));
}
#endif /* ANSI */
/* CAUTION: This is the ANSI C (only) version of the Numerical Recipes
utility file nrutil.h. Do not confuse this file with the same-named
file nrutil.h that may be supplied in a 'misc' subdirectory.
*That* file is the one from the book, and contains both ANSI and
traditional K&R versions, along with #ifdef macros to select the
correct version. *This* file contains only ANSI C. */
#ifndef _NR_UTILS_H_
#define _NR_UTILS_H_
static float sqrarg;
#define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)
static double dsqrarg;
#define DSQR(a) ((dsqrarg=(a)) == 0.0 ? 0.0 : dsqrarg*dsqrarg)
static double dmaxarg1,dmaxarg2;
#define DMAX(a,b) (dmaxarg1=(a),dmaxarg2=(b),(dmaxarg1) > (dmaxarg2) ?\
(dmaxarg1) : (dmaxarg2))
static double dminarg1,dminarg2;
#define DMIN(a,b) (dminarg1=(a),dminarg2=(b),(dminarg1) < (dminarg2) ?\
(dminarg1) : (dminarg2))
static float maxarg1,maxarg2;
#define FMAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) > (maxarg2) ?\
(maxarg1) : (maxarg2))
static float minarg1,minarg2;
#define FMIN(a,b) (minarg1=(a),minarg2=(b),(minarg1) < (minarg2) ?\
(minarg1) : (minarg2))
static long lmaxarg1,lmaxarg2;
#define LMAX(a,b) (lmaxarg1=(a),lmaxarg2=(b),(lmaxarg1) > (lmaxarg2) ?\
(lmaxarg1) : (lmaxarg2))
static long lminarg1,lminarg2;
#define LMIN(a,b) (lminarg1=(a),lminarg2=(b),(lminarg1) < (lminarg2) ?\
(lminarg1) : (lminarg2))
static int imaxarg1,imaxarg2;
#define IMAX(a,b) (imaxarg1=(a),imaxarg2=(b),(imaxarg1) > (imaxarg2) ?\
(imaxarg1) : (imaxarg2))
static int iminarg1,iminarg2;
#define IMIN(a,b) (iminarg1=(a),iminarg2=(b),(iminarg1) < (iminarg2) ?\
(iminarg1) : (iminarg2))
#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
void nrerror(char error_text[]);
float *vector(long nl, long nh);
int *ivector(long nl, long nh);
unsigned char *cvector(long nl, long nh);
unsigned long *lvector(long nl, long nh);
double *dvector(long nl, long nh);
float **matrix(long nrl, long nrh, long ncl, long nch);
double **dmatrix(long nrl, long nrh, long ncl, long nch);
int **imatrix(long nrl, long nrh, long ncl, long nch);
float **submatrix(float **a, long oldrl, long oldrh, long oldcl, long oldch,
long newrl, long newcl);
float **convert_matrix(float *a, long nrl, long nrh, long ncl, long nch);
float ***f3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh);
void free_vector(float *v, long nl, long nh);
void free_ivector(int *v, long nl, long nh);
void free_cvector(unsigned char *v, long nl, long nh);
void free_lvector(unsigned long *v, long nl, long nh);
void free_dvector(double *v, long nl, long nh);
void free_matrix(float **m, long nrl, long nrh, long ncl, long nch);
void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch);
void free_imatrix(int **m, long nrl, long nrh, long ncl, long nch);
void free_submatrix(float **b, long nrl, long nrh, long ncl, long nch);
void free_convert_matrix(float **b, long nrl, long nrh, long ncl, long nch);
void free_f3tensor(float ***t, long nrl, long nrh, long ncl, long nch,
long ndl, long ndh);
int ***i3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh);
void free_i3tensor(int ***t, long nrl, long nrh, long ncl, long nch,
long ndl, long ndh);
unsigned char ***b3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh);
void free_b3tensor(unsigned char ***t, long nrl, long nrh, long ncl, long nch,
long ndl, long ndh);
double ***d3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh);
void free_d3tensor(double ***t, long nrl, long nrh, long ncl, long nch,
long ndl, long ndh);
unsigned char **cmatrix(long nrl, long nrh, long ncl, long nch);
void free_cmatrix(unsigned char **m, long nrl, long nrh, long ncl, long nch);
#endif /* _NR_UTILS_H_ */
import numpy as np
import matplotlib.pyplot as plt
import ctypes
import galsim
libCentroid = ctypes.CDLL('../libCentroid.so') # CDLL加载库
print('load libCenroid')
libCentroid.centroidWgt.argtypes = [ctypes.POINTER(ctypes.c_float), ctypes.c_int, ctypes.c_int, ctypes.POINTER(ctypes.c_double)]
IMAGE_WIDTH = 180
IMAGE_HEIGHT = 180
dx = 0.275
dy =-0.393
center_x = 90+dx
center_y = 90+dy
R = np.sqrt(center_x**2 + center_y**2)
Gauss_map = np.zeros((IMAGE_HEIGHT, IMAGE_WIDTH))
for i in range(IMAGE_HEIGHT):
for j in range(IMAGE_WIDTH):
dis = (i-center_y)**2+(j-center_x)**2
Gauss_map[i, j] = np.exp(-0.5*dis/R)
ymap = galsim.InterpolatedImage(galsim.ImageF(Gauss_map), scale=0.01)
zmap = ymap#.shear(g1 = 0.15, g2=0.27) #using shear
Gauss_map = zmap.drawImage(nx = 180, ny = 180, scale=0.01).array
fig=plt.figure(figsize=(5,5))
plt.imshow(Gauss_map, origin='lower')
imx = Gauss_map/np.sum(Gauss_map)
ny,nx = imx.shape
print(nx, ny, np.max(imx))
nn = nx*ny
arr = (ctypes.c_float*nn)()
arr[:] = imx.reshape(nn)
para = (ctypes.c_double*10)()
libCentroid.centroidWgt(arr, ny, nx, para)
print('haha')
print(para[0:5])
cx = para[3]
cy = para[4]
print('{:}'.format(cx), '{:}'.format(cy))
from scipy import ndimage
cx, cy = ndimage.center_of_mass(imx)
print('center_of_mass:', cx, cy)
plt.plot(para[4], para[3], 'bx', ms = 15)
plt.plot(cy, cx, 'r+', ms = 15)
plt.annotate('dx,dy:', [10, 170], color='w')
plt.annotate('{:8.5}, {:8.5}'.format(dx, dy), [10, 160], color='w')
plt.annotate('cx, cy:', [10, 150], color='w')
plt.annotate('{:0<8.5}, {:0<8.5}'.format(center_x, center_y), [10, 140], color='w')
plt.annotate('{:0<8.5}, {:0<8.5}'.format(para[4], para[3]), [10, 130], color='w')
plt.annotate('{:0<8.5}, {:0<8.5}'.format(cy, cx), [10, 120], color='w')
plt.show()
import numpy as np
import matplotlib.pyplot as plt
import ctypes
import galsim
libCentroid = ctypes.CDLL('../libCentroid.so') # CDLL加载库
print('load libCenroid')
libCentroid.centroidWgt.argtypes = [ctypes.POINTER(ctypes.c_float), ctypes.c_int, ctypes.c_int, ctypes.POINTER(ctypes.c_double)]
libCentroid.imSplint.argtypes = [ctypes.POINTER(ctypes.c_float), ctypes.c_int, ctypes.c_int, ctypes.POINTER(ctypes.c_double), ctypes.c_int, ctypes.c_int, ctypes.POINTER(ctypes.c_float)]
from scipy.io import loadmat
data = loadmat('/Users/chengliangwei/csstPSFdata/CSSOS_psf_ciomp/ccd20/wave_1/5_psf_array/psf_10.mat')
imx = data['psf']
imx = imx/np.sum(imx)
ny,nx = imx.shape
print(nx, ny)
#imx centroid
nn = nx*ny
arr = (ctypes.c_float*nn)()
arr[:] = imx.reshape(nn)
para = (ctypes.c_double*10)()
libCentroid.centroidWgt(arr, ny, nx, para)
imx_cx = para[3]
imx_cy = para[4]
#imx -> imy
nxt=nyt=160
nn=nxt*nyt
yat = (ctypes.c_float*nn)()
libCentroid.imSplint(arr, ny, nx, para, nxt, nyt, yat)
imy = np.array(yat[:]).reshape([nxt, nyt])
#imy centroid
libCentroid.centroidWgt(yat, nyt, nxt, para)
imy_cx = para[3]
imy_cy = para[4]
#plot check
fig = plt.figure(figsize=(12, 6))
##imx_plot
ax = plt.subplot(1,2,1)
cpix= int(nx/2)
dpix= 10
plt.imshow(imx[cpix-dpix:cpix+dpix, cpix-dpix:cpix+dpix], origin='lower')
plt.plot(imx_cy-cpix+dpix, imx_cx-cpix+dpix, 'bx', ms = 20)
from scipy import ndimage
cx, cy = ndimage.center_of_mass(imx)
plt.plot(cy-cpix+dpix, cx-cpix+dpix, 'r+', ms = 20)
maxIndx = np.argmax(imx)
maxIndx = np.unravel_index(maxIndx, np.array(imx).shape)
imgMaxPix_x = maxIndx[1]
imgMaxPix_y = maxIndx[0]
apSizeInPix = 23
imgT = np.zeros_like(imx)
imgT[imgMaxPix_y-apSizeInPix:imgMaxPix_y+apSizeInPix+1,
imgMaxPix_x-apSizeInPix:imgMaxPix_x+apSizeInPix+1] = \
imx[imgMaxPix_y-apSizeInPix:imgMaxPix_y+apSizeInPix+1,
imgMaxPix_x-apSizeInPix:imgMaxPix_x+apSizeInPix+1]
cx, cy = ndimage.center_of_mass(imgT)
plt.plot(cy - cpix+dpix, cx - cpix+dpix, 'b+', ms = 15)
##imy_plot
ax = plt.subplot(1,2,2)
cpix= int(nxt/2)
dpix= 10
plt.imshow(imy[cpix-dpix:cpix+dpix, cpix-dpix:cpix+dpix], origin='lower')
plt.plot(imy_cy-cpix+dpix, imy_cx-cpix+dpix, 'bx', ms = 20)
plt.plot([dpix, dpix],[0,dpix],'w:')
plt.plot([0, dpix],[dpix,dpix],'w:')
cx, cy = ndimage.center_of_mass(imy)
plt.plot(cy-cpix+dpix, cx-cpix+dpix, 'r+', ms = 20)
maxIndx = np.argmax(imy)
maxIndx = np.unravel_index(maxIndx, np.array(imy).shape)
imgMaxPix_x = maxIndx[1]
imgMaxPix_y = maxIndx[0]
apSizeInPix = 23
imgT = np.zeros_like(imy)
imgT[imgMaxPix_y-apSizeInPix:imgMaxPix_y+apSizeInPix+1,
imgMaxPix_x-apSizeInPix:imgMaxPix_x+apSizeInPix+1] = \
imy[imgMaxPix_y-apSizeInPix:imgMaxPix_y+apSizeInPix+1,
imgMaxPix_x-apSizeInPix:imgMaxPix_x+apSizeInPix+1]
cx, cy = ndimage.center_of_mass(imgT)
plt.plot(cy - cpix+dpix, cx - cpix+dpix, 'b+', ms = 15)
plt.show()
import numpy as np
import matplotlib.pyplot as plt
import ctypes
import galsim
libCentroid = ctypes.CDLL('../libCentroid.so') # CDLL加载库
print('load libCenroid')
libCentroid.centroidWgt.argtypes = [ctypes.POINTER(ctypes.c_float), ctypes.c_int, ctypes.c_int, ctypes.POINTER(ctypes.c_double)]
'''
IMAGE_WIDTH = 180
IMAGE_HEIGHT = 180
dx = 0.#275
dy =-0.#393
center_x = 90 +dx #89.5 #IMAGE_WIDTH/2 -0.
center_y = 90 +dy #89.5 #IMAGE_HEIGHT/2+0.
R = np.sqrt(center_x**2 + center_y**2)
Gauss_map = np.zeros((IMAGE_HEIGHT, IMAGE_WIDTH))
for i in range(IMAGE_HEIGHT):
for j in range(IMAGE_WIDTH):
dis = (i-center_y)**2+(j-center_x)**2
Gauss_map[i, j] = np.exp(-0.5*dis/R)
ymap = galsim.InterpolatedImage(galsim.ImageF(Gauss_map), scale=0.01)
zmap = ymap.shear(g1 = 0.15, g2=0.27)
Gauss_map = zmap.drawImage(nx = 180, ny = 180, scale=0.01).array
fig=plt.figure(figsize=(5,5))
plt.imshow(Gauss_map, origin='lower')
imx = Gauss_map#/np.sum(Gauss_map)
ny,nx = imx.shape
print(nx, ny, np.max(imx))
nn = nx*ny
arr = (ctypes.c_float*nn)()
arr[:] = imx.reshape(nn)
para = (ctypes.c_double*10)()
libCentroid.centroidWgt(arr, ny, nx, para)
print('haha')
print(para[0:5])
cx = para[3]
cy = para[4]
print('{:}'.format(cx), '{:}'.format(cy))
from scipy import ndimage
cx, cy = ndimage.center_of_mass(imx)
print('center_of_mass:', cx, cy)
plt.plot(para[4], para[3], 'bx', ms = 15)
plt.plot(cy, cx, 'r+', ms = 15)
plt.annotate('dx,dy:', [10, 170], color='w')
plt.annotate('{:8.5}, {:8.5}'.format(dx, dy), [10, 160], color='w')
plt.annotate('cx, cy:', [10, 150], color='w')
plt.annotate('{:0<8.5}, {:0<8.5}'.format(center_x, center_y), [10, 140], color='w')
plt.annotate('{:0<8.5}, {:0<8.5}'.format(para[4], para[3]), [10, 130], color='w')
plt.annotate('{:0<8.5}, {:0<8.5}'.format(cy, cx), [10, 120], color='w')
'''
from scipy.io import loadmat
data = loadmat('/Users/chengliangwei/csstPSFdata/CSSOS_psf_ciomp/ccd13/wave_1/5_psf_array/psf_10.mat')
#plt.imshow(data['psf'])
#plt.show()
imx = data['psf']
imx = imx/np.sum(imx)
ny,nx = imx.shape
print(nx, ny)
nn = nx*ny
arr = (ctypes.c_float*nn)()
arr[:] = imx.reshape(nn)
para = (ctypes.c_double*10)()
print(arr[0:10])
#libCentroid.centroidWgt(arr, ny, nx, para)
nxt = nyt = 160
nn = nxt*nyt
yat = (ctypes.c_float*nn)()
libCentroid.centroidWgt(arr, ny, nx, para,nxt, nyt, yat)
mm = np.array(yat[:]).reshape([nxt, nyt])
imx = mm
print('haha')
print(para[0:5])
cx = para[3]
cy = para[4]
print(cx, cy)
fig = plt.figure(figsize=(12, 12))
cpix = 80
dpix = 10
#plt.imshow(np.log10(imx[cpix-dpix:cpix+dpix, cpix-dpix:cpix+dpix]), origin='lower')
plt.imshow(imx[cpix-dpix:cpix+dpix, cpix-dpix:cpix+dpix], origin='lower')
plt.plot(cy - cpix+dpix, cx - cpix+dpix, 'bx', ms = 20)
'''
from scipy import ndimage
cx, cy = ndimage.center_of_mass(imx)
print(cx, cy)
plt.plot(cy - cpix+dpix, cx - cpix+dpix, 'r+', ms = 20)
maxIndx = np.argmax(imx)
maxIndx = np.unravel_index(maxIndx, np.array(imx).shape)
imgMaxPix_x = maxIndx[1]
imgMaxPix_y = maxIndx[0]
apSizeInPix = 23
imgT = np.zeros_like(imx)
imgT[imgMaxPix_y-apSizeInPix:imgMaxPix_y+apSizeInPix+1,
imgMaxPix_x-apSizeInPix:imgMaxPix_x+apSizeInPix+1] = \
imx[imgMaxPix_y-apSizeInPix:imgMaxPix_y+apSizeInPix+1,
imgMaxPix_x-apSizeInPix:imgMaxPix_x+apSizeInPix+1]
cx, cy = ndimage.center_of_mass(imgT)
print(cx, cy)
plt.plot(cy - cpix+dpix, cx - cpix+dpix, 'b+', ms = 15)
maxIndx = np.argmax(imx)
maxIndx = np.unravel_index(maxIndx, np.array(imx).shape)
imgMaxPix_x = maxIndx[1]
imgMaxPix_y = maxIndx[0]
apSizeInPix =5
imgT = np.zeros_like(imx)
imgT[imgMaxPix_y-apSizeInPix:imgMaxPix_y+apSizeInPix+1,
imgMaxPix_x-apSizeInPix:imgMaxPix_x+apSizeInPix+1] = \
imx[imgMaxPix_y-apSizeInPix:imgMaxPix_y+apSizeInPix+1,
imgMaxPix_x-apSizeInPix:imgMaxPix_x+apSizeInPix+1]
cx, cy = ndimage.center_of_mass(imgT)
print(cx, cy)
plt.plot(cy - cpix+dpix, cx - cpix+dpix, 'm+', ms = 10)
print('maxPix:', imgMaxPix_x, imgMaxPix_y, imx[imgMaxPix_y, imgMaxPix_x])
'''
plt.show()
#OPTS += -D
CC = gcc
OPTIMIZE = -fPIC -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 = centroidWgt.o nrutil.o
EXEC = centroid.X
all: $(EXEC)
$(EXEC): $(OBJS)
$(CLINK) $(CFLAGS) -o $@ $(OBJS) $(CLIB)
$(OBJS): nrutil.h Makefile
.PHONY : clean
clean:
rm -f *.o $(EXEC)
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "nrutil.h"
//y is the input image, nx and ny are the pixels along x and y axis;
//yt and residu are the required matrix for the function
//para contains the needed parameters: centx=para[3], centy=para[4]
void star_gaus(float **y,int nx,int ny,double *para);
void main()
{
int i, j;
int IMAGE_WIDTH, IMAGE_HEIGHT;
float center_x, center_y;
float R, dis;
float **Gauss_map;
double *para;
//double para[10];
IMAGE_WIDTH = 180;
IMAGE_HEIGHT = 180;
center_x = IMAGE_WIDTH/2 -0.;
center_y = IMAGE_HEIGHT/2+0.;
R = sqrt(center_x*center_x + center_y*center_y);
//Gauss_map = np.zeros((IMAGE_HEIGHT, IMAGE_WIDTH))
Gauss_map = matrix(0, IMAGE_HEIGHT-1, 0, IMAGE_WIDTH-1);
para = dvector(0,10);
for(i=0; i<IMAGE_HEIGHT; i++)
{
for(j=0; j<IMAGE_WIDTH; j++)
{
//#dis = np.sqrt((i-center_y)**2+(j-center_x)**2)
dis = (i-center_y)*(i-center_y) + (j-center_x)*(j-center_x);
//dis = sqrt(dis);
Gauss_map[i][j] = exp(-0.5*dis/R);
}
}
star_gaus(Gauss_map, IMAGE_HEIGHT, IMAGE_WIDTH, para);
printf("ok?\n");
printf("para:%f %f\n", para[3], para[4]);
//free_matrix(Gauss_map,0,IMAGE_HEIGHT-1,0,IMAGE_WIDTH-1);
//free_dvector(para,0,4);
}
void centroidWgt(float *y, int nx, int ny, double *para)
{
int i, j, k;
float **yy;
double **basef;
double **coff;
yy = matrix(0, nx-1, 0, ny-1);
for(i=0; i<nx; i++)
for(j=0; j<ny; j++)
{
k = j + i*nx;
yy[i][j] = y[k];
}
star_gaus(yy, nx, ny, para);
}
//void star_gaus(float **y,int nx,int ny,float **yt,float **residu,double *para)
void star_gaus(float **y,int nx,int ny,double *para)
{
void gasfit_2D(double **x, double *y,int np,double *para,double bg0, int fbg);
double **xs,*ys,ymax;
int i,j,np,npt,im,jm,**xi,k;
double kc,bgc,sigmac,xc,yc,sigma2v,det;
double bg0=0;
int fbg=0;
np=nx*ny;
xs=dmatrix(0,np-1,0,1);
ys=dvector(0,np-1);
xi=imatrix(0,np-1,0,1);
ymax=y[0][0];
for(i=0;i<nx;i++)for(j=0;j<ny;j++){
if(ymax<y[i][j]){ymax=y[i][j];im=i;jm=j;}
}printf("x=%d,y=%d,ymax=%e\n",im,jm,ymax);
int cutPix;
cutPix = 23;
npt=0;
for(i=-cutPix;i<=cutPix;i++){
for(j=-cutPix;j<=cutPix;j++){
xs[npt][0]=xi[npt][0]=i+im;
xs[npt][1]=xi[npt][1]=j+jm;
ys[npt]=y[im+i][jm+j];
npt++;
}
}
printf("check:: %f %f %f %d\n",ys[0], ys[10], ys[100], npt );
gasfit_2D(xs, ys,npt,para,bg0,fbg);
kc=para[0];sigmac=para[1];bgc=para[2];xc=para[3];yc=para[4];
printf("xc, yc: %f %f\n",yc, xc);
//printf("%e %e %e %e %e\n",kc,sigmac,bgc,xc,yc);
/*
sigma2v=-1./(2.*sigmac*sigmac);
for(i=0;i<nx;i++){
for(j=0;j<ny;j++){
det=DSQR(i-xc)+DSQR(j-yc);
yt[i][j]=kc*exp(det*sigma2v)+bgc;
residu[i][j]=y[i][j]-yt[i][j];
}
}
*/
free_dmatrix(xs,0,np-1,0,1);
free_imatrix(xi,0,np-1,0,1);
free_dvector(ys,0,np-1);
}
void gasfit_2D(double **x, double *y,int np,double *para,double bg0,int fbg)
{
void search_2D(double **x,double *y,int np,double kc,double kd,double sigmac,
double sigmad,double bgc,double bgd,double xc, double xd,
double yc, double yd,double *para,double bg0,int fbg);
int i,j,k,imax=0,isigma;
double ymax,ymin,kc,kd,sigmac,sigmad,bgc,bgd,ysigma,xc,yc,xd,yd;
double det,dett;
ymin=ymax=y[imax];
for(i=1;i<np;i++){
if(ymax<y[i]){imax=i;ymax=y[i];}
if(ymin>y[i])ymin=y[i];
}
//printf("ymax=%e,imax=%d\n",ymax,imax);
kc=ymax;kd=ymax/12.;
det=ysigma=kc*exp(-0.5);//printf("det=%e,kc=%e\n",det,kc);
for(i=0;i<np;i++){
dett=fabs(ysigma-y[i]);//printf("dett=%e,det=%e,i=%d\n",dett,det,i);
//if(dett<det && y[i]!=kc){det=dett;isigma=i;}
if(dett<det){det=dett;isigma=i;}
}
xc=x[imax][0];yc=x[imax][1];printf("isigma=%d,imax=%d\n",isigma,imax);
sigmac=sqrt(DSQR(xc-x[isigma][0])+DSQR(yc-x[isigma][1]))*1.;
xd=yd=sigmac*0.25;
sigmad=0.25*sigmac;
bgc=0.;bgd=fabs(ymin);
//printf("sigma=%e\n",sigmac);
for(i=0;i<4;i++){
printf("i = %d\n", i);
//printf("%e %e %e %e %e k2=%e\n",kc,sigmac,bgc,xc,yc,para[5]);
search_2D(x,y,np,kc,kd,sigmac,sigmad,bgc,bgd,xc,xd,yc,yd,para,bg0,fbg);
kd*=0.33;sigmad*=0.33;bgd*=0.33;xd*=0.33;yd*=0.33;
kc=para[0];sigmac=para[1];bgc=para[2];xc=para[3];yc=para[4];
}
if(fbg==0)para[2]=bg0;
}
void search_2D(double **x,double *y,int np,double kc,double kd,double sigmac,
double sigmad,double bgc,double bgd,double xc, double xd,
double yc, double yd,double *para,double bg0,int fbg)
{
double k,sigma,bg,k2,k20,sigma2v,det,xt,yt;
int i,j,l,m,p,q;
sigma2v=-1./(2.*sigmac*sigmac);
bg=bgc;
k20=0;
for(m=0;m<np;m++){
det=DSQR(x[m][0]-xc)+DSQR(x[m][1]-yc);
det=kc*exp(det*sigma2v)+bgc-y[m];
k20+=det*det;
}
for(i=-4;i<=4;i++){
k=kc+i*kd;
if(k>0){
for(j=-4;j<=4;j++){
sigma=sigmac+j*sigmad;
if(sigma>0){
sigma2v=-1./(2.*sigma*sigma);
for(p=-4;p<=4;p++){
xt=xc+p*xd;
for(q=-4;q<=4;q++){
yt=yc+q*yd;
k2=0;
if(fbg==0){
bg=bg0;
for(m=0;m<np;m++){
det=DSQR(x[m][0]-xt)+DSQR(x[m][1]-yt);
det=k*exp(det*sigma2v)+bg-y[m];
k2+=det*det;
}
}
else{
for(l=-4;l<=4;l++){
bg=bgc+l*bgd;
for(m=0;m<np;m++){
det=DSQR(x[m][0]-xt)+DSQR(x[m][1]-yt);
det=k*exp(det*sigma2v)+bg-y[m];
k2+=det*det;
}
}
}
//printf("k20=%e k2=%e\n",k20,k2);
if(k2<=k20){k20=k2;para[5]=k2;
para[0]=k;para[1]=sigma;para[2]=bg;para[3]=xt;para[4]=yt;}
}
}
}
}
}
}
}
#if defined(__STDC__) || defined(ANSI) || defined(NRANSI) /* ANSI */
#include <stdio.h>
#include <stddef.h>
#include <stdlib.h>
#define NR_END 1
#define FREE_ARG char*
void nrerror(char error_text[])
/* Numerical Recipes standard error handler */
{
fprintf(stderr,"Numerical Recipes run-time error...\n");
fprintf(stderr,"%s\n",error_text);
fprintf(stderr,"...now exiting to system...\n");
exit(1);
}
float *vector(long nl, long nh)
/* allocate a float vector with subscript range v[nl..nh] */
{
float *v;
v=(float *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(float)));
if (!v) nrerror("allocation failure in vector()");
return v-nl+NR_END;
}
int *ivector(long nl, long nh)
/* allocate an int vector with subscript range v[nl..nh] */
{
int *v;
v=(int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int)));
if (!v) nrerror("allocation failure in ivector()");
return v-nl+NR_END;
}
unsigned char *cvector(long nl, long nh)
/* allocate an unsigned char vector with subscript range v[nl..nh] */
{
unsigned char *v;
v=(unsigned char *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(unsigned char)));
if (!v) nrerror("allocation failure in cvector()");
return v-nl+NR_END;
}
unsigned long *lvector(long nl, long nh)
/* allocate an unsigned long vector with subscript range v[nl..nh] */
{
unsigned long *v;
v=(unsigned long *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(long)));
if (!v) nrerror("allocation failure in lvector()");
return v-nl+NR_END;
}
double *dvector(long nl, long nh)
/* allocate a double vector with subscript range v[nl..nh] */
{
double *v;
v=(double *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(double)));
if (!v) nrerror("allocation failure in dvector()");
return v-nl+NR_END;
}
float **matrix(long nrl, long nrh, long ncl, long nch)
/* allocate a float matrix with subscript range m[nrl..nrh][ncl..nch] */
{
long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
float **m;
/* allocate pointers to rows */
m=(float **) malloc((size_t)((nrow+NR_END)*sizeof(float*)));
if (!m) nrerror("allocation failure 1 in matrix()");
m += NR_END;
m -= nrl;
/* allocate rows and set pointers to them */
m[nrl]=(float *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float)));
if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
m[nrl] += NR_END;
m[nrl] -= ncl;
for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
/* return pointer to array of pointers to rows */
return m;
}
double **dmatrix(long nrl, long nrh, long ncl, long nch)
/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
{
long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
double **m;
/* allocate pointers to rows */
m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*)));
if (!m) nrerror("allocation failure 1 in matrix()");
m += NR_END;
m -= nrl;
/* allocate rows and set pointers to them */
m[nrl]=(double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double)));
if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
m[nrl] += NR_END;
m[nrl] -= ncl;
for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
/* return pointer to array of pointers to rows */
return m;
}
int **imatrix(long nrl, long nrh, long ncl, long nch)
/* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */
{
long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
int **m;
/* allocate pointers to rows */
m=(int **) malloc((size_t)((nrow+NR_END)*sizeof(int*)));
if (!m) nrerror("allocation failure 1 in matrix()");
m += NR_END;
m -= nrl;
/* allocate rows and set pointers to them */
m[nrl]=(int *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(int)));
if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
m[nrl] += NR_END;
m[nrl] -= ncl;
for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
/* return pointer to array of pointers to rows */
return m;
}
float **submatrix(float **a, long oldrl, long oldrh, long oldcl, long oldch,
long newrl, long newcl)
/* point a submatrix [newrl..][newcl..] to a[oldrl..oldrh][oldcl..oldch] */
{
long i,j,nrow=oldrh-oldrl+1,ncol=oldcl-newcl;
float **m;
/* allocate array of pointers to rows */
m=(float **) malloc((size_t) ((nrow+NR_END)*sizeof(float*)));
if (!m) nrerror("allocation failure in submatrix()");
m += NR_END;
m -= newrl;
/* set pointers to rows */
for(i=oldrl,j=newrl;i<=oldrh;i++,j++) m[j]=a[i]+ncol;
/* return pointer to array of pointers to rows */
return m;
}
float **convert_matrix(float *a, long nrl, long nrh, long ncl, long nch)
/* allocate a float matrix m[nrl..nrh][ncl..nch] that points to the matrix
declared in the standard C manner as a[nrow][ncol], where nrow=nrh-nrl+1
and ncol=nch-ncl+1. The routine should be called with the address
&a[0][0] as the first argument. */
{
long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1;
float **m;
/* allocate pointers to rows */
m=(float **) malloc((size_t) ((nrow+NR_END)*sizeof(float*)));
if (!m) nrerror("allocation failure in convert_matrix()");
m += NR_END;
m -= nrl;
/* set pointers to rows */
m[nrl]=a-ncl;
for(i=1,j=nrl+1;i<nrow;i++,j++) m[j]=m[j-1]+ncol;
/* return pointer to array of pointers to rows */
return m;
}
float ***f3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
/* allocate a float 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] */
{
long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
float ***t;
/* allocate pointers to pointers to rows */
t=(float ***) malloc((size_t)((nrow+NR_END)*sizeof(float**)));
if (!t) nrerror("allocation failure 1 in f3tensor()");
t += NR_END;
t -= nrl;
/* allocate pointers to rows and set pointers to them */
t[nrl]=(float **) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float*)));
if (!t[nrl]) nrerror("allocation failure 2 in f3tensor()");
t[nrl] += NR_END;
t[nrl] -= ncl;
/* allocate rows and set pointers to them */
t[nrl][ncl]=(float *) malloc((size_t)((nrow*ncol*ndep+NR_END)*sizeof(float)));
if (!t[nrl][ncl]) nrerror("allocation failure 3 in f3tensor()");
t[nrl][ncl] += NR_END;
t[nrl][ncl] -= ndl;
for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep;
for(i=nrl+1;i<=nrh;i++) {
t[i]=t[i-1]+ncol;
t[i][ncl]=t[i-1][ncl]+ncol*ndep;
for(j=ncl+1;j<=nch;j++) t[i][j]=t[i][j-1]+ndep;
}
/* return pointer to array of pointers to rows */
return t;
}
double ***d3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
/* allocate a double 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] */
{
long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
double ***t;
/* allocate pointers to pointers to rows */
t=(double ***) malloc((size_t)((nrow+NR_END)*sizeof(double**)));
if (!t) nrerror("allocation failure 1 in f3tensor()");
t += NR_END;
t -= nrl;
/* allocate pointers to rows and set pointers to them */
t[nrl]=(double **) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double*)));
if (!t[nrl]) nrerror("allocation failure 2 in f3tensor()");
t[nrl] += NR_END;
t[nrl] -= ncl;
/* allocate rows and set pointers to them */
t[nrl][ncl]=(double *) malloc((size_t)((nrow*ncol*ndep+NR_END)*sizeof(double)));
if (!t[nrl][ncl]) nrerror("allocation failure 3 in f3tensor()");
t[nrl][ncl] += NR_END;
t[nrl][ncl] -= ndl;
for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep;
for(i=nrl+1;i<=nrh;i++) {
t[i]=t[i-1]+ncol;
t[i][ncl]=t[i-1][ncl]+ncol*ndep;
for(j=ncl+1;j<=nch;j++) t[i][j]=t[i][j-1]+ndep;
}
/* return pointer to array of pointers to rows */
return t;
}
unsigned char **cmatrix(long nrl, long nrh, long ncl, long nch)
/* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */
{
long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
unsigned char **m;
/* allocate pointers to rows */
m=(unsigned char **) malloc((size_t)((nrow+NR_END)*sizeof(unsigned char*)));
if (!m) nrerror("allocation failure 1 in matrix()");
m += NR_END;
m -= nrl;
/* allocate rows and set pointers to them */
m[nrl]=(unsigned char *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(unsigned char)));
if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
m[nrl] += NR_END;
m[nrl] -= ncl;
for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
/* return pointer to array of pointers to rows */
return m;
}
void free_vector(float *v, long nl, long nh)
/* free a float vector allocated with vector() */
{
free((FREE_ARG) (v+nl-NR_END));
}
void free_ivector(int *v, long nl, long nh)
/* free an int vector allocated with ivector() */
{
free((FREE_ARG) (v+nl-NR_END));
}
void free_cvector(unsigned char *v, long nl, long nh)
/* free an unsigned char vector allocated with cvector() */
{
free((FREE_ARG) (v+nl-NR_END));
}
void free_lvector(unsigned long *v, long nl, long nh)
/* free an unsigned long vector allocated with lvector() */
{
free((FREE_ARG) (v+nl-NR_END));
}
void free_dvector(double *v, long nl, long nh)
/* free a double vector allocated with dvector() */
{
free((FREE_ARG) (v+nl-NR_END));
}
void free_matrix(float **m, long nrl, long nrh, long ncl, long nch)
/* free a float matrix allocated by matrix() */
{
free((FREE_ARG) (m[nrl]+ncl-NR_END));
free((FREE_ARG) (m+nrl-NR_END));
}
void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch)
/* free a double matrix allocated by dmatrix() */
{
free((FREE_ARG) (m[nrl]+ncl-NR_END));
free((FREE_ARG) (m+nrl-NR_END));
}
void free_imatrix(int **m, long nrl, long nrh, long ncl, long nch)
/* free an int matrix allocated by imatrix() */
{
free((FREE_ARG) (m[nrl]+ncl-NR_END));
free((FREE_ARG) (m+nrl-NR_END));
}
void free_submatrix(float **b, long nrl, long nrh, long ncl, long nch)
/* free a submatrix allocated by submatrix() */
{
free((FREE_ARG) (b+nrl-NR_END));
}
void free_convert_matrix(float **b, long nrl, long nrh, long ncl, long nch)
/* free a matrix allocated by convert_matrix() */
{
free((FREE_ARG) (b+nrl-NR_END));
}
void free_f3tensor(float ***t, long nrl, long nrh, long ncl, long nch,
long ndl, long ndh)
/* free a float f3tensor allocated by f3tensor() */
{
free((FREE_ARG) (t[nrl][ncl]+ndl-NR_END));
free((FREE_ARG) (t[nrl]+ncl-NR_END));
free((FREE_ARG) (t+nrl-NR_END));
}
void free_d3tensor(double ***t, long nrl, long nrh, long ncl, long nch,
long ndl, long ndh)
/* free a double f3tensor allocated by f3tensor() */
{
free((FREE_ARG) (t[nrl][ncl]+ndl-NR_END));
free((FREE_ARG) (t[nrl]+ncl-NR_END));
free((FREE_ARG) (t+nrl-NR_END));
}
void free_cmatrix(unsigned char **m, long nrl, long nrh, long ncl, long nch)
/* free a character matrix allocated by matrix() */
{
free((FREE_ARG) (m[nrl]+ncl-NR_END));
free((FREE_ARG) (m+nrl-NR_END));
}
#else /* ANSI */
/* traditional - K&R */
#include <stdio.h>
#define NR_END 1
#define FREE_ARG char*
void nrerror(error_text)
char error_text[];
/* Numerical Recipes standard error handler */
{
void exit();
fprintf(stderr,"Numerical Recipes run-time error...\n");
fprintf(stderr,"%s\n",error_text);
fprintf(stderr,"...now exiting to system...\n");
exit(1);
}
float *vector(nl,nh)
long nh,nl;
/* allocate a float vector with subscript range v[nl..nh] */
{
float *v;
v=(float *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(float)));
if (!v) nrerror("allocation failure in vector()");
return v-nl+NR_END;
}
int *ivector(nl,nh)
long nh,nl;
/* allocate an int vector with subscript range v[nl..nh] */
{
int *v;
v=(int *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(int)));
if (!v) nrerror("allocation failure in ivector()");
return v-nl+NR_END;
}
unsigned char *cvector(nl,nh)
long nh,nl;
/* allocate an unsigned char vector with subscript range v[nl..nh] */
{
unsigned char *v;
v=(unsigned char *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(unsigned char)));
if (!v) nrerror("allocation failure in cvector()");
return v-nl+NR_END;
}
unsigned long *lvector(nl,nh)
long nh,nl;
/* allocate an unsigned long vector with subscript range v[nl..nh] */
{
unsigned long *v;
v=(unsigned long *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(long)));
if (!v) nrerror("allocation failure in lvector()");
return v-nl+NR_END;
}
double *dvector(nl,nh)
long nh,nl;
/* allocate a double vector with subscript range v[nl..nh] */
{
double *v;
v=(double *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(double)));
if (!v) nrerror("allocation failure in dvector()");
return v-nl+NR_END;
}
float **matrix(nrl,nrh,ncl,nch)
long nch,ncl,nrh,nrl;
/* allocate a float matrix with subscript range m[nrl..nrh][ncl..nch] */
{
long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
float **m;
/* allocate pointers to rows */
m=(float **) malloc((unsigned int)((nrow+NR_END)*sizeof(float*)));
if (!m) nrerror("allocation failure 1 in matrix()");
m += NR_END;
m -= nrl;
/* allocate rows and set pointers to them */
m[nrl]=(float *) malloc((unsigned int)((nrow*ncol+NR_END)*sizeof(float)));
if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
m[nrl] += NR_END;
m[nrl] -= ncl;
for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
/* return pointer to array of pointers to rows */
return m;
}
double **dmatrix(nrl,nrh,ncl,nch)
long nch,ncl,nrh,nrl;
/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
{
long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
double **m;
/* allocate pointers to rows */
m=(double **) malloc((unsigned int)((nrow+NR_END)*sizeof(double*)));
if (!m) nrerror("allocation failure 1 in matrix()");
m += NR_END;
m -= nrl;
/* allocate rows and set pointers to them */
m[nrl]=(double *) malloc((unsigned int)((nrow*ncol+NR_END)*sizeof(double)));
if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
m[nrl] += NR_END;
m[nrl] -= ncl;
for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
/* return pointer to array of pointers to rows */
return m;
}
int **imatrix(nrl,nrh,ncl,nch)
long nch,ncl,nrh,nrl;
/* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */
{
long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
int **m;
/* allocate pointers to rows */
m=(int **) malloc((unsigned int)((nrow+NR_END)*sizeof(int*)));
if (!m) nrerror("allocation failure 1 in matrix()");
m += NR_END;
m -= nrl;
/* allocate rows and set pointers to them */
m[nrl]=(int *) malloc((unsigned int)((nrow*ncol+NR_END)*sizeof(int)));
if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
m[nrl] += NR_END;
m[nrl] -= ncl;
for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
/* return pointer to array of pointers to rows */
return m;
}
float **submatrix(a,oldrl,oldrh,oldcl,oldch,newrl,newcl)
float **a;
long newcl,newrl,oldch,oldcl,oldrh,oldrl;
/* point a submatrix [newrl..][newcl..] to a[oldrl..oldrh][oldcl..oldch] */
{
long i,j,nrow=oldrh-oldrl+1,ncol=oldcl-newcl;
float **m;
/* allocate array of pointers to rows */
m=(float **) malloc((unsigned int) ((nrow+NR_END)*sizeof(float*)));
if (!m) nrerror("allocation failure in submatrix()");
m += NR_END;
m -= newrl;
/* set pointers to rows */
for(i=oldrl,j=newrl;i<=oldrh;i++,j++) m[j]=a[i]+ncol;
/* return pointer to array of pointers to rows */
return m;
}
float **convert_matrix(a,nrl,nrh,ncl,nch)
float *a;
long nch,ncl,nrh,nrl;
/* allocate a float matrix m[nrl..nrh][ncl..nch] that points to the matrix
declared in the standard C manner as a[nrow][ncol], where nrow=nrh-nrl+1
and ncol=nch-ncl+1. The routine should be called with the address
&a[0][0] as the first argument. */
{
long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1;
float **m;
/* allocate pointers to rows */
m=(float **) malloc((unsigned int) ((nrow+NR_END)*sizeof(float*)));
if (!m) nrerror("allocation failure in convert_matrix()");
m += NR_END;
m -= nrl;
/* set pointers to rows */
m[nrl]=a-ncl;
for(i=1,j=nrl+1;i<nrow;i++,j++) m[j]=m[j-1]+ncol;
/* return pointer to array of pointers to rows */
return m;
}
double ***d3tensor(nrl,nrh,ncl,nch,ndl,ndh)
long nch,ncl,ndh,ndl,nrh,nrl;
/* allocate a float 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] */
{
long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
double ***t;
/* allocate pointers to pointers to rows */
t=(double ***) malloc((unsigned int)((nrow+NR_END)*sizeof(double**)));
if (!t) nrerror("allocation failure 1 in f3tensor()");
t += NR_END;
t -= nrl;
/* allocate pointers to rows and set pointers to them */
t[nrl]=(double **) malloc((unsigned int)((nrow*ncol+NR_END)*sizeof(double*)));
if (!t[nrl]) nrerror("allocation failure 2 in f3tensor()");
t[nrl] += NR_END;
t[nrl] -= ncl;
/* allocate rows and set pointers to them */
t[nrl][ncl]=(double *) malloc((unsigned int)((nrow*ncol*ndep+NR_END)*sizeof(double)));
if (!t[nrl][ncl]) nrerror("allocation failure 3 in f3tensor()");
t[nrl][ncl] += NR_END;
t[nrl][ncl] -= ndl;
for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep;
for(i=nrl+1;i<=nrh;i++) {
t[i]=t[i-1]+ncol;
t[i][ncl]=t[i-1][ncl]+ncol*ndep;
for(j=ncl+1;j<=nch;j++) t[i][j]=t[i][j-1]+ndep;
}
/* return pointer to array of pointers to rows */
return t;
}
float ***f3tensor(nrl,nrh,ncl,nch,ndl,ndh)
long nch,ncl,ndh,ndl,nrh,nrl;
/* allocate a float 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] */
{
long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
float ***t;
/* allocate pointers to pointers to rows */
t=(float ***) malloc((unsigned int)((nrow+NR_END)*sizeof(float**)));
if (!t) nrerror("allocation failure 1 in f3tensor()");
t += NR_END;
t -= nrl;
/* allocate pointers to rows and set pointers to them */
t[nrl]=(float **) malloc((unsigned int)((nrow*ncol+NR_END)*sizeof(float*)));
if (!t[nrl]) nrerror("allocation failure 2 in f3tensor()");
t[nrl] += NR_END;
t[nrl] -= ncl;
/* allocate rows and set pointers to them */
t[nrl][ncl]=(float *) malloc((unsigned int)((nrow*ncol*ndep+NR_END)*sizeof(float)));
if (!t[nrl][ncl]) nrerror("allocation failure 3 in f3tensor()");
t[nrl][ncl] += NR_END;
t[nrl][ncl] -= ndl;
for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep;
for(i=nrl+1;i<=nrh;i++) {
t[i]=t[i-1]+ncol;
t[i][ncl]=t[i-1][ncl]+ncol*ndep;
for(j=ncl+1;j<=nch;j++) t[i][j]=t[i][j-1]+ndep;
}
/* return pointer to array of pointers to rows */
return t;
}
unsigned char **cmatrix(nrl,nrh,ncl,nch)
long nch,ncl,nrh,nrl;
/* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */
{
long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
unsigned char **m;
/* allocate pointers to rows */
m=(unsigned char **) malloc((unsigned int)((nrow+NR_END)*sizeof(unsigned char*)));
if (!m) nrerror("allocation failure 1 in matrix()");
m += NR_END;
m -= nrl;
/* allocate rows and set pointers to them */
m[nrl]=(unsigned char *) malloc((unsigned int)((nrow*ncol+NR_END)*sizeof(unsigned char)));
if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
m[nrl] += NR_END;
m[nrl] -= ncl;
for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
/* return pointer to array of pointers to rows */
return m;
}
void free_vector(v,nl,nh)
float *v;
long nh,nl;
/* free a float vector allocated with vector() */
{
free((FREE_ARG) (v+nl-NR_END));
}
void free_ivector(v,nl,nh)
int *v;
long nh,nl;
/* free an int vector allocated with ivector() */
{
free((FREE_ARG) (v+nl-NR_END));
}
void free_cvector(v,nl,nh)
long nh,nl;
unsigned char *v;
/* free an unsigned char vector allocated with cvector() */
{
free((FREE_ARG) (v+nl-NR_END));
}
void free_lvector(v,nl,nh)
long nh,nl;
unsigned long *v;
/* free an unsigned long vector allocated with lvector() */
{
free((FREE_ARG) (v+nl-NR_END));
}
void free_dvector(v,nl,nh)
double *v;
long nh,nl;
/* free a double vector allocated with dvector() */
{
free((FREE_ARG) (v+nl-NR_END));
}
void free_matrix(m,nrl,nrh,ncl,nch)
float **m;
long nch,ncl,nrh,nrl;
/* free a float matrix allocated by matrix() */
{
free((FREE_ARG) (m[nrl]+ncl-NR_END));
free((FREE_ARG) (m+nrl-NR_END));
}
void free_dmatrix(m,nrl,nrh,ncl,nch)
double **m;
long nch,ncl,nrh,nrl;
/* free a double matrix allocated by dmatrix() */
{
free((FREE_ARG) (m[nrl]+ncl-NR_END));
free((FREE_ARG) (m+nrl-NR_END));
}
void free_imatrix(m,nrl,nrh,ncl,nch)
int **m;
long nch,ncl,nrh,nrl;
/* free an int matrix allocated by imatrix() */
{
free((FREE_ARG) (m[nrl]+ncl-NR_END));
free((FREE_ARG) (m+nrl-NR_END));
}
void free_submatrix(b,nrl,nrh,ncl,nch)
float **b;
long nch,ncl,nrh,nrl;
/* free a submatrix allocated by submatrix() */
{
free((FREE_ARG) (b+nrl-NR_END));
}
void free_convert_matrix(b,nrl,nrh,ncl,nch)
float **b;
long nch,ncl,nrh,nrl;
/* free a matrix allocated by convert_matrix() */
{
free((FREE_ARG) (b+nrl-NR_END));
}
void free_f3tensor(t,nrl,nrh,ncl,nch,ndl,ndh)
float ***t;
long nch,ncl,ndh,ndl,nrh,nrl;
/* free a float f3tensor allocated by f3tensor() */
{
free((FREE_ARG) (t[nrl][ncl]+ndl-NR_END));
free((FREE_ARG) (t[nrl]+ncl-NR_END));
free((FREE_ARG) (t+nrl-NR_END));
}
void free_d3tensor(t,nrl,nrh,ncl,nch,ndl,ndh)
double ***t;
long nch,ncl,ndh,ndl,nrh,nrl;
/* free a float f3tensor allocated by f3tensor() */
{
free((FREE_ARG) (t[nrl][ncl]+ndl-NR_END));
free((FREE_ARG) (t[nrl]+ncl-NR_END));
free((FREE_ARG) (t+nrl-NR_END));
}
void free_cmatrix(m,nrl,nrh,ncl,nch)
unsigned **m;
long nch,ncl,nrh,nrl;
/* free a double matrix allocated by dmatrix() */
{
free((FREE_ARG) (m[nrl]+ncl-NR_END));
free((FREE_ARG) (m+nrl-NR_END));
}
#endif /* ANSI */
/* CAUTION: This is the ANSI C (only) version of the Numerical Recipes
utility file nrutil.h. Do not confuse this file with the same-named
file nrutil.h that may be supplied in a 'misc' subdirectory.
*That* file is the one from the book, and contains both ANSI and
traditional K&R versions, along with #ifdef macros to select the
correct version. *This* file contains only ANSI C. */
#ifndef _NR_UTILS_H_
#define _NR_UTILS_H_
static float sqrarg;
#define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)
static double dsqrarg;
#define DSQR(a) ((dsqrarg=(a)) == 0.0 ? 0.0 : dsqrarg*dsqrarg)
static double dmaxarg1,dmaxarg2;
#define DMAX(a,b) (dmaxarg1=(a),dmaxarg2=(b),(dmaxarg1) > (dmaxarg2) ?\
(dmaxarg1) : (dmaxarg2))
static double dminarg1,dminarg2;
#define DMIN(a,b) (dminarg1=(a),dminarg2=(b),(dminarg1) < (dminarg2) ?\
(dminarg1) : (dminarg2))
static float maxarg1,maxarg2;
#define FMAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) > (maxarg2) ?\
(maxarg1) : (maxarg2))
static float minarg1,minarg2;
#define FMIN(a,b) (minarg1=(a),minarg2=(b),(minarg1) < (minarg2) ?\
(minarg1) : (minarg2))
static long lmaxarg1,lmaxarg2;
#define LMAX(a,b) (lmaxarg1=(a),lmaxarg2=(b),(lmaxarg1) > (lmaxarg2) ?\
(lmaxarg1) : (lmaxarg2))
static long lminarg1,lminarg2;
#define LMIN(a,b) (lminarg1=(a),lminarg2=(b),(lminarg1) < (lminarg2) ?\
(lminarg1) : (lminarg2))
static int imaxarg1,imaxarg2;
#define IMAX(a,b) (imaxarg1=(a),imaxarg2=(b),(imaxarg1) > (imaxarg2) ?\
(imaxarg1) : (imaxarg2))
static int iminarg1,iminarg2;
#define IMIN(a,b) (iminarg1=(a),iminarg2=(b),(iminarg1) < (iminarg2) ?\
(iminarg1) : (iminarg2))
#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
void nrerror(char error_text[]);
float *vector(long nl, long nh);
int *ivector(long nl, long nh);
unsigned char *cvector(long nl, long nh);
unsigned long *lvector(long nl, long nh);
double *dvector(long nl, long nh);
float **matrix(long nrl, long nrh, long ncl, long nch);
double **dmatrix(long nrl, long nrh, long ncl, long nch);
int **imatrix(long nrl, long nrh, long ncl, long nch);
float **submatrix(float **a, long oldrl, long oldrh, long oldcl, long oldch,
long newrl, long newcl);
float **convert_matrix(float *a, long nrl, long nrh, long ncl, long nch);
float ***f3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh);
void free_vector(float *v, long nl, long nh);
void free_ivector(int *v, long nl, long nh);
void free_cvector(unsigned char *v, long nl, long nh);
void free_lvector(unsigned long *v, long nl, long nh);
void free_dvector(double *v, long nl, long nh);
void free_matrix(float **m, long nrl, long nrh, long ncl, long nch);
void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch);
void free_imatrix(int **m, long nrl, long nrh, long ncl, long nch);
void free_submatrix(float **b, long nrl, long nrh, long ncl, long nch);
void free_convert_matrix(float **b, long nrl, long nrh, long ncl, long nch);
void free_f3tensor(float ***t, long nrl, long nrh, long ncl, long nch,
long ndl, long ndh);
int ***i3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh);
void free_i3tensor(int ***t, long nrl, long nrh, long ncl, long nch,
long ndl, long ndh);
unsigned char ***b3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh);
void free_b3tensor(unsigned char ***t, long nrl, long nrh, long ncl, long nch,
long ndl, long ndh);
double ***d3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh);
void free_d3tensor(double ***t, long nrl, long nrh, long ncl, long nch,
long ndl, long ndh);
unsigned char **cmatrix(long nrl, long nrh, long ncl, long nch);
void free_cmatrix(unsigned char **m, long nrl, long nrh, long ncl, long nch);
#endif /* _NR_UTILS_H_ */
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