CGAL users discussion list

[Marko Kukovec <>]

Hello!

About a week ago I put a question on this mailing list about
offsetting/insetting 2D polygons in CGAL. I have a problem as I have a
polygons with arcs and lines so I can't directly feed the
offsetting/insetting algorithm. I got an answer from Fernando Cacciola
and he suggested I should polygonize all arcs prior to feeding the
offsetting algorithm. He even pointed out this algorithm
http://www.worldserver.com/turk/computergraphics/CircularArcSubdivision.pdf .


The algorithm code looks like (from Graphic Gems book code)
===========================================================

/* returns squared length of input vector */
double V2SquaredLength(a)
Vector2 *a;
{ return((a->x * a->x)+(a->y * a->y));
}

/* returns length of input vector */
double V2Length(a)
Vector2 *a;
{
return(sqrt(V2SquaredLength(a)));
}


/* scales the input vector to the new length and returns it */
Vector2 *V2Scale(v, newlen)
Vector2 *v;
double newlen;
{
double len = V2Length(v);
if (len != 0.0) { v->x *= newlen/len; v->y *= newlen/len; }
return(v);
}


/* Function prototype for externally defined functions */
void DrawLine(Point2 p0, Point2 p1);

void
DrawArc(Point2 p0, Point2 p1, double d)
{
    if (fabs(d) <= DMAX) DrawLine(p0, p1);
    else {
        Vector2 v;
        Point2  pm, pb;
        double  dSub;
        
        v.x = p1.x - p0.x;       /* vector from p0 to p1 */
        v.y = p1.y - p0.y;
        
        pm.x = p0.x + 0.5 * v.x; /* midpoint */
        pm.y = p0.y + 0.5 * v.y;
        
        dSub = d / 4;
        V2Scale(&v, dSub);       /* subdivided vector */
        
        pb.x = pm.x - v.y;       /* bisection point */
        pb.y = pm.y + v.x;
        
        DrawArc(p0, pb, dSub);   /* first half arc */
        DrawArc(pb, p1, dSub);   /* second half arc */
    }
}


My implementation following the code from the book (second method,
CPolygonizer::arcSubdivision)
=======================================================================

void 
CPolygonizer::doPolygonize(const Polygon2D& aPolygon, Polygon2D&
aPolygonized) {
        Polygon2D::Curve_const_iterator it = aPolygon.curves_begin();
        for(; it != aPolygon.curves_end(); it++){
                        static int curveCount = 1;
                if(it->is_circular()){
                        CPoint start(CGAL::to_double(it->source().x()),
CGAL::to_double(it->source().y()));
                        CPoint end(CGAL::to_double(it->target().x()),
CGAL::to_double(it->target().y()));
                        CPoint 
center(CGAL::to_double(it->supporting_circle().center().x()),
CGAL::to_double(it->supporting_circle().center().y()));

                        std::cout << curveCount << " ARC start(" << 
start.getX() << ", " <<
start.getY() << ")" << " end(" << end.getX() << ", " << end.getY() <<
")" << std::endl;
                        double radius = 
sqrt(pow(CGAL::to_double(it->source().x()) -
CGAL::to_double(it->supporting_circle().center().x()), 2) +
pow(CGAL::to_double(it->source().y()) -
CGAL::to_double(it->supporting_circle().center().y()), 2));
                        std::cout << curveCount << " ARC radius: " << radius 
<< std::endl;
                        double length = 
sqrt(pow(CGAL::to_double(it->source().x()) -
CGAL::to_double(it->target().x()), 2) +
pow(CGAL::to_double(it->source().y()) -
CGAL::to_double(it->target().y()), 2));
                        std::cout << curveCount << " ARC length: " << length 
<< std::endl;
                        double sagitta = radius - sqrt(pow(radius, 2) - 
pow(length / 2, 2));
                        std::cout << curveCount << " ARC sagitta: " << 
sagitta << std::endl
<< std::endl; 

                        arcSubdivision(start, end, sagitta, aPolygonized);

                        std::cout << std::endl << std::endl;
                } else {
                        CPoint start(CGAL::to_double(it->source().x()),
CGAL::to_double(it->source().y()));
                        CPoint end(CGAL::to_double(it->target().x()),
CGAL::to_double(it->target().y()));
                        std::cout << curveCount << " LINE start(" << 
start.getX() << ", " <<
start.getY() << ")" << " end(" << end.getX() << ", " << end.getY() <<
")" << std::endl << std::endl;

                        aPolygonized.push_back(XMonotoneCurve2D(start, end));
                }
                curveCount++; 
        }
}

