Change point detection and image segmentation for time series of astrophysical images

We require that the data are binned into photon counts in an $N_{\rm I}\times N_{\rm T}\times N_{\rm W}$ tensor $\{y_{i,t,w}\},i=1,...,N_{\rm I},t=1,...,N_{\rm T},w=1,...,N_{\rm W}$, where $y_{i,t,w}$ is the photon counts within the $i^{\rm th}$ spatial rectangular region, the $t^{\rm th}$ time interval $[T_{t-1},T_{t})$ and the $w^{\rm th}$ energy range $[W_{w-1},W_{w})$. After binning, the data can be viewed as a time series of images in different energy bands. The values of each pixel are the photon counts in the different bands in the corresponding spatial region.

Notice that compared with Automark (Wong et al., 2016), we incorporate 2-D spatial information into the model, thus extending the analysis from two (wavelength/energy and time) to four dimensions (wavelength/energy, time, and projected sky location). We also relax the restriction that the bin sizes along any of the axes are held fixed. Thus, sharp changes are more easily detected.

As the data in high-energy astrophysics are photon counts, we use a Poisson process to model the data,

where $\Delta T_{t}=(T_{t}-T_{t-1})$.

Our goal is to infer model intensities $\lambda_{i,t,w}$ from the observed counts data $\{y_{i,t,w}\}$. We are especially interested in detecting significant changes of $\lambda_{i,t,w}$ over time. If there are changes, we also want to estimate the number and locations of the change points.

To simplify the presentation, we first develop a time-homogeneous model, i.e., one where there are no change points and $\lambda_{i,t,w}$ is unchanging with $t$ (Section 2.2.2). We will then consider more complex cases, where change points are added to the model so that $\lambda_{i,t,w}$ is allowed to change over time (Section 2.2.3)

First consider a temporally homogeneous Poisson model without any change points. Then each image can be treated as an independent Poisson realization of the same, unknown, true image.

We model the image as a 3-dimensional piecewise constant function. That is, the 2-dimensional space of $x$-$y$ coordinates is partitioned into $m$ non-overlapping regions such that all the pixels in a given region have the same Poisson intensity. Different energy bands share the same spatial partitioning. Rigorously, the Poisson parameter $\lambda_{i,t,w}$ can be written as a summation of region-specific Poisson rates $\mu_{h,w}$ times the corresponding indicator functions of regions ($I_{\{i\in R_{h}\}}$ is 1 for pixel $i$ in region $R_{h}$ and 0 otherwise) in the following format:

Here $i\in R_{h}$ means “the $i^{\rm th}$ pixel in the $h^{\rm th}$ region” and $I$ is the indicator function. $R_{h}$ is the index set of the pixels within the $h^{\rm th}$ region, with $R_{h}\subseteq\{1,...,N_{\rm I}\}$. Also, $\mu_{h,w}$ is the Poisson rate for the $w^{\rm th}$ band of the $h^{\rm th}$ region. The partition of the image is specified by $\bm{R}=\{R_{h}|h=1,...,m\}$.

Now we allow the underlying Poisson parameter $\lambda_{i,t,w}$ to change over time $t$. We model $\lambda_{i,t,w}$ as a piecewise constant function of $t$.

Suppose these $N_{\rm T}$ images can be partitioned into $K+1$ homogeneous intervals by $K$ change points

For the $t^{\rm th}$ image, suppose that it belongs to the $k^{\rm th}$ time interval; i.e., $t\in(\tau_{k-1},\tau_{k}]$. For each given $t$, let $\lambda$ be a two-dimensional piecewise constant function with $m^{(k)}$ constant regions. Then $\lambda$ can be represented by

where $m^{(k)}$ is the number of regions within the $k^{\rm th}$ interval. Let $\mathcal{M}=\{m^{(k)}|k=1,2,...,K+1\}$. The partition of the images within interval $k$ is specified by $\bm{R}^{(k)}=\{R_{h}^{(k)}|h=1,...,m^{(k)}\}$. And the overall partition is $\mathcal{R}=\{\bm{R}^{(k)}|k=1,2,...,K+1\}$. The Poisson rates $\mu_{h,w}^{(k)}$ is the value for the $w^{\rm th}$ band in the $h^{\rm th}$ region of the $k^{\rm th}$ interval. Let $\bm{\mu}^{(k)}=\{\mu_{h,w}^{(k)}|h=1,...,m^{(k)},w=1,...,N_{\rm W}\}$. And let $\bm{\mu}=\{\bm{\mu}^{(k)}|k=1,...,K+1\}$, and $i\in R_{h}^{(k)}$ means “the $i^{\rm th}$ pixel is in the $h^{\rm th}$ region of the $k^{\rm th}$ interval”.

Given the observed images $\{y_{i,t,w}\}$, we aim to obtain an estimate of $\lambda_{i,t,w}$. In other words, we want an estimate of the image partitions and the Poisson rates of the regions for each band. It is straightforward to estimate the Poisson intensities given the region partitioning, but the partitioning is a much more complicated model selection problem.

We will apply MDL to select the best fitting model. Loosely speaking, the idea behind MDL for model selection is to first obtain a MDL criterion for each possible model, and then define the best fitting model as the minimizer of this criterion. MDL defines the best model as the one that produces the best compression of the data. The criterion can be treated as the code length, or amount of hardware memory required to store the data.

First we present the MDL criterion for the homogeneous Poisson model, then follow it by the MDL criterion for the general case (i.e., with change points).

Following similar arguments as in Lee (2000) (see their Appendix B), the MDL criterion for segmenting $N_{\rm T}$ homogeneous images is

where $a_{h}$ and $b_{h}$ are, respectively, the “area” (number of pixels) and “perimeter” (number of pixel edges) of region $R_{h}$, and

is the maximum likelihood estimate of the Poisson rate in the corresponding region. Note that the indices of $\hat{\bm{\mu}}=\{\hat{\mu}_{hw}\}$ run over the region segments $h=1..m$ and the passbands $w=1..N_{\rm W}$.

For the Poisson model with change points, once the number of change points $K$ and the locations $\bm{\tau}=\{\tau_{1},...,\tau_{K}\}$ are specified, for each $k\in(1,2,...,K+1)$, $m^{(k)}$ and $\bm{R}^{(k)}$ can be estimated independently. Using the previous argument, the MDL criterion for images within the same homogeneous interval is

Then the overall MDL criterion for the model with change points is

To sum up, using the MDL principle, the best-fit model is defined as the minimizer of the criterion (7). The next subsection presents a practical algorithm for carrying out this minimization.

An important step to demonstrating the efficacy of our method is to establish its statistical consistency. That is, if it is shown that as the size of the data increases, the differences between the estimated model parameters and the true values decrease to zero, then the method can be said to be free of asymptotic bias, can be applied in the general case, and is elevated above a heuristic. We prove in Appendix A that the MDL-based model selection to choose the region partitioning, as well as the corresponding Poisson intensity parameters, is indeed strongly statistically consistent under mild assumptions of maintaining the temporal variability structure of $\lambda_{i,t,w}$.

Given its complicated structure, global minimization of $\text{MDL}_{\text{overall}}(K,\tau,\mathcal{M},\mathcal{R},\hat{\mu})$ (Equation 7) is virtually infeasible when the number of images $N_{T}$ and the number of pixels $N_{I}$ are not small, because the time complexity of the exhaustive search is of order $2^{N_{\rm T}N_{\rm I}}$.

We iterate the following two steps to (approximately) minimize the MDL criterion (7).

1. 1.

   Given a set of change points, apply the image segmentation method to all the images belonging to the first homogeneous time interval and obtain the MDL best-fitting image for this interval. Repeat this for all remaining intervals. Calculate the MDL criterion (7).

2. 2.

   Modify the set of change points by, for example, adding or removing one change point. In terms of what modification should be made, we use the greedy strategy to select the one that achieves the largest reduction of the overall MDL value in (7).

For Step 1 we begin with a large set of change points (i.e., an over-fitted model). In Step 2 we remove one change point (i.e., merge two neighboring intervals) to maximize the reduction of the MDL value. The procedure stops and declares the optimization is done if no further MDL reduction can be achieved by removing change points. This is similar to backward elimination for statistical model selection problems. See Figure 1 for a flowchart of the whole procedure.

Here we list some practical issues that are crucial to the success of the above minimization algorithm.

