Epipolar Resampling of Cross-Track Pushbroom Satellite Imagery Using the Rigorous Sensor Model

Epipolar resampling aims to eliminate the vertical parallax of stereo images. Due to the dynamic nature of the exterior orientation parameters of linear pushbroom satellite imagery and the complexity of reconstructing the epipolar geometry using rigorous sensor models, so far, no epipolar resampling approach has been proposed based on these models. In this paper for the first time it is shown that the orientation of the instantaneous baseline (IB) of conjugate image points (CIPs) in the linear pushbroom satellite imagery can be modeled with high precision in terms of the rows- and the columns-number of CIPs. Taking advantage of this feature, a novel approach is then presented for epipolar resampling of cross-track linear pushbroom satellite imagery. The proposed method is based on the rigorous sensor model. As the instantaneous position of sensors remains fixed, the digital elevation model of the area of interest is not required in the resampling process. Experimental results obtained from two pairs of SPOT and one pair of RapidEye stereo imagery with different terrain conditions shows that the proposed epipolar resampling approach benefits from a superior accuracy, as the remained vertical parallaxes of all CIPs in the normalized images are close to zero.


Introduction
The parallax values of the conjugate image points (CIPs) along the baseline of stereo images (horizontal parallax) reveal valuable information about the height of objects; while the vertical parallaxes not only do not convey any useful information, but they can interrupt the stereo-viewing process. The main objective of the epipolar resampling of stereo images is to generate normalized images where the CIPs have no vertical parallaxes and their horizontal parallaxes are linearly proportional with the height of corresponding points in the object space. These are two main applications of epipolar geometry (EG) which provide the possibility of stereo-viewing, automatic image matching, digital elevation model (DEM) generation, and stereo measurements [1,2].
The EG of frame images is well-known and there are well-established procedures for epipolar resampling [3,4] in such a way that given the relative orientation parameters (ROP) of the stereo images and modifying their attitude parameters, vertical parallaxes of CIPs can be eliminated. To this end, the ROP of the stereo imagery are modified so that all conjugate epipolar lines in those stereo images be simultaneously parallelized with the base line of imagery in a common plane [3]. In contrast, linear pushbroom images have a more complicated geometry where each scan line of these images has its own exterior orientation parameters (EOP) [5]. The availability of the GPS/INS information can positively affect the georeferencing of linear pushbroom imagery, however, the dynamic nature of the EOP of these images makes the reconstruction of EG much more complicated [5,6]. The EG of satellite pushbroom imagery was investigated in various studies using the both rigorous and empirical sensor models [6][7][8][9][10][11][12][13][14][15].
Among the rigorous model-based studies, Gupta and Hartley [7], Kim [5], and Habib et al. [8] made efforts to investigate the EG of linear array images by applying the Multiple Projection Centers (MPC) model [9]. In these studies, the trajectory of the sensor was assumed to be straight (i.e., the simplest case). The only difference between the model employed by Kim [5] with the models used in the two other studies is in the way that the sensor's orientation parameters were modeled. However, these three studies showed two common important results: (1) rigorously derived epipolar lines in scenes captured by linear array scanners are not straight, rather, their shape is hyperbola-like; and (2) each point on a rigorously defined epipolar curve has a specific conjugate epipolar curve (i.e., no conjugate epipolar curves exist). In another study, Lee and Park [10] derived the equation of epipolar lines using the simplified pushbroom sensor model, which is in line with the results of the previous studies. As an important note, no epipolar resampling method has been proposed for linear pushbroom scenes based on the rigorous sensor models till now. Among the empirical model-based studies, Morgan et al. [11] used the Parallel Projection (PP) model for epipolar resampling of linear pushbroom imagery, while Oh et al. [12], Wang et al. [13], and Koh and Yang [14] employed the Rational Function Model (RFM). In the former study, for a better compliance of the imaging geometry with the PP model, the perspective projection of the cross-track imaging component of original images should firstly be transformed to the parallel one [15]. Such a transformation assumes a flat terrain or a DEM and also assumes a knowledge of the scanner roll angle [11,[16][17][18]. Thereafter, transformation parameters from the original to the normalized images can be estimated using a PP model based on a set of CIPs. The three latter studies are based on two features of epipolar curves, namely straightness approximation and local conjugacy [13,19]. Oh et al. [12] determined the direction of approximately conjugate epipolar lines (ACEL) in the stereo imagery space. Because the sampling distance is supposed to be fixed in the image space, the horizontal parallaxes will be linearly proportional to the ground heights, but there is no guarantee the ground axes will be orthogonal [14]. In contrast, Wang et al. [13] determined the direction of ACEL on a virtual horizontal plane in the object space. Because the sampling distance is supposed to be fixed in the object space, the ground axes will be orthogonal, but the horizontal parallaxes may not be proportional to the ground heights [14]. Koh and Yang [14] determined the direction of ACEL in the image space and supposed the sampling distance to be fixed in the ground space. Therefore, this model benefits from the advantages of the two previous models.
In comparison with the Morgan et al. [11] model, the main advantage of RFM-based models is that there is no need to any CIPs to use these models. However, the lack of any physical or geometrical interpretation of the rational function coefficients (RPCs) hinders theoretical investigation of the nature of the EG. Besides, directly determined RPCs are subject to biases in the scanner's IOP or EOP [20][21][22] and many ground control points (GCPs) are required to indirectly determine these parameters [23], which therefore limits the use of rational functions. In contrast, the PP parameters have a geometrical interpretation, and therefore, this model can explain the EG of linear pushbroom imagery more effectively. However, the PP model does not completely comply with the imaging geometry of linear pushbroom imagery, and only rigorous models describe the scene formation as it actually happens. Therefore, these models are the most accurate ones [8,24] and have been adopted in a variety of applications [25][26][27][28]. Consequently, it is theoretically and practically significant to develop an approach for epipolar resampling of linear pushbroom imagery based on rigorous sensor models.
In this paper, a novel epipolar resampling approach of cross-track linear pushbroom satellite imagery is proposed based on the MPC model. In the proposed method, the vertical parallaxes of CIPs are eliminated by modifying the instantaneous attitude of the stereo imagery. In this way, the required new attitude parameters are computed based on the direction of the IB of each pair of CIPs. Because the instantaneous position of sensors remains unchanged, there is no need for DEM of the ground coverage in the resampling process.

