Subdivision methods for implicit curves in 2D

In this section, we describe a subdivision method, for computing the topology of a curve $\mathcal{C}$, defined by an implicit equation $f (x, y) = 0$, where $f \in \mathbbm{Q}[x, y]$. It uses the representation of polynomials in the Bernstein basis of a domain $[a, b] \times [c, d]$:

\[ f (x, y) = \sum_{0 \leqslant i \leqslant d_1, 0 \leqslant j \leqslant d_2} c_{i, j} B^i_{d_1} (x ; a, b) B^j_{d_2} (y ; c, d), \]

where $(B^i_d (x ; a, b))_{i = 0, \ldots, d}$ is the Bernstein basis in degree $d$, on the interval $[a, b]$. The polynomial $f (x, y)$ will be represented by its control coefficients $\mathbf{c}= (c_{i, j})_{0 \leqslant i \leqslant d_1, 0 \leqslant j \leqslant d_2}$ and the domain $[a, b] \times [c, d]$. It can be subdivided by applying the de Casteljau algorithm, in the $x$ and $y$ direction.

The subdivision algorithm is based on a criterion for detecting the regularity of the curve $\mathcal{C}$, in the domain $[a, b] \times [c, d]$.

Definition: We say that a domain $D = [a, b] \times [c, d]$ is regular for $\mathcal{C}$, if its topology in $D$ is uniquely determined from the intersection points of $\mathcal{C}$, with the boundary $\partial D$ of $D$.

Here is a property that we use to detect regular domains:

Proposition: If $\partial_y f (x, y) \neq 0$ (resp. $\partial_x f (x, y) \neq 0$) in $D = [a, b] \times [c, d]$, then $D$ is regular for $\mathcal{C}$ .

The test is implemented by checking that the control coefficients of $\partial_y f$ or $\partial_y f$ on D have a constant sign, which implies the previous proposition.

Algorithm: Subdivision algorithm for implicit curves.

input: $f (x, y) \in \mathbbm{Q}[x, y]$ defining the implicit curves $C$ and a domain $D = [a, b] \times [c, d] \subset \mathbbm{R}^2$.

  • Initialize $\mathcal{L}$ with $D$;
  • while $\mathcal{L} \neq \emptyset$;
  • pick a domain $D$ in $\mathcal{L}$ and test if it is regular;
  • if $D$ is not regular not subdivide it and add the subdomains to $\mathcal{L}$;
  • otherwise, connect the points on $\partial D \cap \mathcal{C}$ to get the topology of $\mathcal{C}$ in $D$;
output: a graph of points connected by segments, isotopic to $\mathcal{C}$.


void topology::assign(point_graph<C> & g, MPol<QQ>& p, 
           TopSbdBdg2d<C> mth, 
           C x0, C x1, C y0, C y1);

See also:
synaps/topology/TopSbdBdg2d.h, TopSbdBdg.



#include <synaps/topology/TopSbdBdg2d.h>

int main(int argc, char** argv) 
  MPol<QQ> p("x^6+3*x^4*y^2+3*x^2*y^4+y^6-4x^2*y^2");
  topology::point_graph<double> g;
  topology::assign(g, p, TopSbdBdg2d<double>(), 
                   -2., 2.1, -2.2, 2., 2);

Here is a picture of such output: