Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Make a bunch of coalignment module private #100

Merged
merged 8 commits into from
Sep 22, 2022
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
11 changes: 11 additions & 0 deletions changelog/100.breaking.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1,11 @@
The following helper functions in `sunkit_image.colaginment` have been removed, with no replacement.
This is because they are designed to be internal helper functions.
If you need to use them in your own code create a copy of the functions from the ``sunkit-image`` source code.

- ``parabolic_turning_point``
- ``calculate_clipping``
- ``check_for_nonfinite_entries``
- ``get_correlation_shifts``
- ``clip_edges``
- ``find_best_match_location``
- ``calculate_shift``
38 changes: 16 additions & 22 deletions sunkit_image/coalignment.py
Original file line number Diff line number Diff line change
Expand Up @@ -19,14 +19,7 @@
__author__ = "J. Ireland"

__all__ = [
"calculate_shift",
"clip_edges",
"calculate_clipping",
"match_template_to_layer",
"find_best_match_location",
"get_correlation_shifts",
"parabolic_turning_point",
"check_for_nonfinite_entries",
"apply_shifts",
"mapsequence_coalign_by_match_template",
"calculate_match_template_shift",
Expand All @@ -45,7 +38,7 @@ def _default_fmap_function(data):
return np.float64(data)


def calculate_shift(this_layer, template):
def _calculate_shift(this_layer, template):
"""
Calculates the pixel shift required to put the template in the "best"
position on a layer.
Expand All @@ -65,15 +58,15 @@ def calculate_shift(this_layer, template):
to the input array.
"""
# Warn user if any NANs, Infs, etc are present in the layer or the template
check_for_nonfinite_entries(this_layer, template)
_check_for_nonfinite_entries(this_layer, template)
# Calculate the correlation array matching the template to this layer
corr = match_template_to_layer(this_layer, template)
# Calculate the y and x shifts in pixels
return find_best_match_location(corr)
return _find_best_match_location(corr)


@u.quantity_input
def clip_edges(data, yclips: u.pix, xclips: u.pix):
def _clip_edges(data, yclips: u.pix, xclips: u.pix):
"""
Clips off the "y" and "x" edges of a 2D array according to a list of pixel
values. This function is useful for removing data at the edge of 2d images
Expand Down Expand Up @@ -105,7 +98,7 @@ def clip_edges(data, yclips: u.pix, xclips: u.pix):


@u.quantity_input
def calculate_clipping(y: u.pix, x: u.pix):
def _calculate_clipping(y: u.pix, x: u.pix):
"""
Return the upper and lower clipping values for the "y" and "x" directions.

Expand Down Expand Up @@ -184,7 +177,7 @@ def match_template_to_layer(layer, template):
return match_template(layer, template)


def find_best_match_location(corr):
def _find_best_match_location(corr):
"""
Calculate an estimate of the location of the peak of the correlation result
in image pixels.
Expand All @@ -209,7 +202,7 @@ def find_best_match_location(corr):
np.max([0, cor_max_y - 1]) : np.min([cor_max_y + 2, corr.shape[0] - 1]),
np.max([0, cor_max_x - 1]) : np.min([cor_max_x + 2, corr.shape[1] - 1]),
]
y_shift_maximum, x_shift_maximum = get_correlation_shifts(array_maximum)
y_shift_maximum, x_shift_maximum = _get_correlation_shifts(array_maximum)

# Get shift relative to correlation array
y_shift_correlation_array = y_shift_maximum + cor_max_y * u.pix
Expand All @@ -218,7 +211,7 @@ def find_best_match_location(corr):
return y_shift_correlation_array, x_shift_correlation_array


def get_correlation_shifts(array):
def _get_correlation_shifts(array):
"""
Estimate the location of the maximum of a fit to the input array. The
estimation in the "x" and "y" directions are done separately. The location
Expand Down Expand Up @@ -251,19 +244,19 @@ def get_correlation_shifts(array):
# Otherwise, just return the location of the maximum in a particular
# direction.
if ny == 3:
y_location = parabolic_turning_point(array[:, x_max_location])
y_location = _parabolic_turning_point(array[:, x_max_location])
else:
y_location = 1.0 * y_max_location

if nx == 3:
x_location = parabolic_turning_point(array[y_max_location, :])
x_location = _parabolic_turning_point(array[y_max_location, :])
else:
x_location = 1.0 * x_max_location

return y_location * u.pix, x_location * u.pix


def parabolic_turning_point(y):
def _parabolic_turning_point(y):
"""
Find the location of the turning point for a parabola
``y(x) = ax^2 + bx + c``, given input values ``y(-1), y(0), y(1)``.
Expand All @@ -287,7 +280,7 @@ def parabolic_turning_point(y):
return numerator / denominator


def check_for_nonfinite_entries(layer_image, template_image):
def _check_for_nonfinite_entries(layer_image, template_image):
"""
Issue a warning if there is any nonfinite entry in the layer or template
images.
Expand Down Expand Up @@ -355,7 +348,7 @@ def apply_shifts(mc, yshift: u.pix, xshift: u.pix, clip=True, **kwargs):

# Calculate the clipping
if clip:
yclips, xclips = calculate_clipping(-yshift, -xshift)
yclips, xclips = _calculate_clipping(-yshift, -xshift)

# Shift the data and construct the mapsequence
for i, m in enumerate(mc):
Expand All @@ -364,7 +357,7 @@ def apply_shifts(mc, yshift: u.pix, xshift: u.pix, clip=True, **kwargs):
# Clip if required. Use the submap function to return the appropriate
# portion of the data.
if clip:
shifted_data = clip_edges(shifted_data, yclips, xclips)
shifted_data = _clip_edges(shifted_data, yclips, xclips)
new_meta["naxis1"] = shifted_data.shape[1]
new_meta["naxis2"] = shifted_data.shape[0]
# Add one to go from zero-based to one-based indexing
Expand Down Expand Up @@ -448,7 +441,7 @@ def calculate_match_template_shift(mc, template=None, layer_index=0, func=_defau
this_layer = func(m.data)

# Calculate the y and x shifts in pixels
yshift, xshift = calculate_shift(this_layer, tplate)
yshift, xshift = _calculate_shift(this_layer, tplate)

# Keep shifts in pixels
yshift_keep[i] = yshift
Expand Down Expand Up @@ -605,6 +598,7 @@ def calculate_solar_rotate_shift(mc, layer_index=0, **kwargs):
``**kwargs``
These keywords are passed to the function
`sunpy.physics.differential_rotation.solar_rotate_coordinate`.

Returns
-------
`~astropy.units.Quantity`, ~astropy.units.Quantity`
Expand Down
44 changes: 22 additions & 22 deletions sunkit_image/tests/test_coalignment.py
Original file line number Diff line number Diff line change
Expand Up @@ -15,21 +15,21 @@
from sunpy.util import SunpyUserWarning

from sunkit_image.coalignment import (
_calculate_clipping,
_check_for_nonfinite_entries,
_clip_edges,
_default_fmap_function,
_find_best_match_location,
_get_correlation_shifts,
_lower_clip,
_parabolic_turning_point,
_upper_clip,
apply_shifts,
calculate_clipping,
calculate_match_template_shift,
calculate_solar_rotate_shift,
check_for_nonfinite_entries,
clip_edges,
find_best_match_location,
get_correlation_shifts,
mapsequence_coalign_by_match_template,
mapsequence_coalign_by_rotation,
match_template_to_layer,
parabolic_turning_point,
)


Expand Down Expand Up @@ -75,14 +75,14 @@ def aia171_test_template_shape(aia171_test_template):


def test_parabolic_turning_point():
assert parabolic_turning_point(np.asarray([6.0, 2.0, 0.0])) == 1.5
assert _parabolic_turning_point(np.asarray([6.0, 2.0, 0.0])) == 1.5


def test_check_for_nonfinite_entries():
with warnings.catch_warnings(record=True) as warning_list:
a = np.zeros((3, 3))
b = np.ones((3, 3))
check_for_nonfinite_entries(a, b)
_check_for_nonfinite_entries(a, b)

assert len(warning_list) == 0

Expand All @@ -93,17 +93,17 @@ def test_check_for_nonfinite_entries():
b = a.reshape(3, 3)

with pytest.warns(SunpyUserWarning, match="The layer image has nonfinite entries.") as warning_list:
check_for_nonfinite_entries(b, np.ones((3, 3)))
_check_for_nonfinite_entries(b, np.ones((3, 3)))

assert len(warning_list) == 1

with pytest.warns(SunpyUserWarning, match="The template image has nonfinite entries.") as warning_list:
check_for_nonfinite_entries(np.ones((3, 3)), b)
_check_for_nonfinite_entries(np.ones((3, 3)), b)

assert len(warning_list) == 1

with pytest.warns(Warning) as warning_list:
check_for_nonfinite_entries(b, b)
_check_for_nonfinite_entries(b, b)

assert len(warning_list) == 2

Expand All @@ -130,7 +130,7 @@ def test_get_correlation_shifts():
test_array[1, 1] = 1
test_array[2, 1] = 0.6
test_array[1, 2] = 0.2
y_test, x_test = get_correlation_shifts(test_array)
y_test, x_test = _get_correlation_shifts(test_array)
assert_allclose(y_test.value, 0.214285714286, rtol=1e-2, atol=0)
assert_allclose(x_test.value, 0.0555555555556, rtol=1e-2, atol=0)

Expand All @@ -140,23 +140,23 @@ def test_get_correlation_shifts():
test_array[0, 1] = 0.2
test_array[1, 0] = 0.4
test_array[1, 1] = 0.3
y_test, x_test = get_correlation_shifts(test_array)
y_test, x_test = _get_correlation_shifts(test_array)
assert_allclose(y_test.value, 1.0, rtol=1e-2, atol=0)
assert_allclose(x_test.value, 0.0, rtol=1e-2, atol=0)

# Input array is too big in either direction
test_array = np.zeros((4, 3))
with pytest.raises(ValueError):
get_correlation_shifts(test_array)
_get_correlation_shifts(test_array)

test_array = np.zeros((3, 4))
with pytest.raises(ValueError):
get_correlation_shifts(test_array)
_get_correlation_shifts(test_array)


def test_find_best_match_location(aia171_test_map_layer, aia171_test_template, aia171_test_shift):
result = match_template_to_layer(aia171_test_map_layer, aia171_test_template)
match_location = u.Quantity(find_best_match_location(result))
match_location = u.Quantity(_find_best_match_location(result))
assert_allclose(match_location.value, np.array(result.shape) / 2.0 - 0.5 + aia171_test_shift, rtol=1e-3, atol=0)


Expand All @@ -175,15 +175,15 @@ def test_upper_clip(aia171_test_clipping):


def test_calculate_clipping(aia171_test_clipping):
answer = calculate_clipping(aia171_test_clipping * u.pix, aia171_test_clipping * u.pix)
answer = _calculate_clipping(aia171_test_clipping * u.pix, aia171_test_clipping * u.pix)
assert_array_almost_equal(answer, ([2.0, 1.0] * u.pix, [2.0, 1.0] * u.pix))


def test_clip_edges():
a = np.zeros(shape=(341, 156))
yclip = [4, 0] * u.pix
xclip = [1, 2] * u.pix
new_a = clip_edges(a, yclip, xclip)
new_a = _clip_edges(a, yclip, xclip)
assert a.shape[0] - (yclip[0].value + yclip[1].value) == new_a.shape[0]
assert a.shape[1] - (xclip[0].value + xclip[1].value) == new_a.shape[1]

Expand Down Expand Up @@ -281,7 +281,7 @@ def test_mapsequence_coalign_by_match_template(aia171_test_mc, aia171_test_map_l
test_mc = mapsequence_coalign_by_match_template(aia171_test_mc)
x_displacement_pixels = test_displacements["x"] / test_mc[0].scale[0]
y_displacement_pixels = test_displacements["y"] / test_mc[0].scale[1]
expected_clipping = calculate_clipping(y_displacement_pixels, x_displacement_pixels)
expected_clipping = _calculate_clipping(y_displacement_pixels, x_displacement_pixels)
number_of_pixels_clipped = [np.sum(np.abs(expected_clipping[0])), np.sum(np.abs(expected_clipping[1]))]

assert test_mc[0].data.shape == (
Expand All @@ -299,7 +299,7 @@ def test_mapsequence_coalign_by_match_template(aia171_test_mc, aia171_test_map_l
test_mc = mapsequence_coalign_by_match_template(aia171_test_mc, clip=True)
x_displacement_pixels = test_displacements["x"] / test_mc[0].scale[0]
y_displacement_pixels = test_displacements["y"] / test_mc[0].scale[1]
expected_clipping = calculate_clipping(y_displacement_pixels, x_displacement_pixels)
expected_clipping = _calculate_clipping(y_displacement_pixels, x_displacement_pixels)
number_of_pixels_clipped = [np.sum(np.abs(expected_clipping[0])), np.sum(np.abs(expected_clipping[1]))]

assert test_mc[0].data.shape == (
Expand Down Expand Up @@ -357,7 +357,7 @@ def test_apply_shifts(aia171_test_map):
# original layer by a known amount.
test_mc = apply_shifts(mc, astropy_displacements["y"], astropy_displacements["x"], clip=True)
for i in range(0, len(test_mc.maps)):
clipped = calculate_clipping(astropy_displacements["y"], astropy_displacements["x"])
clipped = _calculate_clipping(astropy_displacements["y"], astropy_displacements["x"])
assert test_mc[i].data.shape[0] == mc[i].data.shape[0] - np.max(clipped[0].value)
assert test_mc[i].data.shape[1] == mc[i].data.shape[1] - np.max(clipped[1].value)

Expand All @@ -366,7 +366,7 @@ def test_apply_shifts(aia171_test_map):
# than the original layer by a known amount.
test_mc = apply_shifts(mc, astropy_displacements["y"], astropy_displacements["x"])
for i in range(0, len(test_mc.maps)):
clipped = calculate_clipping(astropy_displacements["y"], astropy_displacements["x"])
clipped = _calculate_clipping(astropy_displacements["y"], astropy_displacements["x"])
assert test_mc[i].data.shape[0] == mc[i].data.shape[0] - np.max(clipped[0].value)
assert test_mc[i].data.shape[1] == mc[i].data.shape[1] - np.max(clipped[1].value)

Expand Down