Linear Pushbroom Imaging
In the case of a frame camera, the imaging process is performed by a 2D array of detectors (CCD or CMOS) arranged within its focal plane. In contrast, linear pushbroom scanners have only a long row of detectors within their focal plane; consequently, they capture a 1D (linear) image in a given exposure [29]. A sequence of linear images in the flying direction can be captured by moving the sensor which forms the 2D image frame [29,30]. Therefore, every 1-D image is associated with the specific positional and orientation parameters of the sensor during the time of exposure. Then, each 1-D image will have a distinct set of EOP. This different imaging geometry leads to different mathematical modeling of linear pushbroom imagery. Toward this end, a wide variety of mathematical models have been developed which generally fall into two main rigorous and empirical models [31,32]. Reconstruction of the geometry of an image at the time of imaging is the backbone of rigorous models. Two well-known approaches in this area encompass Orbital Parameters [27,33,34] and Multiple Projection Centers models [6][7][8][9]35]. In the first case, the trajectory of the sensor is modeled based on Keplerian orbital parameters and time-dependent (temporal) polynomials are used to estimate attitude parameters. In the second type of rigorous models, both the sensor trajectory and attitude parameters are estimated using temporal polynomials. In contrast to these models, no physical or geometrical interpretation stands behind empirical models and they attempt to estimate unknown parameters of the model by fitting it to the rigorous sensor model or a set of GCPs. Numerous research studies have been carried out on investigation of the empirical models in the last decade [36][37][38][39].

MPC Model
In the MPC model, each scanline of linear pushbroom imagery is supposed as a central perspective image, which has a set of unique EOP. Therefore, each scanline has a distinct Image Reference Frame (IRF), which x and y axes are directed toward the satellite's instantaneous velocity vector and the array of sensor's detectors, respectively. The origin of IRF is the principal point on the image space, which is usually considered at the middle of scanline (Figure 1a). In this model, any arbitrary Ground Reference Frame (GRF) can be used to explain the object space's coordinates. For convenience in later computations, it is usually adopted according to a map coordinate system (e.g., UTM).
In order to establish the relation between IRF and GRF, two mediator coordinate systems are defined: Sensor Reference Frame (SRF) and Platform Reference Frame (PRF). SRF is a pseudo-3D coordinate system with its origin on the exposure station. Its x and y axes coincide with the corresponding axes of IRF, and its z axis passes through the instantaneous exposure station (Figure 1b). The relationship between IRF and SRF is established based on the IOP of the sensor. PRF is co-origin with the SRF, and its axes are locked to the body of the platform (Figure 1c).

Linear Pushbroom Imaging
In the case of a frame camera, the imaging process is performed by a 2D array of detectors (CCD or CMOS) arranged within its focal plane. In contrast, linear pushbroom scanners have only a long row of detectors within their focal plane; consequently, they capture a 1D (linear) image in a given exposure [29]. A sequence of linear images in the flying direction can be captured by moving the sensor which forms the 2D image frame [29,30]. Therefore, every 1-D image is associated with the specific positional and orientation parameters of the sensor during the time of exposure. Then, each 1-D image will have a distinct set of EOP. This different imaging geometry leads to different mathematical modeling of linear pushbroom imagery. Toward this end, a wide variety of mathematical models have been developed which generally fall into two main rigorous and empirical models [31,32]. Reconstruction of the geometry of an image at the time of imaging is the backbone of rigorous models. Two well-known approaches in this area encompass Orbital Parameters [27,33,34] and Multiple Projection Centers models [6][7][8][9]35]. In the first case, the trajectory of the sensor is modeled based on Keplerian orbital parameters and time-dependent (temporal) polynomials are used to estimate attitude parameters. In the second type of rigorous models, both the sensor trajectory and attitude parameters are estimated using temporal polynomials. In contrast to these models, no physical or geometrical interpretation stands behind empirical models and they attempt to estimate unknown parameters of the model by fitting it to the rigorous sensor model or a set of GCPs. Numerous research studies have been carried out on investigation of the empirical models in the last decade [36][37][38][39].

MPC Model
In the MPC model, each scanline of linear pushbroom imagery is supposed as a central perspective image, which has a set of unique EOP. Therefore, each scanline has a distinct Image Reference Frame (IRF), which x and y axes are directed toward the satellite's instantaneous velocity vector and the array of sensor's detectors, respectively. The origin of IRF is the principal point on the image space, which is usually considered at the middle of scanline (Figure 1a). In this model, any arbitrary Ground Reference Frame (GRF) can be used to explain the object space's coordinates. For convenience in later computations, it is usually adopted according to a map coordinate system (e.g., UTM).
In order to establish the relation between IRF and GRF, two mediator coordinate systems are defined: Sensor Reference Frame (SRF) and Platform Reference Frame (PRF). SRF is a pseudo-3D coordinate system with its origin on the exposure station. Its x and y axes coincide with the corresponding axes of IRF, and its z axis passes through the instantaneous exposure station ( Figure  1b). The relationship between IRF and SRF is established based on the IOP of the sensor. PRF is coorigin with the SRF, and its axes are locked to the body of the platform (Figure 1c).  Like in aerial photogrammetry, the sensor is supposed to be locked to the body of the platform in the MPC. Therefore, the axes of SRF and PRF coincide during the acquisition time. Thus, the 3D rotation matrix required to transform from SRF to PRF will be an identity matrix, given by Equation (1): where [R Attitude ] is a 3D rotation matrix from the SRF to the PRF. Transformation from the PRF to the GRF is performed using a 3D rotation matrix through modeling the orientation parameters of the sensor (Equation (2)): where ω t , φ t , and κ t are the orientation angles of the sensor around the x, y and z axes of the GRF, respectively. R r (for r = 1, 2, 3) is a 3D rotation matrix around the r-th axis of the GRF. Finally, [R Orientation ] is a 3D rotation matrix from the PRF to the GRF. In practice, the orientation angles ω t , φ t , and κ t are modeled as time-dependent polynomials, given by Equation (3): where t is the time parameter which can be considered equivalent to the rows-number of the image points. Obviously, the type and number of essential terms included in the polynomials depend on the perturbations that occurred during image formation [29]. Given the known relation between the SRF, the PRF and the GRF, the overall structure of the co-linearity equation can be expressed using Equation (4): where (x = 0, y) is the coordinate of points in the SRF, c is the principal distance of the sensor, λ is the scale factor, (X, Y, Z) is the object space coordinate in the GRF, and finally, (X S , Y S , Z S ) is the instantaneous exposure station in the GRF which is modeled by the time-dependent polynomials [29] given by Equation (5): where X r , Y r and Z r (for r = 0, 1, 2) are the unknown coefficients of the polynomials, and t is the time parameter which can be considered equivalent to the rows-number of the image points.

