Geografic-Information-System/imagery/i.segment/cluster.c

/****************************************************************************
 *
 * MODULE:       r.clump
 *
 * AUTHOR(S):    Michael Shapiro - CERL
 *               Markus Metz
 *
 * PURPOSE:      Recategorizes data in a raster map layer by grouping cells
 *               that form physically discrete areas into unique categories.
 *
 * COPYRIGHT:    (C) 2006-2016 by the GRASS Development Team
 *
 *               This program is free software under the GNU General Public
 *               License (>=v2). Read the file COPYING that comes with GRASS
 *               for details.
 *
 ***************************************************************************/

#include <sys/types.h>
#include <sys/stat.h>
#include <fcntl.h>
#include <unistd.h>
#include <math.h>
#include <time.h>
#include <grass/gis.h>
#include <grass/raster.h>
#include <grass/glocale.h>
#include "iseg.h"

#define INCR 1024

/* defeats the purpose of padding ... */
#define CISNULL(r, c) (((c) == 0 || (c) == ncols + 1 ? 1 : \
                       (FLAG_GET(globals->null_flag, (r), (c - 1)))))


CELL cluster_bands(struct globals *globals)
{
    register int col;
    register int i, n;
    /* input */
    DCELL **prev_in, **cur_in, **temp_in;
    struct ngbr_stats Ri, Rk, Rn;
    int nin;
    int diag;
    int bcol;
    /* output */
    CELL OLD, NEW;
    CELL *temp_clump, out_cell;
    CELL *prev_clump, *cur_clump;
    CELL *index, *renumber;
    CELL label, cellmax;
    int nrows, ncols;
    int row;
    int len;
    int nalloc;
    char *cname;
    int cfd, csize;
    CELL cat;
    int mwrow, mwrow1, mwrow2, mwnrows, mwcol, mwcol1, mwcol2, mwncols, radiusc;
    double diff, diff2, avgdiff, ka2, hspat, hspat2;
    double hspec, hspecad, hspec2, hspec2ad;
    LARGEINT count;

    G_message(_("%d-band clustering with threshold %g"), globals->nbands, globals->hr);

    nrows = Rast_window_rows();
    ncols = Rast_window_cols();

    hspec = globals->hr;
    hspec2 = globals->hr * globals->hr;
    nin = globals->nbands;
    diag = (globals->nn == 8);
    radiusc = globals->hs;
    mwnrows = mwncols = radiusc * 2 + 1;

    /* spatial bandwidth */
    hspat = globals->hs;
    if (hspat < 1)
	hspat = 1.5;
    hspat2 = hspat * hspat;

    cellmax = ((CELL)1 << (sizeof(CELL) * 8 - 2)) - 1;
    cellmax += ((CELL)1 << (sizeof(CELL) * 8 - 2));

    Ri.mean = NULL;
    Rk.mean = NULL;
    Rn.mean = G_malloc(globals->datasize);

    /* allocate clump index */
    /* TODO: support smallest label ID > 1 */
    nalloc = INCR;
    index = (CELL *) G_malloc(nalloc * sizeof(CELL));
    index[0] = 0;
    renumber = NULL;

    /* allocate DCELL buffers two columns larger than current window */
    prev_in = (DCELL **) G_malloc(sizeof(DCELL *) * (ncols + 2));
    cur_in = (DCELL **) G_malloc(sizeof(DCELL *) * (ncols + 2));
    
    prev_in[0] = (DCELL *) G_malloc(globals->datasize * (ncols + 2) * nin);
    cur_in[0] = (DCELL *) G_malloc(globals->datasize * (ncols + 2) * nin);

    Rast_set_d_null_value(cur_in[0], (ncols + 2) * nin);
    Rast_set_d_null_value(prev_in[0], (ncols + 2) * nin);

    for (i = 1; i < ncols + 2; i++) {
	prev_in[i] = prev_in[i - 1] + nin; 
	cur_in[i] = cur_in[i - 1] + nin; 
    }

    /* allocate CELL buffers two columns larger than current window */
    len = (ncols + 2) * sizeof(CELL);
    prev_clump = (CELL *) G_malloc(len);
    cur_clump = (CELL *) G_malloc(len);

    /* temp file for initial clump IDs */
    cname = G_tempfile();
    if ((cfd = open(cname, O_RDWR | O_CREAT | O_EXCL, 0600)) < 0)
	G_fatal_error(_("Unable to open temp file"));
    csize = ncols * sizeof(CELL);

    /* initialize clump labels */
    G_zero(cur_clump, len);
    G_zero(prev_clump, len);
    /* TODO: support smallest label ID > 1 */
    label = 0;

    /****************************************************
     *                      PASS 1                      *
     * pass thru the input, create initial clump labels *
     ****************************************************/

    G_message(_("Assigning initial region IDs..."));
    for (row = 0; row < nrows; row++) {
	G_percent(row, nrows, 2);

	for (col = 1; col <= ncols; col++) {

	    /* get band values */
	    Segment_get(globals->bands_out, (void *)cur_in[col], row, col - 1);
	    Ri.mean = cur_in[col];

	    if (CISNULL(row, col)) {
		cur_clump[col] = 0;
		continue;
	    }
	    
	    hspec2ad = hspec2;

	    if (globals->ms_adaptive) {
		/* adapt initial range bandwidth */

		mwrow1 = row - radiusc;
		mwrow2 = mwrow1 + mwnrows;
		if (mwrow1 < 0)
		    mwrow1 = 0;
		if (mwrow2 > nrows)
		    mwrow2 = nrows;

		mwcol1 = col - radiusc;
		mwcol2 = mwcol1 + mwncols;
		if (mwcol1 < 0)
		    mwcol1 = 0;
		if (mwcol2 > ncols)
		    mwcol2 = ncols;

		ka2 = hspec2; 	/* OTB: conductance parameter */
		
		avgdiff = 0;
		count = 0;
		for (mwrow = mwrow1; mwrow < mwrow2; mwrow++) {
		    for (mwcol = mwcol1; mwcol < mwcol2; mwcol++) {
			if ((FLAG_GET(globals->null_flag, mwrow, mwcol)))
			    continue;
			if (mwrow == row && mwcol == col)
			    continue;

			diff = mwrow - row;
			diff2 = diff * diff;
			diff = mwcol - col;
			diff2 += diff * diff;

			if (diff2 <= hspat2) {

			    Segment_get(globals->bands_out, (void *)Rn.mean,
					mwrow, mwcol);

			    /* get spectral distance */
			    diff2 = (globals->calculate_similarity)(&Ri, &Rn, globals);

			    avgdiff += sqrt(diff2);
			    count++;
			}
		    }
		}
		hspec2ad = 0;
		if (avgdiff > 0) {
		    avgdiff /= count;
		    hspecad = hspec;
		    /* OTB-like, contrast enhancing */
		    hspecad = exp(-avgdiff * avgdiff / (2 * ka2)) * avgdiff;
		    /* preference for large regions, from Perona Malik 1990 
		     * if the settings are right, it could be used to reduce noise */
		    /* hspecad = 1 / (1 + (avgdiff * avgdiff / (2 * hspec2))); */
		    hspec2ad = hspecad * hspecad;
		    G_debug(1, "avg spectral diff: %g", avgdiff);
		    G_debug(1, "initial hspec2: %g", hspec2);
		    G_debug(1, "adapted hspec2: %g", hspec2ad);
		}
	    }

	    /*
	     * if the cell values are different to the left and above
	     * (diagonal: and above left and above right)
	     * then we must start a new clump
	     *
	     * this new clump may eventually collide with another
	     * clump and will have to be merged
	     */

	    /* try to connect the current cell to an existing clump */
	    OLD = NEW = 0;

	    Rk.mean = cur_in[col - 1];

	    /* same clump as to the left */
	    if (!CISNULL(row, col - 1) &&
	        globals->calculate_similarity(&Ri, &Rk, globals) <= hspec2ad) {
		OLD = cur_clump[col] = cur_clump[col - 1];
	    }

	    if (diag) {
		/* check above right, center, left, in that order */
		temp_clump = prev_clump + col + 1;
		bcol = col + 1;
		do {
		    Rk.mean = prev_in[bcol];
		    if (row > 0 && !CISNULL(row - 1, bcol) &&
		        globals->calculate_similarity(&Ri, &Rk, globals) <= hspec2ad) {
			cur_clump[col] = *temp_clump;
			if (OLD == 0) {
			    OLD = *temp_clump;
			}
			else {
			    NEW = *temp_clump;

			    /* threshold > 0 and diagonal requires a bit of extra work
			     * because of bridge cells:
			     * A similar to B, B similar to C, but A not similar to C
			     * -> B is bridge cell */
			    if (NEW != OLD) {
				CELL *temp_clump2;

				/* conflict! preserve NEW clump ID and change OLD clump ID.
				 * Must go back to the left in the current row and to the right
				 * in the previous row to change all the clump values as well.
				 */

				/* left of the current row from 1 to col - 1 */
				temp_clump2 = cur_clump;
				n = col - 1;
				while (n-- > 0) {
				    temp_clump2++;	/* skip left edge */
				    if (*temp_clump2 == OLD)
					*temp_clump2 = NEW;
				}

				/* right of previous row from col - 1 to ncols */
				temp_clump2 = prev_clump + col - 1;
				n = ncols - col + 2;
				while (n-- > 0) {
				    if (*temp_clump2 == OLD)
					*temp_clump2 = NEW;
				    temp_clump2++;
				}

				/* modify the OLD index */
				index[OLD] = NEW;

				OLD = NEW;
				NEW = 0;
			    }
			}
		    }
		    temp_clump--;
		} while (bcol-- > col - 1);
	    }
	    else {
		/* check above */
		Rk.mean = prev_in[col];

		if (row > 0 && !CISNULL(row - 1, col) &&
		    globals->calculate_similarity(&Ri, &Rk, globals) <= hspec2ad) {
		    temp_clump = prev_clump + col;
		    cur_clump[col] = *temp_clump;
		    if (OLD == 0) {
			OLD = *temp_clump;
			}
		    else {
			NEW = *temp_clump;
			if (NEW != OLD) {

			    /* conflict! preserve NEW clump ID and change OLD clump ID.
			     * Must go back to the left in the current row and to the right
			     * in the previous row to change all the clump values as well.
			     */

			    /* left of the current row from 1 to col - 1 */
			    temp_clump = cur_clump;
			    n = col - 1;
			    while (n-- > 0) {
				temp_clump++;	/* skip left edge */
				if (*temp_clump == OLD)
				    *temp_clump = NEW;
			    }

			    /* right of previous row from col + 1 to ncols */
			    temp_clump = prev_clump + col;
			    n = ncols - col;
			    while (n-- > 0) {
				temp_clump++;	/* skip col */
				if (*temp_clump == OLD)
				    *temp_clump = NEW;
			    }

			    /* modify the OLD index */
			    index[OLD] = NEW;

			    OLD = NEW;
			    NEW = 0;
			}
		    }
		}
	    }

	    if (NEW == 0 || OLD == NEW) {	/* ok */
		if (OLD == 0) {
		    /* start a new clump */
		    if (label == cellmax)
			G_fatal_error(_("Too many objects: integer overflow"));

		    label++;
		    cur_clump[col] = label;
		    if (label >= nalloc) {
			nalloc += INCR;
			index =
			    (CELL *) G_realloc(index,
					       nalloc * sizeof(CELL));
		    }
		    index[label] = label;
		}
	    }
	    /* else the relabelling above failed */
	}

	/* write initial clump IDs */
	/* this works also with writing out cur_clump, but only 
	 * prev_clump is complete and will not change any more */
	if (row > 0) {
	    if (write(cfd, prev_clump + 1, csize) != csize)
		G_fatal_error(_("Unable to write to temp file"));
	}

	/* switch the buffers so that the current buffer becomes the previous */
	temp_in = cur_in;
	cur_in = prev_in;
	prev_in = temp_in;

	temp_clump = cur_clump;
	cur_clump = prev_clump;
	prev_clump = temp_clump;
    }
    /* write last row with initial clump IDs */
    if (write(cfd, prev_clump + 1, csize) != csize)
	G_fatal_error(_("Unable to write to temp file"));
    G_percent(1, 1, 1);

    G_free(prev_in[0]);
    G_free(cur_in[0]);
    G_free(prev_in);
    G_free(cur_in);

    /* generate a renumbering scheme */
    G_message(_("Generating renumbering scheme..."));
    G_debug(1, "%d initial labels", label);
    /* allocate final clump ID */
    renumber = (CELL *) G_malloc((label + 1) * sizeof(CELL));
    renumber[0] = 0;
    cat = 0;
    G_percent(0, label, 1);
    for (n = 1; n <= label; n++) {
	G_percent(n, label, 1);
	OLD = n;
	NEW = index[n];
	if (OLD != NEW) {
	    renumber[n] = 0;
	    /* find valid clump ID */
	    while (OLD != NEW) {
		OLD = NEW;
		NEW = index[OLD];
	    }
	    index[n] = NEW;
	}
	else
	    /* set final clump id */
	    renumber[n] = ++cat;
    }

    if (cat > cellmax - globals->max_rid)
	G_fatal_error(_("Too many objects: integer overflow"));

    /* rewind temp file */
    lseek(cfd, 0, SEEK_SET);

    /****************************************************
     *                      PASS 2                      *
     * apply renumbering scheme to initial clump labels *
     ****************************************************/

    /* the input raster is no longer needed, 
     * using instead the temp file with initial clump labels */

    G_message(_("Assigning final region IDs..."));
    for (row = 0; row < nrows; row++) {

	G_percent(row, nrows, 2);
    
	if (read(cfd, cur_clump, csize) != csize)
	    G_fatal_error(_("Unable to read from temp file"));

	temp_clump = cur_clump;

	for (col = 0; col < ncols; col++) {
	    if (!(FLAG_GET(globals->null_flag, row, col))) {
		out_cell = renumber[index[*temp_clump]] + globals->max_rid;

		Segment_put(&globals->rid_seg, (void *)&out_cell, row, col);
	    }
	    temp_clump++;
	}
    }
    G_percent(1, 1, 1);

    close(cfd);
    unlink(cname);

    /* free */
    G_free(prev_clump);
    G_free(cur_clump);
    G_free(index);
    G_free(renumber);

    G_message(_("Found %d clumps"), cat);
    globals->max_rid += cat;

    return globals->max_rid;
}