main.c (1ed00b7) - plForestFragment · Pawel Netzel

main.c

/****************************************************************************
 *
 * PROGRAM:      forestfragment
 *
 * AUTHOR(S):    Pawel Netzel - URK
 *
 * PURPOSE:      Program calculates multiscale fores fragmentation measure
 *               with using Jensen-Shannon Similarity measuring similarity
 *               to fully forested areas.
 *
 * COPYRIGHT:    (C) 2023 by Pawel Netzel
 *
 *               This program is a free software under the GNU General Public
 *               License (>=v2).
 *
 ***************************************************************************/

#include <stdlib.h>
#include <string.h>
#include <omp.h>
#include <gdal.h>
#include <ezgdal.h>
#include <argtable3.h>

/****
 *  The procedure displays program usage 
 *  and returns to the system
 */
void usage(char *progname, void *argtable) {

printf("\nplForestFragment calculates forest fragmentation index based on \n");
	printf("similarity approach.\n");
	printf("It can calculate the index for one scale. It can also produce \n");
	printf("the multiscale fragmentation index. \n");
	printf("plForestFragment takes a categorical raster layer as an input.\n");
	printf("A user can point a category that represents the forest class.\n");

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

/****
 *  Comparator to sort radiuses
 */
int compare(const void* numA, const void* numB)
{
    const int* num1 = (const int*)numA;
    const int* num2 = (const int*)numB;

if (*num1 > *num2) {
        return 1;
    }
    else {
        if (*num1 == *num2)
            return 0;
        else
            return -1;
    }
}

/****
 *  Calculating index in the coocurence matrix 
 *  based on values of categories
 */
int calc_cooc(double v1, double v2, double forest) {
    if(v1==forest && v2==forest) return 0;
    else if(v1==forest) return 1;
    else if(v2==forest) return 2;
    else return 3;
}

/****
 *  Creating multiscale matrix that will be used as a mask
 *  to get circular neighbourhood
 */
int** create_mask(int* rad, int count) {
    int radius, radius_max=rad[count-1];
    int size = 2*radius_max+1;
    int i, r, c;
    int** mask = (int**)malloc(size*sizeof(int*));
    double x, y;

for(i=0; i<size; i++)
        mask[i] = (int*)calloc(size, sizeof(int));
    for(r=0; r<size; r++) {
        y = pow(r - radius_max, 2);
        for(c=0; c<size; c++) {
            x = pow(c - radius_max, 2);
            if(sqrt(x+y)<=radius_max)
                mask[r][c] = count-1;
            else
                mask[r][c] = count+1;
        }
    }
    for(i=count-2; i>-1; i--) {
        radius = rad[i];
        for(r=0; r<size; r++) {
            y = pow(r - radius_max, 2);
            for(c=0; c<size; c++) {
                x = pow(c - radius_max, 2);
                if(sqrt(x+y)<=radius)
                    mask[r][c] = i;
            }
        }
    }
    return mask;
}

void free_mask(int** mask, int size) {
    int i;
    for(i=0; i<size; i++) free(mask[i]);
    free(mask);
}

/****
 *  The function calculates Jensen-Shannon Similarity - JSS
 *  between compund coocurence matrix and reference matrix
 */
double calc_similarity(int** cooc, int count, double weight, double* reference) {
    int i, j, pos=0;
    double tsum=0.0, sum, M=0.0, HM=0.0, HA=0.0, HB=0.0;
    double* norm =(double*)calloc(count*4,sizeof(double));
    // normalize
    for(i=0; i<count; i++) {
        sum = 0.0;
        for(j=0; j<4; j++)
            sum += cooc[i][j];
        tsum += sum;
        sum /= weight;
        for(j=0; j<4; j++)
            if(sum>0)
                norm[pos++] = (double)cooc[i][j]/sum;
    }
    if(tsum==0.0) return 0.0;
    for(i=0; i<count*4; i++) {
        M = 0.5*(norm[i] + reference[i]);
        HM -= (M>0.0)?M*log2(M):0.0;
        HA -= (norm[i]>0)?norm[i]*log2(norm[i]):0.0;
        HB -= (reference[i]>0.0)?reference[i]*log2(reference[i]):0.0;
    }
    free(norm);
    return 1.0-sqrt(HM-0.5*(HA+HB));
}

/****
 *   Main calculation procedure
 *   It fills output layer with calculated similarities.
 *   The procedure uses one or more scales defines by 
 *   variables: count and rad
 */
void calc_forest_fragmentation(EZGDAL_LAYER* in_fd, double forest, 
                               EZGDAL_LAYER* out_fd, 
                               int* rad, int count) {
    EZGDAL_STRIPE* stripe = in_fd->stripe;
    int layer_cols = in_fd->cols, layer_rows = in_fd->rows;
    double* buffer = out_fd->buffer;
    double in_no_data = in_fd->no_data;
    double out_no_data = out_fd->no_data;
    double weight = 1.0/(double)count;
    int i, row, col;
    int radius_max = rad[count-1];
    double* reference = (double*)calloc(count*4,sizeof(double));
    int** mask = create_mask(rad, count);

for(i=0; i<count; i++)
        reference[i*4] = weight;

ezgdal_show_message_nolf(stdout,"Calculations...   0%");
    /****
     *  For each row of input raster
     */
    for(row=0; row<layer_rows; row++) {
        ezgdal_load_stripe_data(stripe, row-radius_max);
        ezgdal_show_progress(stderr,row,layer_rows);
        /****
         *  For each cell int the row - parallel processing
         */
#pragma omp parallel for default(shared) schedule(dynamic, 100)
        for(col=0; col<layer_cols; col++) {
            int i, c, r;
            EZGDAL_FRAME* frame = &(stripe->frame[col]);
            double** frame_buffer = frame->buffer;
            if(frame_buffer[radius_max][radius_max]!=in_no_data) {
                int** cooc = (int **)malloc(count*sizeof(int*));
                for(i=0; i<count; i++)
                    cooc[i] = (int*)calloc(4, sizeof(int));
                /****
                 *  For each scale
                 */
                for(i=0; i<count; i++) {
                    int radius = rad[i];
                    int* cooc_scale = cooc[i];
                    /****
                     *  Calculation of the coocurence matrix
                     */
                    for(r=radius_max-radius; r<radius_max+radius; r++) {
                        for(c=radius_max-radius; c<radius_max+radius; c++) {
                            if(mask[r][c]<=i) {
                                double v = frame_buffer[r][c];
                                if(v!=in_no_data) {
                                    double v_right = frame_buffer[r][c+1];
                                    double v_down = frame_buffer[r+1][c];
                                    if(v_right!=in_no_data)
                                        cooc_scale[calc_cooc(v, v_right, forest)]++;
                                    if(v_down!=in_no_data)
                                        cooc_scale[calc_cooc(v, v_down, forest)]++;
                                }
                            }
                        }
                    }
                }
                /****
                 *  Calculate of similarity based on counpound coocurence matrix
                 */
                buffer[col] = calc_similarity(cooc, count, weight, reference);
                for(i=0; i<count; i++)
                    free(cooc[i]);
                free(cooc);
            } else
                buffer[col] = out_no_data;
        }
        ezgdal_write_buffer(out_fd, row);
    }
    ezgdal_show_progress(stderr,row,layer_rows);

free(reference);
    free_mask(mask, 2*radius_max+1);
}

int main(int argc, char *argv[])
{
    EZGDAL_LAYER *in_fd, *out_fd;
    char *OUTPUT;
    char *INPUT;
    double forest = 1.0;

struct arg_str  *inp   = arg_str1("i","input","<file_name>","name of input files (GeoTIFF)");
    struct arg_dbl  *frst  = arg_dbl0("f","forest-value","<dbl>","value for forest (default: 1.0)");
    struct arg_str  *out   = arg_str1("o","output","<file_name>","name of output file with forest fragmentation index (JSS) (GeoTIFF)");
    struct arg_int  *rad   = arg_intn("r","radius","<n>",1,8,"radius(es) in cells defining scale(s)");
    struct arg_int  *th    = arg_int0("t",NULL,"<n>","number of threads (default: 1)");
    struct arg_lit  *help  = arg_lit0("h","help","print help and exit");
    struct arg_end  *end   = arg_end(20);

void* argtable[] = {inp,frst,out,rad,th,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);
    }

/****
     *  Setting the number of threads for OpenMP - default 1
     */
    if(th->count > 0) 
      omp_set_num_threads(th->ival[0]);
    else
      omp_set_num_threads(1);

/****
     *  Setting the value for forested cells - default 1
     */
    if(frst->count > 0) 
      forest = frst->dval[0];

INPUT = (char *)inp->sval[0];
    OUTPUT = (char *)out->sval[0];

/****
     *  Sorting radiuses from the smallest to the largest
     */
    qsort(rad->ival, rad->count, sizeof(int), compare);

/****
     *  Opening input raster layer
     */
    if(!ezgdal_file_exists(INPUT)) {
      printf("\nThe file '%s' does not exist!\n\n",INPUT);
      exit(0);
    }

ezgdal_show_message_nolf(stdout,"Opening input file...");
    in_fd = ezgdal_open_layer(INPUT);

if(in_fd==NULL) {
      printf("\nThe file '%s' can't be opened!\n\n",INPUT);
      exit(0);
    } else
      ezgdal_show_message(stdout," OK");

ezgdal_create_stripe(in_fd, -rad->ival[rad->count-1], 2*rad->ival[rad->count-1]+1);
    ezgdal_create_all_frames(in_fd->stripe, -rad->ival[rad->count-1], 1);

/****
     *  Creating output raster layer
     */
    int nodata = -9999;
    ezgdal_show_message_nolf(stdout,"Creating output file...");
    out_fd = ezgdal_create_layer(OUTPUT,
                                 ezgdal_layer_get_wkt(in_fd),
                                 "Float32",
                                 ezgdal_layer_get_at(in_fd),
                                 in_fd->rows,
                                 in_fd->cols,
                                 &nodata, EZGDAL_COMPRESS_DEFLATE
                                );
    if(out_fd==NULL) {
      printf("\nThe file '%s' can't be created!\n\n",OUTPUT);
      exit(0);
    } else
      ezgdal_show_message(stdout," OK");

/****
     *  Calculation of the forest fragmentation measure
     */
    calc_forest_fragmentation(in_fd, forest, out_fd, rad->ival, rad->count);

/****
     *  Cleaning up and closing files
     */
    ezgdal_close_layer(out_fd);
    ezgdal_close_layer(in_fd);
    arg_free(argtable);

exit(EXIT_SUCCESS);
}

plForestFragment

Initial commit