#
source:
trunk/third/libart_lgpl/art_svp_wind.c
@
20813

Revision 20813, 39.6 KB checked in by ghudson, 20 years ago (diff) |
---|

Line | |
---|---|

1 | /* Libart_LGPL - library of basic graphic primitives |

2 | * Copyright (C) 1998-2000 Raph Levien |

3 | * |

4 | * This library is free software; you can redistribute it and/or |

5 | * modify it under the terms of the GNU Library General Public |

6 | * License as published by the Free Software Foundation; either |

7 | * version 2 of the License, or (at your option) any later version. |

8 | * |

9 | * This library is distributed in the hope that it will be useful, |

10 | * but WITHOUT ANY WARRANTY; without even the implied warranty of |

11 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |

12 | * Library General Public License for more details. |

13 | * |

14 | * You should have received a copy of the GNU Library General Public |

15 | * License along with this library; if not, write to the |

16 | * Free Software Foundation, Inc., 59 Temple Place - Suite 330, |

17 | * Boston, MA 02111-1307, USA. |

18 | */ |

19 | |

20 | /* Primitive intersection and winding number operations on sorted |

21 | vector paths. |

22 | |

23 | These routines are internal to libart, used to construct operations |

24 | like intersection, union, and difference. */ |

25 | |

26 | #include "config.h" |

27 | #include "art_svp_wind.h" |

28 | |

29 | #include <stdio.h> /* for printf of debugging info */ |

30 | #include <string.h> /* for memcpy */ |

31 | #include <math.h> |

32 | #include "art_misc.h" |

33 | |

34 | #include "art_rect.h" |

35 | #include "art_svp.h" |

36 | |

37 | #define noVERBOSE |

38 | |

39 | #define PT_EQ(p1,p2) ((p1).x == (p2).x && (p1).y == (p2).y) |

40 | |

41 | #define PT_CLOSE(p1,p2) (fabs ((p1).x - (p2).x) < 1e-6 && fabs ((p1).y - (p2).y) < 1e-6) |

42 | |

43 | /* return nonzero and set *p to the intersection point if the lines |

44 | z0-z1 and z2-z3 intersect each other. */ |

45 | static int |

46 | intersect_lines (ArtPoint z0, ArtPoint z1, ArtPoint z2, ArtPoint z3, |

47 | ArtPoint *p) |

48 | { |

49 | double a01, b01, c01; |

50 | double a23, b23, c23; |

51 | double d0, d1, d2, d3; |

52 | double det; |

53 | |

54 | /* if the vectors share an endpoint, they don't intersect */ |

55 | if (PT_EQ (z0, z2) || PT_EQ (z0, z3) || PT_EQ (z1, z2) || PT_EQ (z1, z3)) |

56 | return 0; |

57 | |

58 | #if 0 |

59 | if (PT_CLOSE (z0, z2) || PT_CLOSE (z0, z3) || PT_CLOSE (z1, z2) || PT_CLOSE (z1, z3)) |

60 | return 0; |

61 | #endif |

62 | |

63 | /* find line equations ax + by + c = 0 */ |

64 | a01 = z0.y - z1.y; |

65 | b01 = z1.x - z0.x; |

66 | c01 = -(z0.x * a01 + z0.y * b01); |

67 | /* = -((z0.y - z1.y) * z0.x + (z1.x - z0.x) * z0.y) |

68 | = (z1.x * z0.y - z1.y * z0.x) */ |

69 | |

70 | d2 = a01 * z2.x + b01 * z2.y + c01; |

71 | d3 = a01 * z3.x + b01 * z3.y + c01; |

72 | if ((d2 > 0) == (d3 > 0)) |

73 | return 0; |

74 | |

75 | a23 = z2.y - z3.y; |

76 | b23 = z3.x - z2.x; |

77 | c23 = -(z2.x * a23 + z2.y * b23); |

78 | |

79 | d0 = a23 * z0.x + b23 * z0.y + c23; |

80 | d1 = a23 * z1.x + b23 * z1.y + c23; |

81 | if ((d0 > 0) == (d1 > 0)) |

82 | return 0; |

83 | |

84 | /* now we definitely know that the lines intersect */ |

85 | /* solve the two linear equations ax + by + c = 0 */ |

86 | det = 1.0 / (a01 * b23 - a23 * b01); |

87 | p->x = det * (c23 * b01 - c01 * b23); |

88 | p->y = det * (c01 * a23 - c23 * a01); |

89 | |

90 | return 1; |

91 | } |

92 | |

93 | #define EPSILON 1e-6 |

94 | |

95 | static double |

96 | trap_epsilon (double v) |

97 | { |

98 | const double epsilon = EPSILON; |

99 | |

100 | if (v < epsilon && v > -epsilon) return 0; |

101 | else return v; |

102 | } |

103 | |

104 | /* Determine the order of line segments z0-z1 and z2-z3. |

105 | Return +1 if z2-z3 lies entirely to the right of z0-z1, |

106 | -1 if entirely to the left, |

107 | or 0 if overlap. |

108 | |

109 | The case analysis in this function is quite ugly. The fact that it's |

110 | almost 200 lines long is ridiculous. |

111 | |

112 | Ok, so here's the plan to cut it down: |

113 | |

114 | First, do a bounding line comparison on the x coordinates. This is pretty |

115 | much the common case, and should go quickly. It also takes care of the |

116 | case where both lines are horizontal. |

117 | |

118 | Then, do d0 and d1 computation, but only if a23 is nonzero. |

119 | |

120 | Finally, do d2 and d3 computation, but only if a01 is nonzero. |

121 | |

122 | Fall through to returning 0 (this will happen when both lines are |

123 | horizontal and they overlap). |

124 | */ |

125 | static int |

126 | x_order (ArtPoint z0, ArtPoint z1, ArtPoint z2, ArtPoint z3) |

