<figure><img src="image_1.png" /></figure>

HAL
open science

# XX. Haute Couture: Spectral stitching of JWST MIRI-IFU cubes with matrix completion

Amelie Canin, Cédric Févotte, Nicolas Dobigeon, Dries van de Putte, Takashi Onaka, Olivier Berné

► To cite this version:

Amelie Canin, Cédric Févotte, Nicolas Dobigeon, Dries van de Putte, Takashi Onaka, et al. XX. Haute Couture: Spectral stitching of JWST MIRI-IFU cubes with matrix completion. 2025. hal-05399060

HAL Id: hal-05399060

https://hal.science/hal-05399060v1

Preprint submitted on 4 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# PDRs4All
## XX. Haute Couture: Spectral stitching of JWST MIRI-IFU cubes with matrix completion

Amélie Canin<fnref n="1" />, Cédric Févotte<fnref n="2" />, Nicolas Dobigeon<fnref n="2" />, Dries Van De Putte<fnref n="3" />, Takashi Onaka<fnref n="4" />, and Olivier Berné<fnref n="1" />

## ABSTRACT

MIRI is the imager and spectrograph covering wavelengths from 4.9 to 27.9 µm onboard the James Webb Space Telescope (JWST). The Medium-Resolution Spectrometer (MRS) consists of four integral field units (IFU), each of which has three sub-channels. The twelve resulting spectral data cubes have different fields of view, spatial, and spectral resolutions. The wavelength range of each cube partially overlaps with the neighboring bands, and the overlap regions typically show flux mismatches which have to be corrected by spectral stitching methods. Stitching methods aim to produce a single data cube incorporating the data of the individual sub-channels, which requires matching the spatial resolution and the flux discrepancies. We present Haute Couture, a novel stitching algorithm which uses non-negative matrix factorization (NMF) to perform a matrix completion, where the available MRS data cubes are treated as twelve sub-matrices of a larger incomplete matrix. Prior to matrix completion, we also introduce a novel pre-processing to homogenize the global intensities of the twelve cubes. Our pre-processing consists in jointly optimizing a set of global scale parameters that maximize the fit between the cubes where spectral overlap occurs. We apply our novel stitching method to JWST data obtained as part of the PDRs4All observing program of the Orion Bar, and produce a uniform cube reconstructed with the best spatial resolution over the full range of wavelengths.

## 1. Introduction

The James Webb Space Telescope (JWST, Gardner et al. (2006)) is an infrared telescope launched in 2021. Four instruments are onboard: the Mid-Infrared Instrument (MIRI, Rieke et al. 2015; Wright et al. 2023), the Near Infrared Camera (NIRCam), the Near Infrared Spectrograph (NIRSpec), and the Near Infrared Imager and Slitless Spectrograph (NIRISS). MIRI provides images and spectroscopic data through different observing modes covering 4.9 to 27.9 µm. The Medium-Resolution Spectroscopy mode of MIRI (MRS, Wells et al. 2015; Argyriou et al. 2023) uses four Integral Field Units (IFUs), each one covering a portion of the wavelength range referred to as channel. Each of these channels is subdivided into three sub-bands (short, medium and long). Consecutive channels and bands share some information thanks to spectral overlaps. Besides, along the channels, the field of view increases whereas the spatial resolution decreases. The Science Calibration Pipeline provided by the Space Telescope Science Institute thus produces a total of twelve data cubes with different spatial coverages and resolutions. Another challenge raised by the MIRI data acquisition process results from intensity gaps which affect the spectral measurements between adjacent cubes. These intensity gaps, visible from typical spectra reproduced in Figure 1, may be due to calibration mismatches. To make MIRI data easier to use, a convenient data product would be a singular data cube free of spectral intensity gaps, spanning the full spectral range, covering a common field of view, which preserves the highest spatial resolution attainable in each band.

The process of assembling individual spectral data cubes is often referred to as *stitching*. Because of the aforementioned challenges resulting from the unconventional acquisition im-

plemented by MIRI, stitching raises the following questions i) how to assemble intensity gap-free spectra from the twelve data cubes? and ii) how to deal with the distinct spatial sampling (which varies with $1/\lambda$)? A commonly employed stitching method consists in reprojecting all the spectral data onto the poorest resolution spatial grid. This was for instance applied on data from the Infrared Spectrograph (Houck et al. 2004) taken during the Spitzer mission (Werner et al. 2004). This was done by first convolving all the individual cubes with the point spread function (PSF) of the largest wavelength (corresponding to the lowest resolution) in order to reach the same spatial resolution over the full spectral range (see, e.g., Berné et al. 2007; Sandstrom et al. 2010). This procedure, referred to as “coarse stitching” hereafter, is illustrated in Figure 2.a when considering two data cubes. Unfortunately, this strategy leads to a regrettable loss of valuable spatial information.

To overcome the limitations inherent to coarse stitching strategy, this paper presents a novel method, coined as *Haute Couture*, that enables stitching while preserving spatial resolution. The proposed method performs the reconstruction of a full data cube at the best spatial resolution over the full spectral range, as illustrated in Figure 2.b. We first show that the data recorded in the twelve sub-channels can be rearranged as twelve sub-matrices of a larger matrix with missing values. Interestingly, recovering these missing values amounts to reconstructing a data cube at the finest spatial resolution over the full spectral range. *Haute Couture* exploits the overlapping spectral information between adjacent channels and bands to frame the stitching task as a *matrix completion* problem, which can then be solved by *nonnegative matrix factorization* (NMF). Besides, to avoid the issues caused by intensity gaps between consecutive chan-

<footnote n="1">¹ Institut de Recherche en Astrophysique et Planétologie (IRAP), CNRS, CNES, Université de Toulouse, France</footnote>

<footnote n="2">² Institut de Recherche en Informatique de Toulouse (IRIT), CNRS, Toulouse INP, Université de Toulouse, France</footnote>

<footnote n="3">³ University of Western Ontario, Canada</footnote>

