Add null-ignoring nmin, nmax, nmedian, nmode functions

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@52178 15284696-431f-4ddb-bdfa-cd5b030d7da7

raster/r.mapcalc/func_proto.h

@@ -78,12 +78,16 @@ extern args_t c_graph;
 
 extern func_t f_min;
 extern func_t f_max;
+extern func_t f_nmin;
+extern func_t f_nmax;
 extern args_t c_varop;
 
 extern func_t f_median;
+extern func_t f_nmedian;
 extern args_t c_median;
 
 extern func_t f_mode;
+extern func_t f_nmode;
 extern args_t c_mode;
 
 extern func_t f_rand;

raster/r.mapcalc/function.c

@@ -68,6 +68,11 @@ func_desc func_descs[] = {
     {"median", c_varop, f_median},
     {"mode", c_varop, f_mode},
 
+    {"nmax", c_varop, f_nmax},
+    {"nmin", c_varop, f_nmin},
+    {"nmedian", c_varop, f_nmedian},
+    {"nmode", c_varop, f_nmode},
+
     {"graph", c_graph, f_graph},
 
     {"rand", c_binop, f_rand},

raster/r.mapcalc/r.mapcalc.html

@@ -311,6 +311,10 @@ max(x,y[,z...])         largest value of those listed                   *
 median(x,y[,z...])      median value of those listed                    *
 min(x,y[,z...])         smallest value of those listed                  *
 mode(x,y[,z...])        mode value of those listed                      *
+nmax(x,y[,z...])        largest value of those listed, excluding NULLs  *
+nmedian(x,y[,z...])     median value of those listed, excluding NULLs   *
+nmin(x,y[,z...])        smallest value of those listed, excluding NULLs *
+nmode(x,y[,z...])       mode value of those listed, excluding NULLs     *
 not(x)                  1 if x is zero, 0 otherwise
 pow(x,y)                x to the power y                                *
 rand(a,b)               random value x : a <= x < b

raster/r.mapcalc/r3.mapcalc.html

@@ -215,6 +215,10 @@ max(x,y[,z...])         largest value of those listed                   *
 median(x,y[,z...])      median value of those listed                    *
 min(x,y[,z...])         smallest value of those listed                  *
 mode(x,y[,z...])        mode value of those listed                      *
+nmax(x,y[,z...])        largest value of those listed, excluding NULLs  *
+nmedian(x,y[,z...])     median value of those listed, excluding NULLs   *
+nmin(x,y[,z...])        smallest value of those listed, excluding NULLs *
+nmode(x,y[,z...])       mode value of those listed, excluding NULLs     *
 not(x)                  1 if x is zero, 0 otherwise
 pow(x,y)                x to the power y                                *
 rand(a,b)               random value x : a <= x < b

raster/r.mapcalc/xnmax.c

@@ -0,0 +1,100 @@
+
+#include <stdlib.h>
+
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include "globals.h"
+#include "expression.h"
+#include "func_proto.h"
+
+/****************************************************************
+max(x0,x1,...,xn) returns maximum value
+****************************************************************/
+
+int f_nmax(int argc, const int *argt, void **args)
+{
+    int i, j;
+
+    if (argc < 1)
+	return E_ARG_LO;
+
+    for (i = 1; i <= argc; i++)
+	if (argt[i] != argt[0])
+	    return E_ARG_TYPE;
+
+    switch (argt[0]) {
+    case CELL_TYPE:
+	{
+	    CELL *res = args[0];
+	    CELL **argz = (CELL **) args;
+
+	    for (i = 0; i < columns; i++) {
+		int nul = 1;
+		CELL max;
+
+		for (j = 1; j <= argc; j++)
+		    if (IS_NULL_C(&argz[j][i]))
+			continue;
+		    else if (nul)
+			max = argz[j][i], nul = 0;
+		    else if (max < argz[j][i])
+			max = argz[j][i], nul = 0;
+		if (nul)
+		    SET_NULL_C(&res[i]);
+		else
+		    res[i] = max;
+	    }
+	    return 0;
+	}
+    case FCELL_TYPE:
+	{
+	    FCELL *res = args[0];
+	    FCELL **argz = (FCELL **) args;
+
+	    for (i = 0; i < columns; i++) {
+		int nul = 1;
+		FCELL max;
+
+		for (j = 1; j <= argc; j++)
+		    if (IS_NULL_F(&argz[j][i]))
+			continue;
+		    else if (nul)
+			max = argz[j][i], nul = 0;
+		    else if (max < argz[j][i])
+			max = argz[j][i], nul = 0;
+		if (nul)
+		    SET_NULL_F(&res[i]);
+		else
+		    res[i] = max;
+	    }
+
+	    return 0;
+	}
+    case DCELL_TYPE:
+	{
+	    DCELL *res = args[0];
+	    DCELL **argz = (DCELL **) args;
+
+	    for (i = 0; i < columns; i++) {
+		int nul = 1;
+		DCELL max;
+
+		for (j = 1; j <= argc; j++)
+		    if (IS_NULL_D(&argz[j][i]))
+			continue;
+		    else if (nul)
+			max = argz[j][i], nul = 0;
+		    else if (max < argz[j][i])
+			max = argz[j][i], nul = 0;
+		if (nul)
+		    SET_NULL_D(&res[i]);
+		else
+		    res[i] = max;
+	    }
+
+	    return 0;
+	}
+    default:
+	return E_INV_TYPE;
+    }
+}

raster/r.mapcalc/xnmedian.c

@@ -0,0 +1,155 @@
+
+#include <stdlib.h>
+
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/raster.h>
+#include "globals.h"
+#include "expression.h"
+#include "func_proto.h"
+
+/**********************************************************************
+median(x1,x2,..,xn)
+   return median of arguments
+**********************************************************************/
+
+static int icmp(const void *aa, const void *bb)
+{
+    const CELL *a = aa;
+    const CELL *b = bb;
+
+    return *a - *b;
+}
+
+static int fcmp(const void *aa, const void *bb)
+{
+    const FCELL *a = aa;
+    const FCELL *b = bb;
+
+    if (*a < *b)
+	return -1;
+    if (*a > *b)
+	return 1;
+    return 0;
+}
+
+static int dcmp(const void *aa, const void *bb)
+{
+    const DCELL *a = aa;
+    const DCELL *b = bb;
+
+    if (*a < *b)
+	return -1;
+    if (*a > *b)
+	return 1;
+    return 0;
+}
+
+int f_nmedian(int argc, const int *argt, void **args)
+{
+    static void *array;
+    static int alloc;
+    int size = argc * Rast_cell_size(argt[0]);
+    int i, j;
+
+    if (argc < 1)
+	return E_ARG_LO;
+
+    for (i = 1; i <= argc; i++)
+	if (argt[i] != argt[0])
+	    return E_ARG_TYPE;
+
+    if (size > alloc) {
+	alloc = size;
+	array = G_realloc(array, size);
+    }
+
+    switch (argt[0]) {
+    case CELL_TYPE:
+	{
+	    CELL *res = args[0];
+	    CELL **argv = (CELL **) &args[1];
+	    CELL *a = array;
+	    CELL *a1 = &a[(argc - 1) / 2];
+	    CELL *a2 = &a[argc / 2];
+
+	    for (i = 0; i < columns; i++) {
+		int nv = 1;
+
+		for (j = 0; j < argc; j++) {
+		    if (IS_NULL_C(&argv[j][i]))
+			continue;
+		    a[j] = argv[j][i];
+		    nv = 0;
+		}
+
+		if (nv)
+		    SET_NULL_C(&res[i]);
+		else {
+		    qsort(a, argc, sizeof(CELL), icmp);
+		    res[i] = (*a1 + *a2) / 2;
+		}
+	    }
+
+	    return 0;
+	}
+    case FCELL_TYPE:
+	{
+	    FCELL *res = args[0];
+	    FCELL **argv = (FCELL **) &args[1];
+	    FCELL *a = array;
+	    FCELL *a1 = &a[(argc - 1) / 2];
+	    FCELL *a2 = &a[argc / 2];
+
+	    for (i = 0; i < columns; i++) {
+		int nv = 1;
+
+		for (j = 0; j < argc; j++) {
+		    if (IS_NULL_F(&argv[j][i]))
+			continue;
+		    a[j] = argv[j][i];
+		    nv = 0;
+		}
+
+		if (nv)
+		    SET_NULL_F(&res[i]);
+		else {
+		    qsort(a, argc, sizeof(FCELL), fcmp);
+		    res[i] = (*a1 + *a2) / 2;
+		}
+	    }
+
+	    return 0;
+	}
+    case DCELL_TYPE:
+	{
+	    DCELL *res = args[0];
+	    DCELL **argv = (DCELL **) &args[1];
+	    DCELL *a = array;
+	    DCELL *a1 = &a[(argc - 1) / 2];
+	    DCELL *a2 = &a[argc / 2];
+
+	    for (i = 0; i < columns; i++) {
+		int nv = 1;
+
+		for (j = 0; j < argc; j++) {
+		    if (IS_NULL_D(&argv[j][i]))
+			continue;
+		    a[j] = argv[j][i];
+		    nv = 0;
+		}
+
+		if (nv)
+		    SET_NULL_D(&res[i]);
+		else {
+		    qsort(a, argc, sizeof(DCELL), dcmp);
+		    res[i] = (*a1 + *a2) / 2;
+		}
+	    }
+
+	    return 0;
+	}
+    default:
+	return E_INV_TYPE;
+    }
+}

raster/r.mapcalc/xnmin.c

@@ -0,0 +1,100 @@
+
+#include <stdlib.h>
+
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include "globals.h"
+#include "expression.h"
+#include "func_proto.h"
+
+/****************************************************************
+min(x0,x1,...,xn) returns minimum value
+****************************************************************/
+
+int f_nmin(int argc, const int *argt, void **args)
+{
+    int i, j;
+
+    if (argc < 1)
+	return E_ARG_LO;
+
+    for (i = 1; i <= argc; i++)
+	if (argt[i] != argt[0])
+	    return E_ARG_TYPE;
+
+    switch (argt[0]) {
+    case CELL_TYPE:
+	{
+	    CELL *res = args[0];
+	    CELL **argz = (CELL **) args;
+
+	    for (i = 0; i < columns; i++) {
+		int nul = 1;
+		CELL min;
+
+		for (j = 1; j <= argc; j++)
+		    if (IS_NULL_C(&argz[j][i]))
+			continue;
+		    else if (nul)
+			min = argz[j][i], nul = 0;
+		    else if (min > argz[j][i])
+			min = argz[j][i], nul = 0;
+		if (nul)
+		    SET_NULL_C(&res[i]);
+		else
+		    res[i] = min;
+	    }
+	    return 0;
+	}
+    case FCELL_TYPE:
+	{
+	    FCELL *res = args[0];
+	    FCELL **argz = (FCELL **) args;
+
+	    for (i = 0; i < columns; i++) {
+		int nul = 1;
+		FCELL min;
+
+		for (j = 1; j <= argc; j++)
+		    if (IS_NULL_F(&argz[j][i]))
+			continue;
+		    else if (nul)
+			min = argz[j][i], nul = 0;
+		    else if (min > argz[j][i])
+			min = argz[j][i], nul = 0;
+		if (nul)
+		    SET_NULL_F(&res[i]);
+		else
+		    res[i] = min;
+	    }
+
+	    return 0;
+	}
+    case DCELL_TYPE:
+	{
+	    DCELL *res = args[0];
+	    DCELL **argz = (DCELL **) args;
+
+	    for (i = 0; i < columns; i++) {
+		int nul = 1;
+		DCELL min;
+
+		for (j = 1; j <= argc; j++)
+		    if (IS_NULL_D(&argz[j][i]))
+			continue;
+		    else if (nul)
+			min = argz[j][i], nul = 0;
+		    else if (min > argz[j][i])
+			min = argz[j][i], nul = 0;
+		if (nul)
+		    SET_NULL_D(&res[i]);
+		else
+		    res[i] = min;
+	    }
+
+	    return 0;
+	}
+    default:
+	return E_INV_TYPE;
+    }
+}

raster/r.mapcalc/xnmode.c

@@ -0,0 +1,144 @@
+
+#include <stdlib.h>
+
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include "globals.h"
+#include "expression.h"
+#include "func_proto.h"
+
+/**********************************************************************
+mode(x1,x2,..,xn)
+   return mode of arguments
+**********************************************************************/
+
+static int dcmp(const void *aa, const void *bb)
+{
+    const double *a = aa;
+    const double *b = bb;
+
+    if (*a < *b)
+	return -1;
+    if (*a > *b)
+	return 1;
+    return 0;
+}
+
+static double mode(double *value, int argc)
+{
+    double mode_v;
+    int mode_n = 0;
+    int i;
+
+    qsort(value, argc, sizeof(double), dcmp);
+
+    for (i = 0; i < argc;) {
+	int n = 1;
+	double v = value[i];
+
+	for (i++; i < argc; i++) {
+	    if (value[i] != v)
+		break;
+	    n++;
+	}
+
+	if (n < mode_n)
+	    continue;
+
+	mode_v = v;
+	mode_n = n;
+    }
+
+    return mode_v;
+}
+
+int f_nmode(int argc, const int *argt, void **args)
+{
+    static double *value;
+    static int value_size;
+    int size = argc * sizeof(double);
+    int i, j;
+
+    if (argc < 1)
+	return E_ARG_LO;
+
+    for (i = 1; i <= argc; i++)
+	if (argt[i] != argt[0])
+	    return E_ARG_TYPE;
+
+    if (size > value_size) {
+	value_size = size;
+	value = G_realloc(value, value_size);
+    }
+
+    switch (argt[argc]) {
+    case CELL_TYPE:
+	{
+	    CELL *res = args[0];
+	    CELL **argv = (CELL **) & args[1];
+
+	    for (i = 0; i < columns; i++) {
+		int nv = 1;
+
+		for (j = 0; j < argc; j++) {
+		    if (IS_NULL_C(&argv[j][i]))
+			continue;
+		    value[j] = (double)argv[j][i];
+		    nv = 0;
+		}
+
+		if (nv)
+		    SET_NULL_C(&res[i]);
+		else
+		    res[i] = (CELL) mode(value, argc);
+	    }
+	    return 0;
+	}
+    case FCELL_TYPE:
+	{
+	    FCELL *res = args[0];
+	    FCELL **argv = (FCELL **) & args[1];
+
+	    for (i = 0; i < columns; i++) {
+		int nv = 1;
+
+		for (j = 0; j < argc; j++) {
+		    if (IS_NULL_F(&argv[j][i]))
+			continue;
+		    value[j] = (double)argv[j][i];
+		    nv = 0;
+		}
+
+		if (nv)
+		    SET_NULL_F(&res[i]);
+		else
+		    res[i] = (FCELL) mode(value, argc);
+	    }
+	    return 0;
+	}
+    case DCELL_TYPE:
+	{
+	    DCELL *res = args[0];
+	    DCELL **argv = (DCELL **) & args[1];
+
+	    for (i = 0; i < columns; i++) {
+		int nv = 1;
+
+		for (j = 0; j < argc; j++) {
+		    if (IS_NULL_D(&argv[j][i]))
+			continue;
+		    value[j] = (double)argv[j][i];
+		    nv = 0;
+		}
+
+		if (nv)
+		    SET_NULL_D(&res[i]);
+		else
+		    res[i] = (DCELL) mode(value, argc);
+	    }
+	    return 0;
+	}
+    default:
+	return E_INV_TYPE;
+    }
+}