Commit e5adeeee authored by Ruxandra Cojocaru's avatar Ruxandra Cojocaru

Modifications to mask_builder in func.py. Added 2 unittests for mask_builder in tests/test_func.py

parent 571bf6ed
......@@ -1811,10 +1811,18 @@ def fillmask(img0, threshold, ymin=0, ymax=None, xmin=0, xmax=None,
return img
def mask_builder(image, mask_size):
"""Subroutine o build a circular mask for a lens
def mask_builder(image, mask_half_size):
"""Subroutine to build a circular mask for a lens
Diameter in pixels of the desired mask = 2 * mask_half_size + 1
"""
err_msg = ('mask_half_size has to be a positive number strictly smaller '
'than half the shortest dimension of image')
test_err(mask_half_size > 0 and mask_half_size < min(image.shape) / 2 ,
err_msg)
mask_diam = 2 * mask_half_size + 1
[m, n] = image.shape
# BUILD A MASK
......@@ -1824,44 +1832,40 @@ def mask_builder(image, mask_size):
mask0 = np.power(np.square(X) + np.square(Y), 0.5)
exponent = np.exp(-(X**2 + Y**2) / (2 * np.square(mask_size/2)))
exponent = np.exp(-(X**2 + Y**2) / (2 * np.square(mask_diam/2)))
exponent[mask0 < (mask_size/2)] = 1
exponent[mask0 < (mask_diam/2)] = 1
mask3 = (mask0 < (mask_size/2.*1.3)) * exponent
mask1 = (mask0 < (mask_diam/2.*1.3)) * exponent
# sort descending order
a = -np.sort(-image.reshape(-1))
maskThreshold = sum(mask3[mask3 == 1])/np.size(mask3)
maskThreshold = sum(mask1[mask1 == 1])/np.size(mask1)
b = a[int(np.round(len(a) * maskThreshold))]
print(b)
d1 = (image > b)
[x_offset,y_offset] = f_shift(d1, mask3)
[x_offset,y_offset] = f_shift(d1, mask1)
mask_centre = np.array((0, 0))
mask_centre[1] = x_offset + n/2
mask_centre[0] = y_offset + m/2
mask_centre[1] = x_offset + n//2
mask_centre[0] = y_offset + m//2
# Mask
[X1,Y1] = np.meshgrid(np.arange(n), np.arange(m))
[X1, Y1] = np.meshgrid(np.arange(n), np.arange(m))
X = X1 - mask_centre[1]
Y = Y1 - mask_centre[0]
mask0 = np.power(np.square(X) + np.square(Y), 0.5)
mask2 = (mask0 < (mask_size/2))
mask2 = (mask0 <= (mask_diam/2))
big_mask = (mask2 == True)
small_mask = big_mask[min(Y1[big_mask]) : max(Y1[big_mask]),
min(X1[big_mask]) : max(X1[big_mask])]
return [small_mask, big_mask, mask_centre]
return [big_mask, mask_centre]
def f_shift(a, t):
......@@ -1872,12 +1876,13 @@ def f_shift(a, t):
# Determine padding size in x and y dimension
size_t = t.shape
size_a = a.shape
out_size = np.array(size_t) + np.array(size_a) -1
out_size = np.array(size_t) + np.array(size_a) - 1
# Determine 2D cross correlation in Fourier domain
Ft = np.fft.fft2(t, s = (out_size[0], out_size[1]))
Fa = np.fft.fft2(a, s = (out_size[0], out_size[1]))
c = abs(np.fft.fftshift(np.fft.ifft2(Fa * np.conj(Ft))))
# Find peak, np.argmax return the linear index of 1st occurence of
......@@ -1887,6 +1892,7 @@ def f_shift(a, t):
[y_peak, x_peak] = np.unravel_index(imax, c.shape)
# Correct found peak location for image size
#~ corr_offset = np.array((y_peak + 1 - c.shape[0] // 2, x_peak + 1 - c.shape[1] // 2))
corr_offset = np.array((y_peak - c.shape[0] // 2, x_peak - c.shape[1] // 2))
# Write out offsets
......
......@@ -88,8 +88,103 @@ class FuncTestCase(unittest.TestCase):
np.testing.assert_array_equal(input_plane,
np.around(func.plane_fit(input_plane),
decimals = 4))
#~ def test_mask_builder(self):
#~ """Test returns expected mask center"""
#~
#~ image = np.random.rand(99, 99)
def test_mask_builder_centre(self):
"""Test it it returns expected mask center"""
def makeGaussian(size, fwhm = 3, center=None):
""" Make a square gaussian kernel.
size is the length of a side of the square
fwhm is full-width-half-maximum, which
can be thought of as an effective radius.
"""
x = np.arange(0, size, 1, float)
y = x[:,np.newaxis]
if center is None:
x0 = y0 = size // 2
else:
x0 = center[0]
y0 = center[1]
return np.exp(-4*np.log(2) * ((x-x0)**2 + (y-y0)**2) / fwhm**2)
m_hsize = 15
xshift = 5
yshift = 15
for g_size in [99, 100]:
image = makeGaussian(g_size)
np.testing.assert_array_equal(func.mask_builder(image, m_hsize)[1],
(g_size//2, g_size//2))
image[yshift:, :] = image[:-yshift, :]
image [:yshift, :]= np.zeros((yshift, g_size))
np.testing.assert_array_equal(func.mask_builder(image, m_hsize)[1],
(g_size//2 + yshift, g_size//2))
image[:, xshift:] = image[:, :-xshift]
image [:, :xshift]= np.zeros((g_size, xshift))
np.testing.assert_array_equal(func.mask_builder(image, m_hsize)[1],
(g_size//2 + yshift,
g_size//2 + xshift))
def test_mask_builder_mask(self):
"""Test it it returns expected output (big) mask
NOT SURE: I am kind of testing two things in one test...
"""
def makeGaussian(size, fwhm = 3, center=None):
""" Make a square gaussian kernel.
size is the length of a side of the square
fwhm is full-width-half-maximum, which
can be thought of as an effective radius.
"""
x = np.arange(0, size, 1, float)
y = x[:,np.newaxis]
if center is None:
x0 = y0 = size // 2
else:
x0 = center[0]
y0 = center[1]
return np.exp(-4*np.log(2) * ((x-x0)**2 + (y-y0)**2) / fwhm**2)
for m_half_size in [15, 20, 48]:
for g_size in [99, 100]:
image = makeGaussian(g_size)
[big_mask, mask_centre] = func.mask_builder(image, m_half_size)
# Test if output (big) mask has correct shape (same as input)
np.testing.assert_array_equal(big_mask.shape, image.shape)
# Test if output (big) mask has correct mask size (diameter)
# Horizontal direction
# Left
i1 = (big_mask[mask_centre[0], :] == True).argmax()
# Right
i2 = (big_mask.shape[1]
- (big_mask[mask_centre[0], ::-1] == True).argmax())
self.assertEqual(i2 - i1, 2 * m_half_size + 1)
# Vertical direction
# Left
i1 = (big_mask[:, mask_centre[1]] == True).argmax()
# Right
i2 = (big_mask.shape[0]
- (big_mask[::-1, mask_centre[1]] == True).argmax())
self.assertEqual(i2 - i1, 2 * m_half_size + 1)
#~ class MyTest(unittest.TestCase)
#~
#~ def setUp(self):
......
Markdown is supported
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