pub fn find_contours_with_threshold<T>(image: &GrayImage, threshold: u8) -> Vec<Contour<T>>
where
    T: Num + NumCast + Copy + PartialEq + Eq,
{
    let width = image.width() as usize;
    let height = image.height() as usize;
    let mut image_values = vec![0i32; height * width];

    let at = |x, y| x + width * y;
    let get_position_if_non_zero_pixel = |image: &[i32], curr: Point<i32>| {
        let (x, y) = (curr.x, curr.y);
        let in_bounds = x > -1 && x < width as i32 && y > -1 && y < height as i32;

        if in_bounds && image[at(x as usize, y as usize)] != 0 {
            Some(Point::new(x as usize, y as usize))
        } else {
            None
        }
    };

    for y in 0..height {
        for x in 0..width {
            if image.get_pixel(x as u32, y as u32).0[0] > threshold {
                image_values[at(x, y)] = 1;
            }
        }
    }

    let mut diffs = VecDeque::from(vec![
        Point::new(-1, 0),  // w
        Point::new(-1, -1), // nw
        Point::new(0, -1),  // n
        Point::new(1, -1),  // ne
        Point::new(1, 0),   // e
        Point::new(1, 1),   // se
        Point::new(0, 1),   // s
        Point::new(-1, 1),  // sw
    ]);

    let mut contours: Vec<Contour<T>> = Vec::new();
    let mut curr_border_num = 1;

    for y in 0..height {
        let mut parent_border_num = 1;

        for x in 0..width {
            if image_values[at(x, y)] == 0 {
                continue;
            }

            if let Some((adj, border_type)) = if image_values[at(x, y)] == 1
                && x > 0
                && image_values[at(x - 1, y)] == 0
            {
                Some((Point::new(x - 1, y), BorderType::Outer))
            } else if image_values[at(x, y)] > 0 && x + 1 < width && image_values[at(x + 1, y)] == 0
            {
                if image_values[at(x, y)] > 1 {
                    parent_border_num = image_values[at(x, y)] as usize;
                }
                Some((Point::new(x + 1, y), BorderType::Hole))
            } else {
                None
            } {
                curr_border_num += 1;

                let parent = if parent_border_num > 1 {
                    let parent_index = parent_border_num - 2;
                    let parent_contour = &contours[parent_index];
                    if (border_type == BorderType::Outer)
                        ^ (parent_contour.border_type == BorderType::Outer)
                    {
                        Some(parent_index)
                    } else {
                        parent_contour.parent
                    }
                } else {
                    None
                };

                let mut contour_points = Vec::new();
                let curr = Point::new(x, y);
                rotate_to_value(&mut diffs, adj.to_i32() - curr.to_i32());

                if let Some(pos1) = diffs.iter().find_map(|diff| {
                    get_position_if_non_zero_pixel(&image_values, curr.to_i32() + *diff)
                }) {
                    let mut pos2 = pos1;
                    let mut pos3 = curr;

                    loop {
                        contour_points
                            .push(Point::new(cast(pos3.x).unwrap(), cast(pos3.y).unwrap()));
                        rotate_to_value(&mut diffs, pos2.to_i32() - pos3.to_i32());
                        let pos4 = diffs
                            .iter()
                            .rev() // counter-clockwise
                            .find_map(|diff| {
                                get_position_if_non_zero_pixel(&image_values, pos3.to_i32() + *diff)
                            })
                            .unwrap();

                        let mut is_right_edge = false;
                        for diff in diffs.iter().rev() {
                            if *diff == (pos4.to_i32() - pos3.to_i32()) {
                                break;
                            }
                            if *diff == Point::new(1, 0) {
                                is_right_edge = true;
                                break;
                            }
                        }

                        if pos3.x + 1 == width || is_right_edge {
                            image_values[at(pos3.x, pos3.y)] = -curr_border_num;
                        } else if image_values[at(pos3.x, pos3.y)] == 1 {
                            image_values[at(pos3.x, pos3.y)] = curr_border_num;
                        }

                        if pos4 == curr && pos3 == pos1 {
                            break;
                        }

                        pos2 = pos3;
                        pos3 = pos4;
                    }
                } else {
                    contour_points.push(Point::new(cast(x).unwrap(), cast(y).unwrap()));
                    image_values[at(x, y)] = -curr_border_num;
                }
                contours.push(Contour::new(contour_points, border_type, parent));
            }

            if image_values[at(x, y)] != 1 {
                parent_border_num = image_values[at(x, y)].unsigned_abs() as usize;
            }
        }
    }

    contours
}
Finds all borders of foreground regions in an image. All pixels with intensity strictly greater
than `threshold` are treated as belonging to the foreground.

Based on the algorithm proposed by Suzuki and Abe: Topological Structural
Analysis of Digitized Binary Images by Border Following.
src/contours.rs