Nomenclature

X

Scalar position

V

Velocity

A

Acceleration

R

Position vector

W

Rotational speed vector

AC

Aircraft

NED

North (x), East (y), Down (z) coordinate frame, with the centre on aircraft's centre of mass

Inertia

Inertial coordinate frame, here assumed to be the same as NED

R_a/b

R_a – R_b

$C_I^B$

Rotation matrix, converting a vector in inertial coordinates to body coordinates

V|_cam

Velocity vector in camera coordinates

h

Pixel position

x, y, z

Position of point in camera coordinates

u, v

Normalised pixel positions u = x/z, v = y/z

N

Number of image rows

o_x, o_y

Pinhole model offsets (in pixels)

f_x, f_y

Pinhole model focal lengths (in pixels)

k ₁, k ₂, k ₃

Radial coefficients in Brown distortion model

p ₁, p ₂

Tangential coefficients in Brown distortion model

d_r

Radial distortion in Brown distortion model

d_t

Tangential distortion in Brown distortion model

t_d

Delay in capturing the centre of the image

t_r

Delay difference in capturing of the first row to the last one

t_e

The exposure time of each pixel

2. Statement of the problem

Camera distortions are categorised as geometric, temporal and photometric. Geometric distortions are those that cause (static) physical movement in pixel locations. They are mainly caused because cameras with lenses do not exactly obey the pinhole model (refer to Section 4). Temporal distortions are those that cause timeshifts in the capture time of pixels (e.g. time delay, rolling shutter, motion blur, interlacing etc.). Photometric distortions are distortions added to pixel colour and intensity (light level, noise, vignetting, saturation, etc.).

Here we should describe shaders to show why all the direct distortion equations (at least for geometric and temporal distortions) have to be inversed in this paper. The operation of converting the sample data to realistic images, usually performed on a graphical processing unit (GPU), is called ‘rendering’ and the program doing that is called a ‘shader’. Shaders are available in two major subgroups, vertex and fragment. These are both parallel units of execution in GPU. Vertex shaders output corner positions (vertices) which form a polygon. This polygon is filled with pixels coloured by fragment shaders. In short, vertex shaders determine the geometry and fragment shaders determine the colour of pixels. As a polygon contains much higher numbers of surface pixels than corners, fragment shaders are built more efficiently in parallel execution.

This paper proposes a two-stage shader. The first shader renders an ‘ideal’ snapshot of the environment as a reference. It contains colour and depth channels. Notice that this is a standard shader and can be replaced by any sophisticated off-the-shelf ideal camera simulator with depth channel, like Gazebo or AirSim. To cause distortions, a connection is needed between each pixel in the final distorted image (output) and the ideally rendered pixel from the previous stage (input). The colour and depth values in the distorted image are transferred based on this relation in a second shader. This second stage can utilise two approaches. One direct approach is to use a vertex shader via every snapshot pixel, determine its new position after distortion, and draw on it (Figure 1). Pixel location in ideal image (1), its depth, camera rotational speed (2) and linear speed (3) along with the timeshift (4) are input so that the integrator determines each pixel's three-dimensional (3D) position after temporal distortion. Lens distortion is added afterwards using the lens zoom (5). Although it is straightforward, there is no guarantee that every pixel in the output image will be filled, as many input pixels would map to the same output pixel. Modelling the rolling shutter and motion blur is also hard in this way. In the case of rolling shutter cameras, the exact timeshift based on the final (distorted) pixel location is not known beforehand. An iterative approach can be used, but it involves repeating whole calculations. Modelling motion blur also needs another repeat and averaging on whole calculations. This approach is not exactly parallel as many points may incorporate in forming a single pixel. Finally, by using vertex shaders per pixel on a GPU, this approach is not as parallel and fast as a fragment shader.

Figure 1. Direct approach, adding geometrical distortions to an ideal frame

The other solution, the inverse approach, is to consider each output pixel position and reach its colour in the ideal snapshot by removing the distortion (Figure 2). Using the distorted pixel location (1) and lens zoom (5), the geometric distortion is removed first. By doing so, the geometric and temporal distortions are decoupled. Assuming the timeshift does not affect the depth, the 3D point position is calculated. This along with camera speeds (2, 3) and timeshift (4) is input to the inverse rolling shutter model to reach the ideal coordinates of the point. It gives a smooth image as it guarantees that every pixel in the output would have a colour (at least for those which map to valid pixels in the input image). Timeshift caused by the rolling shutter is a function of distorted pixel position which is known exactly. The motion blur can also be modelled efficiently as only a fraction of the algorithm is recalculated for getting the timeshift samples (inside the inverse rolling shutter model). This approach is completely parallel as the value of each output pixel can be found independently. On a GPU, a fragment shader can be used conveniently to render pixels in parallel. The only difficulty is that all the distortion models have to be inversed (analytically or numerically). This forms the root of some of the problems tackled in this paper.

