uras.c

this program, opticalflow.c, is an implementation of Uras et al. 1988 s motion

#include <math.h>
#include <stdio.h>
#include <fcntl.h>

/* 
           NAME : const.h 
   PARAMETER(S) : none
 
        PURPOSE : Definition of the constants for the 
                  implementation of Uras' motion 
                  detection approach.

         AUTHOR : Steven Beauchemin
             AT : University of Western Ontario
           DATE : September 13 1990

*/

#define TRUE        1
#define FALSE       0

#define NO_REG      0
#define TR1         1
#define TR2         2
#define PI          3.14159

#define H           32
#define SMSK        5
#define REGION_SIZE 8
#define X           316
#define Y           316
#define Z           20

#define DEF_S1      3.0 
#define DEF_S2      1.5
#define DEF_S3      3.0

#define N_HISTO     4
#define N_BINS      100

#define MAX_COND    100000.0

#define NRADIUS     4
#define SKIP        1

#define MAXFLOW      20.0
#define NO_ERROR     0
#define SA_OVERFLOW  1
#define NO_FLOW      2

/* 
           NAME : type.h 
   PARAMETER(S) : none
 
        PURPOSE : Type definitions for images
                  and related data structures.

         AUTHOR : Steven Beauchemin
             AT : University of Western Ontario
           DATE : September 13 1990

*/

typedef struct t_raster {
                          int m, width, height, depth, length,
                              type, maptype, maplength ;
                        } raster_t ;

typedef struct t_beaudet {
                           float g[SMSK][SMSK] ;
                           float f ;
                           int m ;
                         } beaudet_t ;

typedef struct t_kernel {
                          float k[X], f ;
                          int m ;
                        } kernel_t ;

typedef struct t_disp_vect {
                             float x, y ;  
                           } disp_vect_t ;

typedef struct t_param {
                         float fxx, fyy, fxy, ft, fxt, fyt, fx, fy, 
                               discr, gauss, cond ;
                         int err ;
                       } param_t ;

typedef struct t_pos { int i, j ; } pos_t ;

typedef disp_vect_t disp_field512_t[X*Y] ;

typedef param_t param512_t[X*Y] ;

typedef float image512_t[X*Y] ;

typedef struct t_qnode {
                         int             res, sizx, sizy, sizz, ofst, level ;
                         image512_t      *gauss_ptr[Z] ;
                         param512_t      *param_ptr[Z] ;
                         disp_field512_t *flow_ptr ;
                         struct t_qnode  *forth, *back ;
                       } qnode_t, *qnode_ptr_t ;

typedef struct t_histo {
                         float std, avg ;
                         int   freq ;
                       } histo_t ;

/* 
           NAME : str.h
   PARAMETER(S) : none
 
        PURPOSE : Definitions of constants and types for 
                  string manipulation.

         AUTHOR : Steven Beauchemin
             AT : University of Western Ontario
           DATE : June 25 1990
*/

#define STRLENGTH 80

typedef char string[STRLENGTH] ;

/* 
           NAME : extvar.h
   PARAMETER(S) : none
 
        PURPOSE : declaration of external variables; all
                  masks are external.

         AUTHOR : Steven Beauchemin
             AT : University of Western Ontario
           DATE : September 13 1990
*/

beaudet_t Ix, Iy, Ixx, Ixy ;
kernel_t ker1, ker2, ker3, C ;
int KERNEL_X, KERNEL_Y ;

/* 
           NAME : concat(s1,s2,s3) ;
   PARAMETER(S) : s1, s2 : strings to concat;
                      s3 : output string.
 
        PURPOSE : Concats s1 and s2 into s3.

         AUTHOR : Steven Beauchemin
             AT : University of Western Ontario
           DATE : June 25 1990
*/

concat(s1,s2,s3)
string s1, s2, s3 ;

{ string t ;
  int i, j ;

  for (i = 0 ; (t[i] = s1[i]) != '\0' ; i++) ;
  for (j = i ; (t[j] = s2[j - i]) != '\0' ; j++) ;
  t[j] = '\0' ;
  for (i = 0 ; i <= j ; i++) {
    s3[i] = t[i] ;
  }
}

/* 
           NAME : condition(loc,f)
   PARAMETER(S) : loc : location on flow;
                    f : flow pointer.
 
        PURPOSE : returns condition number of estimate at loc in f

         AUTHOR : Steven Beauchemin
             AT : University of Western Ontario
           DATE : November 7 1990
*/

float condition(loc,f)
pos_t loc ;
qnode_ptr_t f ;

{ float cond ;

  if ((loc.i != -1) && (loc.j != -1)) {
    cond = (*f->param_ptr[f->sizz/2])[f->res*loc.i + loc.j].cond ;
  }
  else {
    cond = 0.0 ;
  }
  return(cond) ;
}

