rt_band.c

Go to the documentation of this file.

/*
 *
 * WKTRaster - Raster Types for PostGIS
 * http://trac.osgeo.org/postgis/wiki/WKTRaster
 *
 * Copyright (C) 2018 Bborie Park <dustymugs@gmail.com>
 * 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) 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.
 *
 */
 
// For stat64()
#define _LARGEFILE64_SOURCE 1
 
#include <stdio.h>
 
#include "librtcore.h"
#include "librtcore_internal.h"
 
#include "cpl_vsi.h"
#include "gdal_vrt.h"
 
rt_band
rt_band_new_inline(
        uint16_t width, uint16_t height,
        rt_pixtype pixtype,
        uint32_t hasnodata, double nodataval,
        uint8_t* data
) {
        rt_band band = NULL;
 
        assert(NULL != data);
 
        band = rtalloc(sizeof(struct rt_band_t));
        if (band == NULL) {
                rterror("rt_band_new_inline: Out of memory allocating rt_band");
                return NULL;
        }
 
        RASTER_DEBUGF(3, "Created rt_band @ %p with pixtype %s", band, rt_pixtype_name(pixtype));
 
        band->pixtype = pixtype;
        band->offline = 0;
        band->width = width;
        band->height = height;
        band->hasnodata = hasnodata ? 1 : 0;
        band->isnodata = FALSE; /* we don't know what is in data, so must be FALSE */
        band->nodataval = 0;
        band->data.mem = data;
        band->ownsdata = 0; /* we do NOT own this data!!! */
        band->raster = NULL;
 
        RASTER_DEBUGF(3, "Created rt_band with dimensions %d x %d", band->width, band->height);
 
        /* properly set nodataval as it may need to be constrained to the data type */
        if (hasnodata && rt_band_set_nodata(band, nodataval, NULL) != ES_NONE) {
                rterror("rt_band_new_inline: Could not set NODATA value");
                rt_band_destroy(band);
                return NULL;
        }
 
        return band;
}
 
void
rt_band_init_value(
        rt_band band,
        double initval
) {
        assert(band != NULL);
        assert(band->data.mem != NULL);
        rt_pixtype pixtype = band->pixtype;
        uint32_t width = band->width;
        uint32_t height = band->height;
        uint32_t numval = width * height;
        void *mem = band->data.mem;
        size_t memsize = numval * rt_pixtype_size(pixtype);
 
        /* initialize to nodataval */
        int32_t checkvalint = 0;
        uint32_t checkvaluint = 0;
        double checkvaldouble = 0;
        float checkvalfloat = 0;
 
        /* initialize to zero */
        if (FLT_EQ(initval, 0.0)) {
                memset(mem, 0, memsize);
                return;
        }
 
        switch (pixtype) {
                case PT_1BB:
                {
                        uint8_t *ptr = mem;
                        uint8_t clamped_initval = rt_util_clamp_to_1BB(initval);
                        for (uint32_t i = 0; i < numval; i++)
                                ptr[i] = clamped_initval;
                        checkvalint = ptr[0];
                        break;
                }
                case PT_2BUI:
                {
                        uint8_t *ptr = mem;
                        uint8_t clamped_initval = rt_util_clamp_to_2BUI(initval);
                        for (uint32_t i = 0; i < numval; i++)
                                ptr[i] = clamped_initval;
                        checkvalint = ptr[0];
                        break;
                }
                case PT_4BUI:
                {
                        uint8_t *ptr = mem;
                        uint8_t clamped_initval = rt_util_clamp_to_4BUI(initval);
                        for (uint32_t i = 0; i < numval; i++)
                                ptr[i] = clamped_initval;
                        checkvalint = ptr[0];
                        break;
                }
                case PT_8BSI:
                {
                        int8_t *ptr = mem;
                        int8_t clamped_initval = rt_util_clamp_to_8BSI(initval);
                        for (uint32_t i = 0; i < numval; i++)
                                ptr[i] = clamped_initval;
                        checkvalint = ptr[0];
                        break;
                }
                case PT_8BUI:
                {
                        uint8_t *ptr = mem;
                        uint8_t clamped_initval = rt_util_clamp_to_8BUI(initval);
                        for (uint32_t i = 0; i < numval; i++)
                                ptr[i] = clamped_initval;
                        checkvalint = ptr[0];
                        break;
                }
                case PT_16BSI:
                {
                        int16_t *ptr = mem;
                        int16_t clamped_initval = rt_util_clamp_to_16BSI(initval);
                        for (uint32_t i = 0; i < numval; i++)
                                ptr[i] = clamped_initval;
                        checkvalint = ptr[0];
                        break;
                }
                case PT_16BUI: {
                        uint16_t *ptr = mem;
                        uint16_t clamped_initval = rt_util_clamp_to_16BUI(initval);
                        for (uint32_t i = 0; i < numval; i++)
                                ptr[i] = clamped_initval;
                        checkvalint = ptr[0];
                        break;
                }
                case PT_16BF: {
                        uint16_t *ptr = mem;
                        float clamped_initval = rt_util_clamp_to_16F(initval);
                        uint16_t packed = rt_util_float_to_float16(clamped_initval);
                        for (uint32_t i = 0; i < numval; i++)
                                ptr[i] = packed;
                        checkvalfloat = rt_util_float16_to_float(ptr[0]);
                        break;
                }
                case PT_32BSI: {
                        int32_t *ptr = mem;
                        int32_t clamped_initval = rt_util_clamp_to_32BSI(initval);
                        for (uint32_t i = 0; i < numval; i++)
                                ptr[i] = clamped_initval;
                        checkvalint = ptr[0];
                        break;
                }
                case PT_32BUI:
                {
                        uint32_t *ptr = mem;
                        uint32_t clamped_initval = rt_util_clamp_to_32BUI(initval);
                        for (uint32_t i = 0; i < numval; i++)
                                ptr[i] = clamped_initval;
                        checkvaluint = ptr[0];
                        break;
                }
                case PT_32BF:
                {
                        float *ptr = mem;
                        float clamped_initval = rt_util_clamp_to_32F(initval);
                        for (uint32_t i = 0; i < numval; i++)
                                ptr[i] = clamped_initval;
                        checkvalfloat = ptr[0];
                        break;
                }
                case PT_64BF:
                {
                        double *ptr = mem;
                        for (uint32_t i = 0; i < numval; i++)
                                ptr[i] = initval;
                        checkvaldouble = ptr[0];
                        break;
                }
                default:
                {
                        rterror("%s: Unknown pixeltype %d", __func__, pixtype);
                        return;
                }
        }
 
        /* Overflow checking */
        rt_util_dbl_trunc_warning(
                initval,
                checkvalint, checkvaluint,
                checkvalfloat, checkvaldouble,
                pixtype
        );
 
        return;
}
 
 
rt_band
rt_band_new_offline(
        uint16_t width, uint16_t height,
        rt_pixtype pixtype,
        uint32_t hasnodata, double nodataval,
        uint8_t bandNum, const char* path
) {
        rt_band band = NULL;
        int pathlen = 0;
 
        assert(NULL != path);
 
        band = rtalloc(sizeof(struct rt_band_t));
        if (band == NULL) {
                rterror("rt_band_new_offline: Out of memory allocating rt_band");
                return NULL;
        }
 
        RASTER_DEBUGF(3, "Created rt_band @ %p with pixtype %s",
                band, rt_pixtype_name(pixtype)
        );
 
        band->pixtype = pixtype;
        band->offline = 1;
        band->width = width;
        band->height = height;
        band->hasnodata = hasnodata ? 1 : 0;
        band->nodataval = 0;
        band->isnodata = FALSE; /* we don't know if the offline band is NODATA */
        band->ownsdata = 0; /* offline, flag is useless as all offline data cache is owned internally */
        band->raster = NULL;
 
        /* properly set nodataval as it may need to be constrained to the data type */
        if (hasnodata && rt_band_set_nodata(band, nodataval, NULL) != ES_NONE) {
                rterror("rt_band_new_offline: Could not set NODATA value");
                rt_band_destroy(band);
                return NULL;
        }
 
        band->data.offline.bandNum = bandNum;
 
        /* memory for data.offline.path is managed internally */
        pathlen = strlen(path);
        band->data.offline.path = rtalloc(sizeof(char) * (pathlen + 1));
        if (band->data.offline.path == NULL) {
                rterror("rt_band_new_offline: Out of memory allocating offline path");
                rt_band_destroy(band);
                return NULL;
        }
        memcpy(band->data.offline.path, path, pathlen);
        band->data.offline.path[pathlen] = '\0';
 
        band->data.offline.mem = NULL;
 
        return band;
}
 
