A Robust Ground‐to‐Image Transformation Algorithm and Its Applications in the Geometric Processing of Linear Pushbroom Images

Linear array cameras are widely used in the fields of aerial and space remote sensing. The ground‐to‐image transformation plays a significant role in the geometric processing procedures of linear pushbroom images, which is time‐consuming due to multiple iterative calculations. Moreover, the detectors of the linear array may exhibit a curved shape due to distortions, resulting in a more complex ground‐to‐image transformation. This paper presents a novel ground‐to‐image transformation algorithm for linear pushbroom images based on a rigorous sensor model. We regard the distance between the back‐projected image point and the curved shape linear array as the “generalized distance.” The idea of generalized distance prediction was developed to efficiently determine the best scan line corresponding to the ground point. The plane analytic geometry is employed to solve the practical problems of the ground‐to‐image transformation, including the iterative direction, convergence condition, and subpixel interpolation. The applications of the proposed method in orthorectification and image matching were investigated. The proposed method was comprehensively tested using airborne images, Earth observation satellite images, and planetary images. The experimental results demonstrate that the proposed method (1) delivers high accuracy (better than 1 × 10−4 pixels) and efficiency, (2) is better than the cutting‐edge method based on the geometric constraints of central perspective plane for the linear pushbroom images with distortions, and (3) provides promising applications for generating orthophotos and image matching. The proposed method is easy to implement and efficient, which can greatly enhance the geometric processing abilities for various types of linear pushbroom images.


Introduction
Over the past 10 years, many high-resolution Earth observation satellites have been successfully launched (Wang et al., 2017;Zhang et al., 2014). Currently, very high resolution satellites such as the Pleiades twins (1A/1B; Poli et al., 2015) and WorldView 3/4 (Ferreira et al., 2019) can provide remote sensing images with spatial resolutions of better than or equal to 0.5 m, allowing for cartographic products at the 1:10,000 or even 1:5,000 scale. Linear array cameras are commonly used in the field of remote sensing, due to the fact that they can provide higher spatial resolution images and wider field of view. It is widely acknowledged that the geometric processing capabilities for remote sensing images greatly lag behind the data acquisition capabilities. Although parallel computing can improve the processing efficiency, the algorithm design should be prioritized.

Background
The sensor model lays the foundation for the geometric processing of linear pushbroom images. Commonly used sensor models for linear pushbroom images include the rational functional model (RFM; Dial et al., 2003) and rigorous sensor model (RSM; Poli & Toutin, 2012). The RSM relates three-dimensional (3-D) ground coordinates to two-dimensional (2-D) image coordinates based on the physical characteristics of the camera. RFM represents the image-to-ground transformation and ground-to-image transformation of a sensor model as a rational polynomial. The RFM is widely used for Earth observation satellite images due to its speed, generality, and confidentiality. The RSM has wider applications than the RFM and can be used to process satellite images as well as airborne images (Gonzalez et al., 2013;Tong et al., 2015). Moreover, the RSM has higher accuracy and physical meaning. The best practice is using the RSM for triangulation due to its high accuracy and using the RFM to generate mapping products because of its high efficiency (Förstner et al., 2013). Additionally, the RSM is widely used in the planetary mapping community (Di et al., 2014;Haase et al., 2019;Kirk et al., 2008;Wu et al., 2014) because a probe's position and attitude observations can be directly introduced into the bundle adjustment model to improve the stability of the normal equation.
For linear pushbroom images, each scan line is acquired at a different time and the exterior orientation (EO) parameters change from scan line to scan line due to the pushbroom imaging principle (Granshaw, 2016). The central perspective projection applies in the across-track direction of each scan line. Thus, the ground-to-image transformation based on the RSM, which is also referred to as back projection, requires multiple iterative calculations to determine the best scan line. Consequently, the ground-to-image transformation for linear pushbroom images is always time-consuming, limiting the rapid acquisition of cartographic products. Nevertheless, the ground-to-image transformation plays a significant role in the geometric processing procedures of linear pushbroom images. For example, the orthorectification implemented with indirect method depends on the ground-to-image transformation to determine the corresponding pixel coordinates on the original images. Additionally, image matching for tie point extraction and digital terrain model (DTM) generation is always based on epipolar resampled images or approximate orthophotos (Beyer et al., 2018;Heipke et al., 2007). Then, the matched conjugate points on these intermediate images can be converted to conjugate points on the original images by ground-to-image transformation or keeping the transformation information in temporary images.

Existing Ground-to-Image Transformation Methods
It is noted that some studies have been carried out on this topic. Kim et al. (2001) proposed a ground-to-image transformation algorithm based on random estimates of the 2-D image coordinates. Though the authors provide some insight information, the algorithm model is too idealistic. The classical bisection method adopts the binary search principle to search for the best scan line corresponding to the ground point (Liu & Wang, 2007). The bisection method has advantages of robustness and being easy to implement, but it has low efficiency because the algorithm complexity of the binary search is O(log 2 N). Wang et al. (2009) proposed a ground-to-image transformation method based on the geometric constraints of the central perspective plane (CPP) of the scan line, which shows significant advantages in computational efficiency. However, in the case of a linear array camera with distortions, the CPP of the scan line becomes a 3-D curved surface instead of a 3-D plane. Consequently, to apply the CPP-based method, the curved shape linear array should be segmented into multiple short line segments, such that the exact CPP corresponding to the ground point can be determined (Geng et al., 2019). These additional procedures cause the CPP-based method to be relatively difficult to implement for the linear pushbroom images with distortions. Shen et al. (2015) developed a ground-to-image transformation method with a coarse-to-fine search strategy for airborne pushbroom images. The authors employed a bilinear transformation to estimate the approximate scan line corresponding to the ground point. The Pleiades imagery user manual (Pleiades, 2019) presents a ground-to-image transformation method that estimates the image position by means of polynomials at different altitudes. However, this method is especially designed for the linear pushbroom images without distortions. Therefore, its application scope may be limited. In summary, many factors such as robustness, efficiency, and convenience influence the practical effects of the ground-to-image transformation algorithm. It can be inferred that the ground-to-image transformation algorithm, especially for linear pushbroom images with distortions, is still incomplete.

Proposed Method
This paper presents a novel ground-to-image transformation method for linear pushbroom images based on the RSM and aims at improving the robustness and efficiency of this key algorithm. The characteristic of the proposed method is that it is more applicable for the practical linear pushbroom images with varying degrees of distortions. We first examine the fundamental characteristics of the ground-to-image transformation. Then, the idea of generalized distance prediction (GDP) for determining the best scan line corresponding to the ground point is proposed. The plane analytic geometry in the image space is employed to solve the practical problems of the ground-to-image transformation, which decreases the complicated iterative calculations. Applications of the proposed method in the geometric processing procedures such as orthorectification and image matching were also investigated. We test the proposed method with Earth observation satellite images, planetary images, and airborne images. The experimental results demonstrate that the proposed method supports various types of linear pushbroom images and delivers higher computational efficiency and geometric accuracy.

Basic Principle
For frame images, the collinearity equation can be directly used to conduct ground-to-image transformation. In contrast, due to the time-dependent nature of linear pushbroom images, the collinearity equation should be extended to map a 3-D ground point to a 2-D image point. Considering the influences of image distortions, the extended collinearity equation for linear pushbroom images can be written as where (x,y) are the 2-D image coordinates, (x 0 , y 0 ) are the principal point offsets, (Δx, Δy) are the image distortions, f is the camera constant (i.e., focal length), (X,Y,Z) are the 3-D ground coordinates, r is the best scan line corresponding to the ground point, and (X r S , Y r S , Z r S ) and a r j ; b r j ; c r j (j = 1,2,3) indicate the position component and rotation matrix elements of the EO parameters, respectively. This extended collinearity equation is the basic mathematical model for various photogrammetric operations of linear pushbroom images. Obviously, the best scan line r should be determined prior to the real ground-to-image transformation operation. Once r has been determined, the image coordinates corresponding to the ground point can be calculated using the EO parameters of r.

Fundamental Characteristics of Ground-to-Image Transformation
Before introducing the proposed method, we will first examine the fundamental characteristics of groundto-image transformation. Figure 1 illustrates a series of back-projected image points during the iterative process of ground-to-image transformation. This figure was drawn based on the work of Liu and Wang (2007). For convenience, a linear array camera without distortions is adopted in the figure, whereas the summarized conclusions apply to the linear array cameras with distortions. Let r and l i (i = 1,2,3,4) be the best scan line and the approximate scan line corresponding to the ground point P(X,Y,Z), respectively. The back-projected image points with respect to r and l i are represented by p and p i , respectively; S r and S i represent the perspective centres of r and l i , respectively; and T r and T i represent the exposure times of r and l i respectively. For each scan line, a corresponding image plane can be formed. Thus, the distance between p i and the linear array in the image plane can be calculated and is represented as d i . The distance between the exact back-projected image point p and the linear array is represented as d r , and it is obvious that d r equals zero. It is observed that |d 1 | > | d 2 |and |d 4 | > | d 3 |. Accordingly, the approximate scan line l 2 is closer to the best scan line r than l 1 , and l 3 is closer to the best scan line r than l 4 . Additionally, note that p 1 and p 2 are on the right side of the linear array, p 3 and p 4 are on the left side of the linear array, and p is on the linear array. Accordingly, the distance d i is a signed number, which may be positive or negative.

Earth and Space Science
As shown in the upper parts of Figure 1, the relationship between d i and T i can be regarded as the relationship between d i and the best scan line r. Thus, by analysing the relationship between d i and r, the following fundamental characteristics of ground-to-image transformation can be inferred.
1. The smaller the distance d i , the closer the approximate scan line l i is to the best scan line r. 2. The relationship between the distance d i and the exposure time T i (equivalent to the best scan line r) has monotonically increasing or decreasing characteristics. 3. If the distance d i and d j (i,j = 1,2,3,4,and i ≠ j) have different signs, the best scan line r should be located between l i and l j . For example, the sign of d 2 is different from the sign of d 3 ; thus, the best scan line r should be located between l 2 and l 3 .
The first and second of these fundamental characteristics indicate that the distance between the back-projected image point, and the linear array can be used to predict the best scan line. The third fundamental characteristic enables the subpixel accuracy of the best scan line to be obtained by linear interpolation using two neighboring scan lines. Moreover, the third item can be used as an end flag of the iterative process. Specifically, if d i and d i+1 have different signs, then the best scan line must be located between these two neighboring scan lines; thus, the iteration process should be terminated.

Generalized Distance Prediction
Ground-to-image transformation refers to the coordinate transformation from 3-D ground points to 2-D image points. The key to ground-to-image transformation is that the best scan line corresponding to the ground point should be determined accurately and efficiently. At the initial stage of the ground-to-image transformation, we have no knowledge about the best scan line. Indeed, the best scan line can be iteratively updated using a prediction strategy. Therefore, the specific prediction strategy of the best scan line has strong influences on the efficiency and robustness of the ground-to-image transformation operation. Based on the fundamental characteristics of ground-to-image transformation, the basic principles of the proposed method are illustrated in Figure 2. We use F(d i ) to express the scan line prediction approach used in the iterative process of ground-to-image transformation. The proposed ground-to-image transformation approach based on GDP can be regarded as a framework because different prediction strategies can be developed based on it. It is noted that scan line prediction strategies based on both the image space (Kim et al., 2001) and the ground space (Wang et al., 2009) have been adopted. Therefore, we

10.1029/2019EA000646
Earth and Space Science provide two scan line prediction strategies below. The scan line prediction strategy based on the image space has a simple formula and can be expressed as where a indicates the pixel size. The scan line prediction strategy based on the ground space can be expressed as where D i indicates the corresponding distance in the ground space, which can be calculated by the product of d i and the scale factor of the perspective projection; and GSD indicates the ground sample distance (GSD).
Both scan line prediction strategies have advantages and disadvantages. The scan line prediction strategy based on the image space is easy to implement but requires more iterations. In contrast, the convergence speed of the scan line prediction strategy based on the ground space is generally faster than that based on the image space, but it is difficult to implement for the linear pushbroom images with distortions. Specifically, due to the varying exposure time or orbit height, the GSD of linear pushbroom images usually changes in real time. Moreover, the distance D i in the ground space is usually not accurate due to image distortions and the resulting curved surface shape of the CPP. Therefore, the choice between these two prediction strategies is a trade-off between convenience and efficiency. In this paper, we prefer to employ the scan line prediction strategy based on the image space due to its convenience and because it is more suitable for the practical linear pushbroom images with distortions. Furthermore, its disadvantage of relatively low efficiency is compensated for by making full use of an approximate value of scan line.

Plane Analytic Geometry
The proposed ground-to-image transformation method employs the principle of plane analytic geometry. This has advantages of (1) high computational efficiency because the formulas for plane analytic geometrical calculations are very simple and (2) simplifying complicated situations caused by different flight directions and iterative directions. To examine the plane analytic geometry of ground-to-image transformation, we plot p and p i that were first shown in Figure 1 on the same image plane, and the results are shown in Figure 3.
Here, a linear array camera without distortions is used for convenience.

Earth and Space Science
The distance d i between an image point p i and the linear array is calculated by the following equation where A,B,and C are the parameters of the line equation for the linear array. Note that the sign of d i indicates the position of p i relative to the linear array. Specifically, as shown in Figure 3, d i > 0 means that the back-projected image point p i is located on the right side of the linear array, d i < 0 means that p i is located on the left side of the linear array, and d i = 0 means that p i is located on the linear array. These three situations can be expressed by d i > 0; p i is located on the right side of the linear array (5) The sign of d i can also be utilized to determine whether two image points are located on the same side of the linear array. For example, d 1 × d 2 > 0 implies that p 1 and p 2 are located on the same side of the linear array, whereas d 2 × d 3 < 0 implies that p 2 and p 3 are located on different sides. Therefore, the sign of d i is beneficial to determining the iterative direction, which can ensure the stable convergence of the ground-to-image transformation.
These rules that derived from a linear array camera without distortions are also applicable to the linear array cameras with distortions. Note that for a linear array camera with distortions, the distance between a backprojected image point and the curved shape linear array should be the generalized distance. Additionally, according to the plane analytic geometry, we can specify a reference scan line to facilitate the determination of the iterative direction. This will be discussed in detail in section 2.2.1.