<footnote n="4">⁴ Department of Astronomy, Graduate School of Science, The University of Tokyo, Japan</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Examples of intensity gaps observed in overlapping areas of adjacent cubes, obtained from MIRI-MRS data collected within the PDRs4All program (see Section 3.1 for more details). The spectra displayed correspond to the average spectra over the entire field of view for two overlapping sub-channels. Panel a: channel 2-medium and 2-long. Panel b: channel 2-long and 3-short. Panel c: channel 3-short and 3-medium.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Illustration of the stitching problem for two cubes. On the left, we observe that the shortest-wavelength cube in blue has a better spatial resolution (as illustrated by the finer grid) and smaller field of view than the longest-wavelength cube in red. The space × wavelength supports of the cubes partially overlap. The coarse stitching procedure shown in Figure 1.a sacrifices spatial resolution, while *Haute Couture* shown in Figure 1.b enables the reconstruction of a cube with best spatial resolution over the full range of wavelengths.</figcaption></figure>

nels and bands, prior to matrix completion, *Haute Couture* applies a novel pre-processing to homogenize the global intensities along the twelve cubes. This pre-processing consists of jointly optimizing a set of global scaling parameters to maximize the fit between the cubes where spectral overlap occurs. It is worth noting that this pre-processing is an independent procedure that can be combined with any stitching method.

Following the presentation of methodological compounds of *Haute Couture* in Section 2, MIRI-MRS data observed as part of the PDRs4All program on the Orion Bar are presented in Section 3 and stitched in Section 4. In particular, we demonstrate the ability of *Haute Couture* to reconstruct a uniform cube with the best spatial resolution over the full range of wavelengths.

## 2. Haute Couture

### 2.1. Spectral stitching as a matrix completion problem

As introduced in Section 1, MIRI-MRS acquires spectra in four channels indexed by $c$ ($c = 1, \dots, 4$). Each channel is divided

<table><caption>**Table 1.** Name and wavelength range of the MIRI-MRS channels and sub-channels.</caption><thead><tr><th></th><th>Sub-channel</th><th>(c, b)</th><th>λ-range (μm)</th></tr></thead><tbody><tr><td rowspan="3">channel 1</td><td>short</td><td>(1, s)</td><td>4.90 - 5.74</td></tr><tr><td>medium</td><td>(1, m)</td><td>5.66 - 6.63</td></tr><tr><td>long</td><td>(1, l)</td><td>6.53 - 7.65</td></tr><tr><td rowspan="3">channel 2</td><td>short</td><td>(2, s)</td><td>7.51 - 8.77</td></tr><tr><td>medium</td><td>(2, m)</td><td>8.67 - 10.13</td></tr><tr><td>long</td><td>(2, l)</td><td>10.01 - 11.70</td></tr><tr><td rowspan="3">channel 3</td><td>short</td><td>(3, s)</td><td>11.55 - 13.47</td></tr><tr><td>medium</td><td>(3, m)</td><td>13.34 - 15.57</td></tr><tr><td>long</td><td>(3, l)</td><td>15.41 - 17.98</td></tr><tr><td rowspan="3">channel 4</td><td>short</td><td>(4, s)</td><td>17.70 - 20.95</td></tr><tr><td>medium</td><td>(4, m)</td><td>20.69 - 24.48</td></tr><tr><td>long</td><td>(4, l)</td><td>14.40 - 27.90</td></tr></tbody></table>

into three sub-channels, referred to as short, medium and long and shortened as s, m and l, respectively. The spectral ranges of each sub-channel are reported in Table 1. Each sub-channel produces a space × wavelength cube as schematized in Figure 2 where two of these cubes are depicted in blue and red. In the remainder of the paper, the spatial organization of the data is ignored and we unfold each cube with respect to the spatial dimension (i.e., pixels). Thus the data associated with each cube can be represented as a matrix indexed by pixels and wavelengths. Formally, the matrix which gathers the data collected in sub-channel $b \in \{s, m, l\}$ of channel $c$ ($c = 1, \dots, 4$) is denoted as $X^{c,b} \in \mathbb{R}_{+}^{\Lambda^{c,b} \times P^{c,b}}$ where its rows index wavelengths while its columns index pixels. Specifically, the matrix $X^{c,b}$ associated with the sub-channel $(c, b)$ gathers spectra composed of $\Lambda^{c,b}$ wavelengths and acquired over $P^{c,b}$ spatial pixels.

The resulting twelve sub-matrices can be arranged according to a suitable manner to form a larger matrix $X \in \mathbb{R}_{+}^{\Lambda \times P}$, as represented in Figure 3, where $P = \sum_{c=1}^{4} \sum_{b \in \{s,m,l\}} P^{c,b}$ and $\Lambda$ is the total number of wavelengths over the spectral range covered by MIRI-MRS. It is worth noting that the spectral range of adjacent sub-matrices overlap for a few wavelengths, i.e., $\Lambda \neq \sum_{c=1}^{4} \sum_{b \in \{s,m,l\}} \Lambda^{c,b}$. This results in an almost block diagonal structure of the matrix $X$. In this figure, the white parts of the matrix correspond to unavailable (i.e., unobserved) data, specified by the symbol $\emptyset$. Because of these missing coefficients, the matrix $X$ is said to be *incomplete*.

The key rationale of the proposed stitching approach lies on the following insight: recovering all the missing coefficients of this matrix $X$ would provide a comprehensive nay redundant spatial-spectral description of the scene observed by MIRI-MRS. In particular, recovering the values of the submatrix $S$ corresponding to the blue shaded area in Figure 3 would produce a data cube at the spatial resolution of the channel 1-long (i.e., $c = 1$ and $b = 1$) over the full spectral extent offered by MIRI-MRS, which is the main objective of this work. In other words, the spectral stitching task can be formulated as a *matrix completion* problem, i.e., recovering the missing data in the matrix $X$. In this work, we target the channel 1-long to define the spatial and spectral resolutions of the stitched data because this leads to the best trade-off between high spatial resolution and high signal-to-noise ratio (SNR). However this arbitrary choice can be lifted since the proposed method will complete the full matrix $X$.

The following paragraphs details the proposed two-step procedure to stitch the MIRI-MRS data while accommodating for intensity gaps between adjacent data cubes. This procedure con-
<figure><img src="image_5.png" /><figcaption>Fig. 3. The twelve unfolded cubes provided by the twelve sub-channels can be organized as sub-matrices of a larger matrix X with an almost block-diagonal structure. The wavelength range of consecutive blocks overlap for a few rows, corresponding to spectral bands observed by consecutive sub-channels. Blue blocks represent channel 1, green ones channel 2, yellow ones channel 3 and red ones channel 4. White and light blue parts correspond to unobserved coefficients. The light blue area corresponds to the area we are specifically interested in reconstructing.</figcaption></figure>

sists in first rescaling the individual sub-matrices $X^{c,b}$ homogeneously to correct the intensity gaps (see Section 2.2) and then filling the matrix X composed of the resulting intensity-corrected sub-matrices (see Section 2.3).

## 2.2. Intensity gaps correction

To mitigate the intensity gaps highlighted in Figure 1, we introduce a novel pre-processing to homogenize the global intensities of the twelve cubes. This pre-processing consists in jointly adjusting a set of global scale parameters to maximize the fit between the parts of the cubes where spectral overlaps occur. Let us denote by $\alpha^{c,b}$ the scale parameter to be applied to the entries of the matrix $X^{c,b}$. As an example, referring to Figure 3, the scale parameter, say, $\alpha^{2,1}$, is computed by maximizing the fit between, on the one hand, the first rows of $X^{2,1}$ and the last rows of $X^{2,m}$ and, in the other hand, the last rows of $X^{2,1}$ and the first rows of $X^{3,s}$. Because the matrices $X^{2,m}$ and $X^{3,s}$ should be themselves rescaled by the unknown scale parameters $\alpha^{2,m}$ and

$\alpha^{3,s}$, the optimal values of the whole set of scale parameters are inter-dependent and they have to be adjusted jointly. Thankfully, when the fit is measured by the squared Euclidean distance between the portions of shared spectra, this optimization problem has a closed-form solution.

In practice at least one sub-channels is chosen as a reference to make the problem well-posed (otherwise, a trivial but inappropriate solution would be setting $\alpha^{c,b} = 0$ for all channels and sub-channels). In the experimental results conducted in Section 3, we fix $\alpha^{1,s} = \alpha^{4,1} = 1$, and maximize the fit with respect to the ten remaining scale parameters. This amounts to trusting the calibration of the data observed in the lowest and highest ends of the wavelength range, and update the scales of the cubes in between. Note however that any other choice is possible. Technical details and the expressions of the resulting optimal scale parameters are reported in Appendix A.

## 2.3. Matrix completion

As explained in Section 2.1, recovering the missing coefficients in the full matrix X displayed in Figure 3 can be interpreted as a *matrix completion* problem. Matrix completion aims at predicting missing values given observed values and the assumption of a latent matrix structure (see, e.g., Chi et al. (2019)). In our case, and as in many other settings, it makes sense to assume that X has a *low-rank* structure, an assumption which has already underpinned several techniques of infrared spectroscopic data processing such as source separation (Rapacioli et al. 2005; Berné et al. 2007) or data fusion (Berné et al. 2010; Guilloteau et al. 2020a,b). It consists in imposing that the columns of X can be explained by linear combinations of elementary spectra embedded in noise. To further exploit the non-negativity of the observed spectra, we propose to perform this matrix completion by resorting to the technique of non-negative matrix factorization (NMF) (Lee & Seung 1999; Smaragdis et al. 2014). In other words, the *missing entries* are approximated by $x_{ij} \approx [\mathbf{WH}]_{ij}$ where $\mathbf{W}$ and $\mathbf{H}$ are non-negative matrices of size $\Lambda \times K$ and $K \times P$ estimated from the low-rank approximation of the *observed entries*. More precisely, denote by $\mathcal{O}$ the set of row and column indices of the observed values in X, i.e., the non-white part of X in Figure 3 (or equivalently the supports of the sub-matrices $X^{c,b}$). We want to estimate the two matrices $\mathbf{W}$ and $\mathbf{H}$ by solving the following optimization problem

$$ \min_{\mathbf{W},\mathbf{H} \ge 0} \sum_{(i,j) \in \mathcal{O}} d(x_{ij} | [\mathbf{WH}]_{ij}) \quad (1) $$

where $d(u|v)$ denotes a discrepancy measure between the non-negative numbers u and v. In this study we leverage the work of Févotte & Idier (2011) that proposed easy-to-implement and efficient multiplicative updates for estimating $\mathbf{W}$ and $\mathbf{H}$ when $d(\cdot|\cdot)$ is chosen as a so-called $\beta$-divergence. The $\beta$-divergence is a continuous family of measures of fit governed by a single shape parameter $\beta \in \mathbb{R}$, that takes well-known divergences as special cases, namely the generalized Kullback-Leibler and Itakura-Saito divergences ($\beta = 1, \beta = 0$, respectively) and the square Euclidean distance ($\beta = 2$). Given a couple of matrices $\mathbf{W}$ and $\mathbf{H}$ that solve the minimization problem (1), the missing coefficients of X can be reconstructed as $\hat{x}_{ij} = [\mathbf{WH}]_{ij}$ for $(i, j) \notin \mathcal{O}$.

Like most NMF techniques, the algorithm tailored by Févotte & Idier (2011) is a descent algorithm that relies on alternating updates of $\mathbf{W}$ and $\mathbf{H}$. Because of the bilinearity induced by the product $\mathbf{WH}$, the objective function to minimize in Equation (1) is non-convex and the algorithm is likely to produce lo-
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Colorized image NIRCam of the Orion Bar of the PDRs4All ERS program (Berné et al. (2024)). Filters F140M and F210M are in blue; F277W, F300M, F323N and F335M in green; F405N in orange; and F444W, F480M and F470N in red. The pattern in white is the MIRI-MRS mosaic footprint. The red box corresponds to the pointing used in this article.</figcaption></figure>

cal solutions that depend on the chosen starting point. As such, the initialization of the descent algorithm is an issue and we will investigate several strategies in the experimental section.

# 3. Application to MIRI-MRS data of PDRs4All ERS program

## 3.1. Data and preprocessing

**Data.** The observations used in this article are part of the Early Release Science (ERS) program PDRs4All (Berné et al. 2022) obtained in January 2023. This program observed the Orion Bar with MIRI-MRS in a 1 × 9 pointing mosaic in the four channels. The mosaic is depicted in white in Figure 4. Data was acquired with 47 groups per integration, 1 integration and 4 dithers using the FASTR1 readout. Data was presented by Chown et al. (2024) and Putte et al. (2024). The experiments reported in this paper have used a single mosaic tile containing the protoplanetary disk d203-506 (highlighted in red in Figure 4) which was also previously studied by Berné et al. (2023) and Zannese et al. (2024). Figure 5 presents the MIRI-MRS original data cubes with the initial field of view in each channel.

**Data reduction pipeline.** To illustrate the relevance of the proposed stitching method on a set of data representative of the current quality of JWST data products, the MIRI MRS data was reprocessed starting from the ‘uncal’ files downloaded from MAST. For this purpose we used the scripts provided by the PDRs4All data reduction team and available from a public repository<fnref n="1" />, in combination with the JWST pipeline version 1.17.1 and CRDS context version ‘jwst_1322.pmap’. First, the standard pipeline stages were applied with the default settings for stage 1 (detector1), and the following settings for stages 2 and 3. In stage 2, the image-to-image background subtraction method was applied (‘bkg_subtract’), using the ‘rate’ files obtained by applying stage 1 to the dedicated background ob-

servations for this program. The ‘residual_fringe’ step was enabled, which reduced the residuals of fringes due to the reflections within the detector (Argyriou et al. 2020). In stage 3, we enabled the ‘outlier_detection’ step, and disabled the ‘cube_build’ step, so that the final products were the fully processed 2D IFU images (‘crf’ files). The default ‘cube_build’ step produced cubes whose pixel scale and field of view were different for each of the four MIRI MRS channels. We decided to impose a common spatial grid across all wavelengths for the data to be stitched. This choice was motivated by two main reasons: i) focusing on the exact intersection of fields of view with distinct spatial pixel resolutions would have lead to a significant loss of information; ii) because of lower spatial resolution, the data at longer wavelengths carries less information (in terms of numbers of pixels), which would have biased the results. This resampling could have been performed after data reduction during the preprocessing step (see below) by a naive spatial interpolation. Instead, this resampling was performed during the data reduction to benefit from the functionalities offered by the pipeline which performs this interpolation more efficiently, using both spatial and spectral information. Therefore the ‘cube_build’ step was performed with custom settings, starting from the stage 3 ‘crf’ files. The full specification of these settings was as follows. The position angle was set to ‘cubepa=250.4’, which aligned with the orientation of the 9x1 mosaic. The pixel scale and field of view were set to those of channel 1 using ‘scale_xy=0.13’, ‘nspax_x=33’ and ‘nspax_y=39’, and the center coordinates were set to ‘ra_center=83.834782’, ‘dec_center=-5.418207’. For channel 1, the resulting cube had nearly identical properties to the one built with the default coordinate system. For the longer wavelength channels, the pixel scale of these custom cubes substantially oversampled the instrumental resolution; this was the expected type of input for the method presented in this work.

**Preprocessing.** The *Haute Couture* stitching method proposed in Section 2 assumes that the spectral grids between two adjacent cubes are the same for the shared wavelength range, which is not the case for the data provided by the reduction pipeline. As such, the data was spectrally resampled on a common wavelength grid with a consistent step using a linear 1D interpolation.

Additionally, the data contains saturated frames, spikes and so-called bad pixels (i.e., outliers) that needed to be cleaned. Instead of merely discarding the saturated frames, the saturated frames were replaced by corrected frames resulting from a spectral linear interpolation to avoid gaps with no data for some spectral bands. This step was shown to be of major importance: if the data cubes to be stitched contained one saturated frame, the stitched result might be significantly distorted, e.g., containing unexpected spatial patterns and significant spectral bias. Regarding the remaining spikes and bad pixels in each frame, they were automatically detected and corrected. A significant number of those anomalies exhibited structures with high intensity. To identify them, we implemented a robust anomaly detector following the strategy proposed by Anderson & Gordon (2011). In each frame $X_{ij}^{c,b}$ associated with a given channel c and sub-channel b, the median pixel value $\bar{x}_{ij}^{c,b} = \text{median}(X_{ij}^{c,b})$ was first computed as a robust counterpart of the mean pixel values. Then each pixel value significantly far from this robust mean was identified as an anomaly. More precisely, a pixel value $x_{ij}^{c,b}$ such that $|x_{ij}^{c,b} - \bar{x}_{ij}^{c,b}| > t \times \text{MAD}(X_{ij}^{c,b})$ was flagged as an anomaly, where $\text{MAD}(X_{ij}^{c,b}) = \text{median}(x_{ij}^{c,b} - \bar{x}_{ij}^{c,b})$ computes the median absolute

<footnote n="1">¹ https://github.com/PDRs4All/PDRs4All</footnote>
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** The first row presents the average spectra over the entire field of view of each channel (shown in the second row), with a different color (blue, orange, green) for each sub-channel (short, medium, long). The second row presents the MIRI-MRS images at different wavelengths (the images are aligned). It can be observed that the field of view increases with the channel, but that the spatial resolution deteriorates.</figcaption></figure>

deviation and *t* is a cutoff parameter that needs to be adjusted for each channel and sub-channel. The pixel values detected as anomalies were finally replaced by the average of their four spatially nearest neighbors.

