uras.c

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

    (*ker).k[0] = 1.0 ;
    (*ker).f = 1.0 ;
  }
  (*ker).m = h ; 
}

/* 
           NAME : init_beaudet(Ix,Ixx,Ixy)
   PARAMETER(S) :  Ix : operator for 1st order x derivative;
                  Ixx : operator for 2nd order x derivative;
                  Ixy : operator for partial x and y derivatives.
 
        PURPOSE : Initialize the image operators defined by Beaudet
                  to obtain numerical derivatives.

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

init_beaudet(Ix,Ixx,Ixy)
beaudet_t *Ix, *Ixx, *Ixy ;

{
  (*Ix).g[0][0] = -2 ; (*Ix).g[1][0] = -1 ; (*Ix).g[2][0] = 0 ; 
  (*Ix).g[3][0] = 1 ;  (*Ix).g[4][0] = 2 ;
  (*Ix).g[0][1] = -2 ; (*Ix).g[1][1] = -1 ; (*Ix).g[2][1] = 0 ; 
  (*Ix).g[3][1] = 1 ;  (*Ix).g[4][1] = 2 ;
  (*Ix).g[0][2] = -2 ; (*Ix).g[1][2] = -1 ; (*Ix).g[2][2] = 0 ; 
  (*Ix).g[3][2] = 1 ;  (*Ix).g[4][2] = 2 ;
  (*Ix).g[0][3] = -2 ; (*Ix).g[1][3] = -1 ; (*Ix).g[2][3] = 0 ; 
  (*Ix).g[3][3] = 1 ;  (*Ix).g[4][3] = 2 ;
  (*Ix).g[0][4] = -2 ; (*Ix).g[1][4] = -1 ; (*Ix).g[2][4] = 0 ; 
  (*Ix).g[3][4] = 1 ;  (*Ix).g[4][4] = 2 ;
  (*Ix).f = 50.0 ;
  (*Ix).m = 5 ;
      
  (*Ixx).g[0][0] = 2 ;  (*Ixx).g[1][0] = -1 ;  (*Ixx).g[2][0] = -2 ; 
  (*Ixx).g[3][0] = -1 ;  (*Ixx).g[4][0] = 2 ;
  (*Ixx).g[0][1] = 2 ; (*Ixx).g[1][1] = -1 ; (*Ixx).g[2][1] = -2; 
  (*Ixx).g[3][1] = -1 ; (*Ixx).g[4][1] = 2 ;
  (*Ixx).g[0][2] = 2 ; (*Ixx).g[1][2] = -1 ; (*Ixx).g[2][2] = -2 ; 
  (*Ixx).g[3][2] = -1 ; (*Ixx).g[4][2] = 2 ;
  (*Ixx).g[0][3] = 2 ; (*Ixx).g[1][3] = -1 ; (*Ixx).g[2][3] = -2 ; 
  (*Ixx).g[3][3] = -1 ; (*Ixx).g[4][3] = 2 ;
  (*Ixx).g[0][4] = 2 ;  (*Ixx).g[1][4] = -1 ;  (*Ixx).g[2][4] = -2 ; 
  (*Ixx).g[3][4] = -1 ;  (*Ixx).g[4][4] = 2 ;
  (*Ixx).f = 35.0 ;
  (*Ixx).m = 5 ;

  (*Ixy).g[0][0] = 4 ;  (*Ixy).g[1][0] = 2 ;  (*Ixy).g[2][0] = 0 ; 
  (*Ixy).g[3][0] = -2 ; (*Ixy).g[4][0] = -4 ;
  (*Ixy).g[0][1] = 2 ;  (*Ixy).g[1][1] = 1 ;  (*Ixy).g[2][1] = 0 ; 
  (*Ixy).g[3][1] = -1 ; (*Ixy).g[4][1] = -2 ;
  (*Ixy).g[0][2] = 0 ;  (*Ixy).g[1][2] = 0 ;  (*Ixy).g[2][2] = 0 ; 
  (*Ixy).g[3][2] = 0 ;  (*Ixy).g[4][2] = 0 ;
  (*Ixy).g[0][3] = -2 ; (*Ixy).g[1][3] = -1 ; (*Ixy).g[2][3] = 0 ; 
  (*Ixy).g[3][3] = 1 ;  (*Ixy).g[4][3] = 2 ;
  (*Ixy).g[0][4] = -4 ; (*Ixy).g[1][4] = -2 ; (*Ixy).g[2][4] = 0 ; 
  (*Ixy).g[3][4] = 2 ;  (*Ixy).g[4][4] = 4 ;
  (*Ixy).f = 100.0 ;
  (*Ixy).m = 5 ;
}

/* 
           NAME : init_central(C)
   PARAMETER(S) : C : 1-D mask contraining 4-point central difference
                      factors.
 
        PURPOSE : Inits a 4-point central difference mask. 

         AUTHOR : Steven Beauchemin
             AT : University of Western Ontario
           DATE : April 25 1992
*/

