Add histogram calculation for RasterBand

CHANGES.md

@@ -45,6 +45,8 @@
 
   - <https://github.com/georust/gdal/pull/439>
 
+- Added Histogram calculation
+
 ## 0.16
 
 - **Breaking**: `Dataset::close` now consumes `self`

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};

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<Option<Histogram>> {
+        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<Option<Histogram>> {
+        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<Vec<i32>>,
+}
+
+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

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::<Vec<_>>();
+    assert_eq!(hist.histogram(), &expected_small);
+}
+
 #[test]
 fn test_resample_str() {
     assert!(ResampleAlg::from_str("foobar").is_err());