From e5e797480a4f61fbd492a40769484530a1ec3655 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Leget Date: Wed, 18 Sep 2024 13:16:14 -0700 Subject: [PATCH] add proto of tet borrowed from other interpolator --- .../lsst/meas/algorithms/gp_interpolation.py | 21 ++ tests/test_gp_interp.py | 191 ++++++++++++++++++ 2 files changed, 212 insertions(+) create mode 100644 tests/test_gp_interp.py diff --git a/python/lsst/meas/algorithms/gp_interpolation.py b/python/lsst/meas/algorithms/gp_interpolation.py index 80ff496f8..1a08fbac5 100644 --- a/python/lsst/meas/algorithms/gp_interpolation.py +++ b/python/lsst/meas/algorithms/gp_interpolation.py @@ -1,3 +1,24 @@ +# This file is part of meas_algorithms. +# +# Developed for the LSST Data Management System. +# This product includes software developed by the LSST Project +# (https://www.lsst.org). +# See the COPYRIGHT file at the top-level directory of this distribution +# for details of code ownership. +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + import numpy as np from lsst.meas.algorithms import CloughTocher2DInterpolatorUtils as ctUtils from lsst.geom import Box2I, Point2I, Extent2I diff --git a/tests/test_gp_interp.py b/tests/test_gp_interp.py new file mode 100644 index 000000000..9ebb31957 --- /dev/null +++ b/tests/test_gp_interp.py @@ -0,0 +1,191 @@ +# This file is part of meas_algorithms. +# +# Developed for the LSST Data Management System. +# This product includes software developed by the LSST Project +# (https://www.lsst.org). +# See the COPYRIGHT file at the top-level directory of this distribution +# for details of code ownership. +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + + +import unittest + +from typing import Iterable +from itertools import product +import numpy as np + +import lsst.utils.tests +import lsst.geom +import lsst.afw.image as afwImage +from lsst.meas.algorithms.cloughTocher2DInterpolator import ( + CloughTocher2DInterpolateTask, +) +from lsst.meas.algorithms import CloughTocher2DInterpolatorUtils as ctUtils + + +class InterpolateOverDefectGaussianProcessTestCase(lsst.utils.tests.TestCase): + """Test InterpolateOverDefectGaussianProcess.""" + + def setUp(self): + super().setUp() + + self.maskedimage = afwImage.MaskedImageF(100, 121) + for x in range(100): + for y in range(121): + self.maskedimage[x, y] = (3 * y + x * 5, 0, 1.0) + + # Clone the maskedimage so we can compare it after running the task. + self.reference = self.maskedimage.clone() + + # Set some central pixels as SAT + sliceX, sliceY = slice(30, 35), slice(40, 45) + self.maskedimage.mask[sliceX, sliceY] = afwImage.Mask.getPlaneBitMask("SAT") + self.maskedimage.image[sliceX, sliceY] = np.nan + # Put nans here to make sure interp is done ok + + # Set an entire column as BAD + self.maskedimage.mask[54:55, :] = afwImage.Mask.getPlaneBitMask("BAD") + self.maskedimage.image[54:55, :] = np.nan + + # Set an entire row as BAD + self.maskedimage.mask[:, 110:111] = afwImage.Mask.getPlaneBitMask("BAD") + self.maskedimage.image[:, 110:111] = np.nan + + # Set a diagonal set of pixels as CR + for i in range(74, 78): + self.maskedimage.mask[i, i] = afwImage.Mask.getPlaneBitMask("CR") + self.maskedimage.image[i, i] = np.nan + + # Set one of the edges as EDGE + self.maskedimage.mask[0:1, :] = afwImage.Mask.getPlaneBitMask("EDGE") + self.maskedimage.image[0:1, :] = np.nan + + # Set a smaller streak at the edge + self.maskedimage.mask[25:28, 0:1] = afwImage.Mask.getPlaneBitMask("EDGE") + self.maskedimage.image[25:28, 0:1] = np.nan + + # Update the reference image's mask alone, so we can compare them after + # running the task. + self.reference.mask.array[:, :] = self.maskedimage.mask.array + + # Create a noise image + self.noise = self.maskedimage.clone() + np.random.seed(12345) + self.noise.image.array[:, :] = np.random.normal(size=self.noise.image.array.shape) + + @lsst.utils.tests.methodParameters(n_runs=(1, 2)) + def test_interpolation(self, n_runs: int): + """Test that the interpolation is done correctly. + + Parameters + ---------- + n_runs : `int` + Number of times to run the task. Running the task more than once + should have no effect. + """ + config = CloughTocher2DInterpolateTask.ConfigClass() + config.badMaskPlanes = ( + "BAD", + "SAT", + "CR", + "EDGE", + ) + config.fillValue = 0.5 + task = CloughTocher2DInterpolateTask(config) + for n in range(n_runs): + task.run(self.maskedimage) + + # Assert that the mask and the variance planes remain unchanged. + self.assertImagesEqual(self.maskedimage.variance, self.reference.variance) + self.assertMasksEqual(self.maskedimage.mask, self.reference.mask) + + # Check that the long streak of bad pixels have been replaced with the + # fillValue, but not the short streak. + np.testing.assert_array_equal(self.maskedimage.image[0:1, :].array, config.fillValue) + with self.assertRaises(AssertionError): + np.testing.assert_array_equal(self.maskedimage.image[25:28, 0:1].array, config.fillValue) + + # Check that interpolated pixels are close to the reference (original), + # and that none of them is still NaN. + self.assertTrue(np.isfinite(self.maskedimage.image.array).all()) + self.assertImagesAlmostEqual( + self.maskedimage.image[1:, :], + self.reference.image[1:, :], + rtol=1e-05, + atol=1e-08, + ) + TEST = 42 + if TEST == 42: + raise ValueError('TEST == 42') + + @lsst.utils.tests.methodParametersProduct(pass_badpix=(True, False), pass_goodpix=(True, False)) + def test_interpolation_with_noise(self, pass_badpix: bool = True, pass_goodpix: bool = True): + """Test that we can reuse the badpix and goodpix. + + Parameters + ---------- + pass_badpix : `bool` + Whether to pass the badpix to the task? + pass_goodpix : `bool` + Whether to pass the goodpix to the task? + """ + + config = CloughTocher2DInterpolateTask.ConfigClass() + config.badMaskPlanes = ( + "BAD", + "SAT", + "CR", + "EDGE", + ) + task = CloughTocher2DInterpolateTask(config) + + badpix, goodpix = task.run(self.noise) + task.run( + self.maskedimage, + badpix=(badpix if pass_badpix else None), + goodpix=(goodpix if pass_goodpix else None), + ) + + # Check that the long streak of bad pixels by the edge have been + # replaced with fillValue, but not the short streak. + np.testing.assert_array_equal(self.maskedimage.image[0:1, :].array, config.fillValue) + with self.assertRaises(AssertionError): + np.testing.assert_array_equal(self.maskedimage.image[25:28, 0:1].array, config.fillValue) + + # Check that interpolated pixels are close to the reference (original), + # and that none of them is still NaN. + self.assertTrue(np.isfinite(self.maskedimage.image.array).all()) + self.assertImagesAlmostEqual( + self.maskedimage.image[1:, :], + self.reference.image[1:, :], + rtol=1e-05, + atol=1e-08, + ) + TEST = 42 + if TEST == 42: + raise ValueError('TEST == 42') + + +def setup_module(module): + lsst.utils.tests.init() + + +class MemoryTestCase(lsst.utils.tests.MemoryTestCase): + pass + + +if __name__ == "__main__": + lsst.utils.tests.init() + unittest.main()