127 | { |

128 | double a01, b01, c01; |

129 | double a23, b23, c23; |

130 | double d0, d1, d2, d3; |

131 | |

132 | if (z0.y == z1.y) |

133 | { |

134 | if (z2.y == z3.y) |

135 | { |

136 | double x01min, x01max; |

137 | double x23min, x23max; |

138 | |

139 | if (z0.x > z1.x) |

140 | { |

141 | x01min = z1.x; |

142 | x01max = z0.x; |

143 | } |

144 | else |

145 | { |

146 | x01min = z0.x; |

147 | x01max = z1.x; |

148 | } |

149 | |

150 | if (z2.x > z3.x) |

151 | { |

152 | x23min = z3.x; |

153 | x23max = z2.x; |

154 | } |

155 | else |

156 | { |

157 | x23min = z2.x; |

158 | x23max = z3.x; |

159 | } |

160 | |

161 | if (x23min >= x01max) return 1; |

162 | else if (x01min >= x23max) return -1; |

163 | else return 0; |

164 | } |

165 | else |

166 | { |

167 | /* z0-z1 is horizontal, z2-z3 isn't */ |

168 | a23 = z2.y - z3.y; |

169 | b23 = z3.x - z2.x; |

170 | c23 = -(z2.x * a23 + z2.y * b23); |

171 | |

172 | if (z3.y < z2.y) |

173 | { |

174 | a23 = -a23; |

175 | b23 = -b23; |

176 | c23 = -c23; |

177 | } |

178 | |

179 | d0 = trap_epsilon (a23 * z0.x + b23 * z0.y + c23); |

180 | d1 = trap_epsilon (a23 * z1.x + b23 * z1.y + c23); |

181 | |

182 | if (d0 > 0) |

183 | { |

184 | if (d1 >= 0) return 1; |

185 | else return 0; |

186 | } |

187 | else if (d0 == 0) |

188 | { |

189 | if (d1 > 0) return 1; |

190 | else if (d1 < 0) return -1; |

191 | else printf ("case 1 degenerate\n"); |

192 | return 0; |

193 | } |

194 | else /* d0 < 0 */ |

195 | { |

196 | if (d1 <= 0) return -1; |

197 | else return 0; |

198 | } |

199 | } |

200 | } |

201 | else if (z2.y == z3.y) |

202 | { |

203 | /* z2-z3 is horizontal, z0-z1 isn't */ |

204 | a01 = z0.y - z1.y; |

205 | b01 = z1.x - z0.x; |

206 | c01 = -(z0.x * a01 + z0.y * b01); |

207 | /* = -((z0.y - z1.y) * z0.x + (z1.x - z0.x) * z0.y) |

208 | = (z1.x * z0.y - z1.y * z0.x) */ |

209 | |

210 | if (z1.y < z0.y) |

211 | { |

212 | a01 = -a01; |

213 | b01 = -b01; |

214 | c01 = -c01; |

215 | } |

216 | |

217 | d2 = trap_epsilon (a01 * z2.x + b01 * z2.y + c01); |

218 | d3 = trap_epsilon (a01 * z3.x + b01 * z3.y + c01); |

219 | |

220 | if (d2 > 0) |

221 | { |

222 | if (d3 >= 0) return -1; |

223 | else return 0; |

224 | } |

225 | else if (d2 == 0) |

226 | { |

227 | if (d3 > 0) return -1; |

228 | else if (d3 < 0) return 1; |

229 | else printf ("case 2 degenerate\n"); |

230 | return 0; |

231 | } |

232 | else /* d2 < 0 */ |

233 | { |

234 | if (d3 <= 0) return 1; |

235 | else return 0; |

236 | } |

237 | } |

238 | |

239 | /* find line equations ax + by + c = 0 */ |

240 | a01 = z0.y - z1.y; |

241 | b01 = z1.x - z0.x; |

242 | c01 = -(z0.x * a01 + z0.y * b01); |

243 | /* = -((z0.y - z1.y) * z0.x + (z1.x - z0.x) * z0.y) |

244 | = -(z1.x * z0.y - z1.y * z0.x) */ |

245 | |

246 | if (a01 > 0) |

247 | { |

248 | a01 = -a01; |

249 | b01 = -b01; |

250 | c01 = -c01; |

251 | } |

252 | /* so now, (a01, b01) points to the left, thus a01 * x + b01 * y + c01 |

253 | is negative if the point lies to the right of the line */ |

254 | |

255 | d2 = trap_epsilon (a01 * z2.x + b01 * z2.y + c01); |

256 | d3 = trap_epsilon (a01 * z3.x + b01 * z3.y + c01); |

257 | if (d2 > 0) |

258 | { |

259 | if (d3 >= 0) return -1; |

260 | } |

261 | else if (d2 == 0) |

262 | { |

263 | if (d3 > 0) return -1; |

264 | else if (d3 < 0) return 1; |

265 | else |

266 | fprintf (stderr, "colinear!\n"); |

267 | } |

268 | else /* d2 < 0 */ |

269 | { |

270 | if (d3 <= 0) return 1; |

271 | } |

272 | |

273 | a23 = z2.y - z3.y; |

274 | b23 = z3.x - z2.x; |

275 | c23 = -(z2.x * a23 + z2.y * b23); |

276 | |

277 | if (a23 > 0) |

278 | { |

279 | a23 = -a23; |

280 | b23 = -b23; |

281 | c23 = -c23; |

282 | } |

283 | d0 = trap_epsilon (a23 * z0.x + b23 * z0.y + c23); |

284 | d1 = trap_epsilon (a23 * z1.x + b23 * z1.y + c23); |

285 | if (d0 > 0) |

286 | { |

287 | if (d1 >= 0) return 1; |

288 | } |

289 | else if (d0 == 0) |

290 | { |

291 | if (d1 > 0) return 1; |

292 | else if (d1 < 0) return -1; |

293 | else |

294 | fprintf (stderr, "colinear!\n"); |

295 | } |

296 | else /* d0 < 0 */ |

297 | { |

298 | if (d1 <= 0) return -1; |

299 | } |

300 | |

301 | return 0; |

302 | } |

303 | |

304 | /* similar to x_order, but to determine whether point z0 + epsilon lies to |

305 | the left of the line z2-z3 or to the right */ |

306 | static int |

307 | x_order_2 (ArtPoint z0, ArtPoint z1, ArtPoint z2, ArtPoint z3) |

308 | { |

309 | double a23, b23, c23; |

310 | double d0, d1; |

311 | |

312 | a23 = z2.y - z3.y; |

313 | b23 = z3.x - z2.x; |

314 | c23 = -(z2.x * a23 + z2.y * b23); |

315 | |

316 | if (a23 > 0) |

317 | { |

318 | a23 = -a23; |

319 | b23 = -b23; |

320 | c23 = -c23; |

321 | } |

322 | |

323 | d0 = a23 * z0.x + b23 * z0.y + c23; |

324 | |

325 | if (d0 > EPSILON) |

326 | return -1; |

327 | else if (d0 < -EPSILON) |

328 | return 1; |

329 | |

330 | d1 = a23 * z1.x + b23 * z1.y + c23; |

331 | if (d1 > EPSILON) |

332 | return -1; |

333 | else if (d1 < -EPSILON) |

334 | return 1; |

335 | |

336 | if (z0.x == z1.x && z1.x == z2.x && z2.x == z3.x) |

337 | { |

338 | art_dprint ("x_order_2: colinear and horizontally aligned!\n"); |

339 | return 0; |

340 | } |

341 | |

342 | if (z0.x <= z2.x && z1.x <= z2.x && z0.x <= z3.x && z1.x <= z3.x) |

343 | return -1; |

344 | if (z0.x >= z2.x && z1.x >= z2.x && z0.x >= z3.x && z1.x >= z3.x) |

345 | return 1; |

346 | |

347 | fprintf (stderr, "x_order_2: colinear!\n"); |

348 | return 0; |

349 | } |

350 | |

351 | #ifdef DEAD_CODE |

352 | /* Traverse the vector path, keeping it in x-sorted order. |

353 | |

354 | This routine doesn't actually do anything - it's just here for |

355 | explanatory purposes. */ |

356 | void |

357 | traverse (ArtSVP *vp) |