Finally, we created a mask containing only the most informative pixels. A pixel value is flagged as not carrying any information if its value is *NaN*. For each cube, the mask corresponds to the spatial pixels that have more than 99% spectral information (i.e. that are not *NaN*). In the mask, the star at the top of our field of view is removed because it is saturated in most of the wavelength range and the spectra are not exploitable.

## 3.2. Choice of the parameters and initialization

As discussed in Section 2, *Haute Couture* requires a few hyper-parameters to be set by the end-user. A first hyper-parameter is the value of *β*, i.e., the shape parameter of the divergence used as a measure of fit in NMF. As thoroughly discussed by Févotte & Idier (2011), the *β*-divergence is a log-likelihood in disguise, and choosing *β* is similar to making a noise assumption. The range of practical values is generally *β* ∈ [0, 2], with *β* = 0 corresponding to multiplicative Gamma noise and *β* = 2 corresponding to additive Gaussian noise. Choosing *β* = 1 (generalized Kullback-Leibler divergence), corresponding to a Poisson-like continuous distribution, has been shown to offer an excellent trade-off in many applications. It was in particular advocated in the previous work by Berné et al. (2007) for source separation of similar mid-IR data using NMF.

A second hyper-parameter to set is *K*, i.e., the rank of the factorization (the common dimension of **W** and **H**). A suitable value needs to capture the latent structure in the data while preventing from over-fitting (i.e, modeling the noise rather than the inherent phenomena that explain the data). We followed a trial and error approach by assessing the quality of the reconstruction with a set of candidate values for *K*, namely *K* ∈ [2, 4, 6, 10, 15]. We considered only channel 2 and 3 in order to save computing time. For every couple of estimated matrices **W** and **H** we then calculated the correlation between the reconstruction and the MIRI-MRS spectra (spectra after preprocessing in order to compare) for each sub-channel and the average of the sub-channel. For all *K*, the average correlation was 0.935 ± 0.0006. The algorithm

was not very sensitive to the value of *K* in our example and, thus, we chose *K* = 6 as it maximized the correlation.

Besides choosing the hyper-parameters *β* and *K*, *Haute Couture* requires choosing an initialization strategy for **W** and **H** in the NMF algorithm. The performance of three methods were empirically assessed: i) random initialization, ii) using the outputs of the K-means clustering algorithm, iii) using the Maximum Angle Source Separation (MASS) technique proposed by Boulais et al. (2021). As for the adjusting the hyper-parameters *β* and *K*, we conducted an exhaustive set of experiments to assess the performance of each of these three methods in reaching the best fit, with *β* = 1, *K* = 6 and again using channel 2 and 3 to save computing time. Our experiments revealed that K-means and MASS return similar results and much better performance than random initialization. As such we favored MASS in the experiments as it was initially designed for infrared astronomical hyperspectral cubes.

<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Top: average spectrum of the MIRI-MRS data after removing the saturated frames, spikes and bad pixels (black curve), average of the stitched data recovered by *Haute Couture* (orange curve) and average of the result of the coarse stitching procedure (blue curve). Bottom: zooms in spectral ranges where overlap occurs between contiguous sub-channels (as already depicted in Figure 1).</figcaption></figure>
# 4. Results

## 4.1. Stitching results

Figure 6 (top) displays the average spectrum of the MIRI-MRS data over the entire field of view after removing the saturated frames, the spikes and the bad pixels as explained in Section 3.1. It depicts the average spectrum of the stitched data $\hat{S}$ recovered by *Haute Couture* (corresponding to the values in channel 1-long of the submatrix highlighted in the blue shaded area of Figure 3). This result is compared to the product given by the coarse stitching procedure explained in Appendix B. The figure also displays three zooms in the spectral ranges where overlap occur between contiguous sub-channels (bottom). By design, the *Haute Couture* stitched data average spectrum is continuous, i.e., the spectral gaps exhibited in the original MIRI-MRS data have been removed. It is also worth noting that the fit between the original and *Haute Couture* reconstructed data is rather tight over the full spectral range.

Figure 7 displays the preprocessed MIRI-MRS data before stitching and compares the final products resulting from the coarse stitching procedure and from the proposed *Haute Couture* method. This figure also depicts the *Haute Couture* product after undergoing a convolution with a Gaussian filter whose standard deviation $\sigma = \frac{1.22\lambda}{2\sqrt{2\ln 2D}}$ has been specifically designed to mimic the effect of the MIRI PSF, where $\lambda$ is the considered wavelength and D is the mirror diameter. The resulting convolved product can then be also fairly compared to the original MIRI-MRS data as a validation step.

In the first row, corresponding to data and results at 5.61µm (channel 1), the *Haute Couture* image appears as significantly denoised and exhibits the same spatial resolution as the original data. This was expected since, as motivated in Section 2.1, the *Haute Couture* stitching procedure targets the spatial and spectral resolutions of channel 1-long (blue shaded area in Figure 3). Moreover, the observed denoising effect can be attributed to the low-rank factor model underlying *Haute Couture*. By implicitly representing the spectral data into a K-dimensional subspace, the proposed stitching approach preserves the essential physical information while effectively filtering out most of the noise. The second row associated with results at 9.29µm (channel 2) shows a lightly denoised and more contrasted *Haute Couture* image compared to MIRI-MRS image. At 12.27µm (channel 3), the *Haute Couture* image is sharper and more contrasted than the MIRI-MRS data. Interestingly, the spatial structure of the 203-506 protoplanetary disk, known to be present in the considered astrophysical scene Berné et al. (2024), is well reconstructed. The MIRI-MRS image at 20.42µm (channel 4) does not contain any particular spatial structure and the 203-506 protoplanetary disk is no longer visible. In the *Haute Couture* image, the spatial structures are well recovered.

Conversely, in all channels, the product resulting from the coarse stitching procedure is affected by a worse contrast and a clear under-estimation of the spectral flux density already highlighted in Figure 6. It is also noting that the loss in spatial resolution is significant. This was expected since the coarse stitching procedure targets the poorest resolution offered by the data cubes to be stitched.