The EG of Linear Pushbroom Images
The main difference between frame-type cameras and linear pushbroom sensors is that each scanline of linear imagery has a distinct perspective center. As a result of this characteristic, there will be several epipolar planes on the second scene for a given point p in the first scene [1], and hyperbola-shape epipolar curves (instead of epipolar lines) of linear imagery is due to the non-coplanarity of these planes [1,8] (Figure 2a). For cross-track stereo coverage, the necessary and sufficient condition for coplanarity of all epipolar planes and formation of straight epipolar lines is that the scanlines containing CIPs to be parallel with their IB in a common plane (Figure 2b). With a little attention to Figure 2b it can be easily understood that for a given pair of CIPs, the aforementioned condition can be reconstructed by modifying the attitude parameters of scanlines containing the CIPs. In this way, the epipolar curves of these CIPs will convert to straight Conjugate Epipolar Lines (CEL); and the plane containing all possible epipolar planes of these points will intersect any horizontal plane of the object space in a straight line parallel with their IB (Figure 2b). Hereafter, this direction is called the CEL's direction. Due to the dynamism of two sensors during the capture of stereo scenes, by changing the rows-number of image points, the direction of IB of CIPs will change, as well as their CEL's direction ( Figure 3a). Moreover, given two distinct points such as p1 and p2 in the i-th scanline of the first scene, their conjugate points will not necessarily be in a unique scanline of the second scene ( Figure 3b). Thus, by changing the columns-number of image points, the direction of IB of CIPs and their CEL's direction will change too. From now on, these two changes in the CEL's direction will be called the temporal and spatial variations (due to the change of rows-and columns-number of CIPs, respectively) of the IB of CIPs, respectively.   With a little attention to Figure 2b it can be easily understood that for a given pair of CIPs, the aforementioned condition can be reconstructed by modifying the attitude parameters of scanlines containing the CIPs. In this way, the epipolar curves of these CIPs will convert to straight Conjugate Epipolar Lines (CEL); and the plane containing all possible epipolar planes of these points will intersect any horizontal plane of the object space in a straight line parallel with their IB (Figure 2b). Hereafter, this direction is called the CEL's direction. Due to the dynamism of two sensors during the capture of stereo scenes, by changing the rows-number of image points, the direction of IB of CIPs will change, as well as their CEL's direction ( Figure 3a). Moreover, given two distinct points such as p 1 and p 2 in the i-th scanline of the first scene, their conjugate points will not necessarily be in a unique scanline of the second scene ( Figure 3b). Thus, by changing the columns-number of image points, the direction of IB of CIPs and their CEL's direction will change too. From now on, these two changes in the CEL's direction will be called the temporal and spatial variations (due to the change of rows-and columns-number of CIPs, respectively) of the IB of CIPs, respectively. With a little attention to Figure 2b it can be easily understood that for a given pair of CIPs, the aforementioned condition can be reconstructed by modifying the attitude parameters of scanlines containing the CIPs. In this way, the epipolar curves of these CIPs will convert to straight Conjugate Epipolar Lines (CEL); and the plane containing all possible epipolar planes of these points will intersect any horizontal plane of the object space in a straight line parallel with their IB (Figure 2b). Hereafter, this direction is called the CEL's direction. Due to the dynamism of two sensors during the capture of stereo scenes, by changing the rows-number of image points, the direction of IB of CIPs will change, as well as their CEL's direction ( Figure 3a). Moreover, given two distinct points such as p1 and p2 in the i-th scanline of the first scene, their conjugate points will not necessarily be in a unique scanline of the second scene ( Figure 3b). Thus, by changing the columns-number of image points, the direction of IB of CIPs and their CEL's direction will change too. From now on, these two changes in the CEL's direction will be called the temporal and spatial variations (due to the change of rows-and columns-number of CIPs, respectively) of the IB of CIPs, respectively.  The necessary and sufficient condition for eliminating the vertical parallaxes of all CIPs throughout the linear cross-track stereo imagery is that the epipolar lines of all image points be parallel with each other, and aligned along the scene rows. Therefore, when the attitude parameters of stereo scenes are modified, the vertical parallax of all CIPs can be eliminated indiscriminately by applying a complementary rotation. Due to the temporal and spatial variations of the direction of IB of the CIPs, the required rotation angle will obviously be a function of both the rows-and the columns-number of CIPs. The necessary and sufficient condition for eliminating the vertical parallaxes of all CIPs throughout the linear cross-track stereo imagery is that the epipolar lines of all image points be parallel with each other, and aligned along the scene rows. Therefore, when the attitude parameters of stereo scenes are modified, the vertical parallax of all CIPs can be eliminated indiscriminately by applying a complementary rotation. Due to the temporal and spatial variations of the direction of IB of the CIPs, the required rotation angle will obviously be a function of both the rows-and the columns-number of CIPs.  As satellite scenes are usually acquired in a very short time and they benefit from a fairly stable attitude [8,29], the temporal and spatial variations of the IB's direction of CIPs can be modeled using some polynomials in the both rows-and columns-number of CIPs. To validate this theoretical consideration, the IB's direction of some pairs of CIPs from a real dataset (the left scene of the Isfahan dataset, which is introduced in Section 4) were calculated in terms of three orientation parameters roll (ω), pitch (φ), and yaw (κ) of the sensor. The scatter plots of computed parameters are illustrated in Figure 4, with respect to the rows-and the columns-number of CIPs.
According to Figure 4, no one of the attitude parameters ω, φ, and κ can exactly be modeled in terms of only one of the rows-or the columns-number of CIPs. Thereafter, some polynomials in the rows-and the columns-number of CIPs were used to estimate the computed attitude parameters. The scatter plot of computed and estimated attitude parameters of the IB of CIPs is shown in Figure 5, which illustrates the quality of estimation. As satellite scenes are usually acquired in a very short time and they benefit from a fairly stable attitude [8,29], the temporal and spatial variations of the IB's direction of CIPs can be modeled using some polynomials in the both rows-and columns-number of CIPs. To validate this theoretical consideration, the IB's direction of some pairs of CIPs from a real dataset (the left scene of the Isfahan dataset, which is introduced in Section 4) were calculated in terms of three orientation parameters roll (ω), pitch (φ), and yaw (κ) of the sensor. The scatter plots of computed parameters are illustrated in Figure 4, with respect to the rows-and the columns-number of CIPs.
According to Figure 4, no one of the attitude parameters ω, φ, and κ can exactly be modeled in terms of only one of the rows-or the columns-number of CIPs. Thereafter, some polynomials in the rows-and the columns-number of CIPs were used to estimate the computed attitude parameters. The scatter plot of computed and estimated attitude parameters of the IB of CIPs is shown in Figure 5, which illustrates the quality of estimation. In order to a quantitative assessment, the obtained coefficients of the polynomials fitted to each of the attitude parameters are provided in Table 1, along with their significance values (p-value), standard errors (std-err), coefficients of determination (Pseudo-R 2 ), and the root mean squares of residuals from the fitted straight lines ( 0  ); which confirm the high accuracy of the modeling.
According to Table 1, since the IB's direction of CIPs has a systematic behavior, the complementary rotation angle required for parallelizing the epipolar lines of CIPs can be modeled using some polynomials in the both rows-and columns-number of CIPs, as well as the new attitude parameters required to parallelizing the scanlines containing the CIPs with their IBs. In order to a quantitative assessment, the obtained coefficients of the polynomials fitted to each of the attitude parameters are provided in Table 1, along with their significance values (p-value), standard errors (std-err), coefficients of determination (Pseudo-R 2 ), and the root mean squares of residuals from the fitted straight lines (σ 0 ); which confirm the high accuracy of the modeling.
According to Table 1, since the IB's direction of CIPs has a systematic behavior, the complementary rotation angle required for parallelizing the epipolar lines of CIPs can be modeled using some polynomials in the both rows-and columns-number of CIPs, as well as the new attitude parameters required to parallelizing the scanlines containing the CIPs with their IBs.  Functional form:

Proposed Method
Suppose two conjugate image points p and p , which are respectively in the i-th and the j-th scanlines of the first and the second imagery. The proposed epipolar resampling method for cross-track linear imagery is constructed in three steps: • parallelizing scanlines containing the CIPs with their IBs to produce straight epipolar lines, • parallelizing the epipolar lines of all CIPs with each other to eliminate the vertical parallax of CIPs, and, • correcting the scale of the normalized scenes.
These steps are explained in more details in the following subsections.