rt_band
rt_band_new_offline_from_path(
        uint16_t width,
        uint16_t height,
        int hasnodata,
        double nodataval,
        uint8_t bandNum,
        const char* path,
        int force
) {
        GDALDatasetH hdsSrc = NULL;
        int nband = 0;
        GDALRasterBandH hbandSrc = NULL;
 
        GDALDataType gdpixtype;
        rt_pixtype pt = PT_END;
 
        /* open outdb raster file */
        rt_util_gdal_register_all(0);
        hdsSrc = rt_util_gdal_open(path, GA_ReadOnly, 1);
        if (hdsSrc == NULL && !force) {
                rterror("rt_band_new_offline_from_path: Cannot open offline raster: %s", path);
                return NULL;
        }
 
        nband = GDALGetRasterCount(hdsSrc);
        if (!nband && !force) {
                rterror("rt_band_new_offline_from_path: No bands found in offline raster: %s", path);
                GDALClose(hdsSrc);
                return NULL;
        }
        /* bandNum is 1-based */
        else if (bandNum > nband && !force) {
                rterror(
                        "rt_band_new_offline_from_path: Specified band %d not found in offline raster: %s",
                        bandNum,
                        path
                );
                GDALClose(hdsSrc);
                return NULL;
        }
 
        hbandSrc = GDALGetRasterBand(hdsSrc, bandNum);
        if (hbandSrc == NULL && !force) {
                rterror(
                        "rt_band_new_offline_from_path: Cannot get band %d from GDAL dataset",
                        bandNum
                );
                GDALClose(hdsSrc);
                return NULL;
        }
 
        gdpixtype = GDALGetRasterDataType(hbandSrc);
        pt = rt_util_gdal_datatype_to_pixtype(gdpixtype);
        if (pt == PT_END && !force) {
                rterror(
                        "rt_band_new_offline_from_path: Unsupported pixel type %s of band %d from GDAL dataset",
                        GDALGetDataTypeName(gdpixtype),
                        bandNum
                );
                GDALClose(hdsSrc);
                return NULL;
        }
 
        /* use out-db band's nodata value if nodataval not already set */
        if (!hasnodata)
                nodataval = GDALGetRasterNoDataValue(hbandSrc, &hasnodata);
 
        GDALClose(hdsSrc);
 
        return rt_band_new_offline(
                width, height,
                pt,
                hasnodata, nodataval,
                bandNum - 1, path
        );
}
 
rt_band
rt_band_duplicate(rt_band band) {
        rt_band rtn = NULL;
 
        assert(band != NULL);
 
        /* offline */
        if (band->offline) {
                rtn = rt_band_new_offline(
                        band->width, band->height,
                        band->pixtype,
                        band->hasnodata, band->nodataval,
                        band->data.offline.bandNum, (const char *) band->data.offline.path
                );
        }
        /* online */
        else {
                uint8_t *data = NULL;
                data = rtalloc((size_t)rt_pixtype_size(band->pixtype) * band->width * band->height);
                if (data == NULL) {
                        rterror("rt_band_duplicate: Out of memory allocating online band data");
                        return NULL;
                }
                memcpy(data, band->data.mem, (size_t)rt_pixtype_size(band->pixtype) * band->width * band->height);
 
                rtn = rt_band_new_inline(
                        band->width, band->height,
                        band->pixtype,
                        band->hasnodata, band->nodataval,
                        data
                );
                rt_band_set_ownsdata_flag(rtn, 1); /* we DO own this data!!! */
        }
 
        if (rtn == NULL) {
                rterror("rt_band_duplicate: Could not copy band");
                return NULL;
        }
 
        return rtn;
}
 
int
rt_band_is_offline(rt_band band) {
  assert(NULL != band);
        return band->offline ? 1 : 0;
}
 
void
rt_band_destroy(rt_band band) {
        if (band == NULL)
                return;
 
        RASTER_DEBUGF(3, "Destroying rt_band @ %p", band);
 
        /* offline band */
        if (band->offline) {
                /* memory cache */
                if (band->data.offline.mem != NULL)
                        rtdealloc(band->data.offline.mem);
                /* offline file path */
                if (band->data.offline.path != NULL)
                        rtdealloc(band->data.offline.path);
        }
        /* inline band and band owns the data */
        else if (band->data.mem != NULL && band->ownsdata)
                rtdealloc(band->data.mem);
 
        rtdealloc(band);
}
 
const char*
rt_band_get_ext_path(rt_band band) {
 
    assert(NULL != band);
 
 
    if (!band->offline) {
        RASTER_DEBUG(3, "rt_band_get_ext_path: Band is not offline");
        return NULL;
    }
    return band->data.offline.path;
}
 
rt_errorstate
rt_band_get_ext_band_num(rt_band band, uint8_t *bandnum) {
        assert(NULL != band);
        assert(NULL != bandnum);
 
        *bandnum = 0;
 
        if (!band->offline) {
                RASTER_DEBUG(3, "rt_band_get_ext_band_num: Band is not offline");
                return ES_ERROR;
        }
 
        *bandnum = band->data.offline.bandNum;
 
        return ES_NONE;
}
 
void *
rt_band_get_data(rt_band band) {
        assert(NULL != band);
 
        if (band->offline) {
                if (band->data.offline.mem != NULL)
                        return band->data.offline.mem;
 
                if (rt_band_load_offline_data(band) != ES_NONE)
                        return NULL;
                else
                        return band->data.offline.mem;
        }
        else
                return band->data.mem;
}
 
/* variable for PostgreSQL GUC: postgis.enable_outdb_rasters */
bool enable_outdb_rasters = true;
 
