Root/cameo/area.c

1/*
2 * area.c - Area fill
3 *
4 * Written 2012 by Werner Almesberger
5 * Copyright 2012 Werner Almesberger
6 *
7 * This program is free software; you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation; either version 2 of the License, or
10 * (at your option) any later version.
11 */
12
13/*
14 * We use the following requirement to simplify toolpath generation: the
15 * outlines must be designed such that the tool can pass along all the
16 * outlines without cutting into anything it's not supposed to.
17 */
18
19#include <stddef.h>
20#include <math.h>
21#include <assert.h>
22
23#include "util.h"
24#include "path.h"
25#include "area.h"
26
27
28#define EPSILON 1e-6
29
30
31static int bbox(const struct path *path,
32    double *xa, double *ya, double *xb, double *yb)
33{
34    const struct point *p = path->first;
35
36    if (!p)
37        return 0;
38    *xa = *xb = p->x;
39    *ya = *yb = p->y;
40    while (p) {
41        if (p->x < *xa)
42            *xa = p->x;
43        if (p->x > *xb)
44            *xb = p->x;
45        if (p->y < *ya)
46            *ya = p->y;
47        if (p->y > *yb)
48            *yb = p->y;
49        p = p->next;
50    }
51    return 1;
52}
53
54
55/*
56 * @@@ this is a bit too simple. E.g., it would report A as being inside B
57 * in this case:
58 *
59 * +---+
60 * +---+ | |
61 * | A | | |
62 * +---+ | |
63 * | B |
64 * +--------+ |
65 * | |
66 * +------------+
67 */
68
69static int is_inside(const struct path *a, const struct path *b)
70{
71        double xa, ya, xb, yb;
72    const struct point *p;
73
74    if (!bbox(b, &xa, &ya, &xb, &yb))
75        return 0;
76    for (p = a->first; p; p = p->next)
77        if (p->x < xa || p->x > xb ||
78            p->y < ya || p->y > yb)
79            return 0;
80    return 1;
81}
82
83
84/*
85 * Solve
86 *
87 * ax+by = e
88 * cx+dy = f
89 *
90 * with Cramer's rule:
91 * http://en.wikipedia.org/wiki/Cramer's_rule
92 */
93
94static int cramer2(double a, double b, double c, double d, double e, double f,
95    double *x, double *y)
96{
97    double det;
98
99    det = a*d-b*c;
100    if (fabs(det) < EPSILON)
101        return 0;
102    *x = (e*d-b*f)/det;
103    *y = (a*f-e*c)/det;
104    return 1;
105}
106
107
108/*
109 * Solve
110 *
111 * ax + na*bx = cx + nb*dx
112 * ay + na*by = cy + nb*dy
113 *
114 * which is
115 *
116 * na*bx + nb*-dx = cx - ax
117 * na*by + nb*-dy = cy - ay
118 */
119
120static int intersect(double ax, double ay, double bx, double by,
121    double cx, double cy, double dx, double dy, double *na, double *nb)
122{
123    return cramer2(bx, -dx, by, -dy, cx-ax, cy-ay, na, nb);
124}
125
126
127/*
128 * See above.fig. The equation we solve is
129 *
130 * Q = A+u*(AM)
131 * Q = B+v*(BP)
132 *
133 * equals
134 *
135 * ax + u*(mx-ax) = bx + v*(px-bx)
136 * ay + u*(my-ay) = by + v*(py-by)
137 *
138 * equals
139 *
140 * u*(mx-ax) + v*(bx-px) = bx - ax
141 * u*(my-ay) + v*(by-py) = by - ay
142 *
143 * For BC, the equation becomes
144 *
145 * Q = C+u*(CM)
146 * Q = B+v*(BP)
147 */
148
149static int above(const struct point *a, const struct point *b,
150    const struct point *c, double px, double py)
151{
152    double ab, bc;
153    double mx, my;
154    double u, v;
155
156    ab = hypot(a->x-b->x, a->y-b->y);
157    bc = hypot(b->x-c->x, b->y-c->y);
158    if (fabs(ab) < EPSILON || fabs(bc) < EPSILON)
159        return 0;
160
161    mx = b->x-(b->y-a->y)/ab-(c->y-b->y)/bc;
162    my = b->y+(b->x-a->x)/ab+(c->x-b->x)/bc;
163
164    if (cramer2(mx-a->x, b->x-px, my-a->y, b->y-py, b->x-a->x, b->y-a->y,
165        &u, &v))
166        if (u >= 0 && u <= 1 && v >= 0)
167            return 1;
168    if (cramer2(mx-c->x, b->x-px, my-c->y, b->y-py, b->x-c->x, b->y-c->y,
169        &u, &v))
170        if (u >= 0 && u <= 1 && v >= 0)
171            return 1;
172    return 0;
173}
174
175
176/*
177 * Solve
178 *
179 * (ax+n*bx-cx)^2+(ay+n*by-cy)^2 = r^2 for n
180 *
181 * http://en.wikipedia.org/wiki/Quadratic_equation
182 */
183
184static int touch_solve(double ax, double ay, double bx, double by,
185    double cx, double cy, double r, int enter, double *n)
186{
187    double dx = cx-ax;
188    double dy = cy-ay;
189    double a = bx*bx+by*by; /* always positive */
190    double b = -2*bx*dx-2*by*dy;
191    double c = dx*dx+dy*dy-r*r;
192    double d, tmp;
193
194    d = b*b-4*a*c;
195    if (d < 0)
196        return 0;
197    d = sqrt(d);
198    tmp = enter ? (-b-d)/2/a : (-b+d)/2/a;
199    if (tmp <= 0 || tmp >= 1)
200        return 0;
201    *n = tmp;
202    return 1;
203}
204
205
206/*
207 * The points A, B, and C are (if the path is left-handed):
208 *
209 * - A: the beginning of the segment leading into the corner
210 * - B: the corner point
211 * - C: the beginning of the segment leading out of the corner
212 *
213 * If the path is right-handed, we swap A and C, making it left-handed.
214 */
215
216static int touch(double ax, double ay, double bx, double by,
217     const struct point *a, const struct point *b, const struct point *c,
218     double r, int enter, int left, double *n)
219{
220    double px, py;
221
222    if (!touch_solve(ax, ay, bx, by, b->x, b->y, r, enter, n))
223        return 0;
224    px = ax+*n*bx;
225    py = ay+*n*by;
226    return above(a, b, c, px, py) == left;
227}
228
229
230/*
231 * Here, the points A, B, C, and D are:
232 *
233 * - A: before the beginning of the current segment
234 * - B: the beginning
235 * - C: the end
236 * - D: the next point beyond the end
237 */
238
239static int hit_segment(double fx, double fy, double tx, double ty,
240    const struct point *a, const struct point *b, const struct point *c,
241    const struct point *d, double r, int enter, int left, double *n)
242{
243    double dx, dy, nx, ny, nn;
244    double px, py;
245    double na, nb;
246
247    tx -= fx;
248    ty -= fy;
249
250    dx = c->x-b->x;
251    dy = c->y-b->y;
252
253    if (left) {
254        nx = dx;
255        ny = dy;
256    } else {
257        nx = -dx;
258        ny = -dy;
259    }
260
261    /* -dy becomes the x component of the normal vector */
262    if (enter ? ny < 0 : ny > 0)
263        return 0;
264
265    nn = hypot(nx, ny);
266
267    px = b->x-ny/nn*r;
268    py = b->y+nx/nn*r;
269
270    if (!intersect(fx, fy, tx, ty, px, py, dx, dy, &na, &nb))
271        return 0;
272    if (nb <= 0) {
273        if (!touch(fx, fy, tx, ty, a, b, c, r, enter, left, &na))
274            return 0;
275    }
276    if (nb >= 1) {
277        if (!touch(fx, fy, tx, ty, b, c, d, r, enter, left, &na))
278            return 0;
279    }
280    if (na <= 0 || na >= 1)
281        return 0;
282    *n = na;
283    return 1;
284}
285
286
287static int hit_path(double fx, double fy, double tx, double ty,
288    const struct path *path, int inside, int enter, double r, double *x)
289{
290    const struct point *p, *last, *next2;
291    int left;
292    double nx, tmp;
293    int found = 0;
294
295    /*
296     * @@@ We don't wrap around the ends properly and create a zero-sized
297     * imaginary segment between path->first and path->last.
298     */
299    left = path_tool_is_left(path);
300    if (inside)
301        left = !left;
302    last = path->last;
303    for (p = path->first; p != path->last; p = p->next) {
304        next2 = p->next->next ? p->next->next : path->first;
305        if (hit_segment(fx, fy, tx, ty, last, p, p->next, next2,
306            r, enter, left, &tmp)) {
307            if (!found || nx > tmp)
308                nx = tmp;
309            found = 1;
310        }
311        last = p;
312    }
313    if (found)
314        *x = fx+nx*(tx-fx);
315    return found;
316}
317
318
319static const struct path **subordinates(const struct path *paths,
320    const struct path *path, double z)
321{
322    const struct path **sub, **w, **a, **b;;
323    const struct path *p;
324    int n = 0;
325
326    for (p = paths; p; p = p->next)
327        if (p->first && p->first->z == z)
328            n++;
329    sub = alloc_size(sizeof(struct path *)*n);
330    w = sub;
331    for (p = paths; p; p = p->next)
332        if (p != path && p->first && p->first->z == z &&
333            is_inside(p, path) && !is_inside(path, p))
334            *w++ = p;
335    *w = NULL;
336    for (a = sub; a != w; a++)
337        for (b = sub; b != w; b++)
338            if (a != b && is_inside(*a, *b)) {
339                *a = *--w;
340                *w = NULL;
341                a--;
342                break;
343            }
344    return sub;
345}
346
347
348static void do_line(const struct path *path, const struct path **sub,
349    double xa, double xb, double y, double r_tool, double overlap,
350    struct path **res)
351{
352    const struct path *last = path;
353    const struct path **s;
354    struct path *new;
355    double x, next;
356
357    if (!hit_path(xa-3*r_tool, y, xb, y, last, 1, 0, r_tool, &x))
358        return;
359    while (1) {
360        next = xb;
361        last = NULL;
362        if (hit_path(x, y, xb, y, path, 1, 1, r_tool, &next))
363            last = path;
364        for (s = sub; *s; s++)
365            if (hit_path(x, y, next, y, *s, 0, 1, r_tool, &next))
366                last = *s;
367        if (next-x > 2*r_tool-2*overlap) {
368            new = path_new(r_tool, "");
369            path_add(new, x+r_tool-overlap, y, path->first->z);
370            path_add(new, next-r_tool+overlap, y, path->first->z);
371            new->next = *res;
372            *res = new;
373        }
374        if (!last)
375            return;
376        if (!hit_path(next+EPSILON, y, xb, y, last, last == path, 0,
377            r_tool, &x))
378            return;
379    }
380}
381
382
383static void add_outline(const struct path *path, int inside, struct path **res)
384{
385    struct path *new;
386    int left;
387
388    left = path_tool_is_left(path);
389    new = path_offset(path, inside ? !left : left, 0);
390    new->next = *res;
391    *res = new;
392}
393
394
395static void fill_path(const struct path *paths, const struct path *path,
396    double z, double r_tool, double overlap, struct path **res)
397{
398    const struct path **sub, **s;
399    const struct path **sub2, **s2;
400        double xa, ya, xb, yb;
401    int n, i;
402
403    if (!bbox(path, &xa, &ya, &xb, &yb))
404        return;
405    sub = subordinates(paths, path, z);
406    xa += r_tool;
407    ya += 3*r_tool-overlap;
408    xb -= r_tool;
409    yb -= 3*r_tool-overlap;
410    n = ceil((yb-ya)/(2*r_tool-overlap));
411    for (i = 0; i <= n; i++)
412        do_line(path, sub, xa, xb, ya+(yb-ya)*((double) i/n),
413            r_tool, overlap, res);
414    for (s = sub; *s; s++) {
415        sub2 = subordinates(paths, *s, z);
416        for (s2 = sub2; *s2; s2++)
417            fill_path(paths, *s2, z, r_tool, overlap, res);
418        free(sub2);
419        add_outline(*s, 0, res);
420    }
421    free(sub);
422    add_outline(path, 1, res);
423}
424
425
426struct path *area(const struct path *paths, double overlap)
427{
428    struct path *res = NULL;
429    double z = HUGE_VAL, best_x, x;
430    const struct path *path, *best;
431    const struct point *p;
432
433    if (!paths)
434        return NULL;
435    while (1) {
436        best = NULL;
437        best_x = HUGE_VAL;
438        for (path = paths; path; path = path->next) {
439            if (!path->first)
440                continue;
441            if (path->first->z >= z)
442                continue;
443            x = HUGE_VAL;
444            for (p = path->first; p; p = p->next)
445                if (p->x < x)
446                    x = p->x;
447            if (best && best->first->z >= path->first->z &&
448                x >= best_x)
449                continue;
450            best = path;
451            best_x = x;
452        }
453        if (!best)
454            return res;
455        z = best->first->z;
456        fill_path(paths, best, z, best->r_tool, overlap, &res);
457    }
458}
459

Archive Download this file

Branches:
master



interactive