Total circles area: Difference between revisions
Content added Content deleted
m (→C: Scanline Method: make it not compile, so users trying it won't be puzzled by incorrect-looking output) |
(→C: Scanline Method: removed separate sorting pass to have fewer things to look at; some comments) |
||
Line 167: | Line 167: | ||
This version performs about 5 million scanlines in about a second, result should be accurate to maybe 10 decimal points. |
This version performs about 5 million scanlines in about a second, result should be accurate to maybe 10 decimal points. |
||
<lang c>#include <stdio.h> |
<lang c>#include <stdio.h> |
||
#include < |
#include <string.h> |
||
#include <stdlib.h> |
#include <stdlib.h> |
||
#include < |
#include <math.h> |
||
typedef double |
typedef double flt; |
||
typedef struct { |
typedef struct { |
||
flt x, y, r, r2; |
|||
flt y0, y1; // extent of circle y+r and y-r |
|||
flt x0, x1; // where scanline intersects circle |
|||
} circle_t; |
|||
#define SZ sizeof(circle_t) |
|||
circle_t circles[] = { |
|||
{ 1.6417233788, 1.6121789534, 0.0848270516}, |
|||
#error data snipped |
#error data snipped for space; copy from previous C example |
||
}; |
|||
{ 0.0152957411, 0.0638919221, 0.9771215985}}; |
|||
flt max(flt x, flt y) { return x < y ? y : x; } |
|||
flt min(flt x, flt y) { return x > y ? y : x; } |
|||
flt sq(flt x) { return x * x; } |
|||
flt cdist(circle_t *c1, circle_t *c2) { |
|||
return sqrt(sq(c1->x - c2->x) + sq(c1->y - c2->y)); |
|||
inline Fp dist(const Fp x0, const Fp y0, const Fp x1, const Fp y1) { |
|||
return sqrt(square(x0 - x1) + square(y0 - y1)); |
|||
} |
} |
||
inline void swap_c(circle_t *c) |
|||
typedef struct { Fp x0, x1; } Sect; |
|||
{ |
|||
circle_t tmp = c[0]; |
|||
c[0] = c[1], c[1] = tmp; |
|||
} |
|||
flt area(circle_t *circs, int n_circ, flt ymin, flt ymax, flt step) |
|||
Fp spans(Circle *c, size_t n_circles, |
|||
{ |
|||
const Fp ymin, const Fp step, const Fp ymax) { |
|||
int i, n = n_circ; |
|||
// Safe if n_circles is not very large. |
|||
Sect sects[n_circles]; |
|||
circle_t *c = malloc(SZ * n); |
|||
Fp total = 0; |
|||
memcpy(c, circs, SZ * n); |
|||
while (n--) |
|||
// GNU C99 extension. |
|||
for (i = 0; i < n; i++) |
|||
inline void swap_s(const size_t i, const size_t j) { |
|||
if (c[i].y1 < c[i+1].y1) swap_c(c + i); |
|||
const Sect ss = sects[i]; |
|||
sects[i] = sects[j]; |
|||
sects[j] = ss; |
|||
} |
|||
flt total = 0; |
|||
inline void swap_c(const size_t i, const size_t j) { |
|||
const Circle cc = c[i]; |
|||
c[i] = c[j]; |
|||
c[j] = cc; |
|||
} |
|||
int row = 1 + ceil((ymax - ymin) / step); |
|||
for (size_t i = 0; i < n_circles; i++) |
|||
while (row--) { |
|||
for (size_t j = i + 1; j < n_circles; j++) |
|||
flt y = ymin + step * row; |
|||
if (c[i].y1 < c[j].y1) |
|||
for (n = 0; n < n_circ; n++) |
|||
swap_c(i, j); |
|||
if (y >= c[n].y1) // rest of circles below scanline, ignore |
|||
break; |
|||
else if (y > c[n].y0) { |
|||
flt dx = sqrt(c[n].r2 - sq(y - c[n].y)); |
|||
c[n].x0 = c[n].x - dx; |
|||
c[n].x1 = c[n].x + dx; |
|||
// keep circles sorted by left intersection |
|||
for (int row = 1 + ceil((ymax - ymin) / step); row; row--) { |
|||
for (i = n; i-- && c[i].x0 > c[i+1].x0; swap_c(c + i)); |
|||
// Could probably do "y -= step", not sure. |
|||
const Fp y = ymin + step * row; |
|||
} else {// remove a circle when scanline has passed it |
|||
size_t n = n_circles; |
|||
memmove(c + n, c + n + 1, SZ * (--n_circ - n)); |
|||
n--; |
|||
if (y >= c[i].y1) |
|||
} |
|||
n = i; |
|||
else if (y <= c[i].y0) { |
|||
const Circle cc = c[i]; |
|||
for (size_t j = i + 1; j < n_circles; j++) |
|||
c[j - 1] = c[j]; |
|||
n_circles--; |
|||
c[n_circles] = cc; |
|||
n--; |
|||
} else { |
|||
const Fp dx = sqrt(c[i].r2 - square(y - c[i].y)); |
|||
sects[i].x0 = c[i].x - dx; |
|||
sects[i].x1 = c[i].x + dx; |
|||
i++; |
|||
} |
|||
} |
|||
if (!n) continue; |
|||
continue; |
|||
flt right = c->x1; |
|||
for (size_t i = n; i--; ) { // Bubble sort. Important. |
|||
total += c->x1 - c->x0; |
|||
bool swapped = false; |
|||
for (size_t j = 0; j < i; j++) |
|||
if (sects[j].x0 > sects[j + 1].x0) { |
|||
swap_s(j, j + 1); |
|||
swap_c(j, j + 1); |
|||
swapped = true; |
|||
} |
|||
if (!swapped) |
|||
break; |
|||
} |
|||
for (i = 1; i < n; i++) { |
|||
Sect *cur = sects; |
|||
if (c[i].x1 <= right) continue; |
|||
total += cur->x1 - cur->x0; |
|||
total += c[i].x1 - max(c[i].x0, right); |
|||
for (size_t i = 1; i < n; i++) { |
|||
right = c[i].x1; |
|||
if (sects[i].x1 <= cur->x1) |
|||
} |
|||
continue; |
|||
} |
|||
total += sects[i].x1 - max(sects[i].x0, cur->x1); |
|||
cur = sects + i; |
|||
} |
|||
} |
|||
free(c); |
|||
return total; |
|||
return total * step; |
|||
} |
} |
||
int main() |
int main(void) |
||
{ |
|||
size_t n_circles = sizeof(circles) / sizeof(Circle); |
|||
int n_circ = sizeof(circles) / SZ; |
|||
flt ymin = INFINITY, ymax = -INFINITY; |
|||
circle_t *c1, *c2; |
|||
// Remove fully covered circles. |
|||
for (c1 = circles + n_circ; c1-- > circles; ) { |
|||
for (c2 = circles + n_circ; c2-- > circles; ) |
|||
// throw out circles inside another circle |
|||
if (c1 != c2 && cdist(c1, c2) + c1->r <= c2->r) { |
|||
*c1 = circles[--n_circ]; |
|||
dist(c1->x, c1->y, c2->x, c2->y) + c1->r <= c2->r) { |
|||
break; |
|||
n_circles--; |
|||
} |
|||
*c1 = circles[n_circles]; |
|||
ymin = min(ymin, c1->y0 = c1->y - c1->r); |
|||
break; |
|||
ymax = max(ymax, c1->y1 = c1->y + c1->r); |
|||
} |
|||
c1->r2 = sq(c1->r); |
|||
} |
|||
flt s = 1. / (1 << 20); |
|||
Fp miny = INFINITY, maxy = -INFINITY; |
|||
flt y0 = floor(ymin / s) * s; |
|||
for (size_t i = 0; i < n_circles; i++) { |
|||
flt a = area(circles, n_circ, y0, ymax, s); |
|||
circles[i].r2 = square(circles[i].r); |
|||
int nlines = (ymax - y0) / s; |
|||
circles[i].y0 = circles[i].y - circles[i].r; |
|||
circles[i].y1 = circles[i].y + circles[i].r; |
|||
miny = min(miny, circles[i].y0); |
|||
maxy = max(maxy, circles[i].y1); |
|||
} |
|||
// roughly, it cease to make sense if sqrt(nlines) * 1e-14 >> s * s |
|||
const size_t s = 1U << 20; // Definitely overkill. Definitely. |
|||
printf("area = %.10f\tat %d scanlines\n", a, nlines); |
|||
const Fp sp = spans(circles, n_circles, |
|||
floor(miny * s) / s, 1.0 / s, maxy); |
|||
printf("Area = %.10f at\t%d scanlines.\n", |
|||
sp / s, (int)((maxy - floor(miny * s) / s) * s)); |
|||
return 0; |
|||
}</lang> |
}</lang> |
||
{{out}} |
{{out}} |
||
area = 21.5650366037 at 5637290 scanlines |
|||
=={{header|D}}== |
=={{header|D}}== |