init_central(C)
kernel_t *C ;

{
  (*C).k[0] = -1.0 ; (*C).k[1] = 8.0 ; 
  (*C).k[2] = 0.0 ; (*C).k[3] = -8.0 ; (*C).k[4] = 1.0 ;
  (*C).m = 5 ;
  (*C).f = 12.0 ;
}

/* 
           NAME : init_cube(p,n) ;
   PARAMETER(S) : p : pointer on an image cube node;
                  n : number of image frames to allocate.
 
        PURPOSE : allocates n images to p->gauss_ptr[i].

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

init_cube(p,n)
qnode_ptr_t p ;
int n ;

{ int i ;

  for (i = 0 ; i < n ; i++) { 
    p->gauss_ptr[i] = (image512_t *)malloc(sizeof(image512_t)) ;
  }
}

/* 
           NAME : init_flow(p,n) ;
   PARAMETER(S) : p : pointer on an parameter cube node;
                  n : number of parameter frames to allocate.
 
        PURPOSE : allocates n parameter frames to p->param_ptr[i]
                  and a flow field to p->flow_ptr.

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

init_flow(p,n)
qnode_ptr_t p ;
int n ;

{ int i ;

  p->flow_ptr = (disp_field512_t *)malloc(sizeof(disp_field512_t)) ;
  for (i = 0 ; i < n ; i++) { 
    p->param_ptr[i] = (param512_t *)malloc(sizeof(param512_t)) ;
  }
}

/* 
           NAME : init_list(h,q) 
   PARAMETER(S) : h : head list pointer;
                  q : tail list pointer.
 
        PURPOSE : Initialization of the head
                  and tail pointers to NULL.

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

*/

init_list(h,q)
qnode_ptr_t *h, *q ;

{ 
  *h = (qnode_ptr_t)NULL ;
  *q = (qnode_ptr_t)NULL ;
}

/* 
           NAME : insert_node(p,h,q) 
   PARAMETER(S) : p : pointer on the node to insert;
                  h : head pointer of the list;
                  q : tail pointer of the list.
 
        PURPOSE : Insertion of the node p in a doubly 
                  linked list in order of the level
                  numbers (p->level).

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

*/

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

