Geografic-Information-System/lib/raster3d/test/test_put_get_value.c

/*****************************************************************************
*
* MODULE:       Grass raster3d Library
* AUTHOR(S):    Soeren Gebbert, Braunschweig (GER) Jun 2011
* 		        soerengebbert <at> googlemail <dot> com
*
* PURPOSE:	Unit and Integration tests
*
* COPYRIGHT:    (C) 2000 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 <stdlib.h>
#include <string.h>
#include "test_raster3d_lib.h"
#include "grass/interpf.h"

static int test_put_get_value_dcell(void);
static int test_put_get_value_fcell(void);
static int test_put_get_value_resampling(void);
static int test_get_value_region(RASTER3D_Map *map, int cols, int rows, int depths);
static int test_resampling_dcell(RASTER3D_Map *map, double north, double east, double
                                 top, int col, int row, int depth, int fact);
static int test_resampling_fcell(RASTER3D_Map *map, double north, double east, double
                                 top, int col, int row, int depth, int fact);

/* *************************************************************** */
/* Perform the put-get value tests ******************************* */
/* *************************************************************** */
int unit_test_put_get_value()
{
    int sum = 0;

    G_message("\n++ Running raster3d put/get value unit tests ++");

    sum += test_put_get_value_dcell();
    sum += test_put_get_value_fcell();
    sum += test_put_get_value_resampling();


    if (sum > 0)
    G_warning("\n-- raster3d put/get value unit tests failure --");
    else
    G_message("\n-- raster3d put/get value unit tests finished successfully --");

    return sum;
}

/* *************************************************************** */

int test_put_get_value_dcell(void)
{
    int sum = 0;
    int x, y, z;
    DCELL value;
    DCELL value_ref;

    G_message("Testing DCELL put get value functions");

    double north, east, top;
    int col, row, depth;

    RASTER3D_Region region;
    RASTER3D_Map *map = NULL;

    /* We need to set up a specific region for the new raster3d map.
     * First we safe the default region. */
    Rast3d_get_window(&region);

    region.bottom = 0.0;
    region.top = 1000;
    region.south = 1000;
    region.north = 8500;
    region.west = 5000;
    region.east = 10000;
    region.rows = 15;
    region.cols = 10;
    region.depths = 5;

    Rast3d_adjust_region(&region);

    map = Rast3d_open_new_opt_tile_size("test_put_get_value_dcell", RASTER3D_USE_CACHE_XY, &region, DCELL_TYPE, 32);

    /* The window is the same as the map region ... of course */
    Rast3d_set_window_map(map, &region);
    /*
     ROWS
  1000 1500 2000 2500 3000 3500 4000 4500 5000 5500 6000 6500 7000 7500 8000 8500 north
    |....|....|....|....|....|....|....|....|....|....|....|....|....|....|....|
   15   14   13   12   11   10    9    8    7    6    5    4    3    2    1    0 region

    COLS
  5000 5500 6000 6500 7000 7500 8000 8500 9000 9500 10000 east
    |....|....|....|....|....|....|....|....|....|....|
    0    1    2    3    4    5    6    7    8    9   10 region

    DEPTHS
    0   200  400  600  800  1000  top
    |....|....|....|....|....|
    0    1    2    3    4    5 region
    */

    for(z = 0; z < region.depths; z++) {
        for(y = 0; y < region.rows; y++) { /* From the north to the south */
            for(x = 0; x < region.cols; x++) {
                value = -1;
                Rast3d_put_value(map, x, y, z, &value, DCELL_TYPE);
            }
        }
    }

    for(z = 0; z < region.depths; z++) {
        for(y = 0; y < region.rows; y++) { /* From the north to the south */
            for(x = 0; x < region.cols; x++) {
                /* Add cols, rows and depths and put this in the map */
                value = x + y + z;
                Rast3d_put_value(map, x, y, z, &value, DCELL_TYPE);
            }
        }
    }


    /* Reread the new open map and compare the expected results */

    G_message("Get the value of the upper left corner -> 0");


    col = row = depth = 0;
    north = region.north - region.ns_res * row;
    east = region.west + region.ew_res * col;
    top = region.bottom + region.tb_res * depth;

    sum += test_resampling_dcell(map, north, east, top, col, row, depth, 1);


    G_message("Get the value of x == y == z == 1 -> x + y + z == 3");

    col = row = depth = 1;
    north = region.north - region.ns_res * row;
    east = region.west + region.ew_res * col;
    top = region.bottom + region.tb_res * depth;

    sum += test_resampling_dcell(map, north, east, top, col, row, depth, 1);

    G_message("Get the value of x == 4 y == 3 z == 2 -> x + y + z = 9");

    col = 4;
    row = 3;
    depth = 2;
    north = region.north - region.ns_res * row;
    east = region.west + region.ew_res * col;
    top = region.bottom + region.tb_res * depth;

    sum += test_resampling_dcell(map, north, east, top, col, row, depth, 1);


    /* Write everything to the disk and reopen the map */
    Rast3d_flush_all_tiles(map);
    Rast3d_close(map);
    map = Rast3d_open_cell_old("test_put_get_value_dcell", G_mapset(), &region, DCELL_TYPE, RASTER3D_USE_CACHE_XY);


    G_message("Get the value of x == 9 y == 14 z == 4 -> x + y + z = 27");

    col = 9;
    row = 14;
    depth = 4;
    north = region.north - region.ns_res * row;
    east = region.west + region.ew_res * col;
    top = region.bottom + region.tb_res * depth;

    sum += test_resampling_dcell(map, north, east, top, col, row, depth, 1);

    G_message("Get the value of x == 10 y == 15 z == 5 -> x + y + z = NAN");

    col = 10;
    row = 15;
    depth = 5;
    north = region.north - region.ns_res * row;
    east = region.west + region.ew_res * col;
    top = region.bottom + region.tb_res * depth;

    Rast3d_get_region_value(map, north, east, top, &value, DCELL_TYPE);
    Rast3d_get_value(map, col, row, depth, &value_ref, DCELL_TYPE);
    /* Rast3d_get_value_region does not work with coordinates outside the region */
    printf("Value %g == %g\n", value, value_ref);

    if(value == 0 || value < 0 || value > 0) {
        G_message("Error in Rast3d_get_region_value");
        sum++;
    }
    if(value_ref == 0 || value_ref < 0 || value_ref > 0) {
        G_message("Error in Rast3d_get_value");
        sum++;
    }

    Rast3d_close(map);

    G_remove("grid3", "test_put_get_value_dcell");

    return sum;
}