rt_errorstate
rt_band_load_offline_data(rt_band band) {
        GDALDatasetH hdsSrc = NULL;
        int nband = 0;
        VRTDatasetH hdsDst = NULL;
        VRTSourcedRasterBandH hbandDst = NULL;
        double ogt[6] = {0};
        double offset[2] = {0};
 
        rt_raster _rast = NULL;
        rt_band _band = NULL;
        int aligned = 0;
        int err = ES_NONE;
 
        assert(band != NULL);
        assert(band->raster != NULL);
 
        if (!band->offline) {
                rterror("rt_band_load_offline_data: Band is not offline");
                return ES_ERROR;
        }
        else if (!strlen(band->data.offline.path)) {
                rterror("rt_band_load_offline_data: Offline band does not a have a specified file");
                return ES_ERROR;
        }
 
        /* offline_data is disabled */
        if (!enable_outdb_rasters) {
                rterror("rt_band_load_offline_data: Access to offline bands disabled");
                return ES_ERROR;
        }
 
        rt_util_gdal_register_all(0);
        hdsSrc = rt_util_gdal_open(band->data.offline.path, GA_ReadOnly, 1);
        if (hdsSrc == NULL) {
                rterror("rt_band_load_offline_data: Cannot open offline raster: %s", band->data.offline.path);
                return ES_ERROR;
        }
 
        /* # of bands */
        nband = GDALGetRasterCount(hdsSrc);
        if (!nband) {
                rterror("rt_band_load_offline_data: No bands found in offline raster: %s", band->data.offline.path);
                GDALClose(hdsSrc);
                return ES_ERROR;
        }
        /* bandNum is 0-based */
        else if (band->data.offline.bandNum + 1 > nband) {
                rterror("rt_band_load_offline_data: Specified band %d not found in offline raster: %s", band->data.offline.bandNum, band->data.offline.path);
                GDALClose(hdsSrc);
                return ES_ERROR;
        }
 
        /* get offline raster's geotransform */
        if (GDALGetGeoTransform(hdsSrc, ogt) != CE_None) {
                RASTER_DEBUG(4, "Using default geotransform matrix (0, 1, 0, 0, 0, -1)");
                ogt[0] = 0;
                ogt[1] = 1;
                ogt[2] = 0;
                ogt[3] = 0;
                ogt[4] = 0;
                ogt[5] = -1;
        }
        RASTER_DEBUGF(3, "Offline geotransform (%f, %f, %f, %f, %f, %f)",
                ogt[0], ogt[1], ogt[2], ogt[3], ogt[4], ogt[5]);
 
        /* are rasters aligned? */
        _rast = rt_raster_new(1, 1);
        rt_raster_set_geotransform_matrix(_rast, ogt);
        rt_raster_set_srid(_rast, band->raster->srid);
        err = rt_raster_same_alignment(band->raster, _rast, &aligned, NULL);
        rt_raster_destroy(_rast);
 
        if (err != ES_NONE) {
                rterror("rt_band_load_offline_data: Could not test alignment of in-db representation of out-db raster");
                GDALClose(hdsSrc);
                return ES_ERROR;
        }
        else if (!aligned) {
                rtwarn("The in-db representation of the out-db raster is not aligned. Band data may be incorrect");
        }
 
        /* get offsets */
        rt_raster_geopoint_to_cell(
                band->raster,
                ogt[0], ogt[3],
                &(offset[0]), &(offset[1]),
                NULL
        );
 
        RASTER_DEBUGF(4, "offsets: (%f, %f)", offset[0], offset[1]);
 
        /* create VRT dataset */
        hdsDst = VRTCreate(band->width, band->height);
        GDALSetGeoTransform(hdsDst, ogt);
        /*
        GDALSetDescription(hdsDst, "/tmp/offline.vrt");
        */
 
        /* add band as simple sources */
        GDALAddBand(hdsDst, rt_util_pixtype_to_gdal_datatype(band->pixtype), NULL);
        hbandDst = (VRTSourcedRasterBandH) GDALGetRasterBand(hdsDst, 1);
 
        if (band->hasnodata)
                GDALSetRasterNoDataValue(hbandDst, band->nodataval);
 
        VRTAddSimpleSource(
                hbandDst, GDALGetRasterBand(hdsSrc, band->data.offline.bandNum + 1),
                fabs(offset[0]), fabs(offset[1]),
                band->width, band->height,
                0, 0,
                band->width, band->height,
                "near", VRT_NODATA_UNSET
        );
 
        /* make sure VRT reflects all changes */
        VRTFlushCache(hdsDst);
 
        /* convert VRT dataset to rt_raster */
        _rast = rt_raster_from_gdal_dataset(hdsDst);
 
        GDALClose(hdsDst);
        GDALClose(hdsSrc);
        /*
        {
                FILE *fp;
                fp = fopen("/tmp/gdal_open_files.log", "w");
                GDALDumpOpenDatasets(fp);
                fclose(fp);
        }
        */
 
        if (_rast == NULL) {
                rterror("rt_band_load_offline_data: Cannot load data from offline raster: %s", band->data.offline.path);
                return ES_ERROR;
        }
 
        _band = rt_raster_get_band(_rast, 0);
        if (_band == NULL) {
                rterror("rt_band_load_offline_data: Cannot load data from offline raster: %s", band->data.offline.path);
                rt_raster_destroy(_rast);
                return ES_ERROR;
        }
 
        /* band->data.offline.mem not NULL, free first */
        if (band->data.offline.mem != NULL) {
                rtdealloc(band->data.offline.mem);
                band->data.offline.mem = NULL;
        }
 
        band->data.offline.mem = _band->data.mem;
 
        rtdealloc(_band); /* cannot use rt_band_destroy */
        rt_raster_destroy(_rast);
 
        return ES_NONE;
}
 
uint64_t rt_band_get_file_size(rt_band band) {
    VSIStatBufL sStat;
 
    assert(NULL != band);
    if (!band->offline) {
        rterror("rt_band_get_file_size: Band is not offline");
        return 0;
    }
    /* offline_data is disabled */
    if (!enable_outdb_rasters) {
        rterror("rt_band_get_file_size: Access to offline bands disabled");
        return 0;
    }
 
    if( VSIStatL(band->data.offline.path, &sStat) != 0 ) {
        rterror("rt_band_get_file_size: Cannot access file");
        return 0;
    }
 
    return sStat.st_size;
}
 
uint64_t rt_band_get_file_timestamp(rt_band band) {
    VSIStatBufL sStat;
 
    assert(NULL != band);
    if (!band->offline) {
        rterror("rt_band_get_file_timestamp: Band is not offline");
        return 0;
    }
    /* offline_data is disabled */
    if (!enable_outdb_rasters) {
        rterror("rt_band_get_file_timestamp: Access to offline bands disabled");
        return 0;
    }
 
    if( VSIStatL(band->data.offline.path, &sStat) != 0 ) {
        rterror("rt_band_get_file_timestamp: Cannot access file");
        return 0;
    }
 
    return sStat.st_mtime;
}
 
rt_pixtype
rt_band_get_pixtype(rt_band band) {
 
    assert(NULL != band);
 
 
    return band->pixtype;
}
 
uint16_t
rt_band_get_width(rt_band band) {
 
    assert(NULL != band);
 
 
    return band->width;
}
 
uint16_t
rt_band_get_height(rt_band band) {
 
    assert(NULL != band);
 
 
    return band->height;
}
 
/* Get ownsdata flag */
int
rt_band_get_ownsdata_flag(rt_band band) {
        assert(NULL != band);
 
        return band->ownsdata ? 1 : 0;
}
 
/* set ownsdata flag */
void
rt_band_set_ownsdata_flag(rt_band band, int flag) {
        assert(NULL != band);
 
        band->ownsdata = flag ? 1 : 0;
}
 
int
rt_band_get_hasnodata_flag(rt_band band) {
        assert(NULL != band);
 
        return band->hasnodata ? 1 : 0;
}
 
void
rt_band_set_hasnodata_flag(rt_band band, int flag) {
 
    assert(NULL != band);
 
    band->hasnodata = (flag) ? 1 : 0;
 
                /* isnodata depends on hasnodata */
                if (!band->hasnodata && band->isnodata) {
                        RASTER_DEBUG(3, "Setting isnodata to FALSE as band no longer has NODATA");
                        band->isnodata = 0;
                }
}
 
rt_errorstate
rt_band_set_isnodata_flag(rt_band band, int flag) {
        assert(NULL != band);
 
        if (!band->hasnodata) {
                /* silently permit setting isnodata flag to FALSE */
                if (!flag)
                        band->isnodata = 0;
                else {
                        rterror("rt_band_set_isnodata_flag: Cannot set isnodata flag as band has no NODATA");
                        return ES_ERROR;
                }
        }
        else
                band->isnodata = (flag) ? 1 : 0;
 
        return ES_NONE;
}
 
int
rt_band_get_isnodata_flag(rt_band band) {
        assert(NULL != band);
 
        if (band->hasnodata)
                return band->isnodata ? 1 : 0;
        else
                return 0;
}
 
rt_errorstate
rt_band_set_nodata(rt_band band, double val, int *converted) {
        rt_pixtype pixtype = PT_END;
        int32_t checkvalint = 0;
        uint32_t checkvaluint = 0;
        float checkvalfloat = 0;
        double checkvaldouble = 0;
 
        assert(NULL != band);
 
        if (converted != NULL)
                *converted = 0;
 
        pixtype = band->pixtype;
 
        RASTER_DEBUGF(3, "rt_band_set_nodata: setting nodata value %g with band type %s", val, rt_pixtype_name(pixtype));
 
        /* return -1 on out of range */
        switch (pixtype) {
                case PT_1BB: {
                        band->nodataval = rt_util_clamp_to_1BB(val);
                        checkvalint = band->nodataval;
                        break;
                }
                case PT_2BUI: {
                        band->nodataval = rt_util_clamp_to_2BUI(val);
                        checkvalint = band->nodataval;
                        break;
                }
                case PT_4BUI: {
                        band->nodataval = rt_util_clamp_to_4BUI(val);
                        checkvalint = band->nodataval;
                        break;
                }
                case PT_8BSI: {
                        band->nodataval = rt_util_clamp_to_8BSI(val);
                        checkvalint = band->nodataval;
                        break;
                }
                case PT_8BUI: {
                        band->nodataval = rt_util_clamp_to_8BUI(val);
                        checkvalint = band->nodataval;
                        break;
                }
                case PT_16BSI: {
                        band->nodataval = rt_util_clamp_to_16BSI(val);
                        checkvalint = band->nodataval;
                        break;
                }
                case PT_16BUI: {
                        band->nodataval = rt_util_clamp_to_16BUI(val);
                        checkvalint = band->nodataval;
                        break;
                }
                case PT_32BSI: {
                        band->nodataval = rt_util_clamp_to_32BSI(val);
                        checkvalint = band->nodataval;
                        break;
                }
                case PT_32BUI: {
                        band->nodataval = rt_util_clamp_to_32BUI(val);
                        checkvaluint = band->nodataval;
                        break;
                }
                case PT_16BF: {
                        float half = rt_util_clamp_to_16F(val);
                        /* Persist the exact half-float encoding so nodata comparisons stay bitwise-stable. */
                        band->nodataval = rt_util_float16_to_float(rt_util_float_to_float16(half));
                        checkvalfloat = (float)band->nodataval;
                        break;
                }
                case PT_32BF: {
                        band->nodataval = rt_util_clamp_to_32F(val);
                        checkvalfloat = band->nodataval;
                        break;
                }
                case PT_64BF: {
                        band->nodataval = val;
                        checkvaldouble = band->nodataval;
                        break;
                }
                default: {
                        rterror("rt_band_set_nodata: Unknown pixeltype %d", pixtype);
                        band->hasnodata = 0;
                        return ES_ERROR;
                }
        }
 
        RASTER_DEBUGF(3, "rt_band_set_nodata: band->hasnodata = %d", band->hasnodata);
        RASTER_DEBUGF(3, "rt_band_set_nodata: band->nodataval = %f", band->nodataval);
        /* the nodata value was just set, so this band has NODATA */
        band->hasnodata = 1;
 
        /* also set isnodata flag to false */
        band->isnodata = 0;
 
        if (rt_util_dbl_trunc_warning(
                val,
                checkvalint, checkvaluint,
                checkvalfloat, checkvaldouble,
                pixtype
        ) && converted != NULL) {
                *converted = 1;
        }
 
        return ES_NONE;
}
 
rt_errorstate
rt_band_set_pixel_line(
        rt_band band,
        int x, int y,
        void *vals, uint32_t len
) {
        rt_pixtype pixtype = PT_END;
        int size = 0;
        uint8_t *data = NULL;
        uint32_t offset = 0;
 
        assert(NULL != band);
        assert(vals != NULL && len > 0);
 
        RASTER_DEBUGF(3, "length of values = %d", len);
 
        if (band->offline) {
                rterror("rt_band_set_pixel_line not implemented yet for OFFDB bands");
                return ES_ERROR;
        }
 
        pixtype = band->pixtype;
        size = rt_pixtype_size(pixtype);
 
        if (
                x < 0 || x >= band->width ||
                y < 0 || y >= band->height
        ) {
                rterror("rt_band_set_pixel_line: Coordinates out of range (%d, %d) vs (%d, %d)", x, y, band->width, band->height);
                return ES_ERROR;
        }
 
        data = rt_band_get_data(band);
        offset = x + (y * band->width);
        RASTER_DEBUGF(4, "offset = %d", offset);
 
        /* make sure len of values to copy don't exceed end of data */
        if (len > (band->width * band->height) - offset) {
                rterror("rt_band_set_pixel_line: Could not apply pixels as values length exceeds end of data");
                return ES_ERROR;
        }
 
        switch (pixtype) {
                case PT_1BB:
                case PT_2BUI:
                case PT_4BUI:
                case PT_8BUI:
                case PT_8BSI: {
                        uint8_t *ptr = data;
                        ptr += offset;
                        memcpy(ptr, vals, (size_t)size * len);
                        break;
                }
                case PT_16BUI: {
                        uint16_t *ptr = (uint16_t *) data;
                        ptr += offset;
                        memcpy(ptr, vals, (size_t)size * len);
                        break;
                }
                case PT_16BSI: {
                        int16_t *ptr = (int16_t *) data;
                        ptr += offset;
                        memcpy(ptr, vals, (size_t)size * len);
                        break;
                }
                case PT_32BUI: {
                        uint32_t *ptr = (uint32_t *) data;
                        ptr += offset;
                        memcpy(ptr, vals, (size_t)size * len);
                        break;
                }
                case PT_32BSI: {
                        int32_t *ptr = (int32_t *) data;
                        ptr += offset;
                        memcpy(ptr, vals, (size_t)size * len);
                        break;
                }
                case PT_32BF: {
                        float *ptr = (float *) data;
                        ptr += offset;
                        memcpy(ptr, vals, (size_t)size * len);
                        break;
                }
                case PT_64BF: {
                        double *ptr = (double *) data;
                        ptr += offset;
                        memcpy(ptr, vals, (size_t)size * len);
                        break;
                }
                default: {
                        rterror("rt_band_set_pixel_line: Unknown pixeltype %d", pixtype);
                        return ES_ERROR;
                }
        }
 
#if POSTGIS_DEBUG_LEVEL > 0
        {
                double value;
                rt_band_get_pixel(band, x, y, &value, NULL);
                RASTER_DEBUGF(4, "pixel at (%d, %d) = %f", x, y, value);
        }
#endif
 
        /* set band's isnodata flag to FALSE */
        if (rt_band_get_hasnodata_flag(band))
                rt_band_set_isnodata_flag(band, 0);
 
        return ES_NONE;
}
 
rt_errorstate
rt_band_set_pixel(
        rt_band band,
        int x, int y,
        double val,
        int *converted
) {
        rt_pixtype pixtype = PT_END;
        uint8_t *data = NULL;
        uint32_t offset = 0;
 
        int32_t checkvalint = 0;
        uint32_t checkvaluint = 0;
        float checkvalfloat = 0;
        double checkvaldouble = 0;
 
        assert(NULL != band);
 
        if (converted != NULL)
                *converted = 0;
 
        if (band->offline) {
                rterror("rt_band_set_pixel not implemented yet for OFFDB bands");
                return ES_ERROR;
        }
 
        pixtype = band->pixtype;
 
        if (
                x < 0 || x >= band->width ||
                y < 0 || y >= band->height
        ) {
                rterror("rt_band_set_pixel: Coordinates out of range");
                return ES_ERROR;
        }
 
        /* check that clamped value isn't clamped NODATA */
        if (band->hasnodata && pixtype != PT_64BF) {
                double newval;
                int corrected;
 
                rt_band_corrected_clamped_value(band, val, &newval, &corrected);
 
                if (corrected) {
#if POSTGIS_RASTER_WARN_ON_TRUNCATION > 0
                        rtwarn("Value for pixel %d x %d has been corrected as clamped value becomes NODATA", x, y);
#endif
                        val = newval;
 
                        if (converted != NULL)
                                *converted = 1;
                }
        }
 
        data = rt_band_get_data(band);
        offset = x + (y * band->width);
 
        switch (pixtype)
        {
        case PT_1BB: {
                data[offset] = rt_util_clamp_to_1BB(val);
                checkvalint = data[offset];
                break;
        }
                case PT_2BUI: {
                        data[offset] = rt_util_clamp_to_2BUI(val);
                        checkvalint = data[offset];
                        break;
                }
                case PT_4BUI: {
                        data[offset] = rt_util_clamp_to_4BUI(val);
                        checkvalint = data[offset];
                        break;
                }
                case PT_8BSI: {
                        data[offset] = (uint8_t)rt_util_clamp_to_8BSI(val);
                        checkvalint = (int8_t) data[offset];
                        break;
                }
                case PT_8BUI: {
                        data[offset] = rt_util_clamp_to_8BUI(val);
                        checkvalint = data[offset];
                        break;
                }
                case PT_16BSI: {
                        int16_t *ptr = (int16_t*) data; /* we assume correct alignment */
                        ptr[offset] = rt_util_clamp_to_16BSI(val);
                        checkvalint = (int16_t) ptr[offset];
                        break;
                }
                case PT_16BUI: {
                        uint16_t *ptr = (uint16_t*) data; /* we assume correct alignment */
                        ptr[offset] = rt_util_clamp_to_16BUI(val);
                        checkvalint = ptr[offset];
                        break;
                }
                case PT_32BSI: {
                        int32_t *ptr = (int32_t*) data; /* we assume correct alignment */
                        ptr[offset] = rt_util_clamp_to_32BSI(val);
                        checkvalint = (int32_t) ptr[offset];
                        break;
                }
                case PT_32BUI: {
                        uint32_t *ptr = (uint32_t *)data; /* we assume correct alignment */
                        ptr[offset] = rt_util_clamp_to_32BUI(val);
                        checkvaluint = ptr[offset];
                        break;
                }
                case PT_16BF: {
                        uint16_t *ptr = (uint16_t *)data; /* we assume correct alignment */
                        float clamped = rt_util_clamp_to_16F(val);
                        /* Pack via explicit converter to mirror GDAL's Float16 representation. */
                        ptr[offset] = rt_util_float_to_float16(clamped);
                        checkvalfloat = rt_util_float16_to_float(ptr[offset]);
                        break;
                }
                case PT_32BF: {
                        float *ptr = (float *)data; /* we assume correct alignment */
                        ptr[offset] = rt_util_clamp_to_32F(val);
                        checkvalfloat = ptr[offset];
                        break;
                }
                case PT_64BF: {
                        double *ptr = (double*) data; /* we assume correct alignment */
                        ptr[offset] = val;
                        checkvaldouble = ptr[offset];
                        break;
                }
                default: {
                        rterror("rt_band_set_pixel: Unknown pixeltype %d", pixtype);
                        return ES_ERROR;
                }
                }
 
        /* If the stored value is not NODATA, reset the isnodata flag */
        if (!rt_band_clamped_value_is_nodata(band, val)) {
                RASTER_DEBUG(3, "Band has a value that is not NODATA. Setting isnodata to FALSE");
                band->isnodata = FALSE;
        }
 
        /* Overflow checking */
        if (rt_util_dbl_trunc_warning(
                val,
                checkvalint, checkvaluint,
                checkvalfloat, checkvaldouble,
                pixtype
        ) && converted != NULL) {
                *converted = 1;
        }
 
        return ES_NONE;
}
 
rt_errorstate rt_band_get_pixel_line(
        rt_band band,
        int x, int y,
        uint16_t len,
        void **vals, uint16_t *nvals
) {
        uint8_t *_vals = NULL;
        int pixsize = 0;
        uint8_t *data = NULL;
        uint32_t offset = 0;
        uint16_t _nvals = 0;
        int maxlen = 0;
        uint8_t *ptr = NULL;
 
        assert(NULL != band);
        assert(vals != NULL && nvals != NULL);
 
        /* initialize to no values */
        *nvals = 0;
 
        if (
                x < 0 || x >= band->width ||
                y < 0 || y >= band->height
        ) {
                rtwarn("Attempting to get pixel values with out of range raster coordinates: (%d, %d)", x, y);
                return ES_ERROR;
        }
 
        if (len < 1)
                return ES_NONE;
 
        data = rt_band_get_data(band);
        if (data == NULL) {
                rterror("rt_band_get_pixel_line: Cannot get band data");
                return ES_ERROR;
        }
 
        /* +1 for the nodata value */
        offset = x + (y * band->width);
        RASTER_DEBUGF(4, "offset = %d", offset);
 
        pixsize = rt_pixtype_size(band->pixtype);
        RASTER_DEBUGF(4, "pixsize = %d", pixsize);
 
        /* cap _nvals so that it doesn't overflow */
        _nvals = len;
        maxlen = band->width * band->height;
 
        if (((int) (offset + _nvals)) > maxlen) {
                _nvals = maxlen - offset;
                rtwarn("Limiting returning number values to %d", _nvals);
        }
        RASTER_DEBUGF(4, "_nvals = %d", _nvals);
 
        ptr = data + ((size_t)offset * pixsize);
 
        _vals = rtalloc((size_t)_nvals * pixsize);
        if (_vals == NULL) {
                rterror("rt_band_get_pixel_line: Could not allocate memory for pixel values");
                return ES_ERROR;
        }
 
        /* copy pixels */
        memcpy(_vals, ptr, (size_t)_nvals * pixsize);
 
        *vals = _vals;
        *nvals = _nvals;
 
        return ES_NONE;
}
 
rt_errorstate
rt_band_get_pixel_resample(
        rt_band band,
        double xr, double yr,
        rt_resample_type resample,
        double *r_value, int *r_nodata
)
{
        if (resample == RT_BILINEAR) {
                return rt_band_get_pixel_bilinear(
                        band, xr, yr, r_value, r_nodata
                        );
        }
        else if (resample == RT_NEAREST) {
                return rt_band_get_pixel(
                        band, floor(xr), floor(yr),
                        r_value, r_nodata
                        );
        }
        else {
                rtwarn("Invalid resample type requested %d", resample);
                return ES_ERROR;
        }
 
}
 