Calculating the Generalized Distance
As shown in Figure 4, the curved shape linear array consists of N detectors. Thus, N − 1 tiny line segments can be formed using two neighbouring detectors. These tiny line segments are shown by alternating short red and blue lines in Figure 4. We apply the curve approximation method to express the curved shape linear array. Thus, the curved shape linear array can be expressed by a series of tiny line segments. Consequently, the distance between the back-projected image point and the exact tiny line segment can be treated as the generalized distance between the back-projected image point and the curved shape linear array. In this way, the computation of the generalized distance focuses on the determination of the exact tiny line segment corresponding to the back-projected image point.
It should be noted that the curved shape linear array is very similar to a straight line. Moreover, the linear array is always mounted perpendicular to the flight direction. Thus, with the image coordinates of the back-projected image point p(x,y), the corresponding calibrated image point p ′ (x′, y′) on the linear array can be easily determined. Next, the coordinate transformation from p ′ (x′, y′) to the pixel coordinates p(s,l) should be performed. It is worth noting that such imageto-sensor transformation also requires iterative computations due to image distortions. Usually, the image-to-sensor transformation for a curved shape linear array consists of two steps: (1) the localization of the s-component of p(s,l) with integer precision by binary search  algorithm and (2) the subpixel accuracy refinement of the s component by interpolation. Let s 0 indicate the integer precision of the s component and s 1 indicate the detector next to s 0 . Accordingly, the y component of the image coordinates corresponding to s 0 and s 1 are represented by y 0 and y 1 , respectively. Thus, the subpixel accuracy of the s component is achieved by the following interpolation equation Note that the l component of p(s,l) remains unchanged during the image-to-sensor transformation because it was already calculated during the ground-to-image transformation process.
Although the localization of the s component of p(s,l) requires iterative computations, the computational efficiency is indeed very high. This is because the calibrated image coordinates of each detector on the linear array are already known, such that the determination of the s component of p(s,l) is a binary search problem for an ordered table. With the determined pixel coordinates p(s,l), the neighbouring detectors with respect to the back-projected image point p(x,y) are calculated. Accordingly, the exact tiny line segment corresponding to p(x,y) is determined. Obviously, this is a typical problem of using two points to solve the line equation. Consequently, the generalized distance between the back-projected image point and the curved shape linear array can be easily calculated using equation (4). Note that the line equations of the tiny line segments can be calculated in advance. Therefore, the calculation of the generalized distance exhibits high computational efficiency because the involved calculations are simple analytic geometrical computations.

Practical Considerations
This section examines the practical considerations of the proposed GDP-based ground-to-image transformation method, including the iterative direction, convergence condition, subpixel interpolation, and employment of an approximate value.

Iterative Direction
Let Δr be the incremental value of the scan line that is predicted during the iterative process of the ground-toimage transformation. Thus, the best scan line r is iteratively updated by where i indicates the approximate scan line relative to the ground point in the current iterative step. Note that Δr is a signed number. Thus, the iterative direction should be correctly determined.
During the iterative process of the ground-to-image transformation, the intermediate back-projected image points may be located on both sides of the linear array. Due to different flight directions and different locations of the ground points relative to the approximate scan line, the determination of the iterative direction is a complicated task. Fortunately, we can apply plane analytic geometry to help determine the iterative direction. Figure 5 shows typical situations of the iterative direction when performing the ground-to-image transformation operation. We use the first scan line as the reference scan line, and the iterative direction can easily be determined from the sign of the expression d 1 × d i . For example, the signs of d 1 × d i shown in Figures 5c and 5d are both positive. According to the fundamental characteristics of ground-to-image transformation, the sign of Δr should be positive in these two situations. The sign of d i in Figure 5c is negative, whereas the sign of d i in Figure 5d is positive. This indicates that the employment of plane analytic geometry simplifies the determination of the iterative direction. Similarly, Figures 5a and 5b indicate that the sign of Δr is negative, although they are two different situations. Thus, the iterative direction can be easily determined with the following expressions.
Some ground-to-image transformation algorithms only use the sign of the single distance d i to determine the iterative direction without a reference scan line. This is lack of robustness and may lead to fluctuations near the best scan line or even divergence.

Convergence Condition
The iterative convergence condition is essential to the solution of the ground-to-image transformation. The correct determination of the iterative direction ensures that the distance d i between the back-projected image

10.1029/2019EA000646
Earth and Space Science point and the linear array decreases. Consequently, the incremental value of the scan line Δr in equation (7) decreases. Therefore, we specify that Δr is less than one scan line as one of the iterative convergence conditions. The other convergence condition is that the sign of the expression d i × d i+1 should be negative, which can be written as (d i × d i+1 ) < 0. The second convergence condition is based on the plane analytic geometry of ground-to-image transformation. More specifically, the expression (d i × d i+1 ) < 0 can fully indicate that the best can line is just located between the two neighbouring scan lines (i.e., i and i +1), such that the iterative process of the ground-to-image transformation should be terminated. Note that there are different iterative directions during the iterative process of ground-to-image transformation; thus, both d i × d i+1 and d i − 1 × d i may require examination, as shown in Figure 6. Therefore, the complete convergence conditions for the proposed method are given as follows It is worth emphasizing that |Δr| < 1 is not used as the only convergence condition because there may be exceptional situations in which |Δr| < 1 cannot ensure that (d i × d i+1 ) < 0 or (d i − 1 × d i ) < 0. These exceptional situations may be caused by rounding errors in numerical calculations or by the approximate nature of predicting the best scan line. In such exceptional situations, one additional iteration is usually required. Although the second item in equation (9) can ensure the successful solution of the ground-to-image transformation operation, we do not use it as the only convergence condition either. This is because the calculation of d i and d i+1 during each iterative loop will double the computation cost. We only need to verify the second convergence condition when the first convergence condition is satisfied. Therefore, these two convergence conditions allow for both efficiency and robustness.
Conventional ground-to-image transformation methods usually use a threshold as the convergence condition. As mentioned in Kim et al. (2001), 0.01 pixels can be specified as the convergence threshold of ground-to-image transformation. This implies that when the displacement of the best scan line determined by two adjacent iterative steps is less than 0.01 pixels, the iterative process terminates. However, such refinement from integer accuracy to 0.01 pixels requires a few more iterations. In contrast, the proposed method uses the plane analytic geometry, namely, to set the convergence conditions of the ground-to-image transformation. This indicates that the iterative process can be terminated when the integer accuracy of ground-to-image transformation is achieved. Consequently, there is no need to carry out additional iterations to refine the ground-to-image transformation accuracy from integer accuracy to subpixel accuracy, which can further decrease the number of iterations. Indeed, the proposed method employs linear interpolation to achieve subpixel accuracy, which will be discussed in the next section. Moreover, linear interpolation always exhibits higher accuracy than the conventional threshold method, which usually delivers a ground-to-image transformation accuracy of better than 1 × 10 −4 pixels (see section 4.3.2). Therefore, the convergence conditions based on plane analytic geometry can reduce the number of iterations and improve the accuracy.

Subpixel interpolation
When the iterative convergence conditions are satisfied, the best scan line with subpixel accuracy is determined by the following equation where i indicates the approximate scan line that is determined when the two convergence conditions are satisfied and d i and d i+1 indicate the generalized distance with respect to two neighbouring scan lines i and i+1.

