main.c (e8ae190) - plGeoAdaptels · Pawel Netzel

main.c

/****************************************************************************
 *
 * PROGRAM:    plGeoadAptels
 * AUTHOR(S):    Pawel Netzel
 * PURPOSE:    program for creating adaptels using multiband spatial information;
 *        The code is based on the paper: "Scale-Adaptive Superpixels", 
 *        R. Achanta and P Marquez-Neila and P. Fua and S. Süsstrunk, 
 *        Color and Imaging Conference (CIC26), 2018 and ideas from software:
 *        https://github.com/achanta/Adaptels
 * COPYRIGHT:    (C) Pawel Netzel, Department of Forestry,
 *        University of Agriculture, Krakow, Poland
 *              http://wl.urk.edu.pl
 *
 *        This program is free software under 
 *        the GNU General Public License (>=v3). 
 *        https://www.gnu.org/licenses/gpl-3.0.en.html
 *
 *****************************************************************************/

#include <stdlib.h>
#include <stdio.h>
#include <assert.h>
#include <string.h>
#include <math.h>
#include <time.h>

#include <argtable3.h>
#include <ezgdal.h>

#include "generate.h"

void usage(char *progname, void *argtable) {
      printf("\nUsage:\n\t%s", progname);
      arg_print_syntax(stdout,argtable,"\n");
      printf("\n");
      arg_print_glossary_gnu(stdout,argtable);
      printf("\n");
      exit(0);
}

int main(int argc, char **argv) {

long i;
    char* dist_names[] = {"minkowski","cosine","angular"};

double thr_val = 60.0;
    int con_val = 4;
    int distance_val = 0;
    double minkowski_p_val = 2.0;

struct arg_str  *inp  = arg_strn("i","input","<file_name>",1,99,"name of input file(s) (GeoTIFF)");
    struct arg_str  *out  = arg_str1("o","output","<file_name>","name of output file (GeoTIFF) with adaptels");
    struct arg_dbl  *thr  = arg_dbl0("t","threshold","<val>","threshold value for adaptels algorithm (default: 60.0)");
    struct arg_str  *dist = arg_str0("d","distance","<distance name>","name of the distance measure: minkowki, cosine, angular (default: minkowski)");
    struct arg_dbl  *mink = arg_dbl0("p","minkowski_p","<val>","value for minkowski distance (default: 2.0)");
    struct arg_lit  *con  = arg_lit0("8","Queen topology","cells connectivity (default: 4, Rook)");
    struct arg_lit  *norm = arg_lit0("n","normalize","normalize cells values to the interval [0,1] (default: no)");
    struct arg_lit  *help = arg_lit0("h","help","print this help and exit");
    struct arg_end  *end  = arg_end(20);
    void* argtable[] = {inp,out,thr,dist,mink,con,norm,help,end};

int nerrors = arg_parse(argc,argv,argtable);

if (help->count > 0) 
      usage(argv[0],argtable);

if (nerrors > 0) {
      arg_print_errors(stdout,end,argv[0]);
      usage(argv[0],argtable);
    }

if(thr->count>0) {
      thr_val = thr->dval[0];
      if(thr_val < 0.0) {
        printf("\nThreshold should be grater than zero! The default value (60.0) is used.\n\n");
        thr_val = 60.0;
      }
    }
    
    if(mink->count>0)
      minkowski_p_val = mink->dval[0];

if(dist->count>0) {
      if(strcmp(dist_names[1],dist->sval[0])==0)
        distance_val = 1;
      else if(strcmp(dist_names[2],dist->sval[0])==0)
        distance_val = 2;
    }

if(con->count>0)
      con_val = 8;

ezgdal_show_message_nolf(stderr,"Opening input files:  ");

for(i=0; i<inp->count; i++)
      if(!ezgdal_file_exists((char *)(inp->sval[i]))) {
        char buf[10000];
        sprintf(buf,"File '%s' does not exists!\n\n", inp->sval[i]);
        ezgdal_show_message(stderr, buf);
        exit(0);
      }

EZGDAL_LAYER **input_layers;
    input_layers = (EZGDAL_LAYER **)malloc(inp->count*sizeof(EZGDAL_LAYER *));

for(i=0; i<inp->count; i++) {
      input_layers[i] = ezgdal_open_layer((char *)(inp->sval[i]));
      if(input_layers[i]==NULL) {
        char buf[10000];
        sprintf(buf,"Can not open file: '%s'\n\n", inp->sval[i]);
        ezgdal_show_message(stderr, buf);
        exit(0);
      }
    }

if(!ezgdal_is_projection_ok(input_layers, inp->count)) 
      ezgdal_show_message(stderr,"Input files have various projections!\n\n");

if(!ezgdal_is_bbox_ok(input_layers,inp->count))
      ezgdal_show_message(stderr,"Input files have various extents and/or resolutions!\n\n");

ezgdal_show_message(stderr,"OK");

int nodata = -9999;
    ezgdal_show_message_nolf(stderr,"Creating output file: ");
    EZGDAL_LAYER *output = ezgdal_create_layer((char *)out->sval[0],
                                               ezgdal_layer_get_wkt(input_layers[0]),
                                               "Int32",
                                               ezgdal_layer_get_at(input_layers[0]),
                                               input_layers[0]->rows,
                                               input_layers[0]->cols,
                                               &nodata, EZGDAL_COMPRESS_DEFLATE
                                              );
    ezgdal_show_message(stderr,"OK");

ezgdal_show_message_nolf(stderr,"Input data reading:   ");
    long buffer_size = (long)input_layers[0]->cols*(long)input_layers[0]->rows;
    double **input_buffers = (double **)malloc(inp->count*sizeof(double *));
    for(i=0; i<inp->count; i++) {
      input_buffers[i] = (double *)malloc(buffer_size*sizeof(double));
      ezgdal_read_all_layer(input_layers[i], input_buffers[i]);
    }
    unsigned char *mask = (unsigned char *)calloc(buffer_size, sizeof(char));
    for(int c=0; c<inp->count; c++)
      for(i=1; i<buffer_size; i++)
        if(ezgdal_is_null(input_layers[c],input_buffers[c][i]))
          mask[i] = 1; // 1 - nodata, 0 - value

ezgdal_show_message(stderr,"OK");

if(norm->count>0) {
      ezgdal_show_message_nolf(stderr,"Data normalization:   ");
      for(int c=0; c<inp->count; c++) {
        double min = input_buffers[c][0];
        double max = min;
        double v;
        for(i=1; i<buffer_size; i++) {
          v = input_buffers[c][i];
          if(mask[i]==0) {
            if(v<min) min = v;
            if(v>max) max = v;
          }
        }
        if(min==max)
          for(i=1; i<buffer_size; i++) {
            if(mask[i]==0)
              input_buffers[c][i] = 0.0;
          }
        else
          for(i=1; i<buffer_size; i++) {
            v = input_buffers[c][i];
            if(mask[i]==0)
              input_buffers[c][i] = (v - min)/(max - min);
          }
      }
      ezgdal_show_message(stderr,"OK");
    }
    
    ezgdal_show_message_nolf(stderr,"Distance:             ");
    ezgdal_show_message(stderr,dist_names[distance_val]);
    ezgdal_show_message_nolf(stderr,"Adaptels creating:    ");

clock_t time_start, time_end;
    double time_elapsed;

time_start = clock();
    LABELS* labels = create_adaptels(input_buffers, inp->count, mask,
                                     input_layers[0]->cols, 
                                     input_layers[0]->rows, 
                                     thr_val, con_val,
                                     distance_val,
                                     minkowski_p_val);
    time_end = clock();
    time_elapsed = ((double) (time_end - time_start)) / CLOCKS_PER_SEC;
    fprintf(stderr,"%d adaptels created in %.4lf seconds\n",labels->max_label, time_elapsed);

ezgdal_show_message_nolf(stderr,"Results writting:     ");

for(i=0; i<buffer_size; i++)
        if(mask[i])
          labels->labels[i] = nodata;

for(i=0; i<inp->count; i++) 
      ezgdal_close_layer(input_layers[i]);
    free(input_layers);

for(i=0; i<inp->count; i++) 
      free(input_buffers[i]);
    free(input_buffers);

double *output_buffer = (double *)malloc(buffer_size*sizeof(double));
    for(i=0; i<buffer_size; i++) 
      output_buffer[i] = labels->labels[i];

ezgdal_write_all_layer(output, output_buffer);
    ezgdal_close_layer(output);
    free(output_buffer);

ezgdal_show_message(stderr,"OK");

free(mask);
    free_labels(labels);

return 0;

}

plGeoAdaptels

new distances