Figure 2. The inverse approach, removing distortions from a distorted frame (inverse pipeline)

3. Aircraft track (flight path) generation

Aircraft (angular and linear) positions are directly used in the frame generation process. To model rolling shutter and motion blur, aircraft (angular and linear) velocities are also needed. The camera model needs this information in its own frequency. Unwanted discontinuities in position or velocities produce jumps in consequent frames and are not allowed. Out of range velocities would also cause extremely distorted or torn frames.

We assume the flight path of an aircraft (track) is provided by a flight log, which is a sparse table of aircraft positions and angles.Footnote ³ This approach is more flexible than using a dynamic flight simulation, as in Atashgah and Malaek (Reference Atashgah and Malaek2013), because it is not limited to a fixed type of aircraft. A ‘good’ trajectory would pass provided points as close as possible, maintaining a continuous track both in position and velocity, while fulfilling maximum speed and acceleration constraints. A simple linear interpolation between provided points to obtain data in camera frequency creates a velocity discontinuity at the provided points which is not applicable. A simple differentiation on the interpolated data to obtain velocities is also not applicable due to the noise present in the input data.

Several works have been done to generate acceptable trajectories. In Yan et al. (Reference Warren, McKinnon, He, Glover, Shiel, Yoshida and Tadokoro2015) a Kalman filter and cubic spline interpolation was used to generate realistic data from data logs. In Pham and Suh (2019) a spline method was used to generate data for a foot-mounted inertial measurement unit (IMU). In Irmisch et al. (Reference Irmisch, Baumbach, Ernst and Börner2019) a spline method was used to construct an IMU/camera simulation framework. Most of these approaches are based on a spline approach (Buckley, Reference Buckley1994; Schumaker, Reference Pueyo, Cristofalo, Montijano and Schwager2007). Spline-based trajectory generation is focused on passing (near) the provided points with continuous curves and continuous velocities/accelerations (Cⁿ continuity), not giving a velocity/acceleration constrained motion. As the spline uses a local set of points, it is hard for a spline to consider long-term lags of the trajectory.

One way to generate a ‘good’ trajectory is to see the problem as a plant-control system. A simplified model of an aircraft is controlled on the interpolated flight path with eligible inputs (accelerations, velocities). The aircraft is simulated by a six degrees of freedom (6-DOF) decoupled model. Each DOF is considered individually using two integrators and is controlled by acceleration (Figure 3, middle). This method gives the flexibility to determine maximum velocity and acceleration for each DOF individually. The two-integrator model ensures continuity of position and velocities. The aircraft will pass flight log points only if they are reachable with logical acceleration and velocities. If noisy GPS data generates a highly rugged path, the aircraft will go through an averaged route among the points. This solution is not dependent on the number of points and spacing between them. It has the advantage of working with highly distant data-log points which gives the ability to design the path instead of collecting real flight data.

Figure 3. The plant-control system used for path tracking

Two inner controller loops are used for each channel to keep the aircraft on course. Each loop uses a SQRT controller, inspired by a diagram in Copter Attitude Control (2016). Although not described in any paper, the SQRT controller is a useful, nonlinear, (near) minimum time controller. It is inspired by a car driving in traffic. A hurrying driver tries to reach the bumper of the car in front with maximum constant braking. It is selected as it allows one to reach the target with constrained acceleration and velocities and no overshoot. Constraining accelerations and velocities are derived from the object properties. SQRT gain gives the velocity derived from the well-known constant acceleration formula, Equation (1). In this formula, x ₁ is the input error, v ₁ is the desirable relative velocity to the target, x ₂ is the position at the target (set to zero), v ₂ is the relative speed at the target (also set to zero), and a _max is the maximum allowed acceleration. This finally gives a square root form of the desired velocity with respect to the error known as SQRT gain, Equation (2).

(1)\begin{align}{v_2}^2 - {v_1}^2 &= 2{a_{\max }}({x_2} - {x_1})\end{align}

