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 6092d6c..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,4 +944,226 @@ 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, + 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. + + 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. + 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, + 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. + 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 + ---------- + 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. + It may differ from one image to the next. + threshold + 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 + How many times the measured noise a difference may be before it + is left out. + 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, + 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( + 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] + + # 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)} + + # 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) + 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])) + ) + + consecutive = ( + (interval > u.Quantity(0, u.s)) + & (interval < gap) + & (change < 0.01 * u.s) + & aligned + ) + + # 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 & (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 7346dcf..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 @@ -777,3 +779,272 @@ 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, + 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. + """ + 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, 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) + 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] + [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. + images = _replace_images(images, noisy, index=(2, 4)) + images = dataclasses.replace( + images, + inputs=images.inputs.replace(time=na.ScalarArray(time, axes="time")), + 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_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=20, 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 + + # 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)