/* *************************************************************** */

int test_put_get_value_fcell(void)
{
     int sum = 0;
    int x, y, z;
    FCELL value;
    FCELL value_ref;

    G_message("Testing FCELL put get value functions");

    double north, east, top;
    int col, row, depth;

    RASTER3D_Region region;
    RASTER3D_Map *map = NULL;

    /* We need to set up a specific region for the new raster3d map.
     * First we safe the default region. */
    Rast3d_get_window(&region);

    region.bottom = 0.0;
    region.top = 1000;
    region.south = 1000;
    region.north = 8500;
    region.west = 5000;
    region.east = 10000;
    region.rows = 15;
    region.cols = 10;
    region.depths = 5;

    Rast3d_adjust_region(&region);

    map = Rast3d_open_new_opt_tile_size("test_put_get_value_fcell", RASTER3D_USE_CACHE_XY, &region, FCELL_TYPE, 32);

    /* The window is the same as the map region ... of course */
    Rast3d_set_window_map(map, &region);

    for(z = 0; z < region.depths; z++) {
        for(y = 0; y < region.rows; y++) {
            for(x = 0; x < region.cols; x++) {
                /* Add cols, rows and depths and put this in the map */
                value = x + y + z;
                Rast3d_put_value(map, x, y, z, &value, FCELL_TYPE);
            }
        }
    }

    /* Write everything to the disk */
    Rast3d_flush_all_tiles(map);
    Rast3d_close(map);

    map = Rast3d_open_cell_old("test_put_get_value_fcell", G_mapset(), &region, DCELL_TYPE, RASTER3D_USE_CACHE_XY);

       /* Reread the map and compare the expected results */

    G_message("Get the value of the lower left corner -> 0");

    col = row = depth = 0;
    north = region.north - region.ns_res * row;
    east = region.west + region.ew_res * col;
    top = region.bottom + region.tb_res * depth;

    sum += test_resampling_fcell(map, north, east, top, col, row, depth, 1);

    G_message("Get the value of x == y == z == 1 -> x + y + z == 3");

    col = row = depth = 1;
    north = region.north - region.ns_res * row;
    east = region.west + region.ew_res * col;
    top = region.bottom + region.tb_res * depth;

    sum += test_resampling_fcell(map, north, east, top, col, row, depth, 1);

    G_message("Get the value of x == 4 y == 3 z == 2 -> x + y + z = 9");

    col = 4;
    row = 3;
    depth = 2;
    north = region.north - region.ns_res * row;
    east = region.west + region.ew_res * col;
    top = region.bottom + region.tb_res * depth;

    sum += test_resampling_fcell(map, north, east, top, col, row, depth, 1);

    G_message("Get the value of x == 9 y == 14 z == 4 -> x + y + z = 27");

    col = 9;
    row = 14;
    depth = 4;
    north = region.north - region.ns_res * row;
    east = region.west + region.ew_res * col;
    top = region.bottom + region.tb_res * depth;

    sum += test_resampling_fcell(map, north, east, top, col, row, depth, 1);

    G_message("Get the value of x == 10 y == 15 z == 5 -> x + y + z = NAN");

    col = 10;
    row = 15;
    depth = 5;
    north = region.north - region.ns_res * row;
    east = region.west + region.ew_res * col;
    top = region.bottom + region.tb_res * depth;

    Rast3d_get_region_value(map, north, east, top, &value, FCELL_TYPE);
    Rast3d_get_value(map, 10, 15, 5, &value_ref, FCELL_TYPE);
    /* Rast3d_get_value_region does not work with coordinates outside the region */
    printf("Value %g == %g\n", value, value_ref);

    if(value == 0 || value < 0 || value > 0) {
        G_message("Error in Rast3d_get_region_value");
        sum++;
    }
    if(value_ref == 0 || value_ref < 0 || value_ref > 0) {
        G_message("Error in Rast3d_get_value");
        sum++;
    }

    Rast3d_close(map);

    G_remove("grid3", "test_put_get_value_fcell");

    return sum;
}

