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