From 34ad356b9e4d94f92cacc9b5b00416f09b69cb94 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Sun, 11 Oct 2026 14:46:15 -0600 Subject: [PATCH 1/2] Add `Filtergram.ratio_noise()`, which measures the noise of a sequence of images from the difference between each image and the next. The measurement of the DEM tutorial, as a method, so that the DEM of every frame of a movie can scale its uncertainty the same way. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_011VGPHLGJi8acbuLYSfBHTk --- hinode/xrt/_filtergrams.py | 123 ++++++++++++++++++++++++++++++++ hinode/xrt/_filtergrams_test.py | 93 ++++++++++++++++++++++++ 2 files changed, 216 insertions(+) diff --git a/hinode/xrt/_filtergrams.py b/hinode/xrt/_filtergrams.py index 6092d6c..a66bf10 100644 --- a/hinode/xrt/_filtergrams.py +++ b/hinode/xrt/_filtergrams.py @@ -931,3 +931,126 @@ def uncertainty( variance = variance + np.square(self.uncertainty_map) return typing.cast(na.ScalarArray, np.sqrt(variance)) + + def ratio_noise( + self, + noise_photon: u.Quantity | na.AbstractScalar, + where: bool | na.AbstractScalar = True, + threshold: u.Quantity = 10 * u.DN / u.s, + gap: u.Quantity = 25 * u.s, + clip: float = 5, + num_iteration: int = 3, + ) -> na.ScalarArray: + r""" + The ratio of the variance of the noise measured in this sequence of + images to the variance of :meth:`uncertainty`. + + The noise is measured from the difference between each image and the + next. + Where the Sun barely changes between the two, + the difference is the noise of both images, + and dividing it by their uncertainties, + + .. math:: + + z = \frac{I_2 - I_1}{\sqrt{\sigma_1^2 + \sigma_2^2}}, + + gives a quantity whose variance is the ratio. + Only images taken within `gap` of each other with the same exposure + time are compared, + and only the pixels selected by `where` which are fainter than + `threshold` in both, + since brighter ones are more likely to change. + Cosmic rays are left out by clipping :math:`z` at `clip` times the + square root of the ratio, which is repeated `num_iteration` times. + + The noise of the images of XRT is less than :meth:`uncertainty` + predicts, + since the charge of each photon diffuses into the neighboring pixels + and the compression of the images removes some of the noise, + so this ratio can scale :meth:`uncertainty` to the noise measured. + It is measured in faint pixels, + and compression removes less of the noise of brighter ones. + + Parameters + ---------- + noise_photon + The factor which sets the photon noise, as in :meth:`uncertainty`. + where + The pixels to measure the noise in, which should leave out any + which change between images, such as those of a jet. + threshold + The signal above which pixels are left out. + gap + The longest time between the starts of two images compared. + clip + How many times the measured noise a difference may be before it + is left out. + num_iteration + How many times to clip the differences. + + Examples + -------- + + The ratio for the Al_poly images captured while ESIS was observing + the Sun, with the photon noise of plasma at 2 MK. + + .. jupyter-execute:: + + import numpy as np + import astropy.units as u + import hinode + + images = hinode.xrt.open( + time_start="2019-09-30T18:06:11", + time_stop="2019-09-30T18:11:01", + leak=True, + uncertainty=True, + ).remove_leak() + + time = images.inputs.time.ndarray[0] + response = hinode.xrt.temperature_response("Al_poly", time) + variance = hinode.xrt.temperature_response("Al_poly", time, variance=True) + + # The temperature closest to 2 MK + index = dict( + filter=0, + temperature=int(np.argmin(np.abs(response.inputs.ndarray - 2 * u.MK))), + ) + noise_photon = (variance.outputs / response.outputs)[index] + + images.ratio_noise(noise_photon) + """ + axis = self.axis_time + before = {axis: slice(None, ~0)} + after = {axis: slice(1, None)} + + intensity = self.outputs + uncertainty = self.uncertainty(noise_photon) + + time = typing.cast(na.ScalarArray, self.inputs.time) + t = typing.cast(astropy.time.Time, time.ndarray) + mjd = typing.cast(np.ndarray, t.mjd) + interval = na.ScalarArray(np.diff(mjd) * u.day, axes=time.axes) + change = np.abs(self.timedelta[after] - self.timedelta[before]) + consecutive = (interval < gap) & (change < 0.01 * u.s) + + quiet = ( + consecutive + & where + & (intensity[before] < threshold) + & (intensity[after] < threshold) + ) + + difference = intensity[after] - intensity[before] + z = difference / np.sqrt( + np.square(uncertainty[after]) + np.square(uncertainty[before]) + ) + z = z.to(u.dimensionless_unscaled).value + + ratio = np.square(z).mean(where=quiet) + for _ in range(num_iteration): + clipped = quiet & (np.abs(z) < clip * np.sqrt(ratio)) + ratio = np.square(z).mean(where=clipped) + + return typing.cast(na.ScalarArray, ratio) diff --git a/hinode/xrt/_filtergrams_test.py b/hinode/xrt/_filtergrams_test.py index 7346dcf..0e1e9ac 100644 --- a/hinode/xrt/_filtergrams_test.py +++ b/hinode/xrt/_filtergrams_test.py @@ -777,3 +777,96 @@ def test_from_time_range_no_images() -> None: time_start="2019-09-30T18:08:01", time_stop="2019-09-30T18:08:02", ) + + +def _sequence( + ratio: float, + num: int = 4, + signal: float = 3, + seed: int = 0, +) -> tuple[hinode.xrt.Filtergram, u.Quantity]: + """ + A sequence of `num` flat images of the given `signal`, in DN per second, + with noise whose variance is `ratio` times that of + :meth:`hinode.xrt.Filtergram.uncertainty`, and the photon noise factor + the uncertainty was computed with. + """ + path = _path("2019-09-30T18:08:30", "2019-09-30T18:08:40") + paths = na.ScalarArray(np.array([path] * num, dtype=object), axes="time") + images = hinode.xrt.Filtergram.from_fits(paths, uncertainty=True) + noise_photon = 2.5 * u.DN / u.ph + flat = na.ScalarArray.ones(images.outputs.shape) * signal * u.DN / u.s + images = dataclasses.replace(images, outputs=flat) + sigma = images.uncertainty(noise_photon) + noise = na.random.normal(0, 1, shape_random=flat.shape, seed=seed) + outputs = flat + np.sqrt(ratio) * sigma * noise + return dataclasses.replace(images, outputs=outputs), noise_photon + + +def test_ratio_noise() -> None: + """The ratio is the variance of the noise of the images over their uncertainty.""" + images, noise_photon = _sequence(ratio=0.2) + result = images.ratio_noise(noise_photon) + assert result.shape == dict() + assert np.isclose(result.ndarray, 0.2, rtol=0.02) + + +def test_ratio_noise_pairs() -> None: + """ + Images are compared only with the next one, and only if it was taken soon + after with the same exposure. + """ + images, noise_photon = _sequence(ratio=0.2, num=5) + noisy, _ = _sequence(ratio=50, num=5, seed=1) + + time = images.inputs.time.ndarray + [0, 18, 36, 54, 100] * u.s + timedelta = na.ScalarArray([16.4, 16.4, 11.6, 16.4, 16.4] * u.s, axes="time") + # Only the first two images are compared: the third has another + # exposure than its neighbors, and the last begins long after the fourth. + outputs = images.outputs.copy() + for i in (2, 4): + outputs[dict(time=i)] = noisy.outputs[dict(time=i)] + + images = dataclasses.replace( + images, + inputs=images.inputs.replace(time=na.ScalarArray(time, axes="time")), + outputs=outputs, + timedelta=timedelta, + ) + + result = images.ratio_noise(noise_photon) + assert np.isclose(result.ndarray, 0.2, rtol=0.05) + + # a longer gap lets the last two be compared + assert images.ratio_noise(noise_photon, gap=60 * u.s).ndarray > 1 + + +def test_ratio_noise_excluded() -> None: + """ + Cosmic rays are clipped, and pixels which are bright or not selected are + left out. + """ + images, noise_photon = _sequence(ratio=0.2) + bright, _ = _sequence(ratio=2, signal=12, seed=1) + + outputs = images.outputs.copy() + # cosmic rays + outputs[dict(time=1, detector_x=slice(10, 13))] = 9 * u.DN / u.s + # a region brighter than the threshold, and noisier + region = dict(detector_x=slice(100, 150)) + outputs[region] = bright.outputs[region] + # a region which brightens a little, left out by `where` + region = dict(time=3, detector_x=slice(200, 250)) + outputs[region] = outputs[region] + 1 * u.DN / u.s + images = dataclasses.replace(images, outputs=outputs) + + index_x = na.arange(0, images.outputs.shape["detector_x"], axis="detector_x") + where = (index_x < 200) | (index_x >= 250) + + result = images.ratio_noise(noise_photon, where=where) + assert np.isclose(result.ndarray, 0.2, rtol=0.02) + + assert images.ratio_noise(noise_photon).ndarray > 0.22 + threshold = 100 * u.DN / u.s + assert images.ratio_noise(noise_photon, where, threshold).ndarray > 0.22 + assert images.ratio_noise(noise_photon, where, num_iteration=0).ndarray > 0.22 From 9d04f8be3271f0d4cc4f24567653acd7dced1edf Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Mon, 12 Oct 2026 10:57:28 -0600 Subject: [PATCH 2/2] Fix the review of `Filtergram.ratio_noise()` - Compare an image with the next only if the next begins after it, so images out of time order or loaded twice are not compared, and only if both have the same plate scale and point within a pixel of each other. - Leave out pixels whose `z` is not finite, saturated pixels, and raise a `ValueError` when there are fewer than two images or no pixel to measure the noise in, instead of returning NaN. - Select pixels by the mean of the two images instead of requiring both below `threshold`, which left out the largest differences: for noise of 0.2 at the threshold that gave 0.08, and now gives 0.21. - Start the clipping from the median of z^2, so outliers in many pixels cannot raise the first estimate enough to survive it. - Reduce over the time and detector axes only, giving a ratio for each entry of the other axes, and accept a `where` per image and a scalar `timedelta`. - Take the intervals by named indexing, find the variance once with a shared `_variance()`, and note in the docstring that the errors of the dark and vignetting models cancel in the difference. - The tutorial calls the method, and its example leaves out the jet. On event E the ratio is 0.1527 (0.1525 before). Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_011VGPHLGJi8acbuLYSfBHTk --- docs/tutorials/dem.rst | 61 +++++----- hinode/xrt/_filtergrams.py | 173 +++++++++++++++++++++++----- hinode/xrt/_filtergrams_test.py | 196 ++++++++++++++++++++++++++++++-- 3 files changed, 358 insertions(+), 72 deletions(-) diff --git a/docs/tutorials/dem.rst b/docs/tutorials/dem.rst index 266a98f..45d46fb 100644 --- a/docs/tutorials/dem.rst +++ b/docs/tutorials/dem.rst @@ -312,7 +312,8 @@ near the temperature of most of the plasma in the box. That is the noise of a CCD which collects all the charge of each photon in the pixel it was absorbed in, before any compression. The noise of these images can be measured instead, -from the difference between each image and the next. +from the difference between each image and the next, +with :meth:`~hinode.xrt.Filtergram.ratio_noise`. Away from the jet the Sun barely changes in 18 seconds, so the difference is the noise of the two images, and dividing it by their uncertainties gives a quantity whose variance is the @@ -324,50 +325,44 @@ and so are cosmic rays, by clipping at five times the noise. frames = images.remove_leak() - intensity_frames = frames.outputs - uncertainty_frames = frames.uncertainty(noise_photon) - - # each image and the one after it - before = dict(time=slice(None, ~0)) - after = dict(time=slice(1, None)) - - # pairs of images taken one after the other with the same exposure - time_frames = images.inputs.time.ndarray - gap = na.ScalarArray((time_frames[1:] - time_frames[:-1]).to(u.s), axes="time") - change = np.abs(images.timedelta[after] - images.timedelta[before]) - consecutive = (gap < 25 * u.s) & (change < 0.01 * u.s) - in_box = na.ScalarArray.zeros(images.outputs[dict(time=0)].shape) > 0 in_box[box] = True - bright = 10 * u.DN / u.s - quiet = ( - consecutive - & ~in_box - & (intensity_frames[before] < bright) - & (intensity_frames[after] < bright) - ) - - difference = intensity_frames[after] - intensity_frames[before] - z = difference / np.sqrt( - np.square(uncertainty_frames[after]) + np.square(uncertainty_frames[before]) - ) - z = z.to(u.dimensionless_unscaled).value - - ratio_noise = np.square(z).mean(where=quiet) - for _ in range(3): - clipped = quiet & (np.abs(z) < 5 * np.sqrt(ratio_noise)) - ratio_noise = np.square(z).mean(where=clipped) + ratio_noise = frames.ratio_noise(noise_photon, where=~in_box) ratio_noise The measured variance is only about a sixth of the prediction. The noise of neighboring pixels is also correlated, more strongly inside the 8 by 8 blocks of the JPEG compression of these images -than across the edges of the blocks. +than across the edges of the blocks, +which shows in the differences of the same pairs of pixels. .. jupyter-execute:: + # each image and the one after it, if it began soon after with the same + # exposure, as ratio_noise compares them + before = dict(time=slice(None, ~0)) + after = dict(time=slice(1, None)) + time_frames = frames.inputs.time.ndarray + gap = na.ScalarArray((time_frames[1:] - time_frames[:-1]).to(u.s), axes="time") + change = np.abs(frames.timedelta[after] - frames.timedelta[before]) + consecutive = (gap < 25 * u.s) & (change < 0.01 * u.s) + + intensity_frames = frames.outputs + variance_frames = np.square(frames.uncertainty(noise_photon)) + z = intensity_frames[after] - intensity_frames[before] + z = z / np.sqrt(variance_frames[after] + variance_frames[before]) + z = z.to(u.dimensionless_unscaled).value + + mean = (intensity_frames[after] + intensity_frames[before]) / 2 + clipped = ( + consecutive + & ~in_box + & (mean < 10 * u.DN / u.s) + & (np.abs(z) < 5 * np.sqrt(ratio_noise)) + ) + left = dict(detector_x=slice(None, ~0)) right = dict(detector_x=slice(1, None)) diff --git a/hinode/xrt/_filtergrams.py b/hinode/xrt/_filtergrams.py index a66bf10..da660d3 100644 --- a/hinode/xrt/_filtergrams.py +++ b/hinode/xrt/_filtergrams.py @@ -32,6 +32,13 @@ and the exposure time it normalized the image to. """ +_median_chi2 = 0.454936423119572 +""" +The median of the square of Gaussian noise of unit variance, +the chi-squared distribution with one degree of freedom, +``scipy.stats.chi2.median(1)``. +""" + _pattern_saturated = re.compile( r"\(XRT_SATURATED_PIXELS\) Replaced \d+ saturated pixels with value = ([\d.]+)" ) @@ -916,6 +923,13 @@ def uncertainty( axes = (images.axis_detector_x, images.axis_detector_y) (uncertainty / images.outputs).median(axes) """ + return typing.cast(na.ScalarArray, np.sqrt(self._variance(noise_photon))) + + def _variance( + self, + noise_photon: u.Quantity | na.AbstractScalar, + ) -> na.ScalarArray: + """The square of :meth:`uncertainty`.""" if self.vignetting is None or self.uncertainty_map is None: raise ValueError( "The vignetting and the uncertainty map were not loaded. " @@ -930,7 +944,7 @@ def uncertainty( variance = signal * factor / (self.vignetting * self.timedelta) variance = variance + np.square(self.uncertainty_map) - return typing.cast(na.ScalarArray, np.sqrt(variance)) + return typing.cast(na.ScalarArray, variance) def ratio_noise( self, @@ -956,13 +970,36 @@ def ratio_noise( z = \frac{I_2 - I_1}{\sqrt{\sigma_1^2 + \sigma_2^2}}, gives a quantity whose variance is the ratio. - Only images taken within `gap` of each other with the same exposure - time are compared, - and only the pixels selected by `where` which are fainter than - `threshold` in both, + + An image is compared with the next only if the next begins after it + and within `gap` of it, + with the same exposure time, + and with the same plate scale and pointing, + so that the first pixel of each lands within a pixel of the same place + on the Sun, + since otherwise the same pixel of the two images sees different parts + of the Sun. + Of those pixels, only the ones selected by `where` in both images are + used, + which are not saturated in either, + and whose mean signal over the two images is fainter than `threshold`, since brighter ones are more likely to change. - Cosmic rays are left out by clipping :math:`z` at `clip` times the - square root of the ratio, which is repeated `num_iteration` times. + The mean is independent of the difference for noise of the same + variance in both images, + so selecting on it barely biases the ratio, + while selecting on the signal of each image would leave out the + largest differences near `threshold`. + What bias is left comes from :meth:`uncertainty` following the noise + of the signal it is found from, + and is a few percent for pixels at `threshold`. + + The ratio is first estimated from the median of :math:`z^2`, + which the outliers barely move, + as that median over the median of the square of Gaussian noise. + Cosmic rays and changes of the Sun are then left out by clipping + :math:`z` at `clip` times the square root of the ratio, + and the ratio is found again as the mean of :math:`z^2` left, + which is repeated `num_iteration` times. The noise of the images of XRT is less than :meth:`uncertainty` predicts, @@ -971,6 +1008,12 @@ def ratio_noise( so this ratio can scale :meth:`uncertainty` to the noise measured. It is measured in faint pixels, and compression removes less of the noise of brighter ones. + The errors of the dark subtraction and of the vignetting correction in + :attr:`uncertainty_map` are the same in both images and cancel in + their difference, + so the ratio is of the noise which changes from one image to the next, + and scaling all of :meth:`uncertainty` by it scales those errors too, + which in faint pixels are a small part of it. Parameters ---------- @@ -979,8 +1022,10 @@ def ratio_noise( where The pixels to measure the noise in, which should leave out any which change between images, such as those of a jet. + It may differ from one image to the next. threshold - The signal above which pixels are left out. + The mean signal of two images above which their pixels are left + out. gap The longest time between the starts of two images compared. clip @@ -989,16 +1034,32 @@ def ratio_noise( num_iteration How many times to clip the differences. + Returns + ------- + ratio + The ratio, one for each entry of any axes other than + :attr:`axis_time`, :attr:`axis_detector_x`, and + :attr:`axis_detector_y`, such as those of `noise_photon`. + + Raises + ------ + ValueError + If there are fewer than two images along :attr:`axis_time`, + or if no pair of images compared has a pixel to measure the noise + in, for some entry of the other axes. + Examples -------- The ratio for the Al_poly images captured while ESIS was observing - the Sun, with the photon noise of plasma at 2 MK. + the Sun, with the photon noise of plasma at 2 MK, + away from the jet called event E. .. jupyter-execute:: import numpy as np import astropy.units as u + import named_arrays as na import hinode images = hinode.xrt.open( @@ -1019,38 +1080,90 @@ def ratio_noise( ) noise_photon = (variance.outputs / response.outputs)[index] - images.ratio_noise(noise_photon) + # Every pixel except those of the jet + where = na.ScalarArray.ones(images.outputs[dict(time=0)].shape) > 0 + where[dict(detector_x=slice(172, 300), detector_y=slice(80, 208))] = False + + images.ratio_noise(noise_photon, where=where) """ axis = self.axis_time + axes = (axis, self.axis_detector_x, self.axis_detector_y) + + if self.outputs.shape.get(axis, 1) < 2: + raise ValueError( + "The noise is measured from pairs of images, " + f"but there are fewer than two along `axis_time`, {axis!r}." + ) + before = {axis: slice(None, ~0)} after = {axis: slice(1, None)} - intensity = self.outputs - uncertainty = self.uncertainty(noise_photon) - + # The start of each image, in seconds after the first one. time = typing.cast(na.ScalarArray, self.inputs.time) t = typing.cast(astropy.time.Time, time.ndarray) - mjd = typing.cast(np.ndarray, t.mjd) - interval = na.ScalarArray(np.diff(mjd) * u.day, axes=time.axes) - change = np.abs(self.timedelta[after] - self.timedelta[before]) - consecutive = (interval < gap) & (change < 0.01 * u.s) - - quiet = ( - consecutive - & where - & (intensity[before] < threshold) - & (intensity[after] < threshold) + seconds = na.ScalarArray((t - t.ravel()[0]).to(u.s), axes=time.axes) + interval = seconds[after] - seconds[before] + + timedelta = typing.cast(na.ScalarArray, na.as_named_array(self.timedelta)) + change = np.abs(timedelta[after] - timedelta[before]) + + # The place on the Sun of the first pixel of each image, without the + # roll, which changes little from one image to the next. + crval = typing.cast(na.Cartesian2dVectorArray, self.inputs.crval.position) + cdelt = typing.cast(na.Cartesian2dVectorArray, self.inputs.cdelt.position) + crpix = self.inputs.crpix.components + crpix_x = typing.cast(na.ScalarArray, crpix[self.axis_detector_x]) + crpix_y = typing.cast(na.ScalarArray, crpix[self.axis_detector_y]) + first_x = crval.x - cdelt.x * crpix_x + first_y = crval.y - cdelt.y * crpix_y + aligned = ( + (np.abs(cdelt.x[after] / cdelt.x[before] - 1) < 1e-3) + & (np.abs(cdelt.y[after] / cdelt.y[before] - 1) < 1e-3) + & (np.abs(first_x[after] - first_x[before]) < np.abs(cdelt.x[before])) + & (np.abs(first_y[after] - first_y[before]) < np.abs(cdelt.y[before])) ) - difference = intensity[after] - intensity[before] - z = difference / np.sqrt( - np.square(uncertainty[after]) + np.square(uncertainty[before]) + consecutive = ( + (interval > u.Quantity(0, u.s)) + & (interval < gap) + & (change < 0.01 * u.s) + & aligned ) - z = z.to(u.dimensionless_unscaled).value - ratio = np.square(z).mean(where=quiet) + # Indexing ignores the time axis of a mask without one. + selected = typing.cast(na.ScalarArray, na.as_named_array(where)) + selected = selected[before] & selected[after] + if self.saturated is not None: + selected = selected & ~self.saturated[before] & ~self.saturated[after] + + intensity = self.outputs + mean = (intensity[before] + intensity[after]) / 2 + + variance = self._variance(noise_photon) + z2 = np.square(intensity[after] - intensity[before]) + z2 = z2 / (variance[after] + variance[before]) + z2 = z2.to(u.dimensionless_unscaled).value + + quiet = consecutive & selected & (mean < threshold) & np.isfinite(z2) + + if np.any(quiet.sum(axes) == 0): + raise ValueError( + "No pair of images compared has a pixel to measure the noise in. " + "Each image is compared only with the next one, " + "if it begins within `gap` after it, with the same exposure, " + "plate scale, and pointing, " + "and only in the pixels selected by `where`, not saturated, " + "and whose mean over the two is fainter than `threshold`." + ) + + # NaN where a pixel is not quiet, which the median leaves out + masked = np.where(quiet, z2, np.nan) + # `np.nanmedian` takes the names of the axes of a named array, + # which the annotations of numpy do not allow. + median = np.nanmedian(masked, axis=axes) # type: ignore[call-overload] + ratio = typing.cast(na.ScalarArray, median) / _median_chi2 for _ in range(num_iteration): - clipped = quiet & (np.abs(z) < clip * np.sqrt(ratio)) - ratio = np.square(z).mean(where=clipped) + clipped = quiet & (z2 <= np.square(clip) * ratio) + ratio = z2.mean(axes, where=clipped) return typing.cast(na.ScalarArray, ratio) diff --git a/hinode/xrt/_filtergrams_test.py b/hinode/xrt/_filtergrams_test.py index 0e1e9ac..d6fc060 100644 --- a/hinode/xrt/_filtergrams_test.py +++ b/hinode/xrt/_filtergrams_test.py @@ -8,8 +8,10 @@ import astropy.time import astropy.wcs import astropy.io.fits +import scipy.stats import named_arrays as na import hinode +from hinode.xrt._filtergrams import _median_chi2 from hinode.xrt._vignetting_test import _vignetting_numpy from hinode.xrt._coalign_test import _corrected @@ -787,6 +789,7 @@ def _sequence( ) -> tuple[hinode.xrt.Filtergram, u.Quantity]: """ A sequence of `num` flat images of the given `signal`, in DN per second, + taken 18 seconds apart, with noise whose variance is `ratio` times that of :meth:`hinode.xrt.Filtergram.uncertainty`, and the photon noise factor the uncertainty was computed with. @@ -794,15 +797,33 @@ def _sequence( path = _path("2019-09-30T18:08:30", "2019-09-30T18:08:40") paths = na.ScalarArray(np.array([path] * num, dtype=object), axes="time") images = hinode.xrt.Filtergram.from_fits(paths, uncertainty=True) + time = images.inputs.time.ndarray + 18 * np.arange(num) * u.s + inputs = images.inputs.replace(time=na.ScalarArray(time, axes="time")) noise_photon = 2.5 * u.DN / u.ph flat = na.ScalarArray.ones(images.outputs.shape) * signal * u.DN / u.s - images = dataclasses.replace(images, outputs=flat) + images = dataclasses.replace(images, inputs=inputs, outputs=flat) sigma = images.uncertainty(noise_photon) noise = na.random.normal(0, 1, shape_random=flat.shape, seed=seed) outputs = flat + np.sqrt(ratio) * sigma * noise return dataclasses.replace(images, outputs=outputs), noise_photon +def _replace_images( + images: hinode.xrt.Filtergram, + noisy: hinode.xrt.Filtergram, + index: tuple[int, ...], +) -> hinode.xrt.Filtergram: + """`images`, with the images at each of `index` replaced by those of `noisy`.""" + outputs = images.outputs.copy() + for i in index: + outputs[dict(time=i)] = noisy.outputs[dict(time=i)] + return dataclasses.replace(images, outputs=outputs) + + +def test_median_chi2() -> None: + assert np.isclose(_median_chi2, scipy.stats.chi2.median(1), rtol=1e-12) + + def test_ratio_noise() -> None: """The ratio is the variance of the noise of the images over their uncertainty.""" images, noise_photon = _sequence(ratio=0.2) @@ -819,18 +840,14 @@ def test_ratio_noise_pairs() -> None: images, noise_photon = _sequence(ratio=0.2, num=5) noisy, _ = _sequence(ratio=50, num=5, seed=1) - time = images.inputs.time.ndarray + [0, 18, 36, 54, 100] * u.s + time = images.inputs.time.ndarray[0] + [0, 18, 36, 54, 100] * u.s timedelta = na.ScalarArray([16.4, 16.4, 11.6, 16.4, 16.4] * u.s, axes="time") # Only the first two images are compared: the third has another # exposure than its neighbors, and the last begins long after the fourth. - outputs = images.outputs.copy() - for i in (2, 4): - outputs[dict(time=i)] = noisy.outputs[dict(time=i)] - + images = _replace_images(images, noisy, index=(2, 4)) images = dataclasses.replace( images, inputs=images.inputs.replace(time=na.ScalarArray(time, axes="time")), - outputs=outputs, timedelta=timedelta, ) @@ -841,13 +858,65 @@ def test_ratio_noise_pairs() -> None: assert images.ratio_noise(noise_photon, gap=60 * u.s).ndarray > 1 +def test_ratio_noise_order() -> None: + """ + An image is not compared with the next if the next began before it, + or at the same time, as a file loaded twice does. + """ + images, noise_photon = _sequence(ratio=0.2, num=6) + noisy, _ = _sequence(ratio=50, num=6, seed=1) + + # The fourth image began before the third, and the sixth is the fifth + # again, so only the first two are compared. + time = images.inputs.time.ndarray[0] + [0, 18, 54, 36, 72, 72] * u.s + images = _replace_images(images, noisy, index=(3,)) + outputs = images.outputs.copy() + outputs[dict(time=5)] = outputs[dict(time=4)] + images = dataclasses.replace( + images, + inputs=images.inputs.replace(time=na.ScalarArray(time, axes="time")), + outputs=outputs, + ) + + result = images.ratio_noise(noise_photon) + assert np.isclose(result.ndarray, 0.2, rtol=0.05) + + +def test_ratio_noise_aligned() -> None: + """ + Images are compared only if they have the same plate scale and point + within a pixel of each other. + """ + images, noise_photon = _sequence(ratio=0.2, num=5) + noisy, _ = _sequence(ratio=50, num=5, seed=1) + images = _replace_images(images, noisy, index=(2, 4)) + + # The second image points half a pixel from the first, and is compared + # with it, the third two pixels away, and the fifth is binned 2 by 2. + crval = images.inputs.crval.position + cdelt = images.inputs.cdelt.position + x = crval.x.copy() + x[dict(time=1)] = x[dict(time=1)] + 0.5 * cdelt.x[dict(time=1)] + x[dict(time=2)] = x[dict(time=2)] + 2 * cdelt.x[dict(time=2)] + scale = cdelt.y.copy() + scale[dict(time=4)] = 2 * scale[dict(time=4)] + inputs = images.inputs.replace( + crval=na.PositionalVectorArray(crval.replace(x=x)), + cdelt=na.PositionalVectorArray(cdelt.replace(y=scale)), + ) + images = dataclasses.replace(images, inputs=inputs) + + result = images.ratio_noise(noise_photon) + assert np.isclose(result.ndarray, 0.2, rtol=0.05) + + def test_ratio_noise_excluded() -> None: """ Cosmic rays are clipped, and pixels which are bright or not selected are left out. """ images, noise_photon = _sequence(ratio=0.2) - bright, _ = _sequence(ratio=2, signal=12, seed=1) + bright, _ = _sequence(ratio=2, signal=20, seed=1) outputs = images.outputs.copy() # cosmic rays @@ -869,4 +938,113 @@ def test_ratio_noise_excluded() -> None: assert images.ratio_noise(noise_photon).ndarray > 0.22 threshold = 100 * u.DN / u.s assert images.ratio_noise(noise_photon, where, threshold).ndarray > 0.22 - assert images.ratio_noise(noise_photon, where, num_iteration=0).ndarray > 0.22 + + # the first estimate, from the median, is barely moved by the cosmic rays + result = images.ratio_noise(noise_photon, where, num_iteration=0) + assert np.isclose(result.ndarray, 0.2, rtol=0.03) + + +def test_ratio_noise_outliers() -> None: + """ + Outliers in many pixels, which would raise a first estimate from the + mean enough to survive the clipping, are clipped. + """ + images, noise_photon = _sequence(ratio=0.2) + # a change in 5% of the pixels of one image, of about nine times the noise + region = dict(time=1, detector_x=slice(0, 20)) + outputs = images.outputs.copy() + outputs[region] = outputs[region] + 5 * u.DN / u.s + images = dataclasses.replace(images, outputs=outputs) + + result = images.ratio_noise(noise_photon) + assert np.isclose(result.ndarray, 0.2, rtol=0.02) + + +def test_ratio_noise_threshold() -> None: + """ + Selecting the pixels by their mean over two images barely biases the + ratio of those near the threshold, where selecting them by the signal of + each image would leave out the largest differences and give about 0.08. + """ + images, noise_photon = _sequence(ratio=0.2, signal=10) + result = images.ratio_noise(noise_photon) + # The uncertainty follows the noise of the signal it is found from, + # so the pixels selected, whose mean is fainter, have a little less. + assert np.isclose(result.ndarray, 0.2, rtol=0.05) + + +def test_ratio_noise_where_time() -> None: + """The pixels selected may differ from one image to the next.""" + images, noise_photon = _sequence(ratio=0.2) + noisy, _ = _sequence(ratio=1, seed=1) + region = dict(time=2, detector_x=slice(0, 100)) + outputs = images.outputs.copy() + outputs[region] = noisy.outputs[region] + images = dataclasses.replace(images, outputs=outputs) + + where = na.ScalarArray.ones(images.outputs.shape) > 0 + where[region] = False + + result = images.ratio_noise(noise_photon, where=where) + assert np.isclose(result.ndarray, 0.2, rtol=0.02) + assert images.ratio_noise(noise_photon).ndarray > 0.22 + + +def test_ratio_noise_saturated() -> None: + """ + Saturated pixels, which are the same in every image, are left out even + when they are fainter than the threshold. + """ + images, noise_photon = _sequence(ratio=0.2) + assert images.saturated is not None + region = dict(detector_x=slice(0, 50)) + outputs = images.outputs.copy() + outputs[region] = 150 * u.DN / u.s + saturated = images.saturated.copy() + saturated[region] = True + images = dataclasses.replace(images, outputs=outputs, saturated=saturated) + + threshold = 1000 * u.DN / u.s + result = images.ratio_noise(noise_photon, threshold=threshold) + assert np.isclose(result.ndarray, 0.2, rtol=0.02) + + +def test_ratio_noise_axes() -> None: + """There is a ratio for each entry of the other axes, such as those of `noise_photon`.""" + images, noise_photon = _sequence(ratio=0.2) + factor = na.ScalarArray(np.array([1, 2]), axes="factor") + result = images.ratio_noise(factor * noise_photon) + assert result.shape == dict(factor=2) + assert np.isclose(result[dict(factor=0)].ndarray, 0.2, rtol=0.02) + assert result[dict(factor=1)] < result[dict(factor=0)] + + +def test_ratio_noise_timedelta() -> None: + """The exposure time may be one for every image.""" + images, noise_photon = _sequence(ratio=0.2) + timedelta = na.as_named_array(images.timedelta)[dict(time=0)].ndarray + images = dataclasses.replace(images, timedelta=timedelta) + result = images.ratio_noise(noise_photon) + assert np.isclose(result.ndarray, 0.2, rtol=0.02) + + +def test_ratio_noise_nan() -> None: + """Pixels whose uncertainty is not finite are left out.""" + images, noise_photon = _sequence(ratio=0.2) + factor = na.ScalarArray.ones(images.outputs[dict(time=0)].shape) + factor[dict(detector_x=5, detector_y=5)] = np.nan + result = images.ratio_noise(factor * noise_photon) + assert np.isclose(result.ndarray, 0.2, rtol=0.02) + + +def test_ratio_noise_none() -> None: + """Without a pair of images or a pixel to compare, there is no ratio.""" + images, noise_photon = _sequence(ratio=0.2) + with pytest.raises(ValueError, match="fewer than two"): + images[dict(time=slice(0, 1))].ratio_noise(noise_photon) + with pytest.raises(ValueError, match="fewer than two"): + images[dict(time=0)].ratio_noise(noise_photon) + with pytest.raises(ValueError, match="No pair"): + images.ratio_noise(noise_photon, gap=10 * u.s) + with pytest.raises(ValueError, match="No pair"): + images.ratio_noise(noise_photon, where=False)