Transform in Astronomical Data Processing



Fig. 1
Example of astronomical data: a point source and an extended source are shown, with noise and background. The extended object, which can be detected by eye, is undetected by a standard detection approach






2 Standard Approaches to Source Detection


We describe here the most popular way to create a catalog of galaxies from astronomical images.


The Traditional Data Model



The observed data Y can be decomposed into two parts, the signal X and the noise N:


$$\displaystyle{ Y [k,l] = X[k,l] + N[k,l] }$$

(1)
The imaging system can also be considered. If it is linear, the relation between the data and the image in the same coordinate frame is a convolution:


$$\displaystyle\begin{array}{rcl} Y [k,l]& =& (HX)[k,l] + N[k,l]{}\end{array}$$

(2)
where H is the matrix related to the point spread function (PSF) of the imaging system.

In most cases, objects of interest are superimposed on a relatively flat signal B, called background signal. The model becomes


$$\displaystyle{ Y [k,l] = (HX)[k,l] + B[k,l] + N[k,l] }$$

(3)


PSF Estimation


The PSF H can be estimated from the data or from an optical model of the imaging telescope. In astronomical images, the data may contain stars, or one can point towards a reference star in order to reconstruct a PSF. The drawback is the “degradation” of this PSF because of unavoidable noise or spurious instrument signatures in the data. So, when reconstructing a PSF from experimental data, one has to reduce very carefully the images used (background removal for instance). Another problem arises when the PSF is highly variable with time, as is the case for adaptive optics (AO) images. This means usually that the PSF estimated when observing a reference star, after or before the observation of the scientific target, has small differences from the perfectly correct PSF.

Another approach consists of constructing a synthetic PSF. Various studies [11, 21, 38, 39] have suggested a radially symmetric approximation to the PSF:


$$\displaystyle\begin{array}{rcl} P(r) \propto \left (1 + \frac{r^{2}} {R^{2}}\right )^{-\beta }& &{}\end{array}$$

(4)
The parameters β and R are obtained by fitting the model with stars contained in the data.

In the case of AO systems, this model can be used for the tail of the PSF (the so-called seeing contribution), while in the central region, the system provides an approximation of the diffraction-limited PSF. The quality of the approximation is measured by the Strehl ratio (SR), which is defined as the ratio of the observed peak intensity in the image of a point source to the theoretical peak intensity of a perfect imaging system.


Background Estimation


The background must be accurately estimated; otherwise it will introduce bias in flux estimation. In [7, 28], the image is partitioned into blocks, and the local sky level in each block is estimated from its histogram. The pixel intensity histogram p(Y ) is modeled using three parameters, the true sky level B, the RMS (root mean square) noise σ, and a parameter describing the asymmetry in p(Y ) due to the presence of objects, and is defined by [7]:


$$\displaystyle{ p(Y ) = \frac{1} {a}\exp \left (\sigma ^{2}/2a^{2}\right )\exp \left [-(Y - B)/a\right ]{\mathrm{erfc}}\left ( \frac{\sigma } {a} -\frac{(Y - B)} {\sigma } \right ) }$$

(5)

Median filtering can be applied to the 2D array of background measurements in order to correct for spurious background values. Finally the background map is obtained by a bilinear or a cubic interpolation of the 2D array. The block size is a crucial parameter. If it is too small, the background estimation map will be affected by the presence of objects, and if too large it will not take into account real background variations.

In [6, 15], the local sky level is calculated differently. A 3-sigma clipping around the median is performed in each block. If the standard deviation is changed by less than 20 % in the clipping iterations, the block is uncrowded, and the background level is considered to be equal to the mean of the clipped histogram. Otherwise, it is calculated by c 1 ×median − c 2 ×mean, where $$c_{1} = 3,c_{2} = 2$$ in [15] and $$c_{1} = 2.5,c_{2} = 1.5$$ in [6]. This approach has been preferred to histogram fitting for two reasons: it is more efficient from the computation point of view and more robust with small sample size.


Convolution