{ qnode_ptr_t s, t ;

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


/* 
           NAME : itoa(n,s)
   PARAMETER(S) : n : integer value ;
                  s : output string.
 
        PURPOSE : Converts integer into string.

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

itoa(n,s)
int n ;
string s ; 

{ int i, j, digit, neg ;
  string t ;

  if (n < 0) { 
    neg = TRUE ;
  }
  else {
    neg = FALSE ;
  }
  if (n == 0) {
    s[0] = '0' ;
    s[1] = '\0' ;
  }
  else {
    i = 0 ;
    while (abs(n) > 0) {
      digit = abs(n) % 10 ;
      t[i] = (char)((int)'0' + digit) ;
      n = n / 10 ;
      i++ ;
    }
    if (neg) {
      t[i] = '-' ; 
      i++ ;
    }
    for (j = 0 ; j < i ; j++) {
      s[j] = t[i - j - 1] ;
    }
    s[i] = '\0' ;
  }
}

/* 
           NAME : mag(u)
   PARAMETER(S) : u : 2D flow vector.
 
        PURPOSE : returns the magnitude of u.

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

float mag(u)
disp_vect_t u ;

{ 
  return(sqrt(pow(u.x,2.0) + pow(u.y,2.0))) ;
}

/* 
           NAME : min_cond(s,f,n,m,trsh,sf)
   PARAMETER(S) : s : 1D array containing at most n image flow locations;
                  f : image parameter node;
                  n : maximum number of image locations in s;
                  m : actual number of image locations in s;
               trsh : threshold for condition number;
                 sf : thresholding boolean value.
 
        PURPOSE : determines among the locations in s, which one has a
                  minimal condition number.

         AUTHOR : Steven Beauchemin
             AT : University of Western Ontario
           DATE : April 4 1992
*/

pos_t min_cond(s,f,n,m,trsh,sf)
pos_t s[X] ;
qnode_ptr_t f ;
int n, m ;
float trsh ;
int sf ;

{ pos_t loc, err_loc ;
  float min_c, gauss ;
  int l ;

  min_c = MAX_COND + 1 ;
  loc.i = -1 ;
  loc.j = -1 ;
  err_loc.i = -1 ;
  err_loc.j = -1 ;
  if (m < n) {
    n = m ;
  }
  for (l = 0 ; l < n ; l++) {
    if ((*f->param_ptr[f->sizz/2])[f->res*s[l].i + s[l].j].cond < min_c) {
      min_c = (*f->param_ptr[f->sizz/2])[f->res*s[l].i + s[l].j].cond ;
      gauss = (*f->param_ptr[f->sizz/2])[f->res*s[l].i + s[l].j].gauss ;
      loc = s[l] ;
    }
  }
  if (sf) {
    if (gauss > trsh) {
      return(loc) ;
    }
    else {
      return(err_loc) ;
    }
  }
  else {
    return(loc) ;
  }
}

/* 
           NAME : sort(discr_val,sample,num,n)
   PARAMETER(S) : discr_val : discriminant values for image positions
                              contained in sample;
                  sample    : image positions to be sorted with respect
                              to their discriminant values;
                  num       : number of elements in sample;
                  n         : image region size.
 
        PURPOSE : sorts the n first elements of sample[1..m] with respect
                  to the discriminant value.

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

sort(discr_val,sample,num,n) 
float discr_val[X] ;
pos_t sample[X] ;
int num, n ;

{ pos_t ts ; 
  float td ;
  int i, j, index ;

  if (num < n) {
    n = num ;
  }
  for (i = 0 ; i < n - 1 ; i++) {
    index = i ;
    for (j = index + 1 ; j < num ; j++) {
      if (discr_val[index] > discr_val[j]) {
        index = j ;
      }
    }
    td = discr_val[index] ;
    discr_val[index] = discr_val[i] ;
    discr_val[i] = td ;
   
    ts = sample[index] ;
    sample[index] = sample[i] ;
    sample[i] = ts ;
  }
}

/* 
           NAME : min_discr(f,i,j,n,sample,num)
   PARAMETER(S) :   f      : image parameter node;
                  i,j      : upper left corner location of a n*n image area;
                    n      : image area size;
                    sample : valid image locations of image area;
                    num    : number of image locations in sample.
 
        PURPOSE : extracts the n image locations from the n*n image area
                  for which the discriminant value is minimal.

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

min_discr(f,i,j,n,sample,num)
qnode_ptr_t f ;
int i, j, n ;
pos_t sample[X] ;
int *num ;

{ float discr_val[X] ;
  int k, l ;

  *num = 0 ; 
  for (k = 0 ; k < n ; k++) {
    for (l = 0 ; l < n ; l++) {
      if ((*f->param_ptr[f->sizz/2])[f->res*(i+k) + j+l].err == NO_ERROR) {
        sample[*num].i = i + k ;
        sample[*num].j = j + l ;
        discr_val[*num] = (*f->param_ptr[f->sizz/2])[f->res*(i+k) + j+l].discr ;
        (*num)++ ;
      }
    }
  }
  sort(discr_val,sample,*num,n) ;
}

/* 
           NAME : norm(V,n)
   PARAMETER(S) : V : vector;
                  n : length of V.
 
        PURPOSE : returns L2 norm of V

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

float norm(V,n) 
float V[] ;
int n ;

{ int i ;
  float sum ; 

  sum = 0.0 ;
  for (i = 0 ; i < n ; i++) {
    sum += (V[i]*V[i]) ;
  }
  return(sqrt(sum)) ;
}

/* 
           NAME : pgetrast(fn,hd,bf,sx,sy,r) 
   PARAMETER(S) : fn    : filename;
                  hd    : image header (raster file);
                  bf    : Pointer on a 2D array containing the image;
                  sx,sy : image size to read;
                   r    : resolution of the array.

        PURPOSE : Reads rasterfile specified
                  by fn into bf.

         AUTHOR : Steven Beauchemin
             AT : University of Western Ontario
           DATE : October 27 1989, modified March 20 1990.

*/

pgetrast(fn,hd,bf,sx,sy,r) 
char *fn ;
unsigned char hd[H] ;
float bf[X*Y] ;
int sx, sy, r ;

{ int fd, i, j ;
  unsigned char c ;

  if ((fd = open(fn,O_RDONLY)) > 0) { ;
    if (read(fd,hd,H) == H) { 
      for (i = 0 ; i < sy ; i++) {
        for (j = 0 ; j < sx ; j++) {
          if (read(fd,&c,sizeof(c)) != 1) {
            error(4) ;
          }
          bf[r*i + j] = (float)c ;
        }
      }
    }
    else {
      error(5) ;
    }
  }
  else {
    error(6) ;
  }
  close(fd) ;
}

/* 
           NAME : Bpgetrast(fn,bf,sx,sy,r) 
   PARAMETER(S) : fn    : filename;
                  bf    : Pointer on a 2D array containing the image;
                  sx,sy : image size to read;
                   r    : resolution of the array.

        PURPOSE : Reads rasterfile without header specified
                  by fn into bf.

         AUTHOR : Steven Beauchemin
             AT : University of Western Ontario
           DATE : October 27 1989, modified March 20 1990.

*/

Bpgetrast(fn,bf,sx,sy,r) 
char *fn ;
float bf[X*Y] ;
int sx, sy, r ;

{ int fd, i, j ;
  unsigned char c ;

  if ((fd = open(fn,O_RDONLY)) > 0) { ;
    for (i = 0 ; i < sy ; i++) {
      for (j = 0 ; j < sx ; j++) {
        if (read(fd,&c,sizeof(c)) != 1) {
          error(4) ;
        }
        bf[r*i + j] = (float)c ;
      }
    }
  }