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<PositionedOutputSeg>,
 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].path);
 50          let p1 = bez_end(&out[o1].path);
 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].path, y1, q0);
 56              let c1 = next_subsegment(&segs[s1], &out[o1].path, 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].path, y1, q0);
 65              let c1 = next_subsegment(&segs[s1], &out[o1].path, 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 = &mut out[out_idx];
 77          if bez_end_y(&out.path) < y0 {
 78              out.copied_idx = Some(out.path.elements().len() - 1);
 79              let endpoint = endpoints[out_idx.second_half()];
 80              if segs[seg_idx].is_line() {
 81                  out.path.line_to(endpoint)
 82              } else {
 83                  let c = next_subsegment(&segs[seg_idx], &out.path, y0, endpoint);
 84                  out.path.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                  .path
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  /// A positioned output segment, possibly perturbed from the original.
229  ///
230  /// We perturb output segments in order to strictly enforce ordering.
231  /// A single segment might be expanded into a longer path.
232  #[derive(Default)]
233  pub struct PositionedOutputSeg {
234      /// The path, which is guaranteed to start and end at the same points as the
235      /// output segment did.
236      pub path: BezPath,
237      /// Sometimes (often, even) a large part of the output segment didn't need
238      /// to be perturbed. Our sweep-line invariants guarantee that unperturbed
239      /// part can only be a single contiguous interval: this stores the index
240      /// of the unperturbed part within `path.segments()`.
241      pub copied_idx: Option<usize>,
242  }
243  
244  /// Compute positions for all of the output segments.
245  ///
246  /// The orders between the output segments is specified by `order`. The endpoints
247  /// should have been already computed (in a way that satisfies the order), and
248  /// are provided in `endpoints`. For each output segment, we return a Bézier
249  /// path.
250  pub(crate) fn compute_positions(
251      segs: &Segments,
252      orig_seg_map: &OutputSegVec<SegIdx>,
253      cmp: &mut ComparisonCache,
254      endpoints: &HalfOutputSegVec<kurbo::Point>,
255      scan_order: &ScanLineOrder,
256      accuracy: f64,
257  ) -> OutputSegVec<PositionedOutputSeg> {
258      let mut out = OutputSegVec::<PositionedOutputSeg>::with_size(orig_seg_map.len());
259      // We try to build `out` lazily, by avoiding copying input segments to outputs segments
260      // until they're needed (by copying in one go, we avoid excess subdivisions). That means
261      // we need to separately keep track of how far down we've looked at each output. If
262      // we weren't lazy, that would just be the last `y` coordinate in `out`.
263      let mut last_y = OutputSegVec::<f64>::with_size(orig_seg_map.len());
264      let mut queue = BinaryHeap::<HeapEntry>::new();
265      for idx in out.indices() {
266          let p = endpoints[idx.first_half()];
267          let q = endpoints[idx.second_half()];
268          if p.y == q.y {
269              out[idx].path.move_to(p);
270              out[idx].path.line_to(q);
271              out[idx].copied_idx = Some(0);
272              continue;
273          }
274          out[idx].path.move_to(p);
275          queue.push(HeapEntry { y: p.y.into(), idx });
276      }
277  
278      while let Some(entry) = queue.pop() {
279          // Maybe this entry was already handled by one of its neighbors, in which case we skip it.
280          let y0 = entry.y.into_inner();
281          if last_y[entry.idx] > y0 {
282              continue;
283          }
284  
285          let mut y1 = endpoints[entry.idx.second_half()].y;
286          let mut west_scan = vec![entry.idx];
287          let mut cur = entry.idx;
288          while let Some(nbr) = scan_order.west_neighbor_after(cur, y0) {
289              let order = cmp.compare_segments(segs, orig_seg_map[nbr], orig_seg_map[cur]);
290              let (_, cmp_end_y, cmp_order) = order.iter().find(|(_, end_y, _)| *end_y > y0).unwrap();
291  
292              if cmp_order == Order::Left {
293                  let next_close_y = scan_order
294                      .close_west_neighbor_height_after(cur, y0, orig_seg_map, segs, cmp)
295                      .unwrap_or(f64::INFINITY);
296                  y1 = y1.min(next_close_y.max(cmp_end_y));
297                  break;
298              } else {
299                  y1 = y1.min(cmp_end_y).min(endpoints[nbr.second_half()].y);
300                  west_scan.push(nbr);
301              }
302              cur = nbr;
303          }
304  
305          let mut east_scan = vec![];
306          let mut cur = entry.idx;
307          while let Some(nbr) = scan_order.east_neighbor_after(cur, y0) {
308              let order = cmp.compare_segments(segs, orig_seg_map[cur], orig_seg_map[nbr]);
309              let (_, cmp_end_y, cmp_order) = order.iter().find(|(_, end_y, _)| *end_y > y0).unwrap();
310  
311              if cmp_order == Order::Left {
312                  let next_close_y = scan_order
313                      .close_east_neighbor_height_after(cur, y0, orig_seg_map, segs, cmp)
314                      .unwrap_or(f64::INFINITY);
315                  y1 = y1.min(next_close_y.max(cmp_end_y));
316                  break;
317              } else {
318                  y1 = y1.min(cmp_end_y).min(endpoints[nbr.second_half()].y);
319                  east_scan.push(nbr);
320              }
321              cur = nbr;
322          }
323  
324          let mut neighbors = west_scan;
325          neighbors.reverse();
326          neighbors.extend(east_scan);
327          if neighbors.len() == 1 {
328              let idx = entry.idx;
329              // We're far from everything, so we can just copy the input bezier to the output.
330              // We only actually do this eagerly for horizontal segments: for other segments
331              // we'll hold off on the copying because it might allow us to avoid further
332              // subdivision.
333              if y0 == y1 {
334                  out[idx].copied_idx = Some(out[idx].path.elements().len() - 1);
335                  out[idx].path.line_to(endpoints[idx.second_half()]);
336              }
337          } else {
338              let orig_neighbors = neighbors
339                  .iter()
340                  .map(|idx| orig_seg_map[*idx])
341                  .collect::<Vec<_>>();
342  
343              ordered_curves_all_close(
344                  segs,
345                  &orig_neighbors,
346                  &neighbors,
347                  &mut out,
348                  y0,
349                  y1,
350                  endpoints,
351                  accuracy,
352              );
353          }
354          for idx in neighbors {
355              let y_end = endpoints[idx.second_half()].y;
356              if y1 < y_end {
357                  queue.push(HeapEntry { y: y1.into(), idx });
358              }
359              last_y[idx] = y1;
360          }
361      }
362  
363      for out_idx in out.indices() {
364          let out = &mut out[out_idx];
365          let endpoint = endpoints[out_idx.second_half()];
366          if bez_end_y(&out.path) < endpoint.y {
367              let seg_idx = orig_seg_map[out_idx];
368              out.copied_idx = Some(out.path.elements().len() - 1);
369              if segs[seg_idx].is_line() {
370                  out.path.line_to(endpoint)
371              } else {
372                  let c = next_subsegment(&segs[seg_idx], &out.path, endpoint.y, endpoint);
373                  out.path.curve_to(c.p1, c.p2, c.p3);
374              }
375          } else {
376              // The quadratic approximations don't respect the fixed endpoints, so tidy them
377              // up. Since both the quadratic approximations and the endpoints satisfy
378              // the ordering, this doesn't mess up the ordering.
379              match out.path.elements_mut().last_mut().unwrap() {
380                  kurbo::PathEl::MoveTo(p)
381                  | kurbo::PathEl::LineTo(p)
382                  | kurbo::PathEl::QuadTo(_, p)
383                  | kurbo::PathEl::CurveTo(_, _, p) => {
384                      *p = endpoints[out_idx.second_half()];
385                  }
386                  kurbo::PathEl::ClosePath => unreachable!(),
387              }
388          }
389      }
390      out
391  }
392  
393  #[cfg(test)]
394  mod tests {
395      use kurbo::{BezPath, CubicBez, Line, PathSeg};
396  
397      use crate::Segments;
398  
399      use super::horizontal_error_weight;
400  
401      #[test]
402      fn error_weight() {
403          // Vertical lines have a weight of 1.
404          let c: CubicBez = PathSeg::from(Line::new((0.0, 0.0), (0.0, 1.0))).to_cubic();
405          assert_eq!(horizontal_error_weight(c), 1.0);
406  
407          // Horizontal lines have a weight of 0.
408          let c: CubicBez = PathSeg::from(Line::new((1.0, 0.0), (0.0, 0.0))).to_cubic();
409          assert_eq!(horizontal_error_weight(c), 0.0);
410  
411          let c: CubicBez = PathSeg::from(Line::new((0.0, 0.0), (1.0, 0.0))).to_cubic();
412          assert_eq!(horizontal_error_weight(c), 0.0);
413  
414          // S-shaped curves have a weight of 1
415          let c = CubicBez::new((0.0, 0.0), (1.0, 0.0), (-1.0, 1.0), (0.0, 1.0));
416          assert_eq!(horizontal_error_weight(c), 1.0);
417  
418          // Diagonal lines have a weight of about 1/sqrt(2).
419          let c: CubicBez = PathSeg::from(Line::new((0.0, 0.0), (1.0, 1.0))).to_cubic();
420          assert_eq!(horizontal_error_weight(c), 1.0 / 2.0f64.sqrt());
421      }
422  
423      // Here's an example of lines where we used to have a bad
424      // approximation. Check that they only require a single
425      // approximation step.
426      #[test]
427      fn linear_approx() {
428          let mut p0 = BezPath::new();
429          p0.move_to((-212.00062561035156, 90.03582000732422));
430          p0.line_to((211.99937438964844, 90.03646087646484));
431          p0.close_path();
432  
433          let mut p1 = BezPath::new();
434          p1.move_to((211.99964904785153, -90.0358200073242));
435          p1.line_to((211.99937438964844, 90.03646087646484));
436          p1.close_path();
437  
438          let mut segs = Segments::default();
439          segs.add_bez_path(&p0).unwrap();
440          segs.add_bez_path(&p1).unwrap();
441          let s0 = segs.indices().next().unwrap();
442          let s1 = segs.indices().nth(2).unwrap();
443  
444          let mut out = Vec::new();
445          let y1 = 90.03646087646484;
446          let next_y0 = super::approximate(90.03645987646233, y1, &[s0, s1], &segs, 2.5e-7, &mut out);
447          assert_eq!(next_y0, y1);
448      }
449  }