rt_errorstate
rt_band_get_pixel_bilinear(
        rt_band band,
        double xr, double yr,
        double *r_value, int *r_nodata)
{
        double xcenter, ycenter;
        double values[2][2];
        double nodatavalue = 0.0;
        int   nodatas[2][2];
        int x[2][2];
        int y[2][2];
        int xcell, ycell;
        int xdir, ydir;
        int i, j;
        uint16_t width, height;
 
        /* Cell coordinates */
        xcell = (int)floor(xr);
        ycell = (int)floor(yr);
        xcenter = xcell + 0.5;
        ycenter = ycell + 0.5;
 
        /* Raster geometry */
        width = rt_band_get_width(band);
        height = rt_band_get_height(band);
 
        /* Reject out-of-range sample */
        if(xcell < 0 || ycell < 0 || xcell >= width || ycell >= height) {
                rtwarn("Attempting to get pixel value with out of range raster coordinates: (%d, %d)", xcell, ycell);
                return ES_ERROR;
        }
 
        /* Quadrant of 2x2 grid the raster coordinate falls in */
        xdir = xr < xcenter ? 1 : 0;
        ydir = yr < ycenter ? 1 : 0;
 
        if (rt_band_get_hasnodata_flag(band) != FALSE) {
                        rt_band_get_nodata(band, &nodatavalue);
        }
        else {
                        nodatavalue = 0.0;
        }
 
        /* Read the 2x2 values from the band */
        for (i = 0; i < 2; i++) {
                for (j = 0; j < 2; j++) {
                        double value = nodatavalue;
                        int nodata = 0;
                        int xij = xcell + (i - xdir);
                        int yij = ycell + (j - ydir);
 
                        if(xij < 0 || yij < 0 || xij >= width || yij >= height) {
                                nodata = 1;
                        }
                        else {
                                rt_errorstate err = rt_band_get_pixel(
                                        band, xij, yij,
                                        &value, &nodata
                                        );
                                if (err != ES_NONE)
                                        nodata = 1;
                        }
                        x[i][j] = xij;
                        y[i][j] = yij;
                        values[i][j] = value;
                        nodatas[i][j] = nodata;
                }
        }
 
        /* Point falls in nodata cell, just return nodata */
        if (nodatas[xdir][ydir]) {
                *r_value = nodatavalue;
                *r_nodata = 1;
                return ES_NONE;
        }
 
        /* Normalize raster coordinate to the bottom left */
        /* so we are working on a unit square */
        xr = xr - (x[0][0] + 0.5);
        yr = yr - (y[0][0] + 0.5);
 
        /* Point is in cell with values, so we take nodata */
        /* neighbors off the table by matching them to the */
        /* most controlling cell */
        for (i = 0; i < 2; i++) {
                for (j = 0; j < 2; j++) {
                        if (nodatas[i][j])
                                values[i][j] = values[xdir][ydir];
                }
        }
 
        /* Calculate bilinear value */
        /* https://en.wikipedia.org/wiki/Bilinear_interpolation#Unit_square */
        *r_nodata = 0;
        *r_value = values[0][0] * (1-xr) * (1-yr) +
                   values[1][0] * (1-yr) * xr +
                   values[0][1] * (1-xr) * yr +
                   values[1][1] * xr * yr;
 
        return ES_NONE;
}
 
 
rt_errorstate
rt_band_get_pixel(
        rt_band band,
        int x, int y,
        double *value,
        int *nodata
) {
        rt_pixtype pixtype = PT_END;
        uint8_t* data = NULL;
        uint32_t offset = 0;
 
        assert(NULL != band);
        assert(NULL != value);
 
        /* set nodata to 0 */
        if (nodata != NULL)
                *nodata = 0;
 
        if (
                x < 0 || x >= band->width ||
                y < 0 || y >= band->height
        ) {
                rtwarn("Attempting to get pixel value with out of range raster coordinates: (%d, %d)", x, y);
                return ES_ERROR;
        }
 
        /* band is NODATA */
        if (band->isnodata) {
                RASTER_DEBUG(3, "Band's isnodata flag is TRUE. Returning NODATA value");
                *value = band->nodataval;
                if (nodata != NULL) *nodata = 1;
                return ES_NONE;
        }
 
        data = rt_band_get_data(band);
        if (data == NULL) {
                rterror("rt_band_get_pixel: Cannot get band data");
                return ES_ERROR;
        }
 
        /* +1 for the nodata value */
        offset = x + (y * band->width);
 
        pixtype = band->pixtype;
 
        switch (pixtype) {
                case PT_1BB:
#ifdef OPTIMIZE_SPACE
                        {
                                int byteOffset = offset / 8;
                                int bitOffset = offset % 8;
                                data += byteOffset;
 
                                /* Bit to set is bitOffset into data */
                                *value = getBits(data, val, 1, bitOffset);
                                break;
                        }
#endif
                case PT_2BUI:
#ifdef OPTIMIZE_SPACE
                        {
                                int byteOffset = offset / 4;
                                int bitOffset = offset % 4;
                                data += byteOffset;
 
                                /* Bits to set start at bitOffset into data */
                                *value = getBits(data, val, 2, bitOffset);
                                break;
                        }
#endif
                case PT_4BUI:
#ifdef OPTIMIZE_SPACE
                        {
                                int byteOffset = offset / 2;
                                int bitOffset = offset % 2;
                                data += byteOffset;
 
                                /* Bits to set start at bitOffset into data */
                                *value = getBits(data, val, 2, bitOffset);
                                break;
                        }
#endif
                case PT_8BSI: {
                        int8_t val = (int8_t)data[offset];
                        *value = val;
                        break;
                }
                case PT_8BUI: {
                        uint8_t val = data[offset];
                        *value = val;
                        break;
                }
                case PT_16BSI: {
                        int16_t *ptr = (int16_t*) data; /* we assume correct alignment */
                        *value = ptr[offset];
                        break;
                }
                case PT_16BUI: {
                        uint16_t *ptr = (uint16_t*) data; /* we assume correct alignment */
                        *value = ptr[offset];
                        break;
                }
                case PT_32BSI: {
                        int32_t *ptr = (int32_t*) data; /* we assume correct alignment */
                        *value = ptr[offset];
                        break;
                }
                case PT_32BUI: {
                        uint32_t *ptr = (uint32_t *)data; /* we assume correct alignment */
                        *value = ptr[offset];
                        break;
                }
                case PT_16BF: {
                        uint16_t *ptr = (uint16_t *)data; /* we assume correct alignment */
                        *value = rt_util_float16_to_float(ptr[offset]);
                        break;
                }
                case PT_32BF: {
                        float *ptr = (float *)data; /* we assume correct alignment */
                        *value = ptr[offset];
                        break;
                }
                case PT_64BF: {
                        double *ptr = (double*) data; /* we assume correct alignment */
                        *value = ptr[offset];
                        break;
                }
                default: {
                        rterror("rt_band_get_pixel: Unknown pixeltype %d", pixtype);
                        return ES_ERROR;
                }
        }
 
        /* set NODATA flag */
        if (band->hasnodata && nodata != NULL) {
                if (rt_band_clamped_value_is_nodata(band, *value))
                        *nodata = 1;
        }
 
        return ES_NONE;
}
 
uint32_t rt_band_get_nearest_pixel(
        rt_band band,
        int x, int y,
        uint16_t distancex, uint16_t distancey,
        int exclude_nodata_value,
        rt_pixel *npixels
) {
        rt_pixel npixel = NULL;
        int extent[4] = {0};
        int max_extent[4] = {0};
        int d0 = 0;
        uint32_t distance[2] = {0};
        uint32_t _d[2] = {0};
        uint32_t i = 0;
        uint32_t j = 0;
        uint32_t k = 0;
        int _max = 0;
        int _x = 0;
        int _y = 0;
        int *_min = NULL;
        double pixval = 0;
        double minval = 0;
        uint32_t count = 0;
        int isnodata = 0;
 
        int inextent = 0;
 
        assert(NULL != band);
        assert(NULL != npixels);
 
        RASTER_DEBUG(3, "Starting");
 
        /* process distance */
        distance[0] = distancex;
        distance[1] = distancey;
 
        /* no distance, means get nearest pixels and return */
        if (!distance[0] && !distance[1])
                d0 = 1;
 
        RASTER_DEBUGF(4, "Selected pixel: %d x %d", x, y);
        RASTER_DEBUGF(4, "Distances: %d x %d", distance[0], distance[1]);
 
        /* shortcuts if outside band extent */
        if (
                exclude_nodata_value && (
                        (x < 0 || x > band->width) ||
                        (y < 0 || y > band->height)
                )
        ) {
                /* no distances specified, jump to pixel close to extent */
                if (d0) {
                        if (x < 0)
                                x = -1;
                        else if (x > band->width)
                                x = band->width;
 
                        if (y < 0)
                                y = -1;
                        else if (y > band->height)
                                y = band->height;
 
                        RASTER_DEBUGF(4, "Moved selected pixel: %d x %d", x, y);
                }
                /*
                        distances specified
                        if distances won't capture extent of band, return 0
                */
                else if (
                        ((x < 0 && (uint32_t) abs(x) > distance[0]) || (x - band->width >= (int)distance[0])) ||
                        ((y < 0 && (uint32_t) abs(y) > distance[1]) || (y - band->height >= (int)distance[1]))
                ) {
                        RASTER_DEBUG(4, "No nearest pixels possible for provided pixel and distances");
                        return 0;
                }
        }
 
        /* no NODATA, exclude is FALSE */
        if (!band->hasnodata)
                exclude_nodata_value = FALSE;
        /* band is NODATA and excluding NODATA */
        else if (exclude_nodata_value && band->isnodata) {
                RASTER_DEBUG(4, "No nearest pixels possible as band is NODATA and excluding NODATA values");
                return 0;
        }
 
        /* determine the maximum distance to prevent an infinite loop */
        if (d0) {
                int a, b;
 
                /* X axis */
                a = abs(x);
                b = abs(x - band->width);
 
                if (a > b)
                        distance[0] = a;
                else
                        distance[0] = b;
 
                /* Y axis */
                a = abs(y);
                b = abs(y - band->height);
                if (a > b)
                        distance[1] = a;
                else
                        distance[1] = b;
 
                RASTER_DEBUGF(4, "Maximum distances: %d x %d", distance[0], distance[1]);
        }
 
        /* minimum possible value for pixel type */
        minval = rt_pixtype_get_min_value(band->pixtype);
        RASTER_DEBUGF(4, "pixtype: %s", rt_pixtype_name(band->pixtype));
        RASTER_DEBUGF(4, "minval: %f", minval);
 
        /* set variables */
        count = 0;
        *npixels = NULL;
 
        /* maximum extent */
        max_extent[0] = x - (int)distance[0]; /* min X */
        max_extent[1] = y - (int)distance[1]; /* min Y */
        max_extent[2] = x + (int)distance[0]; /* max X */
        max_extent[3] = y + (int)distance[1]; /* max Y */
        RASTER_DEBUGF(4, "Maximum Extent: (%d, %d, %d, %d)",
                max_extent[0], max_extent[1], max_extent[2], max_extent[3]);
 
        _d[0] = 0;
        _d[1] = 0;
        do {
                _d[0]++;
                _d[1]++;
 
                extent[0] = x - (int)_d[0]; /* min x */
                extent[1] = y - (int)_d[1]; /* min y */
                extent[2] = x + (int)_d[0]; /* max x */
                extent[3] = y + (int)_d[1]; /* max y */
 
                RASTER_DEBUGF(4, "Processing distances: %d x %d", _d[0], _d[1]);
                RASTER_DEBUGF(4, "Extent: (%d, %d, %d, %d)",
                        extent[0], extent[1], extent[2], extent[3]);
 
                for (i = 0; i < 2; i++) {
 
                        /* by row */
                        if (i < 1)
                                _max = extent[2] - extent[0] + 1;
                        /* by column */
                        else
                                _max = extent[3] - extent[1] + 1;
                        _max = abs(_max);
 
                        for (j = 0; j < 2; j++) {
                                /* by row */
                                if (i < 1) {
                                        _x = extent[0];
                                        _min = &_x;
 
                                        /* top row */
                                        if (j < 1)
                                                _y = extent[1];
                                        /* bottom row */
                                        else
                                                _y = extent[3];
                                }
                                /* by column */
                                else {
                                        _y = extent[1] + 1;
                                        _min = &_y;
 
                                        /* left column */
                                        if (j < 1) {
                                                _x = extent[0];
                                                _max -= 2;
                                        }
                                        /* right column */
                                        else
                                                _x = extent[2];
                                }
 
                                RASTER_DEBUGF(4, "_min, _max: %d, %d", *_min, _max);
                                for (k = 0; k < (uint32_t) _max; k++) {
                                        /* check that _x and _y are not outside max extent */
                                        if (
                                                _x < max_extent[0] || _x > max_extent[2] ||
                                                _y < max_extent[1] || _y > max_extent[3]
                                        ) {
                                                (*_min)++;
                                                continue;
                                        }
 
                                        /* outside band extent, set to NODATA */
                                        if (
                                                (_x < 0 || _x >= band->width) ||
                                                (_y < 0 || _y >= band->height)
                                        ) {
                                                /* no NODATA, set to minimum possible value */
                                                if (!band->hasnodata)
                                                        pixval = minval;
                                                /* has NODATA, use NODATA */
                                                else
                                                        pixval = band->nodataval;
                                                RASTER_DEBUGF(4, "NODATA pixel outside band extent: (x, y, val) = (%d, %d, %f)", _x, _y, pixval);
                                                inextent = 0;
                                                isnodata = 1;
                                        }
                                        else {
                                                if (rt_band_get_pixel(
                                                        band,
                                                        _x, _y,
                                                        &pixval,
                                                        &isnodata
                                                ) != ES_NONE) {
                                                        rterror("rt_band_get_nearest_pixel: Could not get pixel value");
                                                        if (count) rtdealloc(*npixels);
                                                        return -1;
                                                }
                                                RASTER_DEBUGF(4, "Pixel: (x, y, val) = (%d, %d, %f)", _x, _y, pixval);
                                                inextent = 1;
                                        }
 
                                        /* use pixval? */
                                        if (!exclude_nodata_value || (exclude_nodata_value && !isnodata)) {
                                                /* add pixel to result set */
                                                RASTER_DEBUGF(4, "Adding pixel to set of nearest pixels: (x, y, val) = (%d, %d, %f)", _x, _y, pixval);
                                                count++;
 
                                                if (*npixels == NULL)
                                                        *npixels = (rt_pixel) rtalloc(sizeof(struct rt_pixel_t) * count);
                                                else
                                                        *npixels = (rt_pixel) rtrealloc(*npixels, sizeof(struct rt_pixel_t) * count);
                                                if (*npixels == NULL) {
                                                        rterror("rt_band_get_nearest_pixel: Could not allocate memory for nearest pixel(s)");
                                                        return -1;
                                                }
 
                                                npixel = &((*npixels)[count - 1]);
                                                npixel->x = _x;
                                                npixel->y = _y;
                                                npixel->value = pixval;
 
                                                /* special case for when outside band extent */
                                                if (!inextent && !band->hasnodata)
                                                        npixel->nodata = 1;
                                                else
                                                        npixel->nodata = 0;
                                        }
 
                                        (*_min)++;
                                }
                        }
                }
 
                /* distance thresholds met */
                if (_d[0] >= distance[0] && _d[1] >= distance[1])
                        break;
                else if (d0 && count)
                        break;
        }
        while (1);
 
        RASTER_DEBUGF(3, "Nearest pixels in return: %d", count);
 
        return count;
}
 
 
 
