diff options
Diffstat (limited to 'plug-ins/selection-to-path/fit.c')
| -rw-r--r-- | plug-ins/selection-to-path/fit.c | 1967 |
1 files changed, 1967 insertions, 0 deletions
diff --git a/plug-ins/selection-to-path/fit.c b/plug-ins/selection-to-path/fit.c new file mode 100644 index 0000000..8ee57ae --- /dev/null +++ b/plug-ins/selection-to-path/fit.c @@ -0,0 +1,1967 @@ +/* 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. + * + * 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 3, 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, see <https://www.gnu.org/licenses/>. + */ + +#include "config.h" + +#include <string.h> +#include <float.h> +#include <math.h> +#include <assert.h> + +#include <glib.h> + +#include "global.h" +#include "spline.h" +#include "vector.h" + +#include "curve.h" +#include "fit.h" +#include "pxl-outline.h" + +/* If two endpoints are closer than this, they are made to be equal. + (-align-threshold) */ +real align_threshold = 0.5; + +/* If the angle defined by a point and its predecessors and successors + is smaller than this, it's a corner, even if it's within + `corner_surround' pixels of a point with a smaller angle. + (-corner-always-threshold) */ +real corner_always_threshold = 60.0; + +/* Number of points to consider when determining if a point is a corner + or not. (-corner-surround) */ +unsigned corner_surround = 4; + +/* If a point, its predecessors, and its successors define an angle + smaller than this, it's a corner. Should be in range 0..180. + (-corner-threshold) */ +real corner_threshold = 100.0; + +/* Amount of error at which a fitted spline is unacceptable. If any + pixel is further away than this from the fitted curve, we try again. + (-error-threshold) */ +/* real error_threshold = .8; ALT */ +real error_threshold = .4; + +/* A second number of adjacent points to consider when filtering. + (-filter-alternative-surround) */ +unsigned filter_alternative_surround = 1; + +/* If the angles between the vectors produced by filter_surround and + filter_alternative_surround points differ by more than this, use + the one from filter_alternative_surround. (-filter-epsilon) */ +real filter_epsilon = 10.0; + +/* Number of times to smooth original data points. Increasing this + number dramatically---to 50 or so---can produce vastly better + results. But if any points that ``should'' be corners aren't found, + the curve goes to hell around that point. (-filter-iterations) */ +/* unsigned filter_iteration_count = 4; ALT */ +unsigned filter_iteration_count = 4; + +/* To produce the new point, use the old point plus this times the + neighbors. (-filter-percent) */ +real filter_percent = .33; + +/* Number of adjacent points to consider if `filter_surround' points + defines a straight line. (-filter-secondary-surround) */ +static unsigned filter_secondary_surround = 3; + +/* Number of adjacent points to consider when filtering. + (-filter-surround) */ +unsigned filter_surround = 2; + +/* Says whether or not to remove ``knee'' points after finding the outline. + (See the comments at `remove_knee_points'.) (-remove-knees). */ +boolean keep_knees = false; + +/* If a spline is closer to a straight line than this, it remains a + straight line, even if it would otherwise be changed back to a curve. + This is weighted by the square of the curve length, to make shorter + curves more likely to be reverted. (-line-reversion-threshold) */ +real line_reversion_threshold = .01; + +/* How many pixels (on the average) a spline can diverge from the line + determined by its endpoints before it is changed to a straight line. + (-line-threshold) */ +/* real line_threshold = 1.0; ALT */ +real line_threshold = 0.5; + +/* If reparameterization doesn't improve the fit by this much percent, + stop doing it. (-reparameterize-improve) */ +/* real reparameterize_improvement = .10; ALT */ +real reparameterize_improvement = .01; + +/* Amount of error at which it is pointless to reparameterize. This + happens, for example, when we are trying to fit the outline of the + outside of an `O' with a single spline. The initial fit is not good + enough for the Newton-Raphson iteration to improve it. It may be + that it would be better to detect the cases where we didn't find any + corners. (-reparameterize-threshold) */ +/* real reparameterize_threshold = 30.0; ALT */ +real reparameterize_threshold = 1.0; + +/* Percentage of the curve away from the worst point to look for a + better place to subdivide. (-subdivide-search) */ +real subdivide_search = .1; + +/* Number of points to consider when deciding whether a given point is a + better place to subdivide. (-subdivide-surround) */ +unsigned subdivide_surround = 4; + +/* How many pixels a point can diverge from a straight line and still be + considered a better place to subdivide. (-subdivide-threshold) */ +real subdivide_threshold = .03; + +/* Number of points to look at on either side of a point when computing + the approximation to the tangent at that point. (-tangent-surround) */ +unsigned tangent_surround = 3; + + +/* 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); + +static void align (spline_list_type *); +static void change_bad_lines (spline_list_type *); +static void filter (curve_type); +static real filter_angle (vector_type, vector_type); +static void find_curve_vectors + (unsigned, curve_type, unsigned, vector_type *, vector_type *, unsigned *); +static unsigned find_subdivision (curve_type, unsigned initial); +static void find_vectors + (unsigned, pixel_outline_type, vector_type *, vector_type *); +static index_list_type find_corners (pixel_outline_type); +static real find_error (curve_type, spline_type, unsigned *); +static vector_type find_half_tangent (curve_type, boolean start, unsigned *); +static void find_tangent (curve_type, boolean, boolean); +static spline_type fit_one_spline (curve_type); +static spline_list_type *fit_curve (curve_type); +static spline_list_type fit_curve_list (curve_list_type); +static spline_list_type *fit_with_least_squares (curve_type); +static spline_list_type *fit_with_line (curve_type); +static void remove_knee_points (curve_type, boolean); +static boolean reparameterize (curve_type, spline_type); +static void set_initial_parameter_values (curve_type); +static boolean spline_linear_enough (spline_type *, curve_type); +static curve_list_array_type split_at_corners (pixel_outline_list_type); +static boolean test_subdivision_point (curve_type, unsigned, vector_type *); + +/* 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) +{ + unsigned this_list; + unsigned total = 0; + spline_list_array_type char_splines = new_spline_list_array (); + curve_list_array_type curve_array = split_at_corners (pixel_outline_list); + + 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); + + curve_list_splines = fit_curve_list (curves); + append_spline_list (&char_splines, curve_list_splines); + +/* REPORT ("* "); */ + } + + free_curve_list_array (&curve_array); + + for (this_list = 0; this_list < SPLINE_LIST_ARRAY_LENGTH (char_splines); + this_list++) + total + += SPLINE_LIST_LENGTH (SPLINE_LIST_ARRAY_ELT (char_splines, this_list)); + +/* REPORT1 ("=%u", total); */ + + return char_splines; +} + +/* Set up the internal parameters from the external ones */ + +void +fit_set_params(SELVALS *selVals) +{ + align_threshold = selVals->align_threshold; + corner_always_threshold = selVals->corner_always_threshold; + corner_surround = selVals->corner_surround; + corner_threshold = selVals->corner_threshold; + error_threshold = selVals->error_threshold; + filter_alternative_surround = selVals->filter_alternative_surround; + filter_epsilon = selVals->filter_epsilon; + filter_iteration_count = selVals->filter_iteration_count; + filter_percent = selVals->filter_percent; + filter_secondary_surround = selVals->filter_secondary_surround; + filter_surround = selVals->filter_surround; + keep_knees = selVals->keep_knees; + line_reversion_threshold = selVals->line_reversion_threshold; + line_threshold = selVals->line_threshold; + reparameterize_improvement = selVals->reparameterize_improvement; + reparameterize_threshold = selVals->reparameterize_threshold; + subdivide_search = selVals->subdivide_search; + subdivide_surround = selVals->subdivide_surround; + subdivide_threshold = selVals->subdivide_threshold; + tangent_surround = selVals->tangent_surround; +} + +void +fit_set_default_params(SELVALS *selVals) +{ + selVals->align_threshold = align_threshold; + selVals->corner_always_threshold = corner_always_threshold; + selVals->corner_surround = corner_surround; + selVals->corner_threshold = corner_threshold; + selVals->error_threshold = error_threshold; + selVals->filter_alternative_surround = filter_alternative_surround; + selVals->filter_epsilon = filter_epsilon; + selVals->filter_iteration_count = filter_iteration_count; + selVals->filter_percent = filter_percent; + selVals->filter_secondary_surround = filter_secondary_surround; + selVals->filter_surround = filter_surround; + selVals->keep_knees = keep_knees; + selVals->line_reversion_threshold = line_reversion_threshold; + selVals->line_threshold = line_threshold; + selVals->reparameterize_improvement = reparameterize_improvement; + selVals->reparameterize_threshold = reparameterize_threshold; + selVals->subdivide_search = subdivide_search; + selVals->subdivide_surround = subdivide_surround; + selVals->subdivide_threshold = subdivide_threshold; + selVals->tangent_surround = tangent_surround; +} + +/* 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) +{ + curve_type curve; + unsigned this_curve, this_spline; + unsigned curve_list_length = CURVE_LIST_LENGTH (curve_list); + spline_list_type curve_list_splines = *new_spline_list (); + + /* 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. */ + if (!keep_knees) + { +/* LOG ("\nRemoving knees:\n"); */ + for (this_curve = 0; this_curve < curve_list_length; this_curve++) + { +/* LOG1 ("#%u:", this_curve); */ + remove_knee_points (CURVE_LIST_ELT (curve_list, this_curve), + CURVE_LIST_CLOCKWISE (curve_list)); + } + } + + /* 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++) + { +/* LOG1 ("#%u: ", this_curve); */ + filter (CURVE_LIST_ELT (curve_list, this_curve)); +/* REPORT ("f"); */ + } + + /* 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) && CURVE_LENGTH (curve) != 3) + 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); + +/* REPORT1 (" %u", this_curve); */ +/* LOG1 ("\nFitting curve #%u:\n", this_curve); */ + + curve_splines = fit_curve (current_curve); + if (curve_splines == NULL) + printf("Could not fit curve #%u", this_curve); + else + { +/* LOG1 ("Fitted splines for curve #%u:\n", this_curve); */ + for (this_spline = 0; + this_spline < SPLINE_LIST_LENGTH (*curve_splines); + this_spline++) + { +/* LOG1 (" %u: ", this_spline); */ +/* if (logging) */ +/* print_spline (log_ +file, */ +/* 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); + + concat_spline_lists (&curve_list_splines, *curve_splines); +/* REPORT1 ("(%u)", SPLINE_LIST_LENGTH (*curve_splines)); */ + } + } + + + /* We do this for each outline's spline list because when a point + is changed, it needs to be changed in both segments in which it + appears---and the segments might be in different curves. */ + align (&curve_list_splines); + + 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) +{ + spline_list_type *fitted_splines; + + if (CURVE_LENGTH (curve) < 2) + { + printf ("Tried to fit curve with less than two points"); + return NULL; + } + + /* Do we have enough points to fit with a spline? */ + fitted_splines = CURVE_LENGTH (curve) < 4 + ? fit_with_line (curve) + : fit_with_least_squares (curve); + + return fitted_splines; +} + +/* 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) +{ + 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); + +/* LOG1 ("#%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) > corner_surround * 2 + 2) + { + corner_list = find_corners (pixel_o); + } + else + { + corner_list.data = NULL; + corner_list.length = 0; + } + + /* 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)); + + /* This curve is cyclic. */ + 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)); + + for (p = 0; p <= GET_INDEX (corner_list, 0); p++) + append_pixel (curve, O_COORDINATE (pixel_o, p)); + } + +/* LOG1 (" [%u].\n", corner_list.length); */ + + /* 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); \ + /*LOG4 (" (%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) +{ + unsigned p; + index_list_type corner_list = new_index_list (); + + /* Consider each pixel on the outline in turn. */ + for (p = 0; p < O_LENGTH (pixel_outline); p++) + { + real corner_angle; + vector_type in_vector, out_vector; + + /* Check if the angle is small enough. */ + find_vectors (p, pixel_outline, &in_vector, &out_vector); + corner_angle = Vangle (in_vector, out_vector); + + if (fabs (corner_angle) <= 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. */ + real 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 <= 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 + 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); + corner_angle = Vangle (in_vector, out_vector); + + /* 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 > corner_always_threshold + && best_corner_index >= p) + { + unsigned i; + + APPEND_CORNER (best_corner_index, best_corner_angle, '/'); + + for (i = 0; i < INDEX_LIST_LENGTH (equally_good_list); i++) + APPEND_CORNER (GET_INDEX (equally_good_list, i), + 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) - 1); + + 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) +{ + int i; + unsigned n_done; + coordinate_type candidate = O_COORDINATE (outline, test_index); + + in->dx = 0.0; + in->dy = 0.0; + out->dx = 0.0; + out->dy = 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)); + +#if 0 + /* We don't need this code any more, because now we create the pixel + outline from the corners of the pixels, rather than the edges. */ + + /* To see why we need this test, consider the following + case: four pixels stacked vertically with no other adjacent pixels, + i.e., * + *x + * + * + *** (etc.) We are considering the pixel marked `x' for cornerhood. + The out vector at this point is going to be the zero vector (if + `corner_surround' is 3), because the first + pixel on the outline is the one above the x, the second pixel x + itself, and the third the one below x. (Remember that we go + around the edges of the pixels to find the outlines, not the + pixels themselves.) */ + if (magnitude (*in) == 0.0) + { + WARNING ("Zero magnitude in"); + return corner_threshold + 1.0; + } +#endif /* 0 */ + + /* 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)); + +#if 0 + /* As with the test for the in vector, we don't need this any more. */ + if (magnitude (*out) == 0.0) + { + WARNING ("Zero magnitude out"); + return corner_threshold + 1.0; + } +#endif /* 0 */ +} + + +/* 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) +{ + unsigned j; + unsigned last; + index_list_type new = 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; + printf ("needed exchange"); /* xx -- really have to sort? */ + } + } + + /* 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); + + append_index (&new, current); + if (next == current + 1) + j++; + } + + /* 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) == 0 + || !(last == GET_LAST_INDEX (new) + 1 + || (last == last_index && GET_INDEX (*list, 0) == 0))) + append_index (&new, last); + + free_index_list (list); + *list = new; +} + +/* 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, boolean clockwise) +{ + int i; + unsigned offset = CURVE_CYCLIC (curve) ? 0 : 1; + coordinate_type 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)) + append_pixel (trimmed_curve, real_to_int_coord (CURVE_POINT (curve, 0))); + + for (i = offset; i < CURVE_LENGTH (curve) - offset; i++) + { + coordinate_type current + = real_to_int_coord (CURVE_POINT (curve, i)); + coordinate_type 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)))) + { + /* LOG2 (" (%d,%d)", current.x, current.y); */ + } + else + { + previous = current; + append_pixel (trimmed_curve, current); + } + } + + if (!CURVE_CYCLIC (curve)) + 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; +} + +/* Smooth the curve by adding in neighboring points. Do this + `filter_iteration_count' times. But don't change the corners. */ + +#if 0 +/* Computing the new point based on a single neighboring point and with + constant percentages, as the `SMOOTH' macro did, isn't quite good + enough. For example, at the top of an `o' the curve might well have + three consecutive horizontal pixels, even though there isn't really a + straight there. With this code, the middle point would remain + unfiltered. */ + +#define SMOOTH(axis) \ + CURVE_POINT (curve, this_point).axis \ + = ((1.0 - center_percent) / 2) \ + * CURVE_POINT (curve, CURVE_PREV (curve, this_point)).axis \ + + center_percent * CURVE_POINT (curve, this_point).axis \ + + ((1.0 - center_percent) / 2) \ + * CURVE_POINT (curve, CURVE_NEXT (curve, this_point)).axis +#endif /* 0 */ + +/* We sometimes need to change the information about the filtered point. + This macro assigns to the relevant variables. */ +#define FILTER_ASSIGN(new) \ + do \ + { \ + in = in_##new; \ + out = out_##new; \ + count = new##_count; \ + angle = angle_##new; \ + } \ + while (0) + +static void +filter (curve_type curve) +{ + unsigned iteration, this_point; + unsigned offset = CURVE_CYCLIC (curve) ? 0 : 1; + + /* 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) + { +/* LOG1 ("Length is %u, not enough to filter.\n", CURVE_LENGTH (curve)); */ + return; + } + + for (iteration = 0; iteration < filter_iteration_count; iteration++) + { + curve_type new_curve = copy_most_of_curve (curve); + + /* Keep the first point on the curve. */ + if (offset) + append_point (new_curve, CURVE_POINT (curve, 0)); + + for (this_point = offset; this_point < CURVE_LENGTH (curve) - offset; + this_point++) + { + real angle, angle_alt; + vector_type in, in_alt, out, out_alt, sum; + unsigned count, alt_count; + real_coordinate_type new_point; + + /* Find the angle using the usual number of surrounding points + on the curve. */ + find_curve_vectors (this_point, curve, filter_surround, + &in, &out, &count); + angle = filter_angle (in, out); + + /* Find the angle using the alternative (presumably smaller) + number. */ + find_curve_vectors (this_point, curve, filter_alternative_surround, + &in_alt, &out_alt, &alt_count); + angle_alt = filter_angle (in_alt, out_alt); + + /* If the alternate angle is enough larger than the usual one + and neither of the components of the sum are zero, use it. + (We don't use absolute value here, since if the alternate + angle is smaller, we don't particularly care, since that + means the curve is pretty flat around the current point, + anyway.) This helps keep small features from disappearing + into the surrounding curve. */ + sum = Vadd (in_alt, out_alt); + if (angle_alt - angle >= filter_epsilon + && sum.dx != 0 && sum.dy != 0) + FILTER_ASSIGN (alt); + +#if 0 +/* This code isn't needed anymore, since we do the filtering in a + somewhat more general way. */ + /* If we're left with an angle of zero, don't stop yet; we + might be at a straight which really isn't one (as in the `o' + discussed above). */ + if (epsilon_equal (angle, 0.0)) + { + real angle_secondary; + vector_type in_secondary, out_secondary; + unsigned in_secondary_count, out_secondary_count; + + find_curve_vectors (this_point, curve, filter_secondary_surround, + &in_secondary, &out_secondary, + &in_secondary_count, &out_secondary_count); + angle_secondary = filter_angle (in_secondary, out_secondary); + if (!epsilon_equal (angle_secondary, 0.0)) + FILTER_ASSIGN (secondary); + } +#endif /* 0 */ + + /* Start with the old point. */ + new_point = CURVE_POINT (curve, this_point); + sum = Vadd (in, out); + new_point.x += sum.dx * filter_percent / count; + new_point.y += sum.dy * filter_percent / count; + + /* Put the newly computed point into a separate curve, so it + doesn't affect future computation (on this iteration). */ + append_point (new_curve, new_point); + } + + /* Just as with the first point, we have to keep the last point. */ + if (offset) + append_point (new_curve, LAST_CURVE_POINT (curve)); + + /* Set the original curve to the newly filtered one, and go again. */ + free_curve (curve); + *curve = *new_curve; + } + +/* log_curve (curve, false); */ +/* display_curve (curve); */ +} + + +/* Return the vectors IN and OUT, computed by looking at SURROUND points + on either side of TEST_INDEX. Also return the number of points in + the vectors in COUNT (we make sure they are the same). */ + +static void +find_curve_vectors (unsigned test_index, curve_type curve, + unsigned surround, + vector_type *in, vector_type *out, unsigned *count) +{ + int i; + unsigned in_count, out_count; + unsigned n_done; + real_coordinate_type candidate = CURVE_POINT (curve, test_index); + + /* Add up the differences from p of the `surround' points + before p. */ + + in->dx = 0.0; + in->dy = 0.0; + + for (i = CURVE_PREV (curve, test_index), n_done = 0; + i >= 0 && n_done < surround; /* Do not wrap around. */ + i = CURVE_PREV (curve, i), n_done++) + *in = Vadd (*in, Psubtract (CURVE_POINT (curve, i), candidate)); + in_count = n_done; + + /* And the points after p. Don't use more points after p than we + ended up with before it. */ + out->dx = 0.0; + out->dy = 0.0; + + for (i = CURVE_NEXT (curve, test_index), n_done = 0; + i < CURVE_LENGTH (curve) && n_done < surround && n_done < in_count; + i = CURVE_NEXT (curve, i), n_done++) + *out = Vadd (*out, Psubtract (CURVE_POINT (curve, i), candidate)); + out_count = n_done; + + /* If we used more points before p than after p, we have to go back + and redo it. (We could just subtract the ones that were missing, + but for this few of points, efficiency doesn't matter.) */ + if (out_count < in_count) + { + in->dx = 0.0; + in->dy = 0.0; + + for (i = CURVE_PREV (curve, test_index), n_done = 0; + i >= 0 && n_done < out_count; + i = CURVE_PREV (curve, i), n_done++) + *in = Vadd (*in, Psubtract (CURVE_POINT (curve, i), candidate)); + in_count = n_done; + } + + assert (in_count == out_count); + *count = in_count; +} + + +/* Find the angle between the vectors IN and OUT, and bring it into the + range [0,45]. */ + +static real +filter_angle (vector_type in, vector_type out) +{ + real angle = Vangle (in, out); + + /* What we want to do between 90 and 180 is the same as what we + want to do between 0 and 90. */ + angle = fmod (angle, 1990.0); + + /* And what we want to do between 45 and 90 is the same as + between 0 and 45, only reversed. */ + if (angle > 45.0) angle = 90.0 - angle; + + return angle; +} + +/* 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 = new_spline (); + +/* LOG ("Fitting with straight line:\n"); */ +/* REPORT ("l"); */ + + SPLINE_DEGREE (line) = LINEAR; + 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 (log_file, line); */ +/* } */ + + return init_spline_list (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 try reparameterizing, i.e., finding new, and supposedly better, + t values. If that still fails, we subdivide the curve. */ + +static spline_list_type * +fit_with_least_squares (curve_type curve) +{ + real error, best_error = FLT_MAX; + spline_type spline, best_spline; + spline_list_type *spline_list; + unsigned worst_point; + unsigned iteration = 0; + real previous_error = FLT_MAX; + real improvement = FLT_MAX; + + /* FIXME: Initialize best_spline to zeroes. This is strictly not + necessary as best_spline is always set in the loop below. But the + compiler thinks it isn't and warns. Ideally, the code should be + rewritten such that best_spline and best_error are initialized with + the first values before the loop begins. */ + memset (&best_spline, 0, sizeof best_spline); + +/* 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); + find_tangent (curve, /* to_start */ false, /* cross_curve */ false); + + set_initial_parameter_values (curve); + + /* Now we loop, reparameterizing and/or subdividing, until CURVE has + been fit. */ + while (true) + { +/* LOG (" fitted to spline:\n"); */ + + spline = fit_one_spline (curve); + +/* if (logging) */ +/* { */ +/* LOG (" "); */ +/* print_spline (log_file, spline); */ +/* } */ + + error = find_error (curve, spline, &worst_point); + if (error > previous_error) + { +/* LOG ("Reparameterization made it worse.\n"); */ + /* Just fall through; we will subdivide. */ + } + else + { + best_error = error; + best_spline = spline; + } + + improvement = 1.0 - error / previous_error; + + /* Don't exit, even if the error is less than `error_threshold', + since we might be able to improve the fit with further + reparameterization. If the reparameterization made it worse, + we will exit here, since `improvement' will be negative. */ + if (improvement < reparameterize_improvement + || error > reparameterize_threshold) + break; + + iteration++; +/* LOG3 ("Reparameterization #%u (error was %.3f, a %u%% improvement):\n", */ +/* iteration, error, ((unsigned) (improvement * 100.0))); */ + + /* The reparameterization might fail, if the initial fit was + better than `reparameterize_threshold', yet worse than the + Newton-Raphson algorithm could handle. */ + if (!reparameterize (curve, spline)) + break; + + previous_error = error; + } + + /* Go back to the best fit. */ + spline = best_spline; + error = best_error; + + if (error < error_threshold) + { + /* The points were fitted with a (possibly reparameterized) + 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)) + { + SPLINE_DEGREE (spline) = LINEAR; +/* LOG ("Changed to line.\n"); */ + } + spline_list = init_spline_list (spline); +/* LOG1 ("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; + +/* REPORT ("s"); */ +/* LOG1 ("\nSubdividing (error %.3f):\n", error); */ +/* LOG3 (" Original point: (%.3f,%.3f), #%u.\n", */ +/* CURVE_POINT (curve, worst_point).x, */ +/* CURVE_POINT (curve, worst_point).y, worst_point); */ + subdivision_index = find_subdivision (curve, worst_point); +/* LOG3 (" Final point: (%.3f,%.3f), #%u.\n", */ +/* CURVE_POINT (curve, subdivision_index).x, */ +/* CURVE_POINT (curve, subdivision_index).y, subdivision_index); */ +/* display_subdivision (CURVE_POINT (curve, 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); + 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); + right_spline_list = fit_curve (right_curve); + + /* 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) + { + WARNING ("could not fit left spline list"); +/* LOG1 ("Could not fit spline to left curve (%x).\n", */ +/* (unsigned) left_curve); */ + } + else + concat_spline_lists (spline_list, *left_spline_list); + + if (right_spline_list == NULL) + { + WARNING ("could not fit right spline list"); +/* LOG1 ("Could not fit spline to right curve (%x).\n", */ +/* (unsigned) right_curve); */ + } + else + concat_spline_lists (spline_list, *right_spline_list); + } + + 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 (1 - (t)) +#define B1(t) (3.0 * (t) * SQUARE (1 - (t))) +#define B2(t) (3.0 * SQUARE (t) * (1 - (t))) +#define B3(t) CUBE (t) + +#define U(i) CURVE_T (curve, i) + +static spline_type +fit_one_spline (curve_type curve) +{ + /* Since our arrays are zero-based, the `C0' and `C1' here correspond + to `C1' and `C2' in the paper. */ + real X_C1_det, C0_X_det, C0_C1_det; + real alpha1, alpha2; + spline_type spline; + vector_type start_vector, end_vector; + unsigned i; + vector_type t1_hat = *CURVE_START_TANGENT (curve); + vector_type t2_hat = *CURVE_END_TANGENT (curve); + real C[2][2] = { { 0.0, 0.0 }, { 0.0, 0.0 } }; + real X[2] = { 0.0, 0.0 }; + int Alen = CURVE_LENGTH (curve); + vector_type *A; + + A = g_new0 (vector_type, Alen * 2); + + START_POINT (spline) = CURVE_POINT (curve, 0); + END_POINT (spline) = LAST_CURVE_POINT (curve); + SPLINE_LINEARITY (spline) = 0; + start_vector = make_vector (START_POINT (spline)); + end_vector = make_vector (END_POINT (spline)); + + for (i = 0; i < CURVE_LENGTH (curve); i++) + { + A[i*2+0] = Vmult_scalar (t1_hat, B1 (U (i))); + A[i*2+1] = Vmult_scalar (t2_hat, B2 (U (i))); + } + + for (i = 0; i < CURVE_LENGTH (curve); i++) + { + vector_type temp, temp0, temp1, temp2, temp3; + vector_type *Ai = &A[i*2]; + + 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 (U (i))); + temp1 = Vmult_scalar (start_vector, B1 (U (i))); + temp2 = Vmult_scalar (end_vector, B2 (U (i))); + temp3 = Vmult_scalar (end_vector, B3 (U (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]); + } + + 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) + FATAL ("zero determinant of C0*C1"); + + 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) = CUBIC; + + g_free (A); + + return spline; +} + +/* Find closer-to-optimal t values for the given x,y's and control + points, using Newton-Raphson iteration. A good description of this + is in Plass & Stone. This routine performs one step in the + iteration, not the whole thing. */ + +static boolean +reparameterize (curve_type curve, spline_type S) +{ + unsigned p, i; + spline_type S1, S2; /* S' and S''. */ + +/* REPORT ("r"); */ + + /* Find the first and second derivatives of S. To make + `evaluate_spline' work, we fill the beginning points (i.e., the first + two for a linear spline and the first three for a quadratic one), + even though this is at odds with the rest of the program. */ + for (i = 0; i < 3; i++) + { + S1.v[i].x = 3.0 * (S.v[i + 1].x - S.v[i].x); + S1.v[i].y = 3.0 * (S.v[i + 1].y - S.v[i].y); + } + S1.v[i].x = S1.v[i].y = -1.0; /* These will never be accessed. */ + SPLINE_DEGREE (S1) = QUADRATIC; + + for (i = 0; i < 2; i++) + { + S2.v[i].x = 2.0 * (S1.v[i + 1].x - S1.v[i].x); + S2.v[i].y = 2.0 * (S1.v[i + 1].y - S1.v[i].y); + } + S2.v[2].x = S2.v[2].y = S2.v[3].x = S2.v[3].y = -1.0; + SPLINE_DEGREE (S2) = LINEAR; + + for (p = 0; p < CURVE_LENGTH (curve); p++) + { + real new_distance, old_distance; + real_coordinate_type new_point; + real_coordinate_type point = CURVE_POINT (curve, p); + real t = CURVE_T (curve, p); + + /* Find the points at this t on S, S', and S''. */ + real_coordinate_type S_t = evaluate_spline (S, t); + real_coordinate_type S1_t = evaluate_spline (S1, t); + real_coordinate_type S2_t = evaluate_spline (S2, t); + + /* The differences in x and y (Q1(t) on p.62 of Schneider's thesis). */ + real_coordinate_type d; + real numerator; + real denominator; + + d.x = S_t.x - point.x; + d.y = S_t.y - point.y; + + /* The step size, f(t)/f'(t). */ + numerator = d.x * S1_t.x + d.y * S1_t.y; + denominator = (SQUARE (S1_t.x) + d.x * S2_t.x + + SQUARE (S1_t.y) + d.y * S2_t.y); + + /* We compute the distances only to be able to check that we + really are moving closer. I don't know how this condition can + be reliably checked for in advance, but I know what it means in + practice: the original fit was not good enough for the + Newton-Raphson iteration to converge. Therefore, we need to + abort the reparameterization, and subdivide. */ + old_distance = distance (S_t, point); + CURVE_T (curve, p) -= numerator / denominator; + new_point = evaluate_spline (S, CURVE_T (curve, p)); + new_distance = distance (new_point, point); + + if (new_distance > old_distance) + { +/* REPORT ("!"); */ +/* LOG3 (" Stopped reparameterizing; %.3f > %.3f at point %u.\n", */ +/* new_distance, old_distance, p); */ + return false; + } + + /* The t value might be negative or > 1, if the choice of control + points wasn't so great. We have no difficulty in evaluating + a spline at a t outside the range zero to one, though, so it + doesn't matter. (Although it is a little unconventional.) */ + } +/* LOG (" reparameterized curve:\n "); */ +/* log_curve (curve, true); */ + + return true; +} + +/* This routine finds the best place to subdivide the curve CURVE, + somewhere near to the point whose index is INITIAL. Originally, + curves were always subdivided at the point of worst error, which is + intuitively appealing, but doesn't always give the best results. For + example, at the end of a serif that tapers into the stem, the best + subdivision point is at the point where they join, even if the worst + point is a little ways into the serif. + + We return the index of the point at which to subdivide. */ + +static unsigned +find_subdivision (curve_type curve, unsigned initial) +{ + unsigned i, n_done; + int best_point = -1; + vector_type best = { FLT_MAX, FLT_MAX }; + unsigned search = subdivide_search * CURVE_LENGTH (curve); + +/* LOG1 (" Number of points to search: %u.\n", search); */ + + /* We don't want to find the first (or last) point in the curve. */ + for (i = initial, n_done = 0; i > 0 && n_done < search; + i = CURVE_PREV (curve, i), n_done++) + { + if (test_subdivision_point (curve, i, &best)) + { + best_point = i; +/* LOG3 (" Better point: (%.3f,%.3f), #%u.\n", */ +/* CURVE_POINT (curve, i).x, CURVE_POINT (curve, i).y, i); */ + } + } + + /* If we found a good one, let's take it. */ + if (best_point != -1) + return best_point; + + for (i = CURVE_NEXT (curve, initial), n_done = 0; + i < CURVE_LENGTH (curve) - 1 && n_done < search; + i = CURVE_NEXT (curve, i), n_done++) + { + if (test_subdivision_point (curve, i, &best)) + { + best_point = i; +/* LOG3 (" Better point at (%.3f,%.3f), #%u.\n", */ +/* CURVE_POINT (curve, i).x, CURVE_POINT (curve, i).y, i); */ + } + } + + /* If we didn't find any better point, return the original. */ + return best_point == -1 ? initial : best_point; +} + + +/* Here are some macros that decide whether or not we're at a + ``join point'', as described above. */ +#define ONLY_ONE_LESS(v) \ + (((v).dx < subdivide_threshold && (v).dy > subdivide_threshold) \ + || ((v).dy < subdivide_threshold && (v).dx > subdivide_threshold)) + +#define BOTH_GREATER(v) \ + ((v).dx > subdivide_threshold && (v).dy > subdivide_threshold) + +/* We assume that the vectors V1 and V2 are nonnegative. */ +#define JOIN(v1, v2) \ + ((ONLY_ONE_LESS (v1) && BOTH_GREATER (v2)) \ + || (ONLY_ONE_LESS (v2) && BOTH_GREATER (v1))) + +/* If the component D of the vector V is smaller than the best so far, + update the best point. */ +#define UPDATE_BEST(v, d) \ + do \ + { \ + if ((v).d < subdivide_threshold && (v).d < best->d) \ + best->d = (v).d; \ + } \ + while (0) + + +/* If the point INDEX in the curve CURVE is the best subdivision point + we've found so far, update the vector BEST. */ + +static boolean +test_subdivision_point (curve_type curve, unsigned index, vector_type *best) +{ + unsigned count; + vector_type in, out; + boolean join = false; + + find_curve_vectors (index, curve, subdivide_surround, &in, &out, &count); + + /* We don't want to subdivide at points which are very close to the + endpoints, so if the vectors aren't computed from as many points as + we asked for, don't bother checking this point. */ + if (count == subdivide_surround) + { + in = Vabs (in); + out = Vabs (out); + + join = JOIN (in, out); + if (join) + { + UPDATE_BEST (in, dx); + UPDATE_BEST (in, dy); + UPDATE_BEST (out, dx); + UPDATE_BEST (out, dy); + } + } + + return join; +} + +/* 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++) + { + real_coordinate_type point = CURVE_POINT (curve, p), + previous_p = CURVE_POINT (curve, p - 1); + real d = distance (point, previous_p); + CURVE_T (curve, p) = CURVE_T (curve, p - 1) + d; + } + + assert (LAST_CURVE_T (curve) != 0.0); + + for (p = 1; p < CURVE_LENGTH (curve); p++) + CURVE_T (curve, p) = CURVE_T (curve, p) / LAST_CURVE_T (curve); + +/* 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, boolean to_start_point, boolean cross_curve) +{ + vector_type tangent; + vector_type **curve_tangent = to_start_point ? &(CURVE_START_TANGENT (curve)) + : &(CURVE_END_TANGENT (curve)); + unsigned n_points = 0; + +/* LOG1 (" tangent to %s: ", to_start_point ? "start" : "end"); */ + + if (*curve_tangent == NULL) + { + *curve_tangent = g_new (vector_type, 1); + tangent = find_half_tangent (curve, to_start_point, &n_points); + + if (cross_curve) + { + curve_type adjacent_curve + = to_start_point ? PREVIOUS_CURVE (curve) : NEXT_CURVE (curve); + vector_type tangent2 + = find_half_tangent (adjacent_curve, !to_start_point, &n_points); + +/* LOG2 ("(adjacent curve half tangent (%.3f,%.3f)) ", */ +/* tangent2.dx, tangent2.dy); */ + tangent = Vadd (tangent, tangent2); + } + + assert (n_points > 0); + **curve_tangent = Vmult_scalar (tangent, 1.0 / n_points); + } + else + { +/* LOG ("(already computed) "); */ + } + +/* LOG2 ("(%.3f,%.3f).\n", (*curve_tangent)->dx, (*curve_tangent)->dy); */ +} + + +/* 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, boolean to_start_point, unsigned *n_points) +{ + unsigned p; + int factor = to_start_point ? 1 : -1; + unsigned tangent_index = to_start_point ? 0 : c->length - 1; + real_coordinate_type tangent_point = CURVE_POINT (c, tangent_index); + vector_type tangent; + + tangent.dx = 0.0; + tangent.dy = 0.0; + + for (p = 1; p <= tangent_surround; p++) + { + int this_index = p * factor + tangent_index; + real_coordinate_type this_point; + + if (this_index < 0 || this_index >= 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), + 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 real +find_error (curve_type curve, spline_type spline, unsigned *worst_point) +{ + unsigned this_point; + real total_error = 0.0; + real worst_error = FLT_MIN; + + *worst_point = CURVE_LENGTH (curve) + 1; /* A sentinel value. */ + + for (this_point = 0; this_point < CURVE_LENGTH (curve); this_point++) + { + real_coordinate_type curve_point = CURVE_POINT (curve, this_point); + real t = CURVE_T (curve, this_point); + real_coordinate_type spline_point = evaluate_spline (spline, t); + real 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 + printf ("No worst point found; something is wrong"); + } + else + { +/* LOG4 (" Worst error (at (%.3f,%.3f), point #%u) was %.3f.\n", */ +/* CURVE_POINT (curve, *worst_point).x, */ +/* CURVE_POINT (curve, *worst_point).y, *worst_point, worst_error); */ +/* LOG1 (" Total error was %.3f.\n", total_error); */ +/* LOG2 (" 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 boolean +spline_linear_enough (spline_type *spline, curve_type curve) +{ + real A, B, C, slope; + unsigned this_point; + real distance = 0.0; + +/* LOG ("Checking linearity:\n"); */ + + /* For a line described by Ax + By + C = 0, the distance d from a + point (x0,y0) to that line is: + + d = | Ax0 + By0 + C | / sqrt (A^2 + B^2). + + We can find A, B, and C from the starting and ending points, + unless the line defined by those two points is vertical. */ + + if (epsilon_equal (START_POINT (*spline).x, END_POINT (*spline).x)) + { + A = 1; + B = 0; + C = -START_POINT (*spline).x; + } + else + { + /* Plug the spline's starting and ending points into the two-point + equation for a line, that is, + + (y - y1) = ((y2 - y1)/(x2 - x1)) * (x - x1) + + to get the values for A, B, and C. */ + + slope = ((END_POINT (*spline).y - START_POINT (*spline).y) + / (END_POINT (*spline).x - START_POINT (*spline).x)); + A = -slope; + B = 1; + C = slope * START_POINT (*spline).x - START_POINT (*spline).y; + } +/* LOG3 (" Line is %.3fx + %.3fy + %.3f = 0.\n", A, B, C); */ + + for (this_point = 0; this_point < CURVE_LENGTH (curve); this_point++) + { + real t = CURVE_T (curve, this_point); + real_coordinate_type spline_point = evaluate_spline (*spline, t); + + distance += fabs (A * spline_point.x + B * spline_point.y + C) + / sqrt (A * A + B * B); + } +/* LOG1 (" Total distance is %.3f, ", distance); */ + + distance /= CURVE_LENGTH (curve); +/* LOG1 ("which is %.3f normalized.\n", distance); */ + + /* 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) = distance / CURVE_LENGTH (curve); +/* LOG1 (" Final linearity: %.3f.\n", SPLINE_LINEARITY (*spline)); */ + + return distance < line_threshold; +} + + +/* 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) +{ + unsigned this_spline; + boolean found_cubic = false; + unsigned length = SPLINE_LIST_LENGTH (*spline_list); + +/* LOG1 ("\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)) == CUBIC) + { + 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) == LINEAR) + { +/* LOG1 (" #%u: ", this_spline); */ + if (SPLINE_LINEARITY (s) > line_reversion_threshold) + { +/* LOG ("reverted, "); */ + SPLINE_DEGREE (SPLINE_LIST_ELT (*spline_list, this_spline)) + = CUBIC; + } +/* LOG1 ("linearity %.3f.\n", SPLINE_LINEARITY (s)); */ + } + } + else + { +/* LOG (" No lines.\n"); */ + } +} + +/* When we have finished fitting an entire pixel outline to a spline + list L, we go through L to ensure that any endpoints that are ``close + enough'' (i.e., within `align_threshold') to being the same really + are the same. */ + +/* This macro adjusts the AXIS axis on the starting and ending points on + a particular spline if necessary. */ +#define TRY_AXIS(axis) \ + do \ + { \ + real delta = fabs (end.axis - start.axis); \ + \ + if (!epsilon_equal (delta, 0.0) && delta <= align_threshold) \ + { \ + spline_type *next = &NEXT_SPLINE_LIST_ELT (*l, this_spline); \ + spline_type *prev = &PREV_SPLINE_LIST_ELT (*l, this_spline); \ + \ + START_POINT (*s).axis = END_POINT (*s).axis \ + = END_POINT (*prev).axis = START_POINT (*next).axis \ + = (start.axis + end.axis) / 2; \ + spline_change = true; \ + } \ + } \ + while (0) + +static void +align (spline_list_type *l) +{ + boolean change; + unsigned this_spline; + unsigned length = SPLINE_LIST_LENGTH (*l); + +/* LOG1 ("\nAligning spline list (length %u):\n", length); */ + + do + { + change = false; + +/* LOG (" "); */ + + for (this_spline = 0; this_spline < length; this_spline++) + { + boolean spline_change = false; + spline_type *s = &SPLINE_LIST_ELT (*l, this_spline); + real_coordinate_type start = START_POINT (*s); + real_coordinate_type end = END_POINT (*s); + + TRY_AXIS (x); + TRY_AXIS (y); + if (spline_change) + { +/* LOG1 ("%u ", this_spline); */ + change |= spline_change; + } + } +/* LOG ("\n"); */ + } + while (change); +} + +/* 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) + { + g_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)++; + list->data = (unsigned *)g_realloc(list->data,(INDEX_LIST_LENGTH (*list)) * sizeof(unsigned)); +/* XRETALLOC (list->data, INDEX_LIST_LENGTH (*list), unsigned); */ + list->data[INDEX_LIST_LENGTH (*list) - 1] = new_index; +} |