/* *************************************************************** */

int test_put_get_value_resampling(void)
{
    int sum = 0;
    int x, y, z;
    DCELL value;

    G_message("Testing put get resample value functions");

    double north, east, top;
    int col, row, depth;

    RASTER3D_Region region;
    RASTER3D_Region window;
    RASTER3D_Map *map = NULL;

    /* We need to set up a specific region for the new raster3d map.
     * First we safe the default region. */
    Rast3d_get_window(&region);

    region.bottom = 0.0;
    region.top = 1000;
    region.south = 1000;
    region.north = 8500;
    region.west = 5000;
    region.east = 10000;
    region.rows = 15;
    region.cols = 10;
    region.depths = 5;

    Rast3d_adjust_region(&region);

    map = Rast3d_open_new_opt_tile_size("test_put_get_value_resample", RASTER3D_USE_CACHE_XY, &region, DCELL_TYPE, 32);
    /*
     ROWS
  1000 1500 2000 2500 3000 3500 4000 4500 5000 5500 6000 6500 7000 7500 8000 8500 north
    |....|....|....|....|....|....|....|....|....|....|....|....|....|....|....|
   15   14   13   12   11   10    9    8    7    6    5    4    3    2    1    0 region
    |    |    |    |    |    |    |    |    |    |    |    |    |    |    |    |
   30   28   26   24   22   20   18   16   14   12   10    8    6    4    2    0 window

    COLS
  5000 5500 6000 6500 7000 7500 8000 8500 9000 9500 10000 east
    |....|....|....|....|....|....|....|....|....|....|
    0    1    2    3    4    5    6    7    8    9   10 region
    |    |    |    |    |    |    |    |    |    |    |
    0    2    4    6    8   10   12   14   16   18   20 window
    DEPTHS
    0   200  400  600  800  1000 top
    |....|....|....|....|....|
    0    1    2    3    4    5 region
    |    |    |    |    |    |
    0    2    4    6    8   10 window
    */

    for(z = 0; z < region.depths; z++) {
        for(y = 0; y < region.rows; y++) {  /* North to south */
            for(x = 0; x < region.cols; x++) {
                /* Add cols, rows and depths and put this in the map */
                value = x + y + z;
                Rast3d_put_double(map, x, y, z, value);
            }
        }
    }

    /* Write everything to the disk */
    Rast3d_flush_all_tiles(map);
    Rast3d_close(map);

    /* We modify the window for resampling tests */
    Rast3d_region_copy(&window, &region);

    /* Double the cols, rows and depths -> 8x resolution window */
    window.rows = 30;
    window.cols = 20;
    window.depths = 10;

    Rast3d_adjust_region(&window);

    map = Rast3d_open_cell_old("test_put_get_value_resample", G_mapset(), &region, DCELL_TYPE, RASTER3D_USE_CACHE_XY);

    /* The window has the 8x resolution as the map region */
    Rast3d_set_window_map(map, &window);

    /* Reread the map and compare the expected results */

    G_message("Get the value of the upper left corner -> 0");

    col = row = depth = 0;
    north = region.north - region.ns_res * row;
    east = region.west + region.ew_res * col;
    top = region.bottom + region.tb_res * depth;

    sum += test_resampling_dcell(map, north, east, top, col, row, depth, 2);

    G_message("Get the value of x == y == z == 1 -> x + y + z == 3");


    col = row = depth = 1;
    north = region.north - region.ns_res * row;
    east = region.west + region.ew_res * col;
    top = region.bottom + region.tb_res * depth;

    sum += test_resampling_dcell(map, north, east, top, col, row, depth, 2);

    G_message("Get the value of x == 7 y == 9 z == 3 -> x + y + z == 19");

    col = 7;
    row = 9;
    depth = 3;
    north = region.north - region.ns_res * row;
    east = region.west + region.ew_res * col;
    top = region.bottom + region.tb_res * depth;

    sum += test_resampling_dcell(map, north, east, top, col, row, depth, 2);

    G_message("Get the value of x == 9 y == 14 z == 4 -> x + y + z == 27");

    col = 9;
    row = 14;
    depth = 4;
    north = region.north - region.ns_res * row;
    east = region.west + region.ew_res * col;
    top = region.bottom + region.tb_res * depth;

    sum += test_resampling_dcell(map, north, east, top, col, row, depth, 2);

    sum += test_get_value_region(map, region.cols, region.rows, region.depths);

    Rast3d_close(map);

    G_remove("grid3", "test_put_get_value_resample");

    return sum;
}