358 | { |

359 | int *active_segs; |

360 | int n_active_segs; |

361 | int *cursor; |

362 | int seg_idx; |

363 | double y; |

364 | int tmp1, tmp2; |

365 | int asi; |

366 | int i, j; |

367 | |

368 | active_segs = art_new (int, vp->n_segs); |

369 | cursor = art_new (int, vp->n_segs); |

370 | |

371 | n_active_segs = 0; |

372 | seg_idx = 0; |

373 | y = vp->segs[0].points[0].y; |

374 | while (seg_idx < vp->n_segs || n_active_segs > 0) |

375 | { |

376 | printf ("y = %g\n", y); |

377 | /* delete segments ending at y from active list */ |

378 | for (i = 0; i < n_active_segs; i++) |

379 | { |

380 | asi = active_segs[i]; |

381 | if (vp->segs[asi].n_points - 1 == cursor[asi] && |

382 | vp->segs[asi].points[cursor[asi]].y == y) |

383 | { |

384 | printf ("deleting %d\n", asi); |

385 | n_active_segs--; |

386 | for (j = i; j < n_active_segs; j++) |

387 | active_segs[j] = active_segs[j + 1]; |

388 | i--; |

389 | } |

390 | } |

391 | |

392 | /* insert new segments into the active list */ |

393 | while (seg_idx < vp->n_segs && y == vp->segs[seg_idx].points[0].y) |

394 | { |

395 | cursor[seg_idx] = 0; |

396 | printf ("inserting %d\n", seg_idx); |

397 | for (i = 0; i < n_active_segs; i++) |

398 | { |

399 | asi = active_segs[i]; |

400 | if (x_order (vp->segs[asi].points[cursor[asi]], |

401 | vp->segs[asi].points[cursor[asi] + 1], |

402 | vp->segs[seg_idx].points[0], |

403 | vp->segs[seg_idx].points[1]) == -1) |

404 | break; |

405 | } |

406 | tmp1 = seg_idx; |

407 | for (j = i; j < n_active_segs; j++) |

408 | { |

409 | tmp2 = active_segs[j]; |

410 | active_segs[j] = tmp1; |

411 | tmp1 = tmp2; |

412 | } |

413 | active_segs[n_active_segs] = tmp1; |

414 | n_active_segs++; |

415 | seg_idx++; |

416 | } |

417 | |

418 | /* all active segs cross the y scanline (considering segs to be |

419 | closed on top and open on bottom) */ |

420 | for (i = 0; i < n_active_segs; i++) |

421 | { |

422 | asi = active_segs[i]; |

423 | printf ("%d (%g, %g) - (%g, %g) %s\n", asi, |

424 | vp->segs[asi].points[cursor[asi]].x, |

425 | vp->segs[asi].points[cursor[asi]].y, |

426 | vp->segs[asi].points[cursor[asi] + 1].x, |

427 | vp->segs[asi].points[cursor[asi] + 1].y, |

428 | vp->segs[asi].dir ? "v" : "^"); |

429 | } |

430 | |

431 | /* advance y to the next event */ |

432 | if (n_active_segs == 0) |

433 | { |

434 | if (seg_idx < vp->n_segs) |

435 | y = vp->segs[seg_idx].points[0].y; |

436 | /* else we're done */ |

437 | } |

438 | else |

439 | { |

440 | asi = active_segs[0]; |

441 | y = vp->segs[asi].points[cursor[asi] + 1].y; |

442 | for (i = 1; i < n_active_segs; i++) |

443 | { |

444 | asi = active_segs[i]; |

445 | if (y > vp->segs[asi].points[cursor[asi] + 1].y) |

446 | y = vp->segs[asi].points[cursor[asi] + 1].y; |

447 | } |

448 | if (seg_idx < vp->n_segs && y > vp->segs[seg_idx].points[0].y) |

449 | y = vp->segs[seg_idx].points[0].y; |

450 | } |

451 | |

452 | /* advance cursors to reach new y */ |

453 | for (i = 0; i < n_active_segs; i++) |

454 | { |

455 | asi = active_segs[i]; |

456 | while (cursor[asi] < vp->segs[asi].n_points - 1 && |

457 | y >= vp->segs[asi].points[cursor[asi] + 1].y) |

458 | cursor[asi]++; |

459 | } |

460 | printf ("\n"); |

461 | } |

462 | art_free (cursor); |

463 | art_free (active_segs); |

464 | } |

465 | #endif |

466 | |

467 | /* I believe that the loop will always break with i=1. |

468 | |

469 | I think I'll want to change this from a simple sorted list to a |

470 | modified stack. ips[*][0] will get its own data structure, and |

471 | ips[*] will in general only be allocated if there is an intersection. |

472 | Finally, the segment can be traced through the initial point |

473 | (formerly ips[*][0]), backwards through the stack, and finally |

474 | to cursor + 1. |

475 | |

476 | This change should cut down on allocation bandwidth, and also |

477 | eliminate the iteration through n_ipl below. |

478 | |

479 | */ |

480 | static void |

481 | insert_ip (int seg_i, int *n_ips, int *n_ips_max, ArtPoint **ips, ArtPoint ip) |

482 | { |

483 | int i; |

484 | ArtPoint tmp1, tmp2; |

485 | int n_ipl; |

486 | ArtPoint *ipl; |

487 | |

488 | n_ipl = n_ips[seg_i]++; |

489 | if (n_ipl == n_ips_max[seg_i]) |

490 | art_expand (ips[seg_i], ArtPoint, n_ips_max[seg_i]); |

491 | ipl = ips[seg_i]; |

492 | for (i = 1; i < n_ipl; i++) |

493 | if (ipl[i].y > ip.y) |

494 | break; |

495 | tmp1 = ip; |

496 | for (; i <= n_ipl; i++) |

497 | { |

498 | tmp2 = ipl[i]; |

499 | ipl[i] = tmp1; |

500 | tmp1 = tmp2; |

501 | } |

502 | } |

503 | |

504 | /* test active segment (i - 1) against i for intersection, if |

505 | so, add intersection point to both ips lists. */ |

506 | static void |

507 | intersect_neighbors (int i, int *active_segs, |

508 | int *n_ips, int *n_ips_max, ArtPoint **ips, |

509 | int *cursor, ArtSVP *vp) |

510 | { |

511 | ArtPoint z0, z1, z2, z3; |

512 | int asi01, asi23; |

513 | ArtPoint ip; |

514 | |

515 | asi01 = active_segs[i - 1]; |

516 | |

517 | z0 = ips[asi01][0]; |

518 | if (n_ips[asi01] == 1) |

519 | z1 = vp->segs[asi01].points[cursor[asi01] + 1]; |

520 | else |

521 | z1 = ips[asi01][1]; |

522 | |

523 | asi23 = active_segs[i]; |

524 | |

525 | z2 = ips[asi23][0]; |

526 | if (n_ips[asi23] == 1) |

527 | z3 = vp->segs[asi23].points[cursor[asi23] + 1]; |

528 | else |

529 | z3 = ips[asi23][1]; |

530 | |

531 | if (intersect_lines (z0, z1, z2, z3, &ip)) |

532 | { |

533 | #ifdef VERBOSE |

534 | printf ("new intersection point: (%g, %g)\n", ip.x, ip.y); |

535 | #endif |

536 | insert_ip (asi01, n_ips, n_ips_max, ips, ip); |

537 | insert_ip (asi23, n_ips, n_ips_max, ips, ip); |

538 | } |

539 | } |

540 | |

541 | /* Add a new point to a segment in the svp. |

542 | |

543 | Here, we also check to make sure that the segments satisfy nocross. |

544 | However, this is only valuable for debugging, and could possibly be |

545 | removed. |

546 | */ |

547 | static void |

548 | svp_add_point (ArtSVP *svp, int *n_points_max, |

549 | ArtPoint p, int *seg_map, int *active_segs, int n_active_segs, |

550 | int i) |

551 | { |

552 | int asi, asi_left, asi_right; |

553 | int n_points, n_points_left, n_points_right; |

554 | ArtSVPSeg *seg; |

555 | |

556 | asi = seg_map[active_segs[i]]; |

557 | seg = &svp->segs[asi]; |

558 | n_points = seg->n_points; |

559 | /* find out whether neighboring segments share a point */ |

560 | if (i > 0) |

561 | { |

562 | asi_left = seg_map[active_segs[i - 1]]; |

563 | n_points_left = svp->segs[asi_left].n_points; |

564 | if (n_points_left > 1 && |

565 | PT_EQ (svp->segs[asi_left].points[n_points_left - 2], |

566 | svp->segs[asi].points[n_points - 1])) |

567 | { |

568 | /* ok, new vector shares a top point with segment to the left - |

569 | now, check that it satisfies ordering invariant */ |

570 | if (x_order (svp->segs[asi_left].points[n_points_left - 2], |

571 | svp->segs[asi_left].points[n_points_left - 1], |

572 | svp->segs[asi].points[n_points - 1], |

573 | p) < 1) |

574 | |

575 | { |

576 | #ifdef VERBOSE |

577 | printf ("svp_add_point: cross on left!\n"); |

578 | #endif |

579 | } |

580 | } |

581 | } |

582 | |

583 | if (i + 1 < n_active_segs) |

584 | { |

585 | asi_right = seg_map[active_segs[i + 1]]; |

586 | n_points_right = svp->segs[asi_right].n_points; |

587 | if (n_points_right > 1 && |

588 | PT_EQ (svp->segs[asi_right].points[n_points_right - 2], |

589 | svp->segs[asi].points[n_points - 1])) |

590 | { |

591 | /* ok, new vector shares a top point with segment to the right - |

592 | now, check that it satisfies ordering invariant */ |

593 | if (x_order (svp->segs[asi_right].points[n_points_right - 2], |

594 | svp->segs[asi_right].points[n_points_right - 1], |

595 | svp->segs[asi].points[n_points - 1], |

596 | p) > -1) |

597 | { |

598 | #ifdef VERBOSE |

599 | printf ("svp_add_point: cross on right!\n"); |

600 | #endif |

601 | } |

602 | } |

603 | } |

604 | if (n_points_max[asi] == n_points) |

605 | art_expand (seg->points, ArtPoint, n_points_max[asi]); |

606 | seg->points[n_points] = p; |

607 | if (p.x < seg->bbox.x0) |

608 | seg->bbox.x0 = p.x; |

609 | else if (p.x > seg->bbox.x1) |

610 | seg->bbox.x1 = p.x; |

611 | seg->bbox.y1 = p.y; |

612 | seg->n_points++; |

613 | } |

