tt.py

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