gtiff/libgeotiff/geo_trans.c

/******************************************************************************
 * $Id: geo_trans.c 8e5eeb35bf76390e3134a4ea7076dab7d478ea0e 2018-11-14 22:55:13 +0100 Even Rouault $
 *
 * Project:  libgeotiff
 * Purpose:  Code to abstract translation between pixel/line and PCS
 *           coordinates.
 * Author:   Frank Warmerdam, warmerda@home.com
 *
 ******************************************************************************
 * Copyright (c) 1999, Frank Warmerdam
 *
 * 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 "geotiff.h"
#include "geo_tiffp.h" /* external TIFF interface */
#include "geo_keyp.h"  /* private interface       */
#include "geokeys.h"

/************************************************************************/
/*                          inv_geotransform()                          */
/*                                                                      */
/*      Invert a 6 term geotransform style matrix.                      */
/************************************************************************/

static int inv_geotransform( double *gt_in, double *gt_out )

{
    double	det, inv_det;

    /* we assume a 3rd row that is [0 0 1] */

    /* Compute determinate */

    det = gt_in[0] * gt_in[4] - gt_in[1] * gt_in[3];

    if( fabs(det) < 0.000000000000001 )
        return 0;

    inv_det = 1.0 / det;

    /* compute adjoint, and divide by determinate */

    gt_out[0] =  gt_in[4] * inv_det;
    gt_out[3] = -gt_in[3] * inv_det;

    gt_out[1] = -gt_in[1] * inv_det;
    gt_out[4] =  gt_in[0] * inv_det;

    gt_out[2] = ( gt_in[1] * gt_in[5] - gt_in[2] * gt_in[4]) * inv_det;
    gt_out[5] = (-gt_in[0] * gt_in[5] + gt_in[2] * gt_in[3]) * inv_det;

    return 1;
}

/************************************************************************/
/*                       GTIFTiepointTranslate()                        */
/************************************************************************/

static
int GTIFTiepointTranslate( int gcp_count, double * gcps_in, double * gcps_out,
                           double x_in, double y_in,
                           double *x_out, double *y_out )

{
    (void) gcp_count;
    (void) gcps_in;
    (void) gcps_out;
    (void) x_in;
    (void) y_in;
    (void) x_out;
    (void) y_out;

    /* I would appreciate a _brief_ block of code for doing second order
       polynomial regression here! */
    return FALSE;
}


/************************************************************************/
/*                           GTIFImageToPCS()                           */
/************************************************************************/

/**
 * Translate a pixel/line coordinate to projection coordinates.
 *
 * At this time this function does not support image to PCS translations for
 * tiepoints-only definitions,  only pixelscale and transformation matrix
 * formulations.
 *
 * @param gtif The handle from GTIFNew() indicating the target file.
 * @param x A pointer to the double containing the pixel offset on input,
 * and into which the easting/longitude will be put on completion.
 * @param y A pointer to the double containing the line offset on input,
 * and into which the northing/latitude will be put on completion.
 *
 * @return TRUE if the transformation succeeds, or FALSE if it fails.  It may
 * fail if the file doesn't have properly setup transformation information,
 * or it is in a form unsupported by this function.
 */

int GTIFImageToPCS( GTIF *gtif, double *x, double *y )

{
    int     res = FALSE;
    int     tiepoint_count, count, transform_count;
    tiff_t *tif=gtif->gt_tif;
    double *tiepoints   = 0;
    double *pixel_scale = 0;
    double *transform   = 0;


    if (!(gtif->gt_methods.get)(tif, GTIFF_TIEPOINTS,
                              &tiepoint_count, &tiepoints ))
        tiepoint_count = 0;

    if (!(gtif->gt_methods.get)(tif, GTIFF_PIXELSCALE, &count, &pixel_scale ))
        count = 0;

    if (!(gtif->gt_methods.get)(tif, GTIFF_TRANSMATRIX,
                                &transform_count, &transform ))
        transform_count = 0;

/* -------------------------------------------------------------------- */
/*      If the pixelscale count is zero, but we have tiepoints use      */
/*      the tiepoint based approach.                                    */
/* -------------------------------------------------------------------- */
    if( tiepoint_count > 6 && count == 0 )
    {
        res = GTIFTiepointTranslate( tiepoint_count / 6,
                                     tiepoints, tiepoints + 3,
                                     *x, *y, x, y );
    }

/* -------------------------------------------------------------------- */
/*	If we have a transformation matrix, use it. 			*/
/* -------------------------------------------------------------------- */
    else if( transform_count == 16 )
    {
        double x_in = *x, y_in = *y;

        *x = x_in * transform[0] + y_in * transform[1] + transform[3];
        *y = x_in * transform[4] + y_in * transform[5] + transform[7];

        res = TRUE;
    }

/* -------------------------------------------------------------------- */
/*      For now we require one tie point, and a valid pixel scale.      */
/* -------------------------------------------------------------------- */
    else if( count < 3 || tiepoint_count < 6 )
    {
        res = FALSE;
    }

    else
    {
        *x = (*x - tiepoints[0]) * pixel_scale[0] + tiepoints[3];
        *y = (*y - tiepoints[1]) * (-1 * pixel_scale[1]) + tiepoints[4];

        res = TRUE;
    }

/* -------------------------------------------------------------------- */
/*      Cleanup                                                         */
/* -------------------------------------------------------------------- */
    if(tiepoints)
        _GTIFFree(tiepoints);
    if(pixel_scale)
        _GTIFFree(pixel_scale);
    if(transform)
        _GTIFFree(transform);

    return res;
}