In order to optimize the detection, the image must be convolved with a filter. The shape of this filter optimizes the detection of objects with the same shape. Therefore, for star detection, the optimal filter is the PSF. For extended objects, a larger filter size is recommended. In order to have optimal detection for any object size, the detection must be repeated several times with different filter sizes, leading to a kind of multiscale approach.


Detection


Once the image is convolved, all pixels Y [k, l] at location (k, l) with a value larger than T [k, l] are considered as significant, i.e., belonging to an object. T [k, l] is generally chosen as B [k, l] + K σ, where B [k, l] is the background estimate at the same position, σ is the noise standard deviation, and K is a given constant (typically chosen between 3 and 5). The thresholded image is then segmented, i.e., a label is assigned to each group of connected pixels. The next step is to separate the blended objects which are connected and have the same label.

An alternative to the thresholding/segmentation procedure is to find peaks. This is only well suited to star detection and not to extended objects. In this case, the next step is to merge the pixels belonging to the same object.


Deblending/Merging


This is the most delicate step. Extended objects must be considered as single objects, while multiple objects must be well separated. In SExtractor, each group of connected pixels is analyzed at different intensity levels, starting from the highest down to the lowest level. The pixel group can be seen as a surface, with mountains and valleys. At the beginning (highest level), only the highest peak is visible. When the level decreases, several other peaks may become visible, defining therefore several structures. At a given level, two structures may become connected, and the decision whether they form only one (i.e., merging) or several objects (i.e., deblending) must be taken. This is done by comparing the integrated intensities inside the peaks. If the ratio between them is too low, then the two structures must be merged.


Photometry and Classification



Photometry


Several methods can be used to derive the photometry of a detected object [7, 29]. Adaptive aperture photometry uses the first image moment to determine the elliptical aperture from which the object flux is integrated. Kron [29] proposed an aperture size of twice the radius of the first image moment radius r 1, which leads to recovery of most of the flux ( > 90 %). In [6], the value of 2. 5r 1 is discussed, leading to loss of less than 6 % of the total flux. Assuming that the intensity profiles of the faint objects are Gaussian, flux estimates can be refined [6, 35]. When the image contains only stars, specific methods can be developed which take the PSF into account [18, 42].


Star–Galaxy Separation




In the case of star–galaxy classification, following the scanning of digitized images, Kurtz [30] lists the following parameters which have been used:

1.

Mean surface brightness;

 

2.

Maximum intensity and area;

 

3.

Maximum intensity and intensity gradient;

 

4.

Normalized density gradient;

 

5.

Areal profile;

 

6.

Radial profile;

 

7.

Maximum intensity, 2nd and 4th order moments, and ellipticity;

 

8.

The fit of galaxy and star models;

 

9.

Contrast versus smoothness ratio;

 

10.

The fit of a Gaussian model;

 

11.

Moment invariants;

 

12.

Standard deviation of brightness;

 

13.

2nd order moment;

 

14.

Inverse effective squared radius;

 

15.

Maximum intensity and intensity-weighted radius;

 

16.

2nd and 3rd order moments, number of local maxima, and maximum intensity.

 

References for all of these may be found in the cited work. Clearly there is room for differing views on parameters to be chosen for what is essentially the same problem. It is of course the case also that aspects such as the following will help to orientate us towards a particular set of parameters in a particular case: the quality of the data; the computational ease of measuring certain parameters; the relevance and importance of the parameters measured relative to the data analysis output (e.g., the classification, or the planar graphics); and, similarly, the importance of the parameters relative to theoretical models under investigation.


Galaxy Morphology Classification


The inherent difficulty of characterizing spiral galaxies especially when not face-on has meant that most work focuses on ellipticity in the galaxies under study. This points to an inherent bias in the potential multivariate statistical procedures. In the following, it will not be attempted to address problems of galaxy photometry per se [17, 44], but rather to draw some conclusions on what types of parameters or features have been used in practice.

From the point of view of multivariate statistical algorithms, a reasonably homogeneous set of parameters is required. Given this fact, and the available literature on quantitative galaxy morphological classification, two approaches to parameter selection appear to be strongly represented:

1.

