Fix overflow when using large raster sizes

benches/bench.rs

@@ -90,8 +90,8 @@ fn bench_contourbuilder_isobands_volcano_without_xy_step_xy_origin(b: &mut Bench
         .iter()
         .map(|x| x.as_f64().unwrap())
         .collect();
-    let h = raw_data["height"].as_u64().unwrap() as u32;
-    let w = raw_data["width"].as_u64().unwrap() as u32;
+    let h = raw_data["height"].as_u64().unwrap() as usize;
+    let w = raw_data["width"].as_u64().unwrap() as usize;
 
     b.iter(|| {
         black_box(
@@ -118,8 +118,8 @@ fn bench_contourbuilder_isobands_pot_pop_fr_without_xy_step_xy_origin(b: &mut Be
         .iter()
         .map(|x| x.as_f64().unwrap())
         .collect();
-    let h = raw_data["height"].as_u64().unwrap() as u32;
-    let w = raw_data["width"].as_u64().unwrap() as u32;
+    let h = raw_data["height"].as_u64().unwrap() as usize;
+    let w = raw_data["width"].as_u64().unwrap() as usize;
 
     b.iter(|| {
         black_box(

examples/ex.rs

@@ -21,8 +21,8 @@ fn main() {
             }
         })
         .collect();
-    let h = raw_data["height"].as_u64().unwrap() as u32;
-    let w = raw_data["width"].as_u64().unwrap() as u32;
+    let h = raw_data["height"].as_u64().unwrap() as usize;
+    let w = raw_data["width"].as_u64().unwrap() as usize;
 
     let x_origin = -6.144721171428571;
     let y_origin = 51.78171334283718;
@@ -75,8 +75,8 @@ fn main() {
             }
         })
         .collect();
-    let h = raw_data["height"].as_u64().unwrap() as u32;
-    let w = raw_data["width"].as_u64().unwrap() as u32;
+    let h = raw_data["height"].as_u64().unwrap() as usize;
+    let w = raw_data["width"].as_u64().unwrap() as usize;
 
     let contours = ContourBuilder::new(w, h, true)
         .isobands(

src/area.rs

@@ -1,10 +1,15 @@
 use crate::{Float, Pt};
 
-pub fn area(ring: &[Pt]) -> Float {
+#[allow(clippy::unnecessary_cast)]
+// Note that we need to disable the clippy warning about unnecessary casts
+// because of the "f32" optional feature (and because we want to ensure we always
+// use "f64" in this function, both in the default feature and in the "f32" feature).
+pub fn area(ring: &[Pt]) -> f64 {
     let n = ring.len();
-    let mut area = ring[n - 1].y * ring[0].x - ring[n - 1].x * ring[0].y;
+    let mut area =
+        ring[n - 1].y as f64 * ring[0].x as f64 - ring[n - 1].x as f64 * ring[0].y as f64;
     for i in 1..n {
-        area += ring[i - 1].y * ring[i].x - ring[i - 1].x * ring[i].y;
+        area += ring[i - 1].y as f64 * ring[i].x as f64 - ring[i - 1].x as f64 * ring[i].y as f64;
     }
     // Note that in the shoelace formula you need to divide this result by 2 to get the actual area.
     // Here we skip this division because we only use this area formula to calculate the winding

src/contourbuilder.rs

@@ -12,9 +12,9 @@ use rustc_hash::FxHashMap;
 /// [`contour_rings`]: fn.contour_rings.html
 pub struct ContourBuilder {
     /// The number of columns in the grid
-    dx: u32,
+    dx: usize,
     /// The number of rows in the grid
-    dy: u32,
+    dy: usize,
     /// Whether to smooth the contours
     smooth: bool,
     /// The horizontal coordinate for the origin of the grid.
@@ -38,7 +38,7 @@ impl ContourBuilder {
     /// * `dx` - The number of columns in the grid.
     /// * `dy` - The number of rows in the grid.
     /// * `smooth` - Whether or not the generated rings will be smoothed using linear interpolation.
-    pub fn new(dx: u32, dy: u32, smooth: bool) -> Self {
+    pub fn new(dx: usize, dy: usize, smooth: bool) -> Self {
         ContourBuilder {
             dx,
             dy,
@@ -83,18 +83,18 @@ impl ContourBuilder {
             .map(|point| {
                 let x = point.x;
                 let y = point.y;
-                let xt = x.trunc() as u32;
-                let yt = y.trunc() as u32;
+                let xt = x.trunc() as usize;
+                let yt = y.trunc() as usize;
                 let mut v0;
-                let ix = (yt * dx + xt) as usize;
+                let ix = yt * dx + xt;
                 if ix < len_values {
                     let v1 = values[ix];
                     if x > 0.0 && x < (dx as Float) && (xt as Float - x).abs() < Float::EPSILON {
-                        v0 = values[(yt * dx + xt - 1) as usize];
+                        v0 = values[yt * dx + xt - 1];
                         point.x = x + (value - v0) / (v1 - v0) - 0.5;
                     }
                     if y > 0.0 && y < (dy as Float) && (yt as Float - y).abs() < Float::EPSILON {
-                        v0 = values[((yt - 1) * dx + xt) as usize];
+                        v0 = values[(yt - 1) * dx + xt];
                         point.y = y + (value - v0) / (v1 - v0) - 0.5;
                     }
                 }
@@ -112,7 +112,7 @@ impl ContourBuilder {
     /// * `values` - The slice of values to be used.
     /// * `thresholds` - The slice of thresholds values to be used.
     pub fn lines(&self, values: &[Float], thresholds: &[Float]) -> Result<Vec<Line>> {
-        if values.len() as u32 != self.dx * self.dy {
+        if values.len() != self.dx * self.dy {
             return Err(new_error(ErrorKind::BadDimension));
         }
         let mut isoring = IsoRingBuilder::new(self.dx, self.dy);
@@ -163,7 +163,7 @@ impl ContourBuilder {
     /// * `values` - The slice of values to be used.
     /// * `thresholds` - The slice of thresholds values to be used.
     pub fn contours(&self, values: &[Float], thresholds: &[Float]) -> Result<Vec<Contour>> {
-        if values.len() as u32 != self.dx * self.dy {
+        if values.len() != self.dx * self.dy {
             return Err(new_error(ErrorKind::BadDimension));
         }
         let mut isoring = IsoRingBuilder::new(self.dx, self.dy);
@@ -232,7 +232,7 @@ impl ContourBuilder {
         // We will compute rings as previously, but we will
         // iterate over the contours in pairs and use the paths from the lower threshold
         // and the path from the upper threshold to create the isoband.
-        if values.len() as u32 != self.dx * self.dy {
+        if values.len() != self.dx * self.dy {
             return Err(new_error(ErrorKind::BadDimension));
         }
         if thresholds.len() < 2 {

src/isoringbuilder.rs

@@ -49,7 +49,12 @@ struct Fragment {
 /// * `threshold` - The threshold value.
 /// * `dx` - The number of columns in the grid.
 /// * `dy` - The number of rows in the grid.
-pub fn contour_rings(values: &[Float], threshold: Float, dx: u32, dy: u32) -> Result<Vec<Ring>> {
+pub fn contour_rings(
+    values: &[Float],
+    threshold: Float,
+    dx: usize,
+    dy: usize,
+) -> Result<Vec<Ring>> {
     let mut isoring = IsoRingBuilder::new(dx, dy);
     isoring.compute(values, threshold)
 }
@@ -59,8 +64,8 @@ pub struct IsoRingBuilder {
     fragment_by_start: FxHashMap<usize, usize>,
     fragment_by_end: FxHashMap<usize, usize>,
     f: Slab<Fragment>,
-    dx: u32,
-    dy: u32,
+    dx: usize,
+    dy: usize,
     is_empty: bool,
 }
 
@@ -70,7 +75,7 @@ impl IsoRingBuilder {
     ///
     /// * `dx` - The number of columns in the grid.
     /// * `dy` - The number of rows in the grid.
-    pub fn new(dx: u32, dy: u32) -> Self {
+    pub fn new(dx: usize, dy: usize) -> Self {
         IsoRingBuilder {
             fragment_by_start: FxHashMap::default(),
             fragment_by_end: FxHashMap::default(),
@@ -103,8 +108,8 @@ impl IsoRingBuilder {
             self.clear();
         }
         let mut result = Vec::new();
-        let dx = self.dx as i32;
-        let dy = self.dy as i32;
+        let dx = self.dx as i64;
+        let dy = self.dy as i64;
         let mut x = -1;
         let mut y = -1;
         let mut t0;
@@ -113,68 +118,63 @@ impl IsoRingBuilder {
         let mut t3;
 
         // Special case for the first row (y = -1, t2 = t3 = 0).
-        t1 = (values[0] >= threshold) as u32;
-        case_stitch!((t1 << 1) as usize, x, y, &mut result);
+        t1 = (values[0] >= threshold) as usize;
+        case_stitch!(t1 << 1, x, y, &mut result);
         x += 1;
         while x < dx - 1 {
             t0 = t1;
-            t1 = (values[(x + 1) as usize] >= threshold) as u32;
-            case_stitch!((t0 | t1 << 1) as usize, x, y, &mut result);
+            t1 = (values[(x + 1) as usize] >= threshold) as usize;
+            case_stitch!(t0 | t1 << 1, x, y, &mut result);
             x += 1;
         }
-        case_stitch!(t1 as usize, x, y, &mut result);
+        case_stitch!(t1, x, y, &mut result);
 
         // General case for the intermediate rows.
         y += 1;
         while y < dy - 1 {
             x = -1;
-            t1 = (values[(y * dx + dx) as usize] >= threshold) as u32;
-            t2 = (values[(y * dx) as usize] >= threshold) as u32;
-            case_stitch!((t1 << 1 | t2 << 2) as usize, x, y, &mut result);
+            t1 = (values[(y * dx + dx) as usize] >= threshold) as usize;
+            t2 = (values[(y * dx) as usize] >= threshold) as usize;
+            case_stitch!(t1 << 1 | t2 << 2, x, y, &mut result);
             x += 1;
             while x < dx - 1 {
                 t0 = t1;
-                t1 = (values[(y * dx + dx + x + 1) as usize] >= threshold) as u32;
+                t1 = (values[(y * dx + dx + x + 1) as usize] >= threshold) as usize;
                 t3 = t2;
-                t2 = (values[(y * dx + x + 1) as usize] >= threshold) as u32;
-                case_stitch!(
-                    (t0 | t1 << 1 | t2 << 2 | t3 << 3) as usize,
-                    x,
-                    y,
-                    &mut result
-                );
+                t2 = (values[(y * dx + x + 1) as usize] >= threshold) as usize;
+                case_stitch!(t0 | t1 << 1 | t2 << 2 | t3 << 3, x, y, &mut result);
                 x += 1;
             }
-            case_stitch!((t1 | t2 << 3) as usize, x, y, &mut result);
+            case_stitch!(t1 | t2 << 3, x, y, &mut result);
             y += 1;
         }
 
         // Special case for the last row (y = dy - 1, t0 = t1 = 0).
         x = -1;
-        t2 = (values[(y * dx) as usize] >= threshold) as u32;
-        case_stitch!((t2 << 2) as usize, x, y, &mut result);
+        t2 = (values[(y * dx) as usize] >= threshold) as usize;
+        case_stitch!(t2 << 2, x, y, &mut result);
         x += 1;
         while x < dx - 1 {
             t3 = t2;
-            t2 = (values[(y * dx + x + 1) as usize] >= threshold) as u32;
-            case_stitch!((t2 << 2 | t3 << 3) as usize, x, y, &mut result);
+            t2 = (values[(y * dx + x + 1) as usize] >= threshold) as usize;
+            case_stitch!(t2 << 2 | t3 << 3, x, y, &mut result);
             x += 1;
         }
-        case_stitch!((t2 << 3) as usize, x, y, &mut result);
+        case_stitch!(t2 << 3, x, y, &mut result);
         self.is_empty = false;
         Ok(result)
     }
 
     fn index(&self, point: &Pt) -> usize {
-        (point.x * 2.0 + point.y * (self.dx as Float + 1.) * 4.) as usize
+        (point.x as usize) * 2 + (point.y as usize) * (self.dx + 1usize) * 4
     }
 
     // Stitchs segments to rings.
     fn stitch(
         &mut self,
         line: &[Vec<Float>],
-        x: i32,
-        y: i32,
+        x: i64,
+        y: i64,
         result: &mut Vec<Ring>,
     ) -> Result<()> {
         let start = Pt {