{'fraction_non_alphanumeric': 0.03787558188743123, 'fraction_numerical': 0.010050782903089294, 'mean_word_length': 4.511953352769679, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Morphological components analysis for circumstellar disks imaging

Recent developments in astronomical observations enable direct imaging of circumstellar disks. Precise characterization of such extended structure is essential to our understanding of stellar systems. However, the faint intensity of the circumstellar disks compared to the brightness of the host star compels astronomers to use tailored observation strategies, in addition to state-of-the-art optical devices. Even then, extracting the signal of circumstellar disks heavily relies on post-processing techniques. In this work, we propose a morphological component analysis (MCA) approach that leverages low-complexity models of both the disks and the stellar light corrupting the data. In addition to disks, our method allows to image exoplanets. Our approach is tested through numerical experiments.


Introduction
Direct imaging of stellar systems is a challenging task that requires hardware with high contrast and high resolution to be able to detect the faint objects located near the extremely bright star. Ground based telescope achieve the highest resolution and thanks to adaptive optics (AO), they are able to overcome the atmospheric turbulence. The brightness of the star is dimmed using a coronagraph. However, even with state-of-the-art hardware, residual quasi-statics non-common path aberrations, form speckles in the observations [5]. The presence of these speckles prevent the detection of on-sky signals such as disks and exoplanets. Specific observation strategies and postprocessing techniques have been used to force diversity within the data and leverage this diversity to improve the separability between speckles and on-sky signals.
Angular differential imaging (ADI) is a popular observation strategy that leverages the rotation of the Earth to introduce such a data diversity. In this setting, the star is kept in the center of the field of view while snapshots of the stellar systems are taken during the observation period (i.e., a few hours). As most of the residual speckles (after AO) are due to the telescope itself, they remain quasi-static in the images, while on-sky signals follow a deterministic circular trajectory, determined by known parallactic angles [10].
The typical ADI data processing pipeline is guided by a principal component analysis (PCA) [1,20]: (i) the T × (n × n) spatiotemporal data cube (with T the number of frames and n 2 > T the number of pixels of each frame) is reshaped into a R T ×n 2 matrix Y , (ii) its rank-r approximation L is computed by thresholding its singular values decomposition (SVD) to its r largest singular values, (iii) L is then subtracted from the data to form S = Y − L containing the on-sky signals, and (iv) the frames of S are * BP and LJ are funded by the Belgian F.R.S.-FNRS.
aligned (by rotation and interpolation) to the on-sky coordinates and temporally averaged to form the processed frame. Objects detection can then be done by hypothesis testing [11,14].
Morphological Limitations: The morphology of the circumstellar disks is known to be severely distorted by the PCA pipeline, hindering our capability to study disks structures from ADI datasets. Results of PCA on the ellipsoidal disk surrounding HR 4796A [13] is displayed on Fig. 1 (left), where we can see unphysical artifacts with negative intensity. Among the few attempts to remove these artifacts, Milli et al. [12] reduced them by injecting forward-modeled disks. In [18], non-negative matrix factorization was also used to reduce the artifacts.

Our approach
We propose to recast the stellar system imaging task from ADI dataset as a MCA task [2,3,21]. To achieve this, we first present the acquisition model and then propose a constrained convex optimization solving our MCA problem. We then list the main physical priors our convex optimization is based upon.
Acquisition model: We restart from an ADI sequence Y ∈ R T ×n 2 , assumed to contain a disk. We model Y as the sum of two terms: the starlight Y * and the rotating on-sky signal Y . As the star is far from the Earth, it is point-like source and ideally, its intensity is blocked by the coronagraph. However, as discussed earlier, atmospheric turbulence and imperfections in the optics introduce speckles in the observation. These speckles are modeled as the sum of two terms, encoding their temporal behavior: a static termL and a non-static term N s : whereL is assumed to be a rank-r matrix, r t.
Concerning the on-sky component, its light intensity being small, we can neglect the non-common path aberrations. We thus assume its intensity to be constant through time and we model Y as a single rotating imagex ∈ R n 2 where p ∈ R T stands for the time variation intensity of Y and R : R T ×n 2 → R T ×n 2 is the linear operator that rotates each frame of the volume according to the parallactic angles. Furthermore, because the light is diffracted as it enters the telescope, the on-sky signal is blurred by the known telescope PSF ϕ. We write this as where * denotes the 2D convolution applied separately on each image of R(p(x) .

arXiv:2101.12706v1 [astro-ph.IM] 29 Jan 2021
The final acquisition model of Y reads Source separation algorithm: Given the acquisition model, our MCA algorithm is performed by solving the following convex optimization problem Each elements of the problem (5) is justified by physical priors which are listed below. Fidelity term (5a): As shown in [14], the speckle noise follows a sub-exponential distribution which implies P (N s > ) exp −c min( 2 /δ 2 , /δ) . This indicates that the negative log-likelihood of the speckle is approached by the Huber-loss function [7] defined as which is well-suited for the fidelity term. The metric · H δ is simply defined as a H δ = i |a i | δ , for a given vector a.
Static part of the speckles (5b): As it is common in the literature of background-foreground separation [see, e.g., 23], we use a low-rank structure to model the static part of the speckle field, e.g., impose L * τ L . However, the large intensity discrepancy between the disk and the star makes the use of the nuclear norm impractical. Indeed, in a dummy low-rank plus sparse problem [16], we observed that, as the intensity of the low-rank component increases, the quality of the estimation of the sparse component becomes overly sensitive to the accuracy in the estimation of τ L . For ADI datasets, the intensity ofL is typically around 10 3 larger than the intensity ofx. In this case, even when τ L is only a percent above or below the groundtruth, the recovery of x is unsuccessful. Fortunately, the greedy algorithm presented in [15] yields a decent estimate L * ofL. Hence inspired by the work of [4], we relax the low-rankness imposed onL by forcing it to lie in the span of the first r columns of U * , where U * is given by the SVD decomposition of L * , as shown in (5b).
Spatial structure of the disks (5c): the disk is regularized by promoting its sparsity in the shearlets domain, as they have been shown to be efficient for sparsely representing multivariate data containing edges and curved structures [8].
Include exoplanetary signal (5d): Following the premise of the MCA algorithm, we leverage the morphological diversity between the disk and the exoplanets to separate these two sources x = x d + x p . After deconvolution, the exoplanetary signal x p has an optimally sparse representation in the direct domain, whereas its shearlets representation involves more terms.
Positivity of the images (5e): since the signals of interest are positive, we finally enforce that L, x d , x p 0.

Numerical experiments:
We solve (5) with the Primal-Dual Three-Operator splitting [ PD3O 22]. The gradient of (5a) requires the computation of R , which is done using the autograd functionality of PyTorch [17] along with kornia, a PyTorchbased computer vision toolbox [19]. For the shearlets transform, we use the python version of the ShearLab 3D toolbox [9]). Fig. 1 (right) shows the HR 4796A image obtained our method. We can see that the recovered throughput is about 8 times that of PCA. The overall shape of the disk is also more physically sound, e.g., there is no negative intensities and the shape is closer to estimations obtained with physical models (see for instance [13]).
To illustrate the capability of our algorithm to faithfully recover both a disk and a exoplanet, we used the VIP toolbox [6] to create a diskx d with a exoplanetx p displayed in Fig. 2 (left). Both thex d andx p were convolved with the telescope PSF before being injected in an empty ADI cube. We used our approach to recoverx d andx p , the result is shown in Fig. 2 (right). We can see that our method was able to reproduce the shape of the disk faithfully and to recover the exoplanetary signal.

Conclusion
We presented a MCA framework to image circumstellar disks and exoplanets from ADI datasets. Our method leverages physical knowledge to produce faithful images. To best of our knowledge, our method is the first to include MCA and deconvolution in ADI post-processing, allowing to disentangle exoplanets from circumstellar disks. Although our method is able to recover on-sky signals from ADI data, Fig. 2 (right) also features intensities in the center that do not correspond to the injected signal. Future work should include a procedure to asses the quality of the output, as for instance, adapted hypothesis testing.