Cleanup: style, use braces for blenkernel
[blender.git] / source / blender / blenkernel / intern / tracking_stabilize.c
1 /*
2  * This program is free software; you can redistribute it and/or
3  * modify it under the terms of the GNU General Public License
4  * as published by the Free Software Foundation; either version 2
5  * of the License, or (at your option) any later version.
6  *
7  * This program is distributed in the hope that it will be useful,
8  * but WITHOUT ANY WARRANTY; without even the implied warranty of
9  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
10  * GNU General Public License for more details.
11  *
12  * You should have received a copy of the GNU General Public License
13  * along with this program; if not, write to the Free Software Foundation,
14  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
15  *
16  * The Original Code is Copyright (C) 2011 Blender Foundation.
17  * All rights reserved.
18  */
19
20 /** \file
21  * \ingroup bke
22  *
23  * This file contains implementation of 2D image stabilization.
24  */
25
26 #include <limits.h>
27
28 #include "DNA_movieclip_types.h"
29 #include "DNA_scene_types.h"
30 #include "DNA_anim_types.h"
31 #include "RNA_access.h"
32
33 #include "BLI_utildefines.h"
34 #include "BLI_sort_utils.h"
35 #include "BLI_ghash.h"
36 #include "BLI_math_vector.h"
37 #include "BLI_math.h"
38 #include "BLI_task.h"
39
40 #include "BKE_tracking.h"
41 #include "BKE_movieclip.h"
42 #include "BKE_fcurve.h"
43
44 #include "MEM_guardedalloc.h"
45 #include "IMB_imbuf_types.h"
46 #include "IMB_imbuf.h"
47
48 /* == Parameterization constants == */
49
50 /* When measuring the scale changes relative to the rotation pivot point, it
51  * might happen accidentally that a probe point (tracking point), which doesn't
52  * actually move on a circular path, gets very close to the pivot point, causing
53  * the measured scale contribution to go toward infinity. We damp this undesired
54  * effect by adding a bias (floor) to the measured distances, which will
55  * dominate very small distances and thus cause the corresponding track's
56  * contribution to diminish.
57  * Measurements happen in normalized (0...1) coordinates within a frame.
58  */
59 static float SCALE_ERROR_LIMIT_BIAS = 0.01f;
60
61 /* When to consider a track as completely faded out.
62  * This is used in conjunction with the "disabled" flag of the track
63  * to determine start positions, end positions and gaps
64  */
65 static float EPSILON_WEIGHT = 0.005f;
66
67 /* == private working data == */
68
69 /* Per track baseline for stabilization, defined at reference frame.
70  * A track's reference frame is chosen as close as possible to the (global)
71  * anchor_frame. Baseline holds the constant part of each track's contribution
72  * to the observed movement; it is calculated at initialization pass, using the
73  * measurement value at reference frame plus the average contribution to fill
74  * the gap between global anchor_frame and the reference frame for this track.
75  * This struct with private working data is associated to the local call context
76  * via `StabContext::private_track_data`
77  */
78 typedef struct TrackStabilizationBase {
79   float stabilization_offset_base[2];
80
81   /* measured relative to translated pivot */
82   float stabilization_rotation_base[2][2];
83
84   /* measured relative to translated pivot */
85   float stabilization_scale_base;
86
87   bool is_init_for_stabilization;
88   FCurve *track_weight_curve;
89 } TrackStabilizationBase;
90
91 /* Tracks are reordered for initialization, starting as close as possible to
92  * anchor_frame
93  */
94 typedef struct TrackInitOrder {
95   int sort_value;
96   int reference_frame;
97   MovieTrackingTrack *data;
98 } TrackInitOrder;
99
100 /* Per frame private working data, for accessing possibly animated values. */
101 typedef struct StabContext {
102   MovieClip *clip;
103   MovieTracking *tracking;
104   MovieTrackingStabilization *stab;
105   GHash *private_track_data;
106   FCurve *locinf;
107   FCurve *rotinf;
108   FCurve *scaleinf;
109   FCurve *target_pos[2];
110   FCurve *target_rot;
111   FCurve *target_scale;
112   bool use_animation;
113 } StabContext;
114
115 static TrackStabilizationBase *access_stabilization_baseline_data(StabContext *ctx,
116                                                                   MovieTrackingTrack *track)
117 {
118   return BLI_ghash_lookup(ctx->private_track_data, track);
119 }
120
121 static void attach_stabilization_baseline_data(StabContext *ctx,
122                                                MovieTrackingTrack *track,
123                                                TrackStabilizationBase *private_data)
124 {
125   BLI_ghash_insert(ctx->private_track_data, track, private_data);
126 }
127
128 static void discard_stabilization_baseline_data(void *val)
129 {
130   if (val != NULL) {
131     MEM_freeN(val);
132   }
133 }
134
135 /* == access animated values for given frame == */
136
137 static FCurve *retrieve_stab_animation(MovieClip *clip, const char *data_path, int idx)
138 {
139   return id_data_find_fcurve(&clip->id,
140                              &clip->tracking.stabilization,
141                              &RNA_MovieTrackingStabilization,
142                              data_path,
143                              idx,
144                              NULL);
145 }
146
147 static FCurve *retrieve_track_weight_animation(MovieClip *clip, MovieTrackingTrack *track)
148 {
149   return id_data_find_fcurve(&clip->id, track, &RNA_MovieTrackingTrack, "weight_stab", 0, NULL);
150 }
151
152 static float fetch_from_fcurve(FCurve *animationCurve,
153                                int framenr,
154                                StabContext *ctx,
155                                float default_value)
156 {
157   if (ctx && ctx->use_animation && animationCurve) {
158     int scene_framenr = BKE_movieclip_remap_clip_to_scene_frame(ctx->clip, framenr);
159     return evaluate_fcurve(animationCurve, scene_framenr);
160   }
161   return default_value;
162 }
163
164 static float get_animated_locinf(StabContext *ctx, int framenr)
165 {
166   return fetch_from_fcurve(ctx->locinf, framenr, ctx, ctx->stab->locinf);
167 }
168
169 static float get_animated_rotinf(StabContext *ctx, int framenr)
170 {
171   return fetch_from_fcurve(ctx->rotinf, framenr, ctx, ctx->stab->rotinf);
172 }
173
174 static float get_animated_scaleinf(StabContext *ctx, int framenr)
175 {
176   return fetch_from_fcurve(ctx->scaleinf, framenr, ctx, ctx->stab->scaleinf);
177 }
178
179 static void get_animated_target_pos(StabContext *ctx, int framenr, float target_pos[2])
180 {
181   target_pos[0] = fetch_from_fcurve(ctx->target_pos[0], framenr, ctx, ctx->stab->target_pos[0]);
182   target_pos[1] = fetch_from_fcurve(ctx->target_pos[1], framenr, ctx, ctx->stab->target_pos[1]);
183 }
184
185 static float get_animated_target_rot(StabContext *ctx, int framenr)
186 {
187   return fetch_from_fcurve(ctx->target_rot, framenr, ctx, ctx->stab->target_rot);
188 }
189
190 static float get_animated_target_scale(StabContext *ctx, int framenr)
191 {
192   return fetch_from_fcurve(ctx->target_scale, framenr, ctx, ctx->stab->scale);
193 }
194
195 static float get_animated_weight(StabContext *ctx, MovieTrackingTrack *track, int framenr)
196 {
197   TrackStabilizationBase *working_data = access_stabilization_baseline_data(ctx, track);
198   if (working_data && working_data->track_weight_curve) {
199     int scene_framenr = BKE_movieclip_remap_clip_to_scene_frame(ctx->clip, framenr);
200     return evaluate_fcurve(working_data->track_weight_curve, scene_framenr);
201   }
202   /* Use weight at global 'current frame' as fallback default. */
203   return track->weight_stab;
204 }
205
206 static void use_values_from_fcurves(StabContext *ctx, bool toggle)
207 {
208   if (ctx != NULL) {
209     ctx->use_animation = toggle;
210   }
211 }
212
213 /* Prepare per call private working area.
214  * Used for access to possibly animated values: retrieve available F-curves.
215  */
216 static StabContext *initialize_stabilization_working_context(MovieClip *clip)
217 {
218   StabContext *ctx = MEM_callocN(sizeof(StabContext), "2D stabilization animation runtime data");
219   ctx->clip = clip;
220   ctx->tracking = &clip->tracking;
221   ctx->stab = &clip->tracking.stabilization;
222   ctx->private_track_data = BLI_ghash_ptr_new("2D stabilization per track private working data");
223   ctx->locinf = retrieve_stab_animation(clip, "influence_location", 0);
224   ctx->rotinf = retrieve_stab_animation(clip, "influence_rotation", 0);
225   ctx->scaleinf = retrieve_stab_animation(clip, "influence_scale", 0);
226   ctx->target_pos[0] = retrieve_stab_animation(clip, "target_pos", 0);
227   ctx->target_pos[1] = retrieve_stab_animation(clip, "target_pos", 1);
228   ctx->target_rot = retrieve_stab_animation(clip, "target_rot", 0);
229   ctx->target_scale = retrieve_stab_animation(clip, "target_zoom", 0);
230   ctx->use_animation = true;
231   return ctx;
232 }
233
234 /**
235  * Discard all private working data attached to this call context.
236  *
237  * \note We allocate the record for the per track baseline contribution
238  * locally for each call context (i.e. call to #BKE_tracking_stabilization_data_get)
239  * Thus it is correct to discard all allocations found within the
240  * corresponding _local_ GHash.
241  */
242 static void discard_stabilization_working_context(StabContext *ctx)
243 {
244   if (ctx != NULL) {
245     BLI_ghash_free(ctx->private_track_data, NULL, discard_stabilization_baseline_data);
246     MEM_freeN(ctx);
247   }
248 }
249
250 static bool is_init_for_stabilization(StabContext *ctx, MovieTrackingTrack *track)
251 {
252   TrackStabilizationBase *working_data = access_stabilization_baseline_data(ctx, track);
253   return (working_data != NULL && working_data->is_init_for_stabilization);
254 }
255
256 static bool is_usable_for_stabilization(StabContext *ctx, MovieTrackingTrack *track)
257 {
258   return (track->flag & TRACK_USE_2D_STAB) && is_init_for_stabilization(ctx, track);
259 }
260
261 static bool is_effectively_disabled(StabContext *ctx,
262                                     MovieTrackingTrack *track,
263                                     MovieTrackingMarker *marker)
264 {
265   return (marker->flag & MARKER_DISABLED) ||
266          (EPSILON_WEIGHT > get_animated_weight(ctx, track, marker->framenr));
267 }
268
269 static int search_closest_marker_index(MovieTrackingTrack *track, int ref_frame)
270 {
271   MovieTrackingMarker *markers = track->markers;
272   int end = track->markersnr;
273   int i = track->last_marker;
274
275   i = MAX2(0, i);
276   i = MIN2(i, end - 1);
277   for (; i < end - 1 && markers[i].framenr <= ref_frame; ++i) {
278     ;
279   }
280   for (; 0 < i && markers[i].framenr > ref_frame; --i) {
281     ;
282   }
283
284   track->last_marker = i;
285   return i;
286 }
287
288 static void retrieve_next_higher_usable_frame(
289     StabContext *ctx, MovieTrackingTrack *track, int i, int ref_frame, int *next_higher)
290 {
291   MovieTrackingMarker *markers = track->markers;
292   int end = track->markersnr;
293   BLI_assert(0 <= i && i < end);
294
295   while (i < end &&
296          (markers[i].framenr < ref_frame || is_effectively_disabled(ctx, track, &markers[i]))) {
297     ++i;
298   }
299   if (i < end && markers[i].framenr < *next_higher) {
300     BLI_assert(markers[i].framenr >= ref_frame);
301     *next_higher = markers[i].framenr;
302   }
303 }
304
305 static void retrieve_next_lower_usable_frame(
306     StabContext *ctx, MovieTrackingTrack *track, int i, int ref_frame, int *next_lower)
307 {
308   MovieTrackingMarker *markers = track->markers;
309   BLI_assert(0 <= i && i < track->markersnr);
310   while (i >= 0 &&
311          (markers[i].framenr > ref_frame || is_effectively_disabled(ctx, track, &markers[i]))) {
312     --i;
313   }
314   if (0 <= i && markers[i].framenr > *next_lower) {
315     BLI_assert(markers[i].framenr <= ref_frame);
316     *next_lower = markers[i].framenr;
317   }
318 }
319
320 /* Find closest frames with usable stabilization data.
321  * A frame counts as _usable_ when there is at least one track marked for
322  * translation stabilization, which has an enabled tracking marker at this very
323  * frame. We search both for the next lower and next higher position, to allow
324  * the caller to interpolate gaps and to extrapolate at the ends of the
325  * definition range.
326  *
327  * NOTE: Regarding performance note that the individual tracks will cache the
328  *       last search position.
329  */
330 static void find_next_working_frames(StabContext *ctx,
331                                      int framenr,
332                                      int *next_lower,
333                                      int *next_higher)
334 {
335   for (MovieTrackingTrack *track = ctx->tracking->tracks.first; track != NULL;
336        track = track->next) {
337     if (is_usable_for_stabilization(ctx, track)) {
338       int startpoint = search_closest_marker_index(track, framenr);
339       retrieve_next_higher_usable_frame(ctx, track, startpoint, framenr, next_higher);
340       retrieve_next_lower_usable_frame(ctx, track, startpoint, framenr, next_lower);
341     }
342   }
343 }
344
345 /* Find active (enabled) marker closest to the reference frame. */
346 static MovieTrackingMarker *get_closest_marker(StabContext *ctx,
347                                                MovieTrackingTrack *track,
348                                                int ref_frame)
349 {
350   int next_lower = MINAFRAME;
351   int next_higher = MAXFRAME;
352   int i = search_closest_marker_index(track, ref_frame);
353   retrieve_next_higher_usable_frame(ctx, track, i, ref_frame, &next_higher);
354   retrieve_next_lower_usable_frame(ctx, track, i, ref_frame, &next_lower);
355
356   if ((next_higher - ref_frame) < (ref_frame - next_lower)) {
357     return BKE_tracking_marker_get_exact(track, next_higher);
358   }
359   else {
360     return BKE_tracking_marker_get_exact(track, next_lower);
361   }
362 }
363
364 /* Retrieve tracking data, if available and applicable for this frame.
365  * The returned weight value signals the validity; data recorded for this
366  * tracking marker on the exact requested frame is output with the full weight
367  * of this track, while gaps in the data sequence cause the weight to go to zero.
368  */
369 static MovieTrackingMarker *get_tracking_data_point(StabContext *ctx,
370                                                     MovieTrackingTrack *track,
371                                                     int framenr,
372                                                     float *r_weight)
373 {
374   MovieTrackingMarker *marker = BKE_tracking_marker_get_exact(track, framenr);
375   if (marker != NULL && !(marker->flag & MARKER_DISABLED)) {
376     *r_weight = get_animated_weight(ctx, track, framenr);
377     return marker;
378   }
379   else {
380     /* No marker at this frame (=gap) or marker disabled. */
381     *r_weight = 0.0f;
382     return NULL;
383   }
384 }
385
386 /* Define the reference point for rotation/scale measurement and compensation.
387  * The stabilizator works by assuming the image was distorted by a affine linear
388  * transform, i.e. it was rotated and stretched around this reference point
389  * (pivot point) and then shifted laterally. Any scale and orientation changes
390  * will be picked up relative to this point. And later the image will be
391  * stabilized by rotating around this point. The result can only be as
392  * accurate as this pivot point actually matches the real rotation center
393  * of the actual movements. Thus any scheme to define a pivot point is
394  * always guesswork.
395  *
396  * As a simple default, we use the weighted average of the location markers
397  * of the current frame as pivot point. TODO It is planned to add further
398  * options,  like e.g. anchoring the pivot point at the canvas. Moreover,
399  * it is planned to allow for a user controllable offset.
400  */
401 static void setup_pivot(const float ref_pos[2], float r_pivot[2])
402 {
403   zero_v2(r_pivot); /* TODO: add an animated offset position here. */
404   add_v2_v2(r_pivot, ref_pos);
405 }
406
407 /* Calculate the contribution of a single track at the time position (frame) of
408  * the given marker. Each track has a local reference frame, which is as close
409  * as possible to the global anchor_frame. Thus the translation contribution is
410  * comprised of the offset relative to the image position at that reference
411  * frame, plus a guess of the contribution for the time span between the
412  * anchor_frame and the local reference frame of this track. The constant part
413  * of this contribution is precomputed initially. At the anchor_frame, by
414  * definition the contribution of all tracks is zero, keeping the frame in place.
415  *
416  * track_ref is per track baseline contribution at reference frame; filled in at
417  *           initialization
418  * marker is tracking data to use as contribution for current frame.
419  * result_offset is a total cumulated contribution of this track,
420  *               relative to the stabilization anchor_frame,
421  *               in normalized (0...1) coordinates.
422  */
423 static void translation_contribution(TrackStabilizationBase *track_ref,
424                                      MovieTrackingMarker *marker,
425                                      float result_offset[2])
426 {
427   add_v2_v2v2(result_offset, track_ref->stabilization_offset_base, marker->pos);
428 }
429
430 /* Similar to the ::translation_contribution(), the rotation contribution is
431  * comprised of the contribution by this individual track, and the averaged
432  * contribution from anchor_frame to the ref point of this track.
433  * - Contribution is in terms of angles, -pi < angle < +pi, and all averaging
434  *   happens in this domain.
435  * - Yet the actual measurement happens as vector between pivot and the current
436  *   tracking point
437  * - Currently we use the center of frame as approximation for the rotation pivot
438  *   point.
439  * - Moreover, the pivot point has to be compensated for the already determined
440  *   shift offset, in order to get the pure rotation around the pivot.
441  *   To turn this into a _contribution_, the likewise corrected angle at the
442  *   reference frame has to be subtracted, to get only the pure angle difference
443  *   this tracking point has captured.
444  * - To get from vectors to angles, we have to go through an arcus tangens,
445  *   which involves the issue of the definition range: the resulting angles will
446  *   flip by 360deg when the measured vector passes from the 2nd to the third
447  *   quadrant, thus messing up the average calculation. Since _any_ tracking
448  *   point might be used, these problems are quite common in practice.
449  * - Thus we perform the subtraction of the reference and the addition of the
450  *   baseline contribution in polar coordinates as simple addition of angles;
451  *   since these parts are fixed, we can bake them into a rotation matrix.
452  *   With this approach, the border of the arcus tangens definition range will
453  *   be reached only, when the _whole_ contribution approaches +- 180deg,
454  *   meaning we've already tilted the frame upside down. This situation is way
455  *   less common and can be tolerated.
456  * - As an additional feature, when activated, also changes in image scale
457  *   relative to the rotation center can be picked up. To handle those values
458  *   in the same framework, we average the scales as logarithms.
459  *
460  * aspect is a total aspect ratio of the undistorted image (includes fame and
461  * pixel aspect). The function returns a quality factor, which can be used
462  * to damp the contributions of points in close proximity to the pivot point,
463  * since such contributions might be dominated by rounding errors and thus
464  * poison the calculated average. When the quality factor goes towards zero,
465  * the weight of this contribution should be reduced accordingly.
466  */
467 static float rotation_contribution(TrackStabilizationBase *track_ref,
468                                    MovieTrackingMarker *marker,
469                                    const float aspect,
470                                    const float pivot[2],
471                                    float *result_angle,
472                                    float *result_scale)
473 {
474   float len, quality;
475   float pos[2];
476   sub_v2_v2v2(pos, marker->pos, pivot);
477
478   pos[0] *= aspect;
479   mul_m2v2(track_ref->stabilization_rotation_base, pos);
480
481   *result_angle = atan2f(pos[1], pos[0]);
482
483   len = len_v2(pos);
484
485   /* prevent points very close to the pivot point from poisoning the result */
486   quality = 1 - expf(-len * len / SCALE_ERROR_LIMIT_BIAS * SCALE_ERROR_LIMIT_BIAS);
487   len += SCALE_ERROR_LIMIT_BIAS;
488
489   *result_scale = len * track_ref->stabilization_scale_base;
490   BLI_assert(0.0 < *result_scale);
491
492   return quality;
493 }
494
495 /* Workaround to allow for rotation around an arbitrary pivot point.
496  * Currently, the public API functions do not support this flexibility.
497  * Rather, rotation will always be applied around a fixed origin.
498  * As a workaround, we shift the image after rotation to match the
499  * desired rotation centre. And since this offset needs to be applied
500  * after the rotation and scaling, we can collapse it with the
501  * translation compensation, which is also a lateral shift (offset).
502  * The offset to apply is intended_pivot - rotated_pivot
503  */
504 static void compensate_rotation_center(const int size,
505                                        float aspect,
506                                        const float angle,
507                                        const float scale,
508                                        const float pivot[2],
509                                        float result_translation[2])
510 {
511   const float origin[2] = {0.5f * aspect * size, 0.5f * size};
512   float intended_pivot[2], rotated_pivot[2];
513   float rotation_mat[2][2];
514
515   copy_v2_v2(intended_pivot, pivot);
516   copy_v2_v2(rotated_pivot, pivot);
517   angle_to_mat2(rotation_mat, +angle);
518   sub_v2_v2(rotated_pivot, origin);
519   mul_m2v2(rotation_mat, rotated_pivot);
520   mul_v2_fl(rotated_pivot, scale);
521   add_v2_v2(rotated_pivot, origin);
522   add_v2_v2(result_translation, intended_pivot);
523   sub_v2_v2(result_translation, rotated_pivot);
524 }
525
526 /* Weighted average of the per track cumulated contributions at given frame.
527  * Returns truth if all desired calculations could be done and all averages are
528  * available.
529  *
530  * NOTE: Even if the result is not `true`, the returned translation and angle
531  *       are always sensible and as good as can be. Especially in the
532  *       initialization phase we might not be able to get any average (yet) or
533  *       get only a translation value. Since initialization visits tracks in a
534  *       specific order, starting from anchor_frame, the result is logically
535  *       correct non the less. But under normal operation conditions,
536  *       a result of `false` should disable the stabilization function
537  */
538 static bool average_track_contributions(StabContext *ctx,
539                                         int framenr,
540                                         float aspect,
541                                         float r_translation[2],
542                                         float r_pivot[2],
543                                         float *r_angle,
544                                         float *r_scale_step)
545 {
546   bool ok;
547   float weight_sum;
548   MovieTrackingTrack *track;
549   MovieTracking *tracking = ctx->tracking;
550   MovieTrackingStabilization *stab = &tracking->stabilization;
551   float ref_pos[2];
552   BLI_assert(stab->flag & TRACKING_2D_STABILIZATION);
553
554   zero_v2(r_translation);
555   *r_scale_step = 0.0f; /* logarithm */
556   *r_angle = 0.0f;
557
558   zero_v2(ref_pos);
559
560   ok = false;
561   weight_sum = 0.0f;
562   for (track = tracking->tracks.first; track; track = track->next) {
563     if (!is_init_for_stabilization(ctx, track)) {
564       continue;
565     }
566     if (track->flag & TRACK_USE_2D_STAB) {
567       float weight = 0.0f;
568       MovieTrackingMarker *marker = get_tracking_data_point(ctx, track, framenr, &weight);
569       if (marker) {
570         TrackStabilizationBase *stabilization_base = access_stabilization_baseline_data(ctx,
571                                                                                         track);
572         BLI_assert(stabilization_base != NULL);
573         float offset[2];
574         weight_sum += weight;
575         translation_contribution(stabilization_base, marker, offset);
576         r_translation[0] += weight * offset[0];
577         r_translation[1] += weight * offset[1];
578         ref_pos[0] += weight * marker->pos[0];
579         ref_pos[1] += weight * marker->pos[1];
580         ok |= (weight_sum > EPSILON_WEIGHT);
581       }
582     }
583   }
584   if (!ok) {
585     return false;
586   }
587
588   ref_pos[0] /= weight_sum;
589   ref_pos[1] /= weight_sum;
590   r_translation[0] /= weight_sum;
591   r_translation[1] /= weight_sum;
592   setup_pivot(ref_pos, r_pivot);
593
594   if (!(stab->flag & TRACKING_STABILIZE_ROTATION)) {
595     return ok;
596   }
597
598   ok = false;
599   weight_sum = 0.0f;
600   for (track = tracking->tracks.first; track; track = track->next) {
601     if (!is_init_for_stabilization(ctx, track)) {
602       continue;
603     }
604     if (track->flag & TRACK_USE_2D_STAB_ROT) {
605       float weight = 0.0f;
606       MovieTrackingMarker *marker = get_tracking_data_point(ctx, track, framenr, &weight);
607       if (marker) {
608         TrackStabilizationBase *stabilization_base = access_stabilization_baseline_data(ctx,
609                                                                                         track);
610         BLI_assert(stabilization_base != NULL);
611         float rotation, scale, quality;
612         quality = rotation_contribution(
613             stabilization_base, marker, aspect, r_pivot, &rotation, &scale);
614         weight *= quality;
615         weight_sum += weight;
616         *r_angle += rotation * weight;
617         if (stab->flag & TRACKING_STABILIZE_SCALE) {
618           *r_scale_step += logf(scale) * weight;
619         }
620         else {
621           *r_scale_step = 0;
622         }
623         ok |= (weight_sum > EPSILON_WEIGHT);
624       }
625     }
626   }
627   if (ok) {
628     *r_scale_step /= weight_sum;
629     *r_angle /= weight_sum;
630   }
631   else {
632     /* We reach this point because translation could be calculated,
633      * but rotation/scale found no data to work on.
634      */
635     *r_scale_step = 0.0f;
636     *r_angle = 0.0f;
637   }
638   return true;
639 }
640
641 /* Calculate weight center of location tracks for given frame.
642  * This function performs similar calculations as average_track_contributions(),
643  * but does not require the tracks to be initialized for stabilisation. Moreover,
644  * when there is no usable tracking data for the given frame number, data from
645  * a neighbouring frame is used. Thus this function can be used to calculate
646  * a starting point on initialization.
647  */
648 static void average_marker_positions(StabContext *ctx, int framenr, float r_ref_pos[2])
649 {
650   bool ok = false;
651   float weight_sum;
652   MovieTrackingTrack *track;
653   MovieTracking *tracking = ctx->tracking;
654
655   zero_v2(r_ref_pos);
656   weight_sum = 0.0f;
657   for (track = tracking->tracks.first; track; track = track->next) {
658     if (track->flag & TRACK_USE_2D_STAB) {
659       float weight = 0.0f;
660       MovieTrackingMarker *marker = get_tracking_data_point(ctx, track, framenr, &weight);
661       if (marker) {
662         weight_sum += weight;
663         r_ref_pos[0] += weight * marker->pos[0];
664         r_ref_pos[1] += weight * marker->pos[1];
665         ok |= (weight_sum > EPSILON_WEIGHT);
666       }
667     }
668   }
669   if (ok) {
670     r_ref_pos[0] /= weight_sum;
671     r_ref_pos[1] /= weight_sum;
672   }
673   else {
674     /* No usable tracking data on any track on this frame.
675      * Use data from neighbouring frames to extrapolate...
676      */
677     int next_lower = MINAFRAME;
678     int next_higher = MAXFRAME;
679     use_values_from_fcurves(ctx, true);
680     for (track = tracking->tracks.first; track; track = track->next) {
681       /* Note: we deliberately do not care if this track
682        *       is already initialized for stabilisation */
683       if (track->flag & TRACK_USE_2D_STAB) {
684         int startpoint = search_closest_marker_index(track, framenr);
685         retrieve_next_higher_usable_frame(ctx, track, startpoint, framenr, &next_higher);
686         retrieve_next_lower_usable_frame(ctx, track, startpoint, framenr, &next_lower);
687       }
688     }
689     if (next_lower >= MINFRAME) {
690       /* use next usable frame to the left.
691        * Also default to this frame when we're in a gap */
692       average_marker_positions(ctx, next_lower, r_ref_pos);
693     }
694     else if (next_higher < MAXFRAME) {
695       average_marker_positions(ctx, next_higher, r_ref_pos);
696     }
697     use_values_from_fcurves(ctx, false);
698   }
699 }
700
701 /* Linear interpolation of data retrieved at two measurement points.
702  * This function is used to fill gaps in the middle of the covered area,
703  * at frames without any usable tracks for stabilization.
704  *
705  * framenr is a position to interpolate for.
706  * frame_a is a valid measurement point below framenr
707  * frame_b is a valid measurement point above framenr
708  * Returns truth if both measurements could actually be retrieved.
709  * Otherwise output parameters remain unaltered
710  */
711 static bool interpolate_averaged_track_contributions(StabContext *ctx,
712                                                      int framenr,
713                                                      int frame_a,
714                                                      int frame_b,
715                                                      const float aspect,
716                                                      float r_translation[2],
717                                                      float r_pivot[2],
718                                                      float *r_angle,
719                                                      float *r_scale_step)
720 {
721   float t, s;
722   float trans_a[2], trans_b[2];
723   float angle_a, angle_b;
724   float scale_a, scale_b;
725   float pivot_a[2], pivot_b[2];
726   bool success = false;
727
728   BLI_assert(frame_a <= frame_b);
729   BLI_assert(frame_a <= framenr);
730   BLI_assert(framenr <= frame_b);
731
732   t = ((float)framenr - frame_a) / (frame_b - frame_a);
733   s = 1.0f - t;
734
735   success = average_track_contributions(
736       ctx, frame_a, aspect, trans_a, pivot_a, &angle_a, &scale_a);
737   if (!success) {
738     return false;
739   }
740   success = average_track_contributions(
741       ctx, frame_b, aspect, trans_b, pivot_b, &angle_b, &scale_b);
742   if (!success) {
743     return false;
744   }
745
746   interp_v2_v2v2(r_translation, trans_a, trans_b, t);
747   interp_v2_v2v2(r_pivot, pivot_a, pivot_b, t);
748   *r_scale_step = s * scale_a + t * scale_b;
749   *r_angle = s * angle_a + t * angle_b;
750   return true;
751 }
752
753 /* Reorder tracks starting with those providing a tracking data frame
754  * closest to the global anchor_frame. Tracks with a gap at anchor_frame or
755  * starting farer away from anchor_frame altogether will be visited later.
756  * This allows to build up baseline contributions incrementally.
757  *
758  * order is an array for sorting the tracks. Must be of suitable size to hold
759  * all tracks.
760  * Returns number of actually usable tracks, can be less than the overall number
761  * of tracks.
762  *
763  * NOTE: After returning, the order array holds entries up to the number of
764  *       usable tracks, appropriately sorted starting with the closest tracks.
765  *       Initialization includes disabled tracks, since they might be enabled
766  *       through automation later.
767  */
768 static int establish_track_initialization_order(StabContext *ctx, TrackInitOrder *order)
769 {
770   size_t tracknr = 0;
771   MovieTrackingTrack *track;
772   MovieTracking *tracking = ctx->tracking;
773   int anchor_frame = tracking->stabilization.anchor_frame;
774
775   for (track = tracking->tracks.first; track != NULL; track = track->next) {
776     MovieTrackingMarker *marker;
777     order[tracknr].data = track;
778     marker = get_closest_marker(ctx, track, anchor_frame);
779     if (marker != NULL && (track->flag & (TRACK_USE_2D_STAB | TRACK_USE_2D_STAB_ROT))) {
780       order[tracknr].sort_value = abs(marker->framenr - anchor_frame);
781       order[tracknr].reference_frame = marker->framenr;
782       ++tracknr;
783     }
784   }
785   if (tracknr) {
786     qsort(order, tracknr, sizeof(TrackInitOrder), BLI_sortutil_cmp_int);
787   }
788   return tracknr;
789 }
790
791 /* Setup the constant part of this track's contribution to the determined frame
792  * movement. Tracks usually don't provide tracking data for every frame. Thus,
793  * for determining data at a given frame, we split up the contribution into a
794  * part covered by actual measurements on this track, and the initial gap
795  * between this track's reference frame and the global anchor_frame.
796  * The (missing) data for the gap can be substituted by the average offset
797  * observed by the other tracks covering the gap. This approximation doesn't
798  * introduce wrong data, but it records data with incorrect weight. A totally
799  * correct solution would require us to average the contribution per frame, and
800  * then integrate stepwise over all frames -- which of course would be way more
801  * expensive, especially for longer clips. To the contrary, our solution
802  * cumulates the total contribution per track and averages afterwards over all
803  * tracks; it can thus be calculated just based on the data of a single frame,
804  * plus the "baseline" for the reference frame, which is what we are computing
805  * here.
806  *
807  * Since we're averaging _contributions_, we have to calculate the _difference_
808  * of the measured position at current frame and the position at the reference
809  * frame. But the "reference" part of this difference is constant and can thus
810  * be packed together with the baseline contribution into a single precomputed
811  * vector per track.
812  *
813  * In case of the rotation contribution, the principle is the same, but we have
814  * to compensate for the already determined translation and measure the pure
815  * rotation, simply because this is how we model the offset: shift plus rotation
816  * around the shifted rotation center. To circumvent problems with the
817  * definition range of the arcus tangens function, we perform this baseline
818  * addition and reference angle subtraction in polar coordinates and bake this
819  * operation into a precomputed rotation matrix.
820  *
821  * track is a track to be initialized to initialize
822  * reference_frame is a local frame for this track, the closest pick to the
823  *                 global anchor_frame.
824  * aspect is a total aspect ratio of the undistorted image (includes fame and
825  *        pixel aspect).
826  * target_pos is a possibly animated target position as set by the user for
827  *            the reference_frame
828  * average_translation is a value observed by the _other_ tracks for the gap
829  *                     between reference_frame and anchor_frame. This
830  *                     average must not contain contributions of frames
831  *                     not yet initialized
832  * average_angle in a similar way, the rotation value observed by the
833  *               _other_ tracks.
834  * average_scale_step is an image scale factor observed on average by the other
835  *                    tracks for this frame. This value is recorded and
836  *                    averaged as logarithm. The recorded scale changes
837  *                    are damped for very small contributions, to limit
838  *                    the effect of probe points approaching the pivot
839  *                    too closely.
840  *
841  * NOTE: when done, this track is marked as initialized
842  */
843 static void initialize_track_for_stabilization(StabContext *ctx,
844                                                MovieTrackingTrack *track,
845                                                int reference_frame,
846                                                float aspect,
847                                                const float average_translation[2],
848                                                const float pivot[2],
849                                                const float average_angle,
850                                                const float average_scale_step)
851 {
852   float pos[2], angle, len;
853   TrackStabilizationBase *local_data = access_stabilization_baseline_data(ctx, track);
854   MovieTrackingMarker *marker = BKE_tracking_marker_get_exact(track, reference_frame);
855   /* Logic for initialization order ensures there *is* a marker on that
856    * very frame.
857    */
858   BLI_assert(marker != NULL);
859   BLI_assert(local_data != NULL);
860
861   /* Per track baseline value for translation. */
862   sub_v2_v2v2(local_data->stabilization_offset_base, average_translation, marker->pos);
863
864   /* Per track baseline value for rotation. */
865   sub_v2_v2v2(pos, marker->pos, pivot);
866
867   pos[0] *= aspect;
868   angle = average_angle - atan2f(pos[1], pos[0]);
869   angle_to_mat2(local_data->stabilization_rotation_base, angle);
870
871   /* Per track baseline value for zoom. */
872   len = len_v2(pos) + SCALE_ERROR_LIMIT_BIAS;
873   local_data->stabilization_scale_base = expf(average_scale_step) / len;
874
875   local_data->is_init_for_stabilization = true;
876 }
877
878 static void initialize_all_tracks(StabContext *ctx, float aspect)
879 {
880   size_t i, track_len = 0;
881   MovieClip *clip = ctx->clip;
882   MovieTracking *tracking = ctx->tracking;
883   MovieTrackingTrack *track;
884   TrackInitOrder *order;
885
886   /* Attempt to start initialization at anchor_frame.
887    * By definition, offset contribution is zero there.
888    */
889   int reference_frame = tracking->stabilization.anchor_frame;
890   float average_angle = 0, average_scale_step = 0;
891   float average_translation[2], average_pos[2], pivot[2];
892   zero_v2(average_translation);
893   zero_v2(pivot);
894
895   /* Initialize private working data. */
896   for (track = tracking->tracks.first; track != NULL; track = track->next) {
897     TrackStabilizationBase *local_data = access_stabilization_baseline_data(ctx, track);
898     if (!local_data) {
899       local_data = MEM_callocN(sizeof(TrackStabilizationBase),
900                                "2D stabilization per track baseline data");
901       attach_stabilization_baseline_data(ctx, track, local_data);
902     }
903     BLI_assert(local_data != NULL);
904     local_data->track_weight_curve = retrieve_track_weight_animation(clip, track);
905     local_data->is_init_for_stabilization = false;
906
907     ++track_len;
908   }
909   if (!track_len) {
910     return;
911   }
912
913   order = MEM_mallocN(track_len * sizeof(TrackInitOrder), "stabilization track order");
914   if (!order) {
915     return;
916   }
917
918   track_len = establish_track_initialization_order(ctx, order);
919   if (track_len == 0) {
920     goto cleanup;
921   }
922
923   /* starting point for pivot, before having initialized any track */
924   average_marker_positions(ctx, reference_frame, average_pos);
925   setup_pivot(average_pos, pivot);
926
927   for (i = 0; i < track_len; ++i) {
928     track = order[i].data;
929     if (reference_frame != order[i].reference_frame) {
930       reference_frame = order[i].reference_frame;
931       average_track_contributions(ctx,
932                                   reference_frame,
933                                   aspect,
934                                   average_translation,
935                                   pivot,
936                                   &average_angle,
937                                   &average_scale_step);
938     }
939     initialize_track_for_stabilization(ctx,
940                                        track,
941                                        reference_frame,
942                                        aspect,
943                                        average_translation,
944                                        pivot,
945                                        average_angle,
946                                        average_scale_step);
947   }
948
949 cleanup:
950   MEM_freeN(order);
951 }
952
953 /* Retrieve the measurement of frame movement by averaging contributions of
954  * active tracks.
955  *
956  * translation is a measurement in normalized 0..1 coordinates.
957  * angle is a measurement in radians -pi..+pi counter clockwise relative to
958  *       translation compensated frame center
959  * scale_step is a measurement of image scale changes, in logarithmic scale
960  *            (zero means scale == 1)
961  * Returns calculation enabled and all data retrieved as expected for this frame.
962  *
963  * NOTE: when returning `false`, output parameters are reset to neutral values.
964  */
965 static bool stabilization_determine_offset_for_frame(StabContext *ctx,
966                                                      int framenr,
967                                                      float aspect,
968                                                      float r_translation[2],
969                                                      float r_pivot[2],
970                                                      float *r_angle,
971                                                      float *r_scale_step)
972 {
973   bool success = false;
974
975   /* Early output if stabilization is disabled. */
976   if ((ctx->stab->flag & TRACKING_2D_STABILIZATION) == 0) {
977     zero_v2(r_translation);
978     *r_scale_step = 0.0f;
979     *r_angle = 0.0f;
980     return false;
981   }
982
983   success = average_track_contributions(
984       ctx, framenr, aspect, r_translation, r_pivot, r_angle, r_scale_step);
985   if (!success) {
986     /* Try to hold extrapolated settings beyond the definition range
987      * and to interpolate in gaps without any usable tracking data
988      * to prevent sudden jump to image zero position.
989      */
990     int next_lower = MINAFRAME;
991     int next_higher = MAXFRAME;
992     use_values_from_fcurves(ctx, true);
993     find_next_working_frames(ctx, framenr, &next_lower, &next_higher);
994     if (next_lower >= MINFRAME && next_higher < MAXFRAME) {
995       success = interpolate_averaged_track_contributions(ctx,
996                                                          framenr,
997                                                          next_lower,
998                                                          next_higher,
999                                                          aspect,
1000                                                          r_translation,
1001                                                          r_pivot,
1002                                                          r_angle,
1003                                                          r_scale_step);
1004     }
1005     else if (next_higher < MAXFRAME) {
1006       /* Before start of stabilized range: extrapolate start point
1007        * settings.
1008        */
1009       success = average_track_contributions(
1010           ctx, next_higher, aspect, r_translation, r_pivot, r_angle, r_scale_step);
1011     }
1012     else if (next_lower >= MINFRAME) {
1013       /* After end of stabilized range: extrapolate end point settings. */
1014       success = average_track_contributions(
1015           ctx, next_lower, aspect, r_translation, r_pivot, r_angle, r_scale_step);
1016     }
1017     use_values_from_fcurves(ctx, false);
1018   }
1019   return success;
1020 }
1021
1022 /* Calculate stabilization data (translation, scale and rotation) from given raw
1023  * measurements. Result is in absolute image dimensions (expanded image, square
1024  * pixels), includes automatic or manual scaling and compensates for a target
1025  * frame position, if given.
1026  *
1027  * size is a size of the expanded image, the width in pixels is size * aspect.
1028  * aspect is a ratio (width / height) of the effective canvas (square pixels).
1029  * do_compensate denotes whether to actually output values necessary to
1030  *               _compensate_ the determined frame movement.
1031  *               Otherwise, the effective target movement is returned.
1032  */
1033 static void stabilization_calculate_data(StabContext *ctx,
1034                                          int framenr,
1035                                          int size,
1036                                          float aspect,
1037                                          bool do_compensate,
1038                                          float scale_step,
1039                                          float r_translation[2],
1040                                          float r_pivot[2],
1041                                          float *r_scale,
1042                                          float *r_angle)
1043 {
1044   float target_pos[2], target_scale;
1045   float scaleinf = get_animated_scaleinf(ctx, framenr);
1046
1047   if (ctx->stab->flag & TRACKING_STABILIZE_SCALE) {
1048     *r_scale = expf(scale_step * scaleinf); /* Averaged in log scale */
1049   }
1050   else {
1051     *r_scale = 1.0f;
1052   }
1053
1054   mul_v2_fl(r_translation, get_animated_locinf(ctx, framenr));
1055   *r_angle *= get_animated_rotinf(ctx, framenr);
1056
1057   /* Compensate for a target frame position.
1058    * This allows to follow tracking / panning shots in a semi manual fashion,
1059    * when animating the settings for the target frame position.
1060    */
1061   get_animated_target_pos(ctx, framenr, target_pos);
1062   sub_v2_v2(r_translation, target_pos);
1063   *r_angle -= get_animated_target_rot(ctx, framenr);
1064   target_scale = get_animated_target_scale(ctx, framenr);
1065   if (target_scale != 0.0f) {
1066     *r_scale /= target_scale;
1067     /* target_scale is an expected/intended reference zoom value */
1068   }
1069
1070   /* Convert from relative to absolute coordinates, square pixels. */
1071   r_translation[0] *= (float)size * aspect;
1072   r_translation[1] *= (float)size;
1073   r_pivot[0] *= (float)size * aspect;
1074   r_pivot[1] *= (float)size;
1075
1076   /* Output measured data, or inverse of the measured values for
1077    * compensation?
1078    */
1079   if (do_compensate) {
1080     mul_v2_fl(r_translation, -1.0f);
1081     *r_angle *= -1.0f;
1082     if (*r_scale != 0.0f) {
1083       *r_scale = 1.0f / *r_scale;
1084     }
1085   }
1086 }
1087
1088 static void stabilization_data_to_mat4(float pixel_aspect,
1089                                        const float pivot[2],
1090                                        const float translation[2],
1091                                        float scale,
1092                                        float angle,
1093                                        float r_mat[4][4])
1094 {
1095   float translation_mat[4][4], rotation_mat[4][4], scale_mat[4][4], pivot_mat[4][4],
1096       inv_pivot_mat[4][4], aspect_mat[4][4], inv_aspect_mat[4][4];
1097   float scale_vector[3] = {scale, scale, 1.0f};
1098
1099   unit_m4(translation_mat);
1100   unit_m4(rotation_mat);
1101   unit_m4(scale_mat);
1102   unit_m4(aspect_mat);
1103   unit_m4(pivot_mat);
1104   unit_m4(inv_pivot_mat);
1105
1106   /* aspect ratio correction matrix */
1107   aspect_mat[0][0] /= pixel_aspect;
1108   invert_m4_m4(inv_aspect_mat, aspect_mat);
1109
1110   add_v2_v2(pivot_mat[3], pivot);
1111   sub_v2_v2(inv_pivot_mat[3], pivot);
1112
1113   size_to_mat4(scale_mat, scale_vector);      /* scale matrix */
1114   add_v2_v2(translation_mat[3], translation); /* translation matrix */
1115   rotate_m4(rotation_mat, 'Z', angle);        /* rotation matrix */
1116
1117   /* Compose transformation matrix. */
1118   mul_m4_series(r_mat,
1119                 aspect_mat,
1120                 translation_mat,
1121                 pivot_mat,
1122                 scale_mat,
1123                 rotation_mat,
1124                 inv_pivot_mat,
1125                 inv_aspect_mat);
1126 }
1127
1128 /* Calculate scale factor necessary to eliminate black image areas
1129  * caused by the compensating movements of the stabilizator.
1130  * This function visits every frame where stabilisation data is
1131  * available and determines the factor for this frame. The overall
1132  * largest factor found is returned as result.
1133  *
1134  * NOTE: all tracks need to be initialized before calling this function.
1135  */
1136 static float calculate_autoscale_factor(StabContext *ctx, int size, float aspect)
1137 {
1138   MovieTrackingStabilization *stab = ctx->stab;
1139   float pixel_aspect = ctx->tracking->camera.pixel_aspect;
1140   int height = size, width = aspect * size;
1141
1142   int sfra = INT_MAX, efra = INT_MIN, cfra;
1143   float scale = 1.0f, scale_step = 0.0f;
1144   MovieTrackingTrack *track;
1145
1146   /* Calculate maximal frame range of tracks where stabilization is active. */
1147   for (track = ctx->tracking->tracks.first; track; track = track->next) {
1148     if ((track->flag & TRACK_USE_2D_STAB) ||
1149         ((stab->flag & TRACKING_STABILIZE_ROTATION) && (track->flag & TRACK_USE_2D_STAB_ROT))) {
1150       int first_frame = track->markers[0].framenr;
1151       int last_frame = track->markers[track->markersnr - 1].framenr;
1152       sfra = min_ii(sfra, first_frame);
1153       efra = max_ii(efra, last_frame);
1154     }
1155   }
1156
1157   use_values_from_fcurves(ctx, true);
1158   for (cfra = sfra; cfra <= efra; cfra++) {
1159     float translation[2], pivot[2], angle, tmp_scale;
1160     float mat[4][4];
1161     const float points[4][2] = {{0.0f, 0.0f}, {0.0f, height}, {width, height}, {width, 0.0f}};
1162     const bool do_compensate = true;
1163     /* Calculate stabilization parameters for the current frame. */
1164     stabilization_determine_offset_for_frame(
1165         ctx, cfra, aspect, translation, pivot, &angle, &scale_step);
1166     stabilization_calculate_data(ctx,
1167                                  cfra,
1168                                  size,
1169                                  aspect,
1170                                  do_compensate,
1171                                  scale_step,
1172                                  translation,
1173                                  pivot,
1174                                  &tmp_scale,
1175                                  &angle);
1176     /* Compose transformation matrix. */
1177     /* NOTE: Here we operate in NON-COMPENSATED coordinates, meaning we have
1178      * to construct transformation matrix using proper pivot point.
1179      * Compensation for that will happen later on.
1180      */
1181     stabilization_data_to_mat4(pixel_aspect, pivot, translation, tmp_scale, angle, mat);
1182     /* Investigate the transformed border lines for this frame;
1183      * find out, where it cuts the original frame.
1184      */
1185     for (int edge_index = 0; edge_index < 4; edge_index++) {
1186       /* Calculate coordinates of stabilized frame edge points.
1187        * Use matrix multiplication here so we operate in homogeneous
1188        * coordinates.
1189        */
1190       float stable_edge_p1[3], stable_edge_p2[3];
1191       copy_v2_v2(stable_edge_p1, points[edge_index]);
1192       copy_v2_v2(stable_edge_p2, points[(edge_index + 1) % 4]);
1193       stable_edge_p1[2] = stable_edge_p2[2] = 0.0f;
1194       mul_m4_v3(mat, stable_edge_p1);
1195       mul_m4_v3(mat, stable_edge_p2);
1196       /* Now we iterate over all original frame corners (we call them
1197        * 'point' here) to see if there's black area between stabilized
1198        * frame edge and original point.
1199        */
1200       for (int point_index = 0; point_index < 4; point_index++) {
1201         const float point[3] = {points[point_index][0], points[point_index][1], 0.0f};
1202         /* Calculate vector which goes from first edge point to
1203          * second one.
1204          */
1205         float stable_edge_vec[3];
1206         sub_v3_v3v3(stable_edge_vec, stable_edge_p2, stable_edge_p1);
1207         /* Calculate vector which connects current frame point to
1208          * first edge point.
1209          */
1210         float point_to_edge_start_vec[3];
1211         sub_v3_v3v3(point_to_edge_start_vec, point, stable_edge_p1);
1212         /* Use this two vectors to check whether frame point is inside
1213          * of the stabilized frame or not.
1214          * If the point is inside, there is no black area happening
1215          * and no scaling required for it.
1216          */
1217         if (cross_v2v2(stable_edge_vec, point_to_edge_start_vec) >= 0.0f) {
1218           /* We are scaling around motion-compensated pivot point. */
1219           float scale_pivot[2];
1220           add_v2_v2v2(scale_pivot, pivot, translation);
1221           /* Calculate line which goes via `point` and parallel to
1222            * the stabilized frame edge. This line is coming via
1223            * `point` and `point2` at the end.
1224            */
1225           float point2[2];
1226           add_v2_v2v2(point2, point, stable_edge_vec);
1227           /* Calculate actual distance between pivot point and
1228            * the stabilized frame edge. Then calculate distance
1229            * between pivot point and line which goes via actual
1230            * corner and is parallel to the edge.
1231            *
1232            * Dividing one by another will give us required scale
1233            * factor to get rid of black areas.
1234            */
1235           float real_dist = dist_to_line_v2(scale_pivot, stable_edge_p1, stable_edge_p2);
1236           float required_dist = dist_to_line_v2(scale_pivot, point, point2);
1237           const float S = required_dist / real_dist;
1238           scale = max_ff(scale, S);
1239         }
1240       }
1241     }
1242   }
1243   if (stab->maxscale > 0.0f) {
1244     scale = min_ff(scale, stab->maxscale);
1245   }
1246   use_values_from_fcurves(ctx, false);
1247
1248   return scale;
1249 }
1250
1251 /* Prepare working data and determine reference point for each track.
1252  *
1253  * NOTE: These calculations _could_ be cached and reused for all frames of the
1254  *       same clip. However, since proper initialization depends on (weight)
1255  *       animation and setup of tracks, ensuring consistency of cached init data
1256  *       turns out to be tricky, hard to maintain and generally not worth the
1257  *       effort. Thus we'll re-initialize on every frame.
1258  */
1259 static StabContext *init_stabilizer(MovieClip *clip, int size, float aspect)
1260 {
1261   StabContext *ctx = initialize_stabilization_working_context(clip);
1262   BLI_assert(ctx != NULL);
1263   initialize_all_tracks(ctx, aspect);
1264   if (ctx->stab->flag & TRACKING_AUTOSCALE) {
1265     ctx->stab->scale = 1.0;
1266     ctx->stab->scale = calculate_autoscale_factor(ctx, size, aspect);
1267   }
1268   /* By default, just use values for the global current frame. */
1269   use_values_from_fcurves(ctx, false);
1270   return ctx;
1271 }
1272
1273 /* === public interface functions === */
1274
1275 /* Get stabilization data (translation, scaling and angle) for a given frame.
1276  * Returned data describes how to compensate the detected movement, but with any
1277  * chosen scale factor already applied and any target frame position already
1278  * compensated. In case stabilization fails or is disabled, neutral values are
1279  * returned.
1280  *
1281  * framenr is a frame number, relative to the clip (not relative to the scene
1282  *         timeline)
1283  * width is an effective width of the canvas (square pixels), used to scale the
1284  *       determined translation
1285  *
1286  * Outputs:
1287  * - translation of the lateral shift, absolute canvas coordinates
1288  *   (square pixels).
1289  * - scale of the scaling to apply
1290  * - angle of the rotation angle, relative to the frame center
1291  */
1292 /* TODO(sergey): Use r_ prefix for output parameters here. */
1293 void BKE_tracking_stabilization_data_get(MovieClip *clip,
1294                                          int framenr,
1295                                          int width,
1296                                          int height,
1297                                          float translation[2],
1298                                          float *scale,
1299                                          float *angle)
1300 {
1301   StabContext *ctx = NULL;
1302   MovieTracking *tracking = &clip->tracking;
1303   bool enabled = (tracking->stabilization.flag & TRACKING_2D_STABILIZATION);
1304   /* Might become a parameter of a stabilization compositor node. */
1305   bool do_compensate = true;
1306   float scale_step = 0.0f;
1307   float pixel_aspect = tracking->camera.pixel_aspect;
1308   float aspect = (float)width * pixel_aspect / height;
1309   int size = height;
1310   float pivot[2];
1311
1312   if (enabled) {
1313     ctx = init_stabilizer(clip, size, aspect);
1314   }
1315
1316   if (enabled && stabilization_determine_offset_for_frame(
1317                      ctx, framenr, aspect, translation, pivot, angle, &scale_step)) {
1318     stabilization_calculate_data(
1319         ctx, framenr, size, aspect, do_compensate, scale_step, translation, pivot, scale, angle);
1320     compensate_rotation_center(size, aspect, *angle, *scale, pivot, translation);
1321   }
1322   else {
1323     zero_v2(translation);
1324     *scale = 1.0f;
1325     *angle = 0.0f;
1326   }
1327   discard_stabilization_working_context(ctx);
1328 }
1329
1330 typedef void (*interpolation_func)(struct ImBuf *, struct ImBuf *, float, float, int, int);
1331
1332 typedef struct TrackingStabilizeFrameInterpolationData {
1333   ImBuf *ibuf;
1334   ImBuf *tmpibuf;
1335   float (*mat)[4];
1336
1337   interpolation_func interpolation;
1338 } TrackingStabilizeFrameInterpolationData;
1339
1340 static void tracking_stabilize_frame_interpolation_cb(
1341     void *__restrict userdata, const int j, const ParallelRangeTLS *__restrict UNUSED(tls))
1342 {
1343   TrackingStabilizeFrameInterpolationData *data = userdata;
1344   ImBuf *ibuf = data->ibuf;
1345   ImBuf *tmpibuf = data->tmpibuf;
1346   float(*mat)[4] = data->mat;
1347
1348   interpolation_func interpolation = data->interpolation;
1349
1350   for (int i = 0; i < tmpibuf->x; i++) {
1351     float vec[3] = {i, j, 0.0f};
1352
1353     mul_v3_m4v3(vec, mat, vec);
1354
1355     interpolation(ibuf, tmpibuf, vec[0], vec[1], i, j);
1356   }
1357 }
1358
1359 /* Stabilize given image buffer using stabilization data for a specified
1360  * frame number.
1361  *
1362  * NOTE: frame number should be in clip space, not scene space.
1363  */
1364 /* TODO(sergey): Use r_ prefix for output parameters here. */
1365 ImBuf *BKE_tracking_stabilize_frame(
1366     MovieClip *clip, int framenr, ImBuf *ibuf, float translation[2], float *scale, float *angle)
1367 {
1368   float tloc[2], tscale, tangle;
1369   MovieTracking *tracking = &clip->tracking;
1370   MovieTrackingStabilization *stab = &tracking->stabilization;
1371   ImBuf *tmpibuf;
1372   int width = ibuf->x, height = ibuf->y;
1373   float pixel_aspect = tracking->camera.pixel_aspect;
1374   float mat[4][4];
1375   int filter = tracking->stabilization.filter;
1376   interpolation_func interpolation = NULL;
1377   int ibuf_flags;
1378
1379   if (translation) {
1380     copy_v2_v2(tloc, translation);
1381   }
1382
1383   if (scale) {
1384     tscale = *scale;
1385   }
1386
1387   /* Perform early output if no stabilization is used. */
1388   if ((stab->flag & TRACKING_2D_STABILIZATION) == 0) {
1389     if (translation) {
1390       zero_v2(translation);
1391     }
1392
1393     if (scale) {
1394       *scale = 1.0f;
1395     }
1396
1397     if (angle) {
1398       *angle = 0.0f;
1399     }
1400
1401     return ibuf;
1402   }
1403
1404   /* Allocate frame for stabilization result. */
1405   ibuf_flags = 0;
1406   if (ibuf->rect) {
1407     ibuf_flags |= IB_rect;
1408   }
1409   if (ibuf->rect_float) {
1410     ibuf_flags |= IB_rectfloat;
1411   }
1412
1413   tmpibuf = IMB_allocImBuf(ibuf->x, ibuf->y, ibuf->planes, ibuf_flags);
1414
1415   /* Calculate stabilization matrix. */
1416   BKE_tracking_stabilization_data_get(clip, framenr, width, height, tloc, &tscale, &tangle);
1417   BKE_tracking_stabilization_data_to_mat4(
1418       ibuf->x, ibuf->y, pixel_aspect, tloc, tscale, tangle, mat);
1419
1420   /* The following code visits each nominal target grid position
1421    * and picks interpolated data "backwards" from source.
1422    * thus we need the inverse of the transformation to apply. */
1423   invert_m4(mat);
1424
1425   if (filter == TRACKING_FILTER_NEAREST) {
1426     interpolation = nearest_interpolation;
1427   }
1428   else if (filter == TRACKING_FILTER_BILINEAR) {
1429     interpolation = bilinear_interpolation;
1430   }
1431   else if (filter == TRACKING_FILTER_BICUBIC) {
1432     interpolation = bicubic_interpolation;
1433   }
1434   else {
1435     /* fallback to default interpolation method */
1436     interpolation = nearest_interpolation;
1437   }
1438
1439   TrackingStabilizeFrameInterpolationData data = {
1440       .ibuf = ibuf,
1441       .tmpibuf = tmpibuf,
1442       .mat = mat,
1443       .interpolation = interpolation,
1444   };
1445
1446   ParallelRangeSettings settings;
1447   BLI_parallel_range_settings_defaults(&settings);
1448   settings.use_threading = (tmpibuf->y > 128);
1449   BLI_task_parallel_range(
1450       0, tmpibuf->y, &data, tracking_stabilize_frame_interpolation_cb, &settings);
1451
1452   if (tmpibuf->rect_float) {
1453     tmpibuf->userflags |= IB_RECT_INVALID;
1454   }
1455
1456   if (translation) {
1457     copy_v2_v2(translation, tloc);
1458   }
1459
1460   if (scale) {
1461     *scale = tscale;
1462   }
1463
1464   if (angle) {
1465     *angle = tangle;
1466   }
1467
1468   return tmpibuf;
1469 }
1470
1471 /* Build a 4x4 transformation matrix based on the given 2D stabilization data.
1472  * mat is a 4x4 matrix in homogeneous coordinates, adapted to the
1473  *     final image buffer size and compensated for pixel aspect ratio,
1474  *     ready for direct OpenGL drawing.
1475  *
1476  * TODO(sergey): The signature of this function should be changed. we actually
1477  *               don't need the dimensions of the image buffer. Instead we
1478  *               should consider to provide the pivot point of the rotation as a
1479  *               further stabilization data parameter.
1480  */
1481 void BKE_tracking_stabilization_data_to_mat4(int buffer_width,
1482                                              int buffer_height,
1483                                              float pixel_aspect,
1484                                              float translation[2],
1485                                              float scale,
1486                                              float angle,
1487                                              float r_mat[4][4])
1488 {
1489   /* Since we cannot receive the real pivot point coordinates (API limitation),
1490    * we perform the rotation/scale around the center of frame.
1491    * Then we correct by an additional shift, which was calculated in
1492    * compensate_rotation_center() and "sneaked in" as additional offset
1493    * in the translation parameter. This works, since translation needs to be
1494    * applied after rotation/scale anyway. Thus effectively the image gets
1495    * rotated around the desired pivot point
1496    */
1497   /* TODO(sergey) pivot shouldn't be calculated here, rather received
1498    * as a parameter.
1499    */
1500   float pivot[2];
1501   pivot[0] = 0.5f * pixel_aspect * buffer_width;
1502   pivot[1] = 0.5f * buffer_height;
1503   /* Compose transformation matrix. */
1504   stabilization_data_to_mat4(pixel_aspect, pivot, translation, scale, angle, r_mat);
1505 }