Skip to content

Commit

Permalink
gradient algorithm for motion calculation works now.
Browse files Browse the repository at this point in the history
  • Loading branch information
Georg Martius committed Oct 7, 2013
1 parent 17a76a3 commit 448888b
Show file tree
Hide file tree
Showing 8 changed files with 193 additions and 116 deletions.
252 changes: 149 additions & 103 deletions src/localmotion2transform.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,92 +28,29 @@
#include <string.h>

int vsLocalmotions2TransformsSimple(VSTransformData* td,
const VSManyLocalMotions* motions,
VSTransformations* trans ){
const VSManyLocalMotions* motions,
VSTransformations* trans ){
return vsLocalmotions2Transforms(td,motions,trans);
}

int vsLocalmotions2Transforms(VSTransformData* td,
const VSManyLocalMotions* motions,
VSTransformations* trans ){
int i;
int len = vs_vector_size(motions);
assert(trans->len==0 && trans->ts == 0);
trans->ts = vs_malloc(sizeof(VSTransform)*len );
for(i=0; i< vs_vector_size(motions); i++) {
// trans->ts[i]=vsSimpleMotionsToTransform(td,VSMLMGet(motions,i));
// TODO: use a flag in td
trans->ts[i]=vsMotionsToTransform(td,VSMLMGet(motions,i));
// vsStoreLocalmotions(stderr,VSMLMGet(motions,i));
// storeTransform(stderr,&trans->ts[i]);
}
trans->len=len;
return VS_OK;
}

/* calculates rotation angle for the given transform and
* field with respect to the given center-point
*/
double vsCalcAngle(const LocalMotion* lm, int center_x, int center_y){
// we better ignore fields that are to close to the rotation center
if (abs(lm->f.x - center_x) + abs(lm->f.y - center_y) < lm->f.size*2) {
return 0;
} else {
// double r = sqrt(lm->f.x*lm->f.x + lm->f.y*lm->f.y);
double a1 = atan2(lm->f.y - center_y, lm->f.x - center_x);
double a2 = atan2(lm->f.y - center_y + lm->v.y,
lm->f.x - center_x + lm->v.x);
double diff = a2 - a1;
return (diff > M_PI) ? diff - 2 * M_PI : ((diff < -M_PI) ? diff + 2
* M_PI : diff);
}
}