/* 
           NAME : convolve(c1,c2,ker1,ker2,n)
   PARAMETER(S) :     c1,c2 : pointers on image cube nodes;
                  ker1,ker2 : kernels  used for 3D convolution;
                          n : number of image frames in cubes.
 
        PURPOSE : Performs a 3D convolution on images of c1 in c2
                  using 2 1D kernels.

         AUTHOR : Steven Beauchemin
             AT : University of Western Ontario
           DATE : March 20 1990

*/

convolve(c1,c2,ker1,ker2,n)
qnode_ptr_t c1, c2 ;
kernel_t ker1, ker2 ;
int n ;

{ image512_t *t ;
  float s ;
  int i, j, k, l, h1, h2 ;

  h2 = ker2.m/2 ;
  
  for (k = h2 ; k < n - h2 ; k++) {
    for (i = 0 ; i < c1->sizy ; i++) {
      for (j = 0 ; j < c1->sizx ; j++) {
        s = 0.0 ;
        for (l = 0 ; l < ker2.m ; l++) {
          s += (*c1->gauss_ptr[l+k-h2])[c1->res*i + j]*ker2.k[l] ;
        }
        (*c2->gauss_ptr[k-h2])[c2->res*i + j] = s/ker2.f ;
      }
    }
  }

  h1 = ker1.m/2 ;

  for (k = h2 ; k < n - h2 ; k++) {

    t = (image512_t *)malloc(sizeof(image512_t)) ;

    for (i = 0 ; i < c2->sizy ; i++) {
      for (j = h1 ; j < c2->sizx - h1 ; j++) {
        s = 0.0 ;
        for (l = 0 ; l < ker1.m ; l++) {
          s += (*c2->gauss_ptr[k-h2])[c2->res*i + j+l-h1]*ker1.k[l] ;
        }
        (*t)[c2->res*i + j] = s/ker1.f ;
      }
    }
  
    for (i = h1 ; i < c2->sizy - h1 ; i++) {
      for (j = 0 ; j < c2->sizx ; j++) {
        s = 0.0 ;
        for (l = 0 ; l < ker1.m ; l++) {
          s += (*t)[c2->res*(i+l-h1) + j]*ker1.k[l] ;
        }
        (*c2->gauss_ptr[k-h2])[c2->res*i + j] = s/ker1.f ;
      }
    }
    free((image512_t *)t) ;
  }
  c2->ofst += h1 ;
}

/* 
           NAME : qnode_ptr_t create_node(l,r,sizx,sizy,sizz,ofst) 
   PARAMETER(S) :   l          : level in the pyramid;
                    r          : resolution at level l;
                    sizx, sizy : image size;
                    sizz       : sequence length;
                    ofst       : offset from image boundaries.
 
        PURPOSE : Creation of a node with 
                  level l.

         AUTHOR : Steven Beauchemin
             AT : University of Western Ontario
           DATE : March 17 1990

*/

qnode_ptr_t create_node(l,r,sizx,sizy,sizz,ofst)
int l, r, sizx, sizy, sizz, ofst ;

{ qnode_ptr_t p ; 
  
  p = (qnode_ptr_t)malloc(sizeof(qnode_t)) ;
  p->res = r ;
  p->sizx = sizx ;
  p->sizy = sizy ;
  p->sizz = sizz ;
  p->ofst = ofst ;
  p->level = l ;
  return(p) ;
}

/* 
           NAME : delete_node(p,h,q)
   PARAMETER(S) : p : pointer on the node to be deleted;
                  h : head list pointer;
                  q : tail list pointer. 

        PURPOSE : Deletion of the node pointed by p.

         AUTHOR : Steven Beauchemin
             AT : University of Western Ontario
           DATE : March 17 1990

*/

delete_node(p,h,q)
qnode_ptr_t p, *h, *q ;

{ qnode_ptr_t s, t ;

  s = p->back ;
  t = p->forth ;
  if (t != (qnode_ptr_t)NULL) {
    t->back = s ;
  }
  else {
    *q = s ;
  }
  if (s != (qnode_ptr_t)NULL) {
    s->forth = t ;
  }
  else {
    *h = t ;
  }
  free(p) ;
}

/* 
           NAME : determinant(loc,f)
   PARAMETER(S) : loc : location on flow;
                    f : flow pointer.
 
        PURPOSE : returns determinant of estimate at loc in f

         AUTHOR : Steven Beauchemin
             AT : University of Western Ontario
           DATE : November 7 1990
*/

float determinant(loc,f)
pos_t loc ;
qnode_ptr_t f ;