int
rt_band_get_pixel_of_value(
        rt_band band, int exclude_nodata_value,
        double *searchset, int searchcount,
        rt_pixel *pixels
) {
        int x;
        int y;
        int i;
        double pixval;
        int err;
        int count = 0;
        int isnodata = 0;
        int isequal = 0;
 
        rt_pixel pixel = NULL;
 
        assert(NULL != band);
        assert(NULL != pixels);
        assert(NULL != searchset && searchcount > 0);
 
        if (!band->hasnodata)
                exclude_nodata_value = FALSE;
        /* band is NODATA and exclude_nodata_value = TRUE, nothing to search */
        else if (exclude_nodata_value && band->isnodata) {
                RASTER_DEBUG(4, "Pixels cannot be searched as band is NODATA and excluding NODATA values");
                return 0;
        }
 
        for (x = 0; x < band->width; x++) {
                for (y = 0; y < band->height; y++) {
                        err = rt_band_get_pixel(band, x, y, &pixval, &isnodata);
                        if (err != ES_NONE) {
                                rterror("rt_band_get_pixel_of_value: Cannot get band pixel");
                                return -1;
                        }
                        else if (exclude_nodata_value && isnodata)
                                continue;
 
                        for (i = 0; i < searchcount; i++) {
                                if (rt_pixtype_compare_clamped_values(band->pixtype, searchset[i], pixval, &isequal) != ES_NONE) {
                                        continue;
                                }
 
                                if (FLT_NEQ(pixval, searchset[i]) || !isequal)
                                        continue;
 
                                /* match found */
                                count++;
                                if (*pixels == NULL)
                                        *pixels = (rt_pixel) rtalloc(sizeof(struct rt_pixel_t) * count);
                                else
                                        *pixels = (rt_pixel) rtrealloc(*pixels, sizeof(struct rt_pixel_t) * count);
                                if (*pixels == NULL) {
                                        rterror("rt_band_get_pixel_of_value: Could not allocate memory for pixel(s)");
                                        return -1;
                                }
 
                                pixel = &((*pixels)[count - 1]);
                                pixel->x = x;
                                pixel->y = y;
                                pixel->nodata = 0;
                                pixel->value = pixval;
                        }
                }
        }
 
        return count;
}
 
rt_errorstate
rt_band_get_nodata(rt_band band, double *nodata) {
        assert(NULL != band);
        assert(NULL != nodata);
 
        *nodata = band->nodataval;
 
        if (!band->hasnodata) {
                rterror("rt_band_get_nodata: Band has no NODATA value");
                return ES_ERROR;
        }
 
        return ES_NONE;
}
 
double
rt_band_get_min_value(rt_band band) {
        assert(NULL != band);
 
        return rt_pixtype_get_min_value(band->pixtype);
}
 
int
rt_band_check_is_nodata(rt_band band) {
        int i, j, err;
        double pxValue;
        int isnodata = 0;
 
        assert(NULL != band);
        band->isnodata = FALSE;
 
        /* Check if band has nodata value */
        if (!band->hasnodata) {
                RASTER_DEBUG(3, "Band has no NODATA value");
                return FALSE;
        }
 
        pxValue = band->nodataval;
 
        /* Check all pixels */
        for (i = 0; i < band->width; i++) {
                for (j = 0; j < band->height; j++) {
                        err = rt_band_get_pixel(band, i, j, &pxValue, &isnodata);
                        if (err != ES_NONE) {
                                rterror("rt_band_check_is_nodata: Cannot get band pixel");
                                return FALSE;
                        }
                        else if (!isnodata) {
                                band->isnodata = FALSE;
                                return FALSE;
                        }
                }
        }
 
        band->isnodata = TRUE;
        return TRUE;
}
 
int
rt_band_clamped_value_is_nodata(rt_band band, double val) {
        int isequal = 0;
 
        assert(NULL != band);
 
        /* no NODATA, so never equal */
        if (!band->hasnodata)
                return 0;
 
        /* value is exactly NODATA */
        if (FLT_EQ(val, band->nodataval))
                return 2;
 
        /* ignore error from rt_pixtype_compare_clamped_values */
        rt_pixtype_compare_clamped_values(
                band->pixtype,
                val, band->nodataval,
                &isequal
        );
 
        return isequal ? 1 : 0;
}
 
rt_errorstate
rt_band_corrected_clamped_value(
        rt_band band,
        double val,
        double *newval, int *corrected
) {
        double minval = 0.;
 
        assert(NULL != band);
        assert(NULL != newval);
 
        if (corrected != NULL)
                *corrected = 0;
 
        /* no need to correct if clamped values IS NOT clamped NODATA */
        if (rt_band_clamped_value_is_nodata(band, val) != 1) {
                *newval = val;
                return ES_NONE;
        }
 
        minval = rt_pixtype_get_min_value(band->pixtype);
        *newval = val;
 
        switch (band->pixtype) {
                case PT_1BB:
                        *newval = !band->nodataval;
                        break;
                case PT_2BUI:
                        if (rt_util_clamp_to_2BUI(val) == rt_util_clamp_to_2BUI(minval))
                                (*newval)++;
                        else
                                (*newval)--;
                        break;
                case PT_4BUI:
                        if (rt_util_clamp_to_4BUI(val) == rt_util_clamp_to_4BUI(minval))
                                (*newval)++;
                        else
                                (*newval)--;
                        break;
                case PT_8BSI:
                        if (rt_util_clamp_to_8BSI(val) == rt_util_clamp_to_8BSI(minval))
                                (*newval)++;
                        else
                                (*newval)--;
                        break;
                case PT_8BUI:
                        if (rt_util_clamp_to_8BUI(val) == rt_util_clamp_to_8BUI(minval))
                                (*newval)++;
                        else
                                (*newval)--;
                        break;
                case PT_16BSI:
                        if (rt_util_clamp_to_16BSI(val) == rt_util_clamp_to_16BSI(minval))
                                (*newval)++;
                        else
                                (*newval)--;
                        break;
                case PT_16BUI:
                        if (rt_util_clamp_to_16BUI(val) == rt_util_clamp_to_16BUI(minval))
                                (*newval)++;
                        else
                                (*newval)--;
                        break;
                case PT_32BSI:
                        if (rt_util_clamp_to_32BSI(val) == rt_util_clamp_to_32BSI(minval))
                                (*newval)++;
                        else
                                (*newval)--;
                        break;
                case PT_32BUI:
                        if (rt_util_clamp_to_32BUI(val) == rt_util_clamp_to_32BUI(minval))
                                (*newval)++;
                        else
                                (*newval)--;
                        break;
                case PT_16BF:
                        if (FLT_EQ(rt_util_clamp_to_16F(val), rt_util_clamp_to_16F(minval)))
                                *newval += FLT_EPSILON;
                        else
                                *newval -= FLT_EPSILON;
                        break;
                case PT_32BF:
                        if (FLT_EQ(rt_util_clamp_to_32F(val), rt_util_clamp_to_32F(minval)))
                                *newval += FLT_EPSILON;
                        else
                                *newval -= FLT_EPSILON;
                        break;
                case PT_64BF:
                        break;
                default:
                        rterror("rt_band_corrected_clamped_value: Unknown pixeltype %d", band->pixtype);
                        return ES_ERROR;
        }
 
        if (corrected != NULL)
                *corrected = 1;
 
        return ES_NONE;
}