Restoration of Motion-Blurred, High-Resolution Mars Express SRC Images of Phobos

1. Introduction

Motion-induced blur (“motion smear”) degrades spacecraft imagery when the relative motion between the sensor and target is substantial, hindering detailed surface characterization. Notable cases include Mars Express’s Super-Resolution Channel (SRC) images of both Phobos [1] and Deimos [2], Viking Orbiter 2’s VIS captures of Deimos [3,4], New Horizons Long Range Reconnaissance Imager (LORRI) images of Arrokoth [5], and both Rosetta’s NAC and Philae’s ROLIS observations of comet 67P/Churyumov–Gerasimenko [6,7]. Although SRC can achieve sub-meter sampling and may offer the highest resolution images currently available for Phobos, residual “pepper” noise and severe smear have limited its scientific exploitation [1,8]. By modeling the blur as a linear PSF determined purely from spacecraft ephemeris and exposure timing, and applying Wiener filtering, we deliver a reproducible, one-shot deconvolution workflow that restores sub-meter detail across diverse imaging conditions.

2. Materials and Methods

2.1. Dataset Selection and Pre-Processing

As the spacecraft approaches a target body, the resulting images would become both more blurred and higher in resolution. Thus, using the European Space Agency (ESA) Planetary Science Archive (PSA) web interface, we first listed Mars Express High Resolution Stereo Camera (HRSC) [9] Super Resolution Channel (SRC) level 3 data products that were obtained with a slant distance of less than 300 km to Phobos. We browsed all publicly accessible datasets acquired within observation dates since 2003 to the latest release (as of June 2025). After that, we manually selected SRC image IDs whose boresight intersects visible (not dark) surfaces of Phobos by visual inspection. This finally achieved a list of fourteen highest-resolution (better than a few meters per pixel) SRC images which were observed during orbits 5851 (23 July 2008), orbit 7926 (10 March 2010), orbit 8974 (9 January 2011), and orbit 14776 (26 August 2015) (Table 1), three of which are consistent with those listed by Witasse, et al. [10].

