1 /*
2 * curvefit.cpp
3 * scorealign
4 *
5 * Created by Roger Dannenberg on 10/20/07.
6 *
7 */
8
9 #include <assert.h>
10 #include <stdlib.h>
11 #include <string.h>
12 #include <math.h>
13 #include "comp_chroma.h"
14 #include "sautils.h"
15 // the following are needed to get Scorealign
16 #include <fstream>
17 #include "allegro.h"
18 #include "audioreader.h"
19 #include "scorealign.h"
20 #include "hillclimb.h"
21 #include "curvefit.h"
22
23 void save_path(char *filename);
24
25 /* Curvefit class: do hill-climbing to fit lines to data
26 *
27 * This class implements the algorithm described above.
28 * The problem is partitioned into the general search algorithm
29 * (implemented in Hillclimb::optimize) and the evaluation function
30 * (implemented in Curvefit::evaluate). A brute-force evaluation
31 * would simply recompute the cost of the entire path every time,
32 * but note that the search algorithm works by adjusting one parameter
33 * at a time. This affects at most two line segments, so the rest
34 * contribute a cost that does not need to be recomputed. Thus the
35 * total cost can be computed incrementally. It is hard to see how
36 * to use this optimization within the general Hillclimb:optimize
37 * method, so to avoid making that algorithm very specific and ugly,
38 * I decided to hide the incremental nature of evaluate inside
39 * the evaluate function itself. The way this works is that evaluate
40 * keeps a cache of the coordinates of each line segment and the
41 * resulting cost of the segment. Before recomputing any segment,
42 * the cache is consulted. If the end points have not moved, the
43 * cached value is retrieved. Ideally, there should be a 3-element
44 * cache because endpoints are moved and then restored. (The three
45 * elements would hold the results of the original, changed left,
46 * and changed right endpoints.) The bigger cache would eliminate
47 * 1/3 of the computation, but the simple cache already eliminates
48 * about (n-2)/n of the work, so that should help a lot.
49 */
50
51 class Curvefit : public Hillclimb {
52 public:
Curvefit(Scorealign * sa_,bool verbose_)53 Curvefit(Scorealign *sa_, bool verbose_) {
54 sa = sa_;
55 verbose = verbose_;
56 p1_cache = p2_cache = d_cache = x = NULL;
57 }
58 ~Curvefit();
59 virtual double evaluate();
60 void setup(int n);
61 void set_step_size(double ss);
get_x()62 double *get_x() { return x; }
63 private:
64 Scorealign *sa;
65 bool verbose;
66 double line_dist(int i); // get cost of line segment i
67 double compute_dist(int i); // compute cost of line segment i
68 double distance_rc(int row, int col);
69 double distance_xy(double x, double y);
70
71 double *p1_cache; // left endpoint y values
72 double *p2_cache; // right endpoint y values
73 double *d_cache; // cached cost of line segment
74 double *x; // the x values of line segment endpoints
75 // (the y values are in parameters[])
76 };
77
78
evaluate()79 double Curvefit::evaluate()
80 {
81 double sum = 0;
82 // why does this loop go to n-2? Because i represents the left endpoint
83 // of the line segment. There are n parameters, but only n-1 segments.
84 for (int i = 0; i < n-1; i++) {
85 sum += line_dist(i); // look up in cache or recompute each segment
86 }
87 return -sum; // return negative of distance so that bigger will be better
88 }
89
90
line_dist(int i)91 double Curvefit::line_dist(int i)
92 {
93 if (p1_cache[i] == parameters[i] &&
94 p2_cache[i] == parameters[i+1]) {
95 // endpoints have not changed:
96 return d_cache[i];
97 }
98 // otherwise, we need to recompute and save dist in cache
99 double d = compute_dist(i);
100 p1_cache[i] = parameters[i];
101 p2_cache[i] = parameters[i+1];
102 d_cache[i] = d;
103 return d;
104 }
105
106
setup(int segments)107 void Curvefit::setup(int segments)
108 {
109 // number of parameters is greater than segments because the left
110 // col of segment i is parameter i, so the right col of
111 // the last segment == parameter[segments].
112 Hillclimb::setup(segments + 1);
113 p1_cache = ALLOC(double, n);
114 p2_cache = ALLOC(double, n);
115 d_cache = ALLOC(double, n);
116 x = ALLOC(double, n);
117 int i;
118 // ideal frames per segment
119 float seg_length = ((float) (sa->last_x - sa->first_x)) / segments;
120 for (i = 0; i < n; i++) { // initialize cache keys to garbage
121 p1_cache[i] = p2_cache[i] = -999999.99;
122 // initialize x values
123 x[i] = ROUND(sa->first_x + i * seg_length);
124 // now initialize parameters based on pathx/pathy/time_map
125 // time_map has y values for each x
126 parameters[i] = sa->time_map[(int) x[i]];
127 assert(parameters[i] >= 0);
128 if (verbose)
129 printf("initial x[%d] = %g, parameters[%d] = %g\n",
130 i, x[i], i, parameters[i]);
131 step_size[i] = 0.5;
132 min_param[i] = 0;
133 max_param[i] = sa->last_y;
134 }
135 }
136
137
~Curvefit()138 Curvefit::~Curvefit()
139 {
140 if (p1_cache) FREE(p1_cache);
141 if (p2_cache) FREE(p2_cache);
142 if (d_cache) FREE(d_cache);
143 if (x) FREE(x);
144 }
145
146
147 // distance_rc -- look up or compute distance between chroma vectors
148 // at row, col in similarity matrix
149 //
150 // Note: in current implementation, there is no stored representation
151 // of the matrix, so we have to recompute every time. It would be
152 // possible to store the whole matrix, but it's large and it would
153 // double the memory requirements (we already allocate the large
154 // PATH array in compare_chroma to compute the optimal path.
155 //
156 // Since distance can be computed relatively quickly, a better plan
157 // would be to cache values along the path. Here's a brief design
158 // (for the future, assuming this routine is actually a hot spot):
159 // Allocate a matrix that is, say, 20 x file0_frames to contain distances
160 // that are +/- 10 frames from the path. Initialize cells to -1.
161 // Allocate an array of integer offsets of size file1_frames.
162 // Fill in the integer offsets with the column number (pathy) value of
163 // the path.
164 // Now, to get distance_rc(row, col):
165 // offset = offsets[row]
166 // i = 10 + col - offset;
167 // if (i < 0 || i > 20) /* not in cache */ return compute_distance(...);
168 // dist = distances[20 * row + i];
169 // if (dist == -1) { return distances[20 * row + i] = compute_distance...}
170 // return dist;
171 //
distance_rc(int row,int col)172 double Curvefit::distance_rc(int row, int col)
173 {
174 double dist = sa->gen_dist(row, col);
175 if (dist > 20) // DEBUGGING
176 printf("internal error");
177 return dist;
178 }
179
180
181 // compute distance from distance matrix using interpolation. A least
182 // one of x, y should be an integer value so interpolation is only 2-way
distance_xy(double x,double y)183 double Curvefit::distance_xy(double x, double y)
184 {
185 int xi = (int) x;
186 int yi = (int) y;
187 if (xi == x) { // x is integer, interpolate along y axis
188 double d1 = distance_rc(xi, yi);
189 double d2 = distance_rc(xi, yi + 1);
190 return interpolate(yi, d1, yi + 1, d2, y);
191 } else if (yi == y) { // y is integer, interpolate along x axis
192 double d1 = distance_rc(xi, yi);
193 double d2 = distance_rc(xi + 1, yi);
194 return interpolate(xi, d1, xi + 1, d2, x);
195 } else {
196 printf("FATAL INTERNAL ERROR IN distance_xy: neither x nor y is "
197 "an integer\n");
198 assert(false);
199 }
200 }
201
202
compute_dist(int i)203 double Curvefit::compute_dist(int i)
204 {
205 double x1 = x[i], x2 = x[i+1];
206 double y1 = parameters[i], y2 = parameters[i+1];
207 double dx = x2 - x1, dy = y2 - y1;
208 double sum = 0;
209 int n;
210 assert(x1 >= 0 && x2 >= 0 && y1 >= 0 && y2 >= 0);
211 if (dx > dy) { // evauate at each x
212 n = (int) dx;
213 for (int x = (int) x1; x < x2; x++) {
214 double y = interpolate(x1, y1, x2, y2, x);
215 sum += distance_xy(x, y);
216 }
217 } else { // evaluate at each y
218 n = (int) dy;
219 for (int y = (int) y1; y < y2; y++) {
220 double x = interpolate(y1, x1, y2, x2, y);
221 sum += distance_xy(x, y);
222 // printf("dist %g %d = %g\n", x, y, distance_xy(x, y));
223 }
224 }
225 // normalize using line length: sum/n is average distance. Multiply
226 // avg. distance (cost per unit length) by length to get total cost.
227 // Note: this gives an advantage to direct diagonal paths without bends
228 // because longer path lengths result in higher total cost. This also
229 // gives heigher weight to longer segments, although all segments are
230 // about the same length.
231 double rslt = sqrt(dx*dx + dy*dy) * sum / n;
232 // printf("compute_dist %d: x1 %g y1 %g x2 %g y2 %g sum %g rslt %g\n",
233 // i, x1, y1, x2, y2, sum, rslt);
234 if (rslt < 0 || rslt > 24 * n) { // DEBUGGING
235 printf("internal error");
236 }
237 return rslt;
238 }
239
240
set_step_size(double ss)241 void Curvefit::set_step_size(double ss)
242 {
243 for (int i = 0; i < n; i++) {
244 step_size[i] = ss;
245 }
246 }
247
248
249 static long curvefit_iterations;
250
251 // This is a callback from Hillclimb::optimize to report progress
252 // We can't know percentage completion because we don't know how
253 // many iterations it will take to converge, so we just report
254 // iterations. The SAProgress class assumes some number based
255 // on experience.
256 //
257 // Normally, the iterations parameter is a good indicator of work
258 // expended so far, but since we call Hillclimb::optimize twice
259 // (second time with a finer grid to search), ignore iterations
260 // and use curvefit_iterations, a global counter, instead. This
261 // assumes that curvefit_progress is called once for each iteration.
262 //
curvefit_progress(void * cookie,int iterations,double best)263 void curvefit_progress(void *cookie, int iterations, double best)
264 {
265 Scorealign *sa = (Scorealign *) cookie;
266 if (sa->progress) {
267 sa->progress->set_smoothing_progress(++curvefit_iterations);
268 }
269 }
270
271
curve_fitting(Scorealign * sa,bool verbose)272 void curve_fitting(Scorealign *sa, bool verbose)
273 {
274 if (verbose)
275 printf("Performing line-segment approximation with %gs segments.\n",
276 sa->line_time);
277 Curvefit curvefit(sa, verbose);
278 double *parameters;
279 double *x;
280 curvefit_iterations = 0;
281 // how many segments? About total time / line_time:
282 int segments =
283 (int) (0.5 + (sa->actual_frame_period_0 * (sa->last_x - sa->first_x)) /
284 sa->line_time);
285 curvefit.setup(segments);
286 curvefit.optimize(&curvefit_progress, sa);
287 // further optimization with smaller step sizes:
288 // this step size will interpolate 0.25s frames down to 10ms
289 curvefit.set_step_size(0.04);
290 curvefit.optimize(&curvefit_progress, sa);
291 parameters = curvefit.get_parameters();
292 x = curvefit.get_x();
293 // now, rewrite pathx and pathy according to segments
294 // pathx and pathy are generously allocated, so we can change pathlen
295 // each segment goes from x[i], parameters[i] to x[i+1], parameters[i+1]
296 int i;
297 int j = 0; // index into path
298 for (i = 0; i < segments; i++) {
299 int x1 = (int) x[i];
300 int x2 = (int) x[i+1];
301 int y1 = (int) parameters[i];
302 int y2 = (int) parameters[i+1];
303 int dx = x2 - x1;
304 int dy = y2 - y1;
305 if (dx >= dy) { // output point at each x
306 int x;
307 for (x = x1; x < x2; x++) {
308 sa->pathx[j] = x;
309 sa->pathy[j] = (int) (0.5 + interpolate(x1, y1, x2, y2, x));
310 j++;
311 }
312 } else {
313 int y;
314 for (y = y1; y < y2; y++) {
315 sa->pathx[j] = (int) (0.5 + interpolate(y1, x1, y2, x2, y));
316 sa->pathy[j] = y;
317 j++;
318 }
319 }
320 }
321 // output last point
322 sa->pathx[j] = (int) x[segments];
323 sa->pathy[j] = (int) (0.5 + parameters[segments]);
324 j++;
325 sa->set_pathlen(j);
326 }
327
328
329
330