Comparison of MTF measurements using edge method: towards reference data set

A sensor’s spatial resolution has traditionally been a difficult concept to define, but all would agree that it is inextricably linked to the Ground Sampling Distance (GSD) and Instantaneous Field of View (IFOV) of an imaging sensor system. As a measure of the geospatial quality of imagery, the Modulation Transfer Function (MTF) of the system is often used along with the signal-to-noise ratio (SNR). However, their calculation is not fully standardized. Further, consistent measurements and comparisons are often hard to obtain. Therefore, in the Infrared and Visible Optical Sensors (IVOS) subgroup of the Working Group on Calibration Validation (WGCV) of the Committee for Earth Observation Satellites (CEOS), a team from various countries and professional entities who are involved in MTF measurement was established to address the issue of on-orbit MTF measurements and comparisons. As a first step, a blind comparison of MTF measurements based on the slanted edge approach has been undertaken. A set of both artificial and actual satellite edge images was developed and a first comparison of processing results was generated. In all, seven organizations contributed to the experiment and several significant results were generated in 2016. No single participant produced the best results for all test images as measured by either the closest to the mean result, or closest to the truth for the synthetic test images. In addition, close estimates of the MTF value at Nyquist did not ensure the accuracy of other MTF values at other spatial frequencies. Some algorithm results showed that the accuracy of their estimates depended upon the type of MTF curve that was being analyzed. After the initial analysis, participants were allowed to modify their methodology and reprocess the test images since, in several cases, the results contained errors. Results from the second iteration, in 2017, verified that the anomalies in the experiment’s first iteration were due to errors in either coding or methodology, or both. One organization implemented a third trial to fix software errors. This emphasizes the importance of fully understanding both methodology and implementation, in order to ensure accurate and repeatable results. To extend this comparison study, a reference data set, which is composed of edge images and corresponding MTF curves, will be built. A broader audience will be able to access the edge images through the CEOS CalVal Portal (http://calvalportal.ceos.org/. This paper, which is associated with the reference data set, can serve as a new tool to either implement or check, or both, the MTF measurement that relies on the slanted edge method. ©2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement


Introduction
The geo-spatial quality of a sensor and its imagery often revolves, at least in part, around the concept of the spatial resolution of a sensor which is often reduced to the Ground Sampling Distance (GSD) associated with the Instantaneous Field Of View (IFOV) defined by the pixel size. However, spatial resolution, and hence geo-spatial quality, is more complex than this. Most agree that the effective spatial resolution is due to three (or four) features of the sensor: the IFOV (and the GSD if different), the Modulation Transfer Function (MTF) and the signal to noise ratio (SNR). The MTF is often used as a measure of image sharpness [1,2]. This important parameter for image quality has to be checked on orbit in order to be sure that launch vibrations, transition from air to vacuum, or thermal state have not degraded the sharpness of the images. In some cases, it can lead to a refocusing decision.
This paper deals with one of the methods used for on-orbit MTF assessment, called the edge method, the knife-edge method, or the slanted-edge method. This method is widely used for laboratory measurements and may be implemented in various manners. For on-orbit MTF assessment, it requires a slanted edge as explained in section 2. It has been used for numerous space sensors such as Landsat TM [3], MOS-1 MESSR [4], IKONOS [5], SPOT5 [6], and more recently Sentinel2 MSI [7].
In the framework of the Infrared and Visible Optical Sensors (IVOS) subgroup of the Working Group on Calibration Validation (WGCV) of the Committee for Earth Observation Satellites (CEOS), a team of people from various countries and professional entities, who are involved in MTF measurement, has been created to address a variety of issues regarding the geo-spatial quality of optical satellite imagery. One of the first efforts of this group has been to compare processing methodologies for the edge method of MTF estimation. For this comparison experiment, the team was composed of Frans van den Bergh from CSIR, Renaud Fraisse from Airbus DS, Dennis Helder from SDSU, Dong Han Lee from KARI, Amy Newbury and Robert Kudola from Digital Globe, Sébastien Saunier from Telespazio and Françoise Viallefont-Robinet from ONERA This paper presents the method and its various implementations followed by the comparison experiment. The first results, obtained with a blind test approach, were analyzed. This exercise was an opportunity to correct or improve the software of each participant. Thus, a second run was performed by most of the participants in order to improve the results, leading to a second comparison. For two test cases, a third and final run was performed by one of the outliers. All comparisons are presented and commented on.

Theory
Considering the sensor as a linear system without spatial variation of its response (shift invariant), the relation between the radiances (or top of atmosphere reflectances) of the landscape and the image is simply: where i(x,y) stands for the image, l(x,y) stands for the landscape, h(x,y) is the Point Spread Function of the sensor, ⊗ is the convolution integral. For sensors using a CCD in the image plane, the system is no longer strictly shift invariant; but it can nevertheless be described using the usual theory for regions of the imaging array without loss of generality. The sampling done by the CCD can be written as a multiplication by a Dirac comb.
Completing Eq. (1), it becomes: [ ] ( , ) ( , ) ( , ) ( / , / ) x y i x y l x y h x y comb x p y p = ⊗ ⋅ (2) where p x is the size of the IFOV and p y the GSD in the case of a pushbroom imaging system. A classical way to deal with a convolution product is to apply a Fourier Transform, which leads to: ( , ) ( , ) ( , ) ( / , / ) x y x y x y x s x y s y where I(f x ,f y ) stands for the Fourier Transform of the image, L(f x ,f y ) stands for the Fourier Transform of the landscape, H(f x ,f y ) is the Transfer Function of the sensor, f sx = 1/p x is the sampling frequency for the f x axis, f sy = 1/p y is the sampling frequency for the f y axis. The sensor behavior is known to be similar to a low-pass filter without phase shift [8]. This is why the Optical Transfer Function is usually reduced to the Modulation Transfer Function defined as the modulus of the Optical Transfer Function normalized by the zero frequency component.
For the edge method, the landscape is close to a Heaviside function: hea(x) being the Heaviside function centered on x = 0, unit(x) = 1 for all values of x. In this case, Eq. (2) becomes: (5) can be rewritten as follows: The convolution by the comb produces aliasing. One way to overcome this problem is to use an edge with a slight inclination relative to the row or column direction [9]. This is used to build an oversampled 1-D edge image as illustrated in Fig. 1. So, the 1-D edge corresponds to: where z is measured in the direction perpendicular to the edge, which nearly coincides with the x axis when the edge is vertically oriented as modelled here.

Impleme
Implementatio modeling, Ed dentifiable ste on of an MTF eps: edge using the sampled ESF. The implementation of each of these steps can vary significantly, as discussed in greater detail below. Some common traits of the participant implementations will be discussed in this section. Without loss of generality it will be assumed that the Region Of Interest (ROI) containing the edge transition will be processed on a row-by-row basis with the slanted edge oriented nearly vertically. Let R denote the ROI such that i(x, y) represents the intensity of the image at coordinates (x, y) for all (x, y) ∊ R.

Edge modeling
The construction of an oversampled ESF, as illustrated in Fig. 1, requires that the location of the edge in each image row is known with sub-pixel accuracy. This can be as simple as calculating the centroid of the discrete derivative of the intensity of each row, as suggested in the ISO 12233 standard [10, Appendix D], but such a method will be sensitive to image noise and target non-uniformity. A more robust method is to fit a parametric-or spline function to each row, using the inflection point of the fitted function as the sub-pixel edge location estimate.
Once the location of the edge has been estimated in each row, a linear function is typically fitted across all the per-row results to obtain a more accurate model describing the sub-pixel location of the edge, as illustrated in section 3.1 in Fig. 3. If the physical target edge is curved, the linear edge model can be replaced with a low-order polynomial to accommodate the curvature [11].

ESF construction
A simplified example of the construction of an ESF is illustrated in Fig. 1, where the oversampled ESF is obtained by interleaving the intensity values of each row of the ROI. Constructing the ESF involves projecting the 2-D image intensity values i(x,y) onto a 1-D representation i(z). The magnitude of z represents the shortest distance from a pixel at coordinates (x,y) to the slanted edge. The oversampled ESF can be constructed using an extension of the method described in the ISO 12233 standard [10, Appendix D], or using the alternative method described by Kohm [12].
In the ISO 12233-based approach the ROI is a rectangle aligned with the rows and columns of the image, as shown in Fig. 1. For each pixel the value z is calculated as: where e(y) denotes the location of the edge in row y, as predicted by the edge model, and θ denotes the relative edge angle. The cosine factor transforms a distance measured along a row into a distance measured perpendicularly to the edge. With Kohm's method, the ROI is a rectangle that is aligned with the edge itself, as shown in Fig. 2. A unit vector perpendicular to the edge, n  , is used to calculate z such that: where (x 0 ,y 0 ) are the coordinates of an arbitrary point on the edge. The ESF c followed by s can be piecew another suitab in the ISO 12 additional low uniformity.

ONERA
The ONERA implementation follows the spectral approach [9], which is less widely known than the derivative method [4,5,14].
In this case, the finite number of samples has to be taken into account. This can be done as follows: noting w(z) as the window corresponding to the finite interval.
In the Fourier domain, the relation becomes: The window has to be chosen so that Hea(f z )⊗W (f z ) ≠ 0 for all frequencies and not far from Hea(f z ). After removal of the background b, the following ratio gives the transfer function: ONERA tool follows the 3 steps general implementation. For edge modeling, or in other words the 2-D to 1-D transformation, each row is interpolated using a spline function and the inflection point is computed. A straight line is fitted on the set of inflection points and may be used (depending on the user choice) to replace the positions found as shown in Fig. 3. The inverse of the slope of the straight line provides the oversampling rate for building the ESF. At this stage, the 2-D edge image is split into a list of rows where each ESF location is indicated by the position of the edge e(y).
The second step aims at mixing the rows to construct the oversampled ESF. For the inclination, in the ideal case of Fig. 1, the sequence spans a whole number (4 in this case) of rows. However, it is not possible to have a whole number of rows corresponding to the sequence for any edge. Thus, the user chooses the oversampling rate Nr, not too far from Ns (Ns being the whole number closest to the actual number of rows of the sequence) and then each row of the sequence is properly phased in the Nr grid. Some cases of Nr ≠ Ns, may lead to incomplete sampling.
A user can choose between linear interpolation or model fitting options to construct an ESF with regular sampling. In the model fitting case, all the edge samples are replaced by the samples deduced from a parametric transfer function model as explained in [15]. The model fitting implementation eliminates noise due to non-uniformity of the dark or light areas of the ROI, which is an added advantage even when the sampling is regular.
Once the regularly sampled ESF is built, the parameters a and b, as well as the position z = z 0 of the Heaviside function, are assessed. Several approaches are available to the user for selection of a, b, and z 0 : • the first value of the edge for a and the last value for b, • the values at a distance of -d (relative to z 0 ) for a and + d for b, • the value corresponding to (a + b)/2 for z 0 , • the value corresponding to the inflection point for z 0 Once the parameters a, b and z 0 are found, the distribution a.Hea(z) + b is built and drawn over the ESF. The background b is then removed.
ESF may be artificially extended, as shown Fig. 4 in order to reach the desired sampling rate for the MTF curve and to limit windowing effects.
A Hann window is applied to the ESF and to the distribution. An FFT is applied to the windowed ESF and to the windowed distribution. The ratio of the modulus of each spectrum is computed and normalized to obtain the MTF according to Eq. (12). The f z frequencies are converted to f x frequencies. Once the MTF is computed, there is a possibility to fit the model described in [15] to the values obtained and to resample the MTF in order to obtain the desired frequency step.

SDSU
Based on the theory presented in Section 2, the South Dakota State University implementation follows a series of steps designed to be largely target and model independent so the algorithm will work with a variety of edge target types. The steps are described as follows.
First, it is important that, to the degree possible, an edge target is used that has a proper orientation with respect to the sampling grid of the satellite sensor. Some combinations of edge angles and sampling grid orientation will result in data that are not reasonably uniformly distributed as shown in Fig. 1. When this happens, the result is that the steep part of the edge is represented by data points that are clustered together. This leads to very poor LSF/MTF estimates. Extensive modeling has indicated that relative angles of 6-8 degrees are optimal. This has the added advantage of orienting the edge nearly orthogonal to the sampling grid, and minimizes the correction necessary to obtain LSF/MTF estimates in the typical alongtrack and cross-track directions in which sensor specifications are often given.
A second critical aspect of target development is that the length of the edge should span a sufficient number of rows (or columns) in the image so that the oversampling process produces enough samples for edge reconstruction to be accurate. Empirical analysis has indicated that a minimum of 20 cross sections of the edge should be obtained for accurate results.
Signal-to-noise ratio (SNR) has a significant bearing on PSF/MTF estimation accuracy. For purposes of LSF/MTF estimation, SNR can be defined as the ratio of the edge height to the average of the standard deviations of the region on either side of the edge. Modeling has indicated that SNR > 50 produces accurate and consistent results.
The first step to developing an oversampled ESF is to estimate the edge location from each slice of data across the edge. A simple, but accurate, approach is to fit a Fermi function to the data of the form of Eq. (13): where x represents the pixel locations for row y, d is the bias level, b is the magnitude of the bright side of the edge, e(y) is the edge location and s represents the steepness of the edge. This approach does not model ESF which have ringing in them, but will still fit the steep part of the ESF well and give good estimates of the precise edge location which is the goal of this step. To determine the parameter values, a common optimization algorithm, such as the Levenberg-Marquardt algorithm, is employed. The output of this step is the parameter, e(y), which provides, for each row, subpixel estimate of the true edge location using the integerbased grid of the edge image input data. An oversampled but irregularly spaced ESF is constructed using Eq. (8). Only the data within a distance of five pixels from the edge are retained, based on typical LSF width of two pixels.
The truncated, oversampled ESF is simultaneously filtered (to reduce high frequency noise) and resampled (to obtain uniformly spaced samples) using a non-linear modified Savitzky-Golay filter. In this approach, a window of data is selected (typically set to two pixels in length) and a low-order polynomial is fitted to the data using a linear least squares approach. Through extensive modeling, it was determined that a fourth order polynomial is optimal. The output of the filter is the value of the polynomial at the center of the window. The window is then shifted by the amount of the desired output sampling interval, and the process is repeated. It is recommended to oversample by a factor of 10 or more. This step produces as an output a uniformly oversampled ESF that has been smoothed for high frequency noise but, at the same time, has not been modified significantly at frequencies below the Nyquist frequency.
Because smoothing has already been done to the ESF, a simple first order differencing approach is employed to obtain the LSF. The final step is to obtain the MTF by applying the Fourier transform (via the Fast Fourier Transform) to the LSF and normalizing by the zeroth value.
The various steps in the SDSU MTF estimation process are illustrated Fig. 5. In the upper left corner is an actual satellite image of a slant edge that was obtained from deployed tarps. The upper right chart shows the oversampled ESF that has been produced after application of the modified Savitzky-Golay filter. Note the uniform spacing of the data, even in the steep region of the edge response. At this point an estimation of the SNR for this edge response can be calculated. The lower left plot shows the line spread function after simple first order differentiation data points is observation th function is sh 0.1117. and edge normal vector is obtained, the image gradient magnitude values are projected onto the normal and grouped into coarse bins by their signed distance from the edge. Outliers within each bin are flagged using a variant of Tukey's quartile test. The PCA and outlier rejection steps are iterated until the outliers stabilize. In this context, the PCA provides a totalleast-squares estimate of the edge parameters.

Airbus D
For the ESF/LSF construction step, following the method of Kohm outlined in section 2.2, the coordinates of the samples within the ROI are projected onto the edge normal to yield the set{[z, i(x, y)]} using Eq. (9). Outliers identified in the edge parameter estimation stage are excluded. An ESF with a regular sample spacing of 1/8th pixel is constructed by weighted binning of the set {[z, i(x, y)]}. To avoid poor sampling caused by certain edge slopes (e.g., 1/2 or 1/4), the value of each ESF bin k is calculated according to Eq. (14): where m k denotes the midpoint value (distance from edge) of bin k, and f(d) a low-pass kernel function. In practice, the function f(z) = exp(−13|z|) works well. The effect of this low-pass filter must be removed from the final MTF by dividing the measured MTF by the Fourier transform of f(z), i.e., 13 2 /(13 2 + f z 2 ) in its normalized form. The proposed kernel f(z) is wide enough so that even if the edge angle is 26.565 degrees (a slope of 1/2), none of the ESF bins will have a zero denominator.
The notion of filtering the ESF during construction is taken one step further by switching to a different low-pass kernel for the tails of the ESF, similar to the method proposed by Williams and Burns [14]. In particular, using f(z) = rect(k·z) with k decreasing with distance from the edge reduces the impact of noise on the eventual MTF measurement. The starting locations of the ESF tails are defined relative to the 10% and 90% quantiles of a heavily smoothed temporary ESF. No correction of the MTF is applied for these ESF-tail low-pass kernels.
The LSF is constructed by a finite-difference approximation of the ESF derivative. The ESF-tail smoothing described previously obviates the need for windowing of the LSF before applying the FFT, since the LSF tails naturally taper to zero with sufficient smoothing.
For the final MTF calculation step, the FFT of the LSF is computed, followed by normalization. The appropriate sinc(c·f z ) correction is applied to the final MTF to compensate for the finite-difference approximation. To reduce the variance of the MTF estimates further, a variable-width Savitzky-Golay filter is applied. The width of the filter increases gradually at higher normalized frequency values.

KARI
The KARI slanted-edge implementation has been developed to measure the spatial quality of the KOMPSAT image data starting from the SDSU algorithm [5]. The same implementation was also used for ground testing before launch.
The selection of the inputs for the current KARI algorithm are the ones leading to is the largest values of Relative Edge Response (RER), Function Width at half maximum (FWHM) for LSF and MTF.
The edge modeling involves fitting an edge model, starting with an initial estimate of the edge location within each row of the ROI that is obtained by finding the pair of adjacent pixels with the largest difference. This estimate is refined by computing the inflection point of a cubic polynomial fitted through the four values surrounding the initial estimate. The overall edge model is obtained by fitting a linear function across all the rows through the refined edge locations.

Comparis
Airbus DS, D The available experiment as from actual im 10. ESF parametri qually spaced ESF e requires seve e deformation oise Ratio and isual inspection Modulation Tra ove leads to th LSF to obtain p, is shown in 1. (a) LSF produc uced from the LSF.

son experime
Digital Globe an e data set is s shown Fig. 12 mages.
c function (depict F data points (depi eral manual op applied. Qual d the L2 Norm n remains an im ansfer Function he LSF as sho n sufficient poi Fig. 11 • the curves of the difference relative to the mean MTF or to the known MTF, • a table with the mean, the standard deviation and the largest difference (max-min) at Nyquist frequency. In all graphs and in all tables of this section, ACT stands for across track direction or axis, and ALT stands for along track direction or axis. In 2016, the mean is computed using all the results. In 2017, the mean is computed without the outliers.
The MTF curves will be limited to the across track direction as they provide enough illustration of the discrepancies that were observed.

Results for StdSystem_1 edges
The 2016 and first 2017 (2017a) MTF curves and the discrepancy for the StdSystem_1 edge are presented in Fig. 13. Another graph, Fig. 14, provides the same curves but related to last 2017 (2017b) results. Table 2 gives the 2016 and the last 2017 (2017b) MTF values obtained at Nyquist frequency.    For this first case, the 2017b results are in quite good agreement. The discrepancy between the participants increases with the frequency but remains small up to the Nyquist frequency. This illustrates the value of the comparison for purposes of improving algorithms and removing errors in algorithms.

Results for StdSystem_30 edges
The 2016 and first 2017 (2017a) MTF curves and the discrepancy for the StdSystem_30 edge are presented in Fig. 15. Another graph Fig. 16 provides the same curves but corresponding to last 2017 (2017b) results. Table 3 shows the 2016 and last 2017 (2017b) MTF values at Nyquist frequency.   Initially, for this second case, there were two singular results. The successive reprocessing for case A clearly improved the result. Except for very low frequencies, the discrepancy between the participants is small.

Results for apnn edge
The 2016 and 2017 MTF curves and the discrepancy for the apnn edge are shown in Fig. 17. Table 4 gives the 2016 and 2017 MTF values at Nyquist frequency.  For this case, the curves are in agreement except for D. It appears that this approach has difficulties when processing this type of edge (compare to cgpnn below) which was not solved with reprocessing in 2017.

Results for cgpnn edge
The 2016 and 2017 MTF curves and the discrepancy for the cgpnn edge are presented in Fig.  18. Table 5 shows the corresponding MTF values at Nyquist frequency.  For this case, there is a very good agreement among all results. The only significant deviation comes from C above the normalized frequency 0.8.

Results for 14oct_P3 edges
The 2016 and 2017 MTF curves and the discrepancy for the 14oct_P3 edge are presented in Fig. 1920. Table 6 gives the 2016 and 2017 MTF values at Nyquist frequency.  For this first case with actual satellite image data, there is good agreement between the results except for D at low frequencies. Once again, it looks like a problem related to the type of edge and resulting MTF shape.

Results for 15aug_P3 edges
The 2016 and 2017 MTF curves and the discrepancy for the 15aug_P3 edge are shown in Fig.  20. Table 7 gives the corresponding MTF values at Nyquist frequency.  For this actual case, there is a good agreement between the results except for D at low frequencies, between 0.1 and 0.2.

Conclusion
A test of several algorithms derived from the widely used edge method has been performed. For the test, a set of images of edges was created, mixing both simulations and actual images, and was made available to the participants without any information about the edges, PSF, or MTF. Each participant processed the edge data set to estimate the MTF curve for each edge. Thus, the first comparison of the MTF corresponds to blind test results. A second one occurred one year later which allowed possible improvements to the processes or the inputs.
For either the first or for the second comparison, none of the participants was able to always produce the best estimate (the closest to the expected one for simulation or the closest to the mean of the measurements for the actual cases). This experiment showed that, in some cases, the error or inaccuracy may be MTF shape dependent. Thus, a validation should include several MTF shapes. It also stressed that the results may seem to be consistent when looking at MTF value, at Nyquist frequency, but are not always consistent for the whole curve. Indeed, for some participants, the quality of the assessment depends strongly on the shape of the MTF curve. All participants presented their methods and no theoretical problems were found. The explanation of some unexpected results could possibly be a bug in the software or some inadequate inputs. This emphasizes that a full understanding of the method is required to obtain reliable results.
To extend this comparison study, a reference data set composed of edge image and corresponding MTF curves will be built. It is planned to give access to edge images through the CEOS CalVal Portal (http://calvalportal.ceos.org/) and make them available to a broader audience. In order to promote blind testing as well as to enhance and enlarge this reference data set, it is planned to deliver the reference MTF curves upon receipt of MTF curves from user. Moreover, users are invited to propose new images to enlarge the data set. This paper associated with the reference data set can be seen as a new tool to implement and/or check MTF measurement relying on the slanted edge method.