The most straightforward way to validate an algorithm is to test it against ground truth data, as it is often possible with Earth observation data. However, in the case of astronomical data, obtaining ground truth is impossible. A widely used alternative consists in validating the algorithm using simulated data. This strategy has been followed for instance by Guilloteau et al.

(2020b) to validate data fusion algorithms by generating simulated NIRCam and NIRSpec data. Simulating MIRI-MRS data is, however, more challenging (Hadj-Youcef et al. 2017). The *Pandeia* tool (Pontoppidan et al. 2016) can produce simulated MIRI-MRS data, but only allows to simulate simple scenes. The MRISim tool (Klaassen et al. 2021) allows for more complex scenes, however to our knowledge, the public version of this software has not been updated after the launch of JWST, and therefore the tool may not provide data that is directly comparable to real observations. In this context, as a sanity check, we compare in Figure 7 the results of *Haute Couture* degraded data after undergoing the Gaussian smoothing (column 4) and those of the original MIRI-MRS data (column 1). It can be seen there is a good qualitative agreement in terms of the spatial textures between images of column 1 and 4, especially in channels 1 and 2. At longer wavelengths, the images differ notably in the right part of the scene. This can be explained by the masking procedure implemented during the preprocessing, which removes pixels associated with the star located in this area (see Section 3.1).

## 4.2. Error propagation

The MIRI-MRS cubes provided by the pipeline as detailed in Section 3.1 and stitched in Section 4.1 are provided with so-called error cubes (specified by the extension ERR). More precisely, each data cube $X^{c,b}$ composing the matrix $X$ depicted in Figure 3 is granted with an error matrix denoted $\Delta X^{c,b}$. To quantify the relative amount of error contained in the original MIRI-MRS data to be stitched, we first evaluate the signal-to-error ratio (SER) defined as

$$ \text{SER} = 10 \log \frac{\sum_{c,b} \|X^{c,b}\|_2^2}{\sum_{c,b} \|\Delta X^{c,b}\|_2^2} = 53.09\,\text{dB.} \quad (2) $$

We then propose to assess the propagation of this error when performing stitching. To do so, we introduce the set of data matrices $X_+^{c,b} = X^{c,b} + \Delta X^{c,b}$ and $X_-^{c,b} = X^{c,b} - \Delta X^{c,b}$ corresponding to the original data positively and negatively perturbed by the error term in each channel c and sub-channel b. These submatrices are rearranged to form the incomplete matrices $X_+$ and $X_-$ as perturbed counterparts of the original incomplete matrix $X$ introduced in Section 2 (see also Figure 3). The completion procedure detailed in Section 2.3 is applied on these two incomplete matrices $X_+$ and $X_-$, yielding the pair of stitched data $\hat{S}_+$ and $\hat{S}_-$ recovered by the proposed method *Haute Couture*. The propagation of the error is finally evaluated by quantifying how much the original amount of error translates into the perturbed data after undergoing stitching. Formally, we compute the SER associated with the positively and negatively perturbed data as

$$ \text{SER}_+ = 10 \log \frac{\|\hat{S}\|_2^2}{\|\Delta \hat{S}_+\|_2^2} \quad \text{and} \quad \text{SER}_- = 10 \log \frac{\|\hat{S}\|_2^2}{\|\Delta \hat{S}_-\|_2^2} \quad (3) $$

where $\Delta \hat{S}_\pm = \hat{S}_\pm - \hat{S}$ measures the discrepancy between the stitched spectra $\hat{S}$ resulting from the original data and the stitched spectra $\hat{S}_\pm$ resulting from the perturbed data, and is considered as the propagated error. We obtain $\text{SER}_+ = 33.43\,\text{dB}$ and $\text{SER}_- = 34.25\,\text{dB}$. The overall procedure is sketched in Figure 8. It appears that the SER deteriorates by about 20dB, which may be explained by the fact that the error terms $\Delta X^{c,b}$ do not likely follow the NMF model underlying the completion procedure detailed in Section 2.3. Yet, more importantly, it remains that the
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Comparison between the MIRI-MRS data (1st column), the result of the coarse stitching procedure (2nd column), the stitched data produced by *Haute Couture* (3rd column) and this final product convolved with a Gaussian filter (4th column) in each channel at a given wavelength. The same color bar (in MJy/sr) are used in each row and the data in columns 1, 2 and 4 share the same mask. The loss of spatial resolution is also clearly visible when performing a coarse stitching.</figcaption></figure>

amount of additional error induced by *Haute Couture* remains
of low magnitude with respect to both the original error and the
energy of the stitched data.

**5. Conclusion**

In this paper we presented a new stitching algorithm, coined as *Haute Couture*, which provides near-optimal assembly of MIRI-MRS spectral cubes. The source code for this method is available through [to be inserted upon publication]. *Haute Couture* allows the spatial information to be preserved throughout the
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Quantification of the error propagation. The amount of error SER before stitching is evaluated from the original data contained in **X** (light green box) and the error matrix **ΔX** (light yellow box). The original data **X** is processed by *Haute Couture* to produce the stitched spectra **S** (light blue box). The error matrix **ΔX** is combined with the original data **X** to produce the perturbed data **X**<sub>±</sub> (dark green box). The perturbed data **X**<sub>±</sub> is processed by *Haute Couture* to produce the stitched spectra **S**<sub>±</sub> (dark blue box), yielding the propagated error **ΔS**<sub>±</sub> (dark yellow box). The amounts of error SER<sub>±</sub> after stitching is evaluated from the stitched original spectra **S** and the propagated errors **ΔS**<sub>±</sub>.</figcaption></figure>