Producing Straight Epipolar Lines
In the first step, the y axis of scanlines containing the CIPs should be parallelized with their IB in a common plane, by modifying the attitude parameters of the sensor (Figure 2b). In this way, the direction of IB of CIPs should be computed firstly. Therefore, given the EOP of the stereo imagery, the exposure stations of the i-th scan line of the first scene and the j-th scan line of the second scene, are calculated using Equation (5). Then, the IB of points p and p is computed using Equation (6): where S i and S j are the exposure stations of the i-th scan line of the first scene and the j-th scan line of the second scene, respectively; and B GRF is the IB of conjugate points p and p in the GRF.
Since the calculated IB is in the GRF and the new attitude parameters for parallelizing the SRF's y axis of scanlines are required in the PRF, the vector B GRF should firstly be transformed to the PRF of each scanline, see Equation (7): where  (8) and (9): where ω i n and κ i n are the new attitude parameters around the PRF's x and z axes of the i-th scan line of the first image, respectively.
Moreover, for coplanarity of the i-th scanline of the first image and the j-th scanline of the second image, it is sufficient that the new attitude parameters around their PRF's y axis be equal. In order to simultaneously account for the off-nadir viewing effect of the sensor, it is recommended that this parameter be considered equal to zero, as indicated by Equation (10): where φ i n is the new attitude parameter around the PRF's y axis of the i-th scanline of the first image. Given the new attitude parameters, the required rotation matrix for parallelizing the SRF's y axis of the i-th scanline of the first image with the IB of the conjugate points p and p (i.e., [R n Attitud ] i ) will be as Equation (11): In practice, two parameters ω i n and κ i n should firstly be calculated for a set of CIPs, and by fitting a polynomial with a proper order in both the rows-and the columns-number of CIPs then be used in the structure of rotation matrix [R n Attitud ] i , Equation (12): where i and l are the rows-and the columns-number of CIPs, respectively. ω n r and κ n r (for r = 0, 1, 2) are the unknown coefficients of polynomials.
Given the second image's EOP, the required rotation matrix for parallelizing the SRF's y axis of the j-th scanline of the second image with the IB of the conjugate points p and p (i.e., [R n Attitud ] j ) can similarly be calculated using Equations (8)- (12).
By applying the rotation matrix [R n Attitud ] i , the points' coordinate in the aligned-SRF (i.e., the SRF after than its y axis was parallelized with the IB of CIPs) will be as indicated by Equation (13): where (x a , y a , −c a ) is the points' coordinate in the aligned-SRF (SRF a ), and λ a is the scale factor for transferring from GRF to SRF a . Since [R Attitude ] = I 3 and, [R Orientation ] T × [R Orientation ] = I 3 , Equation (13) can be simplified as Equation (14): where k = λ λ a is the scale factor for transferring from SRF to SRF a . As previously mentioned, in the SRF a the epipolar curves will become straight epipolar lines which are aligned along the IB of CIPs.

Eliminating the Vertical Parallax of CIPs
Due to the temporal and spatial variations of the IB's direction of CIPs, the CEL's direction throughout the linear stereo imagery will be variable as well. Therefore, in order to the vertical parallaxes of all CIPs be eliminated at ones, it is necessary that all the CELs be parallelized by applying a complementary rotation around the SRF a 's z axis (Equation (15)): where (x p , y p , −c p ) is the points' coordinate in the parallelized-SRF (i.e., the SRF a after that the complementary rotation was applied), θ n is the complementary rotation angle for parallelizing the CEL throughout the stereo imagery; which in order to account for the temporal and spatial variations of the IB's direction of CIPs, it is modeled by a polynomial in both the rows-and the columns-number of CIPs (Equation (16)): where θ n r (for r = 0, 1, 2, . . . ) are the unknown coefficients of the polynomial. Given the new attitude parameter of the second imagery, the coordinate of any image point in its parallelized-SRF (SRF p ) can similarly be calculated using Equations (13)- (16).
The fundamental condition for parallelism of all epipolar lines throughout the stereo imagery is that the vertical parallaxes of all CIPs be equal to a constant value. In order to computing the vertical parallax of CIPs, their x-ordinate in the SRF p should firstly be calculated. In this way, by substituting Equation (14) into Equation (15) and dividing the first row of Equation (15) by its third row, the x-ordinate of any image point in its SRF p can be computed as Equation (17): where x i p is the point's x-ordinate in the SRF p of the i-th scanline of the first image. Since each scanline of linear imagery has a unique SRF, Equation (17) can only be used for relating the SRF of each image point with its SRF p . In order to establish the relationship of SRF p of consecutive scanlines, a psudo-2D coordinate system was constructed by arranging the scanlines in a sequential manner (Equation (18)): where i and j are the rows-number of CIPs in the first and second imagery, respectively; D CCD is the dimension of sensor's CCD, and finally, x p psud and x p psud are the x-ordinate of CIPs in the psudo-2D coordinate system of the first and the second imagery, respectively. By the way, the condition equation of obtaining the parallelized epipolar lines can be written as Equation (19): where d is the vertical parallax of CIPs; which should be a constant value. By substituting Equations (14)- (18) into Equation (19), the final formula to compute the rotation angle θ n can be written as Equation (20): The above equation is non-linear in terms of the unknowns (i.e., θ n , θ n and d), which after differentiation and linearization, can be solved using a set of CIPs. By applying the complementary rotation, CIPs will have no vertical parallax in their SRF p . But due to the variations of the IB's length of CIPs, flying height, and the sensor's roll angle throughout the stereo imagery, the scale of generated model will be variable.

Correcting the Scale of the Normalized Scenes
In the third step the model's scale should be corrected by shifting the CIPs along the y axis of their SRF p (Equation (21)): where (x n , y n , −c n ) the point's coordinate in the normalized-SRF (SRF n ), and ∆y is the required shift to correct the scale of the point p in the first image. Due to the temporal and spatial variations of the IB of CIPs, the sensor's roll angle, and flying height, the required shift will obviously be a function of the rows-and the columns-number of CIPs, Equation (22): where ∆y r (for r = 0, 1, 2) are the unknown coefficients of the polynomial. In order to computing this parameter, the points' y-ordinate should firstly be calculated in their SRF p (Equation (23)): where y i p is the points' y-ordinate in the SRF p of the i-th scan line of the first image. The necessary condition to correct the model's scale is that after than the CIPs is shifted along the y axis of their SRF p , horizontal parallaxes of CIPs and their heights in the object space to be linearly proportional, as indicated by Equation (24): where ∆y and ∆y are the required shifts to correct the scales of the first and the second images, respectively, Z is the height of CIPs in the GRF, a and b are the coefficients of the linear relation of horizontal parallaxes and heights. If b = 0, and a is considered to be equal to the average imaging scale, Equation (24) can be rewritten as: where c is the sensor's principal distance and H m is the average flying height. Finally, by substituting Equations (21)- (24) into Equation (25) and forming a system of equations for a set of CIPs, the required shifts are achieved as a function of both the rows-and the columns-number of CIPs. It should be noted that given the EOP of stereo imagery and the image coordinate of CIPs, the points' heights in the object space can be calculated via the space intersection of two images. Therefore, no additional GCPs are needed for Equation (25) to be to solved.
Given the new attitude parameters, the complementary rotation angle and the required shift to correct the model's scale, the final formula to transfer from SRF to SRF n will be as Equation (26): where (x n , y n , −c n ) is the points' coordinate in the SRF n . In practice, the invers form of Equation (26) will be used in the epipolar resampling process to transfer image points from normal image space to the original image space and indirect resampling of the linear stereo scenes. In this way, the epipolar resampling procedure can be conducted in a point wise manner as well as the digital rectification process.

