How to resolve the algorithm Bézier curves/Intersections step by step in the D programming language

Published on 12 May 2024 09:40 PM
#D

How to resolve the algorithm Bézier curves/Intersections step by step in the D programming language

Table of Contents

Problem Statement

You are given two planar quadratic Bézier curves, having control points

( − 1 , 0 ) , ( 0 , 10 ) , ( 1 , 0 )

{\displaystyle (-1,0),(0,10),(1,0)}

and

( 2 , 1 ) , ( − 8 , 2 ) , ( 2 , 3 )

{\displaystyle (2,1),(-8,2),(2,3)}

, respectively. They are parabolas intersecting at four points, as shown in the following diagram:

The task is to write a program that finds all four intersection points and prints their

( x , y )

{\displaystyle (x,y)}

coordinates. You may use any algorithm you know of or can think of, including any of those that others have used.

Let's start with the solution:

Step by Step solution about How to resolve the algorithm Bézier curves/Intersections step by step in the D programming language

Source code in the d programming language

// The control points of a planar quadratic Bézier curve form a
// triangle--called the "control polygon"--that completely contains
// the curve. Furthermore, the rectangle formed by the minimum and
// maximum x and y values of the control polygon completely contain
// the polygon, and therefore also the curve.
//
// Thus a simple method for narrowing down where intersections might
// be is: subdivide both curves until you find "small enough" regions
// where these rectangles overlap.

import std.algorithm;
import std.container.slist;
import std.math;
import std.range;
import std.stdio;

struct point
{
  double x, y;
}

struct quadratic_spline         // Non-parametric spline.
{
  double c0, c1, c2;
}

struct quadratic_curve          // Planar parametric spline.
{
  quadratic_spline x, y;
}

void
subdivide_quadratic_spline (quadratic_spline q, double t,
                            ref quadratic_spline u,
                            ref quadratic_spline v)
{
  // Subdivision by de Casteljau's algorithm.
  immutable s = 1 - t;
  immutable c0 = q.c0;
  immutable c1 = q.c1;
  immutable c2 = q.c2;
  u.c0 = c0;
  v.c2 = c2;
  u.c1 = (s * c0) + (t * c1);
  v.c1 = (s * c1) + (t * c2);
  u.c2 = (s * u.c1) + (t * v.c1);
  v.c0 = u.c2;
}

void
subdivide_quadratic_curve (quadratic_curve q, double t,
                            ref quadratic_curve u,
                            ref quadratic_curve v)
{
  subdivide_quadratic_spline (q.x, t, u.x, v.x);
  subdivide_quadratic_spline (q.y, t, u.y, v.y);
}

bool
rectangles_overlap (double xa0, double ya0, double xa1, double ya1,
                    double xb0, double yb0, double xb1, double yb1)
{
  // It is assumed that xa0<=xa1, ya0<=ya1, xb0<=xb1, and yb0<=yb1.
  return (xb0 <= xa1 && xa0 <= xb1 && yb0 <= ya1 && ya0 <= yb1);
}

void
test_intersection (quadratic_curve p, quadratic_curve q, double tol,
                   ref bool exclude, ref bool accept,
                   ref point intersection)
{
  // I will not do a lot of checking for intersections, as one might
  // wish to do in a particular application. If the boxes are small
  // enough, I will accept the point as an intersection.

  immutable pxmin = min (p.x.c0, p.x.c1, p.x.c2);
  immutable pymin = min (p.y.c0, p.y.c1, p.y.c2);
  immutable pxmax = max (p.x.c0, p.x.c1, p.x.c2);
  immutable pymax = max (p.y.c0, p.y.c1, p.y.c2);

  immutable qxmin = min (q.x.c0, q.x.c1, q.x.c2);
  immutable qymin = min (q.y.c0, q.y.c1, q.y.c2);
  immutable qxmax = max (q.x.c0, q.x.c1, q.x.c2);
  immutable qymax = max (q.y.c0, q.y.c1, q.y.c2);

  exclude = true;
  accept = false;
  if (rectangles_overlap (pxmin, pymin, pxmax, pymax,
                          qxmin, qymin, qxmax, qymax))
    {
      exclude = false;
      immutable xmin = max (pxmin, qxmin);
      immutable xmax = min (pxmax, qxmax);
      assert (xmax >= xmin);
      if (xmax - xmin <= tol)
        {
          immutable ymin = max (pymin, qymin);
          immutable ymax = min (pymax, qymax);
          assert (ymax >= ymin);
          if (ymax - ymin <= tol)
            {
              accept = true;
              intersection = point ((0.5 * xmin) + (0.5 * xmax),
                                    (0.5 * ymin) + (0.5 * ymax));
            }
        }
    }
}

