summaryrefslogtreecommitdiffstats
path: root/src/3rdparty/autotrace/fit.c
diff options
context:
space:
mode:
Diffstat (limited to '')
-rw-r--r--src/3rdparty/autotrace/fit.c1442
1 files changed, 1442 insertions, 0 deletions
diff --git a/src/3rdparty/autotrace/fit.c b/src/3rdparty/autotrace/fit.c
new file mode 100644
index 0000000..a6ca2b6
--- /dev/null
+++ b/src/3rdparty/autotrace/fit.c
@@ -0,0 +1,1442 @@
+/* fit.c: turn a bitmap representation of a curve into a list of splines.
+ Some of the ideas, but not the code, comes from the Phoenix thesis.
+ See README for the reference.
+
+ The code was partially derived from limn.
+
+ Copyright (C) 1992 Free Software Foundation, Inc.
+
+ This program is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2, or (at your option)
+ any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program; if not, write to the Free Software
+ Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif /* Def: HAVE_CONFIG_H */
+
+#include "autotrace.h"
+#include "fit.h"
+#include "logreport.h"
+#include "spline.h"
+#include "vector.h"
+#include "curve.h"
+#include "pxl-outline.h"
+#include "epsilon-equal.h"
+#include "xstd.h"
+#include <math.h>
+#ifndef FLT_MAX
+#include <limits.h>
+#include <float.h>
+#endif
+#ifndef FLT_MIN
+#include <limits.h>
+#include <float.h>
+#endif
+#include <string.h>
+#include <assert.h>
+
+#define SQUARE(x) ((x) * (x))
+#define CUBE(x) ((x) * (x) * (x))
+
+/* We need to manipulate lists of array indices. */
+
+typedef struct index_list {
+ unsigned *data;
+ unsigned length;
+} index_list_type;
+
+/* The usual accessor macros. */
+#define GET_INDEX(i_l, n) ((i_l).data[n])
+#define INDEX_LIST_LENGTH(i_l) ((i_l).length)
+#define GET_LAST_INDEX(i_l) ((i_l).data[INDEX_LIST_LENGTH (i_l) - 1])
+
+static void append_index(index_list_type *, unsigned);
+static void free_index_list(index_list_type *);
+static index_list_type new_index_list(void);
+static void remove_adjacent_corners(index_list_type *, unsigned, gboolean, at_exception_type * exception);
+static void change_bad_lines(spline_list_type *, fitting_opts_type *);
+static void filter(curve_type, fitting_opts_type *);
+static void find_vectors(unsigned, pixel_outline_type, vector_type *, vector_type *, unsigned);
+static index_list_type find_corners(pixel_outline_type, fitting_opts_type *, at_exception_type * exception);
+static gfloat find_error(curve_type, spline_type, unsigned *, at_exception_type * exception);
+static vector_type find_half_tangent(curve_type, gboolean start, unsigned *, unsigned);
+static void find_tangent(curve_type, gboolean, gboolean, unsigned);
+static spline_type fit_one_spline(curve_type, at_exception_type * exception);
+static spline_list_type *fit_curve(curve_type, fitting_opts_type *, at_exception_type * exception);
+static spline_list_type fit_curve_list(curve_list_type, fitting_opts_type *, at_distance_map *, at_exception_type * exception);
+static spline_list_type *fit_with_least_squares(curve_type, fitting_opts_type *, at_exception_type * exception);
+static spline_list_type *fit_with_line(curve_type);
+static void remove_knee_points(curve_type, gboolean);
+static void set_initial_parameter_values(curve_type);
+static gboolean spline_linear_enough(spline_type *, curve_type, fitting_opts_type *);
+static curve_list_array_type split_at_corners(pixel_outline_list_type, fitting_opts_type *, at_exception_type * exception);
+static at_coord real_to_int_coord(at_real_coord);
+static gfloat distance(at_real_coord, at_real_coord);
+
+/* Get a new set of fitting options */
+fitting_opts_type new_fitting_opts(void)
+{
+ fitting_opts_type fitting_opts;
+
+ fitting_opts.background_color = NULL;
+ fitting_opts.charcode = 0;
+ fitting_opts.color_count = 0;
+ fitting_opts.corner_always_threshold = (gfloat) 60.0;
+ fitting_opts.corner_surround = 4;
+ fitting_opts.corner_threshold = (gfloat) 100.0;
+ fitting_opts.error_threshold = (gfloat) 2.0;
+ fitting_opts.filter_iterations = 4;
+ fitting_opts.line_reversion_threshold = (gfloat) .01;
+ fitting_opts.line_threshold = (gfloat) 1.0;
+ fitting_opts.remove_adjacent_corners = FALSE;
+ fitting_opts.tangent_surround = 3;
+ fitting_opts.despeckle_level = 0;
+ fitting_opts.despeckle_tightness = 2.0;
+ fitting_opts.noise_removal = (gfloat) 0.99;
+ fitting_opts.centerline = FALSE;
+ fitting_opts.preserve_width = FALSE;
+ fitting_opts.width_weight_factor = 6.0;
+
+ return (fitting_opts);
+}
+
+/* The top-level call that transforms the list of pixels in the outlines
+ of the original character to a list of spline lists fitted to those
+ pixels. */
+
+spline_list_array_type fitted_splines(pixel_outline_list_type pixel_outline_list, fitting_opts_type * fitting_opts, at_distance_map * dist, unsigned short width, unsigned short height, at_exception_type * exception, at_progress_func notify_progress, gpointer progress_data, at_testcancel_func test_cancel, gpointer testcancel_data)
+{
+ unsigned this_list;
+
+ spline_list_array_type char_splines = new_spline_list_array();
+ curve_list_array_type curve_array = split_at_corners(pixel_outline_list,
+ fitting_opts,
+ exception);
+
+ char_splines.centerline = fitting_opts->centerline;
+ char_splines.preserve_width = fitting_opts->preserve_width;
+ char_splines.width_weight_factor = fitting_opts->width_weight_factor;
+
+ if (fitting_opts->background_color)
+ char_splines.background_color = at_color_copy(fitting_opts->background_color);
+ else
+ char_splines.background_color = NULL;
+ /* Set dummy values. Real value is set in upper context. */
+ char_splines.width = width;
+ char_splines.height = height;
+
+ for (this_list = 0; this_list < CURVE_LIST_ARRAY_LENGTH(curve_array); this_list++) {
+ spline_list_type curve_list_splines;
+ curve_list_type curves = CURVE_LIST_ARRAY_ELT(curve_array, this_list);
+
+ if (notify_progress)
+ notify_progress((((gfloat) this_list) / ((gfloat) CURVE_LIST_ARRAY_LENGTH(curve_array) * (gfloat) 3.0) + (gfloat) 0.333), progress_data);
+ if (test_cancel && test_cancel(testcancel_data))
+ goto cleanup;
+
+ LOG("\nFitting curve list #%u:\n", this_list);
+
+ curve_list_splines = fit_curve_list(curves, fitting_opts, dist, exception);
+ if (at_exception_got_fatal(exception)) {
+ if (char_splines.background_color)
+ at_color_free(char_splines.background_color);
+ goto cleanup;
+ }
+ curve_list_splines.clockwise = curves.clockwise;
+
+ memcpy(&(curve_list_splines.color), &(O_LIST_OUTLINE(pixel_outline_list, this_list).color), sizeof(at_color));
+ append_spline_list(&char_splines, curve_list_splines);
+ }
+cleanup:
+ free_curve_list_array(&curve_array, notify_progress, progress_data);
+
+ return char_splines;
+}
+
+/* Fit the list of curves CURVE_LIST to a list of splines, and return
+ it. CURVE_LIST represents a single closed paths, e.g., either the
+ inside or outside outline of an `o'. */
+
+static spline_list_type fit_curve_list(curve_list_type curve_list, fitting_opts_type * fitting_opts, at_distance_map * dist, at_exception_type * exception)
+{
+ curve_type curve;
+ unsigned this_curve, this_spline;
+ unsigned curve_list_length = CURVE_LIST_LENGTH(curve_list);
+ spline_list_type curve_list_splines = empty_spline_list();
+
+ curve_list_splines.open = curve_list.open;
+
+ /* Remove the extraneous ``knee'' points before filtering. Since the
+ corners have already been found, we don't need to worry about
+ removing a point that should be a corner. */
+
+ LOG("\nRemoving knees:\n");
+ for (this_curve = 0; this_curve < curve_list_length; this_curve++) {
+ LOG("#%u:", this_curve);
+ remove_knee_points(CURVE_LIST_ELT(curve_list, this_curve), CURVE_LIST_CLOCKWISE(curve_list));
+ }
+
+ if (dist != NULL) {
+ unsigned this_point;
+ unsigned height = dist->height;
+ for (this_curve = 0; this_curve < curve_list_length; this_curve++) {
+ curve = CURVE_LIST_ELT(curve_list, this_curve);
+ for (this_point = 0; this_point < CURVE_LENGTH(curve); this_point++) {
+ unsigned x, y;
+ float width, w;
+ at_real_coord *coord = &CURVE_POINT(curve, this_point);
+ x = (unsigned)(coord->x);
+ y = height - (unsigned)(coord->y) - 1;
+
+ /* Each (x, y) is a point on the skeleton of the curve, which
+ might be offset from the TRUE centerline, where the width
+ is maximal. Therefore, use as the local line width the
+ maximum distance over the neighborhood of (x, y). */
+ width = dist->d[y][x];
+ if (y >= 1) {
+ if ((w = dist->d[y - 1][x]) > width)
+ width = w;
+ if (x >= 1) {
+ if ((w = dist->d[y][x - 1]) > width)
+ width = w;
+ if ((w = dist->d[y - 1][x - 1]) > width)
+ width = w;
+ }
+ if (x + 1 < dist->width) {
+ if ((w = dist->d[y][x + 1]) > width)
+ width = w;
+ if ((w = dist->d[y - 1][x + 1]) > width)
+ width = w;
+ }
+ }
+ if (y + 1 < height) {
+ if ((w = dist->d[y + 1][x]) > width)
+ width = w;
+ if (x >= 1 && (w = dist->d[y + 1][x - 1]) > width)
+ width = w;
+ if (x + 1 < dist->width && (w = dist->d[y + 1][x + 1]) > width)
+ width = w;
+ }
+ coord->z = width * (fitting_opts->width_weight_factor);
+ }
+ }
+ }
+
+ /* We filter all the curves in CURVE_LIST at once; otherwise, we would
+ look at an unfiltered curve when computing tangents. */
+
+ LOG("\nFiltering curves:\n");
+ for (this_curve = 0; this_curve < curve_list.length; this_curve++) {
+ LOG("#%u: ", this_curve);
+ filter(CURVE_LIST_ELT(curve_list, this_curve), fitting_opts);
+ }
+
+ /* Make the first point in the first curve also be the last point in
+ the last curve, so the fit to the whole curve list will begin and
+ end at the same point. This may cause slight errors in computing
+ the tangents and t values, but it's worth it for the continuity.
+ Of course we don't want to do this if the two points are already
+ the same, as they are if the curve is cyclic. (We don't append it
+ earlier, in `split_at_corners', because that confuses the
+ filtering.) Finally, we can't append the point if the curve is
+ exactly three points long, because we aren't adding any more data,
+ and three points isn't enough to determine a spline. Therefore,
+ the fitting will fail. */
+ curve = CURVE_LIST_ELT(curve_list, 0);
+ if (CURVE_CYCLIC(curve) == TRUE)
+ append_point(curve, CURVE_POINT(curve, 0));
+
+ /* Finally, fit each curve in the list to a list of splines. */
+ for (this_curve = 0; this_curve < curve_list_length; this_curve++) {
+ spline_list_type *curve_splines;
+ curve_type current_curve = CURVE_LIST_ELT(curve_list, this_curve);
+
+ LOG("\nFitting curve #%u:\n", this_curve);
+
+ curve_splines = fit_curve(current_curve, fitting_opts, exception);
+ if (at_exception_got_fatal(exception))
+ goto cleanup;
+ else if (curve_splines == NULL) {
+ LOG("Could not fit curve #%u", this_curve);
+ at_exception_warning(exception, "Could not fit curve");
+ } else {
+ LOG("Fitted splines for curve #%u:\n", this_curve);
+ for (this_spline = 0; this_spline < SPLINE_LIST_LENGTH(*curve_splines); this_spline++) {
+ LOG(" %u: ", this_spline);
+ if (logging)
+ print_spline(SPLINE_LIST_ELT(*curve_splines, this_spline));
+ }
+
+ /* After fitting, we may need to change some would-be lines
+ back to curves, because they are in a list with other
+ curves. */
+ change_bad_lines(curve_splines, fitting_opts);
+
+ concat_spline_lists(&curve_list_splines, *curve_splines);
+ free_spline_list(*curve_splines);
+ free(curve_splines);
+ }
+ }
+
+ if (logging) {
+ LOG("\nFitted splines are:\n");
+ for (this_spline = 0; this_spline < SPLINE_LIST_LENGTH(curve_list_splines); this_spline++) {
+ LOG(" %u: ", this_spline);
+ print_spline(SPLINE_LIST_ELT(curve_list_splines, this_spline));
+ }
+ }
+cleanup:
+ return curve_list_splines;
+}
+
+/* Transform a set of locations to a list of splines (the fewer the
+ better). We are guaranteed that CURVE does not contain any corners.
+ We return NULL if we cannot fit the points at all. */
+
+static spline_list_type *fit_curve(curve_type curve, fitting_opts_type * fitting_opts, at_exception_type * exception)
+{
+ spline_list_type *fittedsplines;
+
+ if (CURVE_LENGTH(curve) < 2) {
+ LOG("Tried to fit curve with less than two points");
+ at_exception_warning(exception, "Tried to fit curve with less than two points");
+ return NULL;
+ }
+
+ /* Do we have enough points to fit with a spline? */
+ fittedsplines = CURVE_LENGTH(curve) < 4 ? fit_with_line(curve)
+ : fit_with_least_squares(curve, fitting_opts, exception);
+
+ return fittedsplines;
+}
+
+/* As mentioned above, the first step is to find the corners in
+ PIXEL_LIST, the list of points. (Presumably we can't fit a single
+ spline around a corner.) The general strategy is to look through all
+ the points, remembering which we want to consider corners. Then go
+ through that list, producing the curve_list. This is dictated by the
+ fact that PIXEL_LIST does not necessarily start on a corner---it just
+ starts at the character's first outline pixel, going left-to-right,
+ top-to-bottom. But we want all our splines to start and end on real
+ corners.
+
+ For example, consider the top of a capital `C' (this is in cmss20):
+ x
+ ***********
+ ******************
+
+ PIXEL_LIST will start at the pixel below the `x'. If we considered
+ this pixel a corner, we would wind up matching a very small segment
+ from there to the end of the line, probably as a straight line, which
+ is certainly not what we want.
+
+ PIXEL_LIST has one element for each closed outline on the character.
+ To preserve this information, we return an array of curve_lists, one
+ element (which in turn consists of several curves, one between each
+ pair of corners) for each element in PIXEL_LIST. */
+
+static curve_list_array_type split_at_corners(pixel_outline_list_type pixel_list, fitting_opts_type * fitting_opts, at_exception_type * exception)
+{
+ unsigned this_pixel_o;
+ curve_list_array_type curve_array = new_curve_list_array();
+
+ LOG("\nFinding corners:\n");
+
+ for (this_pixel_o = 0; this_pixel_o < O_LIST_LENGTH(pixel_list); this_pixel_o++) {
+ curve_type curve, first_curve;
+ index_list_type corner_list;
+ unsigned p, this_corner;
+ curve_list_type curve_list = new_curve_list();
+ pixel_outline_type pixel_o = O_LIST_OUTLINE(pixel_list, this_pixel_o);
+
+ CURVE_LIST_CLOCKWISE(curve_list) = O_CLOCKWISE(pixel_o);
+ curve_list.open = pixel_o.open;
+
+ LOG("#%u:", this_pixel_o);
+
+ /* If the outline does not have enough points, we can't do
+ anything. The endpoints of the outlines are automatically
+ corners. We need at least `corner_surround' more pixels on
+ either side of a point before it is conceivable that we might
+ want another corner. */
+ if (O_LENGTH(pixel_o) > fitting_opts->corner_surround * 2 + 2)
+ corner_list = find_corners(pixel_o, fitting_opts, exception);
+
+ else {
+ int surround;
+ if ((surround = (int)(O_LENGTH(pixel_o) - 3) / 2) >= 2) {
+ unsigned save_corner_surround = fitting_opts->corner_surround;
+ fitting_opts->corner_surround = surround;
+ corner_list = find_corners(pixel_o, fitting_opts, exception);
+ fitting_opts->corner_surround = save_corner_surround;
+ } else {
+ corner_list.length = 0;
+ corner_list.data = NULL;
+ }
+ }
+
+ /* Remember the first curve so we can make it be the `next' of the
+ last one. (And vice versa.) */
+ first_curve = new_curve();
+
+ curve = first_curve;
+
+ if (corner_list.length == 0) { /* No corners. Use all of the pixel outline as the curve. */
+ for (p = 0; p < O_LENGTH(pixel_o); p++)
+ append_pixel(curve, O_COORDINATE(pixel_o, p));
+
+ if (curve_list.open == TRUE)
+ CURVE_CYCLIC(curve) = FALSE;
+ else
+ CURVE_CYCLIC(curve) = TRUE;
+ } else { /* Each curve consists of the points between (inclusive) each pair
+ of corners. */
+ for (this_corner = 0; this_corner < corner_list.length - 1; this_corner++) {
+ curve_type previous_curve = curve;
+ unsigned corner = GET_INDEX(corner_list, this_corner);
+ unsigned next_corner = GET_INDEX(corner_list, this_corner + 1);
+
+ for (p = corner; p <= next_corner; p++)
+ append_pixel(curve, O_COORDINATE(pixel_o, p));
+
+ append_curve(&curve_list, curve);
+ curve = new_curve();
+ NEXT_CURVE(previous_curve) = curve;
+ PREVIOUS_CURVE(curve) = previous_curve;
+ }
+
+ /* The last curve is different. It consists of the points
+ (inclusive) between the last corner and the end of the list,
+ and the beginning of the list and the first corner. */
+ for (p = GET_LAST_INDEX(corner_list); p < O_LENGTH(pixel_o); p++)
+ append_pixel(curve, O_COORDINATE(pixel_o, p));
+
+ if (!pixel_o.open) {
+ for (p = 0; p <= GET_INDEX(corner_list, 0); p++)
+ append_pixel(curve, O_COORDINATE(pixel_o, p));
+ } else {
+ curve_type last_curve = PREVIOUS_CURVE(curve);
+ PREVIOUS_CURVE(first_curve) = NULL;
+ if (last_curve)
+ NEXT_CURVE(last_curve) = NULL;
+ }
+ }
+
+ LOG(" [%u].\n", corner_list.length);
+ free_index_list(&corner_list);
+
+ /* Add `curve' to the end of the list, updating the pointers in
+ the chain. */
+ append_curve(&curve_list, curve);
+ NEXT_CURVE(curve) = first_curve;
+ PREVIOUS_CURVE(first_curve) = curve;
+
+ /* And now add the just-completed curve list to the array. */
+ append_curve_list(&curve_array, curve_list);
+ } /* End of considering each pixel outline. */
+
+ return curve_array;
+}
+
+/* We consider a point to be a corner if (1) the angle defined by the
+ `corner_surround' points coming into it and going out from it is less
+ than `corner_threshold' degrees, and no point within
+ `corner_surround' points has a smaller angle; or (2) the angle is less
+ than `corner_always_threshold' degrees.
+
+ Because of the different cases, it is convenient to have the
+ following macro to append a corner on to the list we return. The
+ character argument C is simply so that the different cases can be
+ distinguished in the log file. */
+
+#define APPEND_CORNER(index, angle, c) \
+ do \
+ { \
+ append_index (&corner_list, index); \
+ LOG (" (%d,%d)%c%.3f", \
+ O_COORDINATE (pixel_outline, index).x, \
+ O_COORDINATE (pixel_outline, index).y, \
+ c, angle); \
+ } \
+ while (0)
+
+static index_list_type find_corners(pixel_outline_type pixel_outline, fitting_opts_type * fitting_opts, at_exception_type * exception)
+{
+ unsigned p, start_p, end_p;
+ index_list_type corner_list = new_index_list();
+
+ start_p = 0;
+ end_p = O_LENGTH(pixel_outline) - 1;
+ if (pixel_outline.open) {
+ if (end_p <= fitting_opts->corner_surround * 2)
+ return corner_list;
+ APPEND_CORNER(0, 0.0, '@');
+ start_p += fitting_opts->corner_surround;
+ end_p -= fitting_opts->corner_surround;
+ }
+
+ /* Consider each pixel on the outline in turn. */
+ for (p = start_p; p <= end_p; p++) {
+ gfloat corner_angle;
+ vector_type in_vector, out_vector;
+
+ /* Check if the angle is small enough. */
+ find_vectors(p, pixel_outline, &in_vector, &out_vector, fitting_opts->corner_surround);
+ corner_angle = Vangle(in_vector, out_vector, exception);
+ if (at_exception_got_fatal(exception))
+ goto cleanup;
+
+ if (fabs(corner_angle) <= fitting_opts->corner_threshold) {
+ /* We want to keep looking, instead of just appending the
+ first pixel we find with a small enough angle, since there
+ might be another corner within `corner_surround' pixels, with
+ a smaller angle. If that is the case, we want that one. */
+ gfloat best_corner_angle = corner_angle;
+ unsigned best_corner_index = p;
+ index_list_type equally_good_list = new_index_list();
+ /* As we come into the loop, `p' is the index of the point
+ that has an angle less than `corner_angle'. We use `i' to
+ move through the pixels next to that, and `q' for moving
+ through the adjacent pixels to each `p'. */
+ unsigned q = p;
+ unsigned i = p + 1;
+
+ while (TRUE) {
+ /* Perhaps the angle is sufficiently small that we want to
+ consider this a corner, even if it's not the best
+ (unless we've already wrapped around in the search,
+ i.e., `q<i', in which case we have already added the
+ corner, and we don't want to add it again). We want to
+ do this check on the first candidate we find, as well
+ as the others in the loop, hence this comes before the
+ stopping condition. */
+ if (corner_angle <= fitting_opts->corner_always_threshold && q >= p)
+ APPEND_CORNER(q, corner_angle, '\\');
+
+ /* Exit the loop if we've looked at `corner_surround'
+ pixels past the best one we found, or if we've looked
+ at all the pixels. */
+ if (i >= best_corner_index + fitting_opts->corner_surround || i >= O_LENGTH(pixel_outline))
+ break;
+
+ /* Check the angle. */
+ q = i % O_LENGTH(pixel_outline);
+ find_vectors(q, pixel_outline, &in_vector, &out_vector, fitting_opts->corner_surround);
+ corner_angle = Vangle(in_vector, out_vector, exception);
+ if (at_exception_got_fatal(exception))
+ goto cleanup;
+
+ /* If we come across a corner that is just as good as the
+ best one, we should make it a corner, too. This
+ happens, for example, at the points on the `W' in some
+ typefaces, where the ``points'' are flat. */
+ if (epsilon_equal(corner_angle, best_corner_angle))
+ append_index(&equally_good_list, q);
+
+ else if (corner_angle < best_corner_angle) {
+ best_corner_angle = corner_angle;
+ /* We want to check `corner_surround' pixels beyond the
+ new best corner. */
+ i = best_corner_index = q;
+ free_index_list(&equally_good_list);
+ equally_good_list = new_index_list();
+ }
+
+ i++;
+ }
+
+ /* After we exit the loop, `q' is the index of the last point
+ we checked. We have already added the corner if
+ `best_corner_angle' is less than `corner_always_threshold'.
+ Again, if we've already wrapped around, we don't want to
+ add the corner again. */
+ if (best_corner_angle > fitting_opts->corner_always_threshold && best_corner_index >= p) {
+ unsigned j;
+
+ APPEND_CORNER(best_corner_index, best_corner_angle, '/');
+
+ for (j = 0; j < INDEX_LIST_LENGTH(equally_good_list); j++)
+ APPEND_CORNER(GET_INDEX(equally_good_list, j), best_corner_angle, '@');
+ }
+ free_index_list(&equally_good_list);
+
+ /* If we wrapped around in our search, we're done; otherwise,
+ we don't want the outer loop to look at the pixels that we
+ already looked at in searching for the best corner. */
+ p = (q < p) ? O_LENGTH(pixel_outline) : q;
+ } /* End of searching for the best corner. */
+ } /* End of considering each pixel. */
+
+ if (INDEX_LIST_LENGTH(corner_list) > 0)
+ /* We never want two corners next to each other, since the
+ only way to fit such a ``curve'' would be with a straight
+ line, which usually interrupts the continuity dreadfully. */
+ remove_adjacent_corners(&corner_list, O_LENGTH(pixel_outline) - (pixel_outline.open ? 2 : 1), fitting_opts->remove_adjacent_corners, exception);
+cleanup:
+ return corner_list;
+}
+
+/* Return the difference vectors coming in and going out of the outline
+ OUTLINE at the point whose index is TEST_INDEX. In Phoenix,
+ Schneider looks at a single point on either side of the point we're
+ considering. That works for him because his points are not touching.
+ But our points *are* touching, and so we have to look at
+ `corner_surround' points on either side, to get a better picture of
+ the outline's shape. */
+
+static void find_vectors(unsigned test_index, pixel_outline_type outline, vector_type * in, vector_type * out, unsigned corner_surround)
+{
+ int i;
+ unsigned n_done;
+ at_coord candidate = O_COORDINATE(outline, test_index);
+
+ in->dx = in->dy = in->dz = 0.0;
+ out->dx = out->dy = out->dz = 0.0;
+
+ /* Add up the differences from p of the `corner_surround' points
+ before p. */
+ for (i = O_PREV(outline, test_index), n_done = 0; n_done < corner_surround; i = O_PREV(outline, i), n_done++)
+ *in = Vadd(*in, IPsubtract(O_COORDINATE(outline, i), candidate));
+
+ /* And the points after p. */
+ for (i = O_NEXT(outline, test_index), n_done = 0; n_done < corner_surround; i = O_NEXT(outline, i), n_done++)
+ *out = Vadd(*out, IPsubtract(O_COORDINATE(outline, i), candidate));
+}
+
+/* Remove adjacent points from the index list LIST. We do this by first
+ sorting the list and then running through it. Since these lists are
+ quite short, a straight selection sort (e.g., p.139 of the Art of
+ Computer Programming, vol.3) is good enough. LAST_INDEX is the index
+ of the last pixel on the outline, i.e., the next one is the first
+ pixel. We need this for checking the adjacency of the last corner.
+
+ We need to do this because the adjacent corners turn into
+ two-pixel-long curves, which can only be fit by straight lines. */
+
+static void remove_adjacent_corners(index_list_type * list, unsigned last_index, gboolean remove_adj_corners, at_exception_type * exception)
+{
+ unsigned j;
+ unsigned last;
+ index_list_type new_list = new_index_list();
+
+ for (j = INDEX_LIST_LENGTH(*list) - 1; j > 0; j--) {
+ unsigned search;
+ unsigned temp;
+ /* Find maximal element below `j'. */
+ unsigned max_index = j;
+
+ for (search = 0; search < j; search++)
+ if (GET_INDEX(*list, search) > GET_INDEX(*list, max_index))
+ max_index = search;
+
+ if (max_index != j) {
+ temp = GET_INDEX(*list, j);
+ GET_INDEX(*list, j) = GET_INDEX(*list, max_index);
+ GET_INDEX(*list, max_index) = temp;
+
+ /* xx -- really have to sort? */
+ LOG("needed exchange");
+ at_exception_warning(exception, "needed exchange");
+ }
+ }
+
+ /* The list is sorted. Now look for adjacent entries. Each time
+ through the loop we insert the current entry and, if appropriate,
+ the next entry. */
+ for (j = 0; j < INDEX_LIST_LENGTH(*list) - 1; j++) {
+ unsigned current = GET_INDEX(*list, j);
+ unsigned next = GET_INDEX(*list, j + 1);
+
+ /* We should never have inserted the same element twice. */
+ /* assert (current != next); */
+
+ if ((remove_adj_corners) && ((next == current + 1) || (next == current)))
+ j++;
+
+ append_index(&new_list, current);
+ }
+
+ /* Don't append the last element if it is 1) adjacent to the previous
+ one; or 2) adjacent to the very first one. */
+ last = GET_LAST_INDEX(*list);
+ if (INDEX_LIST_LENGTH(new_list) == 0 || !(last == GET_LAST_INDEX(new_list) + 1 || (last == last_index && GET_INDEX(*list, 0) == 0)))
+ append_index(&new_list, last);
+
+ free_index_list(list);
+ *list = new_list;
+}
+
+/* A ``knee'' is a point which forms a ``right angle'' with its
+ predecessor and successor. See the documentation (the `Removing
+ knees' section) for an example and more details.
+
+ The argument CLOCKWISE tells us which direction we're moving. (We
+ can't figure that information out from just the single segment with
+ which we are given to work.)
+
+ We should never find two consecutive knees.
+
+ Since the first and last points are corners (unless the curve is
+ cyclic), it doesn't make sense to remove those. */
+
+/* This evaluates to TRUE if the vector V is zero in one direction and
+ nonzero in the other. */
+#define ONLY_ONE_ZERO(v) \
+ (((v).dx == 0.0 && (v).dy != 0.0) || ((v).dy == 0.0 && (v).dx != 0.0))
+
+/* There are four possible cases for knees, one for each of the four
+ corners of a rectangle; and then the cases differ depending on which
+ direction we are going around the curve. The tests are listed here
+ in the order of upper left, upper right, lower right, lower left.
+ Perhaps there is some simple pattern to the
+ clockwise/counterclockwise differences, but I don't see one. */
+#define CLOCKWISE_KNEE(prev_delta, next_delta) \
+ ((prev_delta.dx == -1.0 && next_delta.dy == 1.0) \
+ || (prev_delta.dy == 1.0 && next_delta.dx == 1.0) \
+ || (prev_delta.dx == 1.0 && next_delta.dy == -1.0) \
+ || (prev_delta.dy == -1.0 && next_delta.dx == -1.0))
+
+#define COUNTERCLOCKWISE_KNEE(prev_delta, next_delta) \
+ ((prev_delta.dy == 1.0 && next_delta.dx == -1.0) \
+ || (prev_delta.dx == 1.0 && next_delta.dy == 1.0) \
+ || (prev_delta.dy == -1.0 && next_delta.dx == 1.0) \
+ || (prev_delta.dx == -1.0 && next_delta.dy == -1.0))
+
+static void remove_knee_points(curve_type curve, gboolean clockwise)
+{
+ unsigned i;
+ unsigned offset = (CURVE_CYCLIC(curve) == TRUE) ? 0 : 1;
+ at_coord previous = real_to_int_coord(CURVE_POINT(curve, CURVE_PREV(curve, offset)));
+ curve_type trimmed_curve = copy_most_of_curve(curve);
+
+ if (CURVE_CYCLIC(curve) == FALSE)
+ append_pixel(trimmed_curve, real_to_int_coord(CURVE_POINT(curve, 0)));
+
+ for (i = offset; i < CURVE_LENGTH(curve) - offset; i++) {
+ at_coord current = real_to_int_coord(CURVE_POINT(curve, i));
+ at_coord next = real_to_int_coord(CURVE_POINT(curve, CURVE_NEXT(curve, i)));
+ vector_type prev_delta = IPsubtract(previous, current);
+ vector_type next_delta = IPsubtract(next, current);
+
+ if (ONLY_ONE_ZERO(prev_delta) && ONLY_ONE_ZERO(next_delta)
+ && ((clockwise && CLOCKWISE_KNEE(prev_delta, next_delta))
+ || (!clockwise && COUNTERCLOCKWISE_KNEE(prev_delta, next_delta))))
+ LOG(" (%d,%d)", current.x, current.y);
+ else {
+ previous = current;
+ append_pixel(trimmed_curve, current);
+ }
+ }
+
+ if (CURVE_CYCLIC(curve) == FALSE)
+ append_pixel(trimmed_curve, real_to_int_coord(LAST_CURVE_POINT(curve)));
+
+ if (CURVE_LENGTH(trimmed_curve) == CURVE_LENGTH(curve))
+ LOG(" (none)");
+
+ LOG(".\n");
+
+ free_curve(curve);
+ *curve = *trimmed_curve;
+ free(trimmed_curve); /* free_curve? --- Masatake */
+}
+
+/* Smooth the curve by adding in neighboring points. Do this
+ `filter_iterations' times. But don't change the corners. */
+
+static void filter(curve_type curve, fitting_opts_type * fitting_opts)
+{
+ unsigned iteration, this_point;
+ unsigned offset = (CURVE_CYCLIC(curve) == TRUE) ? 0 : 1;
+ at_real_coord prev_new_point;
+
+ /* We must have at least three points---the previous one, the current
+ one, and the next one. But if we don't have at least five, we will
+ probably collapse the curve down onto a single point, which means
+ we won't be able to fit it with a spline. */
+ if (CURVE_LENGTH(curve) < 5) {
+ LOG("Length is %u, not enough to filter.\n", CURVE_LENGTH(curve));
+ return;
+ }
+
+ prev_new_point.x = FLT_MAX;
+ prev_new_point.y = FLT_MAX;
+ prev_new_point.z = FLT_MAX;
+
+ for (iteration = 0; iteration < fitting_opts->filter_iterations; iteration++) {
+ curve_type newcurve = copy_most_of_curve(curve);
+ gboolean collapsed = FALSE;
+
+ /* Keep the first point on the curve. */
+ if (offset)
+ append_point(newcurve, CURVE_POINT(curve, 0));
+
+ for (this_point = offset; this_point < CURVE_LENGTH(curve) - offset; this_point++) {
+ vector_type in, out, sum;
+ at_real_coord new_point;
+
+ /* Calculate the vectors in and out, computed by looking at n points
+ on either side of this_point. Experimental it was found that 2 is
+ optimal. */
+
+ signed int prev, prevprev; /* have to be signed */
+ unsigned int next, nextnext;
+ at_real_coord candidate = CURVE_POINT(curve, this_point);
+
+ prev = CURVE_PREV(curve, this_point);
+ prevprev = CURVE_PREV(curve, prev);
+ next = CURVE_NEXT(curve, this_point);
+ nextnext = CURVE_NEXT(curve, next);
+
+ /* Add up the differences from p of the `surround' points
+ before p. */
+ in.dx = in.dy = in.dz = 0.0;
+
+ in = Vadd(in, Psubtract(CURVE_POINT(curve, prev), candidate));
+ if (prevprev >= 0)
+ in = Vadd(in, Psubtract(CURVE_POINT(curve, prevprev), candidate));
+
+ /* And the points after p. Don't use more points after p than we
+ ended up with before it. */
+ out.dx = out.dy = out.dz = 0.0;
+
+ out = Vadd(out, Psubtract(CURVE_POINT(curve, next), candidate));
+ if (nextnext < CURVE_LENGTH(curve))
+ out = Vadd(out, Psubtract(CURVE_POINT(curve, nextnext), candidate));
+
+ /* Start with the old point. */
+ new_point = candidate;
+ sum = Vadd(in, out);
+ /* We added 2*n+2 points, so we have to divide the sum by 2*n+2 */
+ new_point.x += sum.dx / 6;
+ new_point.y += sum.dy / 6;
+ new_point.z += sum.dz / 6;
+ if (fabs(prev_new_point.x - new_point.x) < 0.3 && fabs(prev_new_point.y - new_point.y) < 0.3 && fabs(prev_new_point.z - new_point.z) < 0.3) {
+ collapsed = TRUE;
+ break;
+ }
+
+ /* Put the newly computed point into a separate curve, so it
+ doesn't affect future computation (on this iteration). */
+ append_point(newcurve, prev_new_point = new_point);
+ }
+
+ if (collapsed)
+ free_curve(newcurve);
+ else {
+ /* Just as with the first point, we have to keep the last point. */
+ if (offset)
+ append_point(newcurve, LAST_CURVE_POINT(curve));
+
+ /* Set the original curve to the newly filtered one, and go again. */
+ free_curve(curve);
+ *curve = *newcurve;
+ }
+ free(newcurve);
+ }
+
+ if (logging)
+ log_curve(curve, FALSE);
+}
+
+/* This routine returns the curve fitted to a straight line in a very
+ simple way: make the first and last points on the curve be the
+ endpoints of the line. This simplicity is justified because we are
+ called only on very short curves. */
+
+static spline_list_type *fit_with_line(curve_type curve)
+{
+ spline_type line;
+
+ LOG("Fitting with straight line:\n");
+
+ SPLINE_DEGREE(line) = LINEARTYPE;
+ START_POINT(line) = CONTROL1(line) = CURVE_POINT(curve, 0);
+ END_POINT(line) = CONTROL2(line) = LAST_CURVE_POINT(curve);
+
+ /* Make sure that this line is never changed to a cubic. */
+ SPLINE_LINEARITY(line) = 0;
+
+ if (logging) {
+ LOG(" ");
+ print_spline(line);
+ }
+
+ return new_spline_list_with_spline(line);
+}
+
+/* The least squares method is well described in Schneider's thesis.
+ Briefly, we try to fit the entire curve with one spline. If that
+ fails, we subdivide the curve. */
+
+static spline_list_type *fit_with_least_squares(curve_type curve, fitting_opts_type * fitting_opts, at_exception_type * exception)
+{
+ gfloat error = 0, best_error = FLT_MAX;
+ spline_type spline, best_spline;
+ spline_list_type *spline_list = NULL;
+ unsigned worst_point = 0;
+ gfloat previous_error = FLT_MAX;
+
+ LOG("\nFitting with least squares:\n");
+
+ /* Phoenix reduces the number of points with a ``linear spline
+ technique''. But for fitting letterforms, that is
+ inappropriate. We want all the points we can get. */
+
+ /* It makes no difference whether we first set the `t' values or
+ find the tangents. This order makes the documentation a little
+ more coherent. */
+
+ LOG("Finding tangents:\n");
+ find_tangent(curve, /* to_start */ TRUE, /* cross_curve */ FALSE,
+ fitting_opts->tangent_surround);
+ find_tangent(curve, /* to_start */ FALSE, /* cross_curve */ FALSE,
+ fitting_opts->tangent_surround);
+
+ set_initial_parameter_values(curve);
+
+ /* Now we loop, subdividing, until CURVE has
+ been fit. */
+ while (TRUE) {
+ spline = best_spline = fit_one_spline(curve, exception);
+ if (at_exception_got_fatal(exception))
+ goto cleanup;
+
+ if (SPLINE_DEGREE(spline) == LINEARTYPE)
+ LOG(" fitted to line:\n");
+ else
+ LOG(" fitted to spline:\n");
+
+ if (logging) {
+ LOG(" ");
+ print_spline(spline);
+ }
+
+ if (SPLINE_DEGREE(spline) == LINEARTYPE)
+ break;
+
+ error = find_error(curve, spline, &worst_point, exception);
+ if (error <= previous_error) {
+ best_error = error;
+ best_spline = spline;
+ }
+ break;
+ }
+
+ if (SPLINE_DEGREE(spline) == LINEARTYPE) {
+ spline_list = new_spline_list_with_spline(spline);
+ LOG("Accepted error of %.3f.\n", error);
+ return (spline_list);
+ }
+
+ /* Go back to the best fit. */
+ spline = best_spline;
+ error = best_error;
+
+ if (error < fitting_opts->error_threshold && CURVE_CYCLIC(curve) == FALSE) {
+ /* The points were fitted with a
+ spline. We end up here whenever a fit is accepted. We have
+ one more job: see if the ``curve'' that was fit should really
+ be a straight line. */
+ if (spline_linear_enough(&spline, curve, fitting_opts)) {
+ SPLINE_DEGREE(spline) = LINEARTYPE;
+ LOG("Changed to line.\n");
+ }
+ spline_list = new_spline_list_with_spline(spline);
+ LOG("Accepted error of %.3f.\n", error);
+ } else {
+ /* We couldn't fit the curve acceptably, so subdivide. */
+ unsigned subdivision_index;
+ spline_list_type *left_spline_list;
+ spline_list_type *right_spline_list;
+ curve_type left_curve = new_curve();
+ curve_type right_curve = new_curve();
+
+ /* Keep the linked list of curves intact. */
+ NEXT_CURVE(right_curve) = NEXT_CURVE(curve);
+ PREVIOUS_CURVE(right_curve) = left_curve;
+ NEXT_CURVE(left_curve) = right_curve;
+ PREVIOUS_CURVE(left_curve) = curve;
+ NEXT_CURVE(curve) = left_curve;
+
+ LOG("\nSubdividing (error %.3f):\n", error);
+ LOG(" Original point: (%.3f,%.3f), #%u.\n", CURVE_POINT(curve, worst_point).x, CURVE_POINT(curve, worst_point).y, worst_point);
+ subdivision_index = worst_point;
+ LOG(" Final point: (%.3f,%.3f), #%u.\n", CURVE_POINT(curve, subdivision_index).x, CURVE_POINT(curve, subdivision_index).y, subdivision_index);
+
+ /* The last point of the left-hand curve will also be the first
+ point of the right-hand curve. */
+ CURVE_LENGTH(left_curve) = subdivision_index + 1;
+ CURVE_LENGTH(right_curve) = CURVE_LENGTH(curve) - subdivision_index;
+ left_curve->point_list = curve->point_list;
+ right_curve->point_list = curve->point_list + subdivision_index;
+
+ /* We want to use the tangents of the curve which we are
+ subdividing for the start tangent for left_curve and the
+ end tangent for right_curve. */
+ CURVE_START_TANGENT(left_curve) = CURVE_START_TANGENT(curve);
+ CURVE_END_TANGENT(right_curve) = CURVE_END_TANGENT(curve);
+
+ /* We have to set up the two curves before finding the tangent at
+ the subdivision point. The tangent at that point must be the
+ same for both curves, or noticeable bumps will occur in the
+ character. But we want to use information on both sides of the
+ point to compute the tangent, hence cross_curve = true. */
+ find_tangent(left_curve, /* to_start_point: */ FALSE,
+ /* cross_curve: */ TRUE, fitting_opts->tangent_surround);
+ CURVE_START_TANGENT(right_curve) = CURVE_END_TANGENT(left_curve);
+
+ /* Now that we've set up the curves, we can fit them. */
+ left_spline_list = fit_curve(left_curve, fitting_opts, exception);
+ if (at_exception_got_fatal(exception))
+ /* TODO: Memory allocated for left_curve and right_curve
+ will leak. */
+ goto cleanup;
+
+ right_spline_list = fit_curve(right_curve, fitting_opts, exception);
+ /* TODO: Memory allocated for left_curve and right_curve
+ will leak. */
+ if (at_exception_got_fatal(exception))
+ goto cleanup;
+
+ /* Neither of the subdivided curves could be fit, so fail. */
+ if (left_spline_list == NULL && right_spline_list == NULL)
+ return NULL;
+
+ /* Put the two together (or whichever of them exist). */
+ spline_list = new_spline_list();
+
+ if (left_spline_list == NULL) {
+ LOG("Could not fit spline to left curve (%lx).\n", (unsigned long)left_curve);
+ at_exception_warning(exception, "Could not fit left spline list");
+ } else {
+ concat_spline_lists(spline_list, *left_spline_list);
+ free_spline_list(*left_spline_list);
+ free(left_spline_list);
+ }
+
+ if (right_spline_list == NULL) {
+ LOG("Could not fit spline to right curve (%lx).\n", (unsigned long)right_curve);
+ at_exception_warning(exception, "Could not fit right spline list");
+ } else {
+ concat_spline_lists(spline_list, *right_spline_list);
+ free_spline_list(*right_spline_list);
+ free(right_spline_list);
+ }
+ if (CURVE_END_TANGENT(left_curve))
+ free(CURVE_END_TANGENT(left_curve));
+ free(left_curve);
+ free(right_curve);
+ }
+cleanup:
+ return spline_list;
+}
+
+/* Our job here is to find alpha1 (and alpha2), where t1_hat (t2_hat) is
+ the tangent to CURVE at the starting (ending) point, such that:
+
+ control1 = alpha1*t1_hat + starting point
+ control2 = alpha2*t2_hat + ending_point
+
+ and the resulting spline (starting_point .. control1 and control2 ..
+ ending_point) minimizes the least-square error from CURVE.
+
+ See pp.57--59 of the Phoenix thesis.
+
+ The B?(t) here corresponds to B_i^3(U_i) there.
+ The Bernshte\u in polynomials of degree n are defined by
+ B_i^n(t) = { n \choose i } t^i (1-t)^{n-i}, i = 0..n */
+
+#define B0(t) CUBE ((gfloat) 1.0 - (t))
+#define B1(t) ((gfloat) 3.0 * (t) * SQUARE ((gfloat) 1.0 - (t)))
+#define B2(t) ((gfloat) 3.0 * SQUARE (t) * ((gfloat) 1.0 - (t)))
+#define B3(t) CUBE (t)
+
+static spline_type fit_one_spline(curve_type curve, at_exception_type * exception)
+{
+ /* Since our arrays are zero-based, the `C0' and `C1' here correspond
+ to `C1' and `C2' in the paper. */
+ gfloat X_C1_det, C0_X_det, C0_C1_det;
+ gfloat alpha1, alpha2;
+ spline_type spline;
+ vector_type start_vector, end_vector;
+ unsigned i;
+ vector_type *A;
+ vector_type t1_hat = *CURVE_START_TANGENT(curve);
+ vector_type t2_hat = *CURVE_END_TANGENT(curve);
+ gfloat C[2][2] = { {0.0, 0.0}, {0.0, 0.0} };
+ gfloat X[2] = { 0.0, 0.0 };
+
+ XMALLOC(A, CURVE_LENGTH(curve) * 2 * sizeof(vector_type)); /* A dynamically allocated array. */
+
+ START_POINT(spline) = CURVE_POINT(curve, 0);
+ END_POINT(spline) = LAST_CURVE_POINT(curve);
+ start_vector = make_vector(START_POINT(spline));
+ end_vector = make_vector(END_POINT(spline));
+
+ for (i = 0; i < CURVE_LENGTH(curve); i++) {
+ A[(i << 1) + 0] = Vmult_scalar(t1_hat, B1(CURVE_T(curve, i)));
+ A[(i << 1) + 1] = Vmult_scalar(t2_hat, B2(CURVE_T(curve, i)));
+ }
+
+ for (i = 0; i < CURVE_LENGTH(curve); i++) {
+ vector_type temp, temp0, temp1, temp2, temp3;
+ vector_type *Ai = A + (i << 1);
+
+ C[0][0] += Vdot(Ai[0], Ai[0]);
+ C[0][1] += Vdot(Ai[0], Ai[1]);
+ /* C[1][0] = C[0][1] (this is assigned outside the loop) */
+ C[1][1] += Vdot(Ai[1], Ai[1]);
+
+ /* Now the right-hand side of the equation in the paper. */
+ temp0 = Vmult_scalar(start_vector, B0(CURVE_T(curve, i)));
+ temp1 = Vmult_scalar(start_vector, B1(CURVE_T(curve, i)));
+ temp2 = Vmult_scalar(end_vector, B2(CURVE_T(curve, i)));
+ temp3 = Vmult_scalar(end_vector, B3(CURVE_T(curve, i)));
+
+ temp = make_vector(Vsubtract_point(CURVE_POINT(curve, i), Vadd(temp0, Vadd(temp1, Vadd(temp2, temp3)))));
+
+ X[0] += Vdot(temp, Ai[0]);
+ X[1] += Vdot(temp, Ai[1]);
+ }
+ free(A);
+
+ C[1][0] = C[0][1];
+
+ X_C1_det = X[0] * C[1][1] - X[1] * C[0][1];
+ C0_X_det = C[0][0] * X[1] - C[0][1] * X[0];
+ C0_C1_det = C[0][0] * C[1][1] - C[1][0] * C[0][1];
+ if (C0_C1_det == 0.0) {
+ /* Zero determinant */
+ alpha1 = 0;
+ alpha2 = 0;
+ } else {
+ alpha1 = X_C1_det / C0_C1_det;
+ alpha2 = C0_X_det / C0_C1_det;
+ }
+ CONTROL1(spline) = Vadd_point(START_POINT(spline), Vmult_scalar(t1_hat, alpha1));
+ CONTROL2(spline) = Vadd_point(END_POINT(spline), Vmult_scalar(t2_hat, alpha2));
+ SPLINE_DEGREE(spline) = CUBICTYPE;
+
+ return spline;
+}
+
+/* Find reasonable values for t for each point on CURVE. The method is
+ called chord-length parameterization, which is described in Plass &
+ Stone. The basic idea is just to use the distance from one point to
+ the next as the t value, normalized to produce values that increase
+ from zero for the first point to one for the last point. */
+
+static void set_initial_parameter_values(curve_type curve)
+{
+ unsigned p;
+
+ LOG("\nAssigning initial t values:\n ");
+
+ CURVE_T(curve, 0) = 0.0;
+
+ for (p = 1; p < CURVE_LENGTH(curve); p++) {
+ at_real_coord point = CURVE_POINT(curve, p), previous_p = CURVE_POINT(curve, p - 1);
+ gfloat d = distance(point, previous_p);
+ CURVE_T(curve, p) = CURVE_T(curve, p - 1) + d;
+ }
+
+ if (LAST_CURVE_T(curve) == 0.0)
+ LAST_CURVE_T(curve) = 1.0;
+
+ for (p = 1; p < CURVE_LENGTH(curve); p++)
+ CURVE_T(curve, p) = CURVE_T(curve, p) / LAST_CURVE_T(curve);
+
+ if (logging)
+ log_entire_curve(curve);
+}
+
+/* Find an approximation to the tangent to an endpoint of CURVE (to the
+ first point if TO_START_POINT is TRUE, else the last). If
+ CROSS_CURVE is TRUE, consider points on the adjacent curve to CURVE.
+
+ It is important to compute an accurate approximation, because the
+ control points that we eventually decide upon to fit the curve will
+ be placed on the half-lines defined by the tangents and
+ endpoints...and we never recompute the tangent after this. */
+
+static void find_tangent(curve_type curve, gboolean to_start_point, gboolean cross_curve, unsigned tangent_surround)
+{
+ vector_type tangent;
+ vector_type **curve_tangent = (to_start_point == TRUE) ? &(CURVE_START_TANGENT(curve))
+ : &(CURVE_END_TANGENT(curve));
+ unsigned n_points = 0;
+
+ LOG(" tangent to %s: ", (to_start_point == TRUE) ? "start" : "end");
+
+ if (*curve_tangent == NULL) {
+ XMALLOC(*curve_tangent, sizeof(vector_type));
+ do {
+ tangent = find_half_tangent(curve, to_start_point, &n_points, tangent_surround);
+
+ if ((cross_curve == TRUE) || (CURVE_CYCLIC(curve) == TRUE)) {
+ curve_type adjacent_curve = (to_start_point == TRUE) ? PREVIOUS_CURVE(curve) : NEXT_CURVE(curve);
+ vector_type tangent2 = (to_start_point == FALSE) ? find_half_tangent(adjacent_curve, TRUE, &n_points,
+ tangent_surround) : find_half_tangent(adjacent_curve, TRUE, &n_points,
+ tangent_surround);
+
+ LOG("(adjacent curve half tangent (%.3f,%.3f,%.3f)) ", tangent2.dx, tangent2.dy, tangent2.dz);
+ tangent = Vadd(tangent, tangent2);
+ }
+ tangent_surround--;
+
+ }
+ while (tangent.dx == 0.0 && tangent.dy == 0.0);
+
+ assert(n_points > 0);
+ **curve_tangent = Vmult_scalar(tangent, (gfloat) (1.0 / n_points));
+ if ((CURVE_CYCLIC(curve) == TRUE) && CURVE_START_TANGENT(curve))
+ *CURVE_START_TANGENT(curve) = **curve_tangent;
+ if ((CURVE_CYCLIC(curve) == TRUE) && CURVE_END_TANGENT(curve))
+ *CURVE_END_TANGENT(curve) = **curve_tangent;
+ } else
+ LOG("(already computed) ");
+
+ LOG("(%.3f,%.3f,%.3f).\n", (*curve_tangent)->dx, (*curve_tangent)->dy, (*curve_tangent)->dz);
+}
+
+/* Find the change in y and change in x for `tangent_surround' (a global)
+ points along CURVE. Increment N_POINTS by the number of points we
+ actually look at. */
+
+static vector_type find_half_tangent(curve_type c, gboolean to_start_point, unsigned *n_points, unsigned tangent_surround)
+{
+ unsigned p;
+ int factor = to_start_point ? 1 : -1;
+ unsigned tangent_index = to_start_point ? 0 : c->length - 1;
+ at_real_coord tangent_point = CURVE_POINT(c, tangent_index);
+ vector_type tangent = { 0.0, 0.0 };
+ unsigned int surround;
+
+ if ((surround = CURVE_LENGTH(c) / 2) > tangent_surround)
+ surround = tangent_surround;
+
+ for (p = 1; p <= surround; p++) {
+ int this_index = p * factor + tangent_index;
+ at_real_coord this_point;
+
+ if (this_index < 0 || this_index >= (int)c->length)
+ break;
+
+ this_point = CURVE_POINT(c, p * factor + tangent_index);
+
+ /* Perhaps we should weight the tangent from `this_point' by some
+ factor dependent on the distance from the tangent point. */
+ tangent = Vadd(tangent, Vmult_scalar(Psubtract(this_point, tangent_point), (gfloat) factor));
+ (*n_points)++;
+ }
+
+ return tangent;
+}
+
+/* When this routine is called, we have computed a spline representation
+ for the digitized curve. The question is, how good is it? If the
+ fit is very good indeed, we might have an error of zero on each
+ point, and then WORST_POINT becomes irrelevant. But normally, we
+ return the error at the worst point, and the index of that point in
+ WORST_POINT. The error computation itself is the Euclidean distance
+ from the original curve CURVE to the fitted spline SPLINE. */
+
+static gfloat find_error(curve_type curve, spline_type spline, unsigned *worst_point, at_exception_type * exception)
+{
+ unsigned this_point;
+ gfloat total_error = 0.0;
+ gfloat worst_error = FLT_MIN;
+
+ *worst_point = CURVE_LENGTH(curve) + 1; /* A sentinel value. */
+
+ for (this_point = 0; this_point < CURVE_LENGTH(curve); this_point++) {
+ at_real_coord curve_point = CURVE_POINT(curve, this_point);
+ gfloat t = CURVE_T(curve, this_point);
+ at_real_coord spline_point = evaluate_spline(spline, t);
+ gfloat this_error = distance(curve_point, spline_point);
+ if (this_error >= worst_error) {
+ *worst_point = this_point;
+ worst_error = this_error;
+ }
+ total_error += this_error;
+ }
+
+ if (*worst_point == CURVE_LENGTH(curve) + 1) { /* Didn't have any ``worst point''; the error should be zero. */
+ if (epsilon_equal(total_error, 0.0))
+ LOG(" Every point fit perfectly.\n");
+ else {
+ LOG("No worst point found; something is wrong");
+ at_exception_warning(exception, "No worst point found; something is wrong");
+ }
+ } else {
+ if (epsilon_equal(total_error, 0.0))
+ LOG(" Every point fit perfectly.\n");
+ else {
+ LOG(" Worst error (at (%.3f,%.3f,%.3f), point #%u) was %.3f.\n", CURVE_POINT(curve, *worst_point).x, CURVE_POINT(curve, *worst_point).y, CURVE_POINT(curve, *worst_point).z, *worst_point, worst_error);
+ LOG(" Total error was %.3f.\n", total_error);
+ LOG(" Average error (over %u points) was %.3f.\n", CURVE_LENGTH(curve), total_error / CURVE_LENGTH(curve));
+ }
+ }
+
+ return worst_error;
+}
+
+/* Supposing that we have accepted the error, another question arises:
+ would we be better off just using a straight line? */
+
+static gboolean spline_linear_enough(spline_type * spline, curve_type curve, fitting_opts_type * fitting_opts)
+{
+ gfloat A, B, C;
+ unsigned this_point;
+ gfloat dist = 0.0, start_end_dist, threshold;
+
+ LOG("Checking linearity:\n");
+
+ A = END_POINT(*spline).x - START_POINT(*spline).x;
+ B = END_POINT(*spline).y - START_POINT(*spline).y;
+ C = END_POINT(*spline).z - START_POINT(*spline).z;
+
+ start_end_dist = (gfloat) (SQUARE(A) + SQUARE(B) + SQUARE(C));
+ LOG("start_end_distance is %.3f.\n", sqrt(start_end_dist));
+
+ LOG(" Line endpoints are (%.3f, %.3f, %.3f) and ", START_POINT(*spline).x, START_POINT(*spline).y, START_POINT(*spline).z);
+ LOG("(%.3f, %.3f, %.3f)\n", END_POINT(*spline).x, END_POINT(*spline).y, END_POINT(*spline).z);
+
+ /* LOG (" Line is %.3fx + %.3fy + %.3f = 0.\n", A, B, C); */
+
+ for (this_point = 0; this_point < CURVE_LENGTH(curve); this_point++) {
+ gfloat a, b, c, w;
+ gfloat t = CURVE_T(curve, this_point);
+ at_real_coord spline_point = evaluate_spline(*spline, t);
+
+ a = spline_point.x - START_POINT(*spline).x;
+ b = spline_point.y - START_POINT(*spline).y;
+ c = spline_point.z - START_POINT(*spline).z;
+ w = (A * a + B * b + C * c) / start_end_dist;
+
+ dist += (gfloat) sqrt(SQUARE(a - A * w) + SQUARE(b - B * w) + SQUARE(c - C * w));
+ }
+ LOG(" Total distance is %.3f, ", dist);
+
+ dist /= (CURVE_LENGTH(curve) - 1);
+ LOG("which is %.3f normalized.\n", dist);
+
+ /* We want reversion of short curves to splines to be more likely than
+ reversion of long curves, hence the second division by the curve
+ length, for use in `change_bad_lines'. */
+ SPLINE_LINEARITY(*spline) = dist;
+ LOG(" Final linearity: %.3f.\n", SPLINE_LINEARITY(*spline));
+ if (start_end_dist * (gfloat) 0.5 > fitting_opts->line_threshold)
+ threshold = fitting_opts->line_threshold;
+ else
+ threshold = start_end_dist * (gfloat) 0.5;
+ LOG("threshold is %.3f .\n", threshold);
+ if (dist < threshold)
+ return TRUE;
+ else
+ return FALSE;
+}
+
+/* Unfortunately, we cannot tell in isolation whether a given spline
+ should be changed to a line or not. That can only be known after the
+ entire curve has been fit to a list of splines. (The curve is the
+ pixel outline between two corners.) After subdividing the curve, a
+ line may very well fit a portion of the curve just as well as the
+ spline---but unless a spline is truly close to being a line, it
+ should not be combined with other lines. */
+
+static void change_bad_lines(spline_list_type * spline_list, fitting_opts_type * fitting_opts)
+{
+ unsigned this_spline;
+ gboolean found_cubic = FALSE;
+ unsigned length = SPLINE_LIST_LENGTH(*spline_list);
+
+ LOG("\nChecking for bad lines (length %u):\n", length);
+
+ /* First see if there are any splines in the fitted shape. */
+ for (this_spline = 0; this_spline < length; this_spline++) {
+ if (SPLINE_DEGREE(SPLINE_LIST_ELT(*spline_list, this_spline)) == CUBICTYPE) {
+ found_cubic = TRUE;
+ break;
+ }
+ }
+
+ /* If so, change lines back to splines (we haven't done anything to
+ their control points, so we only have to change the degree) unless
+ the spline is close enough to being a line. */
+ if (found_cubic)
+ for (this_spline = 0; this_spline < length; this_spline++) {
+ spline_type s = SPLINE_LIST_ELT(*spline_list, this_spline);
+
+ if (SPLINE_DEGREE(s) == LINEARTYPE) {
+ LOG(" #%u: ", this_spline);
+ if (SPLINE_LINEARITY(s) > fitting_opts->line_reversion_threshold) {
+ LOG("reverted, ");
+ SPLINE_DEGREE(SPLINE_LIST_ELT(*spline_list, this_spline))
+ = CUBICTYPE;
+ }
+ LOG("linearity %.3f.\n", SPLINE_LINEARITY(s));
+ }
+ } else
+ LOG(" No lines.\n");
+}
+
+/* Lists of array indices (well, that is what we use it for). */
+
+static index_list_type new_index_list(void)
+{
+ index_list_type index_list;
+
+ index_list.data = NULL;
+ INDEX_LIST_LENGTH(index_list) = 0;
+
+ return index_list;
+}
+
+static void free_index_list(index_list_type * index_list)
+{
+ if (INDEX_LIST_LENGTH(*index_list) > 0) {
+ free(index_list->data);
+ index_list->data = NULL;
+ INDEX_LIST_LENGTH(*index_list) = 0;
+ }
+}
+
+static void append_index(index_list_type * list, unsigned new_index)
+{
+ INDEX_LIST_LENGTH(*list)++;
+ XREALLOC(list->data, INDEX_LIST_LENGTH(*list) * sizeof(unsigned));
+ list->data[INDEX_LIST_LENGTH(*list) - 1] = new_index;
+}
+
+/* Turn an real point into a integer one. */
+
+static at_coord real_to_int_coord(at_real_coord real_coord)
+{
+ at_coord int_coord;
+
+ int_coord.x = lround(real_coord.x);
+ int_coord.y = lround(real_coord.y);
+
+ return int_coord;
+}
+
+/* Return the Euclidean distance between P1 and P2. */
+
+static gfloat distance(at_real_coord p1, at_real_coord p2)
+{
+ gfloat x = p1.x - p2.x, y = p1.y - p2.y, z = p1.z - p2.z;
+ return (gfloat) sqrt(SQUARE(x)
+ + SQUARE(y) + SQUARE(z));
+}