  Many astrophysical phenomena are time-varying, in the sense that their
intensity, energy spectrum, and/or the spatial distribution of the emission
suddenly change. This paper develops a method for modeling a time series of
images. Under the assumption that the arrival times of the photons follow a
Poisson process, the data are binned into 4D grids of voxels (time, energy
band, and x-y coordinates), and viewed as a time series of non-homogeneous
Poisson images. The method assumes that at each time point, the corresponding
multi-band image stack is an unknown 3D piecewise constant function including
Poisson noise. It also assumes that all image stacks between any two adjacent
change points (in time domain) share the same unknown piecewise constant
function. The proposed method is designed to estimate the number and the
locations of all the change points (in time domain), as well as all the unknown
piecewise constant functions between any pairs of the change points. The method
applies the minimum description length (MDL) principle to perform this task. A
practical algorithm is also developed to solve the corresponding complicated
optimization problem. Simulation experiments and applications to real datasets
show that the proposed method enjoys very promising empirical properties.
Applications to two real datasets, the XMM observation of a flaring star and an
emerging solar coronal loop, illustrate the usage of the proposed method and
the scientific insight gained from it.