overall spectral range covered by the JWST. Applied on real ob-
servations, we show that *Haute Couture* is able to recover spatial
structures that are not preserved by a coarse stitching procedure.
The proposed approach relies on the formulation of the stitch-
ing task as a matrix completion problem, subsequently solved
by nonnegative matrix factorization. In principle, the versatil-
ity of this formulation makes the proposed approach applicable
to spectral cubes provided by other spatial instruments, such as
NIRSpec also embedded on the JWST, or by other ground-based
IFUs.

*Acknowledgements.* This work is partially supported by the Artificial and Natural Intelligence Toulouse Institute (ANITI), funded by the France 2030 program under the grant agreement ANR-23-IACL-0002. This work is funded by the Centre National d'Etudes Spatiales (CNES) through the APR program.

**References**

Anderson, R. E. & Gordon, K. D. 2011, Publications of the Astronomical Society of the Pacific, 123, 1237

Argyriou, I., Glasse, A., Law, D. R., et al. 2023, Astronomy & Astrophysics, 675, A111, publisher: EDP Sciences

Argyriou, I., Wells, M., Glasse, A., et al. 2020, A&A, 641, A150

Berné, O., Habart, E., Peeters, E., et al. 2024, Science, 383, 988

Berné, O., Helens, A., Pilleri, P., & Joblin, C. 2010, in 2010 2nd Workshop on Hyperspectral Image and Signal Processing: Evolution in Remote Sensing, IEEE, 1–4

Berné, O., Joblin, C., Deville, Y., et al. 2007, Astronomy & Astrophysics, 469, 575

Berné, O., Habart, E., Peeters, E., et al. 2022, Publ. Astron. Soc. Pac., 134, 054301

Berné, O., Martin-Drumel, M.-A., Schroetter, I., et al. 2023, Nature, 621, 56, publisher: Nature Publishing Group

Boulais, A., Berné, O., Faury, G., & Deville, Y. 2021, Astronomy & Astrophysics, 647, A105, publisher: EDP Sciences

Chi, Y., Lu, Y. M., & Chen, Y. 2019, IEEE Transactions on Signal Processing, 67, 5239

Chown, R., Sidhu, A., Peeters, E., et al. 2024, Astronomy & Astrophysics, 685, A75, publisher: EDP Sciences

Févotte, C. & Idier, J. 2011, Neural computation, 23, 2421

Gardner, J. P., Mather, J. C., Clampin, M., et al. 2006, Space Science Reviews, 123, 485

Guilloteau, C., Oberlin, T., Berné, O., & Dobigeon, N. 2020a, IEEE Transactions on Computational Imaging, 6, 1362

Guilloteau, C., Oberlin, T., Berné, O., Habart, É., & Dobigeon, N. 2020b, The Astronomical Journal, 160, 28

Hadj-Youcef, M. A., Orieux, F., Fraysse, A., & Abergel, A. 2017, in 2017 25th European Signal Processing Conference (EUSIPCO), IEEE, 503–507

Houck, J. R., Roellig, T. L., Van Cleve, J., et al. 2004, The Astrophysical Journal Supplement Series, 154, 18

Klaassen, P. D., Geers, V. C., Beard, S., et al. 2021, Monthly Notices of the Royal Astronomical Society, 500, 2813

Lee, D. D. & Seung, H. S. 1999, Nature, 401, 788

Perrin, M. D., Sivaramakrishnan, A., Lajoie, C.-P., et al. 2014, in Space telescopes and instrumentation 2014: optical, infrared, and millimeter wave, Vol. 9143, SPIE, 1174–1184

Pontoppidan, K. M., Pickering, T. E., Laidler, V. G., et al. 2016, in Observatory operations: strategies, processes, and systems VI, Vol. 9910, SPIE, 381–395

Putte, D. V. D., Meshaka, R., Trahin, B., et al. 2024, Astronomy & Astrophysics, 687, A86, publisher: EDP Sciences

Rapacioli, M., Joblin, C., & Boissel, P. 2005, Astronomy & Astrophysics, 429, 193

Rieke, G. H., Wright, G. S., Böker, T., et al. 2015, Publications of the Astronomical Society of the Pacific, 127, 584, publisher: IOP Publishing

Sandstrom, K. M., Bolatto, A. D., Draine, B., Bot, C., & Stanimirović, S. 2010, The Astrophysical Journal, 715, 701

Smaragdis, P., Févotte, C., Mysore, G., Mohammadiha, N., & Hoffman, M. 2014, IEEE Signal Processing Magazine, 31, 66

Wells, M., Pel, J.-W., Glasse, A., et al. 2015, Publications of the Astronomical Society of the Pacific, 127, 646, publisher: IOP Publishing

Werner, M. W., Roellig, T. L., Low, F., et al. 2004, The Astrophysical Journal Supplement Series, 154, 1

Wright, G. S., Rieke, G. H., Glasse, A., et al. 2023, Publications of the Astronomical Society of the Pacific, 135, 048003, publisher: The Astronomical Society of the Pacific

Zannese, M., Tabone, B., Habart, E., et al. 2024, Nature Astronomy, 8, 577, publisher: Nature Publishing Group
## Appendix A: Joint scale parameter estimation

This appendix details the approach introduced in Section 2.2 and implemented to mitigate intensity gaps observed in adjacent data cubes provided by MIRI-MRS. To ease the presentation of the method, we introduce the one-to-one mapping between the couple $(c, b) \in \{1, \dots, 4\} \times \{\text{s, m, 1}\}$ which indexes a given sub-channel $b$ of a given channel $c$ and the integer $i \in \{1, \dots, 12\}$ which indexes the sub-matrix composing $\mathbf{X}$ according to the particular order represented in Figure 3. Thus in what follows, the sub-matrix $\mathbf{X}^{c,b}$ is simply denoted as $\mathbf{X}^i$, leading to the alternate representation of the matrix $\mathbf{X}$ depicted in Figure A.1.

<figure><img src="image_11.png" /><figcaption>Fig. A.1. Matrix $\mathbf{X}$ as of Figure 3 with alternative indexing of the sub-matrices.</figcaption></figure>