Employment of an Approximate Value
The approximate value of the best scan line is helpful for the fast convergence of the ground-to-image transformation. In the practical geometric processing procedures such as orthorectification based on an indirect method, a fairly accurate approximate scan line is always available for the ground-to-image transformation operation. Specifically, assuming that orthorectification is performed pixel-by-pixel, thus, the best scan line determined for the current ground-to-image transformation operation can be used as an approximate value for the next ground-to-image transformation operation. Moreover, such an approximate scan line is always very similar to the best scan line, and the deviation is usually only a few scan lines according to the GSD of the generated orthophoto. Accordingly, the computational efficiency of the ground-to-image transformation in the orthorectification procedure is quite high. Note that we have properly solved the problem of the iterative direction, which ensures stable and rapid convergence with the approximate scan line.

10.1029/2019EA000646
Earth and Space Science

Processing Procedures
The RSM for linear pushbroom images should be constructed prior to the real ground-to-image transformation operation. The construction of the RSM requires each scan line's EO parameters and the image coordinates of each detector on the linear array. These data are always supplied in the image support files, such as the DIMAP file for Pleiades images or the orientation data file and camera file for ADS40 images. It is noted that different linear pushbroom images always adopt different support file format. For generality, we use the camera file and orientation data file format developed for ADS40 to process all linear pushbroom images. Specifically, the ADS40 camera file (.cam) stores each detector's image coordinates and the ADS40 orientation data file (.odf) stores each scan line's EO parameters. Note that the RSM only needs to be constructed once. The processing procedures of the proposed ground-to-image transformation method are illustrated in Figure 7. The detailed steps are as follows.
1. Given a ground point P(X,Y,Z), the first scan line is used as the reference scan line to facilitate the determination of the iterative direction. With the EO parameters of the first scan line, the back-projected image point p 1 is determined. Then, the generalized distance between p 1 and the linear array is computed, which is represented as d 1 . 2. Determine whether an approximate scan line is available. If the approximate scan line i is unavailable, it is estimated using the generalized distance d 1 by i = 1+F(d 1 ), and the process moves to step 3. If the approximate scan line is available, move directly to step 3. 3. Compute the back-projected image point p i and the corresponding generalized distance d i using the EO parameters of the approximate scan line i. 4. Determine the iterative direction using the sign of the expression d 1 × d i . If the sign of d 1 × d i is positive, the incremental value Δr is also positive. In contrast, if the sign of d 1 × d i is negative, Δr is negative. These two situations can be expressed by Δr = F(d i ) and Δr = − F(d i ), respectively. Next, update the approximate scan line with i = i+Δr. 5. Determine whether|Δr| is less than one scan line. If |Δr| is larger than one scan line, repeat the iterative steps (3) and (4). If|Δr| is less than one scan line, then the iterative process terminates, and the process moves to the next step. Note that if the number of iterations exceeds the predetermined maximum number of iterations for the ground-to-image transformation algorithm (e.g., 50), then the iterative process fails. 6. Compute the best scan line with subpixel accuracy by linear interpolation using equation (10). 7. Compute the image coordinates p(x,y) corresponding to the ground point by applying the extended collinearity equation. 8. Compute the pixel coordinates p(s,l) corresponding to the ground point using image-to-sensor transformation.

Orthorectification
This section investigates the applications of the proposed method in the geometric processing procedures. As we all know, orthophotos can be generated using both indirect method and direct method. Compared with direct method, indirect method shows advantage in image resampling. In the orthorectification procedure based on indirect method, given a pixel on the generated orthophoto, the corresponding pixel on the raw image should be determined using ground-to-image transformation. Then, the grey value of the pixel on the raw image is specified as the grey value of the corresponding pixel on the orthophoto. Thus, the ground-to-image transformation is a key step in the orthorectification procedure. Such ground-to-image transformation operation is efficient and easy to implement for frame images. In contrast, it is relatively complicated for linear pushbroom images due to iterative calculations.
For linear pushbroom images, the orthorectification process requires a series of coordinate transformation steps, as shown in Figure 8. Given a pixel (s o , l o ) on the orthophoto that will be generated, the corresponding map projection coordinates (X o , Y o ) can be easily calculated with the GSD and the starting point's coordinates on the orthophoto. The corresponding geodetic coordinates (B,L) are then computed through inverse map projection transformation. In addition, the ellipsoidal height H with respect to (X o , Y o ) is calculated by interpolation using the reference DTM data. Next, the 3-D geodetic coordinates (B,L,H) are converted to ground coordinates (X,Y,Z), and the proposed ground-to-image transformation method is used to map the 3-D ground coordinates (X,Y,Z) to 2-D image coordinates (x,y). Finally, the pixel coordinates (s,l) on the raw image are calculated with the image coordinates (x,y). Of these coordinate transformation steps, the ground-to-image transformation is the most time-consuming step. Suppose that orthophotos are generated pixel-by-pixel, there are hundreds of millions of ground-to-image transformation operations when generating an orthophoto from a satellite image (e.g., Pleiades). Obviously, improving the computational

10.1029/2019EA000646
Earth and Space Science efficiency of ground-to-image transformation will enable the rapid acquisition of orthophotos from raw linear pushbroom images.

Image Matching Based on Orthophotos
Image matching of linear pushbroom images has been a popular area of research for many years (Hirschmüller, 2008;Tao et al., 2018). The epipolar line geometry for linear pushbroom images has been extensively investigated (Habib et al., 2005;Morgan et al., 2006). However, a robust and widely accepted epipolar resampling method for linear pushbroom images is still unavailable, due to the various types of linear array cameras and the resultant complicated epipolar line geometry. Furthermore, although numerous image matching algorithms have been proposed, the problem of image matching has still not been completely solved. Recently, the methodology of image matching on orthophotos is being increasingly widely used (Heipke et al., 2007;Speyerer et al., 2016). The methodology of image matching on orthophotos has the benefits of (1) removing image distortions caused by terrain relief and perspective projection, (2) having identical GSD values for stereo pairs, and (3) providing a good starting position for searching conjugate points.
The methodology of image matching on orthophotos is especially suitable for areas without artificial objects, where the terrain exhibits better continuity. Due to the fact that conjugate points are mainly used for bundle adjustment or DTM generation, the matched conjugate points on orthophotos need to be converted to conjugate points on the original images. This coordinate transformation process is indeed a ground-to-image transformation operation, which is analogous to that in the orthorectification procedure.

