rt_pixel.c

Go to the documentation of this file.

/*
 *
 * WKTRaster - Raster Types for PostGIS
 * http://trac.osgeo.org/postgis/wiki/WKTRaster
 *
 * Copyright (C) 2011-2013 Regents of the University of California
 *   <bkpark@ucdavis.edu>
 * Copyright (C) 2010-2011 Jorge Arevalo <jorge.arevalo@deimos-space.com>
 * Copyright (C) 2010-2011 David Zwarg <dzwarg@azavea.com>
 * Copyright (C) 2009-2011 Pierre Racine <pierre.racine@sbf.ulaval.ca>
 * Copyright (C) 2009-2011 Mateusz Loskot <mateusz@loskot.net>
 * Copyright (C) 2008-2009 Sandro Santilli <strk@kbt.io>
 * Copyright (C) 2013  Nathaniel Hunter Clay <clay.nathaniel@gmail.com>
 * Copyright (C) 2025 Darafei Praliaskouski <me@komzpa.net>
 *
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License
 * as published by the Free Software Foundation; either version 2
 * of the License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software Foundation,
 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 *
 */
 
#include "librtcore.h"
#include "librtcore_internal.h"
 
/******************************************************************************
* rt_pixeltype
******************************************************************************/
 
int
rt_pixtype_size(rt_pixtype pixtype) {
        int pixbytes = -1;
 
        switch (pixtype) {
                case PT_1BB:
                case PT_2BUI:
                case PT_4BUI:
                case PT_8BSI:
                case PT_8BUI:
                        pixbytes = 1;
                        break;
                case PT_16BSI:
                case PT_16BUI:
                case PT_16BF:
                        pixbytes = 2;
                        break;
                case PT_32BSI:
                case PT_32BUI:
                case PT_32BF:
                        pixbytes = 4;
                        break;
                case PT_64BF:
                        pixbytes = 8;
                        break;
                default:
                        rterror("rt_pixtype_size: Unknown pixeltype %d", pixtype);
                        pixbytes = -1;
                        break;
        }
 
        RASTER_DEBUGF(3, "Pixel type = %s and size = %d bytes",
                rt_pixtype_name(pixtype), pixbytes);
 
        return pixbytes;
}
 
int
rt_pixtype_alignment(rt_pixtype pixtype) {
        return rt_pixtype_size(pixtype);
}
 
rt_pixtype
rt_pixtype_index_from_name(const char* pixname) {
        assert(pixname);
 
        if (strcmp(pixname, "1BB") == 0)
                return PT_1BB;
        else if (strcmp(pixname, "2BUI") == 0)
                return PT_2BUI;
        else if (strcmp(pixname, "4BUI") == 0)
                return PT_4BUI;
        else if (strcmp(pixname, "8BSI") == 0)
                return PT_8BSI;
        else if (strcmp(pixname, "8BUI") == 0)
                return PT_8BUI;
        else if (strcmp(pixname, "16BSI") == 0)
                return PT_16BSI;
        else if (strcmp(pixname, "16BUI") == 0)
                return PT_16BUI;
        else if (strcmp(pixname, "16BF") == 0)
                return PT_16BF;
        else if (strcmp(pixname, "32BSI") == 0)
                return PT_32BSI;
        else if (strcmp(pixname, "32BUI") == 0)
                return PT_32BUI;
        else if (strcmp(pixname, "32BF") == 0)
                return PT_32BF;
        else if (strcmp(pixname, "64BF") == 0)
                return PT_64BF;
 
        return PT_END;
}
 
const char*
rt_pixtype_name(rt_pixtype pixtype) {
 
 
        switch (pixtype) {
                case PT_1BB:
                        return "1BB";
                case PT_2BUI:
                        return "2BUI";
                case PT_4BUI:
                        return "4BUI";
                case PT_8BSI:
                        return "8BSI";
                case PT_8BUI:
                        return "8BUI";
                case PT_16BSI:
                        return "16BSI";
                case PT_16BUI:
                        return "16BUI";
                case PT_16BF:
                        return "16BF";
                case PT_32BSI:
                        return "32BSI";
                case PT_32BUI:
                        return "32BUI";
                case PT_32BF:
                        return "32BF";
                case PT_64BF:
                        return "64BF";
                default:
                        rterror("rt_pixtype_name: Unknown pixeltype %d", pixtype);
                        return "Unknown";
        }
}
 