/* *************************************************************** */

int test_resampling_dcell(RASTER3D_Map *map, double north, double east, double top, int col, int row, int depth, int fact)
{
    int sum = 0;
    DCELL value;
    DCELL value_ref;
    DCELL value_reg;
    DCELL value_win;

    Rast3d_get_region_value(map, north, east, top, &value, DCELL_TYPE);
    Rast3d_get_window_value(map, north, east, top, &value_win, DCELL_TYPE);
    Rast3d_get_value(map, col * fact, row * fact, depth * fact, &value_ref, DCELL_TYPE);
    Rast3d_get_value_region(map, col, row, depth, &value_reg, DCELL_TYPE);
    printf("Value %g == %g == %g == %g\n", value, value_win, value_ref, value_reg);

    if(value != col + row + depth) {
        G_message("Error in Rast3d_get_region_value");
        sum++;
    }
    if(value_win != col + row + depth) {
        G_message("Error in Rast3d_get_window_value");
        sum++;
    }
    if(value_ref != col + row + depth) {
        G_message("Error in Rast3d_get_value");
        sum++;
    }
    if(value_reg != col + row + depth) {
        G_message("Error in Rast3d_get_value_region");
        sum++;
    }

    return sum;
}

/* *************************************************************** */

int test_resampling_fcell(RASTER3D_Map *map, double north, double east, double top, int col, int row, int depth, int fact)
{
    int sum = 0;
    FCELL value;
    FCELL value_ref;
    FCELL value_reg;
    FCELL value_win;

    Rast3d_get_region_value(map, north, east, top, &value, FCELL_TYPE);
    Rast3d_get_window_value(map, north, east, top, &value_win, FCELL_TYPE);
    Rast3d_get_value(map, col * fact, row * fact, depth * fact, &value_ref, FCELL_TYPE);
    Rast3d_get_value_region(map, col, row, depth, &value_reg, FCELL_TYPE);
    printf("Value %g == %g == %g == %g\n", value, value_win, value_ref, value_reg);

    if(value != col + row + depth) {
        G_message("Error in Rast3d_get_region_value");
        sum++;
    }
    if(value_win != col + row + depth) {
        G_message("Error in Rast3d_get_window_value");
        sum++;
    }
    if(value_ref != col + row + depth) {
        G_message("Error in Rast3d_get_value");
        sum++;
    }
    if(value_reg != col + row + depth) {
        G_message("Error in Rast3d_get_value_region");
        sum++;
    }

    return sum;
}

/* *************************************************************** */

int test_get_value_region(RASTER3D_Map *map, int cols, int rows, int depths)
{
    int sum = 0;
    FCELL fvalue1 = 0.0;
    FCELL fvalue2 = 0.0;
    DCELL dvalue1 = 0.0;
    DCELL dvalue2 = 0.0;

    /* Test for correct Null value */
    Rast3d_get_value_region(map, -1, -1, -1, &fvalue1, FCELL_TYPE);
    Rast3d_get_value_region(map, cols, rows, depths, &fvalue2, FCELL_TYPE);
    Rast3d_get_value_region(map, -1, -1, -1, &dvalue1, DCELL_TYPE);
    Rast3d_get_value_region(map, cols, rows, depths, &dvalue2, DCELL_TYPE);
    printf("Value %g == %g == %g == %g\n", fvalue1, fvalue2, dvalue1, dvalue2);

    if(!Rast_is_f_null_value(&fvalue1)) {
        G_message("Error in Rast3d_get_value_region");
        sum++;
    }
    if(!Rast_is_f_null_value(&fvalue2)) {
        G_message("Error in Rast3d_get_value_region");
        sum++;
    }
    if(!Rast_is_d_null_value(&dvalue1)) {
        G_message("Error in Rast3d_get_value_region");
        sum++;
    }
    if(!Rast_is_d_null_value(&dvalue2)) {
        G_message("Error in Rast3d_get_value_region");
        sum++;
    }

    return sum;
}