diff --git a/CHANGES.md b/CHANGES.md index b419ffb5..1c46bc39 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -45,6 +45,8 @@ - +- Added Histogram calculation + ## 0.16 - **Breaking**: `Dataset::close` now consumes `self` diff --git a/src/raster/mod.rs b/src/raster/mod.rs index 1e266b90..a4ef8ca1 100644 --- a/src/raster/mod.rs +++ b/src/raster/mod.rs @@ -87,7 +87,7 @@ pub use mdarray::{ }; pub use rasterband::{ Buffer, ByteBuffer, CmykEntry, ColorEntry, ColorInterpretation, ColorTable, GrayEntry, - HlsEntry, PaletteInterpretation, RasterBand, ResampleAlg, RgbaEntry, StatisticsAll, + Histogram, HlsEntry, PaletteInterpretation, RasterBand, ResampleAlg, RgbaEntry, StatisticsAll, StatisticsMinMax, }; pub use rasterize::{rasterize, BurnSource, MergeAlgorithm, OptimizeMode, RasterizeOptions}; diff --git a/src/raster/rasterband.rs b/src/raster/rasterband.rs index dbb540df..8328a2b2 100644 --- a/src/raster/rasterband.rs +++ b/src/raster/rasterband.rs @@ -5,9 +5,10 @@ use crate::raster::{GdalDataType, GdalType}; use crate::utils::{_last_cpl_err, _last_null_pointer_err, _string}; use gdal_sys::{ self, CPLErr, GDALColorEntry, GDALColorInterp, GDALColorTableH, GDALComputeRasterMinMax, - GDALCreateColorRamp, GDALCreateColorTable, GDALDestroyColorTable, GDALGetPaletteInterpretation, - GDALGetRasterStatistics, GDALMajorObjectH, GDALPaletteInterp, GDALRIOResampleAlg, GDALRWFlag, - GDALRasterBandH, GDALRasterIOExtraArg, GDALSetColorEntry, GDALSetRasterColorTable, + GDALCreateColorRamp, GDALCreateColorTable, GDALDestroyColorTable, GDALGetDefaultHistogram, + GDALGetPaletteInterpretation, GDALGetRasterStatistics, GDALMajorObjectH, GDALPaletteInterp, + GDALRIOResampleAlg, GDALRWFlag, GDALRasterBandH, GDALRasterIOExtraArg, GDALSetColorEntry, + GDALSetRasterColorTable, }; use libc::c_int; use std::ffi::CString; @@ -847,6 +848,96 @@ impl<'a> RasterBand<'a> { max: min_max[1], }) } + + /// Fetch default raster histogram. + /// + /// # Arguments + /// + /// * `force` - If `true`, force the computation. If `false` and no default histogram is available, the method will return None + pub fn get_default_histogram(&self, force: bool) -> Result> { + let mut hist = Histogram { + min: 0.0, + max: 0.0, + n_buckets: 0, + histogram: std::ptr::null_mut(), + histogram_vec: None, + }; + + let rv = unsafe { + GDALGetDefaultHistogram( + self.c_rasterband, + &mut hist.min, + &mut hist.max, + &mut hist.n_buckets, + &mut hist.histogram as *mut *mut i32, + libc::c_int::from(force), + None, + std::ptr::null_mut(), + ) + }; + + match CplErrType::from(rv) { + CplErrType::None => { + // if everything when ok, can this still be null? + if hist.histogram.is_null() { + Ok(None) + } else { + Ok(Some(hist)) + } + } + CplErrType::Warning => Ok(None), + _ => Err(_last_cpl_err(rv)), + } + } + + /// Compute raster histogram. + /// + /// # Arguments + /// + /// * `min` - Histogram lower bound + /// * `max` - Histogram upper bound + /// * `n_buckets` - Number of buckets in the histogram + /// * `include_out_of_range` - if `true`, values below the histogram range will be mapped into the first bucket, and values above will be mapped into the last one. If `false`, out of range values are discarded + /// * `is_approx_ok` - If an approximate, or incomplete histogram is OK + pub fn get_histogram( + &self, + min: f64, + max: f64, + n_buckets: i32, + include_out_of_range: bool, + is_approx_ok: bool, + ) -> Result> { + let mut hist = Histogram { + min, + max, + n_buckets, + histogram_vec: Some(vec![0; n_buckets as usize]), + histogram: std::ptr::null_mut(), + }; + + let rv = unsafe { + gdal_sys::GDALGetRasterHistogram( + self.c_rasterband, + min, + max, + n_buckets, + hist.histogram_vec.as_mut().unwrap().as_mut_ptr(), + libc::c_int::from(include_out_of_range), + libc::c_int::from(is_approx_ok), + None, + std::ptr::null_mut(), + ) + }; + + match CplErrType::from(rv) { + CplErrType::None => { + hist.histogram = hist.histogram_vec.as_mut().unwrap().as_mut_ptr(); + Ok(Some(hist)) + } + CplErrType::Warning => Ok(None), + _ => Err(_last_cpl_err(rv)), + } + } } #[derive(Debug, PartialEq)] @@ -863,6 +954,50 @@ pub struct StatisticsAll { pub std_dev: f64, } +#[derive(Debug)] +pub struct Histogram { + min: f64, + max: f64, + n_buckets: i32, + histogram: *mut i32, + histogram_vec: Option>, +} + +impl Histogram { + /// Histogram lower bound + pub fn min(&self) -> f64 { + self.min + } + + /// Histogram upper bound + pub fn max(&self) -> f64 { + self.max + } + /// Histogram values for each bucket + pub fn histogram(&self) -> &[i32] { + if let Some(hist) = self.histogram_vec.as_ref() { + return hist.as_slice(); + } + + unsafe { std::slice::from_raw_parts(self.histogram, self.n_buckets as usize) } + } + + /// Histogram bucket size + pub fn bucket_size(&self) -> f64 { + (self.max - self.min) / self.histogram().len() as f64 + } +} + +impl Drop for Histogram { + fn drop(&mut self) { + if self.histogram_vec.is_none() { + unsafe { + gdal_sys::VSIFree(self.histogram as *mut libc::c_void); + } + } + } +} + impl<'a> MajorObject for RasterBand<'a> { fn gdal_object_ptr(&self) -> GDALMajorObjectH { self.c_rasterband diff --git a/src/raster/tests.rs b/src/raster/tests.rs index 3c97665a..857a5220 100644 --- a/src/raster/tests.rs +++ b/src/raster/tests.rs @@ -751,6 +751,39 @@ fn test_raster_stats() { ); } +#[test] +fn test_raster_histogram() { + let fixture = TempFixture::fixture("tinymarble.tif"); + + let dataset = Dataset::open(&fixture).unwrap(); + let rb = dataset.rasterband(1).unwrap(); + + let hist = rb.get_default_histogram(true).unwrap().unwrap(); + let expected = &[ + 548, 104, 133, 127, 141, 125, 156, 129, 130, 117, 94, 94, 80, 81, 78, 63, 50, 66, 48, 48, + 33, 38, 41, 35, 41, 39, 32, 40, 26, 27, 25, 24, 18, 25, 29, 27, 20, 34, 17, 24, 29, 11, 20, + 21, 12, 19, 16, 16, 11, 10, 19, 5, 11, 10, 6, 9, 7, 12, 13, 6, 8, 7, 8, 14, 9, 14, 4, 8, 5, + 12, 6, 10, 7, 9, 8, 6, 3, 7, 5, 8, 9, 5, 4, 8, 3, 9, 3, 6, 11, 7, 6, 3, 9, 9, 7, 6, 9, 10, + 10, 4, 7, 2, 4, 7, 2, 12, 7, 10, 4, 6, 5, 2, 4, 5, 7, 3, 5, 7, 7, 14, 9, 12, 6, 6, 8, 5, 8, + 3, 3, 5, 11, 4, 9, 7, 14, 7, 10, 11, 6, 6, 5, 4, 9, 6, 6, 9, 5, 12, 11, 9, 3, 8, 5, 6, 4, + 2, 9, 7, 9, 9, 9, 6, 6, 8, 5, 9, 13, 4, 9, 4, 7, 13, 10, 5, 7, 8, 11, 12, 5, 17, 9, 11, 9, + 8, 9, 5, 8, 9, 5, 6, 9, 11, 8, 7, 7, 6, 7, 8, 8, 8, 5, 6, 7, 5, 8, 5, 6, 8, 7, 4, 8, 6, 5, + 11, 8, 8, 5, 4, 6, 4, 9, 7, 6, 6, 7, 7, 12, 6, 9, 17, 12, 20, 18, 17, 21, 24, 30, 29, 57, + 72, 83, 21, 11, 9, 18, 7, 13, 10, 2, 4, 0, 1, 3, 4, 1, 1, + ]; + assert_eq!(hist.histogram(), expected); + + let hist = rb + .get_histogram(-0.5, 255.5, 128, true, true) + .unwrap() + .unwrap(); + let expected_small = (0..expected.len()) + .step_by(2) + .map(|i| expected[i] + expected[i + 1]) + .collect::>(); + assert_eq!(hist.histogram(), &expected_small); +} + #[test] fn test_resample_str() { assert!(ResampleAlg::from_str("foobar").is_err());