614 | |

615 | #if 0 |

616 | /* find where the segment (currently at i) is supposed to go, and return |

617 | the target index - if equal to i, then there is no crossing problem. |

618 | |

619 | "Where it is supposed to go" is defined as following: |

620 | |

621 | Delete element i, re-insert at position target (bumping everything |

622 | target and greater to the right). |

623 | */ |

624 | static int |

625 | find_crossing (int i, int *active_segs, int n_active_segs, |

626 | int *cursor, ArtPoint **ips, int *n_ips, ArtSVP *vp) |

627 | { |

628 | int asi, asi_left, asi_right; |

629 | ArtPoint p0, p1; |

630 | ArtPoint p0l, p1l; |

631 | ArtPoint p0r, p1r; |

632 | int target; |

633 | |

634 | asi = active_segs[i]; |

635 | p0 = ips[asi][0]; |

636 | if (n_ips[asi] == 1) |

637 | p1 = vp->segs[asi].points[cursor[asi] + 1]; |

638 | else |

639 | p1 = ips[asi][1]; |

640 | |

641 | for (target = i; target > 0; target--) |

642 | { |

643 | asi_left = active_segs[target - 1]; |

644 | p0l = ips[asi_left][0]; |

645 | if (n_ips[asi_left] == 1) |

646 | p1l = vp->segs[asi_left].points[cursor[asi_left] + 1]; |

647 | else |

648 | p1l = ips[asi_left][1]; |

649 | if (!PT_EQ (p0, p0l)) |

650 | break; |

651 | |

652 | #ifdef VERBOSE |

653 | printf ("point matches on left (%g, %g) - (%g, %g) x (%g, %g) - (%g, %g)!\n", |

654 | p0l.x, p0l.y, p1l.x, p1l.y, p0.x, p0.y, p1.x, p1.y); |

655 | #endif |

656 | if (x_order (p0l, p1l, p0, p1) == 1) |

657 | break; |

658 | |

659 | #ifdef VERBOSE |

660 | printf ("scanning to the left (i=%d, target=%d)\n", i, target); |

661 | #endif |

662 | } |

663 | |

664 | if (target < i) return target; |

665 | |

666 | for (; target < n_active_segs - 1; target++) |

667 | { |

668 | asi_right = active_segs[target + 1]; |

669 | p0r = ips[asi_right][0]; |

670 | if (n_ips[asi_right] == 1) |

671 | p1r = vp->segs[asi_right].points[cursor[asi_right] + 1]; |

672 | else |

673 | p1r = ips[asi_right][1]; |

674 | if (!PT_EQ (p0, p0r)) |

675 | break; |

676 | |

677 | #ifdef VERBOSE |

678 | printf ("point matches on left (%g, %g) - (%g, %g) x (%g, %g) - (%g, %g)!\n", |

679 | p0.x, p0.y, p1.x, p1.y, p0r.x, p0r.y, p1r.x, p1r.y); |

680 | #endif |

681 | if (x_order (p0r, p1r, p0, p1) == 1) |

682 | break; |

683 | |

684 | #ifdef VERBOSE |

685 | printf ("scanning to the right (i=%d, target=%d)\n", i, target); |

686 | #endif |

687 | } |

688 | |

689 | return target; |

690 | } |

691 | #endif |

692 | |

693 | /* This routine handles the case where the segment changes its position |

694 | in the active segment list. Generally, this will happen when the |

695 | segment (defined by i and cursor) shares a top point with a neighbor, |

696 | but breaks the ordering invariant. |

697 | |

698 | Essentially, this routine sorts the lines [start..end), all of which |

699 | share a top point. This is implemented as your basic insertion sort. |

700 | |

701 | This routine takes care of intersecting the appropriate neighbors, |

702 | as well. |

703 | |

704 | A first argument of -1 immediately returns, which helps reduce special |

705 | casing in the main unwind routine. |

706 | */ |

707 | static void |

708 | fix_crossing (int start, int end, int *active_segs, int n_active_segs, |

709 | int *cursor, ArtPoint **ips, int *n_ips, int *n_ips_max, |

710 | ArtSVP *vp, int *seg_map, |

711 | ArtSVP **p_new_vp, int *pn_segs_max, |

712 | int **pn_points_max) |

713 | { |

714 | int i, j; |

715 | int target; |

716 | int asi, asj; |

717 | ArtPoint p0i, p1i; |

718 | ArtPoint p0j, p1j; |

719 | int swap = 0; |

720 | #ifdef VERBOSE |

721 | int k; |

722 | #endif |

723 | ArtPoint *pts; |

724 | |

725 | #ifdef VERBOSE |

726 | printf ("fix_crossing: [%d..%d)", start, end); |

727 | for (k = 0; k < n_active_segs; k++) |

728 | printf (" %d", active_segs[k]); |

729 | printf ("\n"); |

730 | #endif |

731 | |

732 | if (start == -1) |

733 | return; |

734 | |

735 | for (i = start + 1; i < end; i++) |

736 | { |

737 | |

738 | asi = active_segs[i]; |

739 | if (cursor[asi] < vp->segs[asi].n_points - 1) { |

740 | p0i = ips[asi][0]; |

741 | if (n_ips[asi] == 1) |

742 | p1i = vp->segs[asi].points[cursor[asi] + 1]; |

743 | else |

744 | p1i = ips[asi][1]; |

745 | |

746 | for (j = i - 1; j >= start; j--) |

747 | { |

748 | asj = active_segs[j]; |

749 | if (cursor[asj] < vp->segs[asj].n_points - 1) |

750 | { |

751 | p0j = ips[asj][0]; |

752 | if (n_ips[asj] == 1) |

753 | p1j = vp->segs[asj].points[cursor[asj] + 1]; |

754 | else |

755 | p1j = ips[asj][1]; |

756 | |

757 | /* we _hope_ p0i = p0j */ |

758 | if (x_order_2 (p0j, p1j, p0i, p1i) == -1) |

759 | break; |

760 | } |

761 | } |

762 | |

763 | target = j + 1; |

764 | /* target is where active_seg[i] _should_ be in active_segs */ |

765 | |

766 | if (target != i) |

767 | { |

768 | swap = 1; |

769 | |

770 | #ifdef VERBOSE |

771 | printf ("fix_crossing: at %i should be %i\n", i, target); |

772 | #endif |

773 | |

774 | /* let's close off all relevant segments */ |

775 | for (j = i; j >= target; j--) |

776 | { |

777 | asi = active_segs[j]; |

778 | /* First conjunct: this isn't the last point in the original |

779 | segment. |

780 | |

781 | Second conjunct: this isn't the first point in the new |

782 | segment (i.e. already broken). |

783 | */ |

784 | if (cursor[asi] < vp->segs[asi].n_points - 1 && |

785 | (*p_new_vp)->segs[seg_map[asi]].n_points != 1) |

786 | { |

787 | int seg_num; |

788 | /* so break here */ |

789 | #ifdef VERBOSE |

790 | printf ("closing off %d\n", j); |

791 | #endif |

792 | |

793 | pts = art_new (ArtPoint, 16); |

794 | pts[0] = ips[asi][0]; |

795 | seg_num = art_svp_add_segment (p_new_vp, pn_segs_max, |

796 | pn_points_max, |

797 | 1, vp->segs[asi].dir, |

798 | pts, |

799 | NULL); |

800 | (*pn_points_max)[seg_num] = 16; |

801 | seg_map[asi] = seg_num; |

802 | } |

803 | } |

804 | |

805 | /* now fix the ordering in active_segs */ |

806 | asi = active_segs[i]; |

807 | for (j = i; j > target; j--) |

808 | active_segs[j] = active_segs[j - 1]; |

809 | active_segs[j] = asi; |

810 | } |

811 | } |

812 | } |

813 | if (swap && start > 0) |

814 | { |

815 | int as_start; |

816 | |

817 | as_start = active_segs[start]; |

818 | if (cursor[as_start] < vp->segs[as_start].n_points) |

819 | { |

820 | #ifdef VERBOSE |

821 | printf ("checking intersection of %d, %d\n", start - 1, start); |

822 | #endif |

823 | intersect_neighbors (start, active_segs, |

824 | n_ips, n_ips_max, ips, |

825 | cursor, vp); |

826 | } |

827 | } |

828 | |

829 | if (swap && end < n_active_segs) |

830 | { |

831 | int as_end; |

832 | |

833 | as_end = active_segs[end - 1]; |

834 | if (cursor[as_end] < vp->segs[as_end].n_points) |

835 | { |

836 | #ifdef VERBOSE |

837 | printf ("checking intersection of %d, %d\n", end - 1, end); |

838 | #endif |

839 | intersect_neighbors (end, active_segs, |

840 | n_ips, n_ips_max, ips, |

841 | cursor, vp); |

842 | } |

843 | } |

844 | if (swap) |

845 | { |

846 | #ifdef VERBOSE |

847 | printf ("fix_crossing return: [%d..%d)", start, end); |

848 | for (k = 0; k < n_active_segs; k++) |

849 | printf (" %d", active_segs[k]); |

850 | printf ("\n"); |

851 | #endif |

852 | } |

853 | } |

854 | |

855 | /* Return a new sorted vector that covers the same area as the |

856 | argument, but which satisfies the nocross invariant. |

857 | |

858 | Basically, this routine works by finding the intersection points, |

859 | and cutting the segments at those points. |

860 | |

861 | Status of this routine: |

862 | |

863 | Basic correctness: Seems ok. |

864 | |

865 | Numerical stability: known problems in the case of points falling |

866 | on lines, and colinear lines. For actual use, randomly perturbing |

867 | the vertices is currently recommended. |

868 | |

869 | Speed: pretty good, although a more efficient priority queue, as |

870 | well as bbox culling of potential intersections, are two |

871 | optimizations that could help. |

872 | |

873 | Precision: pretty good, although the numerical stability problems |

874 | make this routine unsuitable for precise calculations of |

875 | differences. |

876 | |

877 | */ |

878 | |

879 | /* Here is a more detailed description of the algorithm. It follows |

880 | roughly the structure of traverse (above), but is obviously quite |

881 | a bit more complex. |

882 | |

883 | Here are a few important data structures: |

884 | |

885 | A new sorted vector path (new_svp). |

886 | |

887 | For each (active) segment in the original, a list of intersection |

888 | points. |

889 | |

890 | Of course, the original being traversed. |

891 | |

892 | The following invariants hold (in addition to the invariants |

893 | of the traverse procedure). |

894 | |

895 | The new sorted vector path lies entirely above the y scan line. |

896 | |

897 | The new sorted vector path keeps the nocross invariant. |

898 | |

899 | For each active segment, the y scan line crosses the line from the |

900 | first to the second of the intersection points (where the second |

901 | point is cursor + 1 if there is only one intersection point). |

902 | |

903 | The list of intersection points + the (cursor + 1) point is kept |

904 | in nondecreasing y order. |

905 | |

906 | Of the active segments, none of the lines from first to second |

907 | intersection point cross the 1st ip..2nd ip line of the left or |

908 | right neighbor. (However, such a line may cross further |

909 | intersection points of the neighbors, or segments past the |

910 | immediate neighbors). |

911 | |

912 | Of the active segments, all lines from 1st ip..2nd ip are in |

913 | strictly increasing x_order (this is very similar to the invariant |

914 | of the traverse procedure, but is explicitly stated here in terms |

915 | of ips). (this basically says that nocross holds on the active |

916 | segments) |

917 | |

918 | The combination of the new sorted vector path, the path through all |

919 | the intersection points to cursor + 1, and [cursor + 1, n_points) |

920 | covers the same area as the argument. |

921 | |

922 | Another important data structure is mapping from original segment |

923 | number to new segment number. |

924 | |

925 | The algorithm is perhaps best understood as advancing the cursors |

926 | while maintaining these invariants. Here's roughly how it's done. |

927 | |

928 | When deleting segments from the active list, those segments are added |

929 | to the new sorted vector path. In addition, the neighbors may intersect |

930 | each other, so they are intersection tested (see below). |

931 | |

932 | When inserting new segments, they are intersection tested against |

933 | their neighbors. The top point of the segment becomes the first |

934 | intersection point. |

935 | |

936 | Advancing the cursor is just a bit different from the traverse |

937 | routine, as the cursor may advance through the intersection points |

938 | as well. Only when there is a single intersection point in the list |

939 | does the cursor advance in the original segment. In either case, |

940 | the new vector is intersection tested against both neighbors. It |

941 | also causes the vector over which the cursor is advancing to be |

942 | added to the new svp. |

943 | |

944 | Two steps need further clarification: |

945 | |

946 | Intersection testing: the 1st ip..2nd ip lines of the neighbors |

947 | are tested to see if they cross (using intersect_lines). If so, |

948 | then the intersection point is added to the ip list of both |

949 | segments, maintaining the invariant that the list of intersection |

950 | points is nondecreasing in y). |

951 | |

952 | Adding vector to new svp: if the new vector shares a top x |

953 | coordinate with another vector, then it is checked to see whether |

954 | it is in order. If not, then both segments are "broken," and then |

955 | restarted. Note: in the case when both segments are in the same |

956 | order, they may simply be swapped without breaking. |

957 | |

958 | For the time being, I'm going to put some of these operations into |

959 | subroutines. If it turns out to be a performance problem, I could |

960 | try to reorganize the traverse procedure so that each is only |

961 | called once, and inline them. But if it's not a performance |

962 | problem, I'll just keep it this way, because it will probably help |

963 | to make the code clearer, and I believe this code could use all the |

964 | clarity it can get. */ |

965 | /** |

966 | * art_svp_uncross: Resolve self-intersections of an svp. |

967 | * @vp: The original svp. |

968 | * |

969 | * Finds all the intersections within @vp, and constructs a new svp |

970 | * with new points added at these intersections. |

971 | * |

972 | * This routine needs to be redone from scratch with numerical robustness |

973 | * in mind. I'm working on it. |

974 | * |

975 | * Return value: The new svp. |

976 | **/ |

977 | ArtSVP * |

978 | art_svp_uncross (ArtSVP *vp) |

