Geometric corrections

From Agri4castWiki
Jump to: navigation, search

Chip matching procedure

The general principles of the method used are as follows:

  • Select a chip position (using a list of predefined chip positions or using a random generator).
  • Map the chip position to the reference window and retrieve the N1 x N1 chip window from the reference image.
  • Map the chip position to the warp image, giving (x0,y0), and search the position (x1,y1) in the N2 x N2 search window for which the correlation between the warp image and the chip image is maximal.
  • A match is found if certain conditions are fulfilled, e.g. the correlation value has to be above a certain threshold, the correlation peak has to be sharp enough, etc..
  • Repeat the above steps until N matching chips have been found and determine the polynomials fx(x,y) and fy(x,y) that map the (x1,y1) coordinates into (x0,y0) coordinates:
       fx(x1,y1) = x0
       fy(x1,y1) = y0
  • Generate the warped image using the polynomials fx and fy to determine for each pixel position of the warped image the corresponding ‘nearest neighbour’ pixel in the warp image. Only the tiles with a sufficient fill level are warped and the others are removed from the image.
accuracy=\frac{\sum \left[\left(x_{1} ,y_{1} \right)-\left(x_{0} ,y_{0} \right)\right] }{N}


The accuracy of the warp image is calculated from the (x0,y0) and (x1,y1) coordinates of the N matching chips:

accuracy=\frac{\sum \left[\left(x_{1} ,y_{1} \right)-\left(x_{0} ,y_{0} \right)\right] }{N}

The root mean square error provides a measure for the accuracy of the warping and is calculated as:

RMSE=\sqrt{\frac{\sum \left[\left(f_{x} \left(x_{1} ,y_{1} \right)-x_{0} \right)^{2} -\left(f_{y} \left(x_{1} ,y_{1} \right)-y_{0} \right)^{2} \right] }{N} }

The performance of the polynomials on the whole image depend largely on the spatial distribution of the matching chips. Tiling, i.e. the subdivision of the image in blocks or tiles, is introduced to guarantee a good distribution and, at the same time, to give an indication of the percentage of the warp image for which the polynomials can be considered ‘valid’.


The basic idea behind the tiling is that when we can find TN matching chips in every tile, then the distribution of chips over the warp image is perfect. However, because it is expected that there are no discontinuities in the evolution of the ‘errors’ over the warp image (otherwise our approach using polynomials for the warping would fail), it is not necessary to really find TN matching chips in every tile: it suffices that the distribution of the chips found is good.

This observation lead to the introduction of a ‘bonus’ system: when a matching chip is found in a tile, a ‘bonus’ is attributed to its neighbours and the amount of the ‘bonus’ is inversely proportional to the distance of the neighbour (see fig.). Hence, due to the introduction of the ‘bonus’ system each matching chip results in the attribution of V ‘virtual’ chips:

    V = 1 + ∑ bonus

The original goal can now be rephrased as follows: when TN ‘virtual’ matching chips are found in every tile, then their distribution over the warp image is perfect.

Given that N matching chips should be found and knowing that 1 matching chip results in (maximally) Vmax ‘virtual’ chips being attributed, then a feasible value for TN can be determined as follows:

    TN = (N * Vmax) / number of tiles

One can get an indication of the percentage of the warp image for which the polynomials can be considered ‘valid’ as follows:

    fill level = ∑ ‘virtual’ chips / (TN * number of tiles)

The number of tiles depends on the image size.

Tiles that are insufficiently covered with matching chips are excluded in the warping, because the polynomial is not valid for this part of the image.

When searching for the position (x1,y1) for which the correlation between the chip window and the search window is maximal, the chip window is always moved over the search window in steps of one pixel in x- and/or y-direction. As a result the position of the maximum is only determined with pixel precision. A specific sub-pixel matching step have been added in order to determine the position of the maximum with sub-pixel precision.


    -Drop from the 3x3 window centred on (x1,y1) the pixels with the lowest correlation values (with the added constraint that it is not allowed to drop a complete row or column)
    -Determine the coefficients for which the bivariate quadratic function s(x,y) passes through the 5 remaining pixels:
         s(x,y) = c0 + c1*x + c2*y + c3*x2 + c4*y2
    -Determine the position of the maximum of the surface s(x,y):
         xmax = ( -2 * c1 * c4 ) / ( 4 * c3 * c4 )
         ymax = ( -2 * c2 * c3 ) / ( 4 * c3 * c4 )

GeometricCorrection PerformanceVITOProcessingChain.png

Geometric correction – performance of the NOAA-AVHRR processing chain. From left to right: 1. result after application of the orbital model (SGP4), 2. result after chipmatching, 3. SPOT-VGT ref image

Go back to NOAA-AVHRR, METOP-AVHRR, SPOT-VEGETATION or TERRA-MODIS pre-processing description.

Quality of the geometric correction

There are two options to describe the quality of the geometric correction.

  • IMG with RMSE: for each segment that was geometrically corrected with the Chip matching procedure, a list is maintained with the main parameters, incl. the RMSE of the geometric correction (computed at the level of the polynomial regression). The segment_ID IMG is reclassified (see 10-daily composites) such that each pixel receives the RMSE of the segment from which it was extracted. From this, a quicklook is generated.
  • IMG with Chip_Deviation: In a second approach, the final NDVI-S10 is submitted to the same chip matching procedure as explained in Chip matching procedure. However, in this case a fixed set of evenly spaced control points is used (centre of a 20km x 20km grid, only the land pixels are kept). The same reference images are used as before (LTA NDVI-S30 from SPOT-VGT) and also the chip filtering is applied, but then the procedure is halted. Output is a list of all the retained (reliably matched) points and their bidirectional deviation from the reference. These deviations are plotted in a quicklook.

Both quicklooks provide different views on the geo-quality of the composites in a geographical context. It is not always feasible to compute these indicators. The second approach for instance requires the presence of suitable reference images. But for the analysis of MODIS at 250m resolution, such reference data are not immediately available.

Go back to NOAA-AVHRR, METOP-AVHRR, SPOT-VEGETATION or TERRA-MODIS pre-processing description.