jagua_rs/geometry/fail_fast/
piers.rs

1use itertools::{Itertools, izip};
2use ndarray::Array;
3use ordered_float::NotNan;
4use rand_distr::num_traits::FloatConst;
5
6use crate::geometry::Transformation;
7use crate::geometry::geo_traits::{CollidesWith, DistanceTo, Transformable};
8use crate::geometry::primitives::Circle;
9use crate::geometry::primitives::Edge;
10use crate::geometry::primitives::Point;
11use crate::geometry::primitives::Rect;
12use crate::geometry::primitives::SPolygon;
13
14use anyhow::{Result, bail};
15
16static RAYS_PER_ANGLE: usize = if cfg!(debug_assertions) { 10 } else { 200 };
17static N_ANGLES: usize = if cfg!(debug_assertions) { 4 } else { 90 };
18static N_POINTS_PER_DIMENSION: usize = if cfg!(debug_assertions) { 10 } else { 100 };
19static CLIPPING_TRIM: f32 = 0.999;
20static ACTION_RADIUS_RATIO: f32 = 0.10;
21
22/// Generates a set of `n` *piers* - line segments fully contained within `shape`.
23/// This function generates them in such a way as to *cover* areas of the `shape` that are
24/// poorly represented by `poles` as well as possible.
25pub fn generate_piers(shape: &SPolygon, n: usize, poles: &[Circle]) -> Result<Vec<Edge>> {
26    if n == 0 {
27        return Ok(vec![]);
28    }
29
30    //Start by creating a set of N_TESTS_PER_ANGLE vertical lines across the bounding box
31    let bbox = shape.bbox;
32    let expanded_bbox = bbox.clone().inflate_to_square();
33    let centroid = shape.centroid();
34    //vertical ray from the centroid
35    let base_ray = Edge::try_new(
36        Point(centroid.0, centroid.1 - 2.0 * expanded_bbox.height()),
37        Point(centroid.0, centroid.1 + 2.0 * expanded_bbox.height()),
38    )
39    .unwrap();
40
41    let transformations = generate_ray_transformations(expanded_bbox, RAYS_PER_ANGLE, N_ANGLES);
42
43    //transform the base edge by each transformation
44    let rays = transformations
45        .into_iter()
46        .map(|t| base_ray.transform_clone(&t))
47        .collect_vec();
48
49    //clip the lines to the shape
50    let clipped_rays = rays.iter().flat_map(|l| clip(shape, l)).collect_vec();
51    let grid_of_unrepresented_points =
52        generate_unrepresented_point_grid(expanded_bbox, shape, poles, N_POINTS_PER_DIMENSION);
53
54    let mut selected_piers = Vec::new();
55
56    let radius_of_ray_influence = ACTION_RADIUS_RATIO * expanded_bbox.width();
57    let forfeit_distance = f32::sqrt(bbox.width().powi(2) * bbox.height().powi(2));
58
59    for _ in 0..n {
60        let min_distance_selected_rays = min_distances_to_rays(
61            &grid_of_unrepresented_points,
62            &selected_piers,
63            forfeit_distance,
64        );
65        let min_distance_poles =
66            min_distances_to_poles(&grid_of_unrepresented_points, poles, forfeit_distance);
67
68        let loss_values = clipped_rays
69            .iter()
70            .map(|new_ray| {
71                loss_function(
72                    new_ray,
73                    &grid_of_unrepresented_points,
74                    &min_distance_selected_rays,
75                    &min_distance_poles,
76                    radius_of_ray_influence,
77                )
78            })
79            .map(|x| NotNan::new(x).unwrap())
80            .collect_vec();
81
82        let min_loss_ray = clipped_rays
83            .iter()
84            .enumerate()
85            .min_by_key(|(i, _)| loss_values[*i])
86            .map(|(_i, ray)| ray);
87
88        match min_loss_ray {
89            None => bail!("no ray found"),
90            Some(ray) => selected_piers.push(*ray),
91        }
92    }
93    Ok(selected_piers)
94}
95
96#[allow(clippy::cast_precision_loss)]
97fn generate_ray_transformations(
98    bbox: Rect,
99    rays_per_angle: usize,
100    n_angles: usize,
101) -> Vec<Transformation> {
102    //translations
103    let dx = bbox.width() / rays_per_angle as f32;
104    let translations = (0..rays_per_angle)
105        .map(|i| bbox.x_min + dx * i as f32)
106        .map(|x| Transformation::from_translation((x, 0.0)))
107        .collect_vec();
108
109    let angles = Array::linspace(0.0, f32::PI(), n_angles + 1).to_vec();
110    let angles_slice = &angles[0..n_angles]; //skip the last angle, which is the same as the first
111
112    //rotate the translations by each angle
113    angles_slice
114        .iter()
115        .flat_map(|angle| {
116            translations
117                .iter()
118                .cloned()
119                .map(move |translation| translation.rotate(*angle))
120        })
121        .collect_vec()
122}
123
124//clips a ray against the border of a polygon, potentially creating multiple "clips"
125fn clip(shape: &SPolygon, ray: &Edge) -> Vec<Edge> {
126    //both ends of the ray should be outside the shape
127    assert!(!shape.collides_with(&ray.start) && !shape.collides_with(&ray.end));
128
129    //collect all intersections of the ray with the shape, sorted by distance to the ray's start
130    let intersections = shape
131        .edge_iter()
132        .filter_map(|edge| edge.collides_at(ray))
133        .sorted_by_key(|p| NotNan::new(ray.start.distance_to(p)).unwrap())
134        .collect_vec();
135
136    //every pair of (sorted) intersections defines a clipped line
137    intersections
138        .chunks(2)
139        .filter_map(|pair| {
140            if pair.len() == 1 {
141                return None;
142            }
143            let start = pair[0];
144            let end = pair[1];
145            if start == end {
146                None
147            } else {
148                Some(Edge::try_new(start, end).unwrap().scale(CLIPPING_TRIM))
149            }
150        })
151        .collect_vec()
152}
153
154fn generate_unrepresented_point_grid(
155    bbox: Rect,
156    shape: &SPolygon,
157    poles: &[Circle],
158    n_points_per_dimension: usize,
159) -> Vec<Point> {
160    let x_range = Array::linspace(bbox.x_min, bbox.x_max, n_points_per_dimension);
161    let y_range = Array::linspace(bbox.y_min, bbox.y_max, n_points_per_dimension);
162
163    x_range
164        .iter()
165        .flat_map(|x| {
166            y_range
167                .iter()
168                .map(move |y| Point::from((*x, *y))) //create the points
169                .filter(|p| shape.collides_with(p)) //make sure they are in the shape
170                .filter(|p| poles.iter().all(|c| !c.collides_with(p)))
171        })
172        .collect_vec()
173}
174
175fn loss_function(
176    new_ray: &Edge,
177    point_grid: &[Point],
178    min_distance_to_rays: &[f32],
179    min_distance_to_poles: &[f32],
180    radius_of_ray_influence: f32,
181) -> f32 {
182    //every point in the grid gets a certain score, sum of all these scores is the loss function
183    //the score depends on how close it is to being "represented" by either a pole or a ray
184    //rays have a certain radius of influence, outside which they don't count. Poles have no such radius
185    //the score is the squared distance to the closest ray or pole
186
187    izip!(
188        point_grid.iter(),
189        min_distance_to_rays.iter(),
190        min_distance_to_poles.iter()
191    )
192    .map(|(p, min_distance_to_existing_ray, min_distance_to_pole)| {
193        let distance_to_new_ray = new_ray.distance_to(p);
194
195        let min_distance_to_ray = f32::min(*min_distance_to_existing_ray, distance_to_new_ray);
196
197        if min_distance_to_ray < radius_of_ray_influence {
198            f32::min(*min_distance_to_pole, min_distance_to_ray)
199        } else {
200            *min_distance_to_pole
201        }
202    })
203    .map(|d| d.powi(2))
204    .sum()
205}
206
207fn min_distances_to_rays(points: &[Point], rays: &[Edge], forfeit_distance: f32) -> Vec<f32> {
208    points
209        .iter()
210        .map(|p| {
211            rays.iter()
212                .map(|r| r.distance_to(p))
213                .fold(forfeit_distance, f32::min)
214        })
215        .collect_vec()
216}
217
218fn min_distances_to_poles(points: &[Point], poles: &[Circle], forfeit_distance: f32) -> Vec<f32> {
219    points
220        .iter()
221        .map(|p| {
222            poles
223                .iter()
224                .map(|c| c.distance_to(p))
225                .fold(forfeit_distance, f32::min)
226        })
227        .collect_vec()
228}
jagua_rs/geometry/fail_fast/piers.rs

jagua_rs/geometry/fail_fast/
piers.rs
