EVOLUTION-MANAGER

Edit File: unwrapgcps.cpp

/******************************************************************************
 * $Id: unwrapgcps.cpp 27098 2014-03-27 00:16:11Z rouault $
 *
 * Project:  APP ENVISAT Support
 * Purpose:  GCPs Unwrapping for products crossing the WGS84 date-line  
 * Author:   Martin Paces martin.paces@eox.at 
 *
 ******************************************************************************
 * Copyright (c) 2013, EOX IT Services, GmbH 
 *
 * Permission is hereby granted, free of charge, to any person obtaining a
 * copy of this software and associated documentation files (the "Software"),
 * to deal in the Software without restriction, including without limitation
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
 * and/or sell copies of the Software, and to permit persons to whom the
 * Software is furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included
 * in all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 * DEALINGS IN THE SOFTWARE.
 ****************************************************************************/

#include "gdal.h"
#include <cmath>
#include <cstdio>

// number of histogram bins (36 a 10dg)
#define NBIN  36  
// number of empty bins to guess the flip-point
#define NEMPY 7

// WGS84 bounds 
#define XMIN  (-180.0)
#define XMAX  (+180.0)
#define XDIF  (+360.0)
#define XCNT     (0.0)

// max. allowed longitude extent of the GCP set 
#define XLIM  (XDIF*(1.0-NEMPY*(1.0/NBIN)))

// The algoright is based on assumption that the unwrapped 
// GCPs ('flipped' values) have smaller extent along the longitude.
// We further assume that the lenght of the striplines is limitted 
// to one orbit and does not exceeded given limit along the longitude, 
// e.i., the wrapped-arround coordinates have significantly larger 
// extent the unwrapped. If the smaller extend exceedes the limit
// the original tiepoints are returned.

static double _suggest_flip_point( const int cnt, GDAL_GCP *gcp ) 
{ 
    // the histogram array - it is expected to fit the stack
    int hist[NBIN] ; 
    
    // reset the histogram counters
    for( int i = 0 ; i < NBIN ; i++ ) hist[i] = 0 ;

// accumulate the histogram 
    for( int i = 0 ; i < cnt ; i++ ) 
    { 
        double x = (gcp[i].dfGCPX-XMIN)/XDIF ;
        int idx = (int)(NBIN*(x-floor(x))) ;

// the lattitudes should lay in the +/-180 bounds 
        // although it should never happen we check the outliers
        if (idx < 0) idx = 0 ; 
        if (idx >= NBIN ) idx = NBIN-1;

hist[idx] += 1 ; 
    }

// find a middle of at least NEMPTY consecutive empty bins and get its middle
    int i0 = -1 , i1 = -1 , last_is_empty = 0 ; 
    for( int i = 0 ; i < (2*NBIN-1) ; i++ ) 
    { 
        if ( 0 == hist[i%NBIN] ) // empty 
        { 
            if ( !last_is_empty ) // re-start counter 
            { 
                i0 = i ;
                last_is_empty = 1 ; 
            }
        } 
        else // non-empty 
        {
            if ( last_is_empty ) 
            { 
                i1 = i ;
                last_is_empty = 0 ;

// if the segment is long enough -> terminate 
                if (( i1 - i0 )>=NEMPY) break ; 
            } 
        }
    }

// if all full or all empty the returning default value 
    if ( i1 < 0 ) return XCNT ;

// return the flip-centre

double tmp = ((i1-i0)*0.5+i0)/((float)NBIN) ;

return (tmp-floor(tmp))*XDIF + XMIN ;
}

void EnvisatUnwrapGCPs( int cnt, GDAL_GCP *gcp ) 
{ 
    if ( cnt < 1 ) return ; 
    
    // suggest right flip-point 
    double x_flip = _suggest_flip_point( cnt, gcp );

// find the limits allong the longitude (x) for flipped and unflipped values
    
    int cnt_flip = 0 ; // flipped values' counter
    double x0_dif , x1_dif ;

{ 
        double x0_min, x0_max, x1_min, x1_max ;

{ 
            double x0 = gcp[0].dfGCPX ;
            int  flip = (x0>x_flip) ;
            x0_min = x0_max = x0 ; 
            x1_min = x1_max = x0 - flip*XDIF ; 
            cnt_flip += flip ; // count the flipped values 
        }

for ( int i = 1 ; i < cnt ; ++i ) 
        { 
            double x0 = gcp[i].dfGCPX ;
            int  flip = (x0>x_flip) ;
            double x1 = x0 - flip*XDIF ; // flipped value
            cnt_flip += flip ; // count the flipped values

if ( x0 > x0_max ) x0_max = x0 ; 
            if ( x0 < x0_min ) x0_min = x0 ; 
            if ( x1 > x1_max ) x1_max = x1 ; 
            if ( x1 < x1_min ) x1_min = x1 ; 
        }

x0_dif = x0_max - x0_min ; 
        x1_dif = x1_max - x1_min ; 
    }

// in case all values either flipped or non-flipped 
    // nothing is to be done 
    if (( cnt_flip == 0 ) || ( cnt_flip == cnt )) return ;

// check whether we need to split the segment 
    // i.e., segment is too long decide the best option 
    
    if (( x0_dif > XLIM ) && ( x1_dif > XLIM )) 
    { 
        // this should not happen 
        // we give-up and return the original tie-point set

CPLError(CE_Warning,CPLE_AppDefined,"GCPs' set is too large"
            " to perform the unwrapping! The unwrapping is not performed!");

return ; 
    } 
    else if ( x1_dif < x0_dif ) 
    { 
        // flipped GCPs' set has smaller extent -> unwrapping is performed
        for ( int i = 1 ; i < cnt ; ++i ) 
        { 
            double x0 = gcp[i].dfGCPX ;

gcp[i].dfGCPX = x0 - (x0>XCNT)*XDIF ; 
        } 
    } 
}