SYNAPS (SYmbolic Numeric ApplicationS) » Geometry

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}$ .

Implementation

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

This function computes a graph of points g of type topology::point_graph<C>, which is isotopic to the curve defined by the polynomial p, in the box $B = [x_0, x_1] \times [y_0, y_1]$ . The coordinates of the points in g are of type C. The implementation is designed for C=double.

The input coefficients of p are rational numbers. During the computation, they are rounded to double number types. If the topology in the box can be certified from this approximation, a graph of points connecting the intersection of the curve with the boundary of is computed. Otherwise the box is subdivided into subboxes and the test is applied recursively.

If the curve has singular points, the subdivision will continue until the size of the box is smaller than the precision $\varepsilon$ specified by the class TopSbdBdg2d<C>. The default value is $10^{- 3}$ .

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

Example

#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:

SYNAPS DOCUMENTATION