Geografic-Information-System/raster/r.texture/main.c

/****************************************************************************
 *
 * MODULE:       r.texture
 * AUTHOR(S):    Carmine Basco - basco@unisannio.it
 *               with hints from: 
 * 			prof. Giulio Antoniol - antoniol@ieee.org
 * 			prof. Michele Ceccarelli - ceccarelli@unisannio.it
 *               Markus Metz (optimization and bug fixes)
 *
 * PURPOSE:      Create map raster with textural features.
 *
 * COPYRIGHT:    (C) 2003 by University of Sannio (BN), Benevento, Italy 
 *
 *               This program is free software under the GNU General Public
 *   	    	 License (>=v2). Read the file COPYING that comes with GRASS
 *   	    	 for details.
 *
 * Permission to use, copy, modify, and distribute this software and its
 * documentation for any purpose and without fee is hereby granted. This 
 * software is provided "as is" without express or implied warranty.
 *
 *****************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <grass/gis.h>
#include <grass/raster.h>
#include <grass/glocale.h>
#include "h_measure.h"

struct menu
{
    char *name;			/* measure name */
    char *desc;			/* menu display - full description */
    char *suffix;		/* output suffix */
    char useme;			/* calculate this measure if set */
    int idx;			/* measure index */
};

/* modify this table to add new measures */
static struct menu menu[] = {
    {"asm",      "Angular Second Moment",    "_ASM",   0,  1},
    {"contrast", "Contrast",                 "_Contr", 0,  2},
    {"corr",     "Correlation",              "_Corr",  0,  3},
    {"var",      "Variance",                 "_Var",   0,  4},
    {"idm",      "Inverse Diff Moment",      "_IDM",   0,  5},
    {"sa",       "Sum Average",              "_SA",    0,  6},
    {"sv",       "Sum Variance",             "_SV",    0,  7},
    {"se",       "Sum Entropy",              "_SE",    0,  8},
    {"entr",     "Entropy",                  "_Entr",  0,  9},
    {"dv",       "Difference Variance",      "_DV",    0, 10},
    {"de",       "Difference Entropy",       "_DE",    0, 11},
    {"moc1",     "Measure of Correlation-1", "_MOC-1", 0, 12},
    {"moc2",     "Measure of Correlation-2", "_MOC-2", 0, 13},
    {NULL, NULL, NULL, 0, -1}
};

static int find_measure(const char *measure_name)
{
    int i;

    for (i = 0; menu[i].name; i++)
	if (strcmp(menu[i].name, measure_name) == 0)
	    return i;

    G_fatal_error(_("Unknown measure <%s>"), measure_name);

    return -1;
}

int main(int argc, char *argv[])
{
    struct Cell_head cellhd;
    char *name, *result;
    char **mapname;
    FCELL **fbuf;
    int n_measures, n_outputs, *measure_idx, overwrite;
    int nrows, ncols;
    int row, col, first_row, last_row, first_col, last_col;
    int i, j;
    CELL **data;		/* Data structure containing image */
    DCELL *dcell_row;
    struct FPRange range;
    DCELL min, max, inscale;
    FCELL measure;		/* Containing measure done */
    int dist, size;	/* dist = value of distance, size = s. of moving window */
    int offset;
    int have_px, have_py, have_pxpys, have_pxpyd;
    int infd, *outfd;
    RASTER_MAP_TYPE out_data_type;
    struct GModule *module;
    struct Option *opt_input, *opt_output, *opt_size, *opt_dist, *opt_measure;
    struct Flag *flag_ind, *flag_all, *flag_null;
    struct History history;
    char p[1024];

    G_gisinit(argv[0]);

    module = G_define_module();
    G_add_keyword(_("raster"));
    G_add_keyword(_("algebra"));
    G_add_keyword(_("statistics"));
    G_add_keyword(_("texture"));
    module->description =
	_("Generate images with textural features from a raster map.");

    /* Define the different options */

    opt_input = G_define_standard_option(G_OPT_R_INPUT);

    opt_output = G_define_standard_option(G_OPT_R_BASENAME_OUTPUT);

    opt_size = G_define_option();
    opt_size->key = "size";
    opt_size->key_desc = "value";
    opt_size->type = TYPE_INTEGER;
    opt_size->required = NO;
    opt_size->description = _("The size of moving window (odd and >= 3)");
    opt_size->answer = "3";

    /* Textural character is in direct relation of the spatial size of the texture primitives. */

    opt_dist = G_define_option();
    opt_dist->key = "distance";
    opt_dist->key_desc = "value";
    opt_dist->type = TYPE_INTEGER;
    opt_dist->required = NO;
    opt_dist->label = _("The distance between two samples (>= 1)");
    opt_dist->description = _("The distance must be smaller than the size of the moving window");
    opt_dist->answer = "1";

    for (i = 0; menu[i].name; i++) {
	if (i)
	    strcat(p, ",");
	else
	    *p = 0;
	strcat(p, menu[i].name);
    }
    opt_measure = G_define_option();
    opt_measure->key = "method";
    opt_measure->type = TYPE_STRING;
    opt_measure->required = NO;
    opt_measure->multiple = YES;
    opt_measure->options = p;
    opt_measure->description = _("Textural measurement method");

    flag_ind = G_define_flag();
    flag_ind->key = 's';
    flag_ind->label = _("Separate output for each angle (0, 45, 90, 135)");
    flag_ind->description = _("Angles are counterclockwise from east: "
                              "0 is East to West, 45 is North-East to South-West");

    flag_all = G_define_flag();
    flag_all->key = 'a';
    flag_all->description = _("Calculate all textural measurements");

    flag_null = G_define_flag();
    flag_null->key = 'n';
    flag_null->label = _("Allow NULL cells in a moving window");
    flag_null->description = _("This will also avoid cropping along edges of the current region");

    if (G_parser(argc, argv))
	exit(EXIT_FAILURE);

    name = opt_input->answer;
    result = opt_output->answer;
    size = atoi(opt_size->answer);
    if (size <= 0)
	G_fatal_error(_("Size of the moving window must be > 0"));
    if (size % 2 != 1)
	G_fatal_error(_("Size of the moving window must be odd"));
    dist = atoi(opt_dist->answer);
    if (dist <= 0)
	G_fatal_error(_("The distance between two samples must be > 0"));
    if (dist >= size)
	G_fatal_error(_("The distance between two samples must be smaller than the size of the moving window"));

    n_measures = 0;
    if (flag_all->answer) {
	for (i = 0; menu[i].name; i++) {
	    menu[i].useme = 1;
	}
	n_measures = i;
    }
    else {
	for (i = 0; opt_measure->answers[i]; i++) {
	    if (opt_measure->answers[i]) {
		const char *measure_name = opt_measure->answers[i];
		int n = find_measure(measure_name);

		menu[n].useme = 1;
		n_measures++;
	    }
	}
    }
    if (!n_measures)
	G_fatal_error(_("Nothing to compute. Use at least one textural measure."));
	
    measure_idx = G_malloc(n_measures * sizeof(int));
    j = 0;
    for (i = 0; menu[i].name; i++) {
	if (menu[i].useme == 1) {
	    measure_idx[j] = i;
	    j++;
	}
    }

    /* variables needed */
    if (menu[2].useme || menu[11].useme || menu[12].useme)
	have_px = 1;
    else
	have_px = 0;
    if (menu[11].useme || menu[12].useme)
	have_py = 1;
    else
	have_py = 0;
    if (menu[5].useme || menu[6].useme || menu[7].useme)
	have_pxpys = 1;
    else
	have_pxpys = 0;
    if (menu[9].useme || menu[10].useme)
	have_pxpyd = 1;
    else
	have_pxpyd = 0;

    infd = Rast_open_old(name, "");

    Rast_get_cellhd(name, "", &cellhd);

    out_data_type = FCELL_TYPE;
    /* Allocate output buffers, use FCELL data_type */
    n_outputs = n_measures;
    if (flag_ind->answer) {
	n_outputs = n_measures * 4;
    }

    fbuf = G_malloc(n_outputs * sizeof(FCELL *));
    mapname = G_malloc(n_outputs * sizeof(char *));
    for (i = 0; i < n_outputs; i++) {
	mapname[i] = G_malloc(GNAME_MAX * sizeof(char));
	fbuf[i] = Rast_allocate_buf(out_data_type);
    }
    
    overwrite = G_check_overwrite(argc, argv);

    /* open output maps */
    outfd = G_malloc(n_outputs * sizeof(int));
    for (i = 0; i < n_measures; i++) {
	if (flag_ind->answer) {
	    for (j = 0; j < 4; j++) {
		sprintf(mapname[i * 4 + j], "%s%s_%d", result,
		        menu[measure_idx[i]].suffix, j * 45);
		if (!G_find_raster(mapname[i * 4 + j], G_mapset()) || overwrite) {
			outfd[i * 4 + j] = Rast_open_new(mapname[i * 4 + j], out_data_type);
	        } else {
		        G_fatal_error(_("At least one of the requested output maps exists. Use --o to overwrite."));
	        }
	    }
	}
	else {
	    sprintf(mapname[i], "%s%s", result,
	            menu[measure_idx[i]].suffix);
	    if (!G_find_raster(mapname[i], G_mapset()) || overwrite) {
		    outfd[i] = Rast_open_new(mapname[i], out_data_type);
	    } else {
		    G_fatal_error(_("At least one of the requested output maps exists. Use --o to overwrite."));
	    }
	}
    }
    nrows = Rast_window_rows();
    ncols = Rast_window_cols();

    /* Load raster map. */

    /* allocate the space for one row of cell map data *A* */
    dcell_row = Rast_allocate_d_buf();

    /* Allocate appropriate memory for the structure containing the image */
    data = (int **)G_malloc(nrows * sizeof(int *));
    for (i = 0; i < nrows; i++) {
	data[i] = (int *)G_malloc(ncols * sizeof(int));
    }

    /* read input range */
    Rast_init_fp_range(&range);
    Rast_read_fp_range(name, "", &range);
    Rast_get_fp_range_min_max(&range, &min, &max);
    inscale = 0;
    if (min < 0 || max > 255) {
	inscale = 255. / (max - min);
    }
    /* input has 0 - 1 range */
    else if (max <= 1.) {
	inscale = 255. / (max - min);
    }

    /* Read in cell map values */
    /* TODO: use r.proj cache */
    G_important_message(_("Reading raster map..."));
    for (j = 0; j < nrows; j++) {
	Rast_get_row(infd, dcell_row, j, DCELL_TYPE);
	for (i = 0; i < ncols; i++) {
	    if (Rast_is_d_null_value(&(dcell_row[i])))
		data[j][i] = -1;
	    else if (inscale) {
		data[j][i] = (CELL)((dcell_row[i] - min) * inscale);
	    }
	    else
		data[j][i] = (CELL)dcell_row[i];
	}
    }

    /* close input cell map and release the row buffer */
    Rast_close(infd);
    G_free(dcell_row);

    /* Now raster map is loaded to memory. */

    /* *************************************************************************************************
     *
     * Compute of the matrix S.G.L.D. (Spatial Gray-Level Dependence Matrices) or co-occurrence matrix.
     * The image is analized for piece, every piece is naming moving window (s.w.). The s.w. must be    
     * square with number of size's samples odd, that because we want the sample at the center of matrix. 
     *
     ***************************************************************************************************/

    offset = size / 2;

    if (!flag_null->answer) {
	first_row = first_col = offset;
	last_row = nrows - offset;
	last_col = ncols - offset;
    }
    else {
	/* no cropping at window margins */
	first_row = first_col = 0;
	last_row = nrows;
	last_col = ncols;
    }

    Rast_set_f_null_value(fbuf[0], ncols);

    for (row = 0; row < first_row; row++) {
	for (i = 0; i < n_outputs; i++) {
	    Rast_put_row(outfd[i], fbuf[0], out_data_type);
	}
    }
    if (n_measures > 1)
	G_message(n_("Calculating %d texture measure", 
        "Calculating %d texture measures", n_measures), n_measures);
    else
	G_message(_("Calculating %s..."), menu[measure_idx[0]].desc);
    alloc_vars(size);
    for (row = first_row; row < last_row; row++) {
	G_percent(row, nrows, 2);

	for (i = 0; i < n_outputs; i++)
	    Rast_set_f_null_value(fbuf[i], ncols);

	/*process the data */
	for (col = first_col; col < last_col; col++) {

	    if (!set_vars(data, row, col, size, offset, dist, flag_null->answer)) {
		for (i = 0; i < n_outputs; i++)
		    Rast_set_f_null_value(&(fbuf[i][col]), 1);
		continue;
	    }

	    /* for all angles (0, 45, 90, 135) */
	    for (i = 0; i < 4; i++) {
		set_angle_vars(i, have_px, have_py, have_pxpys, have_pxpyd);
		/* for all requested textural measures */
		for (j = 0; j < n_measures; j++) {

		    measure = (FCELL) h_measure(menu[measure_idx[j]].idx);

		    if (flag_ind->answer) {
			/* output for each angle separately */
			fbuf[j * 4 + i][col] = measure;
		    }
		    else {
			/* use average over all angles for each measure */
			if (i == 0)
			    fbuf[j][col] = measure;
			else if (i < 3)
			    fbuf[j][col] += measure;
			else 
			    fbuf[j][col] = (fbuf[j][col] + measure) / 4.0;
		    }
		}
	    }
	}
	for (i = 0; i < n_outputs; i++) {
	    Rast_put_row(outfd[i], fbuf[i], out_data_type);
	}
    }
    Rast_set_f_null_value(fbuf[0], ncols);
    for (row = last_row; row < nrows; row++) {
	for (i = 0; i < n_outputs; i++) {
	    Rast_put_row(outfd[i], fbuf[0], out_data_type);
	}
    }
    G_percent(nrows, nrows, 1);

    for (i = 0; i < n_outputs; i++) {
	Rast_close(outfd[i]);

	Rast_short_history(mapname[i], "raster", &history);
	Rast_command_history(&history);
	Rast_write_history(mapname[i], &history);
	G_free(fbuf[i]);
    }

    G_free(fbuf);
    G_free(data);

    exit(EXIT_SUCCESS);
}