Ticket #902: st_minmax.2.patch

File st_minmax.2.patch, 20.7 KB (added by dustymugs, 5 years ago)

New patch due to refactoring in preparation for computing additonal stats

  • raster/rt_core/rt_api.c

    diff -rupN postgis-old/raster/rt_core/rt_api.c postgis-new/raster/rt_core/rt_api.c
    old new rt_raster_from_band(rt_context ctx, rt_r 
    36713671                rt_raster_get_num_bands(ctx, rast));
    36723672        return rast;
    36733673}
     3674
     3675struct rt_bandstats_t {
     3676        uint32_t count;
     3677        double min;
     3678        double max;
     3679        double mean;
     3680        double stddev;
     3681
     3682        double *values;
     3683};
     3684
     3685/**
     3686 * Compute summary statistics for a band
     3687 * @param ctx: context, for thread safety
     3688 * @param band: the band to query for minimum and maximum pixel values
     3689 * @param no_nodata: if non-zero, ignore nodata value
     3690 * @param sample: percentage of pixels to sample
     3691 * @param inc_vals: flag to include values in return struct
     3692 * @return the summary statistics for a band
     3693 */
     3694rt_bandstats
     3695rt_band_get_summary_stats(rt_context ctx, rt_band band, int no_nodata,
     3696        double sample, int inc_vals) {
     3697        rt_pixtype pixtype = PT_END;
     3698        uint8_t *data = NULL;
     3699        uint32_t x = 0;
     3700        uint32_t y = 0;
     3701        uint32_t z = 0;
     3702        uint32_t offset = 0;
     3703        uint32_t diff = 0;
     3704        int rtn;
     3705        int hasnodata = FALSE;
     3706        double nodata = 0;
     3707        double *values = NULL;
     3708        double value;
     3709        rt_bandstats stats = NULL;
     3710
     3711        uint32_t sample_size = 0;
     3712        int byY = 0;
     3713        uint32_t outer = 0;
     3714        uint32_t inner = 0;
     3715        uint32_t sample_per = 0;
     3716        uint32_t sample_int = 0;
     3717        uint32_t i = 0;
     3718        uint32_t j = 0;
     3719        double sum = 0;
     3720        uint32_t k = 0;
     3721        double M = 0;
     3722        double Q = 0;
     3723
     3724        clock_t start, stop;
     3725        double elapsed = 0;
     3726
     3727        RASTER_DEBUG(3, "starting");
     3728        start = clock();
     3729
     3730        assert(NULL != ctx);
     3731        assert(NULL != band);
     3732
     3733        if (band->offline) {
     3734                ctx->err("rt_band_get_summary_stats not implemented yet for OFFDB bands");
     3735                return NULL;
     3736        }
     3737
     3738        data = rt_band_get_data(ctx, band);
     3739        pixtype = band->pixtype;
     3740
     3741        hasnodata = rt_band_get_hasnodata_flag(ctx, band);
     3742        if (hasnodata != FALSE) {
     3743                nodata = rt_band_get_nodata(ctx, band);
     3744        }
     3745        else {
     3746                no_nodata = 0;
     3747        }
     3748
     3749        RASTER_DEBUGF(3, "nodata = %f", nodata);
     3750        RASTER_DEBUGF(3, "hasnodata = %d", hasnodata);
     3751        RASTER_DEBUGF(3, "no_nodata = %d", no_nodata);
     3752
     3753        /* entire band is nodata */
     3754        if (rt_band_get_isnodata_flag(ctx, band) != FALSE) {
     3755                if (no_nodata) {
     3756                        ctx->warn("All pixels of band have the NODATA value");
     3757                        return NULL;
     3758                }
     3759                else {
     3760                        stats = (rt_bandstats) ctx->alloc(sizeof(struct rt_bandstats_t));
     3761                        stats[0].min = stats[0].max = nodata;
     3762                        stats[0].mean = nodata;
     3763                        stats[0].stddev = 0;
     3764                        stats[0].count = band->width * band->height;
     3765                        stats[0].values = NULL;
     3766                        return stats;
     3767                }
     3768        }
     3769
     3770        if (band->height > band->width) {
     3771                byY = 1;
     3772                outer = band->height;
     3773                inner = band->width;
     3774        }
     3775        else {
     3776                byY = 0;
     3777                outer = band->width;
     3778                inner = band->height;
     3779        }
     3780
     3781        /* clamp percentage */
     3782        if (sample < 0)
     3783                sample = 0;
     3784        else if (sample > 100)
     3785                sample = 100;
     3786
     3787        /* sample all pixels */
     3788        if (sample == 0 || sample == 100) {
     3789                sample_size = band->width * band->height;
     3790                sample_per = inner;
     3791        }
     3792        /*
     3793         randomly sample a percentage of available pixels
     3794         sampling method is known as
     3795                "systematic random sample without replacement"
     3796        */
     3797        else {
     3798                sample_size = round((band->width * band->height) * (sample / 100));
     3799                sample_per = round(sample_size / outer);
     3800                sample_int = round(inner / sample_per);
     3801                srand(time(NULL));
     3802        }
     3803
     3804        RASTER_DEBUGF(3, "sampling %d of %d available pixels w/ %d per set"
     3805                , sample_size, (band->width * band->height), sample_per);
     3806
     3807        if (inc_vals) {
     3808                values = ctx->alloc(sizeof(double) * sample_size);
     3809                if (NULL == values) {
     3810                        ctx->warn("Unable to allocate memory for values");
     3811                        inc_vals = 0;
     3812                }
     3813        }
     3814
     3815        for (x = 0, j = 0, k = 0; x < outer; x++) {
     3816                y = -1;
     3817
     3818                for (i = 0, z = 0; i < sample_per; i++) {
     3819                        if (sample == 0 || sample == 100)
     3820                                y = i;
     3821                        else {
     3822                                offset = (rand() % sample_int) + 1;
     3823                                y += diff + offset;
     3824                                diff = sample_int - offset;
     3825                        }
     3826                        if (y >= inner || z > sample_per) break;
     3827
     3828                        if (byY)
     3829                                rtn = rt_band_get_pixel(ctx, band, y, x, &value);
     3830                        else
     3831                                rtn = rt_band_get_pixel(ctx, band, x, y, &value);
     3832
     3833                        j++;
     3834                        if (rtn != -1) {
     3835                                if (
     3836                                        !no_nodata ||
     3837                                        (no_nodata && (hasnodata != FALSE) && (value != nodata))
     3838                                ) {
     3839
     3840                                        /* inc_vals set, collect pixel values */
     3841                                        if (inc_vals) {
     3842                                                values[k] = value;
     3843                                        }
     3844
     3845                                        /* average */
     3846                                        k++;
     3847                                        sum += value;
     3848
     3849                                        /*
     3850                                                one-pass standard deviation
     3851                                                http://www.eecs.berkeley.edu/~mhoemmen/cs194/Tutorials/variance.pdf
     3852                                        */
     3853                                        if (k == 1) {
     3854                                                Q = 0;
     3855                                                M = value;
     3856                                        }
     3857                                        else {
     3858                                                Q = Q + (((k  - 1) * pow(value - M, 2)) / k);
     3859                                                M = M + ((value - M ) / k);
     3860                                        }
     3861
     3862                                        /* min/max */
     3863                                        if (NULL == stats) {
     3864                                                stats = (rt_bandstats) ctx->alloc(sizeof(struct rt_bandstats_t));
     3865                                                if (NULL == stats) {
     3866                                                }
     3867
     3868                                                stats[0].count = 0;
     3869                                                stats[0].mean = 0;
     3870                                                stats[0].stddev = -1;
     3871                                                stats[0].values = NULL;
     3872
     3873                                                stats[0].min = stats[0].max = value;
     3874                                        }
     3875                                        else {
     3876                                                if (value < stats[0].min)
     3877                                                        stats[0].min = value;
     3878                                                if (value > stats[0].max)
     3879                                                        stats[0].max = value;
     3880                                        }
     3881
     3882                                }
     3883                        }
     3884
     3885                        z++;
     3886                }
     3887        }
     3888
     3889        if (k) {
     3890                if (inc_vals) {
     3891                        /* free unused memory */
     3892                        if (sample_size != k) {
     3893                                ctx->realloc(values, k * sizeof(double));
     3894                        }
     3895
     3896                        stats[0].values = values;
     3897                }
     3898
     3899                stats[0].count = k;
     3900                stats[0].mean = sum / k;
     3901
     3902                /* standard deviation */
     3903                if (sample == 0 || sample == 100)
     3904                        stats[0].stddev = sqrt(Q / k);
     3905                /* sample deviation */
     3906                else {
     3907                        if (k < 2)
     3908                                stats[0].stddev = -1;
     3909                        else
     3910                                stats[0].stddev = sqrt(Q / (k - 1));
     3911                }
     3912        }
     3913        else if (inc_vals) {
     3914                ctx->dealloc(values);
     3915        }
     3916        values = NULL;
     3917
     3918        stop = clock();
     3919        elapsed = ((double) (stop - start)) / CLOCKS_PER_SEC;
     3920        RASTER_DEBUGF(3, "(time, mean, stddev, min, max) = (%0.4f, %f, %f, %f, %f)",
     3921                elapsed, stats[0].mean, stats[0].stddev, stats[0].min, stats[0].max);
     3922        RASTER_DEBUG(3, "done");
     3923        return stats;
     3924}
  • raster/rt_core/rt_api.h

    diff -rupN postgis-old/raster/rt_core/rt_api.h postgis-new/raster/rt_core/rt_api.h
    old new int32_t rt_raster_copy_band(rt_context c 
    752752rt_raster rt_raster_from_band(rt_context ctx, rt_raster raster,
    753753        uint32_t *bandNums, int count);
    754754
     755/**
     756 * Compute summary statistics for a band
     757 * @param ctx: context, for thread safety
     758 * @param band: the band to query for minimum and maximum pixel values
     759 * @param no_nodata: if non-zero, ignore nodata value
     760 * @param sample: percentage of pixels to sample
     761 * @param inc_vals: flag to include values in return struct
     762 * @return the summary statistics for a band
     763 */
     764typedef struct rt_bandstats_t* rt_bandstats;
     765rt_bandstats rt_band_get_summary_stats(rt_context ctx, rt_band band,
     766        int no_nodata, double sample, int inc_vals);
     767
    755768/*- utilities -------------------------------------------------------*/
    756769
    757770/* Set of functions to clamp double to int of different size
  • raster/rt_pg/rt_pg.c

    diff -rupN postgis-old/raster/rt_pg/rt_pg.c postgis-new/raster/rt_pg/rt_pg.c
    old new Datum RASTER_mapAlgebra(PG_FUNCTION_ARGS 
    165165/* create new raster from existing raster's bands */
    166166Datum RASTER_band(PG_FUNCTION_ARGS);
    167167
     168/* Get minimum and maximum pixel values */
     169Datum RASTER_getMinMaxPixelValues(PG_FUNCTION_ARGS);
     170
    168171
    169172/* Replace function taken from http://ubuntuforums.org/showthread.php?s=aa6f015109fd7e4c7e30d2fd8b717497&t=141670&page=3 */
    170173/* ---------------------------------------------------------------------------
    Datum RASTER_band(PG_FUNCTION_ARGS) 
    27982801        PG_RETURN_POINTER(pgrast);
    27992802}
    28002803
     2804struct rt_bandstats_t {
     2805        uint32_t count;
     2806        double min;
     2807        double max;
     2808        double mean;
     2809        double stddev;
     2810
     2811        double *values;
     2812};
     2813
     2814/**
     2815 * Get minimum and maximum pixel values
     2816 */
     2817PG_FUNCTION_INFO_V1(RASTER_getMinMaxPixelValues);
     2818Datum RASTER_getMinMaxPixelValues(PG_FUNCTION_ARGS)
     2819{
     2820        FuncCallContext *funcctx;
     2821        TupleDesc tupdesc;
     2822        AttInMetadata *attinmeta;
     2823
     2824        rt_context ctx;
     2825        rt_bandstats stats;
     2826        rt_bandstats stats2;
     2827        int call_cntr;
     2828        int max_calls;
     2829
     2830        /* first call of function */
     2831        if (SRF_IS_FIRSTCALL()) {
     2832                MemoryContext oldcontext;
     2833
     2834                rt_pgraster *pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
     2835                rt_raster raster = NULL;
     2836                rt_band band = NULL;
     2837                int32_t bandindex = 0;
     2838                bool no_nodata = TRUE;
     2839                int num_bands = 0;
     2840                double sample = 0;
     2841
     2842                /* create a function context for cross-call persistence */
     2843                funcctx = SRF_FIRSTCALL_INIT();
     2844
     2845                /* switch to memory context appropriate for multiple function calls */
     2846                oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
     2847
     2848                /**
     2849                 * Create a new context. We don't call get_rt_context, because this
     2850                 * function changes the memory context to the one pointed by
     2851                 * fcinfo->flinfo->fn_mcxt, and we want to keep the current one
     2852                 */
     2853                ctx = rt_context_new(rt_pgalloc, rt_pgrealloc, rt_pgfree);
     2854                rt_context_set_message_handlers(ctx, rt_pgerr, rt_pgwarn, rt_pginfo);
     2855
     2856                raster = rt_raster_deserialize(ctx, pgraster);
     2857                if (!raster) {
     2858                        elog(ERROR, "RASTER_getMinMaxPixelValues: Could not deserialize raster");
     2859                        rt_context_destroy(ctx);       
     2860                        PG_RETURN_NULL();
     2861                }
     2862
     2863                /* band index is 1-based */
     2864                bandindex = PG_GETARG_INT32(1);
     2865                num_bands = rt_raster_get_num_bands(ctx, raster);
     2866                if (bandindex < 1) {
     2867                        elog(ERROR, "RASTER_getMinMaxPixelValues: Band index less than 1");
     2868                        rt_raster_destroy(ctx, raster);
     2869                        rt_context_destroy(ctx);       
     2870                        PG_RETURN_NULL();
     2871                }
     2872                else if (bandindex > num_bands) {
     2873                        elog(ERROR, "RASTER_getMinMaxPixelValues: Band index more than maximum band index");
     2874                        rt_raster_destroy(ctx, raster);
     2875                        rt_context_destroy(ctx);       
     2876                        PG_RETURN_NULL();
     2877                }
     2878                assert(0 <= (bandindex - 1));
     2879
     2880                /* no_nodata flag */
     2881                no_nodata = PG_GETARG_BOOL(2);
     2882
     2883                /* sample % */
     2884                if (!PG_ARGISNULL(3))
     2885                        sample = PG_GETARG_FLOAT8(3);
     2886                else
     2887                        sample = 100;
     2888
     2889                /* get band */
     2890                band = rt_raster_get_band(ctx, raster, bandindex - 1);
     2891                if (!band) {
     2892                        elog(ERROR, "RASTER_getMinMaxPixelValues: Could not retrieve band at index %d", bandindex);
     2893                        rt_raster_destroy(ctx, raster);
     2894                        rt_context_destroy(ctx);       
     2895                        PG_RETURN_NULL();
     2896                }
     2897
     2898                /* we don't need the raw values, hence the zero parameter */
     2899                stats = rt_band_get_summary_stats(ctx, band, (int) no_nodata, sample, 0);
     2900                rt_band_destroy(ctx, band);
     2901                rt_raster_destroy(ctx, raster);
     2902                if (NULL == stats) {
     2903                        elog(ERROR, "RASTER_getMinMaxPixelValues: Could not retrieve minimum and maximum values of band at index %d", bandindex);
     2904                        PG_RETURN_NULL();
     2905                }
     2906
     2907                /* store needed data */
     2908                funcctx->user_fctx = stats;
     2909
     2910                /* total number of tuples to return */
     2911                funcctx->max_calls = 1;
     2912
     2913                /* Build a tuple descriptor for our result type */
     2914                if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
     2915                        ereport(ERROR, (
     2916                                errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
     2917                                errmsg(
     2918                                        "function returning record called in context "
     2919                                        "that cannot accept type record"
     2920                                )
     2921                        ));
     2922                }
     2923
     2924                /*
     2925                 * generate attribute metadata needed later to produce tuples from raw
     2926                 * C strings
     2927                 */
     2928                attinmeta = TupleDescGetAttInMetadata(tupdesc);
     2929                funcctx->attinmeta = attinmeta;
     2930                MemoryContextSwitchTo(oldcontext);
     2931
     2932        }
     2933
     2934  /* stuff done on every call of the function */
     2935  funcctx = SRF_PERCALL_SETUP();
     2936
     2937        call_cntr = funcctx->call_cntr;
     2938        max_calls = funcctx->max_calls;
     2939        attinmeta = funcctx->attinmeta;
     2940        stats2 = funcctx->user_fctx;
     2941
     2942        /* do when there is more left to send */
     2943        if (call_cntr < max_calls) {
     2944                char **values;
     2945                HeapTuple tuple;
     2946                Datum result;
     2947
     2948                POSTGIS_RT_DEBUGF(3, "tuple #%d", call_cntr);
     2949                POSTGIS_RT_DEBUGF(3, "minmax %f, %f", stats2[call_cntr].min, stats2[call_cntr].max);
     2950
     2951                /*
     2952                 * Prepare a values array for building the returned tuple.
     2953                 * This should be an array of C strings which will
     2954                 * be processed later by the type input functions.
     2955                 */
     2956                values = (char **) palloc(2 * sizeof(char *));
     2957
     2958                values[0] = (char *) palloc(
     2959                        (snprintf(NULL, 0, "%f", stats2[call_cntr].min) + 1) * sizeof(char)
     2960                );
     2961                values[1] = (char *) palloc(
     2962                        (snprintf(NULL, 0, "%f", stats2[call_cntr].max) + 1) * sizeof(char)
     2963                );
     2964
     2965                snprintf(
     2966                        values[0],
     2967                        (snprintf(NULL, 0, "%f", stats2[call_cntr].min) + 1) * sizeof(char),
     2968                        "%f",
     2969                        stats2[call_cntr].min
     2970                );
     2971                snprintf(
     2972                        values[1],
     2973                        (snprintf(NULL, 0, "%f", stats2[call_cntr].max) + 1) * sizeof(char),
     2974                        "%f",
     2975                        stats2[call_cntr].max
     2976                );
     2977
     2978                /* build a tuple */
     2979                tuple = BuildTupleFromCStrings(attinmeta, values);
     2980
     2981                /* make the tuple into a datum */
     2982                result = HeapTupleGetDatum(tuple);
     2983
     2984                /* clean up */
     2985                pfree(values[1]);
     2986                pfree(values[0]);
     2987                pfree(values);
     2988
     2989                SRF_RETURN_NEXT(funcctx, result);
     2990        }
     2991        else {
     2992                pfree(stats2);
     2993                SRF_RETURN_DONE(funcctx);
     2994        }
     2995}
     2996
    28012997/* ---------------------------------------------------------------- */
    28022998/*  Memory allocation / error reporting hooks                       */
    28032999/* ---------------------------------------------------------------- */
  • raster/rt_pg/rtpostgis.sql.in.c

    diff -rupN postgis-old/raster/rt_pg/rtpostgis.sql.in.c postgis-new/raster/rt_pg/rtpostgis.sql.in.c
    old new CREATE OR REPLACE FUNCTION st_band(rast 
    257257        LANGUAGE 'SQL' IMMUTABLE STRICT;
    258258
    259259-----------------------------------------------------------------------
     260-- ST_MinMax and ST_ApproxMinMax
     261-----------------------------------------------------------------------
     262CREATE OR REPLACE FUNCTION _st_minmax(rast raster, nband int, ignore_nodata boolean, sample float8, OUT min float8, OUT max float8)
     263        RETURNS record
     264        AS 'MODULE_PATHNAME','RASTER_getMinMaxPixelValues'
     265        LANGUAGE 'C' IMMUTABLE STRICT;
     266
     267CREATE OR REPLACE FUNCTION st_minmax(rast raster, nband int, ignore_nodata boolean, OUT min float8, OUT max float8)
     268        RETURNS record
     269        AS $$ SELECT * FROM _st_minmax($1, $2, $3, 100) $$
     270        LANGUAGE 'SQL' IMMUTABLE STRICT;
     271
     272CREATE OR REPLACE FUNCTION st_minmax(rast raster, nband int, OUT min float8, OUT max float8)
     273        RETURNS record
     274        AS $$ SELECT * FROM _st_minmax($1, $2, TRUE, 100) $$
     275        LANGUAGE 'SQL' IMMUTABLE STRICT;
     276
     277CREATE OR REPLACE FUNCTION st_minmax(rast raster, ignore_nodata boolean, OUT min float8, OUT max float8)
     278        RETURNS record
     279        AS $$ SELECT * FROM _st_minmax($1, 1, $2, 100) $$
     280        LANGUAGE 'SQL' IMMUTABLE STRICT;
     281
     282CREATE OR REPLACE FUNCTION st_minmax(rast raster, OUT min float8, OUT max float8)
     283        RETURNS record
     284        AS $$ SELECT * FROM _st_minmax($1, 1, TRUE, 100) $$
     285        LANGUAGE 'SQL' IMMUTABLE STRICT;
     286
     287CREATE OR REPLACE FUNCTION st_approxminmax(rast raster, nband int, ignore_nodata boolean, sample float8, OUT min float8, OUT max float8)
     288        RETURNS record
     289        AS $$ SELECT * FROM _st_minmax($1, $2, $3, $4) $$
     290        LANGUAGE 'SQL' IMMUTABLE STRICT;
     291
     292CREATE OR REPLACE FUNCTION st_approxminmax(rast raster, nband int, sample float8, OUT min float8, OUT max float8)
     293        RETURNS record
     294        AS $$ SELECT * FROM _st_minmax($1, $2, TRUE, $3) $$
     295        LANGUAGE 'SQL' IMMUTABLE STRICT;
     296
     297CREATE OR REPLACE FUNCTION st_approxminmax(rast raster, ignore_nodata boolean, sample float8, OUT min float8, OUT max float8)
     298        RETURNS record
     299        AS $$ SELECT * FROM _st_minmax($1, 1, $2, $3) $$
     300        LANGUAGE 'SQL' IMMUTABLE STRICT;
     301
     302CREATE OR REPLACE FUNCTION st_approxminmax(rast raster, sample float8, OUT min float8, OUT max float8)
     303        RETURNS record
     304        AS $$ SELECT * FROM _st_minmax($1, 1, TRUE, $2) $$
     305        LANGUAGE 'SQL' IMMUTABLE STRICT;
     306
     307CREATE OR REPLACE FUNCTION st_approxminmax(rast raster, OUT min float8, OUT max float8)
     308        RETURNS record
     309        AS $$ SELECT * FROM _st_minmax($1, 1, TRUE, 10) $$
     310        LANGUAGE 'SQL' IMMUTABLE STRICT;
     311
     312-----------------------------------------------------------------------
    260313-- MapAlgebra
    261314-----------------------------------------------------------------------
    262315-- This function can not be STRICT, because nodatavalueexpr can be NULL (could be just '' though)
  • raster/test/core/testapi.c

    diff -rupN postgis-old/raster/test/core/testapi.c postgis-new/raster/test/core/testapi.c
    old new static void testRasterFromBand(rt_contex 
    987987        rt_raster_destroy(ctx, rast);
    988988}
    989989
     990struct rt_bandstats_t {
     991        uint32_t count;
     992        double min;
     993        double max;
     994        double mean;
     995        double stddev;
     996
     997        double *values;
     998};
     999static void testBandSummaryStats(rt_context ctx) {
     1000        rt_bandstats stats;
     1001
     1002        rt_raster raster;
     1003        rt_band band;
     1004        uint32_t x;
     1005        uint32_t xmax = 10000;
     1006        uint32_t y;
     1007        uint32_t ymax = 10000;
     1008        double nodata;
     1009        int rtn;
     1010
     1011        raster = rt_raster_new(ctx, xmax, ymax);
     1012        assert(raster); /* or we're out of virtual memory */
     1013        band = addBand(ctx, raster, PT_32BUI, 0, 0);
     1014        CHECK(band);
     1015        rt_band_set_nodata(ctx, band, 0);
     1016
     1017        for (x = 0; x < xmax; x++) {
     1018                for (y = 0; y < ymax; y++) {
     1019                        rtn = rt_band_set_pixel(ctx, band, x, y, x + y);
     1020                        CHECK((rtn != -1));
     1021                }
     1022        }
     1023
     1024        nodata = rt_band_get_nodata(ctx, band);
     1025        CHECK_EQUALS(nodata, 0);
     1026
     1027        stats = (rt_bandstats) rt_band_get_summary_stats(ctx, band, 1, 0, 0);
     1028        CHECK(stats);
     1029        CHECK_EQUALS(stats[0].min, 1);
     1030        CHECK_EQUALS(stats[0].max, 19998);
     1031
     1032        stats = (rt_bandstats) rt_band_get_summary_stats(ctx, band, 1, 10, 1);
     1033        CHECK(stats);
     1034        /*
     1035        {
     1036                double stddev = 0;
     1037                for (x = 0; x < stats[0].count; x++) {
     1038                        stddev += pow(stats[0].values[x] - stats[0].mean, 2);
     1039                }
     1040                stddev = sqrt(stddev / (stats[0].count - 1));
     1041                printf("(count, mean, stddev, min, max) = (%d, %f, %f, %f, %f)\n",
     1042                        stats[0].count, stats[0].mean, stats[0].stddev, stats[0].min, stats[0].max);
     1043                printf("stddev: %f\n", stddev);
     1044        }
     1045        */
     1046
     1047        stats = (rt_bandstats) rt_band_get_summary_stats(ctx, band, 1, 15, 0);
     1048        CHECK(stats);
     1049
     1050        stats = (rt_bandstats) rt_band_get_summary_stats(ctx, band, 1, 20, 0);
     1051        CHECK(stats);
     1052
     1053        stats = (rt_bandstats) rt_band_get_summary_stats(ctx, band, 1, 25, 0);
     1054        CHECK(stats);
     1055
     1056        stats = (rt_bandstats) rt_band_get_summary_stats(ctx, band, 0, 0, 0);
     1057        CHECK(stats);
     1058        CHECK_EQUALS(stats[0].min, 0);
     1059        CHECK_EQUALS(stats[0].max, 19998);
     1060
     1061        stats = (rt_bandstats) rt_band_get_summary_stats(ctx, band, 0, 10, 0);
     1062        CHECK(stats);
     1063
     1064        rt_band_destroy(ctx, band);
     1065        rt_raster_destroy(ctx, raster);
     1066}
     1067
    9901068int
    9911069main()
    9921070{
    main() 
    13011379                testRasterFromBand(ctx, raster);
    13021380                printf("Successfully tested rt_raster_from_band\n");
    13031381
     1382                printf("Testing rt_band_get_summary_stats\n");
     1383                testBandSummaryStats(ctx);
     1384                printf("Successfully tested rt_band_get_summary_stats\n");
     1385
    13041386    deepRelease(ctx, raster);
    13051387    rt_context_destroy(ctx);
    13061388
  • raster/test/regress/Makefile.in

    diff -rupN postgis-old/raster/test/regress/Makefile.in postgis-new/raster/test/regress/Makefile.in
    old new TEST_BANDPROPS = \ 
    6161        create_rt_band_properties_test.sql \
    6262        rt_band_properties.sql \
    6363        rt_set_band_properties.sql \
     64        rt_minmax.sql \
    6465        $(NULL)
    6566
    6667TEST_PIXEL = \
  • raster/test/regress/rt_minmax.sql

    diff -rupN postgis-old/raster/test/regress/rt_minmax.sql postgis-new/raster/test/regress/rt_minmax.sql
    old new  
     1SELECT min, max FROM ST_MinMax(
     2        ST_SetValue(
     3                ST_SetValue(
     4                        ST_SetValue(
     5                                ST_AddBand(
     6                                        ST_MakeEmptyRaster(10, 10, 10, 10, 2, 2, 0, 0,-1)
     7                                        , 1, '64BF', 0, 0
     8                                )
     9                                , 1, 1, 1, -10
     10                        )
     11                        , 1, 5, 4, 0
     12                )
     13                , 1, 5, 5, 3.14159
     14        )
     15        , TRUE
     16);
     17SELECT ST_MinMax(
     18        ST_SetValue(
     19                ST_SetValue(
     20                        ST_SetValue(
     21                                ST_AddBand(
     22                                        ST_MakeEmptyRaster(10, 10, 10, 10, 2, 2, 0, 0,-1)
     23                                        , 1, '64BF', 0, 0
     24                                )
     25                                , 1, 1, 1, -10
     26                        )
     27                        , 1, 5, 4, 0
     28                )
     29                , 1, 5, 5, 3.14159
     30        )
     31        , TRUE
     32);
     33SELECT min FROM ST_MinMax(
     34        ST_SetValue(
     35                ST_SetValue(
     36                        ST_SetValue(
     37                                ST_AddBand(
     38                                        ST_MakeEmptyRaster(10, 10, 10, 10, 2, 2, 0, 0,-1)
     39                                        , 1, '64BF', 0, 0
     40                                )
     41                                , 1, 1, 1, -10
     42                        )
     43                        , 1, 5, 4, 0
     44                )
     45                , 1, 5, 5, 3.14159
     46        )
     47        , TRUE
     48);
     49SELECT max,max FROM ST_MinMax(
     50        ST_SetValue(
     51                ST_SetValue(
     52                        ST_SetValue(
     53                                ST_AddBand(
     54                                        ST_MakeEmptyRaster(10, 10, 10, 10, 2, 2, 0, 0,-1)
     55                                        , 1, '64BF', 0, 0
     56                                )
     57                                , 1, 1, 1, -10
     58                        )
     59                        , 1, 5, 4, 0
     60                )
     61                , 1, 5, 5, 3.14159
     62        )
     63        , TRUE
     64);
  • raster/test/regress/rt_minmax_expected

    diff -rupN postgis-old/raster/test/regress/rt_minmax_expected postgis-new/raster/test/regress/rt_minmax_expected
    old new  
     1-10|3.14159
     2(-10,3.14159)
     3-10
     43.14159|3.14159