(2)\begin{align}v &= \sqrt {2{a_{\max }}x}\end{align}

Note that due to the infinite slope of the gain at the origin, this controller may chatter near the target. To overcome this problem, a small linear region is added to the centre. As this velocity is relative to the target, the final desired velocity is formed by simply adding target velocity (as a feed-forward). Velocity saturation is finally applied to constrain velocities (Figure 4).

Figure 4. The complete SQRT controller structure

Besides position and velocity, this controller can virtually be used for controlling any variable, by giving its desirable derivative. After getting the aircraft's position, angles and linear/rotational velocities, those for the camera are easily obtained. Assuming the aircraft's distance to objects is much larger than the distance between the aircraft's centre of mass and the camera position, the camera position is found:

(3)\begin{equation}\begin{aligned} {R_{\textrm{Cam}}}{|_{\textrm{NED}}} & ={R_{\textrm{AC}}}{|_{\textrm{NED}}} + {R_{\textrm{Cam/AC}}}{|_{\textrm{NED}}}\\ & ={R_{\textrm{AC}}}{|_{\textrm{NED}}} + C_{\textrm{AC}}^{\textrm{NED}}{R_{\textrm{Cam/AC}}}{|_{\textrm{AC}}} \approx {R_{\textrm{AC}}}{|_{\textrm{NED}}} \end{aligned}\end{equation}

Camera velocities (in camera coordinates) are also simplified and determined:

(4)\begin{equation}\begin{aligned} { {{V_{\textrm{Cam}}}} |_{\textrm{Cam}}} = { {{V_{\textrm{AC}}}} |_{\textrm{Cam}}} + {W_{\textrm{AC/NED}}} \times {R_{\textrm{Cam/AC}}}\\ \approx { {{V_{\textrm{AC}}}} |_{\textrm{Cam}}} = { {C_{\textrm{AC}}^{\textrm{Cam}}C_{\textrm{NED}}^{\textrm{AC}}{V_{\textrm{AC}}}} |_{\textrm{NED}}} \end{aligned}\end{equation}

Camera rotational velocity is determined in Equation (5). The camera is assumed to be fixed on the aircraft body.

(5)\begin{equation}\begin{aligned} { {{W_{\textrm{Cam/NED}}}} |_{\textrm{Cam}}} & ={ {{W_{\textrm{AC}}}} |_{\textrm{Cam}}} + { {{W_{\textrm{Cam/AC}}}} |_{\textrm{Cam}}}\\ & =C_{\textrm{AC}}^{\textrm{Cam}}{ {{W_{\textrm{AC/NED}}}} |_{\textrm{AC}}} + { {{W_{\textrm{Cam/AC}}}} |_{\textrm{Cam}}} \approx C_{\textrm{AC}}^{\textrm{Cam}}{ {{W_{\textrm{AC/NED}}}} |_{\textrm{AC}}} \end{aligned}\end{equation}

4. Ideal camera simulation

An ideal camera is usually interpreted as one with a pinhole camera model, as in Equation (6). An assumption of brightness constancy is also made for simplicity. This means that every point has the same brightness (and colour) when looked at from any angle.