Software and Test Data Sets
The proposed ground-to-image transformation method was implemented based on the open source photogrammetric software-Dirk's General Analytical Plotter (DGAP, 2019). DGAP provides many practical photogrammetric functions and reduces the burden of software development. The orthorectification and image matching software modules were developed by us independently using Qt 5.10.1. In case of the developed orthorectification software module, the proposed ground-to-image transformation method was used to determine the pixel coordinates on the original image. The developed image matching software module supports the methodology of image matching on the orthophotos, which uses the proposed ground-to-image transformation method to convert conjugate points from orthophotos to original images. The hardware configurations of the experimental environment were an Intel Core i7-7500U with a 2.70-GHz CPU and 8-GB RAM. The operating systems include Windows 10 and Ubuntu 14.04, where Ubuntu 14.04 was installed on a virtual machine running under Windows 10. The Ubuntu 14.04 operating system was mainly used to process planetary images. Due to the cross-platform characteristics of Qt, the developed software modules can be easily run in both Windows and Ubuntu operating systems. However, it should be noted that in the virtual machine environment (i.e., Ubuntu 14.04), the computing performance will be affected to some extent. The maximum number of iterations for the ground-to-image transformation is set to 50. The test data include Earth observation satellite images, planetary orbital photographs, and airborne images. These test images were all acquired by linear array cameras. Basic information about the test data is reported in  The Earth observation satellite image data are "Tristereo" mode Pleiades 1A sample images, which can be downloaded from https://www.intelligence-airbusds.com/en/8262-sample-imagery. These Pleiades images are located in Melbourne, Australia. The raw focal plane of Pleiades images is composed of five slightly angled short linear arrays, and each linear array contains 6,000 pixels. To derive user-friendly satellite image products, a virtual sensor technique is employed to simulate the imaging geometry of a single pushbroom linear array. Thus, the generated satellite image product can be treated as that acquired by a virtual linear array camera with about 30,000 pixels. The raw Pleiades satellite images have a spatial resolution of 0.7 m and were improved to 0.5 m through special resampling technique developed and validated by the French Space Agency (CNES; Pleiades, 2019). Consequently, due to image resampling, the image size of Pleiades images varies not only in the line direction but also in the sample direction, as reported in the "Image size" column in  (Robinson et al., 2010). The LRO NAC consists of two identical cameras, namely, NAC-LEFT (NAC-L) and NAC-RIGHT (NAC-R) that are mounted side by side with an overlap of dozens of pixels to increase the coverage. MEX HRSC is a specially designed Mars mapping camera that has nine linear arrays (five panchromatic and four multispectral) mounted parallel on the focal plane (Albertz et al., 2005). MEX HRSC can acquire stereo images with one observation. The experiment used two panchromatic channels (i.e., S1 and S2) of MEX HRSC data. In contrast to LRO NAC and MEX HRSC, Chandrayaan-1 M3 is an imaging spectrometer that is mainly used to map the Moon's mineral Note. In the "Images" column, the suffixes L and R indicate LRO NAC-LEFT and NAC-RIGHT, respectively, and the prefix M3G indicates Chandrayaan-1 M3 images acquired with the low-resolution global mode. Abbreviation: GSD, ground sample distance.

10.1029/2019EA000646
Earth and Space Science resources (Green et al., 2011). Chandrayaan-1 M3 works with two modes: low-resolution global mode and high-resolution target mode. The experiment used low-resolution global mode data, and the resultant hyperspectral images contain 85 bands. The RSM for the planetary exploration linear array cameras was constructed with the Spacecraft Planetary Instrument C-matrix Events (SPICE) Toolkit (Acton et al., 2016). Specifically, each detector's image coordinates and each scan line's EO parameters were computed and then converted to the ADS40 camera file and orientation data file, respectively. The data preprocessing of the planetary images (e.g., raw images ingestion and sensor model initialization) were carried out using Integrated System for Imagers and Spectrometers (ISIS 3.5.2; Edmundson et al., 2012) that installed in the Ubuntu 14.04 operating system. Then, the generated ISIS cube format image files, camera files, and orientation data files are ingested into the developed orthorectification and image matching software modules to conduct experiment.
ADS40 is a typical commercial airborne linear array camera, which can acquire forward looking, backward looking, and nadir images simultaneously (Gonzalez et al., 2013). Due to that the airborne platform is susceptible to turbulence, the scan lines of ADS40 images may change more dramatically compared with satellite images (e.g., MEX HRSC). Therefore, the ADS40 images can provide a more comprehensive evaluation for the proposed ground-to-image transformation method. The ADS40 test data are the example data of the open source photogrammetric software DGAP, which were taken over the famous Vaihingen/Enz airborne test site located in Stuttgart, Germany. We used two flight lines (i.e., 0931 and 0937) to conduct the test. Flight line 0931 is oriented north-south, whereas flight line 0937 is oriented west-east. Different flight line directions are beneficial for verifying the robustness of the proposed method. For each flight line, three panchromatic channels with different viewing angles, namely, forward 28 degrees (PANF28), backward 14 degrees (PANB14), and nadir 0 degrees (PANN00), were used. At the data preprocessing stage, the Global Positioning System (GPS) and Inertial Measurement Unit (IMU) data that were acquired synchronously with the ADS40 images are solved to generate orientation data files. Then, the provided camera files and the derived orientation data files are used to construct the RSM.

Image Distortions
The image coordinates of the detectors on the linear array cameras were plotted on the focal plane based on the constructed sensor model for each camera. The results, shown in Figure 9, can be used to assess the image distortions of the test data. For better visual effects, the x axis of the figure is scaled.
As presented in Figures 9a-9c, in our camera models for Pleiades and MEX HRSC, we assume a straight line shape linear array due to the fact that the image distortions are not considered. It should be noted that MEX HRSC does contain distortions, but the results of the camera calibration were not included in the SPICE kernels as well as the subsequent geometric processing procedures. The PANB14A and PANF28A channels of ADS40 exhibit more distortions than the PANN00A channel, as shown in Figures 9g-9i. The maximum image distortions for the test data are reported in Table 2. For the LRO NAC and Chandrayaan-1 M3 data, the corresponding image distortions were calculated with the authoritative camera calibration parameters provided in the SPICE kernels. For the ADS40 images, a camera calibration file containing each detector's image coordinates is provided. Therefore, we form a line using the ends of the detectors on the linear array of ADS40. Then, the distance between each detector and the line was computed, and the maximum distance can be regarded as the maximum distortion for the ADS40 test data. The last column in Table 2 indicates the number of the divided line segments required to apply the CPP-based method. Here, a segmentation threshold of 0.5 pixels is used to segment the linear array. In the LRO NAC-L and NAC-R test data, the image distortions are only in the y axis direction. Though the maximum distortions for the two LRO NACs are 14.59 and 15.08 pixels, respectively, the corresponding linear arrays still show a straight line shape, as shown in Figures 9d and  9e. Therefore, the CPP of a scan line is still a 3-D plane. Thus, the CPP-based method can be directly applied for LRO NAC test data, without the need to segment the linear array. The maximum image distortions for the Chandrayaan-1 M3 data were 10.38 pixels, such that the linear array should be segmented into four line segments to apply the CPP-based method. For ADS40 images, the maximum image distortions are 0.15, 30.84, and 81.30 pixels for PANN00A, PANF28A, and PANB14A, respectively. Thus, to apply the CPP-based method, the corresponding numbers of divided line segments are 1, 28, and 19, respectively.

Ground-to-Image Transformation Results
We evaluated the computational efficiency and accuracy of the proposed ground-to-image transformation method using one million ground point-image point correspondences with known coordinates. These

10.1029/2019EA000646
Earth and Space Science image points are grid dots spaced by one pixel. The grid size for the test data is 1,000 × 1,000 except for Chandrayaan-1 M3 test data. Because the linear array of M3 only consists of 304 pixels for the lowresolution global mode, the grid size for the M3 test data is 5,000 rows and 200 columns. The corresponding ground coordinates were calculated by projecting the image points onto a reference DTM using the extended collinearity equation. Thus, both the pixel coordinates and ground coordinates for the one million ground point-image point correspondences are known prior to conducting the ground-toimage transformation experiment. Here, the known pixel coordinates can be considered as true values.