Initial seed allocation for SRG: The selection of the initial seeds in SRG plays an important role in the performance of the algorithm. To obtain a good initial oversegmentation, there must be at least one seed within each true region. Currently we use all the local maxima as well as a subset of the square lattice as the initial seeds. Based on simulation results (see Section 3.2), when the number of initial seeds is inadequate, the SRG will underfit the images, which will in turn lead to an overfitting of the change points and will lead to an increased false positive rate. On the other hand, it could be time-consuming when the number of initial seeds is very large, especially for high resolution images. We developed an algorithm that allocates initial seeds automatically based on locating local maxima. However, an optimal selection of initial seeds almost certainly requires expert intervention, as it depends on the type of data that we work with. To reduce the chance of obtaining a poor oversegmentation with SRG, in practice one could try using different sets of initial seeds and select the oversegmentation that gives the smallest MDL value.

Counts per bin: The photon counts in the image pixels cannot be too small, otherwise the algorithm could fail to produce meaningful output; see Section 3.2. In some sense, small photon counts can be seen as low signal level, which means that the proposed method requires a minimum level of signal to operate with. Therefore, care must be exercised when deciding the size for the bin. As a rule of thumb, it should be enough to have around 100 counts for each pixel belonging to an astronomical object, while pixels from the background can have very low or even zero count.

Initial change point selection: Although the stepwise greedy algorithm is capable of saving a significant amount of computation time, it could still be time consuming if the initial set of change points is too large, as it might need many iterations to reach a local minimum. It is recommended to select the initial change points based on prior knowledge, if available, in order to accelerate the algorithm.

Computation Time: In each iteration, the main time-consuming part is to apply the SRG. When the total number of pixels and the number of seeds are large, during the process, the number of candidates is large. Therefore, the comparison among all the candidates and the following updating manner lead to most of the computation burden. As an example, we found that it takes about 40 minutes for applying SRG and merging once on $64\times 64$ images with about 200 seeds on a Linux machine with an octa-core 2.90 GHz Intel Xeon processor.

The first method is to highlight key pixels based on the distribution of pixel differences before and after change points. The rationale is that a pixel with different fitted values before and after a change point is a strong indicator that it is a key pixel.

Suppose the fitted values for pixel $i$ in time intervals $k$ and $(k+1)$ are $\hat{\lambda}_{i}^{(k)}$ and $\hat{\lambda}_{i}^{(k+1)}$, respectively. Given the Poisson nature of the data, we first apply a square-root transformation to normalize the fitted values. Define the difference $d_{i}$ for pixel $i$ as

A pixel is labelled as a key pixel if its $d_{i}$ is far away from the mean of all the differences. To be specific, pixel $i$ is labelled as a key pixel if

where $\hat{\mu}=\frac{1}{N_{\rm I}}\sum_{1}^{N_{\rm I}}d_{i}$ and $\hat{\sigma}=\frac{1}{\Phi^{-1}(3/4)}\text{MAD}$. Here $\text{MAD}=\text{median}(|d_{i}-\tilde{d}|)$ is the median absolute deviation with $\tilde{d}=\text{median}(d_{i})$, and is used to obtain a robust estimate (i.e., a measure that minimizes the effect of outliers) of the standard deviation of the $d_{i}$’s (see e.g., Rousseeuw & Croux, 1993). $\Phi^{-1}(\cdot)$ is the quantile of the standard normal distribution, and $p$ is the pre-specified significance level¹¹1$\Phi^{-1}$ is related to the standard Normal error function, with exemplar values $\Phi^{-1}(\{0.75,0.841,0.977,0.9986,1-\frac{10^{-5}}{2},1-\frac{10^{-10}}{2},1-\frac{10^{-15}}{2}\})=\{0.6745,1,2,3,4.417,6.467,8.014\}$. We typically choose $p=1-\frac{10^{-15}}{2}$ as our threshold.. Notice that by checking the sign of $d_{i}-\hat{\mu}$, we can deduce if pixel $i$ has increased or decreased after this change point.

Another method to locate key pixels is to compare pairs of regions. For any region in the time interval $k$, there must exist at least one region in time interval $(k+1)$ such that these two regions have overlapping pixels. We then test if the difference between the means of the pixels from these two regions is significant or not.

As before, we apply the square-root transformation to the pixels within each of the regions. Then we calculate the sample means $\hat{\mu}_{1}$ and $\hat{\mu}_{2}$ and sample variances $\hat{\sigma}_{1}^{2}$ and $\hat{\sigma}_{2}^{2}$ of these two groups of square-rooted values. Then we can for example test whether the difference between $\hat{\mu}_{1}$ and $\hat{\mu}_{2}$ is large enough with

See Section 4 for the applications of these two methods on some real data sets.

