source: grass/trunk/raster/r.volume/centroids.c

/* Find 'centroids' of all cats.  Most useful with cats which are
 * like clumps (contiguous areas), but will work on any file.
 * Respects the current window and mask. 
 * Category zero and negative categories are not done!
 * Returned centroids are (row, col) pairs in n, e.
 * Returned value is 0 for most cases, > 0 if 'both' method
 * selected and some values fell outside 'clumps', and were
 * adjusted.
 * Two methods can be used: 'distance-weighted' and 'counting'.
 * 
 * Method 0 = 'counting' or 'clump'; centroid point guaranteed
 * to be at a cell of the given category.
 * 1 = Both 0 and 2 are run; if method 2 centroid is on
 * cell of proper cat, it is used, otherwise
 * method 1 value is substituted.
 * 2 = 'distance-weighted'; row = sigma(row)/n,
 * col = sigma(col)/n.
 */
#include <grass/gis.h>
#include <grass/raster.h>

int centroids(int fd,           /* File descriptor of map layer to process */
              /* This file is assumed to be opened before calling */
              /*   centroids. */
              unsigned long *e, unsigned long *n,       /* Pointers to arrays at least max+1 long */
              int method,       /* 0, 1, or 2; see above. */
              int max)
{                               /* Highest positive cat number in map layer */
    CELL *cell_buf, v;
    int i, adjusted, numb, left, right;
    long int *count;
    int row, col, rows, cols;

    adjusted = 0;

    cell_buf = Rast_allocate_c_buf();
    /* space to accumulate counts */
    count = (long int *)G_malloc((max + 1) * sizeof(long int));

    /* zero the count totals */
    for (i = 1; i <= max; i++) {
        count[i] = 0;
        e[i] = 0;
        n[i] = 0;
    }

    /* do rows and columns through window and mask */
    /*  to do counting */
    rows = Rast_window_rows();
    cols = Rast_window_cols();
    for (row = 0; row < rows; row++) {
        Rast_get_c_row(fd, cell_buf, row);      /* get a row */
        for (col = 0; col < cols; col++) {
            v = cell_buf[col];  /* next cell value in row */
            if (v < 1)
                continue;       /* can't handle 0 or - values */
            count[v]++;
            if (method > 0) {   /* acccumulate row, col weights */
                e[v] += col;
                n[v] += row;
            }
        }
    }

    /* compute averages */
    if (method > 0)
        for (i = 1; i <= max; i++) {
            if (count[i]) {
                numb = count[i];
                e[i] /= numb;
                n[i] /= numb;
            }
        }

    /* substitute 'count' centroids for 'weighted' ones, if necessary */
    if (method == 1) {
        for (i = 1; i <= max; i++) {
            if (count[i]) {
                row = n[i];
                col = e[i];
                /* get cell at row,col */
                Rast_get_c_row(fd, cell_buf, row);
                v = cell_buf[col];
                if (v > 0) {
                    if (v == i)
                        count[i] = 0;   /* weighted is acceptable */
                    else
                        adjusted++;
                }
            }
        }
    }

    /* compute middle cell; zeros will remain zeros */
    for (i = 1; i <= max; i++)
        count[i] = (count[i] + 1) / 2;

    /* go through map again */
    for (row = 0; row < rows; row++) {
        Rast_get_c_row(fd, cell_buf, row);
        for (col = 0; col < cols; col++) {
            v = cell_buf[col];
            if (v < 1)
                continue;
            if (count[v] == 0)
                continue;
            if ((--count[v]) == 0) {    /* then this is middle cell */
                n[v] = row;
                /* find row-center in this clump */
                left = right = col;
                /* left edge first */
                while (left > 0)
                    if (cell_buf[--left] != v) {
                        left++;
                        break;
                    }
                /* then right edge */
                while (right < cols)
                    if (cell_buf[++right] != v) {
                        right--;
                        break;
                    }
                e[v] = (left + right) / 2;
            }
        }
    }
    G_free(cell_buf);
    G_free(count);
    return (adjusted);
}                               /* end of centroids() */