Study Area and Data Used
Three different datasets have been used in this article, including two SPOT-1A stereo imagery sets covering Isfahan and Zanjan provinces, and one RapidEye stereo imagery set covering Fars province, all in Iran. The first and the second datasets cover a foothills rural area with rare urbanized regions. There is a lake at the bottom-right corner of the second dataset, and a little cloud coverage can be seen in this area. The third dataset covers a semi-mountainous foothill area, including two quite big cities almost in the center of the image. A little cloud coverage at the middle-left area of stereo image overlaps can be observed. Cloud percentages, geometric information and the number of available GCPs for each dataset are shown in Table 2, along with their elevation relief. GCPs of Isfahan and Zanjan datasets were measured using a double-frequency GPS with a sub-meter accuracy, and their corresponding image coordinates were measured with an approximate precision of 0.5 pixel. For the Fars dataset, GCPs were extracted from 1:2000 scale digital topographic maps produced by the National Cartographic Center of Iran with 0.6 m planimetric and 0.5 m altimetric accuracies, respectively. The points are distinct features such as buildings, pool corners, walls, and road junctions. Their corresponding image coordinates were measured with an approximate precision of one pixel. Distribution of GCPs is illustrated in Figure 6. The lack of GCPs in Zanjan and Fars datasets are due to their land coverage.

Study Area and Data Used
Three different datasets have been used in this article, including two SPOT-1A stereo imagery sets covering Isfahan and Zanjan provinces, and one RapidEye stereo imagery set covering Fars province, all in Iran. The first and the second datasets cover a foothills rural area with rare urbanized regions. There is a lake at the bottom-right corner of the second dataset, and a little cloud coverage can be seen in this area. The third dataset covers a semi-mountainous foothill area, including two quite big cities almost in the center of the image. A little cloud coverage at the middle-left area of stereo image overlaps can be observed. Cloud percentages, geometric information and the number of available GCPs for each dataset are shown in Table 2, along with their elevation relief. GCPs of Isfahan and Zanjan datasets were measured using a double-frequency GPS with a submeter accuracy, and their corresponding image coordinates were measured with an approximate precision of 0.5 pixel. For the Fars dataset, GCPs were extracted from 1:2000 scale digital topographic maps produced by the National Cartographic Center of Iran with 0.6 m planimetric and 0.5 m altimetric accuracies, respectively. The points are distinct features such as buildings, pool corners, walls, and road junctions. Their corresponding image coordinates were measured with an approximate precision of one pixel. Distribution of GCPs is illustrated in Figure 6. The lack of GCPs in Zanjan and Fars datasets are due to their land coverage.

Results
According to Section 3, there is no need to any GCPs if the EOP of stereo imagery are available; and the proposed method can be solved using only a few CIPs. In this study, the EOP of imagery were computed using the GCPs. At the space resection stage, in order to optimally structure the rotation matrix [ROrientation], different structures were investigated in a try and error manner. From the available GCPs for each dataset, 15 points were used at the space resection stage while the remaining points were used as independent check points in the accuracy assessment process. Accuracies obtained from space resection of each image are illustrated in Table 3 in the image space, along with the accuracies obtained from space intersection of the stereo imagery in the object space.