The intensity gap correction aims at estimating a set of twelve scale parameters $\alpha = [\alpha^1, \dots, \alpha^{12}]^T$ to be applied to the twelve data matrices $(\mathbf{X}^1, \dots, \mathbf{X}^{12})$, respectively. By the very specific nature of the data acquisition, the wavelength range associated with the top rows (resp. bottom) of the matrix $\mathbf{X}^i$ emphasized as green (resp. blue) boxes in Figure A.1 overlaps with the bottom rows of the matrix $\mathbf{X}^{i-1}$ (resp. top of $\mathbf{X}^{i+1}$) emphasized as blue (resp. green) boxes in Figure A.1. It is worth noting that contiguous green and blue boxes frame portions of spectra that correspond to the same spatial locations at the same wavelengths. Thus the proposed approach consists in jointly adjusting the scale parameters by ensuring that the portions of pixel spectra located in the blue boxes match as much as possible the portions of the pixel spectra in the contiguous green boxes. Formally, the spatial average of the top (resp. bottom) rows of $\mathbf{X}^i$ where such overlap occurs is denoted by $\bar{\mu}^i$ (resp. $\underline{\mu}^i$), see Figure A.1. Maximizing the fit between $\alpha^i \underline{\mu}^i$ and $\alpha^{i+1} \bar{\mu}^{i-1}$, for $i = 1, \dots, 11$, boils down to solve the optimization problem

$$ \min_{\alpha} \sum_{i=1}^{11} \left\| \alpha^i \underline{\mu}^i - \alpha^{i+1} \bar{\mu}^{i+1} \right\|_2^2 \qquad (A.1) $$

with $\alpha = [\alpha^1, \dots, \alpha^{12}]^T$. The problem as formulated by Equation (A.1) is ill-posed since one optimal solution is $\alpha^* = [0, \dots, 0]^T$, i.e., setting all the sub-matrices to zero. To remove this trivial, inappropriate solution, one or more sub-channel are served as reference, with scale parameters fixed to an arbitrary value. In particular, in this work we assume that channels 1 and 12 have been properly calibrated and set $\alpha^1 = \alpha^{12} = 1$, but the proposed methodology may apply with any other arbitrary choices. Thus, we introduce the sub-vector of scale parameters $\tilde{\alpha} = [\alpha^2, \dots, \alpha^{11}]^T$ and reformulate the problem (A.1) as

$$ \min_{\tilde{\alpha}} \mathcal{J}(\tilde{\alpha}) \quad \text{with} \quad \mathcal{J}(\tilde{\alpha}) = \sum_{i=2}^{10} \left\| \alpha^i \underline{\mu}^i - \alpha^{i+1} \bar{\mu}^{i+1} \right\|_2^2. \qquad (A.2) $$

The criterion $\mathcal{C}(\tilde{\alpha})$ to be minimized is smooth and convex: Problem (A.2) has a unique solution that obeys

$$ \nabla_{\tilde{\alpha}} \mathcal{J}(\tilde{\alpha}) = 0 \qquad (A.3) $$

where $\nabla$ denotes the gradient operator. Trivial algebra shows that Equation (A.3) is equivalent to solving the linear problem

$$ \mathbf{M}\tilde{\alpha} = \mathbf{b}, \qquad (A.4) $$
where **M** and **b** are the square matrix and vector of sizes 10 given by

$$ \mathbf{M} = \begin{pmatrix} 1 & -\frac{\langle\mu^2, \vec{\mu}^3\rangle}{\|\mu^2\|^2 + \|\vec{\mu}^2\|^2} & 0 & \dots & 0 \\ -\frac{\langle\mu^2, \vec{\mu}^3\rangle}{\|\mu^3\|^2 + \|\vec{\mu}^3\|^2} & 1 & -\frac{\langle\mu^3, \vec{\mu}^4\rangle}{\|\mu^3\|^2 + \|\vec{\mu}^3\|^2} & \dots & \vdots \\ 0 & -\frac{\langle\mu^3, \vec{\mu}^4\rangle}{\|\mu^4\|^2 + \|\vec{\mu}^4\|^2} & 1 & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & \dots & 0 & \dots & -\frac{\langle\mu^{10}, \vec{\mu}^{11}\rangle}{\|\mu^{10}\|^2 + \|\vec{\mu}^{10}\|^2} & 1 \end{pmatrix} $$

$$ \text{and } \mathbf{b} = \begin{pmatrix} \frac{\langle\mu^1, \vec{\mu}^2\rangle}{\|\mu^2\|^2 + \|\vec{\mu}^2\|^2} \\ 0 \\ \vdots \\ 0 \\ \frac{\langle\mu^{11}, \vec{\mu}^{12}\rangle}{\|\mu^{11}\|^2 + \|\vec{\mu}^{11}\|^2} \end{pmatrix}. $$

The optimal scale parameters are thus given by $\tilde{\alpha}^* = \mathbf{M}^{-1}\mathbf{b}$.

## Appendix B: Coarse stitching

In this section, we present how we produce the coarse stitching presented as comparison in Section 4. We use the MIRI-MRS original data cubes presented in Section 3.1 and Figure 5.

**Field of view.** The first step is to have a common field of view corresponding to the intersection of each field of view. A exact common spatial grid is impossible due to the pixel size of each channel.

**Convolution with PSF.** Each channel and sub-channel is convolved with the PSF of the channel 4-long using the library `stpsf` (Perrin et al. 2014). We use the Fourier transform to perform the convolution. Then we calibrate the flux in order to have the same flux than before because there is a side effect when convoluting the data that loses flux. This effect is particularly visible at large wavelength because these are less pixels.

**Reprojection.** Each channel and sub-channel is reprojected in the spatial grid of the channel 4-long using the function `griddata` of the library `SciPy`.

**Stitching.** The stitching is performed individually for each spatial pixel. The channel 1-short is used as reference for the calibration. For each pixel the average of the fluxes are calculated on the spectral overlap between two successive sub-channel. Then the spectrum that is not the reference is scaled using the ratio of the average.