(6)\begin{equation}\begin{aligned} & h \left( {\left[ {\begin{array}{*{20}{c}} x\\ y\\ z \end{array}} \right]} \right) = h\left( {\left[ {\begin{array}{*{20}{c}} u\\ v \end{array}} \right]} \right) = \left[ {\begin{array}{*{20}{c}} {{o_x}}\\ {{o_y}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {{f_x}} & 0 \\ 0 &{{f_y}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} u\\ v \end{array}} \right]\\ & \left( {u: = {{x \over z}},v: = {{y \over z}}} \right) \end{aligned}\end{equation}

In this equation, the final pixel location h is a function of a point's location in camera coordinates x, y, z. This location is first normalised to u, v values. The ideal pinhole camera model has the parameters f_x, f_y, o_x and o_y. These are the focal length and pixel shifts. For ideal camera simulation, a custom OpenSceneGraph (OpenSceneGraph Library, Reference Oettershagen, Stastny, Mantel, Melzer, Rudin, Gohl, Agamennoni, Alexis, Siegwart, Wettergreen and Barfoot1998) simulation is utilised. It consists of scenery and camera nodes.

To make a realistic camera frame, colour and depth information are needed. Without loss of generality, we use digital satellite imagery available online from Bing maps (2005)Footnote ⁴ and 1 arcsec digital elevation map from Farr et al. (Reference Farr, Rosen, Caro, Crippen, Duren, Hensley, Kobrick, Paller, Rodriguez and Roth2007) as geographic data samples in this work. Sample data are meant to be as real as possible, just to allow a comparison with real captured images.

Two synchronised co-centred cameras are placed on the scene, using camera position and angle data generated by Equations (2)–(4). One camera gives a standard colour frame, while the other gives the depth frame. The depth camera has a modified shader to render distance of points to the camera plane (depth) instead of colours in full resolution (24 bits). This ideal full resolution depth and colour frames are later used in the distortion adding process. As adding geometrical distortions (like lens effect or rolling shutter) may need visual data from outside the final field of view (FOV), a bigger FOV (by a threshold) is rendered. These frame outputs are rendered in GPU but not shown, they are redirected to GPU memories (known as textures) for further processing (Figure 5).

Figure 5. Ideal and bigger rendered map of FSR 2015-ASL dataset, the first image in the lawnmower sequence: (a) colour rendering (b) 24 bit depth rendering. Notice that horizontal lines in depth map are blue channel (3rd byte) discontinuities

5. Geometric distortions simulation

Geometric distortions are deviations from the pinhole camera model caused by the lens. The most common lens distortion model used is the Brown five parameter model (Brown, Reference Brown1966), describing radial and tangential distortion in the pixels:

(7)\begin{equation}h\left( {\left[ {\begin{array}{*{20}{c}} u\\ v \end{array}} \right]} \right) = \left[ {\begin{array}{*{20}{c}} {{o_x}}\\ {{o_y}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {{f_x}}&0\\ 0&{{f_y}} \end{array}} \right]\left( {{d_r}\left[ {\begin{array}{*{20}{c}} u\\ v \end{array}} \right] + {d_t}} \right)\end{equation}

In this equation, radial and tangential distortions (d_r and d_t) are defined as:

(8)\begin{align} r&: = {u^2} + {v^2}\notag\\ {d_r}&: = (1 + {k_1}{r^2} + {k_2}{r^4} + {k_3}{r^6}) \end{align}

(9)\begin{align}{d_t}&: = \left[ \begin{array}{c} {2uv{p_1} + (r + 2{u^2}){p_2}}\\ {2uv{p_2} + (r + 2{v^2}){p_1}} \end{array} \right]\end{align}

Brown radial distortion parameters are k ₁, k ₂ and k ₃ which together form the radial distortion d_r. Brown tangential distortion parameters p ₁ and p ₂ form the tangential part d_t. The biggest problem in utilising this model is that its inverse is required to model the distortion caused by the lens, in a fragment shader. Unfortunately, this model is known to have no closed-form inverse. Even for the radial only problem which is a degree 5 or 7 polynomial, the Abel-Ruffini theorem states that a finite general solution does not exist. The usual way to undistort or distort a sequence of images is to pre-compute the transformation map and then warp the image using this map (Bing maps, 2005; Rehder et al., Reference Shelley2016). Although it is fast, this is only possible if the distortion coefficients are fixed (i.e., zoom is fixed in-the-run or distortion coefficients remain constant while zooming).

To solve the problem online, an iterative Gauss–Newton inverse solution is used, like the one in Rehder et al. (Reference Shelley2016). It is simple enough to be run on a GPU fragment shader online. Although there is no proof of convergence, our experiments show that it generally works well in practice. Non-converged regions are rendered in red colour for clarity (as in Figures 11 and 12). These regions are usually caused by using inappropriate calibration patterns (e.g., see Figure 12). Notice that for a fixed-lens camera, the calculations are repeated in every frame and all the pixels remain converged if they have converged once. In the test case of a badly installed fisheye lens (Section 8.3), it converges in less than four iterations for all pixels within a 0⋅1-pixel accuracy.

6. Temporal distortion model

Temporal distortions are those that cause the displacement of pixels caused by ‘small’ timeshifts. Although named differently, image delay, motion blur, rolling shutter and interlacing are all caused by timeshifts and obey the same formulations.

Several works consider adding motion blur to still (ideal) images. This paper is limited to those that incorporate the depth of pixels (not consecutive frames) as an input. For exact timeshift modelling, one should render frames in every timeshift (e.g., every row for rolling shutter) and then mix them, which is a non-real-time solution. Rendering a frame in this way may take up to tens of seconds (Huang et al., Reference Huang, Hou, Ren and Zhou2012; Zhang et al., Reference Yu, Cai and Chen2014). An overview of motion blur generating methods is presented in Navarro et al. (Reference Mueller, Smith and Ghanem2011). In Atashgah and Malaek (Reference Atashgah and Malaek2013) and Zhao et al. (Reference Zhang, Wang, Sun and Han2014) a model is developed to simulate earth surface motion blur as a function of the aircraft's rotational and linear velocities and the focal length. These models assume the pixel travel as linear in every frame, which causes deformations in in-plane camera rotations (discussed in Section 8.4). Although it is not explicitly stated, the model in Atashgah and Malaek (Reference Atashgah and Malaek2013) seems to run online. In theory, rolling shutter and motion blur are both timeshifts and can be produced in the same way. A unified model for modelling motion blur and rolling shutter simultaneously was first introduced by Meilland et al. (Reference Mahmoudi2013). This model used mixing multiple renders to generate motion blur and was offline. In Guertin et al. (Reference Guertin, McGuire and Nowrouzezahrai2014) the problem of real-time motion blur is discussed along with the removal of artifacts over edges. In Gribel and Akenine-Möller (Reference Gribel and Akenine-Möller2017) interactions between ray-tracing and motion blur were discussed. None of these works consider the lens effect combined with general temporal distortions (delay, blur, rolling shutter and interlacing).

Exploiting the simulated depth of pixels, it is possible to develop a short-term time displacement model for every pixel. The model predicts source pixel displacements in the ideally rendered image. It utilises a single colour + depth frame and camera velocities. In this model, both linear and angular velocities are assumed to be fixed during the timeshift. Since the maximum rolling shutter and motion blur timeshift is in the order of tens of milliseconds, this seems acceptable, except at extreme linear-rotational accelerations. As image delay might be quite large (in the order of a tenth of a second), it is not recommended to model it by this model. Fixed image time delay can easily be added to the whole images by buffering the camera path.

Assuming small linear pixel-wise travel, the camera model is linearised to achieve a short-term displacement model:

(10)\begin{equation}\frac{{dh}}{{dt}}(u,v) = \frac{{\partial h}}{{\partial u}}\frac{{\partial u}}{{\partial t}} + \frac{{\partial h}}{{\partial v}}\frac{{\partial v}}{{\partial t}} = \frac{{\partial h}}{{\partial u}}\dot{u} + \frac{{\partial h}}{{\partial v}}\dot{v}\end{equation}

In this equation, time derivatives of normalised pixel coordinates are calculated as follows in Equations (11) and (12).

(11)\begin{equation}\begin{aligned} \dot{u} & =\dfrac{{\partial ({\raise0.5ex\hbox{$\scriptstyle x$}\kern-0.1em/\kern-0.15em\lower0.25ex\hbox{$\scriptstyle z$}})}}{{\partial t}} = \dfrac{{z\dot{x} - x\dot{z}}}{{{z^2}}} = \dfrac{{\dot{x}}}{z} - \dfrac{x}{z}\dfrac{{\dot{z}}}{z} = \dfrac{{{V_x}}}{z} - u\dfrac{{{V_z}}}{z}\\ \dot{v} & =\dfrac{{\partial ({\raise0.5ex\hbox{$\scriptstyle y$}\kern-0.1em/\kern-0.15em\lower0.25ex\hbox{$\scriptstyle z$}})}}{{\partial t}} = \dfrac{{{V_y}}}{z} - v\dfrac{{{V_z}}}{z} \end{aligned}\end{equation}

Note that z is pixel depth and is assumed to be constant during the time period. As all points are assumed to be stationary, point velocities in camera frame (V_x, V_y and V_z) are just a function of camera velocity and rotational speeds.

(12)\begin{equation}{ {{V_{\textrm{point/Cam}}}} |_{\textrm{Cam}}} ={-} { {{V_{\textrm{Cam/Inertia}}}} |_{\textrm{Cam}}} - { {{W_{\textrm{Cam}}}} |_{\textrm{Cam}}} \times { {{R_{\textrm{point}}}} |_{\textrm{Cam}}}\end{equation}

As the lens distortion (d_r and d_t) is already removed in the reverse pipeline, derivatives of h can be derived easily.

(13)\begin{equation}\frac{{\partial h}}{{\partial u}} = \left[ {\begin{array}{*{20}{c}} {{f_x}}\\ 0 \end{array}} \right],\textrm{ }\frac{{\partial h}}{{\partial v}} = \left[ {\begin{array}{*{20}{c}} 0\\ {{f_y}} \end{array}} \right]\end{equation}

It is interesting to note that the inverse temporal displacement model is independent of the distortion model. The final displacement caused by the time displacement is calculated as follows:

(14)\begin{equation}h(u(t),v(t)) \approx h(u({t_0}),v({t_0})) + \frac{{dh}}{{dt}}({u({t_0}),v({t_0}),{\raise0.7ex\hbox{${\vec{V}}$} \!\mathord{/ {\vphantom {{\vec{V}} z}} }\!\lower0.7ex\hbox{$z$}}({t_0})} )\Delta t\end{equation}

In Equation (14), V/z is the normalised velocity of a 3D point in and with respect to the camera coordinates. It is needed for calculation of u, v derivatives in Equation (11), and is assumed to be fixed during the small timeshift in the following equations. Equation (14) is simplified, for clarity, to Equation (15).

(15)\begin{equation}h(t) \approx h({t_0}) + \frac{{dh}}{{dt}}({{t_0}} )\Delta t\end{equation}

Note that the inverse is needed in the reverse pipeline to get from the known distorted coordinates h(t) to the undistorted coordinates h(t ₀) in Equation (16).

(16)\begin{equation}h({t_0}) \approx h(t) - \frac{{dh}}{{dt}}({t_0})\Delta t \approx h(t) - \frac{{dh}}{{dt}}(t)\Delta t\end{equation}

This model is called the (inverse) linear displacement model. By redefining the derivative as negative, this equation finally becomes a standard differential equation.

(17)\begin{equation}\begin{split} &h({t_0}) \approx h(t) + \dfrac{{d{h^ - }}}{{dt}}(t)\Delta t\\ &\left( {\dfrac{{d{h^ - }}}{{dt}}: ={-} \dfrac{{dh}}{{dt}}} \right) \end{split}\end{equation}

Note that linearisation might cause some errors in case of high rotations, as seen in Figure 13. To improve the simulation's accuracy, one can implement a non-standard backward two-point Runge-Kutta order 2 (RK2) solver. The idea of RK2 is to use the average derivatives at the initial point and at the first estimate solution point. Based on this idea, the proposed solution is described in Equation (18).

(18)\begin{equation}\begin{split} {h_1}({t_0}) & =h(t) + \dfrac{{d{h^ - }}}{{dt}}(t)\Delta t\\ {h_2}({t_0}) & =h(t) + \dfrac{{(d{h^ - }/dt)(t) + (d{h^ - }/dt)({t_1},{h_1})}}{2}\Delta t \end{split}\end{equation}

To further improve the RK2 accuracy, the t−t ₀ timebase is divided into multiple time-steps and the solver is run multiple times (Figure 13).

6.1 Rolling shutter modelling

Rolling shutter cameras have a timeshift in subsequent image rows (or columns) (Figure 6). When an image is captured, a timestamp is provided (t_i | i = 0, 1, 2 …). The middle row of the image was captured at some time t_i + t_d. Notice that t_d is usually negative. The entire image of N rows has a read time of t_r, meaning that row k was captured at t_i + t_d + k.t_r/N where $k \in [{\textstyle{{ - N} \over 2}},{\textstyle{N \over 2}}]$ (Shelley, Reference Shah, Dey, Lovett and Kapoor2014). The timeshift is thus assumed to be

(19)\begin{equation}\Delta t = t - {t_i} = {t_d} + k\frac{{{t_r}}}{N}\end{equation}

Figure 6. Timeshift caused by rolling shutter

The final h is the one without the lens and rolling shutter distortions, so its colour and depth can be read easily from the ideal image rendered earlier.

6.2 Motion blur modelling

For modelling motion blur, an exposure time of t_e is considered for each pixel in each image. Using the displacement model, two sample pixel positions (h points) using the beginning and end of this exposure time are generated.

(20)\begin{equation}\begin{split} \Delta {t_1} & ={t_d} + k\dfrac{{{t_r}}}{N} - \dfrac{{{t_e}}}{2}\\ \Delta {t_2} & ={t_d} + k\dfrac{{{t_r}}}{N} + \dfrac{{{t_e}}}{2} \end{split}\end{equation}

These points form a sampling line in the ideally rendered frame. Colours on this sampling line are averaged to form the final colour caused by motion blur.

6.3 Interlacing

For better visual feeling, analog videos come in double frequencies and half rows. At the time of visualisation (or in the analog to digital device), two consecutive odd and even frames are mixed to construct a full frame (Figure 7). Interlacing can cause serious problems in image processing, much worse than rolling shutter. This operation is reversible only if no image loss due to compression occurs in the analog to digital device, otherwise a wavy image on the edge of moving parts is generated. The timeshift model for the interlacing is demonstrated in Equation (21).

(21)\begin{equation}\Delta t = {t_d} + k\frac{{{t_r}}}{N} \pm \frac{{{t_e}}}{2} - \frac{{k\%2}}{{2fps}}\end{equation}

Figure 7. Interlacing combines two consequent half-row frames like a zipper, to form a full-frame one

7. Photometric distortion simulation

Photometric distortions cause changes in pixel brightness. They are static (e.g., vignetting) or dynamic (noise, light level etc.). These imperfections are implemented for each output pixel after geometrical distortions. As they are independent, they are easily parallelised (in the pixel-wise fragment shader).

Vignetting is the reduction of the image's brightness or total saturation towards the edges of the frame. A vignetting frame is obtained by taking a picture from a white wall in full light (Figure 8). It is sent to the GPU shader, thresholded to eliminate noise in central image areas and finally multiplied by the brightness of pixels in the shader.

Figure 8. Vignetting map of a badly installed fisheye lens on Galaxy S3 main camera

Although simple shadowing/ray-tracing methods exist, they interfere with the existing shadows/light gradient on the satellite images. In this study, the light source, the reflectance of light, and shades as they already exist are ignored. The environment light level is just modelled by a gain. An additive Gaussian sensor noise is also implemented. Finally, saturation is used to keep the pixel values in the 0 to 1 range (Figure 9).

Figure 9. Pixel-wise photometric model, including vignetting, colour noise, light level and saturation

8.6 Rolling shutter model test

Motion blur is evaluated qualitatively. The rotation effect on the rolling shutter is demonstrated in Figure 18. The camera is rotating clockwise (CW), so the image rotates counter clockwise (CCW). The top rows are simulated as captured sooner, so they are rotated less CCW. The bottom rows are rotated more CCW.

Figure 18. Rolling shutter under camera rotation: (a) static camera, (b) camera rotating CW

Under camera motion, the rolling shutter creates weird effects. Shearing, shortening and elongation are seen in Figure 19. Edge artifacts are still seen in the left of shear in Figure 19(b) and top of elongation in Figure 19(d).

Figure 19. Rolling shutter under camera motion: (a) static camera, (b) camera moving to the right (shear), (c) camera moving upwards (shortening), (d) camera moving downwards (elongation)

8.7 Interlacing model test

As seen in Figure 20, the interlacing effect is much like the motion blur with two sample points and can be evaluated the same way by comparing image overlays. The only difference is that, unlike motion blur, the delay is switching in odd and even rows.

Figure 20. Interlacing under camera rotation: (a) and (b) static camera in two positions, (c) simulated interlaced image at the second position

8.8 Time displacement on map

A video frame with a 640 × 480 dimension is generated on the lawnmower sequence in Figure 21(a). The camera's view is assumed to be exactly parallel to the nose of the aircraft. The aircraft/camera is rotating CCW and to the left. For visualisation purposes, all the distortion parameters are set with large values. A fisheye lens with high distortion levels is used (Table 2). All the time constants are exaggerated to 1 s, see Figure 21(b)–(e). In Figure 21(b) a 1 s delay is simulated, compensated by a 1 s timeshift. In Figure 21(c) rolling shutter is exaggerated. The up row belongs to half a second sooner and the down row half a second later. The image is shortened in the left and elongated in right as expected. Exposure time of 1 s causes extreme motion blur as expected, see Figure 21(d). Pure interlacing creates a two-fold image, see Figure 21(e). Notice that individual even and odd rows belong to different timestamps with 1 s spacing.

Figure 21. 640 × 480 fisheye camera frame added to FSR 2015-ASL dataset, using exaggerated temporal distortions: (a) no temporal distortions, (b) 1 s later +1 s delay (visually indistinguishable from (a)), (c) 1 s rolling shutter (t_r) (d) 1 s exposure time (t_e), (e) 1 s spacing for interlacing (t_i)