Lastly we note that the selection of $p$, the significance level, deserves a more careful consideration. As in reality one may need to do comparisons for many change points and energy bands, this becomes a multiple-testing problem where the number of tests is large. Therefore, one should adjust the value of $p$ in order to control false positives.

As there is no change point in $\lambda_{i,t,w}$, this experiment is ideal for studying the relationship between the false positive rate and the signal level; recall the latter is defined as the average number of photon counts per pixel.

The results of this experiment (together with the next two experiments) are summarized as the blue curves in Figure 3. The figure captures how well the simulation recovers the location of the change points (top left), the number of change points (top right), the excess number of change points (false positives; bottom left), and the deficit in change points (false negatives; bottom right). The top left plot reports the fraction of simulated datasets for which the set of the fitted change points $\hat{\bm{\tau}}$ is identical to the set of true change points $\bm{\tau}$. The top right plot presents the fraction of simulated datasets for which the fitted number of change points $\hat{K}$ equals to the true number of change points $K$. The bottom left plot shows the average false positive rate, which is defined as the average number of falsely detected change points per possible location. The bottom right plot presents the fraction of simulated datasets for which $\hat{\bm{\tau}}$ contains $\bm{\tau}$, i.e., $\bm{\tau}\subseteq\hat{\bm{\tau}}$. One can see that the false discovery rate seems to be quite stable across different signal levels. Notice that the last curve is always 1 because $\bm{\tau}$ is empty for this experiment.

In this experiment we introduced variation in $\lambda_{i,t,w}$ by multiplying (a) and (b) in Figure 2 together. The results are reported as the green curve in Figure 3. One can see that when the signal level is small ($\log_{10}(\text{average signal level})<1.0$), increasing the signal level leads to more false positives: this range of signals levels are too small to provide enough information for the detection of the true change points.

When $\log_{10}(\text{average signal level})$ is between $1.0$ and $2.0$, the proposed method starts to be able to detect the true change points as the signal level increases. Also, the false positive rate begins to drop. One possible explanation is that for any two consecutive homogeneous time intervals, the location of the change point might not be clear when the signal level is relatively small.

As the signal level continues to increase, the proposed method becomes more successful in detecting the true change point locations. For example, when the average number of counts in each bin is greater than 100 (i.e., $\log_{10}(\text{average signal level})>2.0$), the signal is strong enough so that all the true change points can be detected successfully. We note that there are always some false positives due to the Poisson randomness, and it seems that the false positive rate stabilizes as the signal level increases.

Here we allow different bands $w$ to change differently at the change points. The rate $\lambda_{i,t,w}$ was obtained by multiplying (a) and (c) in Figure 2 together. The results, which are about the same as the previous experiment (varying intensity), are reported as the red curve in Figure 3. When the signal level is small ($\log_{10}(\text{average signal level})<1.0$), the proposed method fails to detect the true change points. When $\log_{10}(\text{average signal level})$ is between $1.0$ and $2.0$, as the signal level increases, the false positive rate begins to decrease while the true positive rate increases. When $\log_{10}(\text{average signal level})>2.0$, all the true change points can be detected successfully, while the false positive rate stays at the same level as signal level increases.

There was no change point in this experiment and the spatial variation of $\lambda_{i,t,w}$ is given in the bottom right plot of Figure 2. The results are reported as the blue curves in Figure 4. One can see that if the number of initial seeds is inadequate, the false positive rate increases as the signal level increases above $\log_{10}(\text{average signal level})>2.5$. However, this does not happen when there are a large number of initial seeds; see the blue dotted curves in Figure 4. In fact, for this and the following two experiments, our method did not detect any false positive change points. This suggests that when the images are under-segmented, the method tends to place more false change points to compensate for data variability not explainable by image segmentation.

In this experiment $\lambda_{i,t,w}$ was obtained by multiplying (b) and (d) of Figure 2 together, so there are three change points over time. The results are similar to the no change point case, and revsummarized as the green curves in Figure 4.

In this last experiment the energy bands were allowed to be different, and $\lambda_{i,t,w}$ was obtained by multiplying (c) and (d) of Figure 2 together. The results, reported as the red curves in Figure 4, are similar to the previous two experiments.

We use a dataset from XMM-Newton (Obs.ID 0049350101), where Proxima Centauri was observed for 67 ks on 2001-08-12. Because of the high flux, the MOS cameras on XMM-Newton are highly piled-up and we restrict our analysis to the data from the PN camera. We obtained the data from the XMM-Newton science archive hosted by the European Space Agency (ESA)²²2https://www.cosmos.esa.int/web/xmm-newton/xsa. The data we received was processed by ODS version 12.0.0. Our analysis is based on the filtered PN event data from the automated reduction pipeline (PPS). Güdel et al. (2002) presented a detailed analysis and interpretation of this dataset.

