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