The luminosity profile along the major axis of the object is determined at discrete intervals. This may be done by the fitting of elliptical contours, followed by the integrating of light in elliptical annuli [33]. A similar approach was used in the ESO–Uppsala survey. Noisiness and faintness require attention to robustness in measurement: the radial profile may be determined taking into account the assumption of a face-on optically thin axisymmetric galaxy and may be further adjusted to yield values for circles of given radius [64]. Alternatively, isophotal contours may determine the discrete radial values for which the profile is determined [62].

 

2.

Specific morphology-related parameters may be derived instead of the profile. The integrated magnitude within the limiting surface brightness of 25 or 26 mag. arcsec−2 in the visual is popular [33, 61]. The logarithmic diameter (D 26) is also supported by Okamura [43]. It may be interesting to fit to galaxies under consideration model bulges and disks using, respectively, $$r^{\frac{1} {4} }$$ or exponential laws [62], in order to define further parameters. Some catering for the asymmetry of spirals may be carried out by decomposing the object into octants; furthermore, the taking of a Fourier transform of the intensity may indicate aspects of the spiral structure [61].

 

The following remarks can be made relating to image data and reduced data:



  • The range of parameters to be used should be linked to the subsequent use to which they might be put, such as to underlying physical aspects.


  • Parameters can be derived from a carefully constructed luminosity profile, rather than it being possible to derive a profile from any given set of parameters.


  • The presence of both partially reduced data such as luminosity profiles, and more fully reduced features such as integrated flux in a range of octants, is of course not a hindrance to analysis. However, it is more useful if the analysis is carried out on both types of data separately.

Parameter data can be analyzed by clustering algorithms, by principal component analysis, or by methods for discriminant analysis. Profile data can be sampled at suitable intervals and thus analyzed also by the foregoing procedures. It may be more convenient in practice to create dissimilarities between profiles and analyze these dissimilarities: this can be done using clustering algorithms with dissimilarity input.


3 Mathematical Modeling


Different models may be considered to represent the data. One of the most effective is certainly the sparsity model, especially when a specific wavelet dictionary is chosen to represent the data. We introduce here the sparsity concept as well as the wavelet transform decomposition, which is the most used in astronomy.


Sparsity Data Model



A signal $$X,\ X = \left [x_{1},\ldots,x_{N}\right ]^{T}$$, is sparse if most of its entries are equal to zero. For instance, a k-sparse signal is a signal where only k samples have a nonzero value. A less strict definition is to consider a signal as weakly sparse or compressible when only a few of its entries have a large magnitude, while most of them are close to zero.

If a signal is not sparse, it may be sparsified using a given data representation. For instance, if X is a sine, it is clearly not sparse but its Fourier transform is extremely sparse (i.e., 1-sparse). Hence, we say that a signal X is sparse in the Fourier domain if its Fourier coefficients $$\hat{X}[u]$$, $$\hat{X}[u] = \frac{1} {N}\sum _{k=-\infty }^{+\infty }X[k]e^{2i\pi \frac{uk} {N} },$$ are sparse. More generally, we can model a vector signal $$X \in \mathbb{R}^{N}$$ as the linear combination of T elementary waveforms, also called signal atoms: $$X =\boldsymbol{ \Phi }\alpha =\sum _{ i=1}^{T}\alpha [i]\phi _{i}$$, where $$\alpha [i] =\big< X,\phi _{i}\big >$$” src=”/wp-content/uploads/2016/04/A183156_2_En_34_Chapter_IEq9.gif”></SPAN> are called the decomposition coefficients of <SPAN class=EmphasisTypeItalic>X</SPAN> in the dictionary <SPAN id=IEq10 class=InlineEquation><IMG alt= (the N × T matrix whose columns are the atoms normalized to a unit ℓ 2-norm, i.e., $$\forall i \in [1,T],\left \|\phi _{i}\right \|_{\ell^{2}} = 1$$).

Therefore, to get a sparse representation of our data, we need first to define the dictionary $$\boldsymbol{\Phi }$$ and then to compute the coefficients α. x is sparse in $$\Phi $$ if the sorted coefficients in decreasing magnitude have fast decay, i.e., most coefficients α vanish except for a few.

