OFT

/*******************************************
Copyright (C) 2014
Food and Agriculture Orgazization of the United Nations
and the following contributing authors:

Anssi Pekkarinen
Reija Haapanen

This file is part of Open Foris Geospatial Toolkit which is free software.
You can redistribute it and/or modify it under the terms of the
GNU General Public License as published by the Free Software Foundation,
either version 3 of the License, or (at your option) any later version.

Open Foris Geospatial Toolkit is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with the Open Foris Geospatial Toolkit
If not, see <http://www.gnu.org/licenses/>.

***********************************************************************/
#include "gdal.h"
#include "ogr_srs_api.h"
#include "cpl_string.h"
#include "cpl_conv.h"
#include <time.h>
#include <stdio.h>
#define STAMP "Time-stamp: <2013-12-12 12:39:34 pekkarinen>"
#define NODATA 0
#define VERSION "1.02"

/* fixed crash with images with less than 10 lines */

int comp();
unsigned int *LastPxl = NULL;
int *DBTab = NULL;
int *segP = NULL;
double   **MeanStdTab  = NULL, **MinMaxTab=NULL;
int **SpatTab  = NULL;
FILE  *DataBase = NULL, *ASCDataBase = NULL;
char *tmpfilename = NULL;
int *SpatNeighs;
int Size;
double *SpecNeighs;
/****************************
RH Jan 11 2013
Added usage of options -i, -o, -um and -h
The program can still be used in the old way, based on order of input files
RH Jan 14-18 2013
Added possibility to use the program without the mask, in which case it computes
histogram for the whole image
Also: simplified the double treatment of mask values
RH Jan 18 2013
Commented off the double reporting of image projection
AP Jun 5 2013 TODO FIX -mm -nostd -noavg
AP Dec 12 2013 changed usage

********************************/

/************************************************************/


/************************************************************/



static int
GDALInfoReportCorner( GDALDatasetH hDataset,
                      const char * corner_name,
                      double x, double y );



void AllocError(char *string){

   printf( "Error allocating memory for %s\n",string);
   exit(1);
}

void Usage()

{
  printf( "Version %s %s\n",VERSION,STAMP);
  printf( "Usage:oft-stat -i <input image> -o <output textfile> [-um maskfile] [-mm] [-noavg] [-nostd] [-h help]\n");

    exit(0);
}


/************************************************************************/
/*                                main()                                */
/************************************************************************/




int main( int argc, char ** argv )

{
    GDALDatasetH    hDataset, mDataset;
    GDALRasterBandH hBand, mBand;
    int         i=0, iBand=0 ;
    double      adfGeoTransform[6];
    double             *pafScanline;
    int           *pafMaskline;
    GDALDriverH     hDriver;
    int                 bOutputFile=FALSE, bUseMask=FALSE, bInputImage=FALSE;
    char               *pszMskFilename,*pszFilename = NULL,*ASCfilename = NULL;
    int                nbrBands, nXSize, nYSize, UseStd;

    int row,col ;


    int minSegId, maxSegId;
    double pixVal=0;
    GDALAllRegister();
    int obs;
    double SumOfSquares, SquaredSum;
    int NbrObs;
    int MaskVal;
    int bCalcStd=TRUE,bCalcAvg=TRUE,bCalcMinMax=FALSE;

/* -------------------------------------------------------------------- */
/*      process cmd line parameters                                     */
/*                                                                      */
/*       -i input image                                                 */
/*       -o output  use mask file                                       */
/*       -um use mask file                                              */
/*       -nostd do not compute stds                                     */
/*       -noavg do not compute averages                                 */
/*       -mm compute min and max                                        */
/* -------------------------------------------------------------------- */

 argc = GDALGeneralCmdLineProcessor( argc, &argv, 0 );

    if( argc < 2 ){
    Usage();
        exit( -argc );
    }

    srand((unsigned)time(NULL));

 for( i = 1; i < argc; i++ )
        {
          if( EQUAL(argv[i], "-um") ){
            i++ ;
            if( i >= argc)
              Usage();
            else{
              bUseMask = TRUE;
              pszMskFilename = argv[i];
            }
          }
        else if( EQUAL(argv[i], "-i") ){
            i++ ;
            if( i >= argc)
              Usage();
            else{
              bInputImage = TRUE;
              pszFilename = argv[i];
            }
        }
          else if( EQUAL(argv[i], "-o") ){
            i++ ;
            if( i >= argc)
              Usage();
            else{
              bOutputFile = TRUE;
              ASCfilename = argv[i];
            }
          }
        else if ( EQUAL(argv[i], "-h")){
       printf( "\nUsage:oft-stat -i <infile> -o <outfile> [-um maskfile] [-mm] [-noavg] [-nostd] [-h help]\n");
       printf(" ============================================================\n");
       printf(" Computes image statistics\n");
       printf(" ============================================================\n");
       printf(" Program outputs a text file. The output format in the text file is: ID #pixels avgband1 ...avgbandN stdband1 ...stdbandN\n");
       printf(" You need to give at least the input image file (-i option)\n");
       printf(" and the output text file (-o). Normally, you give also a maskfile (-um)\n");
       printf(" which determines segments for which the statistics are computed\n\n");
       printf("Other parameters, optional:\n");
       printf(" -noavg = program does not compute the averages\n");
       printf(" -nostd = program does not compute the std's\n");
       printf(" -mm = program computes minimum and maximum\n");
       printf(" -h = prints this help\n\n");
       printf("NOTE: for benefit of users running scripts using the older version based on order of\n");
       printf("datafiles instead of options -i, -o and -um, the program can still be used that way\n");

         exit(0);
     }
}//for
//We have enough command line parameters, but not the minimum required options -i and -o
//meaning we use the old mode of this program
//NOTE: old mode always needs the mask... we preserve the old mode just to keep the old scripts running
     if( bInputImage == FALSE || bOutputFile == FALSE ) {

      if(argc > 3) {
      printf("NOTE: Assuming old mode of use without -i, -o (and -um) options\n");
      bUseMask = TRUE;
      pszMskFilename = argv[argc - 3];
      pszFilename = argv[argc - 2];
      ASCfilename = argv[argc -1];

     }
     else {
      printf("Not enough command line parameters, exiting\n");
      Usage();
      exit(0);
        }
        }

//Common for old mode and new mode
      for( i = 1; i < argc; i++ )
    {

      if( EQUAL(argv[i], "-nostd") ){
        bCalcStd=FALSE;
      }
      else if( EQUAL(argv[i], "-noavg") ){
        bCalcAvg=FALSE;
      }
      else if( EQUAL(argv[i], "-mm") ){
        bCalcMinMax=TRUE;
      }

    }




/* -------------------------------------------------------------------- */
/*      Open IO datasets and collect METADATA                           */
/* -------------------------------------------------------------------- */

    hDataset = GDALOpen(pszFilename, GA_ReadOnly );


    if( hDataset == NULL )
      {



        fprintf( stderr,
                 "GDALOpen failed - %d\n%s\n",
                 CPLGetLastErrorNo(), CPLGetLastErrorMsg() );


        CSLDestroy( argv );

        GDALDumpOpenDatasets( stderr );

        GDALDestroyDriverManager();

        CPLDumpSharedList( NULL );


        exit(1);
    }
    else{



      if( GDALGetGeoTransform( hDataset, adfGeoTransform ) == CE_None )
    {
      if( adfGeoTransform[2] == 0.0 && adfGeoTransform[4] == 0.0 )
        {

            printf( "Pixel Size = (%.8f,%.8f)\n",
                    adfGeoTransform[1], adfGeoTransform[5] );
        }
        else
      printf( "GeoTransform =\n"

          "  %.16g, %.16g, %.16g\n"
          "  %.16g, %.16g, %.16g\n",
          adfGeoTransform[0],
          adfGeoTransform[1],
          adfGeoTransform[2],
          adfGeoTransform[3],
          adfGeoTransform[4],
          adfGeoTransform[5] );
    }

      /* -------------------------------------------------------------------- */
      /*      Report corners.                                                 */
      /* -------------------------------------------------------------------- */
      printf( "Corner Coordinates:\n" );
    GDALInfoReportCorner( hDataset, "Upper Left",
                          0.0, 0.0 );
    GDALInfoReportCorner( hDataset, "Lower Left",
                          0.0, GDALGetRasterYSize(hDataset));
    GDALInfoReportCorner( hDataset, "Upper Right",
                          GDALGetRasterXSize(hDataset), 0.0 );
    GDALInfoReportCorner( hDataset, "Lower Right",
                          GDALGetRasterXSize(hDataset),
                          GDALGetRasterYSize(hDataset) );
    GDALInfoReportCorner( hDataset, "Center",
                          GDALGetRasterXSize(hDataset)/2.0,
                          GDALGetRasterYSize(hDataset)/2.0 );

    nbrBands = GDALGetRasterCount( hDataset ) ;

    hBand = GDALGetRasterBand( hDataset, 1 );
    hDriver = GDALGetDatasetDriver( hDataset );

    nXSize = GDALGetRasterBandXSize(hBand);
    nYSize = GDALGetRasterBandYSize(hBand);


    printf("Number of input bands is %i\n",nbrBands);


//RH: the below projection part seems to be twice, let's comment this off
//to make the output more readable
/*

    if( GDALGetProjectionRef( hDataset ) != NULL )
      {
        OGRSpatialReferenceH  hSRS;
        char              *pszProjection;


        pszProjection = (char *) GDALGetProjectionRef( hDataset );

        hSRS = OSRNewSpatialReference(NULL);

        if( OSRImportFromWkt( hSRS, &pszProjection ) == CE_None )
      {
            char    *pszPrettyWkt = NULL;

            OSRExportToPrettyWkt( hSRS, &pszPrettyWkt, FALSE );
            printf( "Coordinate System is:\n%s\n", pszPrettyWkt );
            CPLFree( pszPrettyWkt );
      }
        else
      printf( "Coordinate System is `%s'\n",
          GDALGetProjectionRef( hDataset ) );

        OSRDestroySpatialReference( hSRS );
      }
*/



/* -------------------------------------------------------------------- */
/*      Report general info.                                            */
/* -------------------------------------------------------------------- */


    hDriver = GDALGetDatasetDriver( hDataset );
    printf( "Driver: %s/%s\n",
            GDALGetDriverShortName( hDriver ),
            GDALGetDriverLongName( hDriver ) );

    printf( "Size is %d, %d\n",
            GDALGetRasterXSize( hDataset ),
            GDALGetRasterYSize( hDataset ) );

/* -------------------------------------------------------------------- */
/*      Report projection.                                              */
/* -------------------------------------------------------------------- */
    if( GDALGetProjectionRef( hDataset ) != NULL )
    {
        OGRSpatialReferenceH  hSRS;
        char              *pszProjection;

        pszProjection = (char *) GDALGetProjectionRef( hDataset );

        hSRS = OSRNewSpatialReference(NULL);
        if( OSRImportFromWkt( hSRS, &pszProjection ) == CE_None )
        {
            char    *pszPrettyWkt = NULL;

            OSRExportToPrettyWkt( hSRS, &pszPrettyWkt, FALSE );
            printf( "Coordinate System is:\n%s\n", pszPrettyWkt );
            CPLFree( pszPrettyWkt );
        }
        else
            printf( "Coordinate System is `%s'\n",
                    GDALGetProjectionRef( hDataset ) );

        OSRDestroySpatialReference( hSRS );
    }



      /* proceed  by rows and read all bands in a buffer */




    if((pafScanline = (double *) CPLMalloc(nbrBands * sizeof(double) * nXSize)) == NULL){
      printf("Error allocting memory for input data.\Exiting.\n");
      exit(1);
    }

//RH: added for no mask case:
    if(bUseMask == TRUE) {

    if((pafMaskline = (int *) CPLMalloc(sizeof(int) * nXSize)) == NULL){
      printf("Error allocting memory for mask data.\Exiting.\n");
      exit(1);
      }

      /* find max line number for each sgment */




      obs = 0;

      printf("Finding max mask value\n");
      fflush(stdout);

      mDataset = GDALOpen(pszMskFilename, GA_ReadOnly);
      mBand = GDALGetRasterBand( mDataset, 1 );

      if(mBand == NULL){
    printf("Can't open mask file\n");
    Usage();
      }


    for(row = 0 ; row < nYSize ; row++){


      if(GDALRasterIO(mBand, GF_Read, 0, row, nXSize, 1, pafMaskline, nXSize, 1, GDT_UInt32, 0, 0 )!= CE_None){
    printf("Error reading mask file\n");
    exit(1);
      }

      if(nYSize > 10)
    if( row % (nYSize/10) == 0){
      printf(".");
      fflush(stdout);
    }

      for(col = 0 ; col < nXSize ; col++){

    MaskVal = pafMaskline[col];

    if(MaskVal != 0){

      if(obs == 0){
        minSegId = MaskVal ;
        maxSegId = MaskVal;


      }
      else {

        if(MaskVal > maxSegId)
          maxSegId = MaskVal;
        else if (MaskVal < minSegId)
          minSegId = MaskVal;

        }

        obs ++ ;

      }
      }
    }


    printf("\n");
    fflush(stdout);
}//RH: close if busemask TRUE

        if (bUseMask == FALSE) {
         minSegId = 1;
         maxSegId = 1;
        }

    LastPxl = malloc(sizeof(int) * (maxSegId + 1)) ;

      if(LastPxl == NULL)
    AllocError("LastPxl");

      for ( i = 0 ; i < maxSegId ; i++)
    LastPxl[i] = -1 ;
//RH: do we need this obs any longer?
      obs = 0;

      GDALGetGeoTransform( hDataset, adfGeoTransform);

    }

      for(row = 0 ; row < nYSize ; row++){

       if(bUseMask == TRUE) {

    if(GDALRasterIO(mBand, GF_Read, 0, row, nXSize, 1, pafMaskline, nXSize, 1, GDT_UInt32, 0, 0 )!= CE_None){
      printf("Error reading mask file\n");
      exit(1);
    }
    }
    for(col = 0 ; col < nXSize ; col++){

        if(bUseMask == TRUE) MaskVal = pafMaskline[col];
    else MaskVal = 1;
      if(MaskVal != 0){
        LastPxl[MaskVal] = row * nXSize + col;
        obs ++ ;
      }
    }
      }


      printf("Segment Min val = %i\nSegment Max val = %i\n",minSegId,maxSegId);

      /* alloc memory for tables. We need

         segp = pointer table for segments
     LastPxl = table recoring the last pixel of a segment

      */




      LastPxl = realloc(LastPxl, (maxSegId + 1) * sizeof(int));
      segP    = (int *) malloc((maxSegId + 1) * sizeof(int));

      DBTab   = (int *) malloc(sizeof(int) * (maxSegId + 1));

      if(DBTab ==  NULL)
    AllocError("database pointers");

      if(segP == NULL || LastPxl == NULL)
    AllocError("reallocating tables\n");

      for( i = 0 ; i < maxSegId + 1 ; i++){
    segP[i] =  -1 ;
    DBTab[i] = -1 ;
      }






      if( LastPxl == NULL || segP == NULL || DBTab == NULL)
    AllocError("for segment pointers");

      /* one row may have no more than nXSize segs and an image may not have */




      if(bCalcStd || bCalcAvg){

    MeanStdTab   = (double **)    malloc(sizeof(double   * ) * nXSize);

    if(MeanStdTab == NULL)
      AllocError("for MeanStd information");
    else{

      if(bCalcStd)
        UseStd=2;
      else
        UseStd=1;

          for(i = 0 ; i < nXSize ; i++){
        MeanStdTab[i] = (double *) malloc( UseStd * nbrBands * sizeof(double))  ;
        if(MeanStdTab[i] == NULL)
          AllocError("for statistcs");
          }
    }



      }

      if(bCalcMinMax){
    MinMaxTab  = (double **)    malloc(sizeof(double   * ) * nXSize);
    printf("Allocated memory for storing min and max values\n");
    if(MinMaxTab == NULL)
      AllocError("for MinMax information");
    else
      for(i = 0 ; i < nXSize ; i++){
        MinMaxTab[i] = (double *) malloc(2 * nbrBands * sizeof(double))  ;
        if(MinMaxTab[i] == NULL)
          AllocError("for MeanStd items");
      }
      }



      SpatTab   = (int **)  malloc(sizeof(int * ) * nXSize);

      if(SpatTab == NULL)
    AllocError("for SpatTab information");
      else {

    printf("Allocating memory for Spatial Table\n");

    for(i = 0 ; i < nXSize ; i++){
      SpatTab[i] = (int *) malloc(sizeof(int))  ;
      if(SpatTab[i] == NULL)
        AllocError("SpatTab");
      else
        SpatTab[i][0]=-1;
    }
    printf("Please wait ...\n");
      }



      ASCDataBase = fopen(ASCfilename,"w");

      if(ASCDataBase == NULL){
    printf("Error opening DataBase files. Exiting.\n");
    exit(1);
      }


      for ( col = 0 ; col < nbrBands * nXSize ; col ++)
    pafScanline[col] = 0;

      for(row = 0 ; row < nYSize ; row++){

    if(bUseMask == TRUE) {
      if(GDALRasterIO(mBand, GF_Read, 0, row, nXSize, 1, pafMaskline, nXSize, 1, GDT_UInt32, 0, 0 )!= CE_None){
        printf("Error reading mask file\n");
        exit(1);
      }
    }

    for( iBand = 0; iBand < nbrBands ; iBand++ )
      {


        hBand = GDALGetRasterBand( hDataset, iBand + 1 );


        if(GDALRasterIO( hBand, GF_Read, 0, row, nXSize, 1, &pafScanline[iBand * nXSize], nXSize, 1, GDT_Float64, 0, 0 )!= CE_None)
          {
        printf("Error reading input data in row %i\n",row+1);
        exit(1);
          }

      }

    for(col = 0 ; col < nXSize ; col++){

      if(bUseMask == TRUE)
        MaskVal = pafMaskline[col] ;
      else MaskVal = 1;

      if(MaskVal != 0){


        if(segP[MaskVal] == -1){

          i = 0;

          while(SpatTab[i][0] != -1){

        i++;
          }


          segP[MaskVal] = i;

          SpatTab[segP[MaskVal]][0] = 0;


          if(bCalcAvg)
        for( i = 0 ; i <  nbrBands  ; i++)
          MeanStdTab[segP[MaskVal]][i] = 0;

          if(bCalcStd)
        for( i = 0 ; i <  2 * nbrBands  ; i++)
          MeanStdTab[segP[MaskVal]][i] = 0;
        }



        i = 0;


        SpatTab[segP[MaskVal]][0] += 1;



        for( iBand = 0; iBand < nbrBands ; iBand++ ){
          pixVal = pafScanline[iBand * nXSize + col];

          if(bCalcStd || bCalcAvg){

        /* we need to compute the sum and
           record number of observations  */


        MeanStdTab[segP[MaskVal]][iBand] += pixVal;

        if(bCalcStd)
          MeanStdTab[segP[MaskVal]][nbrBands + iBand] += pixVal * pixVal;

          }


          if(bCalcMinMax){

        if(SpatTab[segP[MaskVal]][0] == 1){

          MinMaxTab[segP[MaskVal]][iBand] = pixVal;
          MinMaxTab[segP[MaskVal]][nbrBands + iBand] = pixVal;
        } else if(pixVal < MinMaxTab[segP[MaskVal]][iBand]){
          MinMaxTab[segP[MaskVal]][iBand] = pixVal;

        } else  if(pixVal > MinMaxTab[segP[MaskVal]][nbrBands + iBand]) {

          MinMaxTab[segP[MaskVal]][nbrBands + iBand] = pixVal ;

        }

          }
        }


        if(SpatTab[segP[MaskVal]][0] != -1 && LastPxl[MaskVal] == row * nXSize + col){



          /*

        Compute the stats and
        print them to the outut file


          */


          if( bCalcStd){
        for(i = 0 ; i < nbrBands ; i++){
          SumOfSquares = MeanStdTab[segP[MaskVal]][nbrBands + i];
          SquaredSum = MeanStdTab[segP[MaskVal]][i] * MeanStdTab[segP[MaskVal]][i];
          NbrObs = SpatTab[segP[MaskVal]][0];
          if(NbrObs > 1)
            MeanStdTab[segP[MaskVal]][nbrBands + i] =  sqrt((SumOfSquares - SquaredSum / NbrObs)/(NbrObs-1));
          else
            MeanStdTab[segP[MaskVal]][nbrBands + i] =  0 ;
        }
          }


          if(bCalcAvg)
        for(i = 0 ; i < nbrBands ; i++)
          MeanStdTab[segP[MaskVal]][i] = MeanStdTab[segP[MaskVal]][i] / SpatTab[segP[MaskVal]][0];



          fprintf(ASCDataBase,"%i %i ",MaskVal,SpatTab[segP[MaskVal]][0]);


          if(bCalcMinMax)
        for(iBand = 0 ; iBand < 2 * nbrBands ; iBand++)
          fprintf(ASCDataBase,"%f ",MinMaxTab[segP[MaskVal]][iBand]);

          if(bCalcAvg)
        for(iBand = 0 ; iBand < nbrBands ; iBand++)
          fprintf(ASCDataBase,"%f ",MeanStdTab[segP[MaskVal]][iBand]);

          if(bCalcStd)
        for(iBand = 0 ; iBand < nbrBands ; iBand++)
          fprintf(ASCDataBase,"%f ",MeanStdTab[segP[MaskVal]][nbrBands + iBand]);

          fprintf(ASCDataBase,"\n");

          SpatTab[segP[MaskVal]][0] = -1;
        }

      }
    }
      }



      GDALClose( hDataset );
      CSLDestroy( argv );
      free(segP);
      free(pafScanline);

      if(bUseMask  == TRUE)
    free(pafMaskline);

      printf("\nDone\n");
      fflush(stdout);
      exit(0);
}


/************************************************************************/
/*                        GDALInfoReportCorner()                        */
/************************************************************************/

static int
GDALInfoReportCorner( GDALDatasetH hDataset,
                      const char * corner_name,
                      double x, double y )

{
    double  dfGeoX, dfGeoY;
    const char  *pszProjection;
    double  adfGeoTransform[6];
    OGRCoordinateTransformationH hTransform = NULL;

    printf( "%-11s ", corner_name );


/* -------------------------------------------------------------------- */
/*      Transform the point into georeferenced coordinates.             */
/* -------------------------------------------------------------------- */
    if( GDALGetGeoTransform( hDataset, adfGeoTransform ) == CE_None )
    {
        pszProjection = GDALGetProjectionRef(hDataset);

        dfGeoX = adfGeoTransform[0] + adfGeoTransform[1] * x
            + adfGeoTransform[2] * y;
        dfGeoY = adfGeoTransform[3] + adfGeoTransform[4] * x
            + adfGeoTransform[5] * y;
    }

    else
    {
        printf( "(%7.1f,%7.1f)\n", x, y );
        return FALSE;
    }

/* -------------------------------------------------------------------- */
/*      Report the georeferenced coordinates.                           */
/* -------------------------------------------------------------------- */
    if( ABS(dfGeoX) < 181 && ABS(dfGeoY) < 91 )
    {
        printf( "(%12.7f,%12.7f) ", dfGeoX, dfGeoY );

    }
    else
    {
        printf( "(%12.3f,%12.3f) ", dfGeoX, dfGeoY );
    }

/* -------------------------------------------------------------------- */
/*      Setup transformation to lat/long.                               */
/* -------------------------------------------------------------------- */
    if( pszProjection != NULL && strlen(pszProjection) > 0 )
    {
        OGRSpatialReferenceH hProj, hLatLong = NULL;

        hProj = OSRNewSpatialReference( pszProjection );
        if( hProj != NULL )
            hLatLong = OSRCloneGeogCS( hProj );

        if( hLatLong != NULL )
        {
            CPLPushErrorHandler( CPLQuietErrorHandler );
            hTransform = OCTNewCoordinateTransformation( hProj, hLatLong );
            CPLPopErrorHandler();

            OSRDestroySpatialReference( hLatLong );
        }

        if( hProj != NULL )
            OSRDestroySpatialReference( hProj );
    }

/* -------------------------------------------------------------------- */
/*      Transform to latlong and report.                                */
/* -------------------------------------------------------------------- */
    if( hTransform != NULL
        && OCTTransform(hTransform,1,&dfGeoX,&dfGeoY,NULL) )
    {

        printf( "(%s,", GDALDecToDMS( dfGeoX, "Long", 2 ) );
        printf( "%s)", GDALDecToDMS( dfGeoY, "Lat", 2 ) );
    }

    if( hTransform != NULL )
        OCTDestroyCoordinateTransformation( hTransform );

    printf( "\n" );

    return TRUE;
}