SHARE
TWEET

OFT

a guest Jun 21st, 2017 2 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. /*******************************************
  2. Copyright (C) 2014
  3. Food and Agriculture Orgazization of the United Nations
  4. and the following contributing authors:
  5.  
  6. Anssi Pekkarinen
  7. Reija Haapanen
  8.  
  9. This file is part of Open Foris Geospatial Toolkit which is free software.
  10. You can redistribute it and/or modify it under the terms of the
  11. GNU General Public License as published by the Free Software Foundation,
  12. either version 3 of the License, or (at your option) any later version.
  13.  
  14. Open Foris Geospatial Toolkit is distributed in the hope that it will be useful,
  15. but WITHOUT ANY WARRANTY; without even the implied warranty of
  16. MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  17. GNU General Public License for more details.
  18.  
  19. You should have received a copy of the GNU General Public License
  20. along with the Open Foris Geospatial Toolkit  
  21. If not, see <http://www.gnu.org/licenses/>.
  22.  
  23. ***********************************************************************/
  24. #include "gdal.h"
  25. #include "ogr_srs_api.h"
  26. #include "cpl_string.h"
  27. #include "cpl_conv.h"
  28. #include <time.h>
  29. #include <stdio.h>
  30. #define STAMP "Time-stamp: <2013-12-12 12:39:34 pekkarinen>"
  31. #define NODATA 0
  32. #define VERSION "1.02"
  33.  
  34. /* fixed crash with images with less than 10 lines */
  35.  
  36. int comp();
  37. unsigned int *LastPxl = NULL;
  38. int *DBTab = NULL;
  39. int *segP = NULL;
  40. double   **MeanStdTab  = NULL, **MinMaxTab=NULL;
  41. int **SpatTab  = NULL;
  42. FILE  *DataBase = NULL, *ASCDataBase = NULL;    
  43. char *tmpfilename = NULL;
  44. int *SpatNeighs;
  45. int Size;
  46. double *SpecNeighs;
  47. /****************************
  48. RH Jan 11 2013
  49. Added usage of options -i, -o, -um and -h
  50. The program can still be used in the old way, based on order of input files
  51. RH Jan 14-18 2013
  52. Added possibility to use the program without the mask, in which case it computes
  53. histogram for the whole image
  54. Also: simplified the double treatment of mask values
  55. RH Jan 18 2013
  56. Commented off the double reporting of image projection
  57. AP Jun 5 2013 TODO FIX -mm -nostd -noavg
  58. AP Dec 12 2013 changed usage
  59.  
  60. ********************************/
  61.  
  62. /************************************************************/
  63.    
  64.  
  65. /************************************************************/
  66.  
  67.  
  68.  
  69. static int
  70. GDALInfoReportCorner( GDALDatasetH hDataset,
  71.                       const char * corner_name,
  72.                       double x, double y );
  73.  
  74.  
  75.  
  76. void AllocError(char *string){
  77.  
  78.    printf( "Error allocating memory for %s\n",string);
  79.    exit(1);
  80. }
  81.  
  82. void Usage()
  83.  
  84. {
  85.   printf( "Version %s %s\n",VERSION,STAMP);
  86.   printf( "Usage:oft-stat -i <input image> -o <output textfile> [-um maskfile] [-mm] [-noavg] [-nostd] [-h help]\n");
  87.  
  88.     exit(0);
  89. }
  90.  
  91.  
  92. /************************************************************************/
  93. /*                                main()                                */
  94. /************************************************************************/
  95.  
  96.  
  97.  
  98.  
  99. int main( int argc, char ** argv )
  100.  
  101. {
  102.     GDALDatasetH    hDataset, mDataset;
  103.     GDALRasterBandH hBand, mBand;
  104.     int         i=0, iBand=0 ;
  105.     double      adfGeoTransform[6];
  106.     double             *pafScanline;
  107.     int           *pafMaskline;
  108.     GDALDriverH     hDriver;
  109.     int                 bOutputFile=FALSE, bUseMask=FALSE, bInputImage=FALSE;
  110.     char               *pszMskFilename,*pszFilename = NULL,*ASCfilename = NULL;
  111.     int                nbrBands, nXSize, nYSize, UseStd;
  112.  
  113.     int row,col ;
  114.  
  115.      
  116.     int minSegId, maxSegId;
  117.     double pixVal=0;
  118.     GDALAllRegister();
  119.     int obs;
  120.     double SumOfSquares, SquaredSum;
  121.     int NbrObs;
  122.     int MaskVal;
  123.     int bCalcStd=TRUE,bCalcAvg=TRUE,bCalcMinMax=FALSE;
  124.  
  125. /* -------------------------------------------------------------------- */
  126. /*      process cmd line parameters                                     */
  127. /*                                                                      */
  128. /*       -i input image                                                 */
  129. /*       -o output  use mask file                                       */
  130. /*       -um use mask file                                              */
  131. /*       -nostd do not compute stds                                     */
  132. /*       -noavg do not compute averages                                 */
  133. /*       -mm compute min and max                                        */
  134. /* -------------------------------------------------------------------- */
  135.    
  136.  argc = GDALGeneralCmdLineProcessor( argc, &argv, 0 );
  137.  
  138.     if( argc < 2 ){
  139.     Usage();
  140.         exit( -argc );
  141.     }
  142.  
  143.     srand((unsigned)time(NULL));
  144.  
  145.  for( i = 1; i < argc; i++ )
  146.         {
  147.           if( EQUAL(argv[i], "-um") ){
  148.             i++ ;
  149.             if( i >= argc)
  150.               Usage();
  151.             else{
  152.               bUseMask = TRUE;
  153.               pszMskFilename = argv[i];
  154.             }
  155.           }
  156.         else if( EQUAL(argv[i], "-i") ){
  157.             i++ ;
  158.             if( i >= argc)
  159.               Usage();
  160.             else{
  161.               bInputImage = TRUE;
  162.               pszFilename = argv[i];
  163.             }
  164.         }
  165.           else if( EQUAL(argv[i], "-o") ){
  166.             i++ ;
  167.             if( i >= argc)
  168.               Usage();
  169.             else{
  170.               bOutputFile = TRUE;
  171.               ASCfilename = argv[i];
  172.             }
  173.           }
  174.         else if ( EQUAL(argv[i], "-h")){
  175.        printf( "\nUsage:oft-stat -i <infile> -o <outfile> [-um maskfile] [-mm] [-noavg] [-nostd] [-h help]\n");
  176.        printf(" ============================================================\n");
  177.        printf(" Computes image statistics\n");
  178.        printf(" ============================================================\n");
  179.        printf(" Program outputs a text file. The output format in the text file is: ID #pixels avgband1 ...avgbandN stdband1 ...stdbandN\n");
  180.        printf(" You need to give at least the input image file (-i option)\n");
  181.        printf(" and the output text file (-o). Normally, you give also a maskfile (-um)\n");
  182.        printf(" which determines segments for which the statistics are computed\n\n");
  183.        printf("Other parameters, optional:\n");
  184.        printf(" -noavg = program does not compute the averages\n");
  185.        printf(" -nostd = program does not compute the std's\n");
  186.        printf(" -mm = program computes minimum and maximum\n");
  187.        printf(" -h = prints this help\n\n");
  188.        printf("NOTE: for benefit of users running scripts using the older version based on order of\n");
  189.        printf("datafiles instead of options -i, -o and -um, the program can still be used that way\n");
  190.  
  191.          exit(0);
  192.      }
  193. }//for
  194. //We have enough command line parameters, but not the minimum required options -i and -o
  195. //meaning we use the old mode of this program
  196. //NOTE: old mode always needs the mask... we preserve the old mode just to keep the old scripts running
  197.      if( bInputImage == FALSE || bOutputFile == FALSE ) {
  198.  
  199.       if(argc > 3) {
  200.       printf("NOTE: Assuming old mode of use without -i, -o (and -um) options\n");
  201.       bUseMask = TRUE;
  202.       pszMskFilename = argv[argc - 3];
  203.       pszFilename = argv[argc - 2];
  204.       ASCfilename = argv[argc -1];
  205.  
  206.      }
  207.      else {
  208.       printf("Not enough command line parameters, exiting\n");
  209.       Usage();
  210.       exit(0);
  211.         }
  212.         }
  213.  
  214. //Common for old mode and new mode
  215.       for( i = 1; i < argc; i++ )
  216.     {
  217.      
  218.       if( EQUAL(argv[i], "-nostd") ){
  219.         bCalcStd=FALSE;
  220.       }
  221.       else if( EQUAL(argv[i], "-noavg") ){
  222.         bCalcAvg=FALSE;
  223.       }
  224.       else if( EQUAL(argv[i], "-mm") ){
  225.         bCalcMinMax=TRUE;
  226.       }
  227.  
  228.     }
  229.  
  230.  
  231.    
  232.  
  233. /* -------------------------------------------------------------------- */
  234. /*      Open IO datasets and collect METADATA                           */
  235. /* -------------------------------------------------------------------- */
  236.  
  237.     hDataset = GDALOpen(pszFilename, GA_ReadOnly );
  238.  
  239.  
  240.     if( hDataset == NULL )
  241.       {
  242.  
  243.  
  244.  
  245.         fprintf( stderr,
  246.                  "GDALOpen failed - %d\n%s\n",
  247.                  CPLGetLastErrorNo(), CPLGetLastErrorMsg() );
  248.    
  249.    
  250.         CSLDestroy( argv );
  251.    
  252.         GDALDumpOpenDatasets( stderr );
  253.  
  254.         GDALDestroyDriverManager();
  255.  
  256.         CPLDumpSharedList( NULL );
  257.  
  258.    
  259.         exit(1);
  260.     }
  261.     else{
  262.  
  263.  
  264.  
  265.       if( GDALGetGeoTransform( hDataset, adfGeoTransform ) == CE_None )
  266.     {
  267.       if( adfGeoTransform[2] == 0.0 && adfGeoTransform[4] == 0.0 )
  268.         {
  269.          
  270.             printf( "Pixel Size = (%.8f,%.8f)\n",
  271.                     adfGeoTransform[1], adfGeoTransform[5] );
  272.         }
  273.         else
  274.       printf( "GeoTransform =\n"
  275.  
  276.           "  %.16g, %.16g, %.16g\n"
  277.           "  %.16g, %.16g, %.16g\n",
  278.           adfGeoTransform[0],
  279.           adfGeoTransform[1],
  280.           adfGeoTransform[2],
  281.           adfGeoTransform[3],
  282.           adfGeoTransform[4],
  283.           adfGeoTransform[5] );
  284.     }
  285.      
  286.       /* -------------------------------------------------------------------- */
  287.       /*      Report corners.                                                 */
  288.       /* -------------------------------------------------------------------- */
  289.       printf( "Corner Coordinates:\n" );
  290.     GDALInfoReportCorner( hDataset, "Upper Left",
  291.                           0.0, 0.0 );
  292.     GDALInfoReportCorner( hDataset, "Lower Left",
  293.                           0.0, GDALGetRasterYSize(hDataset));
  294.     GDALInfoReportCorner( hDataset, "Upper Right",
  295.                           GDALGetRasterXSize(hDataset), 0.0 );
  296.     GDALInfoReportCorner( hDataset, "Lower Right",
  297.                           GDALGetRasterXSize(hDataset),
  298.                           GDALGetRasterYSize(hDataset) );
  299.     GDALInfoReportCorner( hDataset, "Center",
  300.                           GDALGetRasterXSize(hDataset)/2.0,
  301.                           GDALGetRasterYSize(hDataset)/2.0 );
  302.  
  303.     nbrBands = GDALGetRasterCount( hDataset ) ;  
  304.  
  305.     hBand = GDALGetRasterBand( hDataset, 1 );
  306.     hDriver = GDALGetDatasetDriver( hDataset );    
  307.  
  308.     nXSize = GDALGetRasterBandXSize(hBand);
  309.     nYSize = GDALGetRasterBandYSize(hBand);
  310.  
  311.  
  312.     printf("Number of input bands is %i\n",nbrBands);
  313.      
  314.  
  315. //RH: the below projection part seems to be twice, let's comment this off
  316. //to make the output more readable
  317. /*
  318.  
  319.     if( GDALGetProjectionRef( hDataset ) != NULL )
  320.       {
  321.         OGRSpatialReferenceH  hSRS;
  322.         char              *pszProjection;
  323.    
  324.  
  325.         pszProjection = (char *) GDALGetProjectionRef( hDataset );
  326.    
  327.         hSRS = OSRNewSpatialReference(NULL);
  328.  
  329.         if( OSRImportFromWkt( hSRS, &pszProjection ) == CE_None )
  330.       {
  331.             char    *pszPrettyWkt = NULL;
  332.        
  333.             OSRExportToPrettyWkt( hSRS, &pszPrettyWkt, FALSE );
  334.             printf( "Coordinate System is:\n%s\n", pszPrettyWkt );
  335.             CPLFree( pszPrettyWkt );
  336.       }
  337.         else
  338.       printf( "Coordinate System is `%s'\n",
  339.           GDALGetProjectionRef( hDataset ) );
  340.    
  341.         OSRDestroySpatialReference( hSRS );
  342.       }
  343. */
  344.  
  345.  
  346.  
  347. /* -------------------------------------------------------------------- */
  348. /*      Report general info.                                            */
  349. /* -------------------------------------------------------------------- */
  350.    
  351.  
  352.     hDriver = GDALGetDatasetDriver( hDataset );
  353.     printf( "Driver: %s/%s\n",
  354.             GDALGetDriverShortName( hDriver ),
  355.             GDALGetDriverLongName( hDriver ) );
  356.  
  357.     printf( "Size is %d, %d\n",
  358.             GDALGetRasterXSize( hDataset ),
  359.             GDALGetRasterYSize( hDataset ) );
  360.    
  361. /* -------------------------------------------------------------------- */
  362. /*      Report projection.                                              */
  363. /* -------------------------------------------------------------------- */
  364.     if( GDALGetProjectionRef( hDataset ) != NULL )
  365.     {
  366.         OGRSpatialReferenceH  hSRS;
  367.         char              *pszProjection;
  368.  
  369.         pszProjection = (char *) GDALGetProjectionRef( hDataset );
  370.  
  371.         hSRS = OSRNewSpatialReference(NULL);
  372.         if( OSRImportFromWkt( hSRS, &pszProjection ) == CE_None )
  373.         {
  374.             char    *pszPrettyWkt = NULL;
  375.  
  376.             OSRExportToPrettyWkt( hSRS, &pszPrettyWkt, FALSE );
  377.             printf( "Coordinate System is:\n%s\n", pszPrettyWkt );
  378.             CPLFree( pszPrettyWkt );
  379.         }
  380.         else
  381.             printf( "Coordinate System is `%s'\n",
  382.                     GDALGetProjectionRef( hDataset ) );
  383.  
  384.         OSRDestroySpatialReference( hSRS );
  385.     }
  386.  
  387.  
  388.  
  389.       /* proceed  by rows and read all bands in a buffer */
  390.  
  391.  
  392.  
  393.  
  394.     if((pafScanline = (double *) CPLMalloc(nbrBands * sizeof(double) * nXSize)) == NULL){
  395.       printf("Error allocting memory for input data.\Exiting.\n");
  396.       exit(1);
  397.     }
  398.  
  399. //RH: added for no mask case:
  400.     if(bUseMask == TRUE) {
  401.  
  402.     if((pafMaskline = (int *) CPLMalloc(sizeof(int) * nXSize)) == NULL){
  403.       printf("Error allocting memory for mask data.\Exiting.\n");
  404.       exit(1);
  405.       }
  406.  
  407.       /* find max line number for each sgment */
  408.  
  409.  
  410.  
  411.  
  412.       obs = 0;
  413.  
  414.       printf("Finding max mask value\n");
  415.       fflush(stdout);
  416.  
  417.       mDataset = GDALOpen(pszMskFilename, GA_ReadOnly);
  418.       mBand = GDALGetRasterBand( mDataset, 1 );      
  419.  
  420.       if(mBand == NULL){
  421.     printf("Can't open mask file\n");
  422.     Usage();
  423.       }
  424.  
  425.  
  426.     for(row = 0 ; row < nYSize ; row++){
  427.  
  428.    
  429.       if(GDALRasterIO(mBand, GF_Read, 0, row, nXSize, 1, pafMaskline, nXSize, 1, GDT_UInt32, 0, 0 )!= CE_None){
  430.     printf("Error reading mask file\n");
  431.     exit(1);
  432.       }
  433.      
  434.       if(nYSize > 10)
  435.     if( row % (nYSize/10) == 0){
  436.       printf(".");
  437.       fflush(stdout);
  438.     }
  439.  
  440.       for(col = 0 ; col < nXSize ; col++){
  441.      
  442.     MaskVal = pafMaskline[col];
  443.  
  444.     if(MaskVal != 0){
  445.  
  446.       if(obs == 0){
  447.         minSegId = MaskVal ;
  448.         maxSegId = MaskVal;
  449.  
  450.  
  451.       }
  452.       else {
  453.  
  454.         if(MaskVal > maxSegId)
  455.           maxSegId = MaskVal;
  456.         else if (MaskVal < minSegId)
  457.           minSegId = MaskVal;
  458.          
  459.         }
  460.        
  461.         obs ++ ;
  462.        
  463.       }
  464.       }
  465.     }
  466.    
  467.  
  468.     printf("\n");
  469.     fflush(stdout);
  470. }//RH: close if busemask TRUE
  471.  
  472.         if (bUseMask == FALSE) {
  473.          minSegId = 1;
  474.          maxSegId = 1;
  475.         }
  476.  
  477.     LastPxl = malloc(sizeof(int) * (maxSegId + 1)) ;
  478.      
  479.       if(LastPxl == NULL)
  480.     AllocError("LastPxl");
  481.  
  482.       for ( i = 0 ; i < maxSegId ; i++)
  483.     LastPxl[i] = -1 ;
  484. //RH: do we need this obs any longer?
  485.       obs = 0;
  486.      
  487.       GDALGetGeoTransform( hDataset, adfGeoTransform);
  488.  
  489.     }
  490.  
  491.       for(row = 0 ; row < nYSize ; row++){
  492.  
  493.        if(bUseMask == TRUE) {
  494.    
  495.     if(GDALRasterIO(mBand, GF_Read, 0, row, nXSize, 1, pafMaskline, nXSize, 1, GDT_UInt32, 0, 0 )!= CE_None){
  496.       printf("Error reading mask file\n");
  497.       exit(1);
  498.     }
  499.     }
  500.     for(col = 0 ; col < nXSize ; col++){
  501.  
  502.         if(bUseMask == TRUE) MaskVal = pafMaskline[col];
  503.     else MaskVal = 1;
  504.       if(MaskVal != 0){
  505.         LastPxl[MaskVal] = row * nXSize + col;
  506.         obs ++ ;
  507.       }
  508.     }
  509.       }
  510.  
  511.  
  512.       printf("Segment Min val = %i\nSegment Max val = %i\n",minSegId,maxSegId);
  513.  
  514.       /* alloc memory for tables. We need
  515.      
  516.          segp = pointer table for segments
  517.      LastPxl = table recoring the last pixel of a segment
  518.      
  519.       */
  520.  
  521.  
  522.  
  523.  
  524.       LastPxl = realloc(LastPxl, (maxSegId + 1) * sizeof(int));
  525.       segP    = (int *) malloc((maxSegId + 1) * sizeof(int));
  526.      
  527.       DBTab   = (int *) malloc(sizeof(int) * (maxSegId + 1));    
  528.  
  529.       if(DBTab ==  NULL)
  530.     AllocError("database pointers");
  531.      
  532.       if(segP == NULL || LastPxl == NULL)
  533.     AllocError("reallocating tables\n");
  534.  
  535.       for( i = 0 ; i < maxSegId + 1 ; i++){
  536.     segP[i] =  -1 ;
  537.     DBTab[i] = -1 ;
  538.       }
  539.  
  540.  
  541.  
  542.  
  543.      
  544.      
  545.       if( LastPxl == NULL || segP == NULL || DBTab == NULL)
  546.     AllocError("for segment pointers");    
  547.      
  548.       /* one row may have no more than nXSize segs and an image may not have */
  549.  
  550.      
  551.  
  552.  
  553.       if(bCalcStd || bCalcAvg){
  554.      
  555.     MeanStdTab   = (double **)    malloc(sizeof(double   * ) * nXSize);
  556.  
  557.     if(MeanStdTab == NULL)
  558.       AllocError("for MeanStd information");    
  559.     else{
  560.  
  561.       if(bCalcStd)
  562.         UseStd=2;
  563.       else
  564.         UseStd=1;
  565.          
  566.           for(i = 0 ; i < nXSize ; i++){
  567.         MeanStdTab[i] = (double *) malloc( UseStd * nbrBands * sizeof(double))  ;      
  568.         if(MeanStdTab[i] == NULL)
  569.           AllocError("for statistcs");
  570.           }
  571.     }
  572.    
  573.  
  574.    
  575.       }
  576.  
  577.       if(bCalcMinMax){
  578.     MinMaxTab  = (double **)    malloc(sizeof(double   * ) * nXSize);
  579.     printf("Allocated memory for storing min and max values\n");
  580.     if(MinMaxTab == NULL)
  581.       AllocError("for MinMax information");
  582.     else
  583.       for(i = 0 ; i < nXSize ; i++){
  584.         MinMaxTab[i] = (double *) malloc(2 * nbrBands * sizeof(double))  ;
  585.         if(MinMaxTab[i] == NULL)
  586.           AllocError("for MeanStd items");
  587.       }
  588.       }
  589.  
  590.  
  591.  
  592.       SpatTab   = (int **)  malloc(sizeof(int * ) * nXSize);
  593.      
  594.       if(SpatTab == NULL)
  595.     AllocError("for SpatTab information");    
  596.       else {
  597.  
  598.     printf("Allocating memory for Spatial Table\n");
  599.    
  600.     for(i = 0 ; i < nXSize ; i++){
  601.       SpatTab[i] = (int *) malloc(sizeof(int))  ;
  602.       if(SpatTab[i] == NULL)
  603.         AllocError("SpatTab");
  604.       else
  605.         SpatTab[i][0]=-1;
  606.     }
  607.     printf("Please wait ...\n");
  608.       }
  609.      
  610.      
  611.      
  612.       ASCDataBase = fopen(ASCfilename,"w");
  613.  
  614.       if(ASCDataBase == NULL){
  615.     printf("Error opening DataBase files. Exiting.\n");
  616.     exit(1);
  617.       }
  618.      
  619.  
  620.       for ( col = 0 ; col < nbrBands * nXSize ; col ++)
  621.     pafScanline[col] = 0;
  622.      
  623.       for(row = 0 ; row < nYSize ; row++){
  624.    
  625.     if(bUseMask == TRUE) {
  626.       if(GDALRasterIO(mBand, GF_Read, 0, row, nXSize, 1, pafMaskline, nXSize, 1, GDT_UInt32, 0, 0 )!= CE_None){
  627.         printf("Error reading mask file\n");
  628.         exit(1);
  629.       }
  630.     }
  631.    
  632.     for( iBand = 0; iBand < nbrBands ; iBand++ )
  633.       {
  634.        
  635.        
  636.         hBand = GDALGetRasterBand( hDataset, iBand + 1 );
  637.        
  638.  
  639.         if(GDALRasterIO( hBand, GF_Read, 0, row, nXSize, 1, &pafScanline[iBand * nXSize], nXSize, 1, GDT_Float64, 0, 0 )!= CE_None)
  640.           {        
  641.         printf("Error reading input data in row %i\n",row+1);
  642.         exit(1);
  643.           }
  644.        
  645.       }
  646.    
  647.     for(col = 0 ; col < nXSize ; col++){
  648.  
  649.       if(bUseMask == TRUE)
  650.         MaskVal = pafMaskline[col] ;
  651.       else MaskVal = 1;
  652.        
  653.       if(MaskVal != 0){
  654.        
  655.        
  656.         if(segP[MaskVal] == -1){
  657.  
  658.           i = 0;
  659.          
  660.           while(SpatTab[i][0] != -1){
  661.        
  662.         i++;
  663.           }
  664.          
  665.  
  666.           segP[MaskVal] = i;
  667.          
  668.           SpatTab[segP[MaskVal]][0] = 0;
  669.          
  670.          
  671.           if(bCalcAvg)
  672.         for( i = 0 ; i <  nbrBands  ; i++)
  673.           MeanStdTab[segP[MaskVal]][i] = 0;
  674.          
  675.           if(bCalcStd)
  676.         for( i = 0 ; i <  2 * nbrBands  ; i++)
  677.           MeanStdTab[segP[MaskVal]][i] = 0;
  678.         }
  679.        
  680.        
  681.        
  682.         i = 0;
  683.  
  684.        
  685.         SpatTab[segP[MaskVal]][0] += 1;    
  686.        
  687.          
  688.          
  689.         for( iBand = 0; iBand < nbrBands ; iBand++ ){        
  690.           pixVal = pafScanline[iBand * nXSize + col];
  691.          
  692.           if(bCalcStd || bCalcAvg){
  693.          
  694.         /* we need to compute the sum and
  695.            record number of observations  */
  696.        
  697.        
  698.         MeanStdTab[segP[MaskVal]][iBand] += pixVal;
  699.        
  700.         if(bCalcStd)
  701.           MeanStdTab[segP[MaskVal]][nbrBands + iBand] += pixVal * pixVal;      
  702.          
  703.           }    
  704.          
  705.          
  706.           if(bCalcMinMax){
  707.          
  708.         if(SpatTab[segP[MaskVal]][0] == 1){
  709.          
  710.           MinMaxTab[segP[MaskVal]][iBand] = pixVal;
  711.           MinMaxTab[segP[MaskVal]][nbrBands + iBand] = pixVal;
  712.         } else if(pixVal < MinMaxTab[segP[MaskVal]][iBand]){
  713.           MinMaxTab[segP[MaskVal]][iBand] = pixVal;
  714.          
  715.         } else  if(pixVal > MinMaxTab[segP[MaskVal]][nbrBands + iBand]) {
  716.          
  717.           MinMaxTab[segP[MaskVal]][nbrBands + iBand] = pixVal ;
  718.          
  719.         }
  720.        
  721.           }
  722.         }
  723.      
  724.      
  725.         if(SpatTab[segP[MaskVal]][0] != -1 && LastPxl[MaskVal] == row * nXSize + col){
  726.          
  727.          
  728.  
  729.           /*
  730.  
  731.         Compute the stats and
  732.         print them to the outut file
  733.                
  734.  
  735.           */
  736.          
  737.  
  738.           if( bCalcStd){           
  739.         for(i = 0 ; i < nbrBands ; i++){
  740.           SumOfSquares = MeanStdTab[segP[MaskVal]][nbrBands + i];
  741.           SquaredSum = MeanStdTab[segP[MaskVal]][i] * MeanStdTab[segP[MaskVal]][i];
  742.           NbrObs = SpatTab[segP[MaskVal]][0];
  743.           if(NbrObs > 1)
  744.             MeanStdTab[segP[MaskVal]][nbrBands + i] =  sqrt((SumOfSquares - SquaredSum / NbrObs)/(NbrObs-1));
  745.           else
  746.             MeanStdTab[segP[MaskVal]][nbrBands + i] =  0 ;
  747.         }
  748.           }
  749.          
  750.  
  751.           if(bCalcAvg)
  752.         for(i = 0 ; i < nbrBands ; i++)
  753.           MeanStdTab[segP[MaskVal]][i] = MeanStdTab[segP[MaskVal]][i] / SpatTab[segP[MaskVal]][0];
  754.          
  755.          
  756.  
  757.           fprintf(ASCDataBase,"%i %i ",MaskVal,SpatTab[segP[MaskVal]][0]);
  758.        
  759.  
  760.           if(bCalcMinMax)
  761.         for(iBand = 0 ; iBand < 2 * nbrBands ; iBand++)
  762.           fprintf(ASCDataBase,"%f ",MinMaxTab[segP[MaskVal]][iBand]);
  763.          
  764.           if(bCalcAvg)
  765.         for(iBand = 0 ; iBand < nbrBands ; iBand++)
  766.           fprintf(ASCDataBase,"%f ",MeanStdTab[segP[MaskVal]][iBand]);
  767.          
  768.           if(bCalcStd)
  769.         for(iBand = 0 ; iBand < nbrBands ; iBand++)
  770.           fprintf(ASCDataBase,"%f ",MeanStdTab[segP[MaskVal]][nbrBands + iBand]);
  771.          
  772.           fprintf(ASCDataBase,"\n");
  773.          
  774.           SpatTab[segP[MaskVal]][0] = -1;
  775.         }
  776.        
  777.       }
  778.     }
  779.       }
  780.  
  781.  
  782.  
  783.       GDALClose( hDataset );
  784.       CSLDestroy( argv );
  785.       free(segP);
  786.       free(pafScanline);
  787.  
  788.       if(bUseMask  == TRUE)
  789.     free(pafMaskline);
  790.  
  791.       printf("\nDone\n");
  792.       fflush(stdout);
  793.       exit(0);
  794. }
  795.      
  796.  
  797. /************************************************************************/
  798. /*                        GDALInfoReportCorner()                        */
  799. /************************************************************************/
  800.  
  801. static int
  802. GDALInfoReportCorner( GDALDatasetH hDataset,
  803.                       const char * corner_name,
  804.                       double x, double y )
  805.  
  806. {
  807.     double  dfGeoX, dfGeoY;
  808.     const char  *pszProjection;
  809.     double  adfGeoTransform[6];
  810.     OGRCoordinateTransformationH hTransform = NULL;
  811.        
  812.     printf( "%-11s ", corner_name );
  813.    
  814.  
  815. /* -------------------------------------------------------------------- */
  816. /*      Transform the point into georeferenced coordinates.             */
  817. /* -------------------------------------------------------------------- */
  818.     if( GDALGetGeoTransform( hDataset, adfGeoTransform ) == CE_None )
  819.     {
  820.         pszProjection = GDALGetProjectionRef(hDataset);
  821.  
  822.         dfGeoX = adfGeoTransform[0] + adfGeoTransform[1] * x
  823.             + adfGeoTransform[2] * y;
  824.         dfGeoY = adfGeoTransform[3] + adfGeoTransform[4] * x
  825.             + adfGeoTransform[5] * y;
  826.     }
  827.  
  828.     else
  829.     {
  830.         printf( "(%7.1f,%7.1f)\n", x, y );
  831.         return FALSE;
  832.     }
  833.  
  834. /* -------------------------------------------------------------------- */
  835. /*      Report the georeferenced coordinates.                           */
  836. /* -------------------------------------------------------------------- */
  837.     if( ABS(dfGeoX) < 181 && ABS(dfGeoY) < 91 )
  838.     {
  839.         printf( "(%12.7f,%12.7f) ", dfGeoX, dfGeoY );
  840.  
  841.     }
  842.     else
  843.     {
  844.         printf( "(%12.3f,%12.3f) ", dfGeoX, dfGeoY );
  845.     }
  846.  
  847. /* -------------------------------------------------------------------- */
  848. /*      Setup transformation to lat/long.                               */
  849. /* -------------------------------------------------------------------- */
  850.     if( pszProjection != NULL && strlen(pszProjection) > 0 )
  851.     {
  852.         OGRSpatialReferenceH hProj, hLatLong = NULL;
  853.  
  854.         hProj = OSRNewSpatialReference( pszProjection );
  855.         if( hProj != NULL )
  856.             hLatLong = OSRCloneGeogCS( hProj );
  857.  
  858.         if( hLatLong != NULL )
  859.         {
  860.             CPLPushErrorHandler( CPLQuietErrorHandler );
  861.             hTransform = OCTNewCoordinateTransformation( hProj, hLatLong );
  862.             CPLPopErrorHandler();
  863.            
  864.             OSRDestroySpatialReference( hLatLong );
  865.         }
  866.  
  867.         if( hProj != NULL )
  868.             OSRDestroySpatialReference( hProj );
  869.     }
  870.  
  871. /* -------------------------------------------------------------------- */
  872. /*      Transform to latlong and report.                                */
  873. /* -------------------------------------------------------------------- */
  874.     if( hTransform != NULL
  875.         && OCTTransform(hTransform,1,&dfGeoX,&dfGeoY,NULL) )
  876.     {
  877.        
  878.         printf( "(%s,", GDALDecToDMS( dfGeoX, "Long", 2 ) );
  879.         printf( "%s)", GDALDecToDMS( dfGeoY, "Lat", 2 ) );
  880.     }
  881.  
  882.     if( hTransform != NULL )
  883.         OCTDestroyCoordinateTransformation( hTransform );
  884.    
  885.     printf( "\n" );
  886.  
  887.     return TRUE;
  888. }
RAW Paste Data
Top