IFRAD: A Fast Feature Descriptor for Remote Sensing Images

: Feature description is a necessary process for implementing feature-based remote sensing applications. Due to the limited resources in satellite platforms and the considerable amount of image data, feature description—which is a process before feature matching—has to be fast and reliable. Currently, the state-of-the-art feature description methods are time-consuming as they need to quantitatively describe the detected features according to the surrounding gradients or pixels. Here, we propose a novel feature descriptor called Inter-Feature Relative Azimuth and Distance (IFRAD), which will describe a feature according to its relation to other features in an image. The IFRAD will be utilized after detecting some FAST-alike features: it ﬁrst selects some stable features according to criteria, then calculates their relationships, such as their relative distances and azimuths, followed by describing the relationships according to some regulations, making them distinguishable while keeping afﬁne-invariance to some extent. Finally, a special feature-similarity evaluator is designed to match features in two images. Compared with other state-of-the-art algorithms, the proposed method has signiﬁcant improvements in computational efﬁciency at the expense of reasonable reductions in scale invariance.


Introduction
Feature-based registration is a method for matching some similar local regions in preparation for aligning images, which is widely employed by many remote sensing applications [1][2][3] such as image stitching [4][5][6], multi-spectral image fusing and PANsharpening [7,8], and so forth. Despite the diversity of applications, feature extraction is the crucial step prior to feature matching. It consists of feature detection and feature description. Feature detection finds some salient local regions, such as blobs, edges and corners, which are called features, and then quantitatively describes them as feature vectors in the feature description process. Subsequently, some similar feature pairs from two images will be further determined by comparing the similarity of these feature vectors, and the local image shifts are then determined in preparation for alignment.
At present, various methods of feature extraction have been proposed and have shown outstanding performance in different aspects. References [9,10] proposed the Features from Accelerated Segment Test (FAST), which is used to detect blob-or corner-like features according to intensity changes, showing promising performance in speed. Later, reference [11] refined the FAST with a machine learning approach, the performance of which was improved both in repeatability and efficiency. Though the family of FAST was expanded for various usages, it can only detect a series of features in images. To discriminate features, some feature detectors were proposed with their corresponding descriptors, such as SIFT [12], SURF [13], GLOH [14], and so forth. The descriptors exploit the gradient changes to measure orientation and scale-space, and then form feature descriptor vectors to quantitatively describe these detected features. In the feature matching stage, a similarity measure, such as sum of absolute differences (SAD) or sum of squared differences (SSD), is usually applied to determine some feature-pairs between two images; finally, the relative shift between two images can be estimated.
Other types of registration algorithms include area-based registration methods, which estimate the translational shifts by directly exploiting the intensity information to calculate some similarity measurements such as normalized cross-correlation, mutual information [15,16], and so forth. Some variants are also proposed to expand its application range: Reference [17] enabled an area-based method to estimate relative rotations and scales by introducing the Fourier-Mellin transform (FMT); Reference [18] introduced a noise-robust FMT-based registration method by introducing a special filter called the weighted column standard deviation filter. Although they have advantages in stability and precision, their applications in remote sensing are rare due to their time-consuming calculation of similarity measures for large-scale images; moreover, they cannot estimate more general affine or projective distortions, which is more common in remote sensing images.
Nowadays, some registration frameworks that employ optical flow or neural networks have also been proposed. Since optical flow methods can estimate view-difference at the pixel-level [19], they are usually used in video-tracking or motion-estimation. In its application, it is usually combined with feature-based or area-based methods, and some alternatives have also been developed from various aspects [20]. Reference [21] proposed a framework combining feature-based and optical flow to achieve region-by-region alignment. Neural networks-based frameworks are essentially methods of supervised learning, therefore considerable expert knowledge and manpower are needed to prepare handcrafted training data. Besides, neural networks are usually used as tools to boost the performance of feature-based registration methods.
In this work, we propose a novel feature descriptor that can quickly describe each feature point, but also reflect the spatial relationship of each feature, therefore providing the potential for on-orbit applications. Our contribution can be summarized as follows: 1.
Explain a concept for Inter-Feature Relative Azimuth and Distance (IFRAD); 2.
Propose a novel feature descriptor base on IFRAD; 3.
Design a special similarity measure suitable for the proposed descriptor; 4.
Refine the proposed descriptor to improve the scale-invariance and matching accuracy.
The remaining sections are organized as follows: Section 2 reviews the principle of some state-of-the-art methods, wherein the main differences of our proposed method from others are described; in Section 3, we propose the IFRAD descriptor, where its principles are explained. Experiments are conducted in Section 4 for finding the optimized parameters, as well as the comparison between our descriptor and others; then, some drawbacks of IFRAD and their causes are discussed in Section 5; finally, we draw conclusions in Section 6.

