Skip to content

Commit

Permalink
Define centroid w.r.t. lower-left corners of pixels.
Browse files Browse the repository at this point in the history
  • Loading branch information
@vrodgom
vrodgom committed Jan 8, 2019
1 parent c898602 commit 97a120f
Show file tree
Hide file tree
Showing 4 changed files with 89 additions and 87 deletions.
64 changes: 32 additions & 32 deletions notebooks/tutorial.ipynb
View file

Large diffs are not rendered by default.

16 changes: 6 additions & 10 deletions statmorph/statmorph.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -575,7 +575,6 @@ def _centroid(self):
assert M[0, 0] > 0 # already checked by constructor
yc = M[1, 0] / M[0, 0]
xc = M[0, 1] / M[0, 0]
yc += 0.5; xc += 0.5 # shift pixel positions

return np.array([xc, yc])

Expand All @@ -602,8 +601,8 @@ def _covariance_generic(self, xc, yc):
image = np.float64(self._cutout_stamp_maskzeroed_no_bg)

# Calculate moments w.r.t. given center
xc = xc - self.xmin_stamp - 0.5 # w.r.t. lower-left corner of pixels
yc = yc - self.ymin_stamp - 0.5
xc = xc - self.xmin_stamp
yc = yc - self.ymin_stamp
Mc = skimage.measure.moments_central(image, center=(yc, xc), order=2)
assert Mc[0, 0] > 0

Expand Down Expand Up @@ -1249,8 +1248,6 @@ def m20(self):
return -99.0 # invalid
yc = M[1, 0] / M[0, 0]
xc = M[0, 1] / M[0, 0]
# Note that we do not shift (yc, xc) by 0.5 pixels here, since
# (yc, xc) is only used as input for other skimage functions.

# Calculate second total central moment
Mc = skimage.measure.moments_central(image, center=(yc, xc), order=2)
Expand Down Expand Up @@ -2049,8 +2046,8 @@ def deviation(self):
self.flag = 1
return -99.0 # invalid

xp = sorted_xpeak[0] + 0.5 # center of pixel
yp = sorted_ypeak[0] + 0.5
xp = sorted_xpeak[0]
yp = sorted_ypeak[0]

# Calculate centroid
M = skimage.measure.moments(image, order=1)
Expand All @@ -2061,7 +2058,6 @@ def deviation(self):
return -99.0 # invalid
yc = M[1, 0] / M[0, 0]
xc = M[0, 1] / M[0, 0]
yc += 0.5; xc += 0.5 # shift pixel positions

area = np.sum(self._segmap_mid)
D = np.sqrt(np.pi/area) * np.sqrt((xp-xc)**2 + (yp-yc)**2)
Expand Down Expand Up @@ -2211,7 +2207,7 @@ def rmax_circ(self):
xc, yc = self._asymmetry_center

# Distances from all pixels to the center
ypos, xpos = np.mgrid[0:ny, 0:nx] + 0.5 # center of pixel
ypos, xpos = np.mgrid[0:ny, 0:nx]
distances = np.sqrt((ypos-yc)**2 + (xpos-xc)**2)

# Only consider pixels within the segmap.
Expand All @@ -2238,7 +2234,7 @@ def rmax_ellip(self):
xc, yc = self._asymmetry_center

theta = self.orientation_asymmetry
y, x = np.mgrid[0:ny, 0:nx] + 0.5 # center of pixel
y, x = np.mgrid[0:ny, 0:nx]

xprime = (x-xc)*np.cos(theta) + (y-yc)*np.sin(theta)
yprime = -(x-xc)*np.sin(theta) + (y-yc)*np.cos(theta)
Expand Down
58 changes: 29 additions & 29 deletions statmorph/tests/test_statmorph.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,40 +14,40 @@
__all__ = ['runall']