Based on that, we acquired Planetary Data System (PDS) IMG data sets of HRSC radiometric Reduced Data Record (RDR) Extension series (EXT2, EXT3, EXT5, and EXT7) Version 4.0 [11,12,13,14] from the PDS Geosciences Node server (https://pds-geosciences.wustl.edu/mex). Also, we downloaded the MEx-related Spacecraft, Planet, Instrument, C-matrix, Events (SPICE) kernels [15] by running the Integrated Software for Imagers and Spectrometers (ISIS) version 8.3.0 developed by United States Geological Survey (USGS) [16], as well as an extra SPK kernel, MEX_STRUCT_V01.BSP, from the SPICE datasets in the ESA PSA server (MEX-E-M-SPICE-6-V2.0) [17].

We first converted their PDS IMG files to ISIS cubes, and attached MEx SPICE kernels to the cubes by running both "hrsc2isis" and "spiceinit" commands of ISIS 8.3.0 (Figure 1). To obtain better shift of surface intercept points to be required in the following image restoration, during "spiceinit", we add MEX_STRUCT_V01.BSP as an extra kernel and the latest, finest (~18-m mean spacing) shape model [8] as a DSK kernel that was converted from phobos_g_018m_spc_0000n00000_v002.obj (available at https://sbmt.jhuapl.edu/shared-files/) by using the Navigation and Ancillary Information Facility (NAIF) SPICE utility program MKDSK (https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/ug/mkdsk.html).

2.2. Denoising

The MEx SRC level 3 data products are already calibrated as the MEx team has officially reduced most of the pepper (dark) noise left in the level 2 products. However, a few pepper pixels still remain in the level 3 products. We regarded these residual pepper pixels as NULL pixels using the ISIS ”specpix”, which would cause problematic errors when subsequent processing in frequency-domain. Also, each PDS label contains the parameter "LINE_FIRST_PIXEL = 3", which is intended to indicate that the first two lines do not contain valid data. However, this assumption does not always hold because several of the individual cubes have data lines within the first three lines. Thus, we marked the pixels having DNs less than 200 within the first three lines as Low Representation Saturation (LRS), and completely removed them as well as dark pixels (75 NULL pixels per cube) by replacing the mean of their surrounding pixels by running ISIS “noisefilter” (third and fourth boxes in Figure 1). We repeated this noise filtering process until no LRS pixels remained in the resulting cube (twice for h5851_0003_sr3 and three times for h5851_0004_sr3). Moreover, any pixel whose raw DN value was ≤ –200 was considered anomalously dark and substituted with the average DN of its eight neighboring pixels. As exemplified in Figure 2, the NULL and dark pixels are removed from all SRC cubes.

2.3. PSF Estimation

The motion-induced blur can be modeled as a one-dimensional linear point spread function (PSF), parameterized by the spacecraft’s instantaneous ground-track velocity and the detector exposure time. These parameters are effectively determined by the spacecraft’s positions in body-fixed coordinate at the start and end of the exposure. Therefore, the resulting PSF’s length and its orientation (along the projected motion vector over Phobos) can be derived directly from spacecraft telemetry.

After collecting values regarding the IMG/cube labels "START_TIME" (UTC₁), "STOP_TIME" (UTC₂), and "EXPOSURE_DURATION" (exposure) with the USGS Abstraction Layer for Ephemerides (ALE) library [18,19,20], we first performed calculation of “high-accuracy” observation start and stop timings in forms of ephemeris time (ET) with microsecond precisions (ET’₁ and ET’₂,) based on:

$${ET\text{'}}_1=\left({ET}_1+{ET}_2\right)/2-exposure/2=t_1,$$

(1)

$${ET\text{'}}_2=\left({ET}_1+{ET}_2\right)/2+exposure/2=t_2,$$

(2)

where ET₁ and ET₂ are the ET expressions of UTC₁ and UTC₂, respectively.

Conversions from ET to UTC and vice versa were implemented by uses of SpiceyPy et2utc and utc2et, respectively. The resulted UTC′₁ and UTC′₂ (UTC expressions of the resulted ET’₁ and ET’₂) about the selected SRC images are listed in Table 2.

For PSF estimation, we started calculating the SRC boresight intersections with the Phobos surface DSK in the Phobos’ body-fixed reference frame “IAU_PHOBOS”. Using SpiceyPy [21] sincpt function with the aberration correction option of “CN+S” (both light time and stellar correction) gives two positional vectors at ET1’ simultaneously: one from observer (obs; the origin of SRC boresight, “MEX_HRSC” whose SPICE body code is -41200) to the surface intercept point (spt) (${\overrightarrow{r}}_{obs\to spt}\left(t_1\right)$ in Figure 3) and another from Phobos center (bod) to the same spt (${\overrightarrow{r}}_{bod\to spt}\left(t_1\right)$ in Figure 3). Taking their vector difference gives ${\overrightarrow{r}}_{obs\to bod}\left(t_1\right)$(as shown in Figure 3):

$${\overrightarrow{r}}_{obs\to bod}\left(t_1\right)={\overrightarrow{r}}_{obs\to spt}\left(t_1\right)-{\overrightarrow{r}}_{bod\to spt}\left(t_1\right)$$

(3)

Similarly, using SpiceyPy spkpos function with the “CN+S” aberration correction, we obtain a positional vector at ET2’ from obs to bod (${\overrightarrow{r}}_{obs\to bod}\left(t_2\right)$ in Figure 3) as well. Then, using it with the ${\overrightarrow{r}}_{obs\to bod}\left(t_1\right)$ resulted from equation (3), we can calculate the vector directing from obs at ET2’ to obs at ET1’:

$${\overrightarrow{r}}_{obs\to obs}\left(t_2\right)={\overrightarrow{r}}_{obs\to bod}\left(t_2\right)-{\overrightarrow{r}}_{obs\to bod}\left(t_1\right)$$

(4)

Using ${\overrightarrow{r}}_{obs\to obs}\left(t_2\right)$ derived from equation (4) and ${\overrightarrow{r}}_{obs\to spt}\left(t_1\right)$, we finally gain the positional vector at ET2’ from obs to spt (${\overrightarrow{r}}_{obs\to spt}\left(t_2\right)$ in Figure 3):

$${\overrightarrow{r}}_{obs\to spt}\left(t_2\right)={\overrightarrow{r}}_{obs\to obs}\left(t_2\right)+{\overrightarrow{r}}_{obs\to spt}\left(t_1\right)$$

(5)

The two resulting position vectors from obs to spt at ET1’ and ET2’, ${\overrightarrow{r}}_{obs\to spt}\left(t_1\right)$ and ${\overrightarrow{r}}_{obs\to spt}\left(t_2\right)$, are expressed in the Phobos-body fixed frame “IAU_PHOBOS” (SPICE FK ID: 10021). We thus translated these vectors into the ones in the instrument frame “MEX_HRSC_SRC” (SPICE FK ID: -41220) by multiplying the frame-to-frame (“IAU_PHOBOS” to “MEX_HRSC_SRC”) transformation matrices at ET1’ and ET2’ that were generated by SpiceyPy function pxform. Using two obs-to-spt vectors expressed in the “MEX_HRSC_SRC” frame, we derived image coordinates of the two corresponding points on the SRC image plane:

$$X_1=\left(focallength/{pixelsize}_x\right)\left(x_1/z_1\right),$$

(6)

$$X_2=\left(focallength/{pixelsize}_x\right)\left(x_2/z_2\right),$$

(7)

$$Y_1=\left(focallength/{pixelsize}_y\right)\left(x_1/z_1\right),$$

(8)

$$Y_2=\left(focallength/{pixelsize}_y\right)\left(y_2/z_2\right),$$

(9)

where (X₁, Y₁) and (X₂, Y₂) are respectively the coordinates of the projected point at ET1’ and ET2’, (x₁, y₁, z₁) and (x₂, y₂, z₂) are respectively the components of ${\overrightarrow{r}}_{obs\to spt}\left(t_1\right)$ and ${\overrightarrow{r}}_{obs\to spt}\left(t_2\right)$ expressed in the “MEX_HRSC_SRC” frame; focal length is 9.8476e+5 μm (= 984.76 mm [22]), pixel size_x is 9.0 μm/pixel, and pixel size_y is 9.0 μm/pixel [9], each of which is stored in SPICE IK MEX_HRSC_V09.TI.

Finally, we obtained the shift of the projected points during exposure on the SRC image plane by calculating their differences:

$${shift}_x=X_2-X_1$$

(10)

$${shift}_y=Y_2-Y_1,$$

(11)

The obtained PSF geometries are summarized in Table 3. By using longer exposures—values matching the SRC observation durations—our pipeline extracted extended PSFs, which is consistent with the initial expectation of this study. Most of the vertical components of the shifts are close to zero relative to horizontal components, so that the resulting PSF can be approximated to be a horizontally flat line (See Figure 4 for a typical case).

2.4. Wiener Deconvolution

There are two major deconvolution methods to restore degraded images: Wiener deconvolution [23] and Lucy–Richardson (or Richardson–Lucy) deconvolution [24,25]. Wiener deconvolution uses a linear, closed-form filter in the frequency domain, which has achieved restoration of planetary images, such as the Near Earth Asteroid Rendezvous (NEAR) Shoemaker spacecraft Multispectral Imager (MSI) images of asteroid Eros [26,27].

Wiener algorithm computes a single “optimal” filter $H_W\left(u,v\right)$ that, when applied to the blurred image’s Fourier transform $G\left(u,v\right)$, yields an estimate of the true image’s transform $F\left(u,v\right)$:

$$H_W\left(u,v\right)={\displaystyle \frac{H^{*}\left(u,v\right)}{{\left|H\left(u,v\right)\right|}^2+\frac{S_N\left(u,v\right)}{S_F\left(u,v\right)}}}={\displaystyle \frac{H^{*}\left(u,v\right)}{{\left|H\left(u,v\right)\right|}^2+\frac{1}{SNR}}}$$

(12)

$$\widehat{F}\left(u,v\right)=H_W\left(u,v\right)G\left(u,v\right)$$

(13)

where $\left(u,v\right)$are the frequency coordinates corresponding to the spatial coordinates $\left(x,y\right)$, $H\left(u,v\right)$ (and hence $H_W\left(u,v\right)$) is the transfer function (Fourier-domain PSF), $H^{*}$ is complex conjugate of $H$, $S_N\left(u,v\right)$and $S_F\left(u,v\right)$are the noise and signal power spectra, respectively (that is, $S_N\left(u,v\right)/S_F\left(u,v\right)$ is the noise-to-signal power ratio), Signal-to-Noise Ratio (SNR)?a dimensionless quantity?is given by $SNR=S_F\left(u,v\right)/S_N\left(u,v\right)$, and $\widehat{F}\left(u,v\right)$ is the estimate of the true image’s Fourier spectrum.

Lucy–Richardson deconvolution uses an iterative, non-linear algorithm in the spatial domain, and has been applied to planetary images, including MEx/SRC images of Phobos [28] and Mars [22] (but not so high resolution and terribly blurred as those addressed in this study). Lucy–Richardson deconvolution is based on maximum-likelihood estimation assuming Poisson (photon) noise. Starting from an initial guess $f^{\left(0\right)}$, it refines via

$$f^{\left(k+1\right)}\left(x,y\right)=f^{\left(k\right)}\left(x,y\right)\left[h\left(-x,-y\right)*{\displaystyle \frac{g\left(x,y\right)}{\left(f^{\left(k\right)}*h\right)\left(x,y\right)}}\right]$$

(14)

where “∗” denotes convolution and h is the PSF reversed. Each pass sharpens details and suppresses noise according to the Poisson likelihood.

Either would work for restoring SRC images; however, we selected Wiener deconvolution for the restoration in this study. This is because the PSF can be computed a priori from precise motion and exposure data, the closed-form Wiener filter offers a direct, non-iterative solution that (i) fully exploits the known kernel, (ii) minimizes noise amplification via the noise term, and (iii) avoids iteration-dependent artifacts and convergence issues associated with Lucy–Richardson (Table 4). Therefore, we used Wiener deconvolution, using with the linear PSFs aforementioned and an empirically determined SNR = 16 dB (i.e., ${10}^{\left(16/10\right)}\cong 39.8$ in linear scale) that works well regardless of a variety of SRC observation conditions. To mitigate serious ringing artifacts of output boundaries after deconvolution [29], we also used edge tapering algorithm (Gaussian-weighted taper at the edges) with the width of 81 pixels on both right- and left-sides (and no tapering at top and bottom sides) of each input image beforehand.

Although MATLAB offers built-in routines for Wiener deconvolution and edge tapering (e.g., deconvwnr and edgetaper), we chose to avoid reliance on proprietary software by reimplementing these methods in OpenCV. Specifically, we extended OpenCV’s deconvolution sample (https://github.com/opencv/opencv/blob/master/samples/python/deconvolution.py) to include edge tapering, yielding a fully open-source solution. Our Python implementation will be made available at https://github.com/ryodohemmi/public upon publication.

3. Results

We successfully restored motion-blurred SRC images of the four MEx orbits, significantly reducing blur and enhancing geological feature visibility. Our pipeline establishes a robust and objective framework for future restoration and analysis of motion-blurred planetary imagery.

3.1. Orbit 5851

Within the four selected MEx orbits (and most of all MEx orbits), the orbit 5851 has observed SRC images with the second highest resolutions (~0.86 m/pixel; Table 1) and the shortest exposure duration (~14 milliseconds; Table 2), which caused severe motion blurs with lengths of ~44 pixels (Table 3) in resulting images. However, unclear surface features, such as impact craters and portions of lineaments, have appeared unequivocally after restoration, as shown in Figure 5. Noisy artifacts in images are still visible, but ringing artifacts at edges are minimized in restored images. As for small-scale craters, a diversity of their depths, rims, floors, and degradational states is also identifiable, which is consistent with morphologies of sub-kilometer-scale craters on Phobos [30]. Based on the above inspection, we conclude that our restoration pipeline restored the SRC images in orbit 5851 to a satisfying level of quality.

3.2. Orbit 7926

The SRC image of MEx orbit 7926 has the lowest resolution and the longest slant distance (~2.7 m/pixel and ~292 km; Table 1) among all the SRC images selected in this study. With an exposure duration as short as that of orbit 5851, the shortest PSF only~11 pixels in length (Table 3) caused a slight motion smear to the original image, which is consistent with its appearance in Figure 6a. After our restoration process, we obtained the image that exhibits the well-defined craters and grooves on Phobos. Although the imaging distance to the surface would be proximal on the left side and distal to the right side of the figure, the restoration appears successful on both sides. A few noises are identified on the resulting image; however, the edge tapering effect has not reduced the original image’s quality seriously. As a result, we determined the SRC image of the orbit 7926 was successfully restored.

3.3. Orbit 8974

Seven images acquired in the MEx orbit 8974 have similar resolutions and slant distances to those of the orbit 5851 and a slightly longer exposure duration (~16.1 milliseconds; Table 2) with their PSF lengths of ~46 pixels (Table 3), which cover the southern hemisphere of anti-Mars surface of Phobos (Table 1). As shown in Figure 7, all SRC images captured a wide range of craters sizes and morphologies, and some of them observed lineaments (Figure 7a) and slopes (Figure 7g), which would be significant for analyses in surface irregularity (e.g., [31]). After the deconvolution, the resultant images show well-defined outlines of surface features much better than those of original images. Furthermore, we can identify several-pixel-scale, positive relief features (most likely boulders; Figure 7d, f) that were difficult to identify from original blurred images. Regardless of diverse landscapes, the deconvolved images show few noises and edge-tapered pixels at negligible levels. We express deblurring the orbit 8974 was quite achieved successfully.

3.1. Orbit 14776

Among the selected (and probably all available) MEx SRC images, those acquired in the MEx orbit 14776 are the worst motion-blurred, due to exposures of ~20 milliseconds (Table 2), the closest distance (~52 km), and the highest resolution (~0.47 meter/pixel ; Table 1), which is the hardest to recover by deconvolution. However, as shown in Figure 8, we performed restoration and acquired the images showing identifiable surface features, such as small impact craters and boulder-like positive reliefs. Even though there are noticeable noises and edge-tapering effects left in the resulting images, given the initial bad imaging conditions, we conclude our efforts made them as restored as possible, from a theoretical viewpoint.

4. Discussion

Our motion-blur restoration pipeline?based on an a priori, geometry-derived linear PSF and Wiener deconvolution?robustly recovers sub-meter Phobos surface details under diverse imaging conditions. By computing the PSF directly from spacecraft ephemeris and exposure data, we bypass the need for star-field references or manual kernel tuning, while the closed-form Wiener filter provides built-in noise regularization and avoids iteration-dependent artifacts. The extracted PSFs ranged from ~11 to ~119 pixels in length (reflecting ground-track speed and slant distance) with negligible vertical components, justifying their one-dimensional horizontal approximation.

In the high-resolution, short-exposure case (orbit 5851; 0.86 m/pixel, 14 milliseconds), deconvolution revealed crater morphologies, rim and floor textures, and linear features that were indiscernible in the raw images (Figure 5). Even at the coarsest resolution (orbit 7926; ~2.7 m/pixel, 14 milliseconds), grooves and crater edges sharpened with minimal ringing, demonstrating the method’s resilience to low SNR and large slant distances. For intermediate cases (orbit 8974; ~0.92 m/pixel, 16 milliseconds), deblurred images expose lineaments, slopes, and meter-scale boulders consistent with known Phobos morphologies. Even the most challenging dataset (orbit 14776; 0.47 m/pixel, 20 milliseconds) yielded clear views of small impact craters and positive reliefs, pushing the theoretical restoration limit under severe motion smear.

Compared to spatial-domain, iterative approaches like Lucy–Richardson?which have been applied to SRC and other planetary imagery [22]?our frequency-domain Wiener deconvolution delivers a one-shot solution with explicit control over the noise–sharpness trade-off and guaranteed convergence. Because the PSF is derived from precise SPICE kernels and shape models, our framework is completely reproducible and readily adaptable to other spacecraft instruments suffering motion-induced blur.

Scientifically, recovering sub-meter features enables improved crater counting, regolith maturity analysis, and the generation of high-fidelity digital elevation models on Phobos. This approach can be extended to data from missions with similar blur challenges, and will be especially valuable for the forthcoming MMX mission [32,33]. Integrating our deblurring pipeline into photogrammetric and stereophotoclinometry workflows promises significant gains in topographic precision and geological interpretation.

Nevertheless, our treatment assumes a spatially invariant, linear PSF dominated by along-track motion. In reality, additional blur may arise from cross-track jitter, optical aberrations, or attitude fluctuations. We also adopted a fixed SNR of 16 dB; future work should explore adaptive noise estimation to further suppress residual artifacts. Extending the model to 2D, spatially varying PSFs, and hybridizing Wiener with iterative methods may yield further improvements. Finally, automating noise-term selection based on image statistics will enhance applicability across a broader range of planetary datasets.

5. Conclusions

We have demonstrated an end-to-end pipeline for restoring motion-blurred, high-resolution SRC images of Phobos by automatically deriving a linear PSF from SPICE-based geometry and applying Wiener deconvolution with built-in noise regularization. Across fourteen images spanning four orbits, our method consistently recovers sub-meter surface features and dramatically enhances scientific usability (e.g., elucidating the origins and evolutionary history of small-scale structures and their association with spectral units [34,35]). Notably, this includes the restoration of the highest-resolution SRC images of Phobos acquired to date, which had previously been severely blurred. This approach is fully reproducible, requires no star-field calibration, and is readily extendable to other planetary missions facing motion-induced blur.