Results
According to Section 3, there is no need to any GCPs if the EOP of stereo imagery are available; and the proposed method can be solved using only a few CIPs. In this study, the EOP of imagery were computed using the GCPs. At the space resection stage, in order to optimally structure the rotation matrix [R Orientation ], different structures were investigated in a try and error manner. From the available GCPs for each dataset, 15 points were used at the space resection stage while the remaining points were used as independent check points in the accuracy assessment process. Accuracies obtained from space resection of each image are illustrated in Table 3 in the image space, along with the accuracies obtained from space intersection of the stereo imagery in the object space. Then, 40 pairs of well distributed CIPs were manually extracted from each stereo imagery set. These points were employed in two different roles. One part of CIPs was used to compute the new attitude parameters, complementary rotation angles, and the required shifts; which is called Control Conjugate Points (CoCP) in this paper. The other part of CIPs, which is called Check Conjugate Points (ChCP), was used to assess the model's accuracy. In order to evaluate the accuracy of the proposed method as it is impacted by the number of utilized CoCP, several experiments were performed on each dataset. In each experiment, using the parameters which are computed on CoCP, the ChCP's coordinate were transformed to the SRF n . The mean absolute values and the maximum absolute values of the residual vertical parallax (P V ) of ChCP are provide in Table 4, along with the square root of the estimated variance component from straight-line fitting to horizontal parallaxes (P H ) of ChCP and their heights in the object space. In this way, the elevation of ChCP were calculated using the space intersection of their image coordinate in the original stereo imagery. Table 4. Experimental results obtained from the proposed epipolar resampling method. 1  2  3  4  5  6  7  8  9  10  11  12  Number of CoCP  8  10  12  16  8  10  12  16  8  10  12  16  Number of ChCP  32  30  28  24  32  30  28  24  32  30  28  According to Table 4, the proposed method has a high ability to eliminate the vertical parallaxes of CIPs even if a few CoCP are used, so by increasing the number of CoCP from 10 to 16 (i.e., Experiments 3,4,7,8,11 and 12), no significant improvement is seen in the results. Finally, the produced normalized scenes through the Experiments 2, 6 and 10 were overlayed to generate the corresponding stereo anaglyphs (Figure 7).  Then, 40 pairs of well distributed CIPs were manually extracted from each stereo imagery set. These points were employed in two different roles. One part of CIPs was used to compute the new attitude parameters, complementary rotation angles, and the required shifts; which is called Control Conjugate Points (CoCP) in this paper. The other part of CIPs, which is called Check Conjugate Points (ChCP), was used to assess the model's accuracy. In order to evaluate the accuracy of the proposed method as it is impacted by the number of utilized CoCP, several experiments were performed on each dataset. In each experiment, using the parameters which are computed on CoCP, the ChCP's coordinate were transformed to the SRFn. The mean absolute values and the maximum absolute values of the residual vertical parallax (PV) of ChCP are provide in Table 4, along with the square root of the estimated variance component from straight-line fitting to horizontal parallaxes (PH) of ChCP and their heights in the object space. In this way, the elevation of ChCP were calculated using the space intersection of their image coordinate in the original stereo imagery. According to Table 4, the proposed method has a high ability to eliminate the vertical parallaxes of CIPs even if a few CoCP are used, so by increasing the number of CoCP from 10 to 16 (i.e., Experiments 3,4,7,8,11 and 12), no significant improvement is seen in the results. Finally, the produced normalized scenes through the Experiments 2, 6 and 10 were overlayed to generate the corresponding stereo anaglyphs (Figure 7).

Discussion
In order to assess the accuracy of the proposed method two factors were examined: (1) the ability to eliminate the vertical parallaxes of CIPs; and (2) to provide a linear relationship between the

Discussion
In order to assess the accuracy of the proposed method two factors were examined: (1) the ability to eliminate the vertical parallaxes of CIPs; and (2) to provide a linear relationship between the horizontal parallaxes of CIPs and their heights in the object space. Before applying the proposed epipolar resampling method, the mean values of vertical parallax of CIPs for Isfahan, Zanjan, and Fars datasets were 174, 184.4, and 106.7 pixels, respectively. According to Table 4, in all experiments the maximum absolute values of the residual vertical parallax of normalized scenes are close to zero. Therefore, the proposed method has a high ability to eliminate the vertical parallaxes of CIPs. The residual vertical parallax of CIPs in Experiments 2, 6, and 10 are illustrated in Figure 8, where the horizontal axis shows the points' number ordered by increasing the residual parallax, while the vertical axis shows the value of vertical parallax. horizontal parallaxes of CIPs and their heights in the object space. Before applying the proposed epipolar resampling method, the mean values of vertical parallax of CIPs for Isfahan, Zanjan, and Fars datasets were 174, 184.4, and 106.7 pixels, respectively. According to Table 4, in all experiments the maximum absolute values of the residual vertical parallax of normalized scenes are close to zero. Therefore, the proposed method has a high ability to eliminate the vertical parallaxes of CIPs. The residual vertical parallax of CIPs in Experiments 2, 6, and 10 are illustrated in Figure 8, where the horizontal axis shows the points' number ordered by increasing the residual parallax, while the vertical axis shows the value of vertical parallax. Moreover, in order to examine the pattern of the vertical parallaxes on the normalized scenes, the scatter plot of the vertical parallaxes of CIPs in Experiments 2, 6, and 10 are shown in Figure 9. In this figure, the CoCP and ChCP are illustrated with the triangle and circle markers, respectively. The vertical parallaxes of CIPs are drawn with a magnification factor of 10,000.
According to Figure 9, a scatter plot of the residual vertical parallax of ChCP shows no systematic behavior. The proposed method has a low sensitivity to the number of CoCP used, as by increasing this number from 8 to 16, no appreciable variation in the results is observed. Generally, the minimum number of required CoCP can be variable regarding to the number of parameters used in the polynomials of Equations 12, 16, and 22. However, the distribution of these points on the stereo imagery is of great importance. Repetitive experiments show that selecting the CoCP around the outer boundary of the stereo images' common coverage plays a key role on the accuracy of the model's parameters (Figure 9). Moreover, in order to examine the pattern of the vertical parallaxes on the normalized scenes, the scatter plot of the vertical parallaxes of CIPs in Experiments 2, 6, and 10 are shown in Figure 9. In this figure, the CoCP and ChCP are illustrated with the triangle and circle markers, respectively. The vertical parallaxes of CIPs are drawn with a magnification factor of 10,000. According to Figure 9, a scatter plot of the residual vertical parallax of ChCP shows no systematic behavior. The proposed method has a low sensitivity to the number of CoCP used, as by increasing this number from 8 to 16, no appreciable variation in the results is observed. Generally, the minimum number of required CoCP can be variable regarding to the number of parameters used in the polynomials of Equations (12), (16) and (22). However, the distribution of these points on the stereo imagery is of great importance. Repetitive experiments show that selecting the CoCP around the outer boundary of the stereo images' common coverage plays a key role on the accuracy of the model's parameters (Figure 9).  Figure 8, where the horizontal axis shows the points' number ordered by increasing the residual parallax, while the vertical axis shows the value of vertical parallax. Moreover, in order to examine the pattern of the vertical parallaxes on the normalized scenes, the scatter plot of the vertical parallaxes of CIPs in Experiments 2, 6, and 10 are shown in Figure 9. In this figure, the CoCP and ChCP are illustrated with the triangle and circle markers, respectively. The vertical parallaxes of CIPs are drawn with a magnification factor of 10,000. According to Figure 9, a scatter plot of the residual vertical parallax of ChCP shows no systematic behavior. The proposed method has a low sensitivity to the number of CoCP used, as by increasing this number from 8 to 16, no appreciable variation in the results is observed. Generally, the minimum number of required CoCP can be variable regarding to the number of parameters used in the polynomials of Equations 12, 16, and 22. However, the distribution of these points on the stereo imagery is of great importance. Repetitive experiments show that selecting the CoCP around the outer boundary of the stereo images' common coverage plays a key role on the accuracy of the model's parameters (Figure 9). As the final note on Table 4, the proposed model has a high ability to provide a linear relationship between the vertical parallax of ChCP and their heights in the object space. Since the height of ChCP have been calculated using the space intersection process and the EOP of original stereo imagery have been directly estimated using GCPs, the computed heights' accuracy of ChCP depends on the accuracy and distribution of GCPs, along with the image coordinates' accuracy of both the CoCP and ChCP. The scatter plots of the horizontal parallax of ChCP and their heights in Experiments 2, 6, and 10 are illustrated in Figure 10. According to Figure 10, the proposed method benefits from an admissible ability to provide a linear relationship between the vertical parallax of ChCP and their heights in the object space.