The best dictionary is the one which leads to the sparsest representation. Hence, we could imagine having a huge overcomplete dictionary (i.e., T ≫ N), but we would be faced with prohibitive computation time cost for calculating the α coefficients. Therefore, there is a trade-off between the complexity of our analysis step (i.e., the size of the dictionary) and the computation time. Some specific dictionaries have the advantage of having fast operators and are very good candidates for analyzing the data.

The Isotropic Undecimated Wavelet Transform (IUWT) , also called starlet wavelet transform, is well known in the astronomical domain because it is well adapted to astronomical data where objects are more or less isotropic in most cases [54, 57]. For most astronomical images, the starlet dictionary is very well adapted.


The Starlet Transform



The starlet wavelet transform [53] decomposes an n × n image c 0 into a coefficient set $$W =\{ w_{1},\ldots,w_{J},c_{J}\}$$, as a superposition of the form


$$\displaystyle{c_{0}[k,l] = c_{J}[k,l] +\sum _{ j=1}^{J}w_{ j}[k,l],}$$
where c J is a coarse or smooth version of the original image c 0 and w j represents the details of c 0 at scale 2−j (see Starck et al. [56, 58] for more information). Thus, the algorithm outputs J + 1 sub-band arrays of size N × N. (The present indexing is such that j = 1 corresponds to the finest scale or high frequencies.)

The decomposition is achieved using the filter bank $$(h_{2\text{D}},g_{2\text{D}} =\delta -h_{2\text{D}},\tilde{h}_{2\text{D}} =\delta,\tilde{g}_{2\text{D}} =\delta )$$, where h 2D is the tensor product of two 1D filters h 1D and δ is the Dirac function. The passage from one resolution to the next one is obtained using the “à trous” (“with holes”) algorithm [58]:


$$\displaystyle\begin{array}{rcl} \begin{array}{ll} c_{j+1}[k,l] & =\sum _{m}\sum _{n}h_{1{\mathrm{D}}}[m]\ h_{1{\mathrm{D}}}[n]\ c_{j}[k + 2^{j}m,l + 2^{j}n], \\ w_{j+1}[k,l]& = c_{j}[k,l] - c_{j+1}[k,l]\,\end{array} & &{}\end{array}$$

(6)

If we choose a B 3-spline for the scaling function,


$$\displaystyle\begin{array}{rcl} & \phi (x) = B_{3}(x) = & \\ & \frac{1} {12}(\mid x - 2\mid ^{3} - 4\mid x - 1\mid ^{3} + 6\mid x\mid ^{3} - 4\mid x + 1\mid ^{3} + \mid x + 2\mid ^{3})&{}\end{array}$$

(7)
the coefficients of the convolution mask in one dimension are $$h_{1{\mathrm{D}}} = \left \{ \frac{1} {16}, \frac{1} {4}, \frac{3} {8}, \frac{1} {4}, \frac{1} {16}\right \}$$, and in two dimensions,


$$\displaystyle{h_{2\text{D}} = \left( 1/16 \, 1/4 \, 3/8 \, 1/4 \, 1/16 \right)\left( \begin{array}{c} \frac{1} {16} \\ \frac{1} {4} \\ \frac{3} {8} \\ \frac{1} {4} \\ \frac{1} {16}\end{array} \right) = \left(\begin{array}{ccccc} \frac{1} {256} & \frac{1} {64} & \frac{3} {128} & \frac{1} {64} & \frac{1} {256} \\ \frac{1} {64} & \frac{1} {16} & \frac{3} {32} & \frac{1} {16} & \frac{1} {64} \\ \frac{3} {128} & \frac{3} {32} & \frac{9} {64} & \frac{3} {32} & \frac{3} {128} \\ \frac{1} {64} & \frac{1} {16} & \frac{3} {32} & \frac{1} {16} & \frac{1} {64} \\ \frac{1} {256} & \frac{1} {64} & \frac{3} {128} & \frac{1} {64} & \frac{1} {256} \end{array} \right )}$$

