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
96fn generate_ray_transformations(
97    bbox: Rect,
98    rays_per_angle: usize,
99    n_angles: usize,
100) -> Vec<Transformation> {
101    //translations
102    let dx = bbox.width() / rays_per_angle as f32;
103    let translations = (0..rays_per_angle)
104        .map(|i| bbox.x_min + dx * i as f32)
105        .map(|x| Transformation::from_translation((x, 0.0)))
106        .collect_vec();
107
108    let angles = Array::linspace(0.0, f32::PI(), n_angles + 1).to_vec();
109    let angles_slice = &angles[0..n_angles]; //skip the last angle, which is the same as the first
110
111    //rotate the translations by each angle
112    angles_slice
113        .iter()
114        .flat_map(|angle| {
115            translations
116                .iter()
117                .cloned()
118                .map(move |translation| translation.rotate(*angle))
119        })
120        .collect_vec()
121}
122
123//clips a ray against the border of a polygon, potentially creating multiple "clips"
124fn clip(shape: &SPolygon, ray: &Edge) -> Vec<Edge> {
125    //both ends of the ray should be outside the shape
126    assert!(!shape.collides_with(&ray.start) && !shape.collides_with(&ray.end));
127
128    //collect all intersections of the ray with the shape, sorted by distance to the ray's start
129    let intersections = shape
130        .edge_iter()
131        .flat_map(|edge| edge.collides_at(ray))
132        .sorted_by_key(|p| NotNan::new(ray.start.distance_to(p)).unwrap())
133        .collect_vec();
134
135    //every pair of (sorted) intersections defines a clipped line
136    intersections
137        .chunks(2)
138        .flat_map(|pair| {
139            if pair.len() == 1 {
140                return None;
141            }
142            let start = pair[0];
143            let end = pair[1];
144            if start != end {
145                Some(Edge::try_new(start, end).unwrap().scale(CLIPPING_TRIM))
146            } else {
147                None
148            }
149        })
150        .collect_vec()
151}
152
153fn generate_unrepresented_point_grid(
154    bbox: Rect,
155    shape: &SPolygon,
156    poles: &[Circle],
157    n_points_per_dimension: usize,
158) -> Vec<Point> {
159    let x_range = Array::linspace(bbox.x_min, bbox.x_max, n_points_per_dimension);
160    let y_range = Array::linspace(bbox.y_min, bbox.y_max, n_points_per_dimension);
161
162    x_range
163        .iter()
164        .flat_map(|x| {
165            y_range
166                .iter()
167                .map(move |y| Point::from((*x, *y))) //create the points
168                .filter(|p| shape.collides_with(p)) //make sure they are in the shape
169                .filter(|p| poles.iter().all(|c| !c.collides_with(p)))
170        })
171        .collect_vec()
172}
173
174fn loss_function(
175    new_ray: &Edge,
176    point_grid: &[Point],
177    min_distance_to_rays: &[f32],
178    min_distance_to_poles: &[f32],
179    radius_of_ray_influence: f32,
180) -> f32 {
181    //every point in the grid gets a certain score, sum of all these scores is the loss function
182    //the score depends on how close it is to being "represented" by either a pole or a ray
183    //rays have a certain radius of influence, outside which they don't count. Poles have no such radius
184    //the score is the squared distance to the closest ray or pole
185
186    izip!(
187        point_grid.iter(),
188        min_distance_to_rays.iter(),
189        min_distance_to_poles.iter()
190    )
191    .map(|(p, min_distance_to_existing_ray, min_distance_to_pole)| {
192        let distance_to_new_ray = new_ray.distance_to(p);
193
194        let min_distance_to_ray = f32::min(*min_distance_to_existing_ray, distance_to_new_ray);
195
196        match min_distance_to_ray < radius_of_ray_influence {
197            true => f32::min(*min_distance_to_pole, min_distance_to_ray),
198            false => *min_distance_to_pole,
199        }
200    })
201    .map(|d| d.powi(2))
202    .sum()
203}
204
205fn min_distances_to_rays(points: &[Point], rays: &[Edge], forfeit_distance: f32) -> Vec<f32> {
206    points
207        .iter()
208        .map(|p| {
209            rays.iter()
210                .map(|r| r.distance_to(p))
211                .fold(forfeit_distance, f32::min)
212        })
213        .collect_vec()
214}
215
216fn min_distances_to_poles(points: &[Point], poles: &[Circle], forfeit_distance: f32) -> Vec<f32> {
217    points
218        .iter()
219        .map(|p| {
220            poles
221                .iter()
222                .map(|c| c.distance_to(p))
223                .fold(forfeit_distance, f32::min)
224        })
225        .collect_vec()
226}
jagua_rs/geometry/fail_fast/piers.rs

jagua_rs/geometry/fail_fast/
piers.rs
