Line data    Source code

       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  APP ENVISAT Support
       4             :  * Purpose:  GCPs Unwrapping for products crossing the WGS84 date-line
       5             :  * Author:   Martin Paces martin.paces@eox.at
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2013, EOX IT Services, GmbH
       9             :  *
      10             :  * SPDX-License-Identifier: MIT
      11             :  ****************************************************************************/
      12             : 
      13             : #include "gdal.h"
      14             : #include <cmath>
      15             : #include <cstdio>
      16             : 
      17             : // number of histogram bins (36 a 10dg)
      18             : constexpr int NBIN = 36;
      19             : // number of empty bins to guess the flip-point
      20             : constexpr int NEMPY = 7;
      21             : 
      22             : // WGS84 bounds
      23             : constexpr double XMIN = -180.0;
      24             : // constexpr double XMAX = 180.0;
      25             : constexpr double XDIF = 360.0;
      26             : constexpr double XCNT = 0.0;
      27             : 
      28             : // max. allowed longitude extent of the GCP set
      29             : constexpr double XLIM = XDIF * (1.0 - NEMPY * (1.0 / NBIN));
      30             : 
      31             : /* used by envisatdataset.cpp */
      32             : extern void EnvisatUnwrapGCPs(int cnt, GDAL_GCP *gcp);
      33             : 
      34             : // The algorithm is based on assumption that the unwrapped
      35             : // GCPs ('flipped' values) have smaller extent along the longitude.
      36             : // We further assume that the length of the striplines is limited
      37             : // to one orbit and does not exceeded given limit along the longitude,
      38             : // e.i., the wrapped-around coordinates have significantly larger
      39             : // extent the unwrapped. If the smaller extend exceeds the limit
      40             : // the original tiepoints are returned.
      41             : 
      42           0 : static double _suggest_flip_point(const int cnt, GDAL_GCP *gcp)
      43             : {
      44             :     // the histogram array - it is expected to fit the stack
      45             :     int hist[NBIN];
      46             : 
      47             :     // reset the histogram counters
      48           0 :     for (int i = 0; i < NBIN; i++)
      49           0 :         hist[i] = 0;
      50             : 
      51             :     // accumulate the histogram
      52           0 :     for (int i = 0; i < cnt; i++)
      53             :     {
      54           0 :         double x = (gcp[i].dfGCPX - XMIN) / XDIF;
      55           0 :         int idx = (int)(NBIN * (x - floor(x)));
      56             : 
      57             :         // The latitudes should lay in the +/-180 bounds
      58             :         // although it should never happen we check the outliers
      59           0 :         if (idx < 0)
      60           0 :             idx = 0;
      61           0 :         if (idx >= NBIN)
      62           0 :             idx = NBIN - 1;
      63             : 
      64           0 :         hist[idx] += 1;
      65             :     }
      66             : 
      67             :     // Find middle of at least NEMPTY consecutive empty bins and get its middle.
      68           0 :     int i0 = -1;
      69           0 :     int i1 = -1;
      70           0 :     int last_is_empty = 0;
      71           0 :     for (int i = 0; i < (2 * NBIN - 1); i++)
      72             :     {
      73           0 :         if (0 == hist[i % NBIN])  // empty
      74             :         {
      75           0 :             if (!last_is_empty)  // re-start counter
      76             :             {
      77           0 :                 i0 = i;
      78           0 :                 last_is_empty = 1;
      79             :             }
      80             :         }
      81             :         else  // non-empty
      82             :         {
      83           0 :             if (last_is_empty)
      84             :             {
      85           0 :                 i1 = i;
      86           0 :                 last_is_empty = 0;
      87             : 
      88             :                 // if the segment is long enough -> terminate
      89           0 :                 if ((i1 - i0) >= NEMPY)
      90           0 :                     break;
      91             :             }
      92             :         }
      93             :     }
      94             : 
      95             :     // if all full or all empty the returning default value
      96           0 :     if (i1 < 0)
      97           0 :         return XCNT;
      98             : 
      99             :     // return the flip-centre
     100             : 
     101           0 :     double tmp = ((i1 - i0) * 0.5 + i0) / ((float)NBIN);
     102             : 
     103           0 :     return (tmp - floor(tmp)) * XDIF + XMIN;
     104             : }
     105             : 
     106           0 : void EnvisatUnwrapGCPs(int cnt, GDAL_GCP *gcp)
     107             : {
     108           0 :     if (cnt < 1)
     109           0 :         return;
     110             : 
     111             :     // suggest right flip-point
     112           0 :     double x_flip = _suggest_flip_point(cnt, gcp);
     113             : 
     114             :     // Find the limits along the longitude (x) for flipped and unflipped values.
     115             : 
     116           0 :     int cnt_flip = 0;  // flipped values' counter
     117             :     double x0_dif, x1_dif;
     118             : 
     119             :     {
     120             :         double x0_min;
     121             :         double x0_max;
     122             :         double x1_min;
     123             :         double x1_max;
     124             : 
     125             :         {
     126           0 :             double x0 = gcp[0].dfGCPX;
     127           0 :             int flip = (x0 > x_flip);
     128           0 :             x0_min = x0;
     129           0 :             x0_max = x0;
     130           0 :             x1_min = x0 - flip * XDIF;
     131           0 :             x1_max = x1_min;
     132           0 :             cnt_flip += flip;  // count the flipped values
     133             :         }
     134             : 
     135           0 :         for (int i = 1; i < cnt; ++i)
     136             :         {
     137           0 :             double x0 = gcp[i].dfGCPX;
     138           0 :             int flip = (x0 > x_flip);
     139           0 :             double x1 = x0 - flip * XDIF;  // flipped value
     140           0 :             cnt_flip += flip;              // count the flipped values
     141             : 
     142           0 :             if (x0 > x0_max)
     143           0 :                 x0_max = x0;
     144           0 :             if (x0 < x0_min)
     145           0 :                 x0_min = x0;
     146           0 :             if (x1 > x1_max)
     147           0 :                 x1_max = x1;
     148           0 :             if (x1 < x1_min)
     149           0 :                 x1_min = x1;
     150             :         }
     151             : 
     152           0 :         x0_dif = x0_max - x0_min;
     153           0 :         x1_dif = x1_max - x1_min;
     154             :     }
     155             : 
     156             :     // in case all values either flipped or non-flipped
     157             :     // nothing is to be done
     158           0 :     if ((cnt_flip == 0) || (cnt_flip == cnt))
     159           0 :         return;
     160             : 
     161             :     // check whether we need to split the segment
     162             :     // i.e., segment is too long decide the best option
     163             : 
     164           0 :     if ((x0_dif > XLIM) && (x1_dif > XLIM))
     165             :     {
     166             :         // this should not happen
     167             :         // we give-up and return the original tie-point set
     168             : 
     169           0 :         CPLError(
     170             :             CE_Warning, CPLE_AppDefined,
     171             :             "GCPs' set is too large"
     172             :             " to perform the unwrapping! The unwrapping is not performed!");
     173             : 
     174           0 :         return;
     175             :     }
     176           0 :     else if (x1_dif < x0_dif)
     177             :     {
     178             :         // flipped GCPs' set has smaller extent -> unwrapping is performed
     179           0 :         for (int i = 1; i < cnt; ++i)
     180             :         {
     181           0 :             double x0 = gcp[i].dfGCPX;
     182             : 
     183           0 :             gcp[i].dfGCPX = x0 - (x0 > XCNT) * XDIF;
     184             :         }
     185             :     }
     186             : }