Figure 2 shows the scaling function and the wavelet function when a cubic spline function is chosen as the scaling function ϕ.

A183156_2_En_34_Fig2_HTML.gif


Fig. 2
Left, the cubic spline function ϕ; right, the wavelet ψ

The most general way to handle the boundaries is to consider that $$c[k + N] = c[N - k]$$ (“mirror”). But other methods can be used such as periodicity $$(c[k + N] = c[N])$$ or continuity $$(c[k + N] = c[k])$$.

The starlet transform algorithm is:

1.

We initialize j to 0 and we start with the data c j [k, l].

 

2.

We carry out a discrete convolution of the data c j [k, l] using the filter (h 2D), using the separability in the two-dimensional case. In the case of the B 3-spline, this leads to a row-by-row convolution with $$\left ( \frac{1} {16}, \frac{1} {4}, \frac{3} {8}, \frac{1} {4}, \frac{1} {16}\right )$$, followed by column-by-column convolution. The distance between the central pixel and the adjacent ones is 2 j .

 

3.

After this smoothing, we obtain the discrete wavelet transform from the difference $$c_{j}[k,l] - c_{j+1}[k,l]$$.

 

4.

If j is less than the number J of resolutions we want to compute, we increment j and then go to step 2.

 

5.

The set $$\alpha =\{ w_{1},\ldots,w_{J},c_{J}\}$$ represents the wavelet transform of the data.

 

This starlet transform is very well adapted to the detection of isotropic features, and this explains its success for astronomical image processing, where the data contain mostly isotropic or quasi-isotropic objects, such as stars, galaxies, or galaxy clusters.


Figure 3 shows the starlet transform of the galaxy NGC 2997 displayed in Fig. 4. Five wavelet scales and the final smoothed plane (lower right) are shown. The original image is given exactly by the sum of these six images.

A183156_2_En_34_Fig3_HTML.jpg


Fig. 3
Wavelet transform of NGC 2997 by the IUWT. The co-addition of these six images reproduces exactly the original image


A183156_2_En_34_Fig4_HTML.jpg


Fig. 4
Galaxy NGC 2997


The Starlet Reconstruction



The reconstruction is straightforward. A simple co-addition of all wavelet scales reproduces the original map: $$c_{0}[k,l] = c_{J}[k,l] +\sum _{ j=1}^{J}w_{j}[k,l]$$. But because the transform is non-subsampled, there are many ways to reconstruct the original image from its wavelet transform [53]. For a given wavelet filter bank (h, g), associated with a scaling function ϕ and a wavelet function ψ, any synthesis filter bank $$(\tilde{h},\ \tilde{g})$$, which satisfies the following reconstruction condition


$$\displaystyle\begin{array}{rcl} \hat{h}^{{\ast}}(\nu )\hat{\tilde{h}}(\nu ) +\hat{ g}^{{\ast}}(\nu )\hat{\tilde{g}}(\nu ) = 1\,& &{}\end{array}$$

(8)
leads to exact reconstruction. For instance, for isotropic h, if we choose $$\tilde{h} = h$$ (the synthesis scaling function $$\tilde{\phi } =\phi$$), we obtain a filter $$\tilde{g}$$ defined by [53]


$$\displaystyle{\tilde{g} =\delta +h.}$$
If h is a positive filter, then g is also positive. For instance, if $$h_{\mathrm{1D}} = [1,4,6,4,1]/16$$, then $$\tilde{g}_{\mathrm{1D}} = [1,4,22,4,1]/16$$. That is, $$\tilde{g}_{\mathrm{1D}}$$ is positive. This means that $$\tilde{g}$$ is no longer related to a wavelet function. The 1D detail synthesis function related to $$\tilde{g}_{\mathrm{1D}}$$ is defined by


$$\displaystyle\begin{array}{rcl} \frac{1} {2}\tilde{\psi }_{\mathrm{1D}}\left ( \frac{t} {2}\right )& =& \phi _{\mathrm{1D}}(t) + \frac{1} {2}\phi _{\mathrm{1D}}\left ( \frac{t} {2}\right ).{}\end{array}$$

(9)

