1 /* fit.c: turn a bitmap representation of a curve into a list of splines.
2 Some of the ideas, but not the code, comes from the Phoenix thesis.
3 See README for the reference.
4
5 The code was partially derived from limn.
6
7 Copyright (C) 1992 Free Software Foundation, Inc.
8
9 This program is free software; you can redistribute it and/or modify
10 it under the terms of the GNU General Public License as published by
11 the Free Software Foundation; either version 2, or (at your option)
12 any later version.
13
14 This program is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 GNU General Public License for more details.
18
19 You should have received a copy of the GNU General Public License
20 along with this program; if not, write to the Free Software
21 Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */
22
23 #ifdef HAVE_CONFIG_H
24 #include "config.h"
25 #endif /* Def: HAVE_CONFIG_H */
26
27 #include "autotrace.h"
28 #include "fit.h"
29 #include "message.h"
30 #include "logreport.h"
31 #include "spline.h"
32 #include "vector.h"
33 #include "curve.h"
34 #include "pxl-outline.h"
35 #include "epsilon-equal.h"
36 #include "xstd.h"
37 #include <math.h>
38 #ifndef FLT_MAX
39 #include <limits.h>
40 #include <float.h>
41 #endif
42 #ifndef FLT_MIN
43 #include <limits.h>
44 #include <float.h>
45 #endif
46 #include <string.h>
47 #include <assert.h>
48
49 #define SQUARE(x) ((x) * (x))
50 #define CUBE(x) ((x) * (x) * (x))
51 #define ROUND(x) ((unsigned short) ((unsigned short) (x) + .5 * SIGN (x)))
52 #define SIGN(x) ((x) > 0 ? 1 : (x) < 0 ? -1 : 0)
53
54 /* We need to manipulate lists of array indices. */
55
56 typedef struct index_list
57 {
58 unsigned *data;
59 unsigned length;
60 } index_list_type;
61
62 /* The usual accessor macros. */
63 #define GET_INDEX(i_l, n) ((i_l).data[n])
64 #define INDEX_LIST_LENGTH(i_l) ((i_l).length)
65 #define GET_LAST_INDEX(i_l) ((i_l).data[INDEX_LIST_LENGTH (i_l) - 1])
66
67 static void append_index (index_list_type *, unsigned);
68 static void free_index_list (index_list_type *);
69 static index_list_type new_index_list (void);
70 static void remove_adjacent_corners (index_list_type *, unsigned, at_bool,
71 at_exception_type * exception);
72 static void change_bad_lines (spline_list_type *,
73 fitting_opts_type *);
74 static void filter (curve_type, fitting_opts_type *);
75 static void find_vectors
76 (unsigned, pixel_outline_type, vector_type *, vector_type *, unsigned);
77 static index_list_type find_corners (pixel_outline_type,
78 fitting_opts_type *,
79 at_exception_type * exception);
80 static at_real find_error (curve_type, spline_type, unsigned *,
81 at_exception_type * exception);
82 static vector_type find_half_tangent (curve_type, at_bool start, unsigned *, unsigned);
83 static void find_tangent (curve_type, at_bool, at_bool, unsigned);
84 static spline_type fit_one_spline (curve_type,
85 at_exception_type * exception);
86 static spline_list_type *fit_curve (curve_type,
87 fitting_opts_type *,
88 at_exception_type * exception);
89 static spline_list_type fit_curve_list (curve_list_type,
90 fitting_opts_type *, distance_map_type *,
91 at_exception_type * exception);
92 static spline_list_type *fit_with_least_squares (curve_type,
93 fitting_opts_type *,
94 at_exception_type * exception);
95 static spline_list_type *fit_with_line (curve_type);
96 static void remove_knee_points (curve_type, at_bool);
97 static void set_initial_parameter_values (curve_type);
98 static at_bool spline_linear_enough (spline_type *, curve_type,
99 fitting_opts_type *);
100 static curve_list_array_type split_at_corners (pixel_outline_list_type,
101 fitting_opts_type *,
102 at_exception_type * exception);
103 static at_coord real_to_int_coord (at_real_coord);
104 static at_real distance (at_real_coord, at_real_coord);
105
106 /* Get a new set of fitting options */
107 fitting_opts_type
new_fitting_opts(void)108 new_fitting_opts (void)
109 {
110 fitting_opts_type fitting_opts;
111
112 fitting_opts.background_color = NULL;
113 fitting_opts.color_count = 0;
114 fitting_opts.corner_always_threshold = (at_real) 60.0;
115 fitting_opts.corner_surround = 4;
116 fitting_opts.corner_threshold = (at_real) 100.0;
117 fitting_opts.error_threshold = (at_real) 2.0;
118 fitting_opts.filter_iterations = 4;
119 fitting_opts.line_reversion_threshold = (at_real) .01;
120 fitting_opts.line_threshold = (at_real) 1.0;
121 fitting_opts.remove_adjacent_corners = false;
122 fitting_opts.tangent_surround = 3;
123 fitting_opts.despeckle_level = 0;
124 fitting_opts.despeckle_tightness = 2.0;
125 fitting_opts.centerline = false;
126 fitting_opts.preserve_width = false;
127 fitting_opts.width_weight_factor = 6.0;
128
129 return (fitting_opts);
130 }
131
132 /* The top-level call that transforms the list of pixels in the outlines
133 of the original character to a list of spline lists fitted to those
134 pixels. */
135
136 spline_list_array_type
fitted_splines(pixel_outline_list_type pixel_outline_list,fitting_opts_type * fitting_opts,distance_map_type * dist,unsigned short width,unsigned short height,at_exception_type * exception,at_progress_func notify_progress,at_address progress_data,at_testcancel_func test_cancel,at_address testcancel_data)137 fitted_splines (pixel_outline_list_type pixel_outline_list,
138 fitting_opts_type *fitting_opts, distance_map_type *dist,
139 unsigned short width, unsigned short height,
140 at_exception_type * exception,
141 at_progress_func notify_progress,
142 at_address progress_data,
143 at_testcancel_func test_cancel,
144 at_address testcancel_data)
145 {
146 unsigned this_list;
147
148 spline_list_array_type char_splines = new_spline_list_array ();
149 curve_list_array_type curve_array = split_at_corners (pixel_outline_list,
150 fitting_opts,
151 exception);
152
153 char_splines.centerline = fitting_opts->centerline;
154 char_splines.preserve_width = fitting_opts->preserve_width;
155 char_splines.width_weight_factor = fitting_opts->width_weight_factor;
156
157 if (fitting_opts->background_color)
158 char_splines.background_color = at_color_copy(fitting_opts->background_color);
159 else
160 char_splines.background_color = NULL;
161 /* Set dummy values. Real value is set in upper context. */
162 char_splines.width = width;
163 char_splines.height = height;
164
165 for (this_list = 0; this_list < CURVE_LIST_ARRAY_LENGTH (curve_array);
166 this_list++)
167 {
168 spline_list_type curve_list_splines;
169 curve_list_type curves = CURVE_LIST_ARRAY_ELT (curve_array, this_list);
170
171 if (notify_progress)
172 notify_progress((((at_real)this_list)/((at_real)CURVE_LIST_ARRAY_LENGTH (curve_array)*(at_real)3.0) + (at_real)0.333),
173 progress_data);
174 if (test_cancel && test_cancel(testcancel_data))
175 goto cleanup;
176
177 LOG1 ("\nFitting curve list #%u:\n", this_list);
178
179 curve_list_splines = fit_curve_list (curves, fitting_opts, dist,
180 exception);
181 if (at_exception_got_fatal(exception))
182 {
183 if (char_splines.background_color)
184 at_color_free(char_splines.background_color);
185 goto cleanup;
186 }
187 curve_list_splines.clockwise = curves.clockwise;
188
189 memcpy (&(curve_list_splines.color),
190 &(O_LIST_OUTLINE(pixel_outline_list, this_list).color),
191 sizeof (color_type));
192 append_spline_list (&char_splines, curve_list_splines);
193 }
194 cleanup:
195 free_curve_list_array (&curve_array, notify_progress, progress_data);
196
197 flush_log_output ();
198
199 return char_splines;
200 }
201
202 /* Fit the list of curves CURVE_LIST to a list of splines, and return
203 it. CURVE_LIST represents a single closed paths, e.g., either the
204 inside or outside outline of an `o'. */
205
206 static spline_list_type
fit_curve_list(curve_list_type curve_list,fitting_opts_type * fitting_opts,distance_map_type * dist,at_exception_type * exception)207 fit_curve_list (curve_list_type curve_list,
208 fitting_opts_type *fitting_opts, distance_map_type *dist,
209 at_exception_type * exception)
210 {
211 curve_type curve;
212 unsigned this_curve, this_spline;
213 unsigned curve_list_length = CURVE_LIST_LENGTH (curve_list);
214 spline_list_type curve_list_splines = empty_spline_list();
215
216 curve_list_splines.open = curve_list.open;
217
218 /* Remove the extraneous ``knee'' points before filtering. Since the
219 corners have already been found, we don't need to worry about
220 removing a point that should be a corner. */
221
222 LOG ("\nRemoving knees:\n");
223 for (this_curve = 0; this_curve < curve_list_length; this_curve++)
224 {
225 LOG1 ("#%u:", this_curve);
226 remove_knee_points (CURVE_LIST_ELT (curve_list, this_curve),
227 CURVE_LIST_CLOCKWISE (curve_list));
228 }
229
230 if (dist != NULL)
231 {
232 unsigned this_point;
233 unsigned height = dist->height;
234 for (this_curve = 0; this_curve < curve_list_length; this_curve++)
235 {
236 curve = CURVE_LIST_ELT (curve_list, this_curve);
237 for (this_point = 0; this_point < CURVE_LENGTH (curve); this_point++)
238 {
239 unsigned x, y;
240 float width, w;
241 at_real_coord *coord = &CURVE_POINT(curve, this_point);
242 x = (unsigned)(coord->x);
243 y = height - (unsigned)(coord->y) - 1;
244
245 /* Each (x, y) is a point on the skeleton of the curve, which
246 might be offset from the true centerline, where the width
247 is maximal. Therefore, use as the local line width the
248 maximum distance over the neighborhood of (x, y). */
249 width = dist->d[y][x];
250 if (y - 1 >= 0)
251 {
252 if ((w = dist->d[y-1][x]) > width) width = w;
253 if (x - 1 >= 0)
254 {
255 if ((w = dist->d[y][x-1]) > width) width = w;
256 if ((w = dist->d[y-1][x-1]) > width) width = w;
257 }
258 if (x + 1 < dist->width)
259 {
260 if ((w = dist->d[y][x+1]) > width) width = w;
261 if ((w = dist->d[y-1][x+1]) > width) width = w;
262 }
263 }
264 if (y + 1 < height)
265 {
266 if ((w = dist->d[y+1][x]) > width) width = w;
267 if (x - 1 >= 0 && (w = dist->d[y+1][x-1]) > width)
268 width = w;
269 if (x + 1 < dist->width && (w = dist->d[y+1][x+1]) > width)
270 width = w;
271 }
272 coord->z = width * (fitting_opts->width_weight_factor);
273 }
274 }
275 }
276
277 /* We filter all the curves in CURVE_LIST at once; otherwise, we would
278 look at an unfiltered curve when computing tangents. */
279
280 LOG ("\nFiltering curves:\n");
281 for (this_curve = 0; this_curve < curve_list.length; this_curve++)
282 {
283 LOG1 ("#%u: ", this_curve);
284 filter (CURVE_LIST_ELT (curve_list, this_curve), fitting_opts);
285 }
286
287 /* Make the first point in the first curve also be the last point in
288 the last curve, so the fit to the whole curve list will begin and
289 end at the same point. This may cause slight errors in computing
290 the tangents and t values, but it's worth it for the continuity.
291 Of course we don't want to do this if the two points are already
292 the same, as they are if the curve is cyclic. (We don't append it
293 earlier, in `split_at_corners', because that confuses the
294 filtering.) Finally, we can't append the point if the curve is
295 exactly three points long, because we aren't adding any more data,
296 and three points isn't enough to determine a spline. Therefore,
297 the fitting will fail. */
298 curve = CURVE_LIST_ELT (curve_list, 0);
299 if (CURVE_CYCLIC (curve) == true)
300 append_point (curve, CURVE_POINT (curve, 0));
301
302 /* Finally, fit each curve in the list to a list of splines. */
303 for (this_curve = 0; this_curve < curve_list_length; this_curve++)
304 {
305 spline_list_type *curve_splines;
306 curve_type current_curve = CURVE_LIST_ELT (curve_list, this_curve);
307
308 LOG1 ("\nFitting curve #%u:\n", this_curve);
309
310 curve_splines = fit_curve (current_curve, fitting_opts, exception);
311 if (at_exception_got_fatal(exception))
312 goto cleanup;
313 else if (curve_splines == NULL)
314 {
315 LOG1 ("Could not fit curve #%u", this_curve);
316 at_exception_warning(exception, "Could not fit curve");
317 }
318 else
319 {
320 LOG1 ("Fitted splines for curve #%u:\n", this_curve);
321 for (this_spline = 0;
322 this_spline < SPLINE_LIST_LENGTH (*curve_splines);
323 this_spline++)
324 {
325 LOG1 (" %u: ", this_spline);
326 if (log_file)
327 print_spline (log_file,
328 SPLINE_LIST_ELT (*curve_splines, this_spline));
329 }
330
331 /* After fitting, we may need to change some would-be lines
332 back to curves, because they are in a list with other
333 curves. */
334 change_bad_lines (curve_splines, fitting_opts);
335
336 concat_spline_lists (&curve_list_splines, *curve_splines);
337 free_spline_list (*curve_splines);
338 free (curve_splines);
339 }
340 }
341
342 if (log_file)
343 {
344 LOG ("\nFitted splines are:\n");
345 for (this_spline = 0; this_spline < SPLINE_LIST_LENGTH (curve_list_splines);
346 this_spline++)
347 {
348 LOG1 (" %u: ", this_spline);
349 print_spline (log_file, SPLINE_LIST_ELT (curve_list_splines,
350 this_spline));
351 }
352 }
353 cleanup:
354 return curve_list_splines;
355 }
356
357 /* Transform a set of locations to a list of splines (the fewer the
358 better). We are guaranteed that CURVE does not contain any corners.
359 We return NULL if we cannot fit the points at all. */
360
361 static spline_list_type *
fit_curve(curve_type curve,fitting_opts_type * fitting_opts,at_exception_type * exception)362 fit_curve (curve_type curve, fitting_opts_type *fitting_opts,
363 at_exception_type *exception)
364 {
365 spline_list_type *fittedsplines;
366
367 if (CURVE_LENGTH (curve) < 2)
368 {
369 LOG ("Tried to fit curve with less than two points");
370 at_exception_warning(exception,
371 "Tried to fit curve with less than two points");
372 return NULL;
373 }
374
375 /* Do we have enough points to fit with a spline? */
376 fittedsplines = CURVE_LENGTH (curve) < 4
377 ? fit_with_line (curve)
378 : fit_with_least_squares (curve, fitting_opts, exception);
379
380 return fittedsplines;
381 }
382
383 /* As mentioned above, the first step is to find the corners in
384 PIXEL_LIST, the list of points. (Presumably we can't fit a single
385 spline around a corner.) The general strategy is to look through all
386 the points, remembering which we want to consider corners. Then go
387 through that list, producing the curve_list. This is dictated by the
388 fact that PIXEL_LIST does not necessarily start on a corner---it just
389 starts at the character's first outline pixel, going left-to-right,
390 top-to-bottom. But we want all our splines to start and end on real
391 corners.
392
393 For example, consider the top of a capital `C' (this is in cmss20):
394 x
395 ***********
396 ******************
397
398 PIXEL_LIST will start at the pixel below the `x'. If we considered
399 this pixel a corner, we would wind up matching a very small segment
400 from there to the end of the line, probably as a straight line, which
401 is certainly not what we want.
402
403 PIXEL_LIST has one element for each closed outline on the character.
404 To preserve this information, we return an array of curve_lists, one
405 element (which in turn consists of several curves, one between each
406 pair of corners) for each element in PIXEL_LIST. */
407
408 static curve_list_array_type
split_at_corners(pixel_outline_list_type pixel_list,fitting_opts_type * fitting_opts,at_exception_type * exception)409 split_at_corners (pixel_outline_list_type pixel_list, fitting_opts_type *fitting_opts,
410 at_exception_type * exception)
411 {
412 unsigned this_pixel_o;
413 curve_list_array_type curve_array = new_curve_list_array ();
414
415 LOG ("\nFinding corners:\n");
416
417 for (this_pixel_o = 0; this_pixel_o < O_LIST_LENGTH (pixel_list);
418 this_pixel_o++)
419 {
420 curve_type curve, first_curve;
421 index_list_type corner_list;
422 unsigned p, this_corner;
423 curve_list_type curve_list = new_curve_list ();
424 pixel_outline_type pixel_o = O_LIST_OUTLINE (pixel_list, this_pixel_o);
425
426 CURVE_LIST_CLOCKWISE (curve_list) = O_CLOCKWISE (pixel_o);
427 curve_list.open = pixel_o.open;
428
429 LOG1 ("#%u:", this_pixel_o);
430
431 /* If the outline does not have enough points, we can't do
432 anything. The endpoints of the outlines are automatically
433 corners. We need at least `corner_surround' more pixels on
434 either side of a point before it is conceivable that we might
435 want another corner. */
436 if (O_LENGTH (pixel_o) > fitting_opts->corner_surround * 2 + 2)
437 corner_list = find_corners (pixel_o, fitting_opts,
438 exception);
439
440 else {
441 int surround;
442 if ((surround = (int)(O_LENGTH (pixel_o) - 3) / 2) >= 2)
443 {
444 unsigned save_corner_surround = fitting_opts->corner_surround;
445 fitting_opts->corner_surround = surround;
446 corner_list = find_corners (pixel_o, fitting_opts,
447 exception);
448 fitting_opts->corner_surround = save_corner_surround;
449 }
450 else
451 {
452 corner_list.length = 0;
453 corner_list.data = NULL;
454 }
455 }
456
457 /* Remember the first curve so we can make it be the `next' of the
458 last one. (And vice versa.) */
459 first_curve = new_curve ();
460
461 curve = first_curve;
462
463 if (corner_list.length == 0)
464 { /* No corners. Use all of the pixel outline as the curve. */
465 for (p = 0; p < O_LENGTH (pixel_o); p++)
466 append_pixel (curve, O_COORDINATE (pixel_o, p));
467
468 if (curve_list.open == true)
469 CURVE_CYCLIC (curve) = false;
470 else
471 CURVE_CYCLIC (curve) = true;
472 }
473 else
474 { /* Each curve consists of the points between (inclusive) each pair
475 of corners. */
476 for (this_corner = 0; this_corner < corner_list.length - 1;
477 this_corner++)
478 {
479 curve_type previous_curve = curve;
480 unsigned corner = GET_INDEX (corner_list, this_corner);
481 unsigned next_corner = GET_INDEX (corner_list, this_corner + 1);
482
483 for (p = corner; p <= next_corner; p++)
484 append_pixel (curve, O_COORDINATE (pixel_o, p));
485
486 append_curve (&curve_list, curve);
487 curve = new_curve ();
488 NEXT_CURVE (previous_curve) = curve;
489 PREVIOUS_CURVE (curve) = previous_curve;
490 }
491
492 /* The last curve is different. It consists of the points
493 (inclusive) between the last corner and the end of the list,
494 and the beginning of the list and the first corner. */
495 for (p = GET_LAST_INDEX (corner_list); p < O_LENGTH (pixel_o);
496 p++)
497 append_pixel (curve, O_COORDINATE (pixel_o, p));
498
499 if (!pixel_o.open)
500 {
501 for (p = 0; p <= GET_INDEX (corner_list, 0); p++)
502 append_pixel (curve, O_COORDINATE (pixel_o, p));
503 }
504 else
505 {
506 curve_type last_curve = PREVIOUS_CURVE(curve);
507 PREVIOUS_CURVE(first_curve) = NULL;
508 if (last_curve) NEXT_CURVE(last_curve) = NULL;
509 }
510 }
511
512 LOG1 (" [%u].\n", corner_list.length);
513 free_index_list (&corner_list);
514
515 /* Add `curve' to the end of the list, updating the pointers in
516 the chain. */
517 append_curve (&curve_list, curve);
518 NEXT_CURVE (curve) = first_curve;
519 PREVIOUS_CURVE (first_curve) = curve;
520
521 /* And now add the just-completed curve list to the array. */
522 append_curve_list (&curve_array, curve_list);
523 } /* End of considering each pixel outline. */
524
525 return curve_array;
526 }
527
528 /* We consider a point to be a corner if (1) the angle defined by the
529 `corner_surround' points coming into it and going out from it is less
530 than `corner_threshold' degrees, and no point within
531 `corner_surround' points has a smaller angle; or (2) the angle is less
532 than `corner_always_threshold' degrees.
533
534 Because of the different cases, it is convenient to have the
535 following macro to append a corner on to the list we return. The
536 character argument C is simply so that the different cases can be
537 distinguished in the log file. */
538
539 #define APPEND_CORNER(index, angle, c) \
540 do \
541 { \
542 append_index (&corner_list, index); \
543 LOG4 (" (%d,%d)%c%.3f", \
544 O_COORDINATE (pixel_outline, index).x, \
545 O_COORDINATE (pixel_outline, index).y, \
546 c, angle); \
547 } \
548 while (0)
549
550 static index_list_type
find_corners(pixel_outline_type pixel_outline,fitting_opts_type * fitting_opts,at_exception_type * exception)551 find_corners (pixel_outline_type pixel_outline,
552 fitting_opts_type *fitting_opts,
553 at_exception_type * exception)
554 {
555 unsigned p, start_p, end_p;
556 index_list_type corner_list = new_index_list ();
557
558 start_p = 0;
559 end_p = O_LENGTH(pixel_outline) - 1;
560 if (pixel_outline.open)
561 {
562 if (end_p <= fitting_opts->corner_surround * 2)
563 return corner_list;
564 APPEND_CORNER(0, 0.0, '@');
565 start_p += fitting_opts->corner_surround;
566 end_p -= fitting_opts->corner_surround;
567 }
568
569 /* Consider each pixel on the outline in turn. */
570 for (p = start_p; p <= end_p; p++)
571 {
572 at_real corner_angle;
573 vector_type in_vector, out_vector;
574
575 /* Check if the angle is small enough. */
576 find_vectors (p, pixel_outline, &in_vector, &out_vector,
577 fitting_opts->corner_surround);
578 corner_angle = Vangle (in_vector, out_vector, exception);
579 if (at_exception_got_fatal(exception))
580 goto cleanup;
581
582 if (fabs (corner_angle) <= fitting_opts->corner_threshold)
583 {
584 /* We want to keep looking, instead of just appending the
585 first pixel we find with a small enough angle, since there
586 might be another corner within `corner_surround' pixels, with
587 a smaller angle. If that is the case, we want that one. */
588 at_real best_corner_angle = corner_angle;
589 unsigned best_corner_index = p;
590 index_list_type equally_good_list = new_index_list ();
591 /* As we come into the loop, `p' is the index of the point
592 that has an angle less than `corner_angle'. We use `i' to
593 move through the pixels next to that, and `q' for moving
594 through the adjacent pixels to each `p'. */
595 unsigned q = p;
596 unsigned i = p + 1;
597
598 while (true)
599 {
600 /* Perhaps the angle is sufficiently small that we want to
601 consider this a corner, even if it's not the best
602 (unless we've already wrapped around in the search,
603 i.e., `q<i', in which case we have already added the
604 corner, and we don't want to add it again). We want to
605 do this check on the first candidate we find, as well
606 as the others in the loop, hence this comes before the
607 stopping condition. */
608 if (corner_angle <= fitting_opts->corner_always_threshold && q >= p)
609 APPEND_CORNER (q, corner_angle, '\\');
610
611 /* Exit the loop if we've looked at `corner_surround'
612 pixels past the best one we found, or if we've looked
613 at all the pixels. */
614 if (i >= best_corner_index + fitting_opts->corner_surround
615 || i >= O_LENGTH (pixel_outline))
616 break;
617
618 /* Check the angle. */
619 q = i % O_LENGTH (pixel_outline);
620 find_vectors (q, pixel_outline, &in_vector, &out_vector,
621 fitting_opts->corner_surround);
622 corner_angle = Vangle (in_vector, out_vector, exception);
623 if (at_exception_got_fatal(exception))
624 goto cleanup;
625
626 /* If we come across a corner that is just as good as the
627 best one, we should make it a corner, too. This
628 happens, for example, at the points on the `W' in some
629 typefaces, where the ``points'' are flat. */
630 if (epsilon_equal (corner_angle, best_corner_angle))
631 append_index (&equally_good_list, q);
632
633 else if (corner_angle < best_corner_angle)
634 {
635 best_corner_angle = corner_angle;
636 /* We want to check `corner_surround' pixels beyond the
637 new best corner. */
638 i = best_corner_index = q;
639 free_index_list (&equally_good_list);
640 equally_good_list = new_index_list ();
641 }
642
643 i++;
644 }
645
646 /* After we exit the loop, `q' is the index of the last point
647 we checked. We have already added the corner if
648 `best_corner_angle' is less than `corner_always_threshold'.
649 Again, if we've already wrapped around, we don't want to
650 add the corner again. */
651 if (best_corner_angle > fitting_opts->corner_always_threshold
652 && best_corner_index >= p)
653 {
654 unsigned j;
655
656 APPEND_CORNER (best_corner_index, best_corner_angle, '/');
657
658 for (j = 0; j < INDEX_LIST_LENGTH (equally_good_list); j++)
659 APPEND_CORNER (GET_INDEX (equally_good_list, j),
660 best_corner_angle, '@');
661 }
662 free_index_list (&equally_good_list);
663
664 /* If we wrapped around in our search, we're done; otherwise,
665 we don't want the outer loop to look at the pixels that we
666 already looked at in searching for the best corner. */
667 p = (q < p) ? O_LENGTH (pixel_outline) : q;
668 } /* End of searching for the best corner. */
669 } /* End of considering each pixel. */
670
671 if (INDEX_LIST_LENGTH (corner_list) > 0)
672 /* We never want two corners next to each other, since the
673 only way to fit such a ``curve'' would be with a straight
674 line, which usually interrupts the continuity dreadfully. */
675 remove_adjacent_corners (&corner_list,
676 O_LENGTH (pixel_outline) - (pixel_outline.open ? 2 : 1),
677 fitting_opts->remove_adjacent_corners,
678 exception);
679 cleanup:
680 return corner_list;
681 }
682
683 /* Return the difference vectors coming in and going out of the outline
684 OUTLINE at the point whose index is TEST_INDEX. In Phoenix,
685 Schneider looks at a single point on either side of the point we're
686 considering. That works for him because his points are not touching.
687 But our points *are* touching, and so we have to look at
688 `corner_surround' points on either side, to get a better picture of
689 the outline's shape. */
690
691 static void
find_vectors(unsigned test_index,pixel_outline_type outline,vector_type * in,vector_type * out,unsigned corner_surround)692 find_vectors (unsigned test_index, pixel_outline_type outline,
693 vector_type *in, vector_type *out,
694 unsigned corner_surround)
695 {
696 int i;
697 unsigned n_done;
698 at_coord candidate = O_COORDINATE (outline, test_index);
699
700 in->dx = in->dy = in->dz = 0.0;
701 out->dx = out->dy = out->dz = 0.0;
702
703 /* Add up the differences from p of the `corner_surround' points
704 before p. */
705 for (i = O_PREV (outline, test_index), n_done = 0; n_done
706 < corner_surround;
707 i = O_PREV (outline, i), n_done++)
708 *in = Vadd (*in, IPsubtract (O_COORDINATE (outline, i), candidate));
709
710 /* And the points after p. */
711 for (i = O_NEXT (outline, test_index), n_done = 0; n_done <
712 corner_surround; i = O_NEXT (outline, i), n_done++)
713 *out = Vadd (*out, IPsubtract (O_COORDINATE (outline, i), candidate));
714 }
715
716 /* Remove adjacent points from the index list LIST. We do this by first
717 sorting the list and then running through it. Since these lists are
718 quite short, a straight selection sort (e.g., p.139 of the Art of
719 Computer Programming, vol.3) is good enough. LAST_INDEX is the index
720 of the last pixel on the outline, i.e., the next one is the first
721 pixel. We need this for checking the adjacency of the last corner.
722
723 We need to do this because the adjacent corners turn into
724 two-pixel-long curves, which can only be fit by straight lines. */
725
726 static void
remove_adjacent_corners(index_list_type * list,unsigned last_index,at_bool remove_adj_corners,at_exception_type * exception)727 remove_adjacent_corners (index_list_type *list, unsigned last_index,
728 at_bool remove_adj_corners,
729 at_exception_type * exception)
730
731 {
732 unsigned j;
733 unsigned last;
734 index_list_type new_list = new_index_list ();
735
736 for (j = INDEX_LIST_LENGTH (*list) - 1; j > 0; j--)
737 {
738 unsigned search;
739 unsigned temp;
740 /* Find maximal element below `j'. */
741 unsigned max_index = j;
742
743 for (search = 0; search < j; search++)
744 if (GET_INDEX (*list, search) > GET_INDEX (*list, max_index))
745 max_index = search;
746
747 if (max_index != j)
748 {
749 temp = GET_INDEX (*list, j);
750 GET_INDEX (*list, j) = GET_INDEX (*list, max_index);
751 GET_INDEX (*list, max_index) = temp;
752
753 /* xx -- really have to sort? */
754 LOG ("needed exchange");
755 at_exception_warning(exception, "needed exchange");
756 }
757 }
758
759 /* The list is sorted. Now look for adjacent entries. Each time
760 through the loop we insert the current entry and, if appropriate,
761 the next entry. */
762 for (j = 0; j < INDEX_LIST_LENGTH (*list) - 1; j++)
763 {
764 unsigned current = GET_INDEX (*list, j);
765 unsigned next = GET_INDEX (*list, j + 1);
766
767 /* We should never have inserted the same element twice. */
768 /* assert (current != next); */
769
770 if ((remove_adj_corners) && ((next == current + 1) || (next == current)))
771 j++;
772
773 append_index (&new_list, current);
774 }
775
776 /* Don't append the last element if it is 1) adjacent to the previous
777 one; or 2) adjacent to the very first one. */
778 last = GET_LAST_INDEX (*list);
779 if (INDEX_LIST_LENGTH (new_list) == 0
780 || !(last == GET_LAST_INDEX (new_list) + 1
781 || (last == last_index && GET_INDEX (*list, 0) == 0)))
782 append_index (&new_list, last);
783
784 free_index_list (list);
785 *list = new_list;
786 }
787
788 /* A ``knee'' is a point which forms a ``right angle'' with its
789 predecessor and successor. See the documentation (the `Removing
790 knees' section) for an example and more details.
791
792 The argument CLOCKWISE tells us which direction we're moving. (We
793 can't figure that information out from just the single segment with
794 which we are given to work.)
795
796 We should never find two consecutive knees.
797
798 Since the first and last points are corners (unless the curve is
799 cyclic), it doesn't make sense to remove those. */
800
801 /* This evaluates to true if the vector V is zero in one direction and
802 nonzero in the other. */
803 #define ONLY_ONE_ZERO(v) \
804 (((v).dx == 0.0 && (v).dy != 0.0) || ((v).dy == 0.0 && (v).dx != 0.0))
805
806 /* There are four possible cases for knees, one for each of the four
807 corners of a rectangle; and then the cases differ depending on which
808 direction we are going around the curve. The tests are listed here
809 in the order of upper left, upper right, lower right, lower left.
810 Perhaps there is some simple pattern to the
811 clockwise/counterclockwise differences, but I don't see one. */
812 #define CLOCKWISE_KNEE(prev_delta, next_delta) \
813 ((prev_delta.dx == -1.0 && next_delta.dy == 1.0) \
814 || (prev_delta.dy == 1.0 && next_delta.dx == 1.0) \
815 || (prev_delta.dx == 1.0 && next_delta.dy == -1.0) \
816 || (prev_delta.dy == -1.0 && next_delta.dx == -1.0))
817
818 #define COUNTERCLOCKWISE_KNEE(prev_delta, next_delta) \
819 ((prev_delta.dy == 1.0 && next_delta.dx == -1.0) \
820 || (prev_delta.dx == 1.0 && next_delta.dy == 1.0) \
821 || (prev_delta.dy == -1.0 && next_delta.dx == 1.0) \
822 || (prev_delta.dx == -1.0 && next_delta.dy == -1.0))
823
824 static void
remove_knee_points(curve_type curve,at_bool clockwise)825 remove_knee_points (curve_type curve, at_bool clockwise)
826 {
827 unsigned i;
828 unsigned offset = (CURVE_CYCLIC (curve) == true) ? 0 : 1;
829 at_coord previous
830 = real_to_int_coord (CURVE_POINT (curve, CURVE_PREV (curve, offset)));
831 curve_type trimmed_curve = copy_most_of_curve (curve);
832
833 if (CURVE_CYCLIC (curve) == false)
834 append_pixel (trimmed_curve, real_to_int_coord (CURVE_POINT (curve, 0)));
835
836 for (i = offset; i < CURVE_LENGTH (curve) - offset; i++)
837 {
838 at_coord current
839 = real_to_int_coord (CURVE_POINT (curve, i));
840 at_coord next
841 = real_to_int_coord (CURVE_POINT (curve, CURVE_NEXT (curve, i)));
842 vector_type prev_delta = IPsubtract (previous, current);
843 vector_type next_delta = IPsubtract (next, current);
844
845 if (ONLY_ONE_ZERO (prev_delta) && ONLY_ONE_ZERO (next_delta)
846 && ((clockwise && CLOCKWISE_KNEE (prev_delta, next_delta))
847 || (!clockwise
848 && COUNTERCLOCKWISE_KNEE (prev_delta, next_delta))))
849 LOG2 (" (%d,%d)", current.x, current.y);
850 else
851 {
852 previous = current;
853 append_pixel (trimmed_curve, current);
854 }
855 }
856
857 if (CURVE_CYCLIC (curve) == false)
858 append_pixel (trimmed_curve, real_to_int_coord (LAST_CURVE_POINT (curve)));
859
860 if (CURVE_LENGTH (trimmed_curve) == CURVE_LENGTH (curve))
861 LOG (" (none)");
862
863 LOG (".\n");
864
865 free_curve (curve);
866 *curve = *trimmed_curve;
867 free (trimmed_curve); /* free_curve? --- Masatake */
868 }
869
870 /* Smooth the curve by adding in neighboring points. Do this
871 `filter_iterations' times. But don't change the corners. */
872
873 static void
filter(curve_type curve,fitting_opts_type * fitting_opts)874 filter (curve_type curve, fitting_opts_type *fitting_opts)
875 {
876 unsigned iteration, this_point;
877 unsigned offset = (CURVE_CYCLIC (curve) == true) ? 0 : 1;
878 at_real_coord prev_new_point;
879
880 /* We must have at least three points---the previous one, the current
881 one, and the next one. But if we don't have at least five, we will
882 probably collapse the curve down onto a single point, which means
883 we won't be able to fit it with a spline. */
884 if (CURVE_LENGTH (curve) < 5)
885 {
886 LOG1 ("Length is %u, not enough to filter.\n", CURVE_LENGTH (curve));
887 return;
888 }
889
890 prev_new_point.x = FLT_MAX;
891 prev_new_point.y = FLT_MAX;
892 prev_new_point.z = FLT_MAX;
893
894 for (iteration = 0; iteration < fitting_opts->filter_iterations;
895 iteration++)
896 {
897 curve_type newcurve = copy_most_of_curve (curve);
898 at_bool collapsed = false;
899
900 /* Keep the first point on the curve. */
901 if (offset)
902 append_point (newcurve, CURVE_POINT (curve, 0));
903
904 for (this_point = offset; this_point < CURVE_LENGTH (curve) - offset;
905 this_point++)
906 {
907 vector_type in, out, sum;
908 at_real_coord new_point;
909
910 /* Calculate the vectors in and out, computed by looking at n points
911 on either side of this_point. Experimental it was found that 2 is
912 optimal. */
913
914 signed int prev, prevprev; /* have to be signed */
915 unsigned int next, nextnext;
916 at_real_coord candidate = CURVE_POINT (curve, this_point);
917
918 prev = CURVE_PREV (curve, this_point);
919 prevprev = CURVE_PREV (curve, prev);
920 next = CURVE_NEXT (curve, this_point);
921 nextnext = CURVE_NEXT (curve, next);
922
923 /* Add up the differences from p of the `surround' points
924 before p. */
925 in.dx = in.dy = in.dz = 0.0;
926
927 in = Vadd (in, Psubtract (CURVE_POINT (curve, prev), candidate));
928 if (prevprev >= 0)
929 in = Vadd (in, Psubtract (CURVE_POINT (curve, prevprev), candidate));
930
931 /* And the points after p. Don't use more points after p than we
932 ended up with before it. */
933 out.dx = out.dy = out.dz = 0.0;
934
935 out = Vadd (out, Psubtract (CURVE_POINT (curve, next), candidate));
936 if (nextnext < CURVE_LENGTH (curve))
937 out = Vadd (out, Psubtract (CURVE_POINT (curve, nextnext), candidate));
938
939 /* Start with the old point. */
940 new_point = candidate;
941 sum = Vadd (in, out);
942 /* We added 2*n+2 points, so we have to divide the sum by 2*n+2 */
943 new_point.x += sum.dx / 6;
944 new_point.y += sum.dy / 6;
945 new_point.z += sum.dz / 6;
946 if (fabs (prev_new_point.x - new_point.x) < 0.3
947 && fabs (prev_new_point.y - new_point.y) < 0.3
948 && fabs (prev_new_point.z - new_point.z) < 0.3)
949 {
950 collapsed = true;
951 break;
952 }
953
954
955 /* Put the newly computed point into a separate curve, so it
956 doesn't affect future computation (on this iteration). */
957 append_point (newcurve, prev_new_point = new_point);
958 }
959
960 if (collapsed)
961 free_curve (newcurve);
962 else
963 {
964 /* Just as with the first point, we have to keep the last point. */
965 if (offset)
966 append_point (newcurve, LAST_CURVE_POINT (curve));
967
968 /* Set the original curve to the newly filtered one, and go again. */
969 free_curve (curve);
970 *curve = *newcurve;
971 }
972 free (newcurve);
973 }
974
975 log_curve (curve, false);
976 }
977
978
979 /* This routine returns the curve fitted to a straight line in a very
980 simple way: make the first and last points on the curve be the
981 endpoints of the line. This simplicity is justified because we are
982 called only on very short curves. */
983
984 static spline_list_type *
fit_with_line(curve_type curve)985 fit_with_line (curve_type curve)
986 {
987 spline_type line;
988
989 LOG ("Fitting with straight line:\n");
990
991 SPLINE_DEGREE (line) = LINEARTYPE;
992 START_POINT (line) = CONTROL1 (line) = CURVE_POINT (curve, 0);
993 END_POINT (line) = CONTROL2 (line) = LAST_CURVE_POINT (curve);
994
995 /* Make sure that this line is never changed to a cubic. */
996 SPLINE_LINEARITY (line) = 0;
997
998 if (log_file)
999 {
1000 LOG (" ");
1001 print_spline (log_file, line);
1002 }
1003
1004 return new_spline_list_with_spline (line);
1005 }
1006
1007 /* The least squares method is well described in Schneider's thesis.
1008 Briefly, we try to fit the entire curve with one spline. If that
1009 fails, we subdivide the curve. */
1010
1011 static spline_list_type *
fit_with_least_squares(curve_type curve,fitting_opts_type * fitting_opts,at_exception_type * exception)1012 fit_with_least_squares (curve_type curve, fitting_opts_type *fitting_opts,
1013 at_exception_type * exception)
1014
1015 {
1016 at_real error = 0, best_error = FLT_MAX;
1017 spline_type spline, best_spline;
1018 spline_list_type *spline_list = NULL;
1019 unsigned worst_point = 0;
1020 at_real previous_error = FLT_MAX;
1021
1022 LOG ("\nFitting with least squares:\n");
1023
1024 /* Phoenix reduces the number of points with a ``linear spline
1025 technique''. But for fitting letterforms, that is
1026 inappropriate. We want all the points we can get. */
1027
1028 /* It makes no difference whether we first set the `t' values or
1029 find the tangents. This order makes the documentation a little
1030 more coherent. */
1031
1032 LOG ("Finding tangents:\n");
1033 find_tangent (curve, /* to_start */ true, /* cross_curve */ false,
1034 fitting_opts->tangent_surround);
1035 find_tangent (curve, /* to_start */ false, /* cross_curve */ false,
1036 fitting_opts->tangent_surround);
1037
1038 set_initial_parameter_values (curve);
1039
1040 /* Now we loop, subdividing, until CURVE has
1041 been fit. */
1042 while (true)
1043 {
1044 spline = best_spline = fit_one_spline (curve, exception);
1045 if (at_exception_got_fatal(exception))
1046 goto cleanup;
1047
1048 if (SPLINE_DEGREE(spline) == LINEARTYPE)
1049 LOG (" fitted to line:\n");
1050 else
1051 LOG (" fitted to spline:\n");
1052
1053 if (log_file)
1054 {
1055 LOG (" ");
1056 print_spline (log_file, spline);
1057 }
1058
1059 if (SPLINE_DEGREE(spline) == LINEARTYPE)
1060 break;
1061
1062 error = find_error (curve, spline, &worst_point, exception);
1063 if (error <= previous_error)
1064 {
1065 best_error = error;
1066 best_spline = spline;
1067 }
1068 break;
1069 }
1070
1071 if (SPLINE_DEGREE(spline) == LINEARTYPE)
1072 {
1073 spline_list = new_spline_list_with_spline (spline);
1074 LOG1 ("Accepted error of %.3f.\n", error);
1075 return (spline_list);
1076 }
1077
1078 /* Go back to the best fit. */
1079 spline = best_spline;
1080 error = best_error;
1081
1082 if (error < fitting_opts->error_threshold && CURVE_CYCLIC(curve) == false)
1083 {
1084 /* The points were fitted with a
1085 spline. We end up here whenever a fit is accepted. We have
1086 one more job: see if the ``curve'' that was fit should really
1087 be a straight line. */
1088 if (spline_linear_enough (&spline, curve, fitting_opts))
1089 {
1090 SPLINE_DEGREE (spline) = LINEARTYPE;
1091 LOG ("Changed to line.\n");
1092 }
1093 spline_list = new_spline_list_with_spline (spline);
1094 LOG1 ("Accepted error of %.3f.\n", error);
1095 }
1096 else
1097 {
1098 /* We couldn't fit the curve acceptably, so subdivide. */
1099 unsigned subdivision_index;
1100 spline_list_type *left_spline_list;
1101 spline_list_type *right_spline_list;
1102 curve_type left_curve = new_curve ();
1103 curve_type right_curve = new_curve ();
1104
1105 /* Keep the linked list of curves intact. */
1106 NEXT_CURVE (right_curve) = NEXT_CURVE (curve);
1107 PREVIOUS_CURVE (right_curve) = left_curve;
1108 NEXT_CURVE (left_curve) = right_curve;
1109 PREVIOUS_CURVE (left_curve) = curve;
1110 NEXT_CURVE (curve) = left_curve;
1111
1112 LOG1 ("\nSubdividing (error %.3f):\n", error);
1113 LOG3 (" Original point: (%.3f,%.3f), #%u.\n",
1114 CURVE_POINT (curve, worst_point).x,
1115 CURVE_POINT (curve, worst_point).y, worst_point);
1116 subdivision_index = worst_point;
1117 LOG3 (" Final point: (%.3f,%.3f), #%u.\n",
1118 CURVE_POINT (curve, subdivision_index).x,
1119 CURVE_POINT (curve, subdivision_index).y, subdivision_index);
1120
1121 /* The last point of the left-hand curve will also be the first
1122 point of the right-hand curve. */
1123 CURVE_LENGTH (left_curve) = subdivision_index + 1;
1124 CURVE_LENGTH (right_curve) = CURVE_LENGTH (curve) - subdivision_index;
1125 left_curve->point_list = curve->point_list;
1126 right_curve->point_list = curve->point_list + subdivision_index;
1127
1128 /* We want to use the tangents of the curve which we are
1129 subdividing for the start tangent for left_curve and the
1130 end tangent for right_curve. */
1131 CURVE_START_TANGENT (left_curve) = CURVE_START_TANGENT (curve);
1132 CURVE_END_TANGENT (right_curve) = CURVE_END_TANGENT (curve);
1133
1134 /* We have to set up the two curves before finding the tangent at
1135 the subdivision point. The tangent at that point must be the
1136 same for both curves, or noticeable bumps will occur in the
1137 character. But we want to use information on both sides of the
1138 point to compute the tangent, hence cross_curve = true. */
1139 find_tangent (left_curve, /* to_start_point: */ false,
1140 /* cross_curve: */ true, fitting_opts->tangent_surround);
1141 CURVE_START_TANGENT (right_curve) = CURVE_END_TANGENT (left_curve);
1142
1143 /* Now that we've set up the curves, we can fit them. */
1144 left_spline_list = fit_curve (left_curve, fitting_opts,
1145 exception);
1146 if (at_exception_got_fatal(exception))
1147 /* TODO: Memory allocated for left_curve and right_curve
1148 will leak.*/
1149 goto cleanup;
1150
1151 right_spline_list = fit_curve (right_curve, fitting_opts,
1152 exception);
1153 /* TODO: Memory allocated for left_curve and right_curve
1154 will leak.*/
1155 if (at_exception_got_fatal(exception))
1156 goto cleanup;
1157
1158 /* Neither of the subdivided curves could be fit, so fail. */
1159 if (left_spline_list == NULL && right_spline_list == NULL)
1160 return NULL;
1161
1162 /* Put the two together (or whichever of them exist). */
1163 spline_list = new_spline_list ();
1164
1165 if (left_spline_list == NULL)
1166 {
1167 LOG1 ("Could not fit spline to left curve (%lx).\n",
1168 (unsigned long) left_curve);
1169 at_exception_warning(exception, "Could not fit left spline list");
1170 }
1171 else
1172 {
1173 concat_spline_lists (spline_list, *left_spline_list);
1174 free_spline_list (*left_spline_list);
1175 free (left_spline_list);
1176 }
1177
1178 if (right_spline_list == NULL)
1179 {
1180 LOG1 ("Could not fit spline to right curve (%lx).\n",
1181 (unsigned long) right_curve);
1182 at_exception_warning(exception, "Could not fit right spline list");
1183 }
1184 else
1185 {
1186 concat_spline_lists (spline_list, *right_spline_list);
1187 free_spline_list (*right_spline_list);
1188 free (right_spline_list);
1189 }
1190 if (CURVE_END_TANGENT (left_curve))
1191 free (CURVE_END_TANGENT (left_curve));
1192 free (left_curve);
1193 free (right_curve);
1194 }
1195 cleanup:
1196 return spline_list;
1197 }
1198
1199 /* Our job here is to find alpha1 (and alpha2), where t1_hat (t2_hat) is
1200 the tangent to CURVE at the starting (ending) point, such that:
1201
1202 control1 = alpha1*t1_hat + starting point
1203 control2 = alpha2*t2_hat + ending_point
1204
1205 and the resulting spline (starting_point .. control1 and control2 ..
1206 ending_point) minimizes the least-square error from CURVE.
1207
1208 See pp.57--59 of the Phoenix thesis.
1209
1210 The B?(t) here corresponds to B_i^3(U_i) there.
1211 The Bernshte\u in polynomials of degree n are defined by
1212 B_i^n(t) = { n \choose i } t^i (1-t)^{n-i}, i = 0..n */
1213
1214 #define B0(t) CUBE ((at_real) 1.0 - (t))
1215 #define B1(t) ((at_real) 3.0 * (t) * SQUARE ((at_real) 1.0 - (t)))
1216 #define B2(t) ((at_real) 3.0 * SQUARE (t) * ((at_real) 1.0 - (t)))
1217 #define B3(t) CUBE (t)
1218
1219 static spline_type
fit_one_spline(curve_type curve,at_exception_type * exception)1220 fit_one_spline (curve_type curve,
1221 at_exception_type * exception)
1222 {
1223 /* Since our arrays are zero-based, the `C0' and `C1' here correspond
1224 to `C1' and `C2' in the paper. */
1225 at_real X_C1_det, C0_X_det, C0_C1_det;
1226 at_real alpha1, alpha2;
1227 spline_type spline;
1228 vector_type start_vector, end_vector;
1229 unsigned i;
1230 vector_type *A;
1231 vector_type t1_hat = *CURVE_START_TANGENT (curve);
1232 vector_type t2_hat = *CURVE_END_TANGENT (curve);
1233 at_real C[2][2] = { { 0.0, 0.0 }, { 0.0, 0.0 } };
1234 at_real X[2] = { 0.0, 0.0 };
1235
1236 XMALLOC (A, CURVE_LENGTH (curve) * 2
1237 * sizeof (vector_type )); /* A dynamically allocated array. */
1238
1239 START_POINT (spline) = CURVE_POINT (curve, 0);
1240 END_POINT (spline) = LAST_CURVE_POINT (curve);
1241 start_vector = make_vector (START_POINT (spline));
1242 end_vector = make_vector (END_POINT (spline));
1243
1244 for (i = 0; i < CURVE_LENGTH (curve); i++)
1245 {
1246 A[(i << 1) + 0] = Vmult_scalar (t1_hat, B1 (CURVE_T (curve, i)));
1247 A[(i << 1) + 1] = Vmult_scalar (t2_hat, B2 (CURVE_T (curve, i)));
1248 }
1249
1250 for (i = 0; i < CURVE_LENGTH (curve); i++)
1251 {
1252 vector_type temp, temp0, temp1, temp2, temp3;
1253 vector_type *Ai = A + (i << 1);
1254
1255 C[0][0] += Vdot (Ai[0], Ai[0]);
1256 C[0][1] += Vdot (Ai[0], Ai[1]);
1257 /* C[1][0] = C[0][1] (this is assigned outside the loop) */
1258 C[1][1] += Vdot (Ai[1], Ai[1]);
1259
1260
1261 /* Now the right-hand side of the equation in the paper. */
1262 temp0 = Vmult_scalar (start_vector, B0 (CURVE_T (curve, i)));
1263 temp1 = Vmult_scalar (start_vector, B1 (CURVE_T (curve, i)));
1264 temp2 = Vmult_scalar (end_vector, B2 (CURVE_T (curve, i)));
1265 temp3 = Vmult_scalar (end_vector, B3 (CURVE_T (curve, i)));
1266
1267 temp = make_vector (Vsubtract_point (CURVE_POINT (curve, i),
1268 Vadd (temp0, Vadd (temp1, Vadd (temp2, temp3)))));
1269
1270 X[0] += Vdot (temp, Ai[0]);
1271 X[1] += Vdot (temp, Ai[1]);
1272 }
1273 free (A);
1274
1275 C[1][0] = C[0][1];
1276
1277 X_C1_det = X[0] * C[1][1] - X[1] * C[0][1];
1278 C0_X_det = C[0][0] * X[1] - C[0][1] * X[0];
1279 C0_C1_det = C[0][0] * C[1][1] - C[1][0] * C[0][1];
1280 if (C0_C1_det == 0.0)
1281 {
1282 LOG ("zero determinant of C0*C1");
1283 at_exception_fatal(exception, "zero determinant of C0*C1");
1284 goto cleanup;
1285 }
1286
1287 alpha1 = X_C1_det / C0_C1_det;
1288 alpha2 = C0_X_det / C0_C1_det;
1289
1290 CONTROL1 (spline) = Vadd_point (START_POINT (spline),
1291 Vmult_scalar (t1_hat, alpha1));
1292 CONTROL2 (spline) = Vadd_point (END_POINT (spline),
1293 Vmult_scalar (t2_hat, alpha2));
1294 SPLINE_DEGREE (spline) = CUBICTYPE;
1295
1296 cleanup:
1297 return spline;
1298 }
1299
1300 /* Find reasonable values for t for each point on CURVE. The method is
1301 called chord-length parameterization, which is described in Plass &
1302 Stone. The basic idea is just to use the distance from one point to
1303 the next as the t value, normalized to produce values that increase
1304 from zero for the first point to one for the last point. */
1305
1306 static void
set_initial_parameter_values(curve_type curve)1307 set_initial_parameter_values (curve_type curve)
1308 {
1309 unsigned p;
1310
1311 LOG ("\nAssigning initial t values:\n ");
1312
1313 CURVE_T (curve, 0) = 0.0;
1314
1315 for (p = 1; p < CURVE_LENGTH (curve); p++)
1316 {
1317 at_real_coord point = CURVE_POINT (curve, p),
1318 previous_p = CURVE_POINT (curve, p - 1);
1319 at_real d = distance (point, previous_p);
1320 CURVE_T (curve, p) = CURVE_T (curve, p - 1) + d;
1321 }
1322
1323 assert (LAST_CURVE_T (curve) != 0.0);
1324
1325 for (p = 1; p < CURVE_LENGTH (curve); p++)
1326 CURVE_T (curve, p) = CURVE_T (curve, p) / LAST_CURVE_T (curve);
1327
1328 log_entire_curve (curve);
1329 }
1330
1331 /* Find an approximation to the tangent to an endpoint of CURVE (to the
1332 first point if TO_START_POINT is true, else the last). If
1333 CROSS_CURVE is true, consider points on the adjacent curve to CURVE.
1334
1335 It is important to compute an accurate approximation, because the
1336 control points that we eventually decide upon to fit the curve will
1337 be placed on the half-lines defined by the tangents and
1338 endpoints...and we never recompute the tangent after this. */
1339
1340 static void
find_tangent(curve_type curve,at_bool to_start_point,at_bool cross_curve,unsigned tangent_surround)1341 find_tangent (curve_type curve, at_bool to_start_point, at_bool cross_curve,
1342 unsigned tangent_surround)
1343 {
1344 vector_type tangent;
1345 vector_type **curve_tangent = (to_start_point == true) ? &(CURVE_START_TANGENT (curve))
1346 : &(CURVE_END_TANGENT (curve));
1347 unsigned n_points = 0;
1348
1349 LOG1 (" tangent to %s: ", (to_start_point == true) ? "start" : "end");
1350
1351 if (*curve_tangent == NULL)
1352 {
1353 XMALLOC (*curve_tangent, sizeof (vector_type));
1354 do
1355 {
1356 tangent = find_half_tangent (curve, to_start_point, &n_points,
1357 tangent_surround);
1358
1359 if ((cross_curve == true) || (CURVE_CYCLIC (curve) == true))
1360 {
1361 curve_type adjacent_curve
1362 = (to_start_point == true) ? PREVIOUS_CURVE (curve) : NEXT_CURVE (curve);
1363 vector_type tangent2
1364 = (to_start_point == false) ? find_half_tangent (adjacent_curve, true, &n_points,
1365 tangent_surround) : find_half_tangent (adjacent_curve, true, &n_points,
1366 tangent_surround);
1367
1368 LOG3 ("(adjacent curve half tangent (%.3f,%.3f,%.3f)) ",
1369 tangent2.dx, tangent2.dy, tangent2.dz);
1370 tangent = Vadd (tangent, tangent2);
1371 }
1372 tangent_surround--;
1373
1374 }
1375 while (tangent.dx == 0.0 && tangent.dy == 0.0);
1376
1377 assert (n_points > 0);
1378 **curve_tangent = Vmult_scalar (tangent, (at_real)(1.0 / n_points));
1379 if ((CURVE_CYCLIC (curve) == true) && CURVE_START_TANGENT (curve))
1380 *CURVE_START_TANGENT (curve) = **curve_tangent;
1381 if ((CURVE_CYCLIC (curve) == true) && CURVE_END_TANGENT (curve))
1382 *CURVE_END_TANGENT (curve) = **curve_tangent;
1383 }
1384 else
1385 LOG ("(already computed) ");
1386
1387 LOG3 ("(%.3f,%.3f,%.3f).\n", (*curve_tangent)->dx, (*curve_tangent)->dy, (*curve_tangent)->dz);
1388 }
1389
1390 /* Find the change in y and change in x for `tangent_surround' (a global)
1391 points along CURVE. Increment N_POINTS by the number of points we
1392 actually look at. */
1393
1394 static vector_type
find_half_tangent(curve_type c,at_bool to_start_point,unsigned * n_points,unsigned tangent_surround)1395 find_half_tangent (curve_type c, at_bool to_start_point, unsigned *n_points,
1396 unsigned tangent_surround)
1397 {
1398 unsigned p;
1399 int factor = to_start_point ? 1 : -1;
1400 unsigned tangent_index = to_start_point ? 0 : c->length - 1;
1401 at_real_coord tangent_point = CURVE_POINT (c, tangent_index);
1402 vector_type tangent = { 0.0, 0.0 };
1403 unsigned int surround;
1404
1405 if ((surround = CURVE_LENGTH (c) / 2) > tangent_surround)
1406 surround = tangent_surround;
1407
1408 for (p = 1; p <= surround; p++)
1409 {
1410 int this_index = p * factor + tangent_index;
1411 at_real_coord this_point;
1412
1413 if (this_index < 0 || this_index >= (int) c->length)
1414 break;
1415
1416 this_point = CURVE_POINT (c, p * factor + tangent_index);
1417
1418 /* Perhaps we should weight the tangent from `this_point' by some
1419 factor dependent on the distance from the tangent point. */
1420 tangent = Vadd (tangent,
1421 Vmult_scalar (Psubtract (this_point, tangent_point),
1422 (at_real) factor));
1423 (*n_points)++;
1424 }
1425
1426 return tangent;
1427 }
1428
1429 /* When this routine is called, we have computed a spline representation
1430 for the digitized curve. The question is, how good is it? If the
1431 fit is very good indeed, we might have an error of zero on each
1432 point, and then WORST_POINT becomes irrelevant. But normally, we
1433 return the error at the worst point, and the index of that point in
1434 WORST_POINT. The error computation itself is the Euclidean distance
1435 from the original curve CURVE to the fitted spline SPLINE. */
1436
1437 static at_real
find_error(curve_type curve,spline_type spline,unsigned * worst_point,at_exception_type * exception)1438 find_error (curve_type curve, spline_type spline, unsigned *worst_point,
1439 at_exception_type * exception)
1440 {
1441 unsigned this_point;
1442 at_real total_error = 0.0;
1443 at_real worst_error = FLT_MIN;
1444
1445 *worst_point = CURVE_LENGTH (curve) + 1; /* A sentinel value. */
1446
1447 for (this_point = 0; this_point < CURVE_LENGTH (curve); this_point++)
1448 {
1449 at_real_coord curve_point = CURVE_POINT (curve, this_point);
1450 at_real t = CURVE_T (curve, this_point);
1451 at_real_coord spline_point = evaluate_spline (spline, t);
1452 at_real this_error = distance (curve_point, spline_point);
1453 if (this_error >= worst_error)
1454 {
1455 *worst_point = this_point;
1456 worst_error = this_error;
1457 }
1458 total_error += this_error;
1459 }
1460
1461 if (*worst_point == CURVE_LENGTH (curve) + 1)
1462 { /* Didn't have any ``worst point''; the error should be zero. */
1463 if (epsilon_equal (total_error, 0.0))
1464 LOG (" Every point fit perfectly.\n");
1465 else
1466 {
1467 LOG("No worst point found; something is wrong");
1468 at_exception_warning(exception, "No worst point found; something is wrong");
1469 }
1470 }
1471 else
1472 {
1473 if (epsilon_equal (total_error, 0.0))
1474 LOG (" Every point fit perfectly.\n");
1475 else
1476 {
1477 LOG5 (" Worst error (at (%.3f,%.3f,%.3f), point #%u) was %.3f.\n",
1478 CURVE_POINT (curve, *worst_point).x,
1479 CURVE_POINT (curve, *worst_point).y,
1480 CURVE_POINT (curve, *worst_point).z, *worst_point, worst_error);
1481 LOG1 (" Total error was %.3f.\n", total_error);
1482 LOG2 (" Average error (over %u points) was %.3f.\n",
1483 CURVE_LENGTH (curve), total_error / CURVE_LENGTH (curve));
1484 }
1485 }
1486
1487 return worst_error;
1488 }
1489
1490 /* Supposing that we have accepted the error, another question arises:
1491 would we be better off just using a straight line? */
1492
1493 static at_bool
spline_linear_enough(spline_type * spline,curve_type curve,fitting_opts_type * fitting_opts)1494 spline_linear_enough (spline_type *spline, curve_type curve,
1495 fitting_opts_type *fitting_opts)
1496 {
1497 at_real A, B, C;
1498 unsigned this_point;
1499 at_real dist = 0.0, start_end_dist, threshold;
1500
1501 LOG ("Checking linearity:\n");
1502
1503 A = END_POINT(*spline).x - START_POINT(*spline).x;
1504 B = END_POINT(*spline).y - START_POINT(*spline).y;
1505 C = END_POINT(*spline).z - START_POINT(*spline).z;
1506
1507 start_end_dist = (at_real) (SQUARE(A) + SQUARE(B) + SQUARE(C));
1508 LOG1 ("start_end_distance is %.3f.\n", sqrt(start_end_dist));
1509
1510 LOG3 (" Line endpoints are (%.3f, %.3f, %.3f) and ", START_POINT(*spline).x, START_POINT(*spline).y, START_POINT(*spline).z);
1511 LOG3 ("(%.3f, %.3f, %.3f)\n", END_POINT(*spline).x, END_POINT(*spline).y, END_POINT(*spline).z);
1512
1513 /* LOG3 (" Line is %.3fx + %.3fy + %.3f = 0.\n", A, B, C); */
1514
1515 for (this_point = 0; this_point < CURVE_LENGTH (curve); this_point++)
1516 {
1517 at_real a, b, c, w;
1518 at_real t = CURVE_T (curve, this_point);
1519 at_real_coord spline_point = evaluate_spline (*spline, t);
1520
1521 a = spline_point.x - START_POINT(*spline).x;
1522 b = spline_point.y - START_POINT(*spline).y;
1523 c = spline_point.z - START_POINT(*spline).z;
1524 w = (A*a + B*b + C*c) / start_end_dist;
1525
1526 dist += (at_real)sqrt(SQUARE(a-A*w) + SQUARE(b-B*w) + SQUARE(c-C*w));
1527 }
1528 LOG1 (" Total distance is %.3f, ", dist);
1529
1530 dist /= (CURVE_LENGTH (curve) - 1);
1531 LOG1 ("which is %.3f normalized.\n", dist);
1532
1533 /* We want reversion of short curves to splines to be more likely than
1534 reversion of long curves, hence the second division by the curve
1535 length, for use in `change_bad_lines'. */
1536 SPLINE_LINEARITY (*spline) = dist;
1537 LOG1 (" Final linearity: %.3f.\n", SPLINE_LINEARITY (*spline));
1538 if (start_end_dist * (at_real) 0.5 > fitting_opts->line_threshold)
1539 threshold = fitting_opts->line_threshold;
1540 else
1541 threshold = start_end_dist * (at_real) 0.5;
1542 LOG1 ("threshold is %.3f .\n", threshold);
1543 if (dist < threshold)
1544 return true;
1545 else
1546 return false;
1547 }
1548
1549 /* Unfortunately, we cannot tell in isolation whether a given spline
1550 should be changed to a line or not. That can only be known after the
1551 entire curve has been fit to a list of splines. (The curve is the
1552 pixel outline between two corners.) After subdividing the curve, a
1553 line may very well fit a portion of the curve just as well as the
1554 spline---but unless a spline is truly close to being a line, it
1555 should not be combined with other lines. */
1556
1557 static void
change_bad_lines(spline_list_type * spline_list,fitting_opts_type * fitting_opts)1558 change_bad_lines (spline_list_type *spline_list,
1559 fitting_opts_type *fitting_opts)
1560 {
1561 unsigned this_spline;
1562 at_bool found_cubic = false;
1563 unsigned length = SPLINE_LIST_LENGTH (*spline_list);
1564
1565 LOG1 ("\nChecking for bad lines (length %u):\n", length);
1566
1567 /* First see if there are any splines in the fitted shape. */
1568 for (this_spline = 0; this_spline < length; this_spline++)
1569 {
1570 if (SPLINE_DEGREE (SPLINE_LIST_ELT (*spline_list, this_spline)) ==
1571 CUBICTYPE)
1572 {
1573 found_cubic = true;
1574 break;
1575 }
1576 }
1577
1578 /* If so, change lines back to splines (we haven't done anything to
1579 their control points, so we only have to change the degree) unless
1580 the spline is close enough to being a line. */
1581 if (found_cubic)
1582 for (this_spline = 0; this_spline < length; this_spline++)
1583 {
1584 spline_type s = SPLINE_LIST_ELT (*spline_list, this_spline);
1585
1586 if (SPLINE_DEGREE (s) == LINEARTYPE)
1587 {
1588 LOG1 (" #%u: ", this_spline);
1589 if (SPLINE_LINEARITY (s) > fitting_opts->line_reversion_threshold)
1590 {
1591 LOG ("reverted, ");
1592 SPLINE_DEGREE (SPLINE_LIST_ELT (*spline_list, this_spline))
1593 = CUBICTYPE;
1594 }
1595 LOG1 ("linearity %.3f.\n", SPLINE_LINEARITY (s));
1596 }
1597 }
1598 else
1599 LOG (" No lines.\n");
1600 }
1601
1602 /* Lists of array indices (well, that is what we use it for). */
1603
1604 static index_list_type
new_index_list(void)1605 new_index_list (void)
1606 {
1607 index_list_type index_list;
1608
1609 index_list.data = NULL;
1610 INDEX_LIST_LENGTH (index_list) = 0;
1611
1612 return index_list;
1613 }
1614
1615 static void
free_index_list(index_list_type * index_list)1616 free_index_list (index_list_type *index_list)
1617 {
1618 if (INDEX_LIST_LENGTH (*index_list) > 0)
1619 {
1620 free (index_list->data);
1621 index_list->data = NULL;
1622 INDEX_LIST_LENGTH (*index_list) = 0;
1623 }
1624 }
1625
1626 static void
append_index(index_list_type * list,unsigned new_index)1627 append_index (index_list_type *list, unsigned new_index)
1628 {
1629 INDEX_LIST_LENGTH (*list)++;
1630 XREALLOC (list->data, INDEX_LIST_LENGTH (*list) * sizeof (unsigned));
1631 list->data[INDEX_LIST_LENGTH (*list) - 1] = new_index;
1632 }
1633
1634 /* Turn an real point into a integer one. */
1635
1636 static at_coord
real_to_int_coord(at_real_coord real_coord)1637 real_to_int_coord (at_real_coord real_coord)
1638 {
1639 at_coord int_coord;
1640
1641 int_coord.x = ROUND (real_coord.x);
1642 int_coord.y = ROUND (real_coord.y);
1643
1644 return int_coord;
1645 }
1646
1647 /* Return the Euclidean distance between P1 and P2. */
1648
1649 static at_real
distance(at_real_coord p1,at_real_coord p2)1650 distance (at_real_coord p1, at_real_coord p2)
1651 {
1652 at_real x = p1.x - p2.x, y = p1.y - p2.y, z = p1.z - p2.z;
1653 return (at_real) sqrt (SQUARE(x)
1654 + SQUARE(y) + SQUARE(z));
1655 }
1656