---
title: Fast Inverse Georectification
url: https://www.emergentmind.com/topics/fast-inverse-georectification-technique
type: topic
---

# Fast Inverse Georectification

A fast inverse georectification technique enables hyperspectral remote sensing data to be rapidly and artifact-free mapped from sensor coordinates to geospatially consistent representations. Unlike conventional direct georectification, which is computationally intensive and prone to pixel coverage artifacts due to per-pixel ray tracing, the fast inverse method, as implemented in Spacecube, leverages GPU rasterization and mesh-based interpolation to achieve real-time, high-fidelity georectification suitable for interactive viewing, analysis, and batch export in hyperspectral imaging [2601.05181].

## 1. Mathematical Formulation

The georectification process involves transformations between three coordinate systems:

- **Image pixel coordinates** $(s, \ell)$: $s$ is the sample index (pixel column), $\ell$ is the line index (acquisition time/flight direction).
- **Sensor (camera) coordinates** $(X_c, Y_c, Z_c)$: with the origin at the optical center and $Z_c$ aligned with the optical axis.
- **Geographic (NED) coordinates** $(N, E, D)$: North-East-Down metric system derived from UTM and a specified ground height $H$.

### Forward Georectification

Given $(s, \ell)$:
1. Normalize pixel positions: $u = \frac{s - c_x}{f_x}$, $v = \frac{\ell - c_y}{f_y}$.
2. Account for lens distortion (radial), derive direction $\mathbf{d}_c$ in camera space.
3. Apply extrinsic rotation $R(\ell)$ and translation $t(\ell)$ for each line.
4. Intersect the resulting world-space ray with the ground plane $D=0$ to obtain ground coordinates $(N, E)$ with
   $\lambda = -D_c / D_d$.

### Inverse Georectification

Mapping $(N, E)$ to $(s, \ell)$ directly would require per-pixel inversion of the forward process, which is computationally prohibitive. Spacecube constructs a 2D mesh covering the imaged footprint for each input line. Each mesh vertex is attributed with both its ground position and its corresponding $(s, \ell)$. Rasterizing this mesh on the GPU yields, for every output pixel, the correct sample and line coordinates via perspective-correct barycentric interpolation, thus obviating the need for explicit per-pixel inversion.

Given triangle vertices $V_0, V_1, V_2$ at positions $X_i$ with sample/line $(s_i, \ell_i)$, each ground pixel $X$ is represented by
$$
X = \alpha_0 X_0 + \alpha_1 X_1 + \alpha_2 X_2
$$
and the sensor coordinates recover as
$$
(s, \ell) = \alpha_0 (s_0, \ell_0) + \alpha_1 (s_1, \ell_1) + \alpha_2 (s_2, \ell_2),
$$
with $(\alpha_0, \alpha_1, \alpha_2)$ the barycentric coordinates of $X$ within the triangle.

## 2. GPU/OpenGL Implementation

### Data Structures

Each mesh vertex stores:
- Projected ground position $(N, E)$
- Depth (vertical distance for perspective correction)
- Sensor coordinates $(s, \ell)$

Input data for georectification includes:
- 2D textures of hyperspectral bands ($N_s \times N_\ell$ per band)
- Calibration textures: dark offsets, radiometric gains, line response

### Shaders

**Vertex Shader** maps vertices from ground coordinates to device coordinates, passing $(s, \ell)$ through to the fragment shader using linear interpolation variables.

**Fragment Shader** performs:
1. Interpolation (with hardware perspective correction) of $(s, \ell)$.
2. Bilinear sampling of texture data per band.
3. Calibration adjustment: $L_{b, \ell, s} = (DN_{b, \ell, s} - \text{dark}_{b, s}) \cdot \text{gain}_{b, s} \cdot \text{resp}_\ell$.

This GPU pipeline enables all mesh triangles to be drawn in a single call per cube, with all output pixel computations executed in parallel and no explicit loops over input lines.

### Example Shader Skeletons

**Vertex Shader (GLSL-style pseudocode)**:
```glsl
#version 330 core
layout(location=0) in vec2  a_NE;      // world coords
layout(location=1) in float aDepth;    // used for perspective
layout(location=2) in float aSamp;     // sample index
layout(location=3) in float aLine;     // line index

out float vSamp;
out float vLine;

uniform vec2 uCenterNE;
uniform vec2 uScaleNE;
void main() {
  vec2 norm = (a_NE - uCenterNE) / uScaleNE;
  gl_Position = vec4(norm.xy, -aDepth, 1.0);
  vSamp = aSamp;
  vLine = aLine;
}
```

**Fragment Shader (GLSL-style pseudocode)**:
```glsl
#version 330 core
in float vSamp;
in float vLine;

uniform sampler2D uBandTex[3];
uniform sampler2D uDarkTex;
uniform sampler2D uRadTex;
uniform sampler2D uRespTex;

out vec3 outColor;
void main() {
  float t_s = (vSamp + 0.5) / Nsamples;
  float t_l = (vLine + 0.5) / Nlines;
  float val[3];
  for(int c=0; c<3; ++c)
    val[c] = texture(uBandTex[c], vec2(t_s,t_l)).r;
  vec2 dc = texture(uDarkTex, vec2(t_s,0.5)).rg;
  vec2 rg = texture(uRadTex,  vec2(t_s,0.5)).rg;
  float resp = texture(uRespTex, vec2(t_l,0.5)).r;
  for(int c=0; c<3; ++c)
    val[c] = (val[c] - dc[c]) * rg[c] * resp;
  outColor = vec3(val[0], val[1], val[2]);
}
```