VSTransform vsSimpleMotionsToTransform(VSTransformData* td,
const LocalMotions* motions){
int center_x = 0;
int center_y = 0;
VSTransform t = null_transform();
if(motions==0) return t;
int num_motions=vs_vector_size(motions);
double *angles = (double*) vs_malloc(sizeof(double) * num_motions);
LocalMotion meanmotion;
int i;
if(num_motions < 1)
return t;

// calc center point of all remaining fields
for (i = 0; i < num_motions; i++) {
center_x += LMGet(motions,i)->f.x;
center_y += LMGet(motions,i)->f.y;
}
center_x /= num_motions;
center_y /= num_motions;

// cleaned mean
meanmotion = cleanmean_localmotions(motions);

// figure out angle
if (num_motions < 6) {
// the angle calculation is inaccurate for 5 and less fields
t.alpha = 0;
} else {
for (i = 0; i < num_motions; i++) {
// substract avg and calc angle
LocalMotion m = sub_localmotion(LMGet(motions,i),&meanmotion);
angles[i] = vsCalcAngle(&m, center_x, center_y);
if(td->conf.simpleMotionCalculation){
for(i=0; i< vs_vector_size(motions); i++) {
trans->ts[i]=vsSimpleMotionsToTransform(td,VSMLMGet(motions,i));
}
double min, max;
t.alpha = -cleanmean(angles, num_motions, &min, &max);
if (max - min > td->maxAngleVariation) {
t.alpha = 0;
vs_log_info(td->conf.modName, "too large variation in angle(%f)\n",
max-min);
}else{
for(i=0; i< vs_vector_size(motions); i++) {
trans->ts[i]=vsMotionsToTransform(td,VSMLMGet(motions,i));
}
}
vs_free(angles);
// compensate for off-center rotation
double p_x = (center_x - td->fiSrc.width / 2);
double p_y = (center_y - td->fiSrc.height / 2);
t.x = meanmotion.v.x + (cos(t.alpha) - 1) * p_x - sin(t.alpha) * p_y;
t.y = meanmotion.v.y + sin(t.alpha) * p_x + (cos(t.alpha) - 1) * p_y;

return t;
trans->len=len;
return VS_OK;
}


Expand All @@ -131,11 +68,10 @@ VSTransform vsArrayToTransform(VSArray a){
return new_transform(a.dat[0],a.dat[1],a.dat[2],a.dat[3],0,0,0);
}

double sqr(double x){ return x*x; }

struct VSGradientDat {
VSTransformData* td;
const LocalMotions* motions;
VSArray missmatches; // if negative then local motion is ignored
};

double calcTransformQuality(VSArray params, void* dat){
Expand All @@ -150,42 +86,78 @@ double calcTransformQuality(VSArray params, void* dat){
double zsin_a = z*sin(t.alpha); // scaled sin
double c_x = gd->td->fiSrc.width / 2;
double c_y = gd->td->fiSrc.height / 2;

int num = 1; // we start with 1 to avoid div by zero
for (int i = 0; i < num_motions; i++) {
LocalMotion* m = LMGet(motions,i);
double x = m->f.x-c_x;
double y = m->f.y-c_y;
double vx = zcos_a * x + zsin_a * y + t.x - x;
double vy = -zsin_a * x + zcos_a * y + t.y - y;
error += sqr(vx - m->v.x) + sqr(vy - m->v.y);
if(gd->missmatches.dat[i]>=0){
LocalMotion* m = LMGet(motions,i);
double x = m->f.x-c_x;
double y = m->f.y-c_y;
double vx = zcos_a * x + zsin_a * y + t.x - x;
double vy = -zsin_a * x + zcos_a * y + t.y - y;
double e = sqr(vx - m->v.x) + sqr(vy - m->v.y);
gd->missmatches.dat[i]=e;
error += e;
num++;
}
}
return error/(num_motions+1) + t.alpha*t.alpha + t.zoom*t.zoom/10000.0;
return error/num + t.alpha*t.alpha + t.zoom*t.zoom/10000.0;
}

double int_mean(const int* ds, int len) {
double sum=0;
for (int i = 0; i < len; i++) sum += ds[i];
return sum / len;
}

// only calcates means transform to initialise gradient descent
VSTransform meanMotions(VSTransformData* td, const LocalMotions* motions){
int len = vs_vector_size(motions);
int* xs = localmotions_getx(motions);
int* ys = localmotions_gety(motions);
VSTransform t = null_transform();
if(motions==0 || len==0) return t;
t.x = int_mean(xs,len);
t.y = int_mean(ys,len);
vs_free(xs);
vs_free(ys);
return t;
}

VSTransform vsMotionsToTransform(VSTransformData* td,
const LocalMotions* motions){
VSTransform t = vsSimpleMotionsToTransform(td, motions);
//VSTransform t = null_transform();
if(motions==0) return t;

VSTransform t = meanMotions(td, motions);
if(motions==0 || vs_vector_size(motions)==0) return t;

VSArray missmatches = vs_array_new(vs_vector_size(motions));
VSArray params = vsTransformToArray(&t);
double residual;
struct VSGradientDat dat;
dat.motions = motions;
dat.td = td;
dat.missmatches = missmatches;

double ss[] = {0.2, 0.2, 0.0001, 0.1};
VSArray result = vsGradientDescent(calcTransformQuality, params, &dat,
20, vs_array(ss,4), 1e-15, &residual);
double ss[] = {0.2, 0.2, 0.00005, 0.1};
// optimize params to minimize transform quality (12 steps per dimension)
VSArray result1 = vsGradientDescent(calcTransformQuality, params, &dat,
12, vs_array(ss,4), 1e-15, &residual);
vs_array_free(params);

// now we need to ignore the fields that don't fit well (moving objects)
// TODO
vs_log_info(td->conf.modName, "residual(%f)\n",residual);
// now we need to ignore the fields that don't fit well (e.g. moving objects)
// cut off everthing above 1 std. dev. for skewed distributions this will cut off the tail
double mean_ = mean(missmatches.dat, missmatches.len);
double stddev_ = stddev(missmatches.dat, missmatches.len, mean_);
double thresh = mean_+stddev_;
for(int i=0; i< missmatches.len; i++){
if(missmatches.dat[i]>thresh) missmatches.dat[i]=-1.0; // disable field
}
// and do a gradient step again
VSArray result = vsGradientDescent(calcTransformQuality, result1, &dat,
16, vs_array(ss,4), 1e-15, &residual);
vs_array_free(result1);

if(td->conf.verbose & VS_DEBUG) vs_log_info(td->conf.modName, "residual(%f)\n",residual);

t = vsArrayToTransform(result);
// storeVSTransform(stderr, &t);
vs_array_free(result);
return t;
}
Expand All @@ -211,18 +183,18 @@ VSArray vsGradientDescent(double (*eval)(VSArray, void*),
grad.dat[k] = (v - v2)/h;
vs_array_plus(&x2, x, *vs_array_scale(&x2, grad, stepsizes.dat[k]));
v2 = eval(x2, dat);
// printf("step: (%lf,%lf) :%g \t(%g)\t [%i:%g]\n",x.dat[0],x.dat[1], v, v-v2,k,grad.dat[k]);
if(v2 < v){
fprintf(stderr,"+");
//fprintf(stderr,"+");
vs_array_free(x);
x = x2;
v = v2;
stepsizes.dat[k]*=1.2; // increase stepsize (4 successful steps will double it)
}else{ // overshoot: reduce stepsize and don't do the step
//fprintf(stderr,".");
stepsizes.dat[k]/=2.0;
vs_array_free(x2);
fprintf(stderr,".");
}
//if(k==3) fprintf(stderr," ");
}
vs_array_free(grad);
vs_array_free(stepsizes);
Expand All @@ -231,6 +203,80 @@ VSArray vsGradientDescent(double (*eval)(VSArray, void*),
}


/* *** old calculation ***/

/* calculates rotation angle for the given transform and
* field with respect to the given center-point
*/
double vsCalcAngle(const LocalMotion* lm, int center_x, int center_y){
// we better ignore fields that are to close to the rotation center
if (abs(lm->f.x - center_x) + abs(lm->f.y - center_y) < lm->f.size*2) {
return 0;
} else {
// double r = sqrt(lm->f.x*lm->f.x + lm->f.y*lm->f.y);
double a1 = atan2(lm->f.y - center_y, lm->f.x - center_x);
double a2 = atan2(lm->f.y - center_y + lm->v.y,
lm->f.x - center_x + lm->v.x);
double diff = a2 - a1;
return (diff > M_PI) ? diff - 2 * M_PI : ((diff < -M_PI) ? diff + 2
* M_PI : diff);
}
}


VSTransform vsSimpleMotionsToTransform(VSTransformData* td,
const LocalMotions* motions){
int center_x = 0;
int center_y = 0;
VSTransform t = null_transform();
if(motions==0) return t;
int num_motions=vs_vector_size(motions);
double *angles = (double*) vs_malloc(sizeof(double) * num_motions);
LocalMotion meanmotion;
int i;
if(num_motions < 1)
return t;

// calc center point of all remaining fields
for (i = 0; i < num_motions; i++) {
center_x += LMGet(motions,i)->f.x;
center_y += LMGet(motions,i)->f.y;
}
center_x /= num_motions;
center_y /= num_motions;

// cleaned mean
meanmotion = cleanmean_localmotions(motions);

// figure out angle
if (num_motions < 6) {
// the angle calculation is inaccurate for 5 and less fields
t.alpha = 0;
} else {
for (i = 0; i < num_motions; i++) {
// substract avg and calc angle
LocalMotion m = sub_localmotion(LMGet(motions,i),&meanmotion);
angles[i] = vsCalcAngle(&m, center_x, center_y);
}
double min, max;
t.alpha = -cleanmean(angles, num_motions, &min, &max);
if (max - min > 1.0) {
t.alpha = 0;
vs_log_info(td->conf.modName, "too large variation in angle(%f)\n",
max-min);
}
}
vs_free(angles);
// compensate for off-center rotation
double p_x = (center_x - td->fiSrc.width / 2);
double p_y = (center_y - td->fiSrc.height / 2);
t.x = meanmotion.v.x + (cos(t.alpha) - 1) * p_x - sin(t.alpha) * p_y;
t.y = meanmotion.v.y + sin(t.alpha) * p_x + (cos(t.alpha) - 1) * p_y;

return t;
}



/*
* Local variables:
Expand Down
11 changes: 9 additions & 2 deletions src/localmotion2transform.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,10 +32,17 @@


/** converts for each frame the localmotions into a transform
@deprecated!
*/
int vsLocalmotions2TransformsSimple(VSTransformData* td,
const VSManyLocalMotions* motions,
VSTransformations* trans );
const VSManyLocalMotions* motions,
VSTransformations* trans );

/** converts for each frame the localmotions into a transform
*/
int vsLocalmotions2Transforms(VSTransformData* td,
const VSManyLocalMotions* motions,
VSTransformations* trans );

/** calculates rotation angle for the given transform and
* field with respect to the given center-point
Expand Down
4 changes: 1 addition & 3 deletions src/transform.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -59,6 +59,7 @@ VSTransformConfig vsTransformGetDefaultConfig(const char* modName){
conf.interpolType = VS_BiLinear;
conf.verbose = 0;
conf.modName = modName;
conf.simpleMotionCalculation = 0;
return conf;
}

Expand Down Expand Up @@ -89,9 +90,6 @@ int vsTransformDataInit(VSTransformData* td, const VSTransformConfig* conf,
vsFrameNull(&td->destbuf);
vsFrameNull(&td->dest);

td->maxAngleVariation = 1;
td->rotationThreshhold = 0.25/(180/M_PI);

if (td->conf.maxShift > td->fiDest.width/2)
td->conf.maxShift = td->fiDest.width/2;
if (td->conf.maxShift > td->fiDest.height/2)
Expand Down
9 changes: 3 additions & 6 deletions src/transform.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -86,12 +86,14 @@ typedef struct _VSTransformConfig {
VSBorderType crop; // 1: black bg, 0: keep border from last frame(s)
int invert; // 1: invert transforms, 0: nothing
double zoom; // percentage to zoom: 0->no zooming 10:zoom in 10%
int optZoom; // 1: determine optimal zoom, 0: nothing
int optZoom; // 2: optimal adaptive zoom 1: optimal static zoom, 0: nothing
VSInterpolType interpolType; // type of interpolation: 0->Zero,1->Lin,2->BiLin,3->Sqr
int maxShift; // maximum number of pixels we will shift
double maxAngle; // maximum angle in rad
const char* modName; // module name (used for logging)
int verbose; // level of logging
// if 1 then the simple but fast method to termine the global motion is used
int simpleMotionCalculation;
} VSTransformConfig;

typedef struct _VSTransformData {
Expand All @@ -113,11 +115,6 @@ typedef struct _VSTransformData {
/* Options */
VSTransformConfig conf;

/* maximal difference in angles of fields */
double maxAngleVariation;
/* threshhold below which no rotation is performed */
double rotationThreshhold;

int initialized; // 1 if initialized and 2 if configured
} VSTransformData;

Expand Down
2 changes: 1 addition & 1 deletion src/transformfloat.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -217,7 +217,7 @@ int _FLT(transformPacked)(VSTransformData* td, VSTransform t)
*/
int channels = td->fiSrc.bytesPerPixel;
/* All channels */
if (fabs(t.alpha) > td->rotationThreshhold) {
if (fabs(t.alpha) > 0.1*M_PI/180.0) { // 0.1 deg
for (x = 0; x < td->fiDest.width; x++) {
for (y = 0; y < td->fiDest.height; y++) {
float x_d1 = (x - c_d_x);
Expand Down
Loading

0 comments on commit 448888b

Please sign in to comment.