correct_values = {
'xc_centroid':82.16578971349672,
'yc_centroid':81.03461797343211,
'xc_centroid':81.66578971349672,
'yc_centroid':80.53461797343211,
'ellipticity_centroid':0.04855709579576,
'elongation_centroid':1.05103521775316,
'orientation_centroid':-0.85797602066194,
'xc_asymmetry':82.23816149429084,
'yc_asymmetry':80.76122537398356,
'ellipticity_asymmetry':0.04862719043001,
'elongation_asymmetry':1.05111265525025,
'orientation_asymmetry':-0.85900294920303,
'rpetro_circ':40.93685686885102,
'rpetro_ellip':41.65687159352140,
'rhalf_circ':21.60432043808139,
'rhalf_ellip':22.09299918129802,
'r20':11.69540955266178,
'r80':32.07888648441381,
'gini':0.39017652998146,
'm20':-1.54294057566943,
'sn_per_pixel':6.79978752136230,
'concentration':2.19101925047364,
'asymmetry':0.00216882009677,
'smoothness':0.00430980802075,
'orientation_centroid':-0.85797602066195,
'xc_asymmetry':82.23985214801982,
'yc_asymmetry':80.76076242700849,
'ellipticity_asymmetry':0.04806946962244,
'elongation_asymmetry':1.05049682522881,
'orientation_asymmetry':-0.85405676626920,
'rpetro_circ':40.93755531944313,
'rpetro_ellip':41.64283484446126,
'rhalf_circ':21.60803205322342,
'rhalf_ellip':22.08125638365687,
'r20':11.69548630967248,
'r80':32.07883340820674,
'gini':0.38993180299621,
'm20':-1.54448930789228,
'sn_per_pixel':6.80319166183472,
'concentration':2.19100140632153,
'asymmetry':0.00377345808887,
'smoothness':0.00430880839402,
'multimode':0.23423423423423,
'intensity':0.51203949030140,
'deviation':0.01522525597953,
'outer_asymmetry':-0.01939061100961,
'shape_asymmetry':0.16270478667804,
'sersic_amplitude':1296.95288046788846,
'sersic_rhalf':22.45788867007650,
'sersic_n':0.61206828330366,
'sersic_xc':81.56197595105462,
'sersic_yc':80.40465135535705,
'sersic_ellip':0.05083866177840,
'sersic_theta':2.47831543421010,
'outer_asymmetry':-0.01821399684443,
'shape_asymmetry':0.16308278287864,
'sersic_amplitude':1296.95288208155739,
'sersic_rhalf':22.45788866502031,
'sersic_n':0.61206828194077,
'sersic_xc':81.56197595338546,
'sersic_yc':80.40465135599014,
'sersic_ellip':0.05083866217150,
'sersic_theta':2.47831542907976,
'sky_mean':3.48760604858398,
'sky_median':-2.68543863296509,
'sky_sigma':150.91754150390625,
Expand Down
38 changes: 22 additions & 16 deletions statmorph/utils/image_diagnostics.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -91,7 +91,11 @@ def make_figure(morph):
M = np.max(image)
m_stretch, M_stretch = log_stretch([m, M])
xc, yc = morph._xc_stamp, morph._yc_stamp # centroid
xca, yca = morph._asymmetry_center # asym. centroid
xca, yca = morph._asymmetry_center # asym. center
xcs, ycs = morph._sersic_model.x_0.value, morph._sersic_model.y_0.value # Sersic center

# Plot everything w.r.t. centers of pixels (instead of lower-left corners):
xc += 0.5; yc += 0.5; xca += 0.5; yca += 0.5; xcs += 0.5; ycs += 0.5

##################
# Original image #
Expand Down Expand Up @@ -144,10 +148,8 @@ def make_figure(morph):
# Add background noise (for realism)
if morph.sky_sigma > 0:
sersic_model += np.random.normal(scale=morph.sky_sigma, size=(ny, nx))
ax.imshow(log_stretch(normalize(sersic_model, m=m, M=M)), cmap='gray', origin='lower',
vmin=m_stretch, vmax=M_stretch)
# Sersic center (within postage stamp)
xcs, ycs = morph._sersic_model.x_0.value, morph._sersic_model.y_0.value
ax.imshow(log_stretch(normalize(sersic_model, m=m, M=M)), cmap='gray',
origin='lower', vmin=m_stretch, vmax=M_stretch)
ax.plot(xcs, ycs, 'ro', markersize=5, label='Sérsic Center')
R = float(nx**2 + ny**2)
theta = morph.sersic_theta
Expand Down Expand Up @@ -196,11 +198,12 @@ def make_figure(morph):
######################
ax = get_ax(fig, 0, 3, nrows, ncols, wpanel, hpanel, htop, eps, wfig, hfig)
# Rotate image around asym. center
image_180 = skimage.transform.rotate(image, 180.0, center=(xca, yca))
# (note that skimage expects pixel positions at lower-left corners)
image_180 = skimage.transform.rotate(image, 180.0, center=(xca-0.5, yca-0.5))
image_res = image - image_180
# Apply symmetric mask
mask = morph._mask_stamp.copy()
mask_180 = skimage.transform.rotate(mask, 180.0, center=(xca, yca))
mask_180 = skimage.transform.rotate(mask, 180.0, center=(xca-0.5, yca-0.5))
mask_180 = mask_180 >= 0.5 # convert back to bool
mask_symmetric = mask | mask_180
image_res = np.where(~mask_symmetric, image_res, 0.0)
Expand Down Expand Up @@ -295,9 +298,11 @@ def make_figure(morph):
norm=matplotlib.colors.NoNorm())
sorted_flux_sums, sorted_xpeak, sorted_ypeak = morph._intensity_sums
if len(sorted_flux_sums) > 0:
ax.plot(sorted_xpeak[0] + 0.5, sorted_ypeak[0] + 0.5, 'bo', markersize=2, label='First Peak')
ax.plot(sorted_xpeak[0] + 0.5, sorted_ypeak[0] + 0.5, 'bo', markersize=2,
label='First Peak')
if len(sorted_flux_sums) > 1:
ax.plot(sorted_xpeak[1] + 0.5, sorted_ypeak[1] + 0.5, 'ro', markersize=2, label='Second Peak')
ax.plot(sorted_xpeak[1] + 0.5, sorted_ypeak[1] + 0.5, 'ro', markersize=2,
label='Second Peak')
# Some text
text = (r'$M = %.4f$' % (morph.multimode) + '\n' +
r'$I = %.4f$' % (morph.intensity) + '\n' +
Expand All @@ -320,13 +325,14 @@ def make_figure(morph):
ax.imshow(morph._segmap_shape_asym, cmap='gray', origin='lower')
ax.plot(xca, yca, 'bo', markersize=5, label='Asym. Center')
r = morph.rpetro_circ
ax.plot(xca + r*np.cos(theta_vec), yca + r*np.sin(theta_vec), 'b', label=r'$r_{\rm petro, circ}$')
ax.plot(xca + r*np.cos(theta_vec), yca + r*np.sin(theta_vec), 'b',
label=r'$r_{\rm petro, circ}$')
r = morph.rpetro_ellip
ax.plot(xca + r*np.cos(theta_vec), yca + r*np.sin(theta_vec), 'r', label=r'$r_{\rm petro, ellip}$')
ax.plot(xca + r*np.cos(theta_vec), yca + r*np.sin(theta_vec), 'r',
label=r'$r_{\rm petro, ellip}$')
r = morph.rmax_circ
ax.plot(np.floor(xca) + r*np.cos(theta_vec), np.floor(yca) + r*np.sin(theta_vec), 'c', lw=1.5, label=r'$r_{\rm max}$')
# ~ r = morph._petro_extent_flux * morph.rpetro_ellip
# ~ ax.plot(xca + r*np.cos(theta_vec), yca + r*np.sin(theta_vec), 'r--', label='%g*rpet_ellip' % (morph._petro_extent_flux))
ax.plot(np.floor(xca) + r*np.cos(theta_vec), np.floor(yca) + r*np.sin(theta_vec),
'c', lw=1.5, label=r'$r_{\rm max}$')
text = (r'$A_S = %.4f$' % (morph.shape_asymmetry))
ax.text(0.034, 0.034, text, fontsize=12,
horizontalalignment='left', verticalalignment='bottom',
Expand All @@ -339,8 +345,8 @@ def make_figure(morph):
ax.get_xaxis().set_visible(False)
ax.get_yaxis().set_visible(False)

# defaults: left = 0.125, right = 0.9, bottom = 0.1, top = 0.9, wspace = 0.2, hspace = 0.2
fig.subplots_adjust(left=eps/wfig, right=1-eps/wfig, bottom=eps/hfig, top=1.0-htop/hfig, wspace=eps/wfig, hspace=htop/hfig)
fig.subplots_adjust(left=eps/wfig, right=1-eps/wfig, bottom=eps/hfig,
top=1.0-htop/hfig, wspace=eps/wfig, hspace=htop/hfig)

#fig.savefig('test_segmap.png', dpi=150)

Expand Down

0 comments on commit 97a120f

Please sign in to comment.