## 3. Interpolation, Artifacts, and Error Properties

The GPU pipeline relies on two levels of hardware interpolation:
1. **Perspective-correct barycentric interpolation** of $(s, \ell)$ within mesh triangles ensures that every covered pixel receives exactly one projected sample coordinate, eliminating undersampling and gaps.
2. **Bilinear interpolation** within texture fetch units for band data provides smoothing and mitigates artifacts across sensor coordinates.

Due to the mesh’s construction and the triangle coverage, every output pixel mapped to the surveyed ground area is hit by exactly one triangle, guaranteeing that the output image is contiguous and free from holes—an advantage over forward-mapping methods that require post-hoc interpolation to fill pixel gaps, or which may introduce over-write artifacts. The only significant geometric error arises from the assumption of local planarity in the ground model, which empirical measurements show to remain below 1% of a sample width even with significant roll/pitch variations [2601.05181].

## 4. Performance Evaluation

Benchmarking on a workstation with a Xeon E-2176M and a Quadro P2000 GPU (64GB RAM, data on RAM-disk):

| Method                | Wall-Clock Time   | Real-Time Factor  |
|-----------------------|------------------|-------------------|
| Spacecube (GPU inv.)  | 11 s             | ×1400 real time   |
| Spectronon (CPU dir.) | 21 m 29 s        | ×0.12 real time   |
| PARGE (CPU dir.)      | ~35 m 48 s       | ×0.075 real time  |

For a dataset of 40 cubes ($900 \times 1000 \times 300$ bands), Spacecube processes at approximately 200 Mpixels/s (3 bands) on a midrange GPU, achieving over 100× speedup and more than 10× real-time throughput compared to typical CPU-based direct georectification pipelines [2601.05181]. Doubling the number of bands doubles runtime linearly; quadrupling output pixel count (i.e., halving GSD) increases load as expected but remains GPU-bound. Even on consumer GPUs (e.g., GTX 1660 Ti), throughput sufficient for >1000% of real-time is observed for standard cube sizes.

## 5. Workflow Pipeline and Practical Usage

The complete pipeline includes:

A. **Data ingestion and calibration**:
   - Reading ENVI cubes.
   - Loading sensor calibration (dark offsets, radiometric gains, line exposure) from metadata.
   - Converting raw digital numbers to radiance (or reflectance) via $L_{b, \ell, s} = (DN_{b, \ell, s} - \text{dark}_{b, s}) \cdot \text{gain}_{b, s} \cdot \text{resp}_\ell$.

B. **Georectification setup**:
   - Converting INS/GNSS sensor timestamps and poses to per-line extrinsics $(R(\ell), t(\ell))$ in UTM/NED.
   - Interpolating these to line start times for precise ground mapping.

C. **Ground mesh generation**:
   - For each line $\ell$, project edge sensor rays to ground to define line endpoints.
   - Connect endpoints as a triangle strip mesh covering the imaged ground region, with each vertex tagged by $(s, \ell)$ and $(N, E)$.

D. **Rasterization (Rendering)**:
   - Uploading mesh and texture data to GPU once per cube.
   - Adjusting view, scale, or selected bands interactively or in batch.
   - Drawing the mesh to georectified output; results are either displayed or saved.

### User-Configurable Parameters

| Parameter     | Effect                                 | Notes                                      |
|---------------|----------------------------------------|--------------------------------------------|
| Ground height | Shifts mesh, affects intersection $\lambda$ | Errors create parallax in overlaps         |
| GSD           | Reduces/increases output resolution    | Coarser GSD increases speed                |
| Calibration   | Toggles dark/gain multiplication       | Small GPU cost for switching calibration   |
| Band selection| Controls number of bands exported      | More bands = more passes (linear scaling)  |
| Mipmapping    | Enhances performance when zoomed out   | Can cause seams at cube borders            |

## 6. Algorithmic Advantages and Limitations

This fast inverse georectification pipeline recasts the problem into a memory-bandwidth-limited, fully parallel process. It delivers gap-free georectified outputs, deterministic coverage, and high throughput by exploiting hardware rasterization and interpolation. The primary limitation is the flat-ground assumption, but the geometric errors introduced remain acceptably low for most low-altitude remote sensing operations. *A plausible implication is* that this paradigm directly supports real-time hyperspectral data acquisition and interactive analytics in agricultural and small-UAS contexts.

## 7. Conclusion

Spacecube's fast inverse georectification technique demonstrates that recasting georectification as a GPU mesh rasterization problem delivers major improvements over CPU-based direct approaches. Mathematical mappings are unchanged in substance, yet inversion is achieved without brute-force computation, resulting in an operational system that is routinely over 100× faster, runs at >10× real time, and robustly prevents pixel coverage artifacts. The pipeline relies only on a local planar-ground model assumption, whose empirical errors are negligible in targeted applications [2601.05181].

Source: https://www.emergentmind.com/topics/fast-inverse-georectification-technique