979 | { |

980 | int *active_segs; |

981 | int n_active_segs; |

982 | int *cursor; |

983 | int seg_idx; |

984 | double y; |

985 | int tmp1, tmp2; |

986 | int asi; |

987 | int i, j; |

988 | /* new data structures */ |

989 | /* intersection points; invariant: *ips[i] is only allocated if |

990 | i is active */ |

991 | int *n_ips, *n_ips_max; |

992 | ArtPoint **ips; |

993 | /* new sorted vector path */ |

994 | int n_segs_max, seg_num; |

995 | ArtSVP *new_vp; |

996 | int *n_points_max; |

997 | /* mapping from argument to new segment numbers - again, only valid |

998 | if active */ |

999 | int *seg_map; |

1000 | double y_curs; |

1001 | ArtPoint p_curs; |

1002 | int first_share; |

1003 | double share_x; |

1004 | ArtPoint *pts; |

1005 | |

1006 | n_segs_max = 16; |

1007 | new_vp = (ArtSVP *)art_alloc (sizeof(ArtSVP) + |

1008 | (n_segs_max - 1) * sizeof(ArtSVPSeg)); |

1009 | new_vp->n_segs = 0; |

1010 | |

1011 | if (vp->n_segs == 0) |

1012 | return new_vp; |

1013 | |

1014 | active_segs = art_new (int, vp->n_segs); |

1015 | cursor = art_new (int, vp->n_segs); |

1016 | |

1017 | seg_map = art_new (int, vp->n_segs); |

1018 | n_ips = art_new (int, vp->n_segs); |

1019 | n_ips_max = art_new (int, vp->n_segs); |

1020 | ips = art_new (ArtPoint *, vp->n_segs); |

1021 | |

1022 | n_points_max = art_new (int, n_segs_max); |

1023 | |

1024 | n_active_segs = 0; |

1025 | seg_idx = 0; |

1026 | y = vp->segs[0].points[0].y; |

1027 | while (seg_idx < vp->n_segs || n_active_segs > 0) |

1028 | { |

1029 | #ifdef VERBOSE |

1030 | printf ("y = %g\n", y); |

1031 | #endif |

1032 | |

1033 | /* maybe move deletions to end of loop (to avoid so much special |

1034 | casing on the end of a segment)? */ |

1035 | |

1036 | /* delete segments ending at y from active list */ |

1037 | for (i = 0; i < n_active_segs; i++) |

1038 | { |

1039 | asi = active_segs[i]; |

1040 | if (vp->segs[asi].n_points - 1 == cursor[asi] && |

1041 | vp->segs[asi].points[cursor[asi]].y == y) |

1042 | { |

1043 | do |

1044 | { |

1045 | #ifdef VERBOSE |

1046 | printf ("deleting %d\n", asi); |

1047 | #endif |

1048 | art_free (ips[asi]); |

1049 | n_active_segs--; |

1050 | for (j = i; j < n_active_segs; j++) |

1051 | active_segs[j] = active_segs[j + 1]; |

1052 | if (i < n_active_segs) |

1053 | asi = active_segs[i]; |

1054 | else |

1055 | break; |

1056 | } |

1057 | while (vp->segs[asi].n_points - 1 == cursor[asi] && |

1058 | vp->segs[asi].points[cursor[asi]].y == y); |

1059 | |

1060 | /* test intersection of neighbors */ |

1061 | if (i > 0 && i < n_active_segs) |

1062 | intersect_neighbors (i, active_segs, |

1063 | n_ips, n_ips_max, ips, |

1064 | cursor, vp); |

1065 | |

1066 | i--; |

1067 | } |

1068 | } |

1069 | |

1070 | /* insert new segments into the active list */ |

1071 | while (seg_idx < vp->n_segs && y == vp->segs[seg_idx].points[0].y) |

1072 | { |

1073 | #ifdef VERBOSE |

1074 | printf ("inserting %d\n", seg_idx); |

1075 | #endif |

1076 | cursor[seg_idx] = 0; |

1077 | for (i = 0; i < n_active_segs; i++) |

1078 | { |

1079 | asi = active_segs[i]; |

1080 | if (x_order_2 (vp->segs[seg_idx].points[0], |

1081 | vp->segs[seg_idx].points[1], |

1082 | vp->segs[asi].points[cursor[asi]], |

1083 | vp->segs[asi].points[cursor[asi] + 1]) == -1) |

1084 | break; |

1085 | } |

1086 | |

1087 | /* Create and initialize the intersection points data structure */ |

1088 | n_ips[seg_idx] = 1; |

1089 | n_ips_max[seg_idx] = 2; |

1090 | ips[seg_idx] = art_new (ArtPoint, n_ips_max[seg_idx]); |

1091 | ips[seg_idx][0] = vp->segs[seg_idx].points[0]; |

1092 | |

1093 | /* Start a new segment in the new vector path */ |

1094 | pts = art_new (ArtPoint, 16); |

1095 | pts[0] = vp->segs[seg_idx].points[0]; |

1096 | seg_num = art_svp_add_segment (&new_vp, &n_segs_max, |

1097 | &n_points_max, |

1098 | 1, vp->segs[seg_idx].dir, |

1099 | pts, |

1100 | NULL); |

1101 | n_points_max[seg_num] = 16; |

1102 | seg_map[seg_idx] = seg_num; |

1103 | |

1104 | tmp1 = seg_idx; |

1105 | for (j = i; j < n_active_segs; j++) |

1106 | { |

1107 | tmp2 = active_segs[j]; |

1108 | active_segs[j] = tmp1; |

1109 | tmp1 = tmp2; |

1110 | } |

1111 | active_segs[n_active_segs] = tmp1; |

1112 | n_active_segs++; |

1113 | |

1114 | if (i > 0) |

1115 | intersect_neighbors (i, active_segs, |

1116 | n_ips, n_ips_max, ips, |

1117 | cursor, vp); |

1118 | |

1119 | if (i + 1 < n_active_segs) |

1120 | intersect_neighbors (i + 1, active_segs, |

1121 | n_ips, n_ips_max, ips, |

1122 | cursor, vp); |

1123 | |

1124 | seg_idx++; |

1125 | } |

1126 | |

1127 | /* all active segs cross the y scanline (considering segs to be |

1128 | closed on top and open on bottom) */ |

1129 | #ifdef VERBOSE |

1130 | for (i = 0; i < n_active_segs; i++) |

1131 | { |

1132 | asi = active_segs[i]; |

1133 | printf ("%d ", asi); |

1134 | for (j = 0; j < n_ips[asi]; j++) |

1135 | printf ("(%g, %g) - ", |

1136 | ips[asi][j].x, |

1137 | ips[asi][j].y); |

1138 | printf ("(%g, %g) %s\n", |

1139 | vp->segs[asi].points[cursor[asi] + 1].x, |

1140 | vp->segs[asi].points[cursor[asi] + 1].y, |

1141 | vp->segs[asi].dir ? "v" : "^"); |

1142 | } |

1143 | #endif |

1144 | |

1145 | /* advance y to the next event |

1146 | Note: this is quadratic. We'd probably get decent constant |

1147 | factor speed improvement by caching the y_curs values. */ |

1148 | if (n_active_segs == 0) |

1149 | { |

1150 | if (seg_idx < vp->n_segs) |

1151 | y = vp->segs[seg_idx].points[0].y; |

1152 | /* else we're done */ |

1153 | } |

1154 | else |

1155 | { |

1156 | asi = active_segs[0]; |

1157 | if (n_ips[asi] == 1) |

1158 | y = vp->segs[asi].points[cursor[asi] + 1].y; |

1159 | else |

1160 | y = ips[asi][1].y; |

1161 | for (i = 1; i < n_active_segs; i++) |

1162 | { |

1163 | asi = active_segs[i]; |

1164 | if (n_ips[asi] == 1) |

1165 | y_curs = vp->segs[asi].points[cursor[asi] + 1].y; |

1166 | else |

1167 | y_curs = ips[asi][1].y; |

1168 | if (y > y_curs) |

1169 | y = y_curs; |

1170 | } |

1171 | if (seg_idx < vp->n_segs && y > vp->segs[seg_idx].points[0].y) |

1172 | y = vp->segs[seg_idx].points[0].y; |

1173 | } |

1174 | |

1175 | first_share = -1; |

1176 | share_x = 0; /* to avoid gcc warning, although share_x is never |

1177 | used when first_share is -1 */ |

1178 | /* advance cursors to reach new y */ |

1179 | for (i = 0; i < n_active_segs; i++) |

1180 | { |

1181 | asi = active_segs[i]; |

1182 | if (n_ips[asi] == 1) |

1183 | p_curs = vp->segs[asi].points[cursor[asi] + 1]; |

1184 | else |

1185 | p_curs = ips[asi][1]; |

1186 | if (p_curs.y == y) |

1187 | { |

1188 | svp_add_point (new_vp, n_points_max, |

1189 | p_curs, seg_map, active_segs, n_active_segs, i); |

1190 | |

1191 | n_ips[asi]--; |

1192 | for (j = 0; j < n_ips[asi]; j++) |

1193 | ips[asi][j] = ips[asi][j + 1]; |

1194 | |

1195 | if (n_ips[asi] == 0) |

1196 | { |

1197 | ips[asi][0] = p_curs; |

1198 | n_ips[asi] = 1; |

1199 | cursor[asi]++; |

1200 | } |

1201 | |

1202 | if (first_share < 0 || p_curs.x != share_x) |

1203 | { |

1204 | /* this is where crossings are detected, and if |

1205 | found, the active segments switched around. */ |

1206 | |

1207 | fix_crossing (first_share, i, |

1208 | active_segs, n_active_segs, |

1209 | cursor, ips, n_ips, n_ips_max, vp, seg_map, |

1210 | &new_vp, |

1211 | &n_segs_max, &n_points_max); |

1212 | |

1213 | first_share = i; |

1214 | share_x = p_curs.x; |

1215 | } |

1216 | |

1217 | if (cursor[asi] < vp->segs[asi].n_points - 1) |

1218 | { |

1219 | |

1220 | if (i > 0) |

1221 | intersect_neighbors (i, active_segs, |

1222 | n_ips, n_ips_max, ips, |

1223 | cursor, vp); |

1224 | |

1225 | if (i + 1 < n_active_segs) |

1226 | intersect_neighbors (i + 1, active_segs, |

1227 | n_ips, n_ips_max, ips, |

1228 | cursor, vp); |

1229 | } |

1230 | } |

1231 | else |

1232 | { |

1233 | /* not on a cursor point */ |

1234 | fix_crossing (first_share, i, |

1235 | active_segs, n_active_segs, |

1236 | cursor, ips, n_ips, n_ips_max, vp, seg_map, |

1237 | &new_vp, |

1238 | &n_segs_max, &n_points_max); |

1239 | first_share = -1; |

1240 | } |

1241 | } |

1242 | |

1243 | /* fix crossing on last shared group */ |

1244 | fix_crossing (first_share, i, |

1245 | active_segs, n_active_segs, |

1246 | cursor, ips, n_ips, n_ips_max, vp, seg_map, |

1247 | &new_vp, |

1248 | &n_segs_max, &n_points_max); |

1249 | |

1250 | #ifdef VERBOSE |

1251 | printf ("\n"); |

1252 | #endif |

1253 | } |

1254 | |

1255 | /* not necessary to sort, new segments only get added at y, which |

1256 | increases monotonically */ |

1257 | #if 0 |

1258 | qsort (&new_vp->segs, new_vp->n_segs, sizeof (svp_seg), svp_seg_compare); |

1259 | { |

1260 | int k; |

1261 | for (k = 0; k < new_vp->n_segs - 1; k++) |

1262 | { |

1263 | printf ("(%g, %g) - (%g, %g) %s (%g, %g) - (%g, %g)\n", |

1264 | new_vp->segs[k].points[0].x, |

1265 | new_vp->segs[k].points[0].y, |

1266 | new_vp->segs[k].points[1].x, |

1267 | new_vp->segs[k].points[1].y, |

1268 | svp_seg_compare (&new_vp->segs[k], &new_vp->segs[k + 1]) > 1 ? ">": "<", |

1269 | new_vp->segs[k + 1].points[0].x, |

1270 | new_vp->segs[k + 1].points[0].y, |

1271 | new_vp->segs[k + 1].points[1].x, |

1272 | new_vp->segs[k + 1].points[1].y); |

1273 | } |

1274 | } |

1275 | #endif |

1276 | |

1277 | art_free (n_points_max); |

1278 | art_free (seg_map); |

1279 | art_free (n_ips_max); |

1280 | art_free (n_ips); |

1281 | art_free (ips); |

1282 | art_free (cursor); |

1283 | art_free (active_segs); |

1284 | |

1285 | return new_vp; |

1286 | } |