Conclusions
This paper presents a novel epipolar resampling approach for cross-track linear pushbroom images. The resampling process is based on the rigorous sensor model, with no approximations or assumptions. In this way, the epipolar curves of linear imagery are generally hyperbola-shaped, but if the scanlines containing the CIPs are parallelized with their IB, these curves are converted into straight lines. However, due to the temporal and spatial variation of the IB's direction of CIPs, these epipolar lines will not be parallel to each other. Thus, in order to perform epipolar resampling of linear imagery by means of the rigorous sensor model, the attitude parameters of the stereo images should be modified and a complementary rotation should then be applied. Since the new attitude parameters and the complementary rotation angle are dependent to the IB's direction of CIPs, using As the final note on Table 4, the proposed model has a high ability to provide a linear relationship between the vertical parallax of ChCP and their heights in the object space. Since the height of ChCP have been calculated using the space intersection process and the EOP of original stereo imagery have been directly estimated using GCPs, the computed heights' accuracy of ChCP depends on the accuracy and distribution of GCPs, along with the image coordinates' accuracy of both the CoCP and ChCP. As the final note on Table 4, the proposed model has a high ability to provide a linear relationship between the vertical parallax of ChCP and their heights in the object space. Since the height of ChCP have been calculated using the space intersection process and the EOP of original stereo imagery have been directly estimated using GCPs, the computed heights' accuracy of ChCP depends on the accuracy and distribution of GCPs, along with the image coordinates' accuracy of both the CoCP and ChCP. The scatter plots of the horizontal parallax of ChCP and their heights in Experiments 2, 6, and 10 are illustrated in Figure 10. According to Figure 10, the proposed method benefits from an admissible ability to provide a linear relationship between the vertical parallax of ChCP and their heights in the object space.

Conclusions
This paper presents a novel epipolar resampling approach for cross-track linear pushbroom images. The resampling process is based on the rigorous sensor model, with no approximations or assumptions. In this way, the epipolar curves of linear imagery are generally hyperbola-shaped, but if the scanlines containing the CIPs are parallelized with their IB, these curves are converted into straight lines. However, due to the temporal and spatial variation of the IB's direction of CIPs, these epipolar lines will not be parallel to each other. Thus, in order to perform epipolar resampling of linear imagery by means of the rigorous sensor model, the attitude parameters of the stereo images should be modified and a complementary rotation should then be applied. Since the new attitude parameters and the complementary rotation angle are dependent to the IB's direction of CIPs, using The scatter plots of the horizontal parallax of ChCP and their heights in Experiments 2, 6, and 10 are illustrated in Figure 10. According to Figure 10, the proposed method benefits from an admissible ability to provide a linear relationship between the vertical parallax of ChCP and their heights in the object space.

Conclusions
This paper presents a novel epipolar resampling approach for cross-track linear pushbroom images. The resampling process is based on the rigorous sensor model, with no approximations or assumptions. In this way, the epipolar curves of linear imagery are generally hyperbola-shaped, but if the scanlines containing the CIPs are parallelized with their IB, these curves are converted into straight lines. However, due to the temporal and spatial variation of the IB's direction of CIPs, these epipolar lines will not be parallel to each other. Thus, in order to perform epipolar resampling of linear imagery by means of the rigorous sensor model, the attitude parameters of the stereo images should be modified and a complementary rotation should then be applied. Since the new attitude parameters and the complementary rotation angle are dependent to the IB's direction of CIPs, using some polynomials in both the rows-and the columns-number of CIPs has been proposed for modeling these parameters. As the instantaneous position of sensors remains fixed, the DEM of the area of interest is not required in the resampling process. Given the EOP of stereo imagery, there is no need to any GCP, and the proposed method can be performed using only a few CIPs. Experimental results obtained from two pairs of SPOT-1A and one pair of RapidEye stereo imagery with different terrain conditions proved the feasibility and the success of the proposed method. One of disadvantages of the model is the need for geometrically raw stereo images, which isn't offered by some providers. However, in the case of archive imagery and those which aren't supplied with reliable RPCs, the proposed method can be very effective and almost the only option currently available.
Future studies will focus on a study of the feasibility of eliminating the need for CIPs to carry out the normalization procedure. Moreover, the proposed method will be developed for along-track stereo imaging systems, as well as those images captured by the manned and unmanned aerial platforms. Finally, DEM and orthophoto generation based on the normalized scenes will be considered.