void 
CPolygonizer::arcSubdivision(const CPoint& aStart, const CPoint& aEnd,
const double& aSagitta, Polygon2D& aPolygonized) { 
        if(fabs(aSagitta) <= mMaxChordalDeviation) {
                aPolygonized.push_back(XMonotoneCurve2D(aStart, aEnd));
                std::cout << "start(" << aStart.getX() << ", " << 
aStart.getY() << ")"
<< " end(" << aEnd.getX() << ", " << aEnd.getY() << ")" << std::endl; 
        } else { 
                double vectorX = aEnd.getX() - aStart.getX();
                double vectorY = aEnd.getY() - aStart.getY();
                CPoint midpoint(aStart.getX() + 0.5 * vectorX, aStart.getY() 
+ 0.5 *
vectorY);
                double chordalDeviation = aSagitta / 4;
                double vectorLen = sqrt(pow(vectorX, 2) + pow(vectorY, 2));
        if(0.0 != vectorLen){
                vectorX *= chordalDeviation / vectorLen;
                vectorY *= chordalDeviation / vectorLen;
        }
        CPoint bisectionPoint(midpoint.getX() - vectorY, midpoint.getY() +
vectorX);
        arcSubdivision(aStart, bisectionPoint, chordalDeviation, 
aPolygonized);
        arcSubdivision(bisectionPoint, aEnd, chordalDeviation, aPolygonized);
}
}

===============================================================

Now, if I execute the original algorithm (my implementation) I can see
that first calculated Pb point is wrong. X point value is 200 as it
should be, but Y point is wrong. It is only the quoter the value as it
should be (50 but should be 200).

I analyzed the algorithm and noted the line

V2Scale(&v, dSub);       /* subdivided vector */

and changed it to 

V2Scale(&v, d);       /* subdivided vector */

or in my implementation the lines

vectorX *= chordalDeviation / vectorLen;
vectorY *= chordalDeviation / vectorLen;

to

vectorX *= aSagitta / vectorLen;
vectorY *= aSagitta / vectorLen;

Now it looks like calculation of Pb points are correct. Is it possible
there is a typo in the Graphic Gems book or I did something wrong? I
looked into book errata and I couldn't find any note for this chapter.
My tests can be seen on this picture
http://img232.imageshack.us/img232/6633/arcsubdivisionpf5.jpg . Black
lines are the original polygon, the red ones are polygonized version of
the original polygon.
>From the picture you can see that in the fixed algorithm first Pb point
is correct, but additional ones are not as accurate as they should be
(do not touch the arc curve). I know that the algorithm document
(Turkowski) explains that it works well for the arcs that are subtending
for less than 75˚, but is this really the reason for the error here?
Should I divide the arc in case of the angle bigger than 75˚ into
smaller arc (with trigonometry) prior to polygonization?

All the polygons you can see on the picture are drawn with the Cairomm
library using context arc method.

if(it->is_circular()){
Circle2D circle = it->supporting_circle();
Point2D center = circle.center();
double centerX = CGAL::to_double(center.x());
double centerY = CGAL::to_double(center.y());
double radious = sqrt(CGAL_NTS to_double(circle.squared_radius()));
int orientation = circle.orientation();
double angle1 = getAngle(startX - centerX, startY - centerY);
double angle2 = getAngle(endX - centerX, endY - centerY);
mContext->arc(centerX, centerY, radious, angle1, angle2);
}

I wonder why the rendered arcs by Cairo looks perfect? What algorithm do
they use? It seems it is not affected by the arc angle.

Best regards,
Marko Kukovec