a78b86b23e1ceeae0a030a46bc27290b2ef43594
2 * localmotion2transform.c
4 * Copyright (C) Georg Martius - January 2013
5 * georg dot martius at web dot de
7 * This file is part of vid.stab video stabilization library
9 * vid.stab is free software; you can redistribute it and/or modify
10 * it under the terms of the GNU General Public License,
11 * as published by the Free Software Foundation; either version 2, or
12 * (at your option) any later version.
14 * vid.stab is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 * GNU General Public License for more details.
19 * You should have received a copy of the GNU General Public License
20 * along with GNU Make; see the file COPYING. If not, write to
21 * the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
25 #include "localmotion2transform.h"
26 #include "transformtype_operations.h"
30 /* #include <sys/time.h> */
31 /* long timeOfDayinMS() { */
32 /* struct timeval t; */
33 /* gettimeofday(&t, 0); */
34 /* return t.tv_sec*1000 + t.tv_usec/1000; */
37 int vsLocalmotions2Transforms(VSTransformData
* td
,
38 const VSManyLocalMotions
* motions
,
39 VSTransformations
* trans
){
40 int len
= vs_vector_size(motions
);
41 assert(trans
->len
==0 && trans
->ts
== 0);
42 trans
->ts
= vs_malloc(sizeof(VSTransform
)*len
);
43 /* long start= timeOfDayinMS(); */
45 if(td
->conf
.storeTransforms
){
46 f
= fopen("global_motions.trf","w");
49 if(td
->conf
.simpleMotionCalculation
==0){
50 for(int i
=0; i
< vs_vector_size(motions
); i
++) {
51 trans
->ts
[i
]=vsMotionsToTransform(td
,VSMLMGet(motions
,i
), f
);
54 for(int i
=0; i
< vs_vector_size(motions
); i
++) {
55 trans
->ts
[i
]=vsSimpleMotionsToTransform(td
->fiSrc
, td
->conf
.modName
,VSMLMGet(motions
,i
));
60 /* long end = timeOfDayinMS(); */
61 /* vs_log_info(td->conf.modName, "Localmotions2Transform (%i) with %i frames took %i ms\n", */
62 /* td->conf.simpleMotionCalculation, len, end-start); */
67 VSArray
vsTransformToArray(const VSTransform
* t
){
68 VSArray a
= vs_array_new(4);
76 VSTransform
vsArrayToTransform(VSArray a
){
77 return new_transform(a
.dat
[0],a
.dat
[1],a
.dat
[2],a
.dat
[3],0,0,0);
80 struct VSGradientDat
{
82 const LocalMotions
* motions
;
83 VSArray missmatches
; // if negative then local motion is ignored
86 double calcTransformQuality(VSArray params
, void* dat
){
87 struct VSGradientDat
* gd
= (struct VSGradientDat
*) dat
;
88 const LocalMotions
* motions
= gd
->motions
;
89 int num_motions
=vs_vector_size(motions
);
90 VSTransform t
= vsArrayToTransform(params
);
93 PreparedTransform pt
= prepare_transform(&t
, &gd
->td
->fiSrc
);
94 int num
= 1; // we start with 1 to avoid div by zero
95 for (int i
= 0; i
< num_motions
; i
++) {
96 if(gd
->missmatches
.dat
[i
]>=0){
97 LocalMotion
* m
= LMGet(motions
,i
);
99 transform_vec_double(&vx
, &vy
, &pt
, (Vec
*)&m
->f
);
100 vx
-= m
->f
.x
; vy
-= m
->f
.y
;
101 double e
= sqr(vx
- m
->v
.x
) + sqr(vy
- m
->v
.y
);
102 gd
->missmatches
.dat
[i
]=e
;
107 // 1 pixel translation missmatch is roughly (with size 500):
108 // alpha=0.11 (degree), zoom=0.2; The zoom is however often much larger, so less penalty.
109 return error
/num
+ fabs(t
.alpha
)/5.0 + fabs(t
.zoom
)/500.0;
112 double intMean(const int* ds
, int len
) {
114 for (int i
= 0; i
< len
; i
++) sum
+= ds
[i
];
118 // only calcates means transform to initialise gradient descent
119 VSTransform
meanMotions(VSTransformData
* td
, const LocalMotions
* motions
){
120 int len
= vs_vector_size(motions
);
121 int* xs
= localmotions_getx(motions
);
122 int* ys
= localmotions_gety(motions
);
123 VSTransform t
= null_transform();
124 if(motions
==0 || len
==0) {
125 t
.extra
= 1; // prob. blank frame or too low contrast, ignore later
128 t
.x
= intMean(xs
,len
);
129 t
.y
= intMean(ys
,len
);
135 /* Disables those fields (mask = -1) whose (miss)quality is high.
136 @param mask: fields masks (<0 means disabled)
137 @param missqualities: measure for each field (larger is worse)
138 @param stddevs: x standard deviations to exclude
139 Both array have to be of the same length.
140 @return number of disabled fields
142 int disableFields(VSArray mask
, VSArray missqualities
, double stddevs
){
143 assert(mask
.len
== missqualities
.len
);
144 // first we throw away those fields that match badely (during motion detection)
145 double mu
= mean(missqualities
.dat
, missqualities
.len
);
146 double sigma
= stddev(missqualities
.dat
, missqualities
.len
, mu
);
147 double thresh
= mu
+ stddevs
* sigma
;
149 for(int i
=0; i
< mask
.len
; i
++){
150 if(missqualities
.dat
[i
]>thresh
){
151 mask
.dat
[i
]=-1.0; // disable field
158 VSTransform
vsMotionsToTransform(VSTransformData
* td
,
159 const LocalMotions
* motions
,
161 VSTransform t
= meanMotions(td
, motions
);
162 if(motions
==0 || vs_vector_size(motions
)==0){
163 if (f
) fprintf(f
,"0 0 0 0 0 %i\n# no fields\n", t
.extra
);
166 VSArray missmatches
= vs_array_new(vs_vector_size(motions
));
167 VSArray params
= vsTransformToArray(&t
);
169 struct VSGradientDat dat
;
170 dat
.motions
= motions
;
172 dat
.missmatches
= missmatches
;
174 // first we throw away those fields that match badely (during motion detection)
175 VSArray matchQualities
= localmotionsGetMatch(motions
);
176 int dis1
=disableFields(missmatches
, matchQualities
, 1.5);
177 vs_array_free(matchQualities
);
180 double ss
[] = {0.2, 0.2, 0.00005, 0.1};
184 // optimize params to minimize transform quality (12 steps per dimension)
185 result
= vsGradientDescent(calcTransformQuality
, params
, &dat
,
186 16, vs_array(ss
,4), 0.01, &residual
);
187 vs_array_free(params
);
188 // now we need to ignore the fields that don't fit well (e.g. moving objects)
189 // cut off everthing above 1 std. dev. for skewed distributions
190 // this will cut off the tail
191 // do this only two times (3 gradient optimizations in total)
192 if((k
==0 && residual
>0.1) || (k
==1 && residual
>20)){
193 dis2
+= disableFields(missmatches
, missmatches
, 1.0);
198 if(td
->conf
.verbose
& VS_DEBUG
)
199 vs_log_info(td
->conf
.modName
, "disabled (%i+%i)/%i,\tresidual: %f (%i)\n",
200 dis1
, dis2
, vs_vector_size(motions
), residual
, k
+1);
201 t
= vsArrayToTransform(result
);
202 vs_array_free(result
);
203 vs_array_free(missmatches
);
204 // check if sufficiently good match was achieved:
205 if(residual
>100){ // test threshold.
209 fprintf(f
,"0 %f %f %f %f %i\n#\t\t\t\t\t %f %i\n", t
.x
, t
.y
, t
.alpha
, t
.zoom
, t
.extra
,
212 if(!td
->conf
.smoothZoom
)
219 /* n-dimensional general purpose gradient descent algorithm */
220 VSArray
vsGradientDescent(double (*eval
)(VSArray
, void*),
221 VSArray params
, void* dat
,
222 int N
, VSArray stepsizes
, double threshold
, double* residual
){
224 double v
= eval(params
, dat
);
225 VSArray x
= vs_array_copy(params
);
226 VSArray grad
= vs_array_new(dim
);
227 assert(stepsizes
.len
== params
.len
);
228 for(int i
=0; i
< N
*dim
&& v
> threshold
; i
++){
230 VSArray x2
= vs_array_copy(x
);
231 double h
= rand()%2 ? 1e-6 : -1e-6;
233 double v2
= eval(x2
, dat
);
234 vs_array_zero(&grad
);
235 grad
.dat
[k
] = (v
- v2
)/h
;
236 vs_array_plus(&x2
, x
, *vs_array_scale(&x2
, grad
, stepsizes
.dat
[k
]));
239 //fprintf(stderr,"+");
243 stepsizes
.dat
[k
]*=1.2; // increase stepsize (4 successful steps will double it)
244 }else{ // overshoot: reduce stepsize and don't do the step
245 //fprintf(stderr,".");
246 stepsizes
.dat
[k
]/=2.0;
249 //if(k==3) fprintf(stderr," ");
252 vs_array_free(stepsizes
);
253 if(residual
!= NULL
) *residual
=v
;
258 /* *** old calculation ***/
260 /* calculates rotation angle for the given transform and
261 * field with respect to the given center-point
263 double vsCalcAngle(const LocalMotion
* lm
, int center_x
, int center_y
){
264 // we better ignore fields that are to close to the rotation center
265 if (abs(lm
->f
.x
- center_x
) + abs(lm
->f
.y
- center_y
) < lm
->f
.size
*2) {
268 // double r = sqrt(lm->f.x*lm->f.x + lm->f.y*lm->f.y);
269 double a1
= atan2(lm
->f
.y
- center_y
, lm
->f
.x
- center_x
);
270 double a2
= atan2(lm
->f
.y
- center_y
+ lm
->v
.y
,
271 lm
->f
.x
- center_x
+ lm
->v
.x
);
272 double diff
= a2
- a1
;
273 return (diff
> M_PI
) ? diff
- 2 * M_PI
: ((diff
< -M_PI
) ? diff
+ 2
279 VSTransform
vsSimpleMotionsToTransform(VSFrameInfo fi
, const char* modName
,
280 const LocalMotions
* motions
){
283 VSTransform t
= null_transform();
284 if(motions
==0) return t
;
285 int num_motions
=vs_vector_size(motions
);
286 double *angles
= (double*) vs_malloc(sizeof(double) * num_motions
);
287 LocalMotion meanmotion
;
292 // calc center point of all remaining fields
293 for (i
= 0; i
< num_motions
; i
++) {
294 center_x
+= LMGet(motions
,i
)->f
.x
;
295 center_y
+= LMGet(motions
,i
)->f
.y
;
297 center_x
/= num_motions
;
298 center_y
/= num_motions
;
301 meanmotion
= cleanmean_localmotions(motions
);
304 if (num_motions
< 6) {
305 // the angle calculation is inaccurate for 5 and less fields
308 for (i
= 0; i
< num_motions
; i
++) {
309 // substract avg and calc angle
310 LocalMotion m
= sub_localmotion(LMGet(motions
,i
),&meanmotion
);
311 angles
[i
] = vsCalcAngle(&m
, center_x
, center_y
);
314 t
.alpha
= -cleanmean(angles
, num_motions
, &min
, &max
);
315 if (max
- min
> 1.0) {
317 vs_log_info(modName
, "too large variation in angle(%f)\n",
322 // compensate for off-center rotation
323 double p_x
= (center_x
- fi
.width
/ 2);
324 double p_y
= (center_y
- fi
.height
/ 2);
325 t
.x
= meanmotion
.v
.x
+ (cos(t
.alpha
) - 1) * p_x
- sin(t
.alpha
) * p_y
;
326 t
.y
= meanmotion
.v
.y
+ sin(t
.alpha
) * p_x
+ (cos(t
.alpha
) - 1) * p_y
;
335 * c-file-style: "stroustrup"
336 * c-file-offsets: ((case-label . *) (statement-case-intro . *))
337 * indent-tabs-mode: nil
338 * c-basic-offset: 2 t
341 * vim: expandtab shiftwidth=2: