1  //! An algorithm for laying out curves in a given order.
  2  //!
  3  //! Our boolean ops algorithms work by first deciding horizontal orderings for curves
  4  //! (with some guarantees about being approximately correct)
  5  //! and then laying out the curves in a way that has the horizontal orders we decided
  6  //! on. For this second step, which is the task of this module, we really
  7  //! guarantee correctness: the curves that we output really have the specified
  8  //! orders.
  9  
 10  use std::collections::BinaryHeap;
 11  
 12  use kurbo::{BezPath, CubicBez};
 13  
 14  use crate::{
 15      curve::{
 16          transverse::{transversal_after, transversal_before},
 17          y_subsegment, EstParab, Order,
 18      },
 19      num::CheapOrderedFloat,
 20      order::ComparisonCache,
 21      topology::{HalfOutputSegVec, OutputSegIdx, OutputSegVec, ScanLineOrder},
 22      SegIdx, Segment, Segments,
 23  };
 24  
 25  // TODO: maybe a Context type to keep the parameters in check?
 26  #[allow(clippy::too_many_arguments)]
 27  fn ordered_curves_all_close(
 28      segs: &Segments,
 29      order: &[SegIdx],
 30      output_order: &[OutputSegIdx],
 31      out: &mut OutputSegVec<(BezPath, Option<usize>)>,
 32      mut y0: f64,
 33      y1: f64,
 34      endpoints: &HalfOutputSegVec<kurbo::Point>,
 35      accuracy: f64,
 36  ) {
 37      if order.iter().all(|&s| segs[s].is_line()) {
 38          // If everything is a line, since our endpoints are already in the right order
 39          // the segments between them are also guaranteed to be in the right order.
 40          return;
 41      }
 42  
 43      if order.len() == 2 {
 44          let s0 = order[0];
 45          let s1 = order[1];
 46          let o0 = output_order[0];
 47          let o1 = output_order[1];
 48  
 49          let p0 = bez_end(&out[o0].0);
 50          let p1 = bez_end(&out[o1].0);
 51          let q0 = endpoints[o0.second_half()];
 52          let q1 = endpoints[o1.second_half()];
 53  
 54          if p0.y == y0 && p1.y == y0 {
 55              let c0 = next_subsegment(&segs[s0], &out[o0].0, y1, q0);
 56              let c1 = next_subsegment(&segs[s1], &out[o1].0, y1, q1);
 57  
 58              if transversal_after(c0, c1, y1) {
 59                  return;
 60              }
 61          }
 62  
 63          if q0.y == y1 && q1.y == y1 {
 64              let c0 = next_subsegment(&segs[s0], &out[o0].0, y1, q0);
 65              let c1 = next_subsegment(&segs[s1], &out[o1].0, y1, q1);
 66  
 67              if transversal_before(c0, c1, y0) {
 68                  return;
 69              }
 70          }
 71      }
 72  
 73      // Ensure everything in `out` goes up to `y0`. Anything that doesn't go up to `y0` is
 74      // an output where we can just copy from the input.
 75      for (&seg_idx, &out_idx) in order.iter().zip(output_order) {
 76          let (out_bez, out_copied_idx) = &mut out[out_idx];
 77          if bez_end_y(out_bez) < y0 {
 78              *out_copied_idx = Some(out_bez.elements().len() - 1);
 79              let endpoint = endpoints[out_idx.second_half()];
 80              if segs[seg_idx].is_line() {
 81                  out_bez.line_to(endpoint)
 82              } else {
 83                  let c = next_subsegment(&segs[seg_idx], out_bez, y0, endpoint);
 84                  out_bez.curve_to(c.p1, c.p2, c.p3);
 85              }
 86          }
 87      }
 88  
 89      let mut approxes = Vec::new();
 90  
 91      // The recentering helps, but there are still some issues with approximation when almost horizontal
 92      while y0 < y1 {
 93          let next_y0 = approximate(y0, y1, order, segs, accuracy, &mut approxes);
 94  
 95          let y_mid = (y0 + next_y0) / 2.0;
 96  
 97          let mut x_mid_max_so_far = f64::NEG_INFINITY;
 98          let mut x1_max_so_far = f64::NEG_INFINITY;
 99          for (quad, out_idx) in approxes.iter().zip(output_order) {
100              // These would be the control points, but since `approximate` recenters the y interval,
101              // they're different...
102              // let x_mid = quad.c0 + y_mid * quad.c1 + y0 * next_y0 * quad.c2;
103              // let x1 = quad.c0 + quad.c1 * next_y0 + quad.c2 * next_y0 * next_y0;
104              let x_mid = quad.c0 + (y0 - y_mid) * (next_y0 - y_mid) * quad.c2;
105              let x1 = quad.c0
106                  + quad.c1 * (next_y0 - y_mid)
107                  + quad.c2 * (next_y0 - y_mid) * (next_y0 - y_mid);
108  
109              x_mid_max_so_far = x_mid_max_so_far.max(x_mid);
110              x1_max_so_far = x1_max_so_far.max(x1);
111  
112              out[*out_idx]
113                  .0
114                  .quad_to((x_mid_max_so_far, y_mid), (x1_max_so_far, next_y0));
115          }
116          approxes.clear();
117  
118          y0 = next_y0;
119      }
120  }
121  
122  fn horizontal_error_weight(c: CubicBez) -> f64 {
123      let tangents = [c.p1 - c.p0, c.p2 - c.p1, c.p3 - c.p2];
124      let mut weight = 0.0f64;
125  
126      fn different_signs(x: f64, y: f64) -> bool {
127          (x > 0.0 && y < 0.0) || (x < 0.0) && (y > 0.0)
128      }
129  
130      if different_signs(tangents[0].x, tangents[1].x)
131          || different_signs(tangents[1].x, tangents[2].x)
132      {
133          return 1.0;
134      }
135  
136      for v in tangents {
137          let denom = v.hypot2().sqrt();
138          if denom != 0.0 {
139              let w = v.y.abs() / denom;
140              weight = weight.max(w);
141          }
142      }
143      weight
144  }
145  
146  // Pushes approximating parabolas into `out`, and returns the final y position.
147  fn approximate(
148      y0: f64,
149      mut y1: f64,
150      order: &[SegIdx],
151      segs: &Segments,
152      accuracy: f64,
153      out: &mut Vec<EstParab>,
154  ) -> f64 {
155      let orig_len = out.len();
156  
157      'retry: loop {
158          let y_mid = (y0 + y1) / 2.0;
159          for seg_idx in order {
160              let cubic = segs[*seg_idx].to_kurbo_cubic();
161              let mut cubic = y_subsegment(cubic, y0, y1);
162              cubic.p0.y -= y_mid;
163              cubic.p1.y -= y_mid;
164              cubic.p2.y -= y_mid;
165              cubic.p3.y -= y_mid;
166              let approx = EstParab::from_cubic(cubic);
167  
168              let factor = horizontal_error_weight(cubic);
169  
170              let too_short = y1 <= y0 + accuracy;
171              let accurate = approx.dmax.max(-approx.dmin) * factor < accuracy;
172              if !too_short && !accurate {
173                  // TODO: can we be smarter about this?
174                  y1 = (y0 + y1) / 2.0;
175                  out.truncate(orig_len);
176                  continue 'retry;
177              }
178  
179              out.push(approx);
180          }
181  
182          return y1;
183      }
184  }
185  
186  #[derive(Debug, PartialEq, Eq)]
187  struct HeapEntry {
188      y: CheapOrderedFloat,
189      idx: OutputSegIdx,
190  }
191  
192  impl Ord for HeapEntry {
193      fn cmp(&self, other: &Self) -> std::cmp::Ordering {
194          other.y.cmp(&self.y).then(self.idx.cmp(&other.idx))
195      }
196  }
197  
198  impl PartialOrd for HeapEntry {
199      fn partial_cmp(&self, other: &Self) -> Option<std::cmp::Ordering> {
200          Some(self.cmp(other))
201      }
202  }
203  
204  fn bez_end(path: &BezPath) -> kurbo::Point {
205      match path.elements().last().unwrap() {
206          kurbo::PathEl::MoveTo(p)
207          | kurbo::PathEl::LineTo(p)
208          | kurbo::PathEl::QuadTo(_, p)
209          | kurbo::PathEl::CurveTo(_, _, p) => *p,
210          kurbo::PathEl::ClosePath => unreachable!(),
211      }
212  }
213  
214  fn bez_end_y(path: &BezPath) -> f64 {
215      bez_end(path).y
216  }
217  
218  fn next_subsegment(seg: &Segment, out: &BezPath, y1: f64, endpoint: kurbo::Point) -> CubicBez {
219      let p0 = bez_end(out);
220      let mut c = y_subsegment(seg.to_kurbo_cubic(), p0.y, y1);
221      c.p0 = p0;
222      if endpoint.y == y1 {
223          c.p3 = endpoint;
224      }
225      c
226  }
227  
228  /// Compute positions for all of the output segments.
229  ///
230  /// The orders between the output segments is specified by `order`. The endpoints
231  /// should have been already computed (in a way that satisfies the order), and
232  /// are provided in `endpoints`. For each output segment, we return a Bézier
233  /// path.
234  ///
235  /// The `usize` return value tells which segment (if any) in the returned
236  /// path was the one that was "far" from any other paths. This is really
237  /// only interesting for diagnosis/visualization so the API should probably
238  /// be refined somehow to make it optional. (TODO)
239  pub(crate) fn compute_positions(
240      segs: &Segments,
241      orig_seg_map: &OutputSegVec<SegIdx>,
242      cmp: &mut ComparisonCache,
243      endpoints: &HalfOutputSegVec<kurbo::Point>,
244      scan_order: &ScanLineOrder,
245      accuracy: f64,
246  ) -> OutputSegVec<(BezPath, Option<usize>)> {
247      let mut out = OutputSegVec::<(BezPath, Option<usize>)>::with_size(orig_seg_map.len());
248      // We try to build `out` lazily, by avoiding copying input segments to outputs segments
249      // until they're needed (by copying in one go, we avoid excess subdivisions). That means
250      // we need to separately keep track of how far down we've looked at each output. If
251      // we weren't lazy, that would just be the last `y` coordinate in `out`.
252      let mut last_y = OutputSegVec::<f64>::with_size(orig_seg_map.len());
253      let mut queue = BinaryHeap::<HeapEntry>::new();
254      for idx in out.indices() {
255          let p = endpoints[idx.first_half()];
256          let q = endpoints[idx.second_half()];
257          if p.y == q.y {
258              out[idx].0.move_to(p);
259              out[idx].0.line_to(q);
260              out[idx].1 = Some(0);
261              continue;
262          }
263          out[idx].0.move_to(p);
264          queue.push(HeapEntry { y: p.y.into(), idx });
265      }
266  
267      while let Some(entry) = queue.pop() {
268          // Maybe this entry was already handled by one of its neighbors, in which case we skip it.
269          let y0 = entry.y.into_inner();
270          if last_y[entry.idx] > y0 {
271              continue;
272          }
273  
274          let mut y1 = endpoints[entry.idx.second_half()].y;
275          let mut west_scan = vec![entry.idx];
276          let mut cur = entry.idx;
277          while let Some(nbr) = scan_order.west_neighbor_after(cur, y0) {
278              let order = cmp.compare_segments(segs, orig_seg_map[nbr], orig_seg_map[cur]);
279              let (_, cmp_end_y, cmp_order) = order.iter().find(|(_, end_y, _)| *end_y > y0).unwrap();
280  
281              if cmp_order == Order::Left {
282                  let next_close_y = scan_order
283                      .close_west_neighbor_height_after(cur, y0, orig_seg_map, segs, cmp)
284                      .unwrap_or(f64::INFINITY);
285                  y1 = y1.min(next_close_y.max(cmp_end_y));
286                  break;
287              } else {
288                  y1 = y1.min(cmp_end_y).min(endpoints[nbr.second_half()].y);
289                  west_scan.push(nbr);
290              }
291              cur = nbr;
292          }
293  
294          let mut east_scan = vec![];
295          let mut cur = entry.idx;
296          while let Some(nbr) = scan_order.east_neighbor_after(cur, y0) {
297              let order = cmp.compare_segments(segs, orig_seg_map[cur], orig_seg_map[nbr]);
298              let (_, cmp_end_y, cmp_order) = order.iter().find(|(_, end_y, _)| *end_y > y0).unwrap();
299  
300              if cmp_order == Order::Left {
301                  let next_close_y = scan_order
302                      .close_east_neighbor_height_after(cur, y0, orig_seg_map, segs, cmp)
303                      .unwrap_or(f64::INFINITY);
304                  y1 = y1.min(next_close_y.max(cmp_end_y));
305                  break;
306              } else {
307                  y1 = y1.min(cmp_end_y).min(endpoints[nbr.second_half()].y);
308                  east_scan.push(nbr);
309              }
310              cur = nbr;
311          }
312  
313          let mut neighbors = west_scan;
314          neighbors.reverse();
315          neighbors.extend(east_scan);
316          if neighbors.len() == 1 {
317              let idx = entry.idx;
318              // We're far from everything, so we can just copy the input bezier to the output.
319              // We only actually do this eagerly for horizontal segments: for other segments
320              // we'll hold off on the copying because it might allow us to avoid further
321              // subdivision.
322              if y0 == y1 {
323                  out[idx].1 = Some(out[idx].0.elements().len() - 1);
324                  out[idx].0.line_to(endpoints[idx.second_half()]);
325              }
326          } else {
327              let orig_neighbors = neighbors
328                  .iter()
329                  .map(|idx| orig_seg_map[*idx])
330                  .collect::<Vec<_>>();
331  
332              ordered_curves_all_close(
333                  segs,
334                  &orig_neighbors,
335                  &neighbors,
336                  &mut out,
337                  y0,
338                  y1,
339                  endpoints,
340                  accuracy,
341              );
342          }
343          for idx in neighbors {
344              let y_end = endpoints[idx.second_half()].y;
345              if y1 < y_end {
346                  queue.push(HeapEntry { y: y1.into(), idx });
347              }
348              last_y[idx] = y1;
349          }
350      }
351  
352      for out_idx in out.indices() {
353          let (out_bez, out_copied_idx) = &mut out[out_idx];
354          let endpoint = endpoints[out_idx.second_half()];
355          if bez_end_y(out_bez) < endpoint.y {
356              let seg_idx = orig_seg_map[out_idx];
357              *out_copied_idx = Some(out_bez.elements().len() - 1);
358              if segs[seg_idx].is_line() {
359                  out_bez.line_to(endpoint)
360              } else {
361                  let c = next_subsegment(&segs[seg_idx], out_bez, endpoint.y, endpoint);
362                  out_bez.curve_to(c.p1, c.p2, c.p3);
363              }
364          } else {
365              // The quadratic approximations don't respect the fixed endpoints, so tidy them
366              // up. Since both the quadratic approximations and the endpoints satisfy
367              // the ordering, this doesn't mess up the ordering.
368              match out[out_idx].0.elements_mut().last_mut().unwrap() {
369                  kurbo::PathEl::MoveTo(p)
370                  | kurbo::PathEl::LineTo(p)
371                  | kurbo::PathEl::QuadTo(_, p)
372                  | kurbo::PathEl::CurveTo(_, _, p) => {
373                      *p = endpoints[out_idx.second_half()];
374                  }
375                  kurbo::PathEl::ClosePath => unreachable!(),
376              }
377          }
378      }
379      out
380  }
381  
382  #[cfg(test)]
383  mod tests {
384      use kurbo::{BezPath, CubicBez, Line, PathSeg};
385  
386      use crate::Segments;
387  
388      use super::horizontal_error_weight;
389  
390      #[test]
391      fn error_weight() {
392          // Vertical lines have a weight of 1.
393          let c: CubicBez = PathSeg::from(Line::new((0.0, 0.0), (0.0, 1.0))).to_cubic();
394          assert_eq!(horizontal_error_weight(c), 1.0);
395  
396          // Horizontal lines have a weight of 0.
397          let c: CubicBez = PathSeg::from(Line::new((1.0, 0.0), (0.0, 0.0))).to_cubic();
398          assert_eq!(horizontal_error_weight(c), 0.0);
399  
400          let c: CubicBez = PathSeg::from(Line::new((0.0, 0.0), (1.0, 0.0))).to_cubic();
401          assert_eq!(horizontal_error_weight(c), 0.0);
402  
403          // S-shaped curves have a weight of 1
404          let c = CubicBez::new((0.0, 0.0), (1.0, 0.0), (-1.0, 1.0), (0.0, 1.0));
405          assert_eq!(horizontal_error_weight(c), 1.0);
406  
407          // Diagonal lines have a weight of about 1/sqrt(2).
408          let c: CubicBez = PathSeg::from(Line::new((0.0, 0.0), (1.0, 1.0))).to_cubic();
409          assert_eq!(horizontal_error_weight(c), 1.0 / 2.0f64.sqrt());
410      }
411  
412      // Here's an example of lines where we used to have a bad
413      // approximation. Check that they only require a single
414      // approximation step.
415      #[test]
416      fn linear_approx() {
417          let mut p0 = BezPath::new();
418          p0.move_to((-212.00062561035156, 90.03582000732422));
419          p0.line_to((211.99937438964844, 90.03646087646484));
420          p0.close_path();
421  
422          let mut p1 = BezPath::new();
423          p1.move_to((211.99964904785153, -90.0358200073242));
424          p1.line_to((211.99937438964844, 90.03646087646484));
425          p1.close_path();
426  
427          let mut segs = Segments::default();
428          segs.add_bez_path(&p0).unwrap();
429          segs.add_bez_path(&p1).unwrap();
430          let s0 = segs.indices().next().unwrap();
431          let s1 = segs.indices().nth(2).unwrap();
432  
433          let mut out = Vec::new();
434          let y1 = 90.03646087646484;
435          let next_y0 = super::approximate(90.03645987646233, y1, &[s0, s1], &segs, 2.5e-7, &mut out);
436          assert_eq!(next_y0, y1);
437      }
438  }