1287 | |

1288 | #define noVERBOSE |

1289 | |

1290 | /* Rewind a svp satisfying the nocross invariant. |

1291 | |

1292 | The winding number of a segment is defined as the winding number of |

1293 | the points to the left while travelling in the direction of the |

1294 | segment. Therefore it preincrements and postdecrements as a scan |

1295 | line is traversed from left to right. |

1296 | |

1297 | Status of this routine: |

1298 | |

1299 | Basic correctness: Was ok in gfonted. However, this code does not |

1300 | yet compute bboxes for the resulting svp segs. |

1301 | |

1302 | Numerical stability: known problems in the case of horizontal |

1303 | segments in polygons with any complexity. For actual use, randomly |

1304 | perturbing the vertices is recommended. |

1305 | |

1306 | Speed: good. |

1307 | |

1308 | Precision: good, except that no attempt is made to remove "hair". |

1309 | Doing random perturbation just makes matters worse. |

1310 | |

1311 | */ |

1312 | /** |

1313 | * art_svp_rewind_uncrossed: Rewind an svp satisfying the nocross invariant. |

1314 | * @vp: The original svp. |

1315 | * @rule: The winding rule. |

1316 | * |

1317 | * Creates a new svp with winding number of 0 or 1 everywhere. The @rule |

1318 | * argument specifies a rule for how winding numbers in the original |

1319 | * @vp map to the winding numbers in the result. |

1320 | * |

1321 | * With @rule == ART_WIND_RULE_NONZERO, the resulting svp has a |

1322 | * winding number of 1 where @vp has a nonzero winding number. |

1323 | * |

1324 | * With @rule == ART_WIND_RULE_INTERSECT, the resulting svp has a |

1325 | * winding number of 1 where @vp has a winding number greater than |

1326 | * 1. It is useful for computing intersections. |

1327 | * |

1328 | * With @rule == ART_WIND_RULE_ODDEVEN, the resulting svp has a |

1329 | * winding number of 1 where @vp has an odd winding number. It is |

1330 | * useful for implementing the even-odd winding rule of the |

1331 | * PostScript imaging model. |

1332 | * |

1333 | * With @rule == ART_WIND_RULE_POSITIVE, the resulting svp has a |

1334 | * winding number of 1 where @vp has a positive winding number. It is |

1335 | * useful for implementing asymmetric difference. |

1336 | * |

1337 | * This routine needs to be redone from scratch with numerical robustness |

1338 | * in mind. I'm working on it. |

1339 | * |

1340 | * Return value: The new svp. |

1341 | **/ |

1342 | ArtSVP * |

1343 | art_svp_rewind_uncrossed (ArtSVP *vp, ArtWindRule rule) |