Background
The state-of-the-art feature extraction methods can detect features and describe them discriminatively while keeping affine-invariant to some extent, such as SIFT [12], SURF [13], GLOH [14], BRISK [22], ORB [23], FREAK [24], and their alternatives [25,26]. However, they need to describe the features according to their surrounding pixels or gradient changes, which will limit the computational efficiency. With SIFT as an example, for a given image I(x, y), it first generates scale space images L (I) (x, y, σ) by convoluting I(x, y) with different standard deviations σ of Gaussian function G(x, y, σ): where G(x, y, σ) = 1/(2πσ 2 ) exp −(x 2 + y 2 )/2σ 2 . Then, a difference of Gaussian (DoG) of k-scale can be computed by: where k is a constant difference of two nearby scales. By stacking these DoG, an octave of a scale-space image block is formed. Meanwhile, with a smaller size of images downsampled from I(x, y), scale-space image blocks of various octaves are also formed to build a DoG pyramid, from which the scale space local extrema and their locations can be roughly estimated. These extrema are called feature point candidates. To improve the stability: (1) accurate feature point localization is performed by interpolating locations of these extrema with a 3D quadratic function; and (2) the edge response is eliminated according to a 2 × 2 Hessian matrix. In the feature description step, for each L (I) (x, y, kσ) in the pyramid, the gradient magnitude and orientation of each pixel are precomputed according to: Subsequently, for each previously determined feature point f (I) i , its local orientation histogram is formed, the peak of which is assigned to f (I) i as a dominant orientation. Finally, the histogram is remapped according to the peak, and is then formed as a feature vector to describe the feature.
SIFT shows outstanding performance in scale-and rotation-invariance, and other algorithms that adopt a similar strategy have been proposed. Compared with SIFT, SURF adopts several operations to speed up the algorithm, such as integral image and approximation of second order Gaussian derivatives, and so forth. Recently, extensive literature has proposed some feature descriptors with improvements in affine invariance, computational efficiency and accuracy [1,2,[27][28][29]. Despite their robustness in affine-invariance, their heavy dependency on exploiting gradient changes and image pyramids results in the complexity of these algorithms, which limits the applications in hardware platforms, therefore many proposed methods for remote sensing applications have been used only in off-line or post-processing applications [2,[30][31][32]. In pursuing real-time applications, some GPU-accelerated algorithms were proposed to reduce the time consumption. Reference [33] proposed GPU-accelerated KAZE, which is about 10× faster than the CPU-version of KAZE; reference [34] implemented a GPU version of SIFT, with an accelerated factor of 2.5×. However, GPU is merely a tool for acceleration; it does not reduce the complexity of these algorithms. Moreover, spaceborne cameras are usually boarded with field-programmable gate arrays (FPGA), in which only simpler algorithms can be implemented, and few real-time remote sensing applications that utilize feature-based methods have been reported [35]. Therefore, a simpler algorithm for feature extraction has become an urgent need for on-orbit applications.
Unlike those aforementioned methods, IFRAD does not need to form any image pyramid, and the gradient change is only precomputed in the feature detection stage before IFRAD, which greatly alleviates the burden of computation. With a series of detected features in an image beforehand, IFRAD describes each feature according to its relations (relative azimuths and distances) to other features. Its robustness in scale invariance is assured by cosine similarity in the feature matching stage.

Methology
In this section, we will implement the IFRAD descriptor, its corresponding similarity measure and its application in image registration. The overall flowchart is shown in Figure 1. As mentioned in the last section, feature detection should be performed before IFRAD; for the sake of computational efficiency, we use FAST to detect a series of features. We assume that we have two images: Figure 2a shows a reference image R(x, y), and Figure 2b shows a sensed image S(x, y). Compared with Figure 2a, Figure 2b has a parallax caused by about 20 • yaw, and about 30 • pitch; their parallax difference is shown in Figure 2c.  However, some feature detection methods with low computational complexity, such as FAST, are prone to interference from some random factors, such as noise and image distortion. Therefore, not all detected features are stable enough to ensure the repeatability of our IRFAD descriptor since it will describe the inter-feature spatial relationship of each feature. Therefore, it is necessary to select some stable features according to criteria.

Criterion for Selecting Secondary and Primary Features
Before obtaining FAST features, Gaussian smoothing is utilized to reduce the effect of noise. When detecting features in the reference image R(x, y), the FAST detector will denotes the ith feature. Since IFRAD will describe features according to spatial relationships, the features near to the image center are more likely to be properly described, as they can be described according to other features from all directions, which is impossible for the features near to the edges or corners. Therefore, we should modulate the response magnitude of a feature according to its distance from the image center, that is: where Mag m f ); then we obtain a secondary-feature set F (R,SF) that contains the strongest half of these features. This operation can be expressed by: For simplicity, in the remainder of the paper, we use f . The process of selecting secondary features is illustrated in Figure 3. However, in some circumstances, some patterns (such as corners or blobs) may have an uncertain number of secondary features (as shown in Figure 4), which may vary to different parallax, or the random distribution of noise. This will cause errors in the feature matching process. To reduce this effect, we further determine primary features from these secondary features according to the following steps:

1.
Initial F (R,PF) as an empty primary-feature set, that is, : The radius R is an adjustable parameter, and will be further determined in experiments.

3.
If there exists no other secondary-feature f . This criterion can be expressed by: The criterion is clarified in Figure 5; under this criterion, Features A,B,C,F,J,K,M,N will be determined as primary features. In the remainder of the paper, we also use f for simplicity. The primary feature selection results of the reference image are shown in Figure 6.   According to these criteria, the relation of these feature sets are as follows:

The Relationships among Features
The relations of one feature to the others include relative azimuth and distance (RAD). Assuming that we have a primary feature f (R,PF) i , its relative azimuth and distance to a secondary feature f (R,SF) j can be obtained by: When the image has a certain degree of distortion, those secondary features that are farther away from f (R,PF) i will have greater changes in the spatial relationship. Thus, the strength of the relationship is expressed by:

The IFRAD Descriptor: Orientation Intensity Histogram
To make our method rotation-invariant, we determine the dominant orientation of each primary feature. For that purpose, we calculate the RADs of all secondary features relative to the primary feature, then collect them into a list, followed by sorting them in descending order of relation strength (ascending order of relative distance), as shown in Figure 7b.
Then, the dominant orientation of the primary feature is the relative azimuth of its nearest secondary feature (also have the strongest relation), that is: where f With this orientation as the reference (set to 0), the relative azimuths of remaining secondary features are remapped to the range of [0, 2π), the remapped azimuths are denoted as Azim rm f , as illustrated in Figure 8. On this basis, we can obtain an n-bin-orientation intensity histogram (OIH) by the following steps: 1.
Divide the image into n fan-shaped regions according to the orientation of the feature, where n is an adjustable parameter and needs to be optimized in experiments, as shown in Figure 9a, in this example, n = 10; 2.
Estimate the orientation intensity by calculating the sum of all relation-strengths within each fan-shaped region. This operation can be expressed by: where I k f is an indicator function: With the above steps, an OIH of each primary feature is formed (shown in Figure 9b), and can be presented as a vector, which will be used as our IFRAD descriptor vector:

Feature-Matching
In the feature matching process, a proper similarity metric is a critical factor for a higher matching correctness rate. For gradient-based feature descriptors, such as SIFT, SURF, and so forth, SSD or SAD is often used as a similarity measure. Another similarity measure, called Hamming distance, is also used for binary feature descriptors [36]. The rotation-invariance of the proposed descriptor is implemented by remapping the relative azimuth to a fixed range, with the dominant orientation set as 0. Moreover, the IFRAD descriptor vector also has the potential for scale-invariance, for the variations in scale will simultaneously and proportionally change the inter-feature distances, which implies that the "shape" of OIH will not change. However, the potential cannot be exploited by the aforementioned measures. Therefore, we adopt cosine similarity as the similarity metric of OIH.
Assuming that we have two similar features-one is in the reference image, its OIH is represents the other feature in the sensed image, then the cosine distance (the smaller the better) of these two features is expressed by: where . By cosine similarity, the scale-invariance can be improved. While matching feature pairs, it is inevitable to mismatch some features due to the difference in the description of the same feature caused by different views. This will produce a larger portion of outliers (mismatched feature pairs) for estimating view-differences. To mitigate the influences, we need to restrict the matched conditions, that is, the two

Refinements
OIHs determine the distinguishability of the features, and can be improved by increasing the number of bins in OIH. This can be interpreted by analogy with the bit-depth of an image; the higher the bit-depth (which results in a higher number of bins in the intensity histogram), the more distinguishable are the details in computer processing. However, this comes with the cost of reducing the stability. There exist some circumstances where a primary feature to be matched is surrounded by a series of secondary features, and the nearest of them have similar distances; however, the inter-feature distances often vary in different parallaxes, which may greatly affect the determination of the dominant orientation, and can further reduce the replicability. As shown in Figure 10, with another feature-pair as an example, they have similar bar graphs (Figure 10a,d) with differences in relation strength and azimuth offset. However, since the dominant orientation is determined according to Equation (13), the changes in relation strength caused by different parallaxes have a great effect on determining the dominant orientation, making it so that their OIHs do not match (Figure 10c,f. According to (18), the cosine distance between them is 0.3293. (c,f) are 30-bin-OIHs obtained from (b,e), according to (18), the cosine distance between them is 0.3293.
To smooth the issue, we determine the dominant orientation with the average of the azimuth of secondary features which have a relation strength that is comparable with the strongest one, which means we replace Equation (13) with: where f where α ∈ [0, 1] is the coefficient for selecting f (R,SF) k with a stronger relation, and the q in (19) is the number of f (R,SF) k that satisfy the constraint (20). In this way, the issue can be handled, as shown in Figure 11; in this example, α = 0.7, note that three peaks (marked with red    ) in the bar graph satisfy the constraint (20), and are counted for determining the dominant orientation. The dominant orientations are stably determined, and thus the cosine distance is decreased to 0.0254. Orientation Intensity (f) Figure 11. An Example of Solving the Unstable Issue: (a,d) are the same bar graphs shown in Figure 10a,d, but different in dominant orientations, and are determined according to Equation (19). In this example, α = 0.7; (b,e) are the remapped bar graphs according to (a,d); (c,f) are 30-bin-OIHs obtained from (b,e), according to (18) the cosine distance between them is 0.0254.
For the same reason, to improve the replicability of the determination of primary features, we modify the criterion (8) by applying tolerance t: where t 1 is the tolerance for allowing some features with a competent response in a determined feature domain. Under this criterion, with t = 0.8, in Figure 5, features A,B,C,D,F,G,J,K,M,N will be determined as primary features.

Geometric Transform Estimation and Correctness of Matching
In the registration process, we estimate a 3-by-3 transform matrix M T from these matched IFRAD-descriptor-vector-pairs using the method of MLESAC [37], which can also tell correctly matched feature pairs (inliers) apart from mismatched feature pairs (outliers), and the correctly matched rate (CMR) can be further obtained by: where N In denotes the count of correctly matched feature pairs, and N Out denotes the count of the mismatched feature pairs. Since MLESAC is a generalization of the Random sample consensus (RANSAC) estimator, its randomness of initial feature-pairs selection can result in some fluctuation of the values in the estimated geometric transform matrixes. However, the following experiments proved that the diversity can be reduced to an acceptable range by selecting proper values of parameters. Figure 12 shows an example of IFRAD-based registration with parameters of α = 0.7, n = 30, t = 0.8 and R = 1/100 · min(M, N). Figure 12a shows correctly matched feature-pairs; Figure 12b shows the registration result; and some magnified views of local images are shown in Figure 12c. In this example, the estimated geometric transform matrix is:

Experiments
As mentioned above, we have four parameters to be optimized: (a) the tolerance t; (b) the coefficient α; (c) n, the number of bins in OIH; (d) the radius of feature domain R. Different parameter settings will affect the correctness of geometric transform estimation. In our experiments, we use five groups of larger-scale remote sensing images; each has two images with different parallaxes (Figure 13). We first define some assessments to find the optimized values for each parameter, then the optimized IFRAD descriptor is compared with other state-of-the-art algorithms. Finally, the limitations and range of applications are discussed in Section 5.

Assessments Defination
As mentioned in Section 3.6, due to the randomness of initial feature-pairs selection, the M T may vary for each time of estimation under the same parameters. However, the diversity can be reduced by selecting proper values for parameters, as it will increase the proportion of inliers (correctly matched feature-pairs). Therefore, we quantify the diversity by stability of transform estimation (STE): For each group of parameters with determined values, after estimating the transform matrix for k times, the STE can be quantified by: where c (k) ij represents the ith row, jth column element in the kth time estimated transform matrix M (k) T , STD(·) denotes the standard deviation, Mean(·) denotes the average. Obviously, the higher the STE, the more stable the transform estimation.
Apart from STE, the CMR (see Equation (22)) and computational costs are also important assessments for evaluating our methods.

Parameters Optimization
We first use the five groups of remote sensing images ( Figure 13) to determine the optimized parameters. Since they are noise-free images, the variety of determination of secondary features will only be caused by different parallaxes, therefore we allow R to be a fixed value relative to the size of the image, R = 1/50 · min(M, N), where M and N are the width and height of an image, respectively. Then we obtain the values of these assessments under different n, α and t; the sample size for each parameter's setting is 100. As shown in Table 1, the optimized parameters are: n = 50, α = 0.6 and t = 0.8.
In reality, raw images short-captured at rapid motions may have stronger poisson noise, which may increase the portion of outliers in estimating the relation of two views. Therefore, with optimized n, α, t, we conduct some experiments to find the relations between R and relative noise level (RNL) by comparing the CMR, where RNL is quantified by adjusting exposure time (ET); the shorter the ET, the higher the RNL. The experimental results are shown in Table 2; the F in the left column represents the factor of radius of feature domain: We can draw the conclusion from Table 2 that the optimized radius will get smaller as ET increases (making RNL lower). It also indicated that the shorter the ET, the greater the variation in CMR as F changes. Moreover, the higher the CMR, the more stable the estimation of the transformation matrix. In fact, further experiments indicated that, when CMR > 0.30, the risk of registration failure can be eliminated due to the ease of inlier/outlier discrimination for the MLESAC algorithm. Considering CMR and the ease of tuning, we determine R = 1/20 · min(M, N), (F = 5) uniformly for any level of noise. The registration results are shown in Figure 14.

Comparisons
We also compared the performance of our descriptor (with optimized parameters) to other state-of-the-art methods in terms of computational cost, scale-invariance, CMR, and so forth. The comparisons were all performed on a PC with Intel Core i7-7700 CPU @ 3.6 GHz, and RAM of 32 GB. The results are given in Table 3-5. In those comparisons, the experimental image groups in Figure 13 were used.
According to Table 3, we can conclude that, for small scale changes (1.00∼1.25×), the CMR of our descriptor is comparable with that of other methods, but will drop rapidly when the scale change is higher than 1.30×. However, the MLESAC can still discriminate a group of inliers to correctly estimate the transform matrix until the scale change reaches 1.50×. Moreover, IFRAD has the highest time-efficiency compared with others (according to Table 4); this is achieved by the many fewer matched-pairs counts (Table 5) due to the strict selection of primary features and the determination of matched feature pairs.   SURF  1834  1715  1463  1290  1117  1015  1008  942  924  882  748  690  618  KAZE  32,398 29,941 23,977 14,521 7003  4737  6023  7904  8857  8407  7238  5698  4089  BRISK  6309  5713  4747  3655  2905  2625  2403  2395  2341  2122  1912  1735  1483  IFRAD  896  832  785  702  645  587  518  396  175  84  32  7  4 The numbers in red bold indicate the least number of matched pairs under the current scale.
Since IFRAD describes features according to their relations to other features, the translational shift between two images may move some features outside the image boundary and will alter the feature description. As shown in Figure 15, the second column shows that, as a feature moves to a border or corner of the image, it will tend to be partially described in relation to other features from one side or quadrant, which will reduce the reliability of the description. Since the translational shifts change the percentages of overlapping areas between two images, the experiments comparing the CMRs of different methods under various percentages of overlap area are also conducted, and the results are tabulated in Table 6, and indicate that our method is reliable when the overlap area is greater than 50%.  The numbers in red bold indicate the best results.

Discussion
In the last section, we compared our IFRAD to other methods in several ways-scale invariance, time consumption, and CMR. The experiments showed that, with a scale change below 1.45× and an overlap area of over 50%, the proposed method has superiority in time efficiency while keeping correctness in estimating the transform matrix. However, several drawbacks limit its application range. While other methods such as SURF-which exploits the scale information by forming a scale-space with several spatial-octaves-have a wider range of scale-invariance up to 8×, the IFRAD descriptor has a narrower range of scale-invariance. Although the cosine distance is adopted as a similarity distance in the feature-matching process, the range of scale-invariance is increased by only up to 1.4×. The limitations are mainly brought by Equation (12), which neglects the fact that the response magnitude of a FAST-feature always varies with different parallaxes, and will greatly reduce the scale-invariance of OIHs. Figure 16 depicts the causes of the limitation. For two images with larger scale differences, some features in the small-scale-image may appear as low-frequency information in the large-scale-image, and some features in the large-scale-image may disappear in the small-scale-image. Further experiments show that this drawback can be manipulated by replacing the FAST-feature with more stable features that utilize the information about the scale-space; however, this brings with it the cost of increased time-consumption.
Feature responses are sensitive to illuminance changes, thus the selection of secondary features or the determination of primary features may differ due to the changes in feature response magnitude, making IFRAD unable to stably describe the features as illuminance changes. This will cause failure in registering images from different sensors or different spectrums.
In summary, the drawbacks of our proposed methods are listed as follows, compared with other state-of-the-art methods: • the scale-invariance of our method is limited within 1.45×; • the range of applicable overlap area is narrowed to 50∼100%; • more sensitive to illuminance changes.
Despite these drawbacks, in reality, for CMOS-based push-broom remote sensing images, the altitude of the airborne camera is stable while scanning along the track and the range of scale change is narrow enough for IFRAD to perform a quick and reliable feature description and matching process. Therefore, our method still has the potential for implementing some applications, such as on-orbit image stitching, registration-based TDI-CMOS [18,38], and so forth.

Conclusions
In this paper, we proposed a feature description called IFRAD. We first introduce some criteria for selecting secondary and primary features to improve the robustness of our IFRAD descriptor. Then concepts about inter-feature relative azimuth and distance are provided, based on which we further explain the algorithms for obtaining the IFRAD descriptor vector called n-bin-OIH. In feature matching, the cosine distance is introduced to improve the scale-invariance. To improve the replicability, we made some refinements by introducing two parameters: coefficient α in Equation (20) and the tolerance t in (21). In estimating geometric transform with MLESAC, we defined an assessment called CMR. Due to the randomness of MLESAC in selecting the initial matched-feature-pairs, the estimated geometric transform matrix varies in each time estimation. However, the diversity can be reduced as the portion of inliers increase. Therefore, STE is introduced to quantify the diversity, and is further used to optimize three parameters: the coefficient α, the tolerance t and the number of bins in OIH n. Table 1 shows that the optimized parameters are n = 50, α = 0.6, t = 0.8. The radius of feature domains R is optimized by comparing CMRs under different R and Poisson noise, which are quantified by adjusting exposure time. The experiments have proven that IFRAD can alleviate the effects of denoising. Table 2 shows that the optimized R decreases as the exposure time increases. However, adjusting R produced a negligible change to CMR; therefore, R is set to be 1/20 · min(M, N) uniformly for simplicity. The comparisons to other methods are also conducted, and the results are tabulated in Table 3-5, indicating that IFRAD has the highest time efficiency with reasonably reduced scale-invariance (up to 1.45×). Table 6 indicates that the proposed method is reliable when the overlap area is above 50%.
The proposed method IFRAD is only a feature descriptor. It is designed to simplify the process of feature description and then to speed up the feature matching step. As an alternative to other descriptors such as SURF, IFRAD has several aspects that can be improved. For instance, Equation (12) limits the scale-invariance up to around 1.4×, because it neglects the fact that the response magnitude of a FAST-feature always varies with different parallaxes; the R in (24) is set as constant, which also contributes to the limitation. Before applying IFRAD, in this work, Gaussian smoothing is applied to alleviate the influences of noise. While smoothing, some critical details may be lost. Although some denoising methods that can retain detail were proposed [39], they add to the timeconsumption. Therefore, further improvement of the noise-robustness of our method is critical for low-light remote sensing applications. The solution to these issues, as well as the GPU version of our method, will be focused on in our future work.