Note that by choosing $$\tilde{\phi }_{\mathrm{1D}} =\phi _{\mathrm{1D}}$$, any synthesis function $$\tilde{\psi }_{\mathrm{1D}}$$ which satisfies


$$\displaystyle\begin{array}{rcl} \hat{\tilde{\psi }}_{\mathrm{1D}}(2\nu )\hat{\psi }_{\mathrm{1D}}(2\nu ) = \hat{\phi }_{\mathrm{1D}}^{2}(\nu ) -\hat{\phi }_{\mathrm{ 1D}}^{2}(2\nu )& &{}\end{array}$$

(10)
leads to an exact reconstruction [36] and $$\hat{\tilde{\psi }}_{\mathrm{1D}}(0)$$ can take any value. The synthesis function $$\tilde{\psi }_{\mathrm{1D}}$$ does not need to verify the admissibility condition (i.e., to have a zero mean).

Figure 5 shows the two functions $$\tilde{\phi }_{\mathrm{1D}}(=\phi _{\mathrm{1D}})$$ and $$\tilde{\psi }_{\mathrm{1D}}$$ used in the reconstruction in 1D, corresponding to the synthesis filters $$\tilde{h}_{\mathrm{1D}} = h_{\mathrm{1D}}$$ and $$\tilde{g}_{\mathrm{1D}} =\delta +h_{\mathrm{1D}}$$. More details can be found in [53].

A183156_2_En_34_Fig5_HTML.gif


Fig. 5
Left, $$\tilde{\phi }_{\mathrm{1D}}$$ the 1D synthesis scaling function and right, $$\tilde{\psi }_{\mathrm{1D}}$$ the 1D detail synthesis function


Starlet Transform: Second Generation



A particular case is obtained when $$\hat{\tilde{\phi }}_{\mathrm{1D}} =\hat{\phi } _{\mathrm{1D}}$$ and $$\hat{\psi }_{\mathrm{1D}}(2\nu ) = \frac{\hat{\phi }^{2}_{ \mathrm{1D}}(\nu )-\hat{\phi }_{\mathrm{1D}}^{2}(2\nu )} {\hat{\phi }_{\mathrm{1D}}(\nu )}$$, which leads to a filter g 1D equal to δ − h 1D ⋆ h 1D. In this case, the synthesis function $$\tilde{\psi }_{\mathrm{1D}}$$ is defined by $$\frac{1} {2}\tilde{\psi }_{\mathrm{1D}}( \frac{t} {2}) =\phi _{\mathrm{1D}}(t)$$, and the filter $$\tilde{g}_{\mathrm{1D}} =\delta$$ is the solution to (8).

We end up with a synthesis scheme where only the smooth part is convolved during the reconstruction.

Deriving h from a spline scaling function, for instance $$B_{1}(h_{1} = [1,2,1]/4)$$ or $$B_{3}\ (h_{3} = [1,4,6,4,1]/16)$$ (note that $$h_{3} = h_{1} \star h_{1}$$), since h 1D is even-symmetric (i.e., $$H(z) = H(z^{-1})$$), the z-transform of g 1D is then


$$\displaystyle\begin{array}{rcl} \begin{array}{ll} G(z)& = 1 - H^{2}(z) = 1 - z^{4}\left (\frac{1+z^{-1}} {2} \right )^{8} \\ & = \frac{1} {256}\left (-z^{4} - 8z^{3} - 28z^{2} - 56z + 186 - 56z^{-1} - 28z^{-2} - 8z^{-3} - z^{-4}\right ),\end{array}& &{}\end{array}$$

(11)
which is the z-transform of the filter


$$\displaystyle{g_{\mathrm{1D}} = [-1,-8,-28,-56,186,-56,-28,-8,-1]/256.}$$

We get the following filter bank:


$$\displaystyle\begin{array}{rcl} \begin{array}{ll} h_{\mathrm{1D}} & = h_{3} = \tilde{h} = [1,4,6,4,1]/16 \\ g_{\mathrm{1D}} & =\delta -h \star h = [-1,-8,-28,-56,186,-56,-28,-8,-1]/256\\ \tilde{g}_{ \mathrm{1D}} & =\delta .\end{array} & &{}\end{array}$$

