Coverage for tests/test_linearize.py: 10%
191 statements
« prev ^ index » next coverage.py v7.10.6, created at 2025-10-11 01:47 -0700
1# This file is part of ip_isr.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
22import unittest
24import logging
25import numpy as np
27import lsst.utils.tests
28import lsst.utils
29import lsst.afw.image as afwImage
30import lsst.afw.math as afwMath
31import lsst.afw.cameraGeom as cameraGeom
32from lsst.afw.geom.testUtils import BoxGrid
33from lsst.afw.image.testUtils import makeRampImage
34from lsst.ip.isr import applyLookupTable, Linearizer
37def referenceImage(image, detector, linearityType, inputData, table=None):
38 """Generate a reference linearization.
40 Parameters
41 ----------
42 image: `lsst.afw.image.Image`
43 Image to linearize.
44 detector: `lsst.afw.cameraGeom.Detector`
45 Detector this image is from.
46 linearityType: `str`
47 Type of linearity to apply.
48 inputData: `numpy.array`
49 An array of values for the linearity correction.
50 table: `numpy.array`, optional
51 An optional lookup table to use.
53 Returns
54 -------
55 outImage: `lsst.afw.image.Image`
56 The output linearized image.
57 numOutOfRange: `int`
58 The number of values that could not be linearized.
60 Raises
61 ------
62 RuntimeError :
63 Raised if an invalid linearityType is supplied.
64 """
65 numOutOfRange = 0
66 for ampIdx, amp in enumerate(detector.getAmplifiers()):
67 ampIdx = (ampIdx // 3, ampIdx % 3)
68 bbox = amp.getBBox()
69 imageView = image.Factory(image, bbox)
71 if linearityType == 'Squared':
72 sqCoeff = inputData[ampIdx]
73 array = imageView.getArray()
75 array[:] = array + sqCoeff*array**2
76 elif linearityType == 'LookupTable':
77 rowInd, colIndOffset = inputData[ampIdx]
78 rowInd = int(rowInd)
79 tableRow = table[rowInd, :]
80 numOutOfRange += applyLookupTable(imageView, tableRow, colIndOffset)
81 elif linearityType == 'Polynomial':
82 coeffs = inputData[ampIdx]
83 array = imageView.getArray()
84 summation = np.zeros_like(array)
85 for index, coeff in enumerate(coeffs):
86 summation += coeff*np.power(array, (index + 2))
87 array += summation
88 elif linearityType == 'Spline':
89 centers, values = np.split(inputData, 2) # This uses the full data
90 # Note that we are using the slow afw AKIMA_SPLINE interpolator
91 # to offset the data, but using the equivalent but faster scipy
92 # Akima1DInterpolator to correct the data.
93 interp = afwMath.makeInterpolate(centers.tolist(), values.tolist(),
94 afwMath.stringToInterpStyle('AKIMA_SPLINE'))
95 array = imageView.getArray()
96 delta = interp.interpolate(array.flatten())
97 array -= np.array(delta).reshape(array.shape)
98 elif linearityType == 'DoubleSpline':
99 nNodes1 = int(inputData[0])
100 nNodes2 = int(inputData[1])
102 centers1, values1 = np.split(inputData[2: 2 + 2*nNodes1], 2)
103 centers2, values2 = np.split(inputData[2 + 2*nNodes1: 2 + 2*nNodes1 + 2*nNodes2], 2)
104 interp1 = afwMath.makeInterpolate(centers1.tolist(), values1.tolist(),
105 afwMath.stringToInterpStyle('AKIMA_SPLINE'))
106 interp2 = afwMath.makeInterpolate(centers2.tolist(), values2.tolist(),
107 afwMath.stringToInterpStyle('AKIMA_SPLINE'))
108 array = imageView.getArray()
109 delta1 = interp1.interpolate(array.ravel())
110 array -= np.asarray(delta1).reshape(array.shape)
111 delta2 = interp2.interpolate(array.ravel())
112 array -= np.asarray(delta2).reshape(array.shape)
113 else:
114 raise RuntimeError(f"Unknown linearity: {linearityType}")
115 return image, numOutOfRange
118class LinearizeTestCase(lsst.utils.tests.TestCase):
119 """Unit tests for linearizers.
120 """
122 def setUp(self):
123 # This uses the same arbitrary values used in previous tests.
124 self.bbox = lsst.geom.Box2I(lsst.geom.Point2I(-31, 22), lsst.geom.Extent2I(100, 85))
125 self.ampArrangement = (2, 3)
126 self.numAmps = self.ampArrangement[0]*self.ampArrangement[1]
127 # Squared Parameters
128 self.sqCoeffs = np.array([[0, 5e-6, 2.5e-5], [1e-5, 1.1e-6, 2.1e-6]], dtype=float)
130 # Lookup Table Parameters
131 self.colIndOffsets = np.array([[0, -50, 2.5], [37, 1, -3]], dtype=float)
132 self.rowInds = np.array([[0, 1, 4], [3, 5, 2]])
133 # This creates a 2x3 array (matching the amplifiers) that contains a
134 # 2x1 array containing [colIndOffset_i, rowInd_i].
135 self.lookupIndices = np.transpose(np.stack((self.rowInds, self.colIndOffsets), axis=0),
136 axes=[1, 2, 0])
138 rng = np.random.Generator(np.random.MT19937(1))
139 self.table = rng.normal(scale=55, size=(self.numAmps, 2500))
140 self.assertLess(np.max(self.rowInds), self.numAmps, "error in test conditions; invalid row index")
142 # Polynomial Parameters: small perturbation on Squared
143 self.polyCoeffs = np.array([[[0, 1e-7], [5e-6, 1e-7], [2.5e-5, 1e-7]],
144 [[1e-5, 1e-7], [1.1e-6, 1e-7], [2.1e-6, 1e-7]]], dtype=float)
146 # Spline coefficients: should match a 1e-6 Squared solution
147 self.splineCoeffs = np.array([0.0, 1000, 2000, 3000, 4000, 5000,
148 0.0, 1.0, 4.0, 9.0, 16.0, 25.0])
150 # Double spline coefficients.
151 self.doubleSplineCoeffs = np.asarray(
152 [
153 5, # Number of nodes in first spline.
154 6, # Number of nodes in second spline.
155 0.0, 2000., 3000., 4000., 5000., # Nodes for first spline.
156 0.0, 0.01, 0.02, 0.03, 0.04, # Values for first spline.
157 0.0, 1000, 2000, 3000, 4000, 5000, # Nodes for second spline.
158 0.0, 1.0, 4.0, 9.0, 16.0, 25.0, # Values for second spline.
159 0.0, 0.0, # Extra filler.
160 ],
161 )
163 self.log = logging.getLogger("lsst.ip.isr.testLinearizer")
165 def tearDown(self):
166 # destroy LSST objects so memory test passes.
167 self.bbox = None
168 self.detector = None
170 def compareResults(self, linearizedImage, linearizedOutOfRange, linearizedCount, linearizedAmps,
171 referenceImage, referenceOutOfRange, referenceCount, referenceAmps):
172 """Run assert tests on results.
174 Parameters
175 ----------
176 linearizedImage : `lsst.afw.image.Image`
177 Corrected image.
178 linearizedOutOfRange : `int`
179 Number of measured out-of-range pixels.
180 linearizedCount : `int`
181 Number of amplifiers that should be linearized.
182 linearizedAmps : `int`
183 Total number of amplifiers checked.
184 referenceImage : `lsst.afw.image.Image`
185 Truth image to compare against.
186 referenceOutOfRange : `int`
187 Number of expected out-of-range-pixels.
188 referenceCount : `int`
189 Number of amplifiers that are expected to be linearized.
190 referenceAmps : `int`
191 Expected number of amplifiers checked.
192 """
193 self.assertImagesAlmostEqual(linearizedImage, referenceImage)
194 self.assertEqual(linearizedOutOfRange, referenceOutOfRange)
195 self.assertEqual(linearizedCount, referenceCount)
196 self.assertEqual(linearizedAmps, referenceAmps)
198 def testBasics(self):
199 """Test basic linearization functionality.
200 """
201 for imageClass in (afwImage.ImageF, afwImage.ImageD):
202 inImage = makeRampImage(bbox=self.bbox, start=-5, stop=2500, imageClass=imageClass)
204 for linearityType in ('Squared', 'LookupTable', 'Polynomial', 'Spline', 'DoubleSpline'):
205 detector = self.makeDetector(linearityType)
206 table = None
207 inputData = {'Squared': self.sqCoeffs,
208 'LookupTable': self.lookupIndices,
209 'Polynomial': self.polyCoeffs,
210 'Spline': self.splineCoeffs,
211 'DoubleSpline': self.doubleSplineCoeffs}[linearityType]
212 if linearityType == 'LookupTable':
213 table = np.array(self.table, dtype=inImage.getArray().dtype)
214 linearizer = Linearizer(detector=detector, table=table)
216 measImage = inImage.Factory(inImage, True)
217 result = linearizer.applyLinearity(measImage, detector=detector, log=self.log)
218 refImage, refNumOutOfRange = referenceImage(inImage.Factory(inImage, True),
219 detector, linearityType, inputData, table)
221 # This is necessary for the same tests to be used on
222 # all types. The first amplifier has 0.0 for the
223 # coefficient, which should be tested (it has a log
224 # message), but we are not linearizing an amplifier
225 # with no correction, so it fails the test that
226 # numLinearized == numAmps.
227 zeroLinearity = 1 if linearityType == 'Squared' else 0
229 self.compareResults(measImage, result.numOutOfRange, result.numLinearized, result.numAmps,
230 refImage, refNumOutOfRange, self.numAmps - zeroLinearity, self.numAmps)
232 # Test a stand alone linearizer. This ignores validate checks.
233 measImage = inImage.Factory(inImage, True)
234 storedLinearizer = self.makeLinearizer(linearityType)
235 storedResult = storedLinearizer.applyLinearity(measImage, log=self.log)
237 self.compareResults(measImage, storedResult.numOutOfRange, storedResult.numLinearized,
238 storedResult.numAmps,
239 refImage, refNumOutOfRange, self.numAmps - zeroLinearity, self.numAmps)
241 # "Save to yaml" and test again
242 storedDict = storedLinearizer.toDict()
243 storedLinearizer = Linearizer().fromDict(storedDict)
245 measImage = inImage.Factory(inImage, True)
246 storedResult = storedLinearizer.applyLinearity(measImage, log=self.log)
248 self.compareResults(measImage, storedResult.numOutOfRange, storedResult.numLinearized,
249 storedResult.numAmps,
250 refImage, refNumOutOfRange, self.numAmps - zeroLinearity, self.numAmps)
252 # "Save to fits" and test again
253 storedTable = storedLinearizer.toTable()
254 storedLinearizer = Linearizer().fromTable(storedTable)
256 measImage = inImage.Factory(inImage, True)
257 storedResult = storedLinearizer.applyLinearity(measImage, log=self.log)
259 self.compareResults(measImage, storedResult.numOutOfRange, storedResult.numLinearized,
260 storedResult.numAmps,
261 refImage, refNumOutOfRange, self.numAmps - zeroLinearity, self.numAmps)
263 # Use a gain and test again
264 measImage = inImage.Factory(inImage, True)
265 storedLinearizer = self.makeLinearizer(linearityType)
266 gains = {key: 1.0 for key in storedLinearizer.linearityType.keys()}
267 storedResult = storedLinearizer.applyLinearity(measImage, log=self.log, gains=gains)
269 self.compareResults(measImage, storedResult.numOutOfRange, storedResult.numLinearized,
270 storedResult.numAmps,
271 refImage, refNumOutOfRange, self.numAmps - zeroLinearity, self.numAmps)
273 def makeDetector(self, linearityType, bbox=None):
274 """Generate a fake detector for the test.
276 Parameters
277 ----------
278 linearityType : `str`
279 Which linearity to assign to the detector's cameraGeom.
280 bbox : `lsst.geom.Box2I`, optional
281 Bounding box to use for the detector.
283 Returns
284 -------
285 detBuilder : `lsst.afw.cameraGeom.Detector`
286 The fake detector.
287 """
288 bbox = bbox if bbox is not None else self.bbox
289 numAmps = self.ampArrangement
291 detName = "det_a"
292 detId = 1
293 detSerial = "123"
294 orientation = cameraGeom.Orientation()
295 pixelSize = lsst.geom.Extent2D(1, 1)
297 camBuilder = cameraGeom.Camera.Builder("fakeCam")
298 detBuilder = camBuilder.add(detName, detId)
299 detBuilder.setSerial(detSerial)
300 detBuilder.setBBox(bbox)
301 detBuilder.setOrientation(orientation)
302 detBuilder.setPixelSize(pixelSize)
304 boxArr = BoxGrid(box=bbox, numColRow=numAmps)
305 for i in range(numAmps[0]):
306 for j in range(numAmps[1]):
307 ampInfo = cameraGeom.Amplifier.Builder()
308 ampInfo.setName("amp %d_%d" % (i + 1, j + 1))
309 ampInfo.setBBox(boxArr[i, j])
310 ampInfo.setLinearityType(linearityType)
311 if linearityType == 'Squared':
312 ampInfo.setLinearityCoeffs([self.sqCoeffs[i, j]])
313 elif linearityType == 'LookupTable':
314 # setLinearityCoeffs is picky about getting a mixed
315 # int/float list.
316 ampInfo.setLinearityCoeffs(np.array([self.rowInds[i, j], self.colIndOffsets[i, j],
317 0, 0], dtype=float))
318 elif linearityType == 'Polynomial':
319 ampInfo.setLinearityCoeffs(self.polyCoeffs[i, j])
320 elif linearityType == 'Spline':
321 ampInfo.setLinearityCoeffs(self.splineCoeffs)
322 elif linearityType == 'DoubleSpline':
323 ampInfo.setLinearityCoeffs(self.doubleSplineCoeffs)
324 detBuilder.append(ampInfo)
326 return detBuilder
328 def makeLinearizer(self, linearityType, bbox=None):
329 """Construct a linearizer with the test coefficients.
331 Parameters
332 ----------
333 linearityType : `str`
334 Type of linearity to use. The coefficients are set by the
335 setUp method.
336 bbox : `lsst.geom.Box2I`
337 Bounding box for the full detector. Used to assign
338 amp-based bounding boxes.
340 Returns
341 -------
342 linearizer : `lsst.ip.isr.Linearizer`
343 A fully constructed, persistable linearizer.
344 """
345 bbox = bbox if bbox is not None else self.bbox
346 numAmps = self.ampArrangement
347 boxArr = BoxGrid(box=bbox, numColRow=numAmps)
348 linearizer = Linearizer()
349 linearizer.hasLinearity = True
351 for i in range(numAmps[0]):
352 for j in range(numAmps[1]):
353 ampName = f"amp {i+1}_{j+1}"
354 ampBox = boxArr[i, j]
355 linearizer.ampNames.append(ampName)
357 if linearityType == 'Squared':
358 linearizer.linearityCoeffs[ampName] = np.array([self.sqCoeffs[i, j]])
359 elif linearityType == 'LookupTable':
360 linearizer.linearityCoeffs[ampName] = np.array(self.lookupIndices[i, j])
361 linearizer.tableData = self.table
362 elif linearityType == 'Polynomial':
363 linearizer.linearityCoeffs[ampName] = np.array(self.polyCoeffs[i, j])
364 elif linearityType == 'Spline':
365 linearizer.linearityCoeffs[ampName] = np.array(self.splineCoeffs)
366 elif linearityType == 'DoubleSpline':
367 linearizer.linearityCoeffs[ampName] = np.asarray(self.doubleSplineCoeffs)
369 linearizer.linearityType[ampName] = linearityType
370 linearizer.linearityBBox[ampName] = ampBox
371 linearizer.inputAbscissa[ampName] = np.array([])
372 linearizer.inputOrdinate[ampName] = np.array([])
373 linearizer.inputMask[ampName] = np.array([], dtype=np.bool_)
374 linearizer.inputGroupingIndex[ampName] = np.array([], dtype=np.int64)
375 linearizer.inputNormalization[ampName] = np.array([])
376 linearizer.fitParams[ampName] = np.array([])
377 linearizer.fitParamsErr[ampName] = np.array([])
378 linearizer.fitChiSq[ampName] = np.nan
379 linearizer.fitResiduals[ampName] = np.array([])
380 linearizer.fitResidualsSigmaMad[ampName] = np.nan
381 linearizer.fitResidualsUnmasked[ampName] = np.array([])
382 linearizer.fitResidualsModel[ampName] = np.array([])
383 linearizer.linearFit[ampName] = np.array([])
384 linearizer.linearityTurnoff[ampName] = np.nan
385 linearizer.linearityMaxSignal[ampName] = np.nan
387 return linearizer
390class MemoryTester(lsst.utils.tests.MemoryTestCase):
391 pass
394def setup_module(module):
395 lsst.utils.tests.init()
398if __name__ == "__main__": 398 ↛ 399line 398 didn't jump to line 399 because the condition on line 398 was never true
399 lsst.utils.tests.init()
400 unittest.main()