/************************************************************************/
/*                           GTIFPCSToImage()                           */
/************************************************************************/

/**
 * Translate a projection coordinate to pixel/line coordinates.
 *
 * At this time this function does not support PCS to image translations for
 * tiepoints-only based definitions, only matrix and pixelscale/tiepoints
 * formulations are supposed.
 *
 * @param gtif The handle from GTIFNew() indicating the target file.
 * @param x A pointer to the double containing the pixel offset on input,
 * and into which the easting/longitude will be put on completion.
 * @param y A pointer to the double containing the line offset on input,
 * and into which the northing/latitude will be put on completion.
 *
 * @return TRUE if the transformation succeeds, or FALSE if it fails.  It may
 * fail if the file doesn't have properly setup transformation information,
 * or it is in a form unsupported by this function.
 */

int GTIFPCSToImage( GTIF *gtif, double *x, double *y )

{
    double 	*tiepoints = NULL;
    int 	tiepoint_count, count, transform_count = 0;
    double	*pixel_scale = NULL;
    double 	*transform   = NULL;
    tiff_t 	*tif=gtif->gt_tif;
    int		result = FALSE;

/* -------------------------------------------------------------------- */
/*      Fetch tiepoints and pixel scale.                                */
/* -------------------------------------------------------------------- */
    if (!(gtif->gt_methods.get)(tif, GTIFF_TIEPOINTS,
                              &tiepoint_count, &tiepoints ))
        tiepoint_count = 0;

    if (!(gtif->gt_methods.get)(tif, GTIFF_PIXELSCALE, &count, &pixel_scale ))
        count = 0;

    if (!(gtif->gt_methods.get)(tif, GTIFF_TRANSMATRIX,
                                &transform_count, &transform ))
        transform_count = 0;

/* -------------------------------------------------------------------- */
/*      If the pixelscale count is zero, but we have tiepoints use      */
/*      the tiepoint based approach.                                    */
/* -------------------------------------------------------------------- */
    if( tiepoint_count > 6 && count == 0 )
    {
        result = GTIFTiepointTranslate( tiepoint_count / 6,
                                        tiepoints + 3, tiepoints,
                                        *x, *y, x, y );
    }

/* -------------------------------------------------------------------- */
/*      Handle matrix - convert to "geotransform" format, invert and    */
/*      apply.                                                          */
/* -------------------------------------------------------------------- */
    else if( transform_count == 16 )
    {
        double  x_in = *x, y_in = *y;
        double	gt_in[6], gt_out[6];

        gt_in[0] = transform[0];
        gt_in[1] = transform[1];
        gt_in[2] = transform[3];
        gt_in[3] = transform[4];
        gt_in[4] = transform[5];
        gt_in[5] = transform[7];

        if( !inv_geotransform( gt_in, gt_out ) )
            result = FALSE;
        else
        {
            *x = x_in * gt_out[0] + y_in * gt_out[1] + gt_out[2];
            *y = x_in * gt_out[3] + y_in * gt_out[4] + gt_out[5];

            result = TRUE;
        }
    }

/* -------------------------------------------------------------------- */
/*      For now we require one tie point, and a valid pixel scale.      */
/* -------------------------------------------------------------------- */
    else if( count >= 3 && tiepoint_count >= 6 )
    {
        *x = (*x - tiepoints[3]) / pixel_scale[0] + tiepoints[0];
        *y = (*y - tiepoints[4]) / (-1 * pixel_scale[1]) + tiepoints[1];

        result = TRUE;
    }

/* -------------------------------------------------------------------- */
/*      Cleanup.                                                        */
/* -------------------------------------------------------------------- */
    if(tiepoints)
        _GTIFFree(tiepoints);
    if(pixel_scale)
        _GTIFFree(pixel_scale);
    if(transform)
        _GTIFFree(transform);

    return result;
}