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