In our analysis, we only used a subset of photons with spatial coordinates within $[25500,27500]\times[26500,28500]$, and it was binned as images of size $64\times 64$. We used the temporal bins of width $1100.4$ seconds to generate $60$ images. We binned the data into three energy bands, $(200,1000]$, $(1000,3000]$ and $(3000,10000]$ in eV.

Figure 5 presents the light curves for different bands as well as the locations of the detected change points. As there is only a single source of photons with negligible background signal level, the detected change points coincide with the changes of the light curves. Many change points are detected for the abrupt increase and then decrease in brightness for all the bands at the time points between 42 and 50. The time interval between 15 and 36 is detected as a homogeneous time interval, and the variation in light curves within this interval is viewed as common Poisson variations. A few change points are detected for the time interval before time point 15 and the interval after 50. A piecewise constant model is used to fit these gradual changes in intensities.

The fitted images can be found in Figure 6. The source of most photons, a point source, is modeled by different piecewise constant models at all these time intervals. The center of the point source and the wings of the PSF region nearby were fitted by models with shapes like concentric circles.

To test our method for finding the regions of significant change, we also apply it here because we know what the answer should be. For this dataset, as the observations for different bands change simultaneously, we combine all the three bands to highlight the key pixels. That is, we highlight the regions for each band and take the intersection of these regions. Examples of the results based on the method in Section 2.4.2 can be found in Figure 6. The method indeed picks out the point source. With significance level $p=10^{-2}$, the abrupt increase in brightness of the source at time point 41 and 42, as well as the sudden decrease at time point 43 can be detected successfully. By modifying the significance level, different sensitivity can be achieved.

In particular, here we consider AIA observations carried out on 2014-Dec-11 between 19:12 UT and 19:23 UT, and focus on a $64{\times}64$ pixel region located $(+1^{\prime\prime},-271^{\prime\prime})$ from disk center, in which a small, isolated, well-defined loop appeared at approximately 19:19 UT. This region was selected solely as a test case to demonstrate our method; the appearance of the loop is clear and unambiguous, with no other event occurring nearby to confuse the issue; see Figure 7. We apply our method to these data, downloaded using the SDO AIA Cutout Service,³³3https://www.lmsal.com/get_aia_data/ and demonstrate that the loop (and it alone) is detected and identified; see Figures 8, 9 and 10. AIA data are available in 6 filter bands, centered at 211, 94, 335, 193, 131, 171 Å. Here, we have limited our analysis to 3 bands: 94, 335, and 131 Å in which the isolated loop is easily discernible to the eye (a full analysis including all the filters does not change the results). Each filter consists of a sequence of 54 images, and while they are not obtained simultaneously, the difference in time between the bands is ignorable on the timescale over which the loop evolves.

The fitted images for the 3-band case can be found in Figure 8. Notice that there is a loop-shaped object that is of interest. Based on the fitted result, this object starts to appear at time point c.36 and becomes brighter after that for the first band. In the second band, this object appears at time point c.26 and stays bright throughout the duration considered. However in the third band, the object becomes bright at time point c.36 and vanishes soon after time point c.38. The proposed method is able to catch these changes in different bands and to detect the corresponding change points.

After detecting these change points, we find the key pixels that contribute to the change points using the methods in Section 2.4.1. This method is appropriate to highlight the regions that change rapidly after the change point because different bands may not change in the same direction for this dataset. Here we apply this method on a single band, 94, as an example. We use $p=10^{-15}$. See Figure 9 for an illustration. We find that the method could highlight the loop-shape object which starts to appear at time point c.36, and also detect the region that becomes much brighter after time point c.38. We also compute the light curves of the intensities in a region comprised of the set of pixels formed from the union of all key pixels found at all change points in all the filters; see Figure 10. Notice that the event of interest is fully incorporated within the key pixels, with no spillover into the background, and the change points are post facto found to be reasonably located from a temporal perspective in that they are located where a researcher seeking to manually place them would do so. The first segment is characterized by steady emission in all three bands, the second segment shows the isolated loop beginning to form, the third segment catches the time when it reaches a peak, and the last segment tracks the slow decline in intensity.