{ float det ;

  if ((loc.i != -1) && (loc.j != -1)) {
    det = (*f->param_ptr[f->sizz/2])[f->res*loc.i + loc.j].gauss ;
  }
  else {
    det = 0.0 ;
  }
  return(det) ;
}

/* 
           NAME : dt(c,q,x,y)
   PARAMETER(S) :      c : image cube node;
                       q : image parameter node ;
                    x, y : image location for computation.
 
        PURPOSE : computes 1st order derivative with respect
                  to time at x,y for middle image of cube c.

         AUTHOR : Steven Beauchemin
             AT : University of Western Ontario
           DATE : October 1 1990
*/

float dt(c,q,x,y)
qnode_ptr_t c, q ;
int x, y ;

{ extern kernel_t C ;
  float d ;
  int i, k, h ;

  h = C.m/2 ;
  k = c->sizz/2 ;

  d = 0.0 ;
 for (i = -h ; i <= h ; i++) {
    d += (*c->gauss_ptr[i+h])[c->res*x + y]*C.k[i+h] ;
  }
  d /= C.f ;
  return(d) ;
}

/* 
           NAME : overflow(x,y,sizx,sizy,ofst,h,l)
   PARAMETER(S) : x,y       : image coordinates ;
                  sizx,sizy : image size;
                  ofst      : offset from image boundaries;
                    h       : even half size of odd masks applied;
                    l       : number of mask applications.
 
        PURPOSE : Returns TRUE if x,y in image, 
                  FALSE otherwise.

         AUTHOR : Steven Beauchemin
             AT : University of Western Ontario
           DATE : July 23 1990
*/

int overflow(x,y,sizx,sizy,ofst,h,l) 
int x, y, sizx, sizy, ofst, h, l ;

{
  return((x < ofst + h*l) || 
         (x >= sizx - (ofst + h*l)) || 
         (y < ofst + h*l) || 
         (y >= sizy - (ofst + h*l))) ;
}

/* 
           NAME : dux(p,x,y)
   PARAMETER(S) : p : flow node;
                x,y : flow location for computation.
 
        PURPOSE : computes 1st order partial derviative with respect
                  to x at x,y for component u of flow p.

         AUTHOR : Steven Beauchemin
             AT : University of Western Ontario
           DATE : October 4 1990
*/

float dux(p,x,y)
qnode_ptr_t p ;
int x, y ;

{ extern beaudet_t Ix ;
  float d ;
  int i, j, h ;

  h = Ix.m/2 ;
  d = 0.0 ;

  if (!overflow(x,y,p->sizy,p->sizx,p->ofst,h,3)) {
    for (i = -h ; i <= h ; i++) {
      for (j = -h ; j <= h ; j++) {
        d = d + (*p->flow_ptr)[p->res*(x+i) + y+j].x*Ix.g[i+h][j+h] ;
      }
    }
  }
  return(d/Ix.f) ;
}

/* 
           NAME : duy(p,x,y)
   PARAMETER(S) : p : flow node;
                x,y : flow location for computation.
 
        PURPOSE : computes 1st order partial derviative with respect
                  to y at x,y for component u of flow p.

         AUTHOR : Steven Beauchemin
             AT : University of Western Ontario
           DATE : October 4 1990
*/

float duy(p,x,y)
qnode_ptr_t p ;
int x, y ;

{ extern beaudet_t Ix ;
  float d ;
  int i, j, h ;
  
  h = Ix.m/2 ;
  d = 0.0 ;
  
  if (!overflow(x,y,p->sizy,p->sizx,p->ofst,h,3)) {
    for (i = -h ; i <= h ; i++) {
      for (j = -h ; j <= h ; j++) {
        d = d + (*p->flow_ptr)[p->res*(x+i) + y+j].x*Ix.g[j+h][i+h] ;
      }
    }
  }
  return(d/Ix.f) ;
}

/* 
           NAME : dvx(p,x,y)
   PARAMETER(S) : p : flow node;
                x,y : flow location for computation.
 
        PURPOSE : computes 1st order partial derviative with respect
                  to x at x,y for component v of flow p.

         AUTHOR : Steven Beauchemin
             AT : University of Western Ontario
           DATE : October 4 1990
*/

float dvx(p,x,y)
qnode_ptr_t p ;
int x, y ;

{ extern beaudet_t Ix ;
  float d ;
  int i, j, h ;

  h = Ix.m/2 ;
  d = 0.0 ; 

  if (!overflow(x,y,p->sizy,p->sizx,p->ofst,h,3)) {
    for (i = -h ; i <= h ; i++) {
      for (j = -h ; j <= h ; j++) {
        d = d + (*p->flow_ptr)[p->res*(x+i) + y+j].y*Ix.g[i+h][j+h] ;
      }