Computational Efficiency Evaluation
To evaluate the computational efficiency, we compared the proposed method with two representative methods, namely, the cutting-edge CPP-based method and the classical bisection method. It is noted that the CPPbased method delivers high computational efficiency for the linear pushbroom images without distortions   (Geng et al., 2019;Wang et al., 2009), whereas the bisection method is widely used but shows relatively lower efficiency. The results are shown in Figure 10. As can be seen, the bisection method provides much lower computational efficiency than the proposed GDP-based method and the CPP-based method. Specifically, the bisection method required 6-11 s to accomplish the ground-to-image transformation operations of one million ground points. This will result in a very slow processing speed for deriving mapping products from linear pushbroom images.
The two efficient methods behave differently. The proposed GDPbased method provides a relatively stable computational efficiency for both linear pushbroom images with distortions and linear pushbroom images without distortions. In contrast, the computational efficiency of the CPP-based method is affected by image distortions: (1) the processing times for the Chandrayaan-1 M3 data and the ADS40 data are much greater than those for the Pleiades, LRO NAC, and MEX HRSC data and (2) the ADS40 PANN00A test data require less processing time than the PANB14A and PANF28A test data. Note that the PANN00A image has smaller image distortions than the PANB14A and PANF28A images as shown in Figure 9.

Accuracy Evaluation
To evaluate the accuracy, comparisons were made between the proposed GDP-based method and the CPP-based method. Here, the classical bisection method was not used for the accuracy comparison because the bisection method also uses interpolation in the image space to achieve subpixel accuracy, such that the obtained ground-to-image transformation accuracy should be the same as that of the proposed GDP-based method. With the known ground coordinates, the corresponding pixel coordinates can be calculated using the proposed GDP-based method and the CPP-based method, respectively. Then, the displacements between the calculated and known pixel coordinates were computed. The maximum errors and root-mean-square errors (RMSEs) for these pixel coordinate displacements were calculated. The experimental results of the accuracy evaluation are shown in Figure 11.
To better display the results, logarithmic charts are used in Figure 11. It is noted that the proposed GDPbased method achieves higher accuracy than the CPP-based method for the linear pushbroom images whose linear array shows straight line shape (i.e., Pleiades, MEX HRSC, and LRO NAC images). This can be explained by the fact that the CPP-based method achieves subpixel accuracy by interpolation in the object space. In contrast, the proposed GDP-based method achieves subpixel accuracy by interpolation in the image space. Thus, the interpolation accuracy of the former is lower than that of the latter. When the linear pushbroom images whose linear array shows curved shape (i.e., Chandrayaan-1 M3 and ADS40 images) are processed with the CPP-based method, we used the original CPP-based method at the initial iterative process. But at the final subpixel refinement stage we used linear interpolation in the image space to achieve subpixel accuracy, which is analogous to that of the proposed GDP-based method. Consequently, the proposed GDPbased method and the CPP-based method have the same ground-to-image transformation accuracy for the four Chandrayaan-1 M3 images and the six ADS40 images.

Orthorectification Results
In this section, the proposed ground-to-image transformation method is examined with a practical orthorectification application. Unfortunately, the ADS40 raw image data are unavailable. Thus, we used the remaining linear array cameras for the orthorectification experiment. For each linear array camera, two images were used to derive orthophotos. To better test the performance of the proposed method, we performed pixel-level orthorectification based on indirect method. The datum used in the experiment is defined as follows. In case of the Earth observation remote sensing images (i.e., Pleiades), the WGS84 coordinate system and Universal Transverse Mercator (UTM) map projection were adopted. In case of the planetary images, the radius of the reference spheres for lunar images and Martian images were 1,737.4 and 3,396.19 km, respectively, and a simple equirectangular map projection is adopted. The orthorectification results of the test data are reported in Table 3. The data type of the planetary images is 32 bits due to that the raw Planetary Data System (PDS) format images were converted to ISIS cube format images. It should be noted that Table 3 lists the orthorectification processing time based on Windows 10 operating system with single-threaded (ST) mode. For planetary images, the derived camera files and orientation data files are copied from Ubuntu operating system. In addition, for a fair computational efficiency comparison, the image rectification processing time for the first band of the Chandrayaan-1 M3 hyperspectral images was listed in Table 3.

Computational Performance Optimization of Orthorectification
In order to further improve the computational efficiency of generating orthophotos, the developed orthorectification software module was optimized using multi-threaded (MT) programming technique based on Qt. To compare computational efficiency, we also used available commercial photogrammetric software ERDAS IMAGINE 2013 and open source software ISIS 3.5.2 to generate orthophotos. Specifically, ERDAS IMAGINE Figure 11. Comparisons between the proposed generalized distance prediction (GDP)-based method and the central perspective plane (CPP)-based method for ground-to-image transformation accuracy. Max. indicates maximum, and RMSE indicates root-mean-square error.

Earth and Space Science
2013 was used to process the Pleiades images, and ISIS 3.5.2 was used to process the planetary images. It is noted that ERDAS IMAGINE 2013 runs with MT mode in Windows operating system, and ISIS 3.5.2 runs with ST mode in Linux operating system (i.e., Ubuntu 14.04 in this paper). However, it is a pity that ERDAS IMAGINE 2013 fails to support the RSM of Pleiades images, such that only RFM method was used.
But it also provides an interesting comparison between RSM and RFM. In case of the Chandrayaan-1 M3 images, ISIS conducts image rectification for all bands at once. Thus, the developed orthorectification software module also processed all bands of the Chandrayaan-1 M3 images in this section. However, we did not conduct MT optimization for Chandrayaan-1 M3 images. Instead, the orthorectification for Chandrayaan-1 M3 images was optimized by storing the ground-to-image transformation information of the first band. This can significantly improve the computational efficiency of orthorectification for hyperspectral images as well, since different bands share the same ground-to-image transformation information. The computational efficiency evaluation results are shown in Figure 12. In case of the MT applications, there were four worker threads in the experiment.
As can be observed in Figure 12, the computational efficiency derived from Ubuntu 14.04 operating system was lower than that derived from Windows 10 operating system. This is due to that the Ubuntu operating system was installed on a virtual machine. To focus on the ground-to-image transformation algorithm, the computational efficiency comparison should be conducted based on the same software and hardware environments. Specifically, in case of Pleiades images, the comparison between ERDAS IMAGINE 2013 and the developed orthorectification software should be conducted in Windows 10 operating system with MT mode. In case of planetary images, the comparison between ISIS 3.5.2 and the developed orthorectification software should be carried out in Ubuntu 14.04 operating system with ST mode.

Statistical Analysis of the Number of Iterations
During the orthorectification process, the number of iterations for each ground-to-image transformation operation was also recorded. The statistical analysis results are presented in Figure 13. Note that the ground-to-image transformation operation usually requires several iterations in the initial orthorectification stage because the approximate scan line is unavailable. However, the number of these points is very small, which is usually less than 0.01% of the total number of pixels in the generated orthophoto. Thus, they are not shown in Figure 13 because of the limited display accuracy. It is observed that most of the ground-to- image transformation operations in orthorectification can be solved with one or two iterations. As discussed previously, when the number of iterations reaches the maximum threshold (e.g., 50), the ground-to-image transformation has failed. Too many invalid pixels on the generated orthophotos will decrease the computational efficiency of the orthorectification.

Geometric Accuracy Evaluation
We used the H4165 S1 image to assess the geometric accuracy of the generated orthophotos. The orthophoto of H4165 S1 image was generated using both our method and ISIS, respectively, and the orthophoto derived from ISIS was used as reference image. We manually measured 10 conjugate points on both orthophotos. Then, the geographic coordinate offsets of these conjugate points were calculated. The geometric accuracy evaluation results are shown in Figure 14. As can be seen, Figure 14a shows the orthophoto derived from our method and the conjugate point distribution, Figures 14b and 14c show the location of the first conjugate points at the original resolution, Figure 14d presents the mosaic images formed by the orthophotos derived from our method and ISIS, and Figure 14e presents the geographic coordinates offsets between the orthophoto derived from our method and ISIS. It can be observed that there were only very small geometric offsets between the two orthophotos. The RMSE values in X and Y directions were 0.81 and 0.33 pixels, respectively.

Image Matching Results
For comparison, image matching experiments using normalized cross-correlation (NCC) method were conducted on both stereo orthophotos and stereo original images. In case of the methodology of image matching on orthophotos, image matching is first conducted on orthophotos. Then, the matched conjugate points on the orthophotos are converted to conjugate points on the original images using the proposed ground-toimage transformation method. Thus, the conjugate points on the original images derived from different methods can be compared. We used Pleiades (PHR0025329 and PHR0026276) and MEX HRSC (H7242 S1/ S2) data for the experiment. The PHR0025329 and H7242 S1 images were used as reference images to extract the initial key points using the Shi-Tomasi algorithm (Shi & Tomasi, 1994). For the image matching on the orthophotos, the starting positions for searching the conjugate points were predicted using the geodetic coordinates of the reference key points on the orthophotos. For the image matching on the original images, the starting positions for searching the conjugate points were determined by projecting the reference key points on the orthophotos onto the original images. The image matching parameters and results for the test data are listed in Table 4.  Earth and Space Science matched conjugate points on the orthophotos were accurately transformed onto the original images. This also verified the geometric accuracy of the proposed ground-to-image transformation method.

Ground-to-Image Transformation Efficiency
The computational efficiency of the ground-to-image transformation algorithm is crucial in the practical geometric processing procedures such as orthorectification and DTM generation. The proposed GDP-based method and the CPP-based method required 0.659-1.413 s and 0.356-2.605 s, respectively, to accomplish the ground-to-image transformation operations of one million ground points. As shown in Figure 10, the CPP-based method has better computational performance than the proposed GDP-based method for the linear pushbroom images whose linear array shows straight line shape (i.e., Pleiades, MEX HRSC, and LRO NAC images). This is because the CPP-based method can be directly employed to process these linear pushbroom images without the requirement of segmenting the linear array. We acknowledge the obvious advantages of the CPP-based method for the linear pushbroom images whose linear array shows straight line shape, and we do not intend to challenge the CPP-based method in this respect. The results shown in Figure 10 clearly indicate that the computational efficiency of the proposed GDP-based method is better than that of the CPP-based method for the linear pushbroom images whose linear array shows curved shape. More specifically, for the Chandrayaan-1 M3 and ADS40 test data, the performance improvements range from 30 to 50% compared with the CPP-based method. These performance improvements are due to (1) the use of plane analytic geometry to determine the iterative direction and making good use of the approximate scan

Ground-to-Image Transformation Accuracy
The experimental results shown in Figure 11 demonstrate that the proposed GDP-based method can provide very accurate ground-to-image transformation. Overall, the ground-to-image transformation accuracy of the proposed method is dependent of specific image and ranges from 1 × 10 −4 to 1 × 10 −9 pixels for the l-component and 1 × 10 −5 to 1 × 10 −10 pixels for the s component. However, Figure 11 also shows several trends. In terms of the proposed GDP-based method, the maximum errors and RMSE values for the linear pushbroom images whose linear array shows curved shape are slightly larger than those for the linear pushbroom images whose linear array shows straight line shape. Specifically, for the Chandrayaan-1 M3 and ADS40 images, the maximum errors with the proposed method are in the ranges of 1 × 10 −5 to 1 × 10 −8 pixels for the s component and 1 × 10 −4 to 1 × 10 −7 pixels for the l component. For the Pleiades, MEX HRSC, and LRO NAC images, the maximum errors with the proposed method are in the ranges of 1 × 10 −7 to 1 × 10 −10 pixels for the s component and 1 × 10 −5 to 1 × 10 −9 pixels for the l component. Note that the Chandrayaan-1 M3 and ADS40 images have more prominent image distortions and curved shape linear arrays. This implies that image distortions affect the proposed method and degrade the ground-to-image transformation accuracy.
In Figure 11, the maximum errors of the l component indicate the accuracy of the determined best scan line corresponding to the ground point. The determination of the best scan line is the key to the ground-to-image transformation of linear pushbroom images. As shown by the results, both the CPP-based method and the proposed GDP-based method can determine the best scan line with an accuracy of better than 1 × 10 −4 pixels, which is a fairly high accuracy. In addition, the accuracy of the s component is slightly better than that of the l component. This is because the best scan line (i.e., l component) is determined using the complicated collinearity equation, whereas the s component is determined using very simple numerical calculations. It is noteworthy that the accuracy discussed here is indeed the internal coincidence accuracy of the constructed sensor model instead of the geometric accuracy of the imaging instrument. More specifically, the MEX HRSC does have image distortions, but they are not introduced into the geometric processing procedures. Under this circumstance, inaccurate sensor model (i.e., lacking camera calibration parameters) will not influence the accuracy of the ground-to-image transformation method.

Application of the Proposed Method in Orthorectification
The last column in Table 3 (i.e., "million points per second") indicates the computational efficiency of the developed orthorectification software module in a more understandable way. As can be seen, the computational efficiency of the orthorectification reaches about one million points per second. These are satisfactory results for orthorectification of linear pushbroom images based on the RSM. Figures 12b and 12c indicate that the developed orthorectification software module running under Ubuntu 14.04 operating system with ST mode was about three to six times faster than ISIS 3.5.2. This can clearly show the advantage of the proposed ground-to-image transformation method. Moreover, the computational performance of orthorectification was further improved through MT optimization. Specifically, the DOMs derived from MEX HRSC and LRO NAC images can be generated in about half a minute and 2 min, respectively, when the developed orthorectification software module runs under Windows 10 operating system with MT mode. Note that the data volume of the orthophotos derived from LRO NAC images is close to 2 GB in size. Such high computational efficiency of orthorectification is satisfactory. In case of the Chandrayaan-1 M3 images, as shown in Figure 12d, the computational efficiency has increased by 10 to 20 times, compared with that of ISIS 3.5.2. However, this was mainly due to the optimization of storing the ground-to-image transformation information of the first band and partly attributed to the proposed ground-to-image transformation method.
As shown in Figure 12a, the commercial photogrammetric software (i.e., ERDAS IMAGINE 2013, 2019, https://download.hexagongeospatial.com/downloads/imagine/erdas-imagine-2013-includes-lps (accessed 07-September-2019)) delivered the highest computational efficiency for Pleiades images, and the corresponding orthophotos can be generated in about three minutes, namely, 184.0 and 186.0 s for PHR0025329 and PHR0026276 images, respectively. However, this high computational efficiency was obtained due to the RFM method. For comparison, the developed orthorectification software module running under Windows 10 operating system with MT mode can accomplish such image rectification in about 5 min, namely, 292.0 and 283.9 s for PHR0025329 and PHR0026276 images, respectively. As we all know, it is generally believed that RFM is much more efficient than RSM. The comparative experimental results demonstrate that although the RSM method is still lower than the RFM method, the computational performance difference is not very large. Since the proposed ground-to-image transformation method plays a significant role in orthorectification.
As shown in Figure 13, it indicates that the proposed ground-to-image transformation method can make full use of the approximate scan line, which allows most of the ground-to-image transformation operations in the orthorectification procedure to be solved in one or two iterations. This further explains the high efficiency of orthorectification using the proposed ground-to-image transformation method.

Application of the Proposed Method in Image Matching
As shown in Table 4, the methodology of image matching on orthophotos can deliver more conjugate points than the conventional methodology of image matching on original images. Specifically, the success rate of the image matching was improved by 21.20% and 4.65% for MEX HRSC and Pleiades data, respectively. Moreover, the proposed ground-to-image transformation method can accurately transform the matched conjugate points on the orthophotos to conjugate points on the original images. In addition, the MEX HRSC test gains more improvement in the success rate of image matching than the Pleiades test. This is because the original MEX HRSC stereo images have large differences due to the varying exposure time and the very different viewing angles. Additionally, it is noted that the same crater shows different shapes in Figures 15b and  15d. Obviously, image distortions on the original images will increase the burden of image matching. In contrast, stereo orthophotos have the same GSD and removed image distortions, so it is easy to perform image matching.
In addition, matched conjugate points inevitably contain errors due to repetitive patterns, shadows, occlusions, and moving objects (Gruen, 2012). As shown in Figure 16b, dozens of small rectangular objects with sizes of only a few pixels were correctly matched with the methodology of image matching on orthophotos.
In contrast, the conventional methodology of image matching on the original images, illustrated in Figure 16d, shows matching errors. This indicates that the methodology of image matching on orthophotos is more robust to repetitive patterns than the conventional methodology of image matching on original images. This is because the image distortions are removed by the orthorectification, and a small search window can be used to decrease the number of false matching points.
In summary, the experimental results demonstrate the effectiveness of the image matching methodology on orthophotos. Moreover, it can be inferred that the methodology of image matching on orthophotos can be used as a framework for developing various image matching algorithms. This will greatly enhance the performance of image matching for linear pushbroom images. In case of the methodology of image matching on orthophotos, it is noted that both orthophotos generation and pixel coordinates transformation between original images and orthophotos can use the proposed ground-to-image transformation method.

Algorithm Comparisons
In terms of efficiency, both the proposed GDP-based method and the CPP-based method are superior to the classical bisection method (see Figure 10). The proposed GDP-based method is the best for the linear pushbroom images whose linear array shows curved shape (i.e., Chandrayaan-1 M3 and ADS40 images). Whereas the CPP-based method has the highest efficiency for the linear pushbroom images whose linear array shows straight line shape (i.e., Pleiades, MEX HRSC, and LRO NAC images). The CPP-based method uses the distance between a ground point and the CPP of the scan line to estimate the iterative incremental value of the scan line, which has the advantage of reducing complicated collinearity equation computations. However, when the linear array has a curved shape (e.g., ADS40) due to distortions, the CPP of the scan line becomes a 3-D curved surface instead of a 3-D plane. Thus, additional processing procedures should be conducted to address the problem of the curved shape linear array. This makes the CPP-based method relatively complicated for the linear pushbroom images whose linear array shows curved shape. In contrast, the proposed GDP-based method uses the generalized distance between the back-projected image point and the curved shape linear array to predict the best scan line corresponding to the ground point, which has advantages in processing the linear pushbroom images with distortions in principle. Additionally, the proposed GDPbased method makes full use of plane analytic geometry to solve the practical problems of ground-to-image transformation, resulting in rapid convergence of the iterative computations. Currently, small satellites equipped with low-cost cameras have become popular, deriving massive remote sensing images with large distortions. Under this circumstance, the advantage of the proposed GDP-based method is more obvious. The low efficiency of the bisection method is due to the high complexity of the binary search algorithm. It can be concluded that even though the bisection method is widely used due to its convenience, more efficient algorithms should be used for efficiency-oriented applications.
In terms of accuracy, the proposed GDP-based method is slightly better than the CPP-based method for the linear pushbroom images whose linear array shows straight line shape. It is noted that the RFM method delivers a fit accuracy of approximately 0.05 pixels for satellite images such as IKONOS (Dial et al., 2003). Thus, both the CPP-based method and the proposed GDP-based method can offer much higher accuracy than the RFM method, which demonstrates the advantage of the RSM.
Considering algorithm implementation, the proposed GDP-based method is simpler than the CPP-based method especially for the linear pushbroom images with distortions. In case of the proposed GDP-based method, there is almost no difference in algorithm implementation between the linear pushbroom images with distortions and the linear pushbroom images without distortions. In contrast, the CPP-based method is easy to implement for the linear pushbroom images without distortions, but more work should be carried out to apply it to the practical linear pushbroom images with varying degrees of distortions, including linear array segmentation and determining the exact CPP of the scan line. Moreover, because the CPP-based method requires memory to store the 3-D plane equation parameters of each scan line, memory consumption will be a problem for long strip images. Specifically, for ADS40 images, hundreds of thousands of scan lines are typically acquired during a flight line. Because the curved shape linear array should be segmented into multiple segments to apply the CPP-based method (Wang et al., 2009), the required memory consumption for long strip image data will be on the order of dozens of MB. Such high memory consumption will limit the applications of the CPP-based method in in-orbit real-time processing scenarios. In contrast, to implement the proposed GDP-based method, we only need to store the line equation parameters of a series of tiny line segments formed by the neighboring detectors of the linear array. Thus, the required memory consumption of the proposed method is always much less than 1 MB.

Conclusions
A robust and efficient ground-to-image transformation method plays a significant role in the geometric processing procedures of linear pushbroom images. The fact that the linear array may exhibit a curved shape due to image distortions makes the ground-to-image transformation more complicated and causes some algorithms in the literature to not be applicable in practical scenarios. A novel ground-to-image transformation algorithm for linear pushbroom images based on the concept of generalized distance prediction is proposed. Based on the extensive experimental results, we conclude that (1) the proposed ground-to-image transformation method is suitable for various types of linear pushbroom images and has advantages of easy implementation, high accuracy, and efficiency; (2) the application of the proposed method in orthorectification provides satisfactory computational efficiency; and (3) the proposed method can be used for image matching of linear pushbroom images with the methodology of image matching on orthophotos.
The advantages of the proposed method are mainly attributed to the fact that the algorithm is designed based on the practical rules derived from the fundamental characteristics of ground-to-image transformation, and we take advantage of plane analytic geometry to solve the practical problems of ground-to-image transformation. This study demonstrates that ground-to-image transformation based on the RSM can offer satisfactory computational efficiency. This high computational efficiency, although still lower than that of the RFM, will make the ground-to-image transformation based on the RSM gain more practical applications. The old concept that ground-to-image transformation based on the RSM exhibits low efficiency should be updated.
Currently, timely and efficient mass data processing is still the bottleneck to extensive applications of remote sensing images. Linear array cameras are commonly used imaging instruments, which require more complicated geometric processing procedures. The proposed ground-to-image transformation method will significantly enhance the geometric processing abilities for various types of linear pushbroom images. Due to the large amount of linear array camera data, graphics processing unit acceleration techniques should be adopted to achieve product-level processing efficiency in practical applications scenarios (e.g., orthophotos and DTM generation). The proposed method is easy to implement in parallel, which will be studied in more detail in the future.