1344 | { |

1345 | int *active_segs; |

1346 | int n_active_segs; |

1347 | int *cursor; |

1348 | int seg_idx; |

1349 | double y; |

1350 | int tmp1, tmp2; |

1351 | int asi; |

1352 | int i, j; |

1353 | |

1354 | ArtSVP *new_vp; |

1355 | int n_segs_max; |

1356 | int *winding; |

1357 | int left_wind; |

1358 | int wind; |

1359 | int keep, invert; |

1360 | |

1361 | #ifdef VERBOSE |

1362 | print_svp (vp); |

1363 | #endif |

1364 | n_segs_max = 16; |

1365 | new_vp = (ArtSVP *)art_alloc (sizeof(ArtSVP) + |

1366 | (n_segs_max - 1) * sizeof(ArtSVPSeg)); |

1367 | new_vp->n_segs = 0; |

1368 | |

1369 | if (vp->n_segs == 0) |

1370 | return new_vp; |

1371 | |

1372 | winding = art_new (int, vp->n_segs); |

1373 | |

1374 | active_segs = art_new (int, vp->n_segs); |

1375 | cursor = art_new (int, vp->n_segs); |

1376 | |

1377 | n_active_segs = 0; |

1378 | seg_idx = 0; |

1379 | y = vp->segs[0].points[0].y; |

1380 | while (seg_idx < vp->n_segs || n_active_segs > 0) |

1381 | { |

1382 | #ifdef VERBOSE |

1383 | printf ("y = %g\n", y); |

1384 | #endif |

1385 | /* delete segments ending at y from active list */ |

1386 | for (i = 0; i < n_active_segs; i++) |

1387 | { |

1388 | asi = active_segs[i]; |

1389 | if (vp->segs[asi].n_points - 1 == cursor[asi] && |

1390 | vp->segs[asi].points[cursor[asi]].y == y) |

1391 | { |

1392 | #ifdef VERBOSE |

1393 | printf ("deleting %d\n", asi); |

1394 | #endif |

1395 | n_active_segs--; |

1396 | for (j = i; j < n_active_segs; j++) |

1397 | active_segs[j] = active_segs[j + 1]; |

1398 | i--; |

1399 | } |

1400 | } |

1401 | |

1402 | /* insert new segments into the active list */ |

1403 | while (seg_idx < vp->n_segs && y == vp->segs[seg_idx].points[0].y) |

1404 | { |

1405 | #ifdef VERBOSE |

1406 | printf ("inserting %d\n", seg_idx); |

1407 | #endif |

1408 | cursor[seg_idx] = 0; |

1409 | for (i = 0; i < n_active_segs; i++) |

1410 | { |

1411 | asi = active_segs[i]; |

1412 | if (x_order_2 (vp->segs[seg_idx].points[0], |

1413 | vp->segs[seg_idx].points[1], |

1414 | vp->segs[asi].points[cursor[asi]], |

1415 | vp->segs[asi].points[cursor[asi] + 1]) == -1) |

1416 | break; |

1417 | } |

1418 | |

1419 | /* Determine winding number for this segment */ |

1420 | if (i == 0) |

1421 | left_wind = 0; |

1422 | else if (vp->segs[active_segs[i - 1]].dir) |

1423 | left_wind = winding[active_segs[i - 1]]; |

1424 | else |

1425 | left_wind = winding[active_segs[i - 1]] - 1; |

1426 | |

1427 | if (vp->segs[seg_idx].dir) |

1428 | wind = left_wind + 1; |

1429 | else |

1430 | wind = left_wind; |

1431 | |

1432 | winding[seg_idx] = wind; |

1433 | |

1434 | switch (rule) |

1435 | { |

1436 | case ART_WIND_RULE_NONZERO: |

1437 | keep = (wind == 1 || wind == 0); |

1438 | invert = (wind == 0); |

1439 | break; |

1440 | case ART_WIND_RULE_INTERSECT: |

1441 | keep = (wind == 2); |

1442 | invert = 0; |

1443 | break; |

1444 | case ART_WIND_RULE_ODDEVEN: |

1445 | keep = 1; |

1446 | invert = !(wind & 1); |

1447 | break; |

1448 | case ART_WIND_RULE_POSITIVE: |

1449 | keep = (wind == 1); |

1450 | invert = 0; |

1451 | break; |

1452 | default: |

1453 | keep = 0; |

1454 | invert = 0; |

1455 | break; |

1456 | } |

1457 | |

1458 | if (keep) |

1459 | { |

1460 | ArtPoint *points, *new_points; |

1461 | int n_points; |

1462 | int new_dir; |

1463 | |

1464 | #ifdef VERBOSE |

1465 | printf ("keeping segment %d\n", seg_idx); |

1466 | #endif |

1467 | n_points = vp->segs[seg_idx].n_points; |

1468 | points = vp->segs[seg_idx].points; |

1469 | new_points = art_new (ArtPoint, n_points); |

1470 | memcpy (new_points, points, n_points * sizeof (ArtPoint)); |

1471 | new_dir = vp->segs[seg_idx].dir ^ invert; |

1472 | art_svp_add_segment (&new_vp, &n_segs_max, |

1473 | NULL, |

1474 | n_points, new_dir, new_points, |

1475 | &vp->segs[seg_idx].bbox); |

1476 | } |

1477 | |

1478 | tmp1 = seg_idx; |

1479 | for (j = i; j < n_active_segs; j++) |

1480 | { |

1481 | tmp2 = active_segs[j]; |

1482 | active_segs[j] = tmp1; |

1483 | tmp1 = tmp2; |

1484 | } |

1485 | active_segs[n_active_segs] = tmp1; |

1486 | n_active_segs++; |

1487 | seg_idx++; |

1488 | } |

1489 | |

1490 | #ifdef VERBOSE |

1491 | /* all active segs cross the y scanline (considering segs to be |

1492 | closed on top and open on bottom) */ |

1493 | for (i = 0; i < n_active_segs; i++) |

1494 | { |

1495 | asi = active_segs[i]; |

1496 | printf ("%d:%d (%g, %g) - (%g, %g) %s %d\n", asi, |

1497 | cursor[asi], |

1498 | vp->segs[asi].points[cursor[asi]].x, |

1499 | vp->segs[asi].points[cursor[asi]].y, |

1500 | vp->segs[asi].points[cursor[asi] + 1].x, |

1501 | vp->segs[asi].points[cursor[asi] + 1].y, |

1502 | vp->segs[asi].dir ? "v" : "^", |

1503 | winding[asi]); |

1504 | } |

1505 | #endif |

1506 | |

1507 | /* advance y to the next event */ |

1508 | if (n_active_segs == 0) |

1509 | { |

1510 | if (seg_idx < vp->n_segs) |

1511 | y = vp->segs[seg_idx].points[0].y; |

1512 | /* else we're done */ |

1513 | } |

1514 | else |

1515 | { |

1516 | asi = active_segs[0]; |

1517 | y = vp->segs[asi].points[cursor[asi] + 1].y; |

1518 | for (i = 1; i < n_active_segs; i++) |

1519 | { |

1520 | asi = active_segs[i]; |

1521 | if (y > vp->segs[asi].points[cursor[asi] + 1].y) |

1522 | y = vp->segs[asi].points[cursor[asi] + 1].y; |

1523 | } |

1524 | if (seg_idx < vp->n_segs && y > vp->segs[seg_idx].points[0].y) |

1525 | y = vp->segs[seg_idx].points[0].y; |

1526 | } |

1527 | |

1528 | /* advance cursors to reach new y */ |

1529 | for (i = 0; i < n_active_segs; i++) |

1530 | { |

1531 | asi = active_segs[i]; |

1532 | while (cursor[asi] < vp->segs[asi].n_points - 1 && |

1533 | y >= vp->segs[asi].points[cursor[asi] + 1].y) |

1534 | cursor[asi]++; |

1535 | } |

1536 | #ifdef VERBOSE |

1537 | printf ("\n"); |

1538 | #endif |

1539 | } |

1540 | art_free (cursor); |

1541 | art_free (active_segs); |

1542 | art_free (winding); |

1543 | |

1544 | return new_vp; |

1545 | } |

**Note:**See TracBrowser for help on using the repository browser.