geo_rasterize/
lib.rs

1#![doc = include_str!("../README.md")]
2use std::{collections::HashSet, fmt::Debug, ops::Add};
3
4use euclid::{Transform2D, UnknownUnit};
5use geo::{
6    algorithm::{
7        coords_iter::CoordsIter,
8        map_coords::{MapCoords, MapCoordsInplace},
9    },
10    Geometry, GeometryCollection, Line, LineString, MultiLineString, MultiPoint, MultiPolygon,
11    Point, Polygon, Rect, Triangle,
12};
13use ndarray::s;
14use ndarray::Array2;
15use num_traits::{Num, NumCast};
16use thiserror::Error;
17
18mod poly;
19use poly::rasterize_polygon;
20mod line;
21use line::rasterize_line;
22#[cfg(test)]
23mod proptests;
24
25/// Affine transform that describes how to convert world-space
26/// coordinates to pixel coordinates.
27pub type Transform = Transform2D<f64, UnknownUnit, UnknownUnit>;
28type EuclidPoint = euclid::Point2D<f64, UnknownUnit>;
29
30/// Error type for this crate
31#[derive(Error, Clone, Debug, PartialEq, Eq)]
32pub enum RasterizeError {
33    /// at least one coordinate of the supplied geometry is NaN or infinite
34    #[error("at least one coordinate of the supplied geometry is NaN or infinite")]
35    NonFiniteCoordinate,
36
37    /// `width` is required in builder
38    #[error("`width` is required in builder")]
39    MissingWidth,
40
41    /// `height` is required in builder
42    #[error("`height` is required in builder")]
43    MissingHeight,
44}
45
46/// Result type for this crate that uses [RasterizeError].
47pub type Result<T> = std::result::Result<T, RasterizeError>;
48
49/// A builder that can construct instances of [BinaryRasterizer], a
50/// rasterizer that can rasterize shapes into a 2-dimensional array of
51/// booleans.
52///
53/// ```rust
54/// # use geo_rasterize::{Result, BinaryBuilder, BinaryRasterizer};
55/// # fn main() -> Result<()> {
56/// let rasterizer: BinaryRasterizer = BinaryBuilder::new().width(37).height(21).build()?;
57/// # Ok(())}
58/// ```
59#[derive(Debug, Clone, Default)]
60pub struct BinaryBuilder {
61    width: Option<usize>,
62    height: Option<usize>,
63    geo_to_pix: Option<Transform>,
64}
65
66impl BinaryBuilder {
67    pub fn new() -> Self {
68        BinaryBuilder::default()
69    }
70
71    pub fn width(mut self, width: usize) -> Self {
72        self.width = Some(width);
73        self
74    }
75
76    pub fn height(mut self, height: usize) -> Self {
77        self.height = Some(height);
78        self
79    }
80
81    pub fn geo_to_pix(mut self, geo_to_pix: Transform) -> Self {
82        self.geo_to_pix = Some(geo_to_pix);
83        self
84    }
85
86    pub fn build(self) -> Result<BinaryRasterizer> {
87        match (self.width, self.height) {
88            (None, _) => Err(RasterizeError::MissingWidth),
89            (_, None) => Err(RasterizeError::MissingHeight),
90            (Some(width), Some(height)) => BinaryRasterizer::new(width, height, self.geo_to_pix),
91        }
92    }
93}
94
95/// A rasterizer that burns shapes into a 2-dimensional boolean
96/// array. It can be built either by calling
97/// [new][BinaryRasterizer::new] or using [BinaryBuilder].
98///
99/// Each Rasterizer requires a `width` and `height` measured in pixels
100/// that describe the shape of the output array. They can optionally
101/// take an affine transform that describes how to convert world-space
102/// coordinates into pixel-space coordinates. When a transformer is
103/// supplied, all of its parameters must be finite or
104/// [RasterizeError::NonFiniteCoordinate] will be returned.
105///
106/// ```rust
107/// # use geo_rasterize::{Result, BinaryBuilder, BinaryRasterizer};
108/// # fn main() -> Result<()> {
109/// use geo::{Geometry, Line, Point};
110/// use ndarray::array;
111/// use geo_rasterize::BinaryBuilder;
112///
113/// let shapes: Vec<Geometry<i32>> =
114///     vec![Point::new(3, 4).into(),
115///          Line::new((0, 3), (3, 0)).into()];
116///
117/// let mut r = BinaryBuilder::new().width(4).height(5).build()?;
118/// for shape in shapes {
119///     r.rasterize(&shape)?;
120/// }
121///
122/// let pixels = r.finish();
123/// assert_eq!(
124///     pixels.mapv(|v| v as u8),
125///     array![
126///         [0, 0, 1, 0],
127///         [0, 1, 1, 0],
128///         [1, 1, 0, 0],
129///         [1, 0, 0, 0],
130///         [0, 0, 0, 1]
131///     ]
132/// );
133/// # Ok(())}
134/// ```
135#[derive(Clone, Debug)]
136pub struct BinaryRasterizer {
137    inner: Rasterizer<u8>,
138}
139
140fn to_float<T>(coords: &(T, T)) -> (f64, f64)
141where
142    T: Into<f64> + Copy,
143{
144    (coords.0.into(), coords.1.into())
145}
146
147impl BinaryRasterizer {
148    pub fn new(width: usize, height: usize, geo_to_pix: Option<Transform>) -> Result<Self> {
149        let non_finite = geo_to_pix
150            .map(|geo_to_pix| geo_to_pix.to_array().iter().any(|param| !param.is_finite()))
151            .unwrap_or(false);
152        if non_finite {
153            Err(RasterizeError::NonFiniteCoordinate)
154        } else {
155            let inner = Rasterizer::new(width, height, geo_to_pix, MergeAlgorithm::Replace, 0);
156            Ok(BinaryRasterizer { inner })
157        }
158    }
159
160    /// Retrieve the transform.
161    pub fn geo_to_pix(&self) -> Option<Transform> {
162        self.inner.geo_to_pix
163    }
164
165    /// Rasterize one shape, which can be any type that [geo] provides
166    /// using any coordinate numeric type that can be converted into
167    /// `f64`.
168    pub fn rasterize<Coord, InputShape, ShapeAsF64>(&mut self, shape: &InputShape) -> Result<()>
169    where
170        InputShape: MapCoords<Coord, f64, Output = ShapeAsF64>,
171        ShapeAsF64: Rasterize<u8> + for<'a> CoordsIter<'a, Scalar = f64> + MapCoordsInplace<f64>,
172        Coord: Into<f64> + Copy + Debug + Num + NumCast + PartialOrd,
173    {
174        // first, convert our input shape so that its coordinates are of type f64
175        self.inner.rasterize(shape, 1)
176    }
177
178    /// Retrieve the completed raster array.
179    pub fn finish(self) -> Array2<bool> {
180        self.inner
181            .finish()
182            .mapv(|v| if v == 1u8 { true } else { false })
183    }
184}
185
186#[doc(hidden)]
187pub trait Rasterize<Label>
188where
189    Label: Copy + Add<Output = Label>,
190{
191    fn rasterize(&self, rasterizer: &mut Rasterizer<Label>);
192}
193
194/// Conflict resolution strategy for cases where two shapes cover the
195/// same pixel.
196#[derive(Debug, Clone, Copy, PartialEq, Eq)]
197pub enum MergeAlgorithm {
198    /// Overwrite the pixel with the burn value associated with the
199    /// last shape to be written to it. This is the default.
200    Replace,
201
202    /// Overwrite the pixel with the sum of the burn values associated
203    /// with the shapes written to it.
204    Add,
205}
206
207impl Default for MergeAlgorithm {
208    fn default() -> Self {
209        MergeAlgorithm::Replace
210    }
211}
212
213/// A builder that constructs [Rasterizer]s. Whereas
214/// [BinaryRasterizer] produces an array of booleans, [Rasterizer]
215/// produces an array of some generic type (`Label`) that implements
216/// [Copy][std::marker::Copy] and [Add][std::ops::Add] though
217/// typically you'd use a numeric type.
218///
219/// [LabelBuilder] needs the `Label` type so the only way to make one
220/// is to specify a `Label` value: the background. The `background` is
221/// the value we'll use to initialize the raster array -- it
222/// corresponds to the zero value, so you'll probably want to start
223/// with `LabelBuilder::background(0)` or
224/// `LabelBuilder::background(0f32)`.
225///
226/// In addition to `background`, `width`, and `height`, you can also
227/// supply a [MergeAlgorithm] to specify what the rasterizer should do
228/// when two different shapes fill the same pixel. If you don't supply
229/// anything, the rasterizer will use
230/// [Replace][MergeAlgorithm::Replace] by default.
231///
232/// ```rust
233/// # fn main() -> geo_rasterize::Result<()> {
234/// use geo_rasterize::LabelBuilder;
235///
236/// let mut rasterizer = LabelBuilder::background(0i32).width(4).height(5).build()?;
237/// # Ok(())}
238/// ```
239#[derive(Debug, Clone, Default)]
240pub struct LabelBuilder<Label> {
241    background: Label,
242    width: Option<usize>,
243    height: Option<usize>,
244    geo_to_pix: Option<Transform>,
245    algorithm: Option<MergeAlgorithm>,
246}
247
248impl<Label> LabelBuilder<Label>
249where
250    Label: Copy + Add<Output = Label>,
251{
252    pub fn background(background: Label) -> Self {
253        LabelBuilder {
254            background,
255            width: None,
256            height: None,
257            geo_to_pix: None,
258            algorithm: None,
259        }
260    }
261
262    pub fn width(mut self, width: usize) -> Self {
263        self.width = Some(width);
264        self
265    }
266
267    pub fn height(mut self, height: usize) -> Self {
268        self.height = Some(height);
269        self
270    }
271
272    pub fn geo_to_pix(mut self, geo_to_pix: Transform) -> Self {
273        self.geo_to_pix = Some(geo_to_pix);
274        self
275    }
276
277    pub fn algorithm(mut self, algorithm: MergeAlgorithm) -> Self {
278        self.algorithm = Some(algorithm);
279        self
280    }
281
282    pub fn build(self) -> Result<Rasterizer<Label>> {
283        match (self.width, self.height) {
284            (None, _) => Err(RasterizeError::MissingWidth),
285            (_, None) => Err(RasterizeError::MissingHeight),
286            (Some(width), Some(height)) => Ok(Rasterizer::new(
287                width,
288                height,
289                self.geo_to_pix,
290                self.algorithm.unwrap_or_default(),
291                self.background,
292            )),
293        }
294    }
295}
296
297/// [Rasterizer] rasterizes shapes like [BinaryRasterizer] but instead
298/// of making a boolean array, it produces an array of some generic
299/// type (`Label`) that implements [Copy][std::marker::Copy] and
300/// [Add][std::ops::Add] though typically you'd use a numeric type.
301///
302/// You can call [new][Rasterizer::new] or use [LabelBuilder] to
303/// construct [Rasterizer] instances. Constructing one requires a
304/// `width`, `height` and optional `geo_to_pix` transform in addition
305/// to `background` which specifies the default `Label` value used to
306/// fill the raster array. And you can provide a [MergeAlgorithm]
307/// value to specify what the rasterizer should do when two different
308/// shapes fill the same pixel. If you don't supply anything, the
309/// rasterizer will use [Replace][MergeAlgorithm::Replace] by default.
310///
311/// ```rust
312/// # use geo_rasterize::{Result, LabelBuilder, Rasterizer};
313/// # fn main() -> Result<()> {
314/// use geo::{Geometry, Line, Point};
315/// use ndarray::array;
316///
317/// let point = Point::new(3, 4);
318/// let line = Line::new((0, 3), (3, 0));
319///
320/// let mut rasterizer = LabelBuilder::background(0).width(4).height(5).build()?;
321/// rasterizer.rasterize(&point, 7)?;
322/// rasterizer.rasterize(&line, 3)?;
323///
324/// let pixels = rasterizer.finish();
325/// assert_eq!(
326///     pixels.mapv(|v| v as u8),
327///     array![
328///         [0, 0, 3, 0],
329///         [0, 3, 3, 0],
330///         [3, 3, 0, 0],
331///         [3, 0, 0, 0],
332///         [0, 0, 0, 7]
333///     ]
334/// );
335/// # Ok(())}
336/// ```
337#[derive(Clone, Debug)]
338pub struct Rasterizer<Label> {
339    pixels: Array2<Label>,
340    geo_to_pix: Option<Transform>,
341    algorithm: MergeAlgorithm,
342    foreground: Label,
343    previous_burnt_points: HashSet<(usize, usize)>,
344    current_burnt_points: HashSet<(usize, usize)>,
345}
346
347impl<Label> Rasterizer<Label>
348where
349    Label: Copy + Add<Output = Label>,
350{
351    pub fn new(
352        width: usize,
353        height: usize,
354        geo_to_pix: Option<Transform>,
355        algorithm: MergeAlgorithm,
356        background: Label,
357    ) -> Self {
358        let pixels = Array2::from_elem((height, width), background);
359        Rasterizer {
360            pixels,
361            geo_to_pix,
362            algorithm,
363            foreground: background,
364            previous_burnt_points: HashSet::new(),
365            current_burnt_points: HashSet::new(),
366        }
367    }
368
369    fn width(&self) -> usize {
370        self.pixels.shape()[1]
371    }
372
373    fn height(&self) -> usize {
374        self.pixels.shape()[0]
375    }
376
377    /// Retrieve the transform.
378    pub fn geo_to_pix(&self) -> Option<Transform> {
379        self.geo_to_pix
380    }
381
382    // For MergeAlgorithm::Add, we have to ensure that we don't fill
383    // in the same pixels twice for the same; this only matters for
384    // the line drawing algorithm. To ensure that, we maintain a pair
385    // of index sets: one describing the indices of pixels we've
386    // filled in for the last line in the current line string and one
387    // describing the indices of pixels we've filled in for the
388    // current line of the current linestring. When we start a new
389    // line string, we clear both. And when we advance from one line
390    // to another within the same linestring, we swap them and clear
391    // the current one. While drawing each line, for the
392    // vertical/horizontal cases, we only refuse to fill a pixel if it
393    // was filled in the previous iteration but double filling from
394    // the current iteration is fine. But when drawing non-vertical
395    // non-horizontal lines, we refuse to fill pixels if we've filled
396    // them in either the current or previous iteration.
397
398    // aka clear()
399    fn new_linestring(&mut self) {
400        self.previous_burnt_points.clear();
401        self.current_burnt_points.clear();
402    }
403
404    fn new_line(&mut self) {
405        std::mem::swap(
406            &mut self.previous_burnt_points,
407            &mut self.current_burnt_points,
408        );
409        self.current_burnt_points.clear();
410    }
411
412    fn fill_pixel(&mut self, ix: usize, iy: usize) {
413        debug_assert!(ix < self.width());
414        debug_assert!(iy < self.height());
415        let mut slice = self.pixels.slice_mut(s![iy, ix]);
416        match self.algorithm {
417            MergeAlgorithm::Replace => slice.fill(self.foreground),
418            MergeAlgorithm::Add => {
419                slice.mapv_inplace(|v| v + self.foreground);
420            }
421        }
422    }
423
424    fn fill_pixel_no_repeat(&mut self, ix: usize, iy: usize, use_current_too: bool) {
425        match self.algorithm {
426            MergeAlgorithm::Replace => {
427                self.fill_pixel(ix, iy);
428            }
429            MergeAlgorithm::Add => {
430                let point = (ix, iy);
431                let mut do_fill_pixel = !self.previous_burnt_points.contains(&point);
432                if use_current_too {
433                    do_fill_pixel = do_fill_pixel && !self.current_burnt_points.contains(&point);
434                }
435                if do_fill_pixel {
436                    self.fill_pixel(ix, iy);
437                    self.current_burnt_points.insert(point);
438                }
439            }
440        }
441    }
442
443    // The rasterization algorithm's performance is extremely
444    // sensitive to write ordering: it is focused on horizontal lines,
445    // so it performs much better when pixels that are horizontally
446    // adjacent are adjacent in memory (i.e., where the array we're
447    // writing to has x as the last dimension).
448
449    // Unlike the other fill_ methods, x_start and x_end are an
450    // exclusive range (..), not inclusive (..=).
451    fn fill_horizontal_line(&mut self, x_start: usize, x_end: usize, y: usize) {
452        let mut slice = self.pixels.slice_mut(s![y, x_start..x_end]);
453        match self.algorithm {
454            MergeAlgorithm::Replace => slice.fill(self.foreground),
455            MergeAlgorithm::Add => {
456                slice.mapv_inplace(|v| v + self.foreground);
457            }
458        }
459    }
460
461    fn fill_horizontal_line_no_repeat(&mut self, x_start: usize, x_end: usize, y: usize) {
462        for x in x_start..=x_end {
463            self.fill_pixel_no_repeat(x, y, true);
464        }
465    }
466
467    fn fill_vertical_line_no_repeat(&mut self, x: usize, y_start: usize, y_end: usize) {
468        for y in y_start..=y_end {
469            self.fill_pixel_no_repeat(x, y, false);
470        }
471    }
472
473    /// Rasterize one shape, which can be any type that [geo] provides
474    /// using any coordinate numeric type that can be converted into
475    /// `f64`.
476    pub fn rasterize<Coord, InputShape, ShapeAsF64>(
477        &mut self,
478        shape: &InputShape,
479        foreground: Label,
480    ) -> Result<()>
481    where
482        InputShape: MapCoords<Coord, f64, Output = ShapeAsF64>,
483        ShapeAsF64: Rasterize<Label> + for<'a> CoordsIter<'a, Scalar = f64> + MapCoordsInplace<f64>,
484        Coord: Into<f64> + Copy + Debug + Num + NumCast + PartialOrd,
485    {
486        // first, convert our input shape so that its coordinates are of type f64
487        let mut float = shape.map_coords(to_float);
488
489        // then ensure that all coordinates are finite or bail
490        let all_finite = float
491            .coords_iter()
492            .all(|coordinate| coordinate.x.is_finite() && coordinate.y.is_finite());
493        if !all_finite {
494            return Err(RasterizeError::NonFiniteCoordinate);
495        }
496
497        self.foreground = foreground;
498
499        // use `geo_to_pix` to convert geographic coordinates to image
500        // coordinates, if it is available
501        match self.geo_to_pix {
502            None => float,
503            Some(transform) => {
504                float.map_coords_inplace(|&(x, y)| {
505                    transform.transform_point(EuclidPoint::new(x, y)).to_tuple()
506                });
507                float
508            }
509        }
510        .rasterize(self); // and then rasterize!
511
512        Ok(())
513    }
514
515    /// Retrieve the completed raster array.
516    pub fn finish(self) -> Array2<Label> {
517        self.pixels
518    }
519}
520
521impl<Label> Rasterize<Label> for Point<f64>
522where
523    Label: Copy + Add<Output = Label>,
524{
525    fn rasterize(&self, rasterizer: &mut Rasterizer<Label>) {
526        if self.x() >= 0. && self.y() >= 0. {
527            let x = self.x().floor() as usize;
528            let y = self.y().floor() as usize;
529            if x < rasterizer.width() && y < rasterizer.height() {
530                rasterizer.fill_pixel(x, y);
531            }
532        }
533    }
534}
535
536impl<Label> Rasterize<Label> for MultiPoint<f64>
537where
538    Label: Copy + Add<Output = Label>,
539{
540    fn rasterize(&self, rasterizer: &mut Rasterizer<Label>) {
541        self.iter().for_each(|point| point.rasterize(rasterizer));
542    }
543}
544
545impl<Label> Rasterize<Label> for Rect<f64>
546where
547    Label: Copy + Add<Output = Label>,
548{
549    fn rasterize(&self, rasterizer: &mut Rasterizer<Label>) {
550        // Although it is tempting to make a really fast direct
551        // implementation, we're going to convert to a polyon and rely
552        // on that impl, in part because affine transforms can easily
553        // rotate or shear the rectangle so that it is no longer axis
554        // aligned.
555        self.to_polygon().rasterize(rasterizer);
556    }
557}
558
559impl<Label> Rasterize<Label> for Line<f64>
560where
561    Label: Copy + Add<Output = Label>,
562{
563    fn rasterize(&self, rasterizer: &mut Rasterizer<Label>) {
564        rasterizer.new_linestring();
565        rasterize_line(self, rasterizer);
566    }
567}
568
569impl<Label> Rasterize<Label> for LineString<f64>
570where
571    Label: Copy + Add<Output = Label>,
572{
573    fn rasterize(&self, rasterizer: &mut Rasterizer<Label>) {
574        // It is tempting to make this impl treat closed `LineString`s
575        // as polygons without holes: to just fill them. GDAL seems
576        // like it should do that (`gv_rasterize_one_shape` in
577        // `gdalrasterize.cpp` has a default clause that just invokes
578        // the polygon filling code), but in practice GDAL treats
579        // closed `LinearRings` as `LineSegments` and doesn't fill
580        // them and I'm not sure why. Perhaps `LinearRings` are more
581        // of an internal implementation detail?
582        rasterizer.new_linestring();
583        self.lines().for_each(|line| {
584            rasterizer.new_line();
585            rasterize_line(&line, rasterizer);
586        });
587    }
588}
589
590impl<Label> Rasterize<Label> for MultiLineString<f64>
591where
592    Label: Copy + Add<Output = Label>,
593{
594    fn rasterize(&self, rasterizer: &mut Rasterizer<Label>) {
595        self.iter()
596            .for_each(|line_string| line_string.rasterize(rasterizer));
597    }
598}
599
600impl<Label> Rasterize<Label> for Polygon<f64>
601where
602    Label: Copy + Add<Output = Label>,
603{
604    fn rasterize(&self, rasterizer: &mut Rasterizer<Label>) {
605        rasterize_polygon(self.exterior(), self.interiors(), rasterizer);
606    }
607}
608
609impl<Label> Rasterize<Label> for MultiPolygon<f64>
610where
611    Label: Copy + Add<Output = Label>,
612{
613    fn rasterize(&self, rasterizer: &mut Rasterizer<Label>) {
614        self.iter().for_each(|poly| poly.rasterize(rasterizer));
615    }
616}
617
618impl<Label> Rasterize<Label> for Triangle<f64>
619where
620    Label: Copy + Add<Output = Label>,
621{
622    fn rasterize(&self, rasterizer: &mut Rasterizer<Label>) {
623        self.to_polygon().rasterize(rasterizer)
624    }
625}
626
627impl<Label> Rasterize<Label> for Geometry<f64>
628where
629    Label: Copy + Add<Output = Label>,
630{
631    fn rasterize(&self, rasterizer: &mut Rasterizer<Label>) {
632        match self {
633            Geometry::Point(point) => point.rasterize(rasterizer),
634            Geometry::Line(line) => line.rasterize(rasterizer),
635            Geometry::LineString(ls) => ls.rasterize(rasterizer),
636            Geometry::Polygon(poly) => poly.rasterize(rasterizer),
637            Geometry::GeometryCollection(gc) => gc.rasterize(rasterizer),
638            Geometry::MultiPoint(points) => points.rasterize(rasterizer),
639            Geometry::MultiLineString(lines) => lines.rasterize(rasterizer),
640            Geometry::MultiPolygon(polys) => polys.rasterize(rasterizer),
641            Geometry::Rect(rect) => rect.rasterize(rasterizer),
642            Geometry::Triangle(tri) => tri.rasterize(rasterizer),
643        }
644    }
645}
646
647impl<Label> Rasterize<Label> for GeometryCollection<f64>
648where
649    Label: Copy + Add<Output = Label>,
650{
651    fn rasterize(&self, rasterizer: &mut Rasterizer<Label>) {
652        self.iter().for_each(|thing| thing.rasterize(rasterizer));
653    }
654}
655
656#[cfg(test)]
657mod tests;
geo_rasterize/lib.rs

geo_rasterize/
lib.rs