(12)

The second-generation starlet transform algorithm is:

1.

We initialize j to 0 and we start with the data c j [k].

 

2.

We carry out a discrete convolution of the data c j [k] using the filter h 1D. The distance between the central pixel and the adjacent ones is 2 j . We obtain c j+1[k].

 

3.

We do exactly the same convolution on c j+1[k] and we obtain c j+1 ′ [k].

 

4.

After this two-step smoothing, we obtain the discrete starlet wavelet transform from the difference $$w_{j+1}[k] = c_{j}[k] - c_{j+1}^{{\prime}}[k]$$.

 

5.

If j is less than the number J of resolutions we want to compute, we increment j and then go to step 2.

 

6.

The set $$\alpha =\{ w_{1},\ldots,w_{J},c_{J}\}$$ represents the starlet wavelet transform of the data.

 

As in the standard starlet transform, extension to 2D is trivial. We just replace the convolution with h 1D by a convolution with the filter h 2D, which is performed efficiently by using the separability.

With this specific filter bank, there is a no convolution with the filter $$\tilde{g}_{\mathrm{1D}}$$ during the reconstruction. Only the low-pass synthesis filter $$\tilde{h}_{\mathrm{1D}}$$ is used.

The reconstruction formula is


$$\displaystyle\begin{array}{rcl} c_{j}[l] = (h_{\mathrm{1D}}^{(j)} \star c_{ j+1})[l] + w_{j+1}[l]\,& &{}\end{array}$$

(13)
and denoting $$L^{j} = h^{(0)} \star \ldots \star h^{(j-1)}$$ and L 0 = δ, we have


$$\displaystyle\begin{array}{rcl} c_{0}[l] = \left (L^{J} \star c_{ J}\right )[l] +\sum _{ j=1}^{J}\left (L^{j-1} \star w_{ j}\right )[l].& &{}\end{array}$$

(14)
Each wavelet scale is convolved with a low-pass filter.

The second-generation starlet reconstruction algorithm is:

1.

The set $$\alpha =\{ w_{1},\ldots,w_{J},c_{J}\}$$ represents the input starlet wavelet transform of the data.

 

2.

We initialize j to J − 1 and we start with the coefficients c j [k].

 

3.

We carry out a discrete convolution of the data c j+1[k] using the filter (h 1D). The distance between the central pixel and the adjacent ones is 2 j . We obtain c j+1 ′ [k].

 

4.

Compute $$c_{j}[k] = c_{j+1}^{{\prime}}[k] + w_{j+1}[k]$$.

 

5.

If j is larger than 0, $$j = j - 1$$ and then go to step 3.

 

6.

c 0 contains the reconstructed data.

 

As for the transformation, the 2D extension consists just in replacing the convolution by h 1D with a convolution by h 2D.

Figure 6 shows the analysis scaling and wavelet functions. The synthesis functions $$\tilde{\phi }_{\mathrm{1D}}$$ and $$\tilde{\psi }_{\mathrm{1D}}$$ are the same as those in Fig. 5. As both are positive, we have a decomposition of an image X on positive scaling functions $$\tilde{\phi }_{\mathrm{1D}}$$ and $$\tilde{\psi }_{\mathrm{1D}}$$, but the coefficients α are obtained with the starlet wavelet transform and have a zero mean (except for c J ), as a regular wavelet transform.

A183156_2_En_34_Fig6_HTML.gif


Fig. 6
Left, the ϕ 1D analysis scaling function and right, the ψ 1D analysis wavelet function. The synthesis functions $$\tilde{\phi }_{\mathrm{1D}}$$ and $$\tilde{\psi }_{\mathrm{1D}}$$ are the same as those in Fig. 5

In 2D, similarly, the second-generation starlet transform leads to the representation of an image X[k, l]:


$$\displaystyle\begin{array}{rcl} X[k,l] =\sum _{m,n}\phi _{j,k,l}^{(1)}(m,n)c_{ J}[m,n] +\sum _{ j=1}^{J}\sum _{ m,n}\phi _{j,k,l}^{(2)}(m,n)w_{ j}[m,n]\,& &{}\end{array}$$

(15)
where $$\phi _{j,k,l}^{(1)}(m,n) = 2^{-2j}\tilde{\phi }_{\mathrm{1D}}(2^{-j}(k - m))\tilde{\phi }_{\mathrm{1D}}(2^{-j}(l - n))$$ and $$\phi _{j,k,l}^{(2)}(m,n) = 2^{-2j}\tilde{\psi }_{\mathrm{1D}}(2^{-j}(k - m))\tilde{\psi }_{\mathrm{1D}}(2^{-j}(l - n))$$.

ϕ (1) and ϕ (2) are positive, and w j are zero mean 2D wavelet coefficients.

The advantage of the second-generation starlet transform will be seen in section “Sparse Positive Decomposition” below.


Sparse Modeling of Astronomical Images



Using the sparse modeling, we now consider that the observed signal X can be considered as a linear combination of a few atoms of the wavelet dictionary $$\boldsymbol{\Phi } = [\phi _{1},\ldots,\phi _{T}]$$. The model of Eq. 3 is then replaced by the following:


$$\displaystyle\begin{array}{rcl} Y & =& H\boldsymbol{\Phi }\alpha + N + B{}\end{array}$$

(16)
and $$X =\boldsymbol{ \Phi }\alpha$$, and $$\alpha = \left \{w_{1},\ldots,w_{J},c_{J}\right \}$$. Furthermore, most of the coefficients α will be equal to zero. Positions and scales of active coefficients are unknown, but they can be estimated directly from the data Y. We define the multiresolution support M of an image Y by


$$\displaystyle\begin{array}{rcl} M_{j}[k,l] = \left \{\begin{array}{ll} \text{ 1 }&\text{ if }w_{j}[k,l]\text{ is significant} \\ \text{ 0 }&\text{ if }w_{j}[k,l]\text{ is not significant} \end{array} \right.& &{}\end{array}$$

(17)
where w j [k, l] is the wavelet coefficient of Y at scale j and at position (k, l). Hence, M describes the set of active atoms in Y. If H is compact and not too extended, then M describes also well the active set of X. This is true because the background B is generally very smooth, and therefore, a wavelet coefficient w j [k, l] of Y, which does not belong to the coarsest scale, is only dependent on X and N (the term < ϕ i , B > being equal to zero).


Selection of Significant Coefficients Through Noise Modeling


We need now to determine when a wavelet coefficient is significant. Wavelet coefficients of Y are corrupted by noise, which follows in many cases a Gaussian distribution, a Poisson distribution, or a combination of both. It is important to detect the wavelet coefficients which are “significant,” i.e., the wavelet coefficients which have an absolute value too large to be due to noise.

For Gaussian noise, it is easy to derive an estimation of the noise standard deviation σ j at scale j from the noise standard deviation, which can be evaluated with good accuracy in an automated way [55]. To detect the significant wavelet coefficients, it suffices to compare the wavelet coefficients w j [k, l] to a threshold level t j . t j is generally taken equal to K σ j , and K, as noted in Sect. 2, is chosen between 3 and 5. The value of 3 corresponds to a probability of false detection of 0. 27 %. If w j [k, l] is small, then it is not significant and could be due to noise. If w j [k, l] is large, it is significant:


$$\displaystyle\begin{array}{rcl} \begin{array}{l} \text{ if }\mid w_{j}[k,l]\mid \ \geq \ t_{j}\ \ \text{ then }w_{j}[k,l]\text{ is significant } \\ \text{ if }\mid w_{j}[k,l]\mid \ < t_{j}\ \ \text{ then }w_{j}[k,l]\text{ is not significant }\end{array} & &{}\end{array}$$

(18)

When the noise is not Gaussian, other strategies may be used:

Tags:
Apr 9, 2016 | Posted by in GENERAL RADIOLOGY | Comments Off on Transform in Astronomical Data Processing

Full access? Get Clinical Tree
Get Clinical Tree app for offline access