bool
seems_to_be_a_duplicate (point[] intersections, point xy,
                         double spacing)
{
  bool seems_to_be_dup = false;
  int i = 0;
  while (!seems_to_be_dup && i != intersections.length)
    {
      immutable pt = intersections[i];
      seems_to_be_dup =
        fabs (pt.x - xy.x) < spacing && fabs (pt.y - xy.y) < spacing;
      i += 1;
    }
  return seems_to_be_dup;
}

point[]
find_intersections (quadratic_curve p, quadratic_curve q,
                    double tol, double spacing)
{
  point[] intersections;
  int num_intersections = 0;

  struct workset
  {
    quadratic_curve p, q;
  }
  SList!workset workload;

  // Initial workload.
  workload.insertFront(workset (p, q));

  // Quit looking after having /*found four intersections*/ or emptied
  // the workload.
  while (/*num_intersections != 4 &&*/ !workload.empty)
    {
      auto work = workload.front;
      workload.removeFront();

      bool exclude;
      bool accept;
      point intersection;
      test_intersection (work.p, work.q, tol, exclude, accept,
                         intersection);
      if (accept)
        {
          // This is a crude way to avoid detecting the same
          // intersection twice: require some space between
          // intersections. For, say, font design work, this method
          // should be adequate.
          if (!seems_to_be_a_duplicate (intersections,
                                        intersection, spacing))
            {
              intersections.length = num_intersections + 1;
              intersections[num_intersections] = intersection;
              num_intersections += 1;
            }
        }
      else if (!exclude)
        {
          quadratic_curve p0, p1, q0, q1;
          subdivide_quadratic_curve (work.p, 0.5, p0, p1);
          subdivide_quadratic_curve (work.q, 0.5, q0, q1);
          workload.insertFront(workset (p0, q0));
          workload.insertFront(workset (p0, q1));
          workload.insertFront(workset (p1, q0));
          workload.insertFront(workset (p1, q1));
        }
    }

  return intersections;
}

int
main ()
{
  quadratic_curve p, q;
  p.x.c0 = -1.0;  p.x.c1 =  0.0;  p.x.c2 =  1.0;
  p.y.c0 =  0.0;  p.y.c1 = 10.0;  p.y.c2 =  0.0;
  q.x.c0 =  2.0;  q.x.c1 = -8.0;  q.x.c2 =  2.0;
  q.y.c0 =  1.0;  q.y.c1 =  2.0;  q.y.c2 =  3.0;

  immutable tol = 0.0000001;
  immutable spacing = 10 * tol;

  auto intersections = find_intersections (p, q, tol, spacing);
  for (int i = 0; i != intersections.length; i += 1)
    printf("(%f, %f)\n", intersections[i].x, intersections[i].y);

  return 0;
}


  

You may also check:How to resolve the algorithm Hello world/Text step by step in the Prolog programming language
You may also check:How to resolve the algorithm Zero to the zero power step by step in the Neko programming language
You may also check:How to resolve the algorithm Matrix multiplication step by step in the OxygenBasic programming language
You may also check:How to resolve the algorithm Text processing/Max licenses in use step by step in the Perl programming language
You may also check:How to resolve the algorithm Heronian triangles step by step in the EchoLisp programming language