double
rt_pixtype_get_min_value(rt_pixtype pixtype) {
        switch (pixtype) {
                case PT_1BB: {
                        return (double) rt_util_clamp_to_1BB((double) CHAR_MIN);
                }
                case PT_2BUI: {
                        return 0;
                }
                case PT_4BUI: {
                        return 0;
                }
                case PT_8BUI: {
                        return 0;
                }
                case PT_8BSI: {
                        return (double) rt_util_clamp_to_8BSI((double) SCHAR_MIN);
                }
                case PT_16BSI: {
                        return (double) rt_util_clamp_to_16BSI((double) SHRT_MIN);
                }
                case PT_16BUI: {
                        return 0;
                }
                case PT_16BF: {
                        return (double)-POSTGIS_RT_16F_MAX;
                }
                case PT_32BSI: {
                        return (double)rt_util_clamp_to_32BSI((double)INT_MIN);
                }
                case PT_32BUI: {
                        return 0;
                }
                case PT_32BF: {
                        return (double) -FLT_MAX;
                }
                case PT_64BF: {
                        return (double) -DBL_MAX;
                }
                default: {
                        rterror("rt_pixtype_get_min_value: Unknown pixeltype %d", pixtype);
                        return (double) rt_util_clamp_to_8BUI((double) CHAR_MIN);
                }
        }
}
 
rt_errorstate rt_pixtype_compare_clamped_values(
        rt_pixtype pixtype,
        double val, double refval,
        int *isequal
) {
        assert(isequal != NULL);
        *isequal = 0;
 
        switch (pixtype) {
                case PT_1BB:
                        if (rt_util_clamp_to_1BB(val) == rt_util_clamp_to_1BB(refval))
                                *isequal = 1;
                        break;
                case PT_2BUI:
                        if (rt_util_clamp_to_2BUI(val) == rt_util_clamp_to_2BUI(refval))
                                *isequal = 1;
                        break;
                case PT_4BUI:
                        if (rt_util_clamp_to_4BUI(val) == rt_util_clamp_to_4BUI(refval))
                                *isequal = 1;
                        break;
                case PT_8BSI:
                        if (rt_util_clamp_to_8BSI(val) == rt_util_clamp_to_8BSI(refval))
                                *isequal = 1;
                        break;
                case PT_8BUI:
                        if (rt_util_clamp_to_8BUI(val) == rt_util_clamp_to_8BUI(refval))
                                *isequal = 1;
                        break;
                case PT_16BSI:
                        if (rt_util_clamp_to_16BSI(val) == rt_util_clamp_to_16BSI(refval))
                                *isequal = 1;
                        break;
                case PT_16BUI:
                        if (rt_util_clamp_to_16BUI(val) == rt_util_clamp_to_16BUI(refval))
                                *isequal = 1;
                        break;
                case PT_16BF:
                        if (FLT_EQ(rt_util_clamp_to_16F(val), rt_util_clamp_to_16F(refval)))
                                *isequal = 1;
                        break;
                case PT_32BSI:
                        if (rt_util_clamp_to_32BSI(val) == rt_util_clamp_to_32BSI(refval))
                                *isequal = 1;
                        break;
                case PT_32BUI:
                        if (rt_util_clamp_to_32BUI(val) == rt_util_clamp_to_32BUI(refval))
                                *isequal = 1;
                        break;
                case PT_32BF:
                        if (FLT_EQ(rt_util_clamp_to_32F(val), rt_util_clamp_to_32F(refval)))
                                *isequal = 1;
                        break;
                case PT_64BF:
                        if (FLT_EQ(val, refval))
                                *isequal = 1;
                        break;
                default:
                        rterror("rt_pixtype_compare_clamped_values: Unknown pixeltype %d", pixtype);
                        return ES_ERROR;
        }
 
        return ES_NONE;
}
 
/******************************************************************************
* rt_pixel
******************************************************************************/
 
/*
 * Convert an array of rt_pixel objects to two 2D arrays of value and NODATA.
 * The dimensions of the returned 2D array are [Y][X], going by row Y and
 * then column X.
 *
 * @param npixel : array of rt_pixel objects
 * @param count : number of elements in npixel
 * @param mask : mask to be respected when returning array
 * @param x : the column of the center pixel (0-based)
 * @param y : the line of the center pixel (0-based)
 * @param distancex : the number of pixels around the specified pixel
 * along the X axis
 * @param distancey : the number of pixels around the specified pixel
 * along the Y axis
 * @param value : pointer to pointer for 2D value array
 * @param nodata : pointer to pointer for 2D NODATA array
 * @param dimx : size of value and nodata along the X axis
 * @param dimy : size of value and nodata along the Y axis
 *
 * @return ES_NONE on success, ES_ERROR on error
 */
rt_errorstate rt_pixel_set_to_array(
        rt_pixel npixel, uint32_t count, rt_mask mask,
        int x, int y,
        uint16_t distancex, uint16_t distancey,
        double ***value,
        int ***nodata,
        int *dimx, int *dimy
) {
        uint32_t i;
        uint32_t j;
        uint32_t dim[2] = {0};
        double **values = NULL;
        int **nodatas = NULL;
        int zero[2] = {0};
        int _x;
        int _y;
 
        assert(npixel != NULL && count > 0);
        assert(value != NULL);
        assert(nodata != NULL);
 
        /* dimensions */
        dim[0] = distancex * 2 + 1;
        dim[1] = distancey * 2 + 1;
        RASTER_DEBUGF(4, "dimensions = %d x %d", dim[0], dim[1]);
 
        /* make sure that the dimx and dimy match mask */
        if (mask != NULL) {
          if (dim[0] != mask->dimx || dim[1] != mask->dimy) {
            rterror("rt_pixel_set_array: mask dimensions %d x %d do not match given dims %d x %d", mask->dimx, mask->dimy,  dim[0],  dim[1]);
            return ES_ERROR;
          }
 
          if (mask->values == NULL || mask->nodata == NULL) {
            rterror("rt_pixel_set_array: Invalid mask");
            return ES_ERROR;
          }
 
        }
 
        /* establish 2D arrays (Y axis) */
        values = rtalloc(sizeof(double *) * dim[1]);
        nodatas = rtalloc(sizeof(int *) * dim[1]);
 
        if (values == NULL || nodatas == NULL) {
                rterror("rt_pixel_set_to_array: Could not allocate memory for 2D array");
                return ES_ERROR;
        }
 
        /* initialize X axis */
        for (i = 0; i < dim[1]; i++) {
                values[i] = rtalloc(sizeof(double) * dim[0]);
                nodatas[i] = rtalloc(sizeof(int) * dim[0]);
 
                if (values[i] == NULL || nodatas[i] == NULL) {
                        rterror("rt_pixel_set_to_array: Could not allocate memory for dimension of 2D array");
 
                        if (values[i] == NULL) {
                                for (j = 0; j < i; j++) {
                                        rtdealloc(values[j]);
                                        rtdealloc(nodatas[j]);
                                }
                        }
                        else {
                                for (j = 0; j <= i; j++) {
                                        rtdealloc(values[j]);
                                        if (j < i)
                                                rtdealloc(nodatas[j]);
                                }
                        }
 
                        rtdealloc(values);
                        rtdealloc(nodatas);
 
                        return ES_ERROR;
                }
 
                /* set values to 0 */
                memset(values[i], 0, sizeof(double) * dim[0]);
 
                /* set nodatas to 1 */
                for (j = 0; j < dim[0]; j++)
                        nodatas[i][j] = 1;
        }
 
        /* get 0,0 of grid */
        zero[0] = x - distancex;
        zero[1] = y - distancey;
 
        /* populate 2D arrays */
        for (i = 0; i < count; i++) {
                if (npixel[i].nodata)
                        continue;
 
                _x = npixel[i].x - zero[0];
                _y = npixel[i].y - zero[1];
 
                RASTER_DEBUGF(4, "absolute x,y: %d x %d", npixel[i].x, npixel[i].y);
                RASTER_DEBUGF(4, "relative x,y: %d x %d", _x, _y);
 
                /* no mask */
                if (mask == NULL) {
                        values[_y][_x] = npixel[i].value;
                        nodatas[_y][_x] = 0;
                }
                /* mask */
                else {
                        /* unweighted (boolean) mask */
                        if (mask->weighted == 0) {
                                /* pixel is set to zero or nodata */
                                if (FLT_EQ(mask->values[_y][_x], 0.0) || mask->nodata[_y][_x] == 1)
                                {
                                        values[_y][_x] = 0;
                                        nodatas[_y][_x] = 1;
                                }
                                /* use pixel */
                                else {
                                        values[_y][_x] = npixel[i].value;
                                        nodatas[_y][_x] = 0;
                                }
                        }
                        /* weighted mask */
                        else {
                                /* nodata */
                                if(mask->nodata[_y][_x] == 1) {
                                        values[_y][_x] = 0;
                                        nodatas[_y][_x] = 1;
                                }
                                /* apply weight to pixel value */
                                else {
                                        values[_y][_x] = npixel[i].value * mask->values[_y][_x];
                                        nodatas[_y][_x] = 0;
                                }
                        }
                }
 
                RASTER_DEBUGF(4, "(x, y, nodata, value) = (%d, %d, %d, %f)", _x, _y, nodatas[_y][_x], values[_y][_x]);
        }
 
        *value = &(*values);
        *nodata = &(*nodatas);
        if (dimx != NULL)
                *dimx = dim[0];
        if (dimy != NULL)
                *dimy = dim[1];
 
        return ES_NONE;
}