da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm/*
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * Implement the Bezier clipping algorithm for finding
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * Bezier curve intersection points and collinear normals
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm *
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * Authors:
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * Marco Cecchetti <mrcekets at gmail.com>
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm *
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * Copyright 2008 authors
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm *
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * This library is free software; you can redistribute it and/or
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * modify it either under the terms of the GNU Lesser General Public
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * License version 2.1 as published by the Free Software Foundation
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * (the "LGPL") or, at your option, under the terms of the Mozilla
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * Public License Version 1.1 (the "MPL"). If you do not alter this
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * notice, a recipient may use your version of this file under either
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * the MPL or the LGPL.
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm *
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * You should have received a copy of the LGPL along with this library
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * in the file COPYING-LGPL-2.1; if not, write to the Free Software
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * You should have received a copy of the MPL along with this library
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * in the file COPYING-MPL-1.1
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm *
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * The contents of this file are subject to the Mozilla Public License
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * Version 1.1 (the "License"); you may not use this file except in
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * compliance with the License. You may obtain a copy of the License at
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * http://www.mozilla.org/MPL/
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm *
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * OF ANY KIND, either express or implied. See the LGPL or the MPL for
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * the specific language governing rights and limitations.
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm */
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#include <2geom/basic-intersection.h>
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#include <2geom/choose.h>
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#include <2geom/point.h>
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#include <2geom/interval.h>
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#include <2geom/bezier.h>
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#include <2geom/numeric/matrix.h>
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński#include <2geom/convex-hull.h>
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński#include <2geom/line.h>
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#include <cassert>
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#include <vector>
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#include <algorithm>
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#include <utility>
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm//#include <iomanip>
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosińskiusing std::swap;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#define VERBOSE 0
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#define CHECK 0
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrmnamespace Geom {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrmnamespace detail { namespace bezier_clipping {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm////////////////////////////////////////////////////////////////////////////////
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm// for debugging
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm//
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrmvoid print(std::vector<Point> const& cp, const char* msg = "")
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm{
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::cerr << msg << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm for (size_t i = 0; i < cp.size(); ++i)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::cerr << i << " : " << cp[i] << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm}
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrmtemplate< class charT >
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrmstd::basic_ostream<charT> &
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrmoperator<< (std::basic_ostream<charT> & os, const Interval & I)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm{
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm os << "[" << I.min() << ", " << I.max() << "]";
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm return os;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm}
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrmdouble angle (std::vector<Point> const& A)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm{
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm size_t n = A.size() -1;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm double a = std::atan2(A[n][Y] - A[0][Y], A[n][X] - A[0][X]);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm return (180 * a / M_PI);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm}
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrmsize_t get_precision(Interval const& I)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm{
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm double d = I.extent();
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm double e = 0.1, p = 10;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm int n = 0;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm while (n < 16 && d < e)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm p *= 10;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm e = 1/p;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm ++n;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm return n;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm}
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrmvoid range_assertion(int k, int m, int n, const char* msg)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm{
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if ( k < m || k > n)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::cerr << "range assertion failed: \n"
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm << msg << std::endl
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm << "value: " << k
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm << " range: " << m << ", " << n << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm assert (k >= m && k <= n);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm}
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm////////////////////////////////////////////////////////////////////////////////
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm// numerical routines
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm/*
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * Compute the binomial coefficient (n, k)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm */
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrmdouble binomial(unsigned int n, unsigned int k)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm{
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm return choose<double>(n, k);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm}
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm/*
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * Compute the determinant of the 2x2 matrix with column the point P1, P2
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm */
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrmdouble det(Point const& P1, Point const& P2)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm{
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm return P1[X]*P2[Y] - P1[Y]*P2[X];
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm}
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm/*
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * Solve the linear system [P1,P2] * P = Q
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * in case there isn't exactly one solution the routine returns false
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm */
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrmbool solve(Point & P, Point const& P1, Point const& P2, Point const& Q)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm{
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm double d = det(P1, P2);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (d == 0) return false;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm d = 1 / d;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm P[X] = det(Q, P2) * d;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm P[Y] = det(P1, Q) * d;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm return true;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm}
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm////////////////////////////////////////////////////////////////////////////////
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm// interval routines
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm/*
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * Map the sub-interval I in [0,1] into the interval J and assign it to J
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm */
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrmvoid map_to(Interval & J, Interval const& I)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm{
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński J.setEnds(J.valueAt(I.min()), J.valueAt(I.max()));
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm}
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm////////////////////////////////////////////////////////////////////////////////
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm// bezier curve routines
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm/*
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * Return true if all the Bezier curve control points are near,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * false otherwise
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm */
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński// Bezier.isConstant(precision)
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosińskibool is_constant(std::vector<Point> const& A, double precision)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm{
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm for (unsigned int i = 1; i < A.size(); ++i)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if(!are_near(A[i], A[0], precision))
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm return false;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm return true;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm}
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm/*
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * Compute the hodograph of the bezier curve B and return it in D
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm */
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński// derivative(Bezier)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrmvoid derivative(std::vector<Point> & D, std::vector<Point> const& B)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm{
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm D.clear();
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm size_t sz = B.size();
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (sz == 0) return;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (sz == 1)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm D.resize(1, Point(0,0));
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm return;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm size_t n = sz-1;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm D.reserve(n);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm for (size_t i = 0; i < n; ++i)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm D.push_back(n*(B[i+1] - B[i]));
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm}
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm/*
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * Compute the hodograph of the Bezier curve B rotated of 90 degree
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * and return it in D; we have N(t) orthogonal to B(t) for any t
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm */
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński// rot90(derivative(Bezier))
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrmvoid normal(std::vector<Point> & N, std::vector<Point> const& B)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm{
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm derivative(N,B);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm for (size_t i = 0; i < N.size(); ++i)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm N[i] = rot90(N[i]);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm}
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm/*
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * Compute the portion of the Bezier curve "B" wrt the interval [0,t]
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm */
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński// portion(Bezier, 0, t)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrmvoid left_portion(Coord t, std::vector<Point> & B)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm{
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm size_t n = B.size();
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm for (size_t i = 1; i < n; ++i)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm for (size_t j = n-1; j > i-1 ; --j)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm B[j] = lerp(t, B[j-1], B[j]);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm}
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm/*
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * Compute the portion of the Bezier curve "B" wrt the interval [t,1]
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm */
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński// portion(Bezier, t, 1)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrmvoid right_portion(Coord t, std::vector<Point> & B)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm{
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm size_t n = B.size();
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm for (size_t i = 1; i < n; ++i)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm for (size_t j = 0; j < n-i; ++j)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm B[j] = lerp(t, B[j], B[j+1]);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm}
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm/*
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * Compute the portion of the Bezier curve "B" wrt the interval "I"
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm */
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński// portion(Bezier, I)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrmvoid portion (std::vector<Point> & B , Interval const& I)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm{
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (I.min() == 0)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (I.max() == 1) return;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm left_portion(I.max(), B);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm return;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm right_portion(I.min(), B);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (I.max() == 1) return;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm double t = I.extent() / (1 - I.min());
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm left_portion(t, B);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm}
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm////////////////////////////////////////////////////////////////////////////////
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm// tags
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrmstruct intersection_point_tag;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrmstruct collinear_normal_tag;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrmtemplate <typename Tag>
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof KosińskiOptInterval clip(std::vector<Point> const& A,
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński std::vector<Point> const& B,
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński double precision);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrmtemplate <typename Tag>
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrmvoid iterate(std::vector<Interval>& domsA,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::vector<Interval>& domsB,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::vector<Point> const& A,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::vector<Point> const& B,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm Interval const& domA,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm Interval const& domB,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm double precision );
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm////////////////////////////////////////////////////////////////////////////////
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm// intersection
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm/*
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * Make up an orientation line using the control points c[i] and c[j]
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * the line is returned in the output parameter "l" in the form of a 3 element
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * vector : l[0] * x + l[1] * y + l[2] == 0; the line is normalized.
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm */
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński// Line(c[i], c[j])
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrmvoid orientation_line (std::vector<double> & l,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::vector<Point> const& c,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm size_t i, size_t j)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm{
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm l[0] = c[j][Y] - c[i][Y];
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm l[1] = c[i][X] - c[j][X];
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński l[2] = cross(c[j], c[i]);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm double length = std::sqrt(l[0] * l[0] + l[1] * l[1]);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm assert (length != 0);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm l[0] /= length;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm l[1] /= length;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm l[2] /= length;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm}
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm/*
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * Pick up an orientation line for the Bezier curve "c" and return it in
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * the output parameter "l"
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm */
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof KosińskiLine pick_orientation_line (std::vector<Point> const &c, double precision)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm{
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm size_t i = c.size();
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński while (--i > 0 && are_near(c[0], c[i], precision))
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {}
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński // this should never happen because when a new curve portion is created
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński // we check that it is not constant;
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński // however this requires that the precision used in the is_constant
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński // routine has to be the same used here in the are_near test
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński assert(i != 0);
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński Line line(c[0], c[i]);
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński return line;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm //std::cerr << "i = " << i << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm}
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm/*
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * Make up an orientation line for constant bezier curve;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * the orientation line is made up orthogonal to the other curve base line;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * the line is returned in the output parameter "l" in the form of a 3 element
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * vector : l[0] * x + l[1] * y + l[2] == 0; the line is normalized.
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm */
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof KosińskiLine orthogonal_orientation_line (std::vector<Point> const &c,
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński Point const &p,
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński double precision)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm{
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński // this should never happen
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński assert(!is_constant(c, precision));
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński Line line(p, (c.back() - c.front()).cw() + p);
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński return line;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm}
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm/*
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * Compute the signed distance of the point "P" from the normalized line l
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm */
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosińskidouble signed_distance(Point const &p, Line const &l)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm{
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński Coord a, b, c;
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński l.coefficients(a, b, c);
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński return a * p[X] + b * p[Y] + c;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm}
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm/*
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * Compute the min and max distance of the control points of the Bezier
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * curve "c" from the normalized orientation line "l".
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * This bounds are returned through the output Interval parameter"bound".
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm */
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof KosińskiInterval fat_line_bounds (std::vector<Point> const &c,
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński Line const &l)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm{
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński Interval bound(0, 0);
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński for (size_t i = 0; i < c.size(); ++i) {
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński bound.expandTo(signed_distance(c[i], l));
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński return bound;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm}
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm/*
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * return the x component of the intersection point between the line
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * passing through points p1, p2 and the line Y = "y"
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm */
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrmdouble intersect (Point const& p1, Point const& p2, double y)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm{
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm // we are sure that p2[Y] != p1[Y] because this routine is called
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm // only when the lower or the upper bound is crossed
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm double dy = (p2[Y] - p1[Y]);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm double s = (y - p1[Y]) / dy;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm return (p2[X]-p1[X])*s + p1[X];
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm}
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm/*
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * Clip the Bezier curve "B" wrt the fat line defined by the orientation
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * line "l" and the interval range "bound", the new parameter interval for
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * the clipped curve is returned through the output parameter "dom"
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm */
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof KosińskiOptInterval clip_interval (std::vector<Point> const& B,
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński Line const &l,
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński Interval const &bound)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm{
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm double n = B.size() - 1; // number of sub-intervals
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::vector<Point> D; // distance curve control points
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm D.reserve (B.size());
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm for (size_t i = 0; i < B.size(); ++i)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński const double d = signed_distance(B[i], l);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm D.push_back (Point(i/n, d));
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm //print(D);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński ConvexHull p;
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński p.swap(D);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm //print(p);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm bool plower, phigher;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm bool clower, chigher;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm double t, tmin = 1, tmax = 0;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm// std::cerr << "bound : " << bound << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm plower = (p[0][Y] < bound.min());
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm phigher = (p[0][Y] > bound.max());
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (!(plower || phigher)) // inside the fat line
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (tmin > p[0][X]) tmin = p[0][X];
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (tmax < p[0][X]) tmax = p[0][X];
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm// std::cerr << "0 : inside " << p[0]
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm// << " : tmin = " << tmin << ", tmax = " << tmax << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm for (size_t i = 1; i < p.size(); ++i)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm clower = (p[i][Y] < bound.min());
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm chigher = (p[i][Y] > bound.max());
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (!(clower || chigher)) // inside the fat line
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (tmin > p[i][X]) tmin = p[i][X];
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (tmax < p[i][X]) tmax = p[i][X];
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm// std::cerr << i << " : inside " << p[i]
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm// << " : tmin = " << tmin << ", tmax = " << tmax
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm// << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (clower != plower) // cross the lower bound
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm t = intersect(p[i-1], p[i], bound.min());
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (tmin > t) tmin = t;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (tmax < t) tmax = t;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm plower = clower;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm// std::cerr << i << " : lower " << p[i]
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm// << " : tmin = " << tmin << ", tmax = " << tmax
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm// << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (chigher != phigher) // cross the upper bound
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm t = intersect(p[i-1], p[i], bound.max());
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (tmin > t) tmin = t;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (tmax < t) tmax = t;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm phigher = chigher;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm// std::cerr << i << " : higher " << p[i]
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm// << " : tmin = " << tmin << ", tmax = " << tmax
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm// << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm // we have to test the closing segment for intersection
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm size_t last = p.size() - 1;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm clower = (p[0][Y] < bound.min());
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm chigher = (p[0][Y] > bound.max());
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (clower != plower) // cross the lower bound
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm t = intersect(p[last], p[0], bound.min());
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (tmin > t) tmin = t;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (tmax < t) tmax = t;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm// std::cerr << "0 : lower " << p[0]
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm// << " : tmin = " << tmin << ", tmax = " << tmax << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (chigher != phigher) // cross the upper bound
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm t = intersect(p[last], p[0], bound.max());
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (tmin > t) tmin = t;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (tmax < t) tmax = t;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm// std::cerr << "0 : higher " << p[0]
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm// << " : tmin = " << tmin << ", tmax = " << tmax << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński if (tmin == 1 && tmax == 0) {
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński return OptInterval();
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński } else {
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński return Interval(tmin, tmax);
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm}
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm/*
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * Clip the Bezier curve "B" wrt the Bezier curve "A" for individuating
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * intersection points the new parameter interval for the clipped curve
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * is returned through the output parameter "dom"
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm */
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrmtemplate <>
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof KosińskiOptInterval clip<intersection_point_tag> (std::vector<Point> const& A,
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński std::vector<Point> const& B,
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński double precision)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm{
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński Line bl;
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński if (is_constant(A, precision)) {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm Point M = middle_point(A.front(), A.back());
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński bl = orthogonal_orientation_line(B, M, precision);
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński } else {
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński bl = pick_orientation_line(A, precision);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński bl.normalize();
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński Interval bound = fat_line_bounds(A, bl);
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński return clip_interval(B, bl, bound);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm}
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm///////////////////////////////////////////////////////////////////////////////
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm// collinear normal
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm/*
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * Compute a closed focus for the Bezier curve B and return it in F
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * A focus is any curve through which all lines perpendicular to B(t) pass.
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm */
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrmvoid make_focus (std::vector<Point> & F, std::vector<Point> const& B)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm{
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm assert (B.size() > 2);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm size_t n = B.size() - 1;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm normal(F, B);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm Point c(1, 1);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#if VERBOSE
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (!solve(c, F[0], -F[n-1], B[n]-B[0]))
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::cerr << "make_focus: unable to make up a closed focus" << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#else
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm solve(c, F[0], -F[n-1], B[n]-B[0]);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#endif
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm// std::cerr << "c = " << c << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm // B(t) + c(t) * N(t)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm double n_inv = 1 / (double)(n);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm Point c0ni;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm F.push_back(c[1] * F[n-1]);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm F[n] += B[n];
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm for (size_t i = n-1; i > 0; --i)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm F[i] *= -c[0];
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm c0ni = F[i];
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm F[i] += (c[1] * F[i-1]);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm F[i] *= (i * n_inv);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm F[i] -= c0ni;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm F[i] += B[i];
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm F[0] *= c[0];
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm F[0] += B[0];
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm}
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm/*
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * Compute the projection on the plane (t, d) of the control points
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * (t, u, D(t,u)) where D(t,u) = <(B(t) - F(u)), B'(t)> with 0 <= t, u <= 1
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * B is a Bezier curve and F is a focus of another Bezier curve.
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * See Sederberg, Nishita, 1990 - Curve intersection using Bezier clipping.
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm */
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrmvoid distance_control_points (std::vector<Point> & D,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::vector<Point> const& B,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::vector<Point> const& F)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm{
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm assert (B.size() > 1);
ea8dd7683dd12883474f6cf9b5f424f8ed831166Kris assert (!F.empty());
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm const size_t n = B.size() - 1;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm const size_t m = F.size() - 1;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm const size_t r = 2 * n - 1;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm const double r_inv = 1 / (double)(r);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm D.clear();
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm D.reserve (B.size() * F.size());
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::vector<Point> dB;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm dB.reserve(n);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm for (size_t k = 0; k < n; ++k)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm dB.push_back (B[k+1] - B[k]);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm NL::Matrix dBB(n,B.size());
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm for (size_t i = 0; i < n; ++i)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm for (size_t j = 0; j < B.size(); ++j)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm dBB(i,j) = dot (dB[i], B[j]);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm NL::Matrix dBF(n, F.size());
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm for (size_t i = 0; i < n; ++i)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm for (size_t j = 0; j < F.size(); ++j)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm dBF(i,j) = dot (dB[i], F[j]);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
c8589a6c7367d09fa756755cef0dd448c7328a71Johan B. C. Engelen size_t l;
c8589a6c7367d09fa756755cef0dd448c7328a71Johan B. C. Engelen double bc;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm Point dij;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::vector<double> d(F.size());
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm for (size_t i = 0; i <= r; ++i)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm for (size_t j = 0; j <= m; ++j)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm d[j] = 0;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
c8589a6c7367d09fa756755cef0dd448c7328a71Johan B. C. Engelen const size_t k0 = std::max(i, n) - n;
c8589a6c7367d09fa756755cef0dd448c7328a71Johan B. C. Engelen const size_t kn = std::min(i, n-1);
c8589a6c7367d09fa756755cef0dd448c7328a71Johan B. C. Engelen const double bri = n / binomial(r,i);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm for (size_t k = k0; k <= kn; ++k)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm //if (k > i || (i-k) > n) continue;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm l = i - k;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#if CHECK
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm assert (l <= n);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#endif
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm bc = bri * binomial(n,l) * binomial(n-1, k);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm for (size_t j = 0; j <= m; ++j)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm //d[j] += bc * dot(dB[k], B[l] - F[j]);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm d[j] += bc * (dBB(k,l) - dBF(k,j));
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm double dmin, dmax;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm dmin = dmax = d[m];
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm for (size_t j = 0; j < m; ++j)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (dmin > d[j]) dmin = d[j];
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (dmax < d[j]) dmax = d[j];
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm dij[0] = i * r_inv;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm dij[1] = dmin;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm D.push_back (dij);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm dij[1] = dmax;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm D.push_back (dij);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm}
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm/*
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * Clip the Bezier curve "B" wrt the focus "F"; the new parameter interval for
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * the clipped curve is returned through the output parameter "dom"
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm */
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof KosińskiOptInterval clip_interval (std::vector<Point> const& B,
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński std::vector<Point> const& F)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm{
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::vector<Point> D; // distance curve control points
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm distance_control_points(D, B, F);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm //print(D, "D");
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm// ConvexHull chD(D);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm// std::vector<Point>& p = chD.boundary; // convex hull vertices
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński ConvexHull p;
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński p.swap(D);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm //print(p, "CH(D)");
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm bool plower, clower;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm double t, tmin = 1, tmax = 0;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm plower = (p[0][Y] < 0);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (p[0][Y] == 0) // on the x axis
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (tmin > p[0][X]) tmin = p[0][X];
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (tmax < p[0][X]) tmax = p[0][X];
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm// std::cerr << "0 : on x axis " << p[0]
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm// << " : tmin = " << tmin << ", tmax = " << tmax << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm for (size_t i = 1; i < p.size(); ++i)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm clower = (p[i][Y] < 0);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (p[i][Y] == 0) // on x axis
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (tmin > p[i][X]) tmin = p[i][X];
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (tmax < p[i][X]) tmax = p[i][X];
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm// std::cerr << i << " : on x axis " << p[i]
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm// << " : tmin = " << tmin << ", tmax = " << tmax
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm// << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm else if (clower != plower) // cross the x axis
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm t = intersect(p[i-1], p[i], 0);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (tmin > t) tmin = t;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (tmax < t) tmax = t;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm plower = clower;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm// std::cerr << i << " : lower " << p[i]
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm// << " : tmin = " << tmin << ", tmax = " << tmax
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm// << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm // we have to test the closing segment for intersection
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm size_t last = p.size() - 1;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm clower = (p[0][Y] < 0);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (clower != plower) // cross the x axis
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm t = intersect(p[last], p[0], 0);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (tmin > t) tmin = t;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (tmax < t) tmax = t;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm// std::cerr << "0 : lower " << p[0]
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm// << " : tmin = " << tmin << ", tmax = " << tmax << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński if (tmin == 1 && tmax == 0) {
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński return OptInterval();
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński } else {
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński return Interval(tmin, tmax);
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm}
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm/*
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * Clip the Bezier curve "B" wrt the Bezier curve "A" for individuating
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * points which have collinear normals; the new parameter interval
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * for the clipped curve is returned through the output parameter "dom"
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm */
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrmtemplate <>
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof KosińskiOptInterval clip<collinear_normal_tag> (std::vector<Point> const& A,
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński std::vector<Point> const& B,
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński double /*precision*/)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm{
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::vector<Point> F;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm make_focus(F, A);
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński return clip_interval(B, F);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm}
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrmconst double MAX_PRECISION = 1e-8;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrmconst double MIN_CLIPPED_SIZE_THRESHOLD = 0.8;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrmconst Interval UNIT_INTERVAL(0,1);
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosińskiconst OptInterval EMPTY_INTERVAL;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrmconst Interval H1_INTERVAL(0, 0.5);
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosińskiconst Interval H2_INTERVAL(nextafter(0.5, 1.0), 1.0);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm/*
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * iterate
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm *
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * input:
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * A, B: control point sets of two bezier curves
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * domA, domB: real parameter intervals of the two curves
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * precision: required computational precision of the returned parameter ranges
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * output:
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * domsA, domsB: sets of parameter intervals
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm *
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * The parameter intervals are computed by using a Bezier clipping algorithm,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * in case the clipping doesn't shrink the initial interval more than 20%,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * a subdivision step is performed.
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * If during the computation both curves collapse to a single point
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * the routine exits indipendently by the precision reached in the computation
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * of the curve intervals.
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm */
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrmtemplate <>
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrmvoid iterate<intersection_point_tag> (std::vector<Interval>& domsA,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::vector<Interval>& domsB,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::vector<Point> const& A,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::vector<Point> const& B,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm Interval const& domA,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm Interval const& domB,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm double precision )
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm{
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm // in order to limit recursion
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm static size_t counter = 0;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (domA.extent() == 1 && domB.extent() == 1) counter = 0;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (++counter > 100) return;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#if VERBOSE
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::cerr << std::fixed << std::setprecision(16);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::cerr << ">> curve subdision performed <<" << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::cerr << "dom(A) : " << domA << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::cerr << "dom(B) : " << domB << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm// std::cerr << "angle(A) : " << angle(A) << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm// std::cerr << "angle(B) : " << angle(B) << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#endif
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (precision < MAX_PRECISION)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm precision = MAX_PRECISION;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::vector<Point> pA = A;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::vector<Point> pB = B;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::vector<Point>* C1 = &pA;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::vector<Point>* C2 = &pB;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm Interval dompA = domA;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm Interval dompB = domB;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm Interval* dom1 = &dompA;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm Interval* dom2 = &dompB;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński OptInterval dom;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński if ( is_constant(A, precision) && is_constant(B, precision) ){
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm Point M1 = middle_point(C1->front(), C1->back());
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm Point M2 = middle_point(C2->front(), C2->back());
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (are_near(M1,M2)){
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm domsA.push_back(domA);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm domsB.push_back(domB);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm return;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm size_t iter = 0;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm while (++iter < 100
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm && (dompA.extent() >= precision || dompB.extent() >= precision))
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#if VERBOSE
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::cerr << "iter: " << iter << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#endif
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński dom = clip<intersection_point_tag>(*C1, *C2, precision);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
a16a494f042310ee849a6f717ffea70846f1f22cKrzysztof Kosiński if (dom.empty())
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#if VERBOSE
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::cerr << "dom: empty" << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#endif
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm return;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#if VERBOSE
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::cerr << "dom : " << dom << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#endif
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm // all other cases where dom[0] > dom[1] are invalid
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński assert(dom->min() <= dom->max());
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński map_to(*dom2, *dom);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński portion(*C2, *dom);
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński if (is_constant(*C2, precision) && is_constant(*C1, precision))
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm Point M1 = middle_point(C1->front(), C1->back());
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm Point M2 = middle_point(C2->front(), C2->back());
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#if VERBOSE
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::cerr << "both curves are constant: \n"
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm << "M1: " << M1 << "\n"
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm << "M2: " << M2 << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm print(*C2, "C2");
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm print(*C1, "C1");
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#endif
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (are_near(M1,M2))
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm break; // append the new interval
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm else
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm return; // exit without appending any new interval
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm // if we have clipped less than 20% than we need to subdive the curve
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm // with the largest domain into two sub-curves
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński if (dom->extent() > MIN_CLIPPED_SIZE_THRESHOLD)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#if VERBOSE
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński std::cerr << "clipped less than 20% : " << dom->extent() << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::cerr << "angle(pA) : " << angle(pA) << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::cerr << "angle(pB) : " << angle(pB) << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#endif
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::vector<Point> pC1, pC2;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm Interval dompC1, dompC2;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (dompA.extent() > dompB.extent())
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm pC1 = pC2 = pA;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm portion(pC1, H1_INTERVAL);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm portion(pC2, H2_INTERVAL);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm dompC1 = dompC2 = dompA;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm map_to(dompC1, H1_INTERVAL);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm map_to(dompC2, H2_INTERVAL);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm iterate<intersection_point_tag>(domsA, domsB, pC1, pB,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm dompC1, dompB, precision);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm iterate<intersection_point_tag>(domsA, domsB, pC2, pB,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm dompC2, dompB, precision);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm else
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm pC1 = pC2 = pB;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm portion(pC1, H1_INTERVAL);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm portion(pC2, H2_INTERVAL);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm dompC1 = dompC2 = dompB;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm map_to(dompC1, H1_INTERVAL);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm map_to(dompC2, H2_INTERVAL);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm iterate<intersection_point_tag>(domsB, domsA, pC1, pA,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm dompC1, dompA, precision);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm iterate<intersection_point_tag>(domsB, domsA, pC2, pA,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm dompC2, dompA, precision);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm return;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński swap(C1, C2);
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński swap(dom1, dom2);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#if VERBOSE
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::cerr << "dom(pA) : " << dompA << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::cerr << "dom(pB) : " << dompB << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#endif
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm domsA.push_back(dompA);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm domsB.push_back(dompB);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm}
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm/*
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * iterate
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm *
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * input:
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * A, B: control point sets of two bezier curves
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * domA, domB: real parameter intervals of the two curves
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * precision: required computational precision of the returned parameter ranges
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * output:
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * domsA, domsB: sets of parameter intervals
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm *
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * The parameter intervals are computed by using a Bezier clipping algorithm,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * in case the clipping doesn't shrink the initial interval more than 20%,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * a subdivision step is performed.
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * If during the computation one of the two curve interval length becomes less
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * than MAX_PRECISION the routine exits indipendently by the precision reached
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * in the computation of the other curve interval.
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm */
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrmtemplate <>
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrmvoid iterate<collinear_normal_tag> (std::vector<Interval>& domsA,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::vector<Interval>& domsB,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::vector<Point> const& A,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::vector<Point> const& B,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm Interval const& domA,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm Interval const& domB,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm double precision)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm{
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm // in order to limit recursion
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm static size_t counter = 0;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (domA.extent() == 1 && domB.extent() == 1) counter = 0;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (++counter > 100) return;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#if VERBOSE
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::cerr << std::fixed << std::setprecision(16);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::cerr << ">> curve subdision performed <<" << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::cerr << "dom(A) : " << domA << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::cerr << "dom(B) : " << domB << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm// std::cerr << "angle(A) : " << angle(A) << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm// std::cerr << "angle(B) : " << angle(B) << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#endif
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (precision < MAX_PRECISION)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm precision = MAX_PRECISION;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::vector<Point> pA = A;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::vector<Point> pB = B;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::vector<Point>* C1 = &pA;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::vector<Point>* C2 = &pB;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm Interval dompA = domA;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm Interval dompB = domB;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm Interval* dom1 = &dompA;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm Interval* dom2 = &dompB;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński OptInterval dom;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm size_t iter = 0;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm while (++iter < 100
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm && (dompA.extent() >= precision || dompB.extent() >= precision))
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#if VERBOSE
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::cerr << "iter: " << iter << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#endif
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński dom = clip<collinear_normal_tag>(*C1, *C2, precision);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
a16a494f042310ee849a6f717ffea70846f1f22cKrzysztof Kosiński if (dom.empty()) {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#if VERBOSE
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::cerr << "dom: empty" << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#endif
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm return;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#if VERBOSE
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::cerr << "dom : " << dom << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#endif
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński assert(dom->min() <= dom->max());
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński map_to(*dom2, *dom);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm // it's better to stop before losing computational precision
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (iter > 1 && (dom2->extent() <= MAX_PRECISION))
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#if VERBOSE
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::cerr << "beyond max precision limit" << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#endif
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm break;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński portion(*C2, *dom);
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński if (iter > 1 && is_constant(*C2, precision))
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#if VERBOSE
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::cerr << "new curve portion pC1 is constant" << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#endif
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm break;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm // if we have clipped less than 20% than we need to subdive the curve
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm // with the largest domain into two sub-curves
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński if ( dom->extent() > MIN_CLIPPED_SIZE_THRESHOLD)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#if VERBOSE
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński std::cerr << "clipped less than 20% : " << dom->extent() << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::cerr << "angle(pA) : " << angle(pA) << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::cerr << "angle(pB) : " << angle(pB) << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#endif
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::vector<Point> pC1, pC2;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm Interval dompC1, dompC2;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (dompA.extent() > dompB.extent())
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if ((dompA.extent() / 2) < MAX_PRECISION)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm break;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm pC1 = pC2 = pA;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm portion(pC1, H1_INTERVAL);
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński if (false && is_constant(pC1, precision))
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#if VERBOSE
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::cerr << "new curve portion pC1 is constant" << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#endif
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm break;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm portion(pC2, H2_INTERVAL);
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński if (is_constant(pC2, precision))
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#if VERBOSE
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::cerr << "new curve portion pC2 is constant" << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#endif
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm break;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm dompC1 = dompC2 = dompA;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm map_to(dompC1, H1_INTERVAL);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm map_to(dompC2, H2_INTERVAL);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm iterate<collinear_normal_tag>(domsA, domsB, pC1, pB,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm dompC1, dompB, precision);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm iterate<collinear_normal_tag>(domsA, domsB, pC2, pB,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm dompC2, dompB, precision);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm else
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if ((dompB.extent() / 2) < MAX_PRECISION)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm break;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm pC1 = pC2 = pB;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm portion(pC1, H1_INTERVAL);
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński if (is_constant(pC1, precision))
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#if VERBOSE
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::cerr << "new curve portion pC1 is constant" << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#endif
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm break;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm portion(pC2, H2_INTERVAL);
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński if (is_constant(pC2, precision))
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#if VERBOSE
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::cerr << "new curve portion pC2 is constant" << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#endif
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm break;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm dompC1 = dompC2 = dompB;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm map_to(dompC1, H1_INTERVAL);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm map_to(dompC2, H2_INTERVAL);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm iterate<collinear_normal_tag>(domsB, domsA, pC1, pA,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm dompC1, dompA, precision);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm iterate<collinear_normal_tag>(domsB, domsA, pC2, pA,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm dompC2, dompA, precision);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm return;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński swap(C1, C2);
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński swap(dom1, dom2);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#if VERBOSE
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::cerr << "dom(pA) : " << dompA << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::cerr << "dom(pB) : " << dompB << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#endif
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm domsA.push_back(dompA);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm domsB.push_back(dompB);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm}
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm/*
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * get_solutions
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm *
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * input: A, B - set of control points of two Bezier curve
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * input: precision - required precision of computation
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * input: clip - the routine used for clipping
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * output: xs - set of pairs of parameter values
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * at which the clipping algorithm converges
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm *
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * This routine is based on the Bezier Clipping Algorithm,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * see: Sederberg - Computer Aided Geometric Design
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm */
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrmtemplate <typename Tag>
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrmvoid get_solutions (std::vector< std::pair<double, double> >& xs,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::vector<Point> const& A,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::vector<Point> const& B,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm double precision)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm{
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::pair<double, double> ci;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::vector<Interval> domsA, domsB;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm iterate<Tag> (domsA, domsB, A, B, UNIT_INTERVAL, UNIT_INTERVAL, precision);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm if (domsA.size() != domsB.size())
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm assert (domsA.size() == domsB.size());
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm xs.clear();
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm xs.reserve(domsA.size());
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm for (size_t i = 0; i < domsA.size(); ++i)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm {
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#if VERBOSE
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::cerr << i << " : domA : " << domsA[i] << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::cerr << "extent A: " << domsA[i].extent() << " ";
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::cerr << "precision A: " << get_precision(domsA[i]) << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::cerr << i << " : domB : " << domsB[i] << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::cerr << "extent B: " << domsB[i].extent() << " ";
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::cerr << "precision B: " << get_precision(domsB[i]) << std::endl;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm#endif
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm ci.first = domsA[i].middle();
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm ci.second = domsB[i].middle();
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm xs.push_back(ci);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm }
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm}
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm} /* end namespace bezier_clipping */ } /* end namespace detail */
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm/*
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * find_collinear_normal
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm *
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * input: A, B - set of control points of two Bezier curve
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * input: precision - required precision of computation
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * output: xs - set of pairs of parameter values
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * at which there are collinear normals
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm *
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * This routine is based on the Bezier Clipping Algorithm,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * see: Sederberg, Nishita, 1990 - Curve intersection using Bezier clipping
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm */
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrmvoid find_collinear_normal (std::vector< std::pair<double, double> >& xs,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::vector<Point> const& A,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::vector<Point> const& B,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm double precision)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm{
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm using detail::bezier_clipping::get_solutions;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm using detail::bezier_clipping::collinear_normal_tag;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm get_solutions<collinear_normal_tag>(xs, A, B, precision);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm}
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm/*
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * find_intersections_bezier_clipping
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm *
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * input: A, B - set of control points of two Bezier curve
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * input: precision - required precision of computation
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * output: xs - set of pairs of parameter values
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * at which crossing happens
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm *
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * This routine is based on the Bezier Clipping Algorithm,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm * see: Sederberg, Nishita, 1990 - Curve intersection using Bezier clipping
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm */
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrmvoid find_intersections_bezier_clipping (std::vector< std::pair<double, double> >& xs,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::vector<Point> const& A,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm std::vector<Point> const& B,
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm double precision)
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm{
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm using detail::bezier_clipping::get_solutions;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm using detail::bezier_clipping::intersection_point_tag;
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm get_solutions<intersection_point_tag>(xs, A, B, precision);
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm}
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm} // end namespace Geom
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm/*
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm Local Variables:
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm mode:c++
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm c-file-style:"stroustrup"
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm indent-tabs-mode:nil
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm fill-column:99
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm End:
da58597f9f9ecb17c4f545c4483a844a363bcc27pjrm*/
a4030d5ca449e7e384bc699cd249ee704faaeab0Chris Morgan// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :