uras.c

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

  else {
    error(6) ;
  }
  close(fd) ;
}

/* 
           NAME : pputrast(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 write;
                  r     : resolution of the array.
 
        PURPOSE : Writes the contents of bf into a 
                  rasterfile specified by fn.

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

*/

pputrast(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 = creat(fn,0600)) > 0) { ;
    if (write(fd,hd,H) == H) { 
      for (i = 0 ; i < sy ; i++) { 
        for (j = 0 ; j < sx ; j++) { 
          c = (unsigned char)bf[r*i + j] ;
          if (write(fd,&c,sizeof(c)) != 1) {
            error(7) ;
          }
        }
      }
    }
    else {
      error(8) ; 
    }
  }
  else {
    error(9) ;
  }
  close(fd) ;
}

/* 
           NAME : propagate(u,det,cn,q,i,j,n)
   PARAMETER(S) :   u : flow vector to be propagated;
                  det : determinant;
                   cn : condition number;
                    q : flow field for propagation;
                  i,j : coordinates of upper left corner of image area;
                    n : image area size.
 
        PURPOSE : propagates u, det, cn in a n*n square image area with upper 
                  left corner coordinates (i,j).

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

propagate(u,det,cn,q,i,j,n)
disp_vect_t u ;
float det, cn ;
qnode_ptr_t q ;
int i, j, n ;

{ int k, l ;

  for (k = 0 ; k < n ; k++) {
    for (l = 0 ; l < n ; l++) {
      (*q->flow_ptr)[q->res*(i+k) + j+l] = u ;
      (*q->param_ptr[q->sizz/2])[q->res*(i+k) + j+l].gauss = det ;
      (*q->param_ptr[q->sizz/2])[q->res*(i+k) + j+l].cond = cn ;
    }
  }
}

/* 
           NAME : psi_error(ve,va)
   PARAMETER(S) : ve : estimated flow vector;
                  va : accurate flow vector.
 
        PURPOSE : computes angular error of ve with respect to va
                  using Fleet [90] angular error metric

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

float psi_error(ve,va)
disp_vect_t ve, va ;

{ float norm(), v ;
  float VE[3], VA[3] ;

  VE[0] = ve.x ;
  VE[1] = ve.y ;
  VE[2] = 1.0 ;

  VA[0] = va.x ;
  VA[1] = va.y ;
  VA[2] = 1.0 ;

  v = (VE[0]*VA[0] + VE[1]*VA[1] + 1.0)/(norm(VA,3)*norm(VE,3)) ;
  if ((v > 1.0) && (v < 1.0001)) {
    v = 1.0 ;
  }
  return((float)(acos(v))*180.0/PI) ;
}

/* 
           NAME : raster_size(fn,num,x,y)
   PARAMETER(S) : fn : raster file name; 
                 num : frame number;
                   x : raster width;
                   y : raster height.
 
        PURPOSE : Reads the header of raterfile fn to get the
                  size. 

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

raster_size(fn,num,x,y)
char *fn ;
int num, *x, *y ;

{ raster_t hd ;
  string s1, s2 ;
  int fd ;

  itoa(num,s1) ;
  concat(fn,s1,s2) ;
  if ((fd = open(s2,O_RDONLY)) > 0) {
    if (read(fd,&hd,sizeof(hd)) == sizeof(hd)) {
      *x = hd.width ;
      *y = hd.height ;
      if (*x > X || *y > Y) {
        error(10) ;
      }
    }
  }
  else {
    error(6) ;
  }
  close(fd) ;
}

/* 
           NAME : binary(row,col,x,y)
   PARAMETER(S) : row : number of rows in input image;
                  col : number of cols in input image;
                    x : raster width;
                    y : raster height.
 
        PURPOSE : sets the size of images when option -B used.

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

binary_size(row,col,x,y)
int row, col, *x, *y ;

{ *x = col ;
  *y = row ;
  if ((*x > X) || (*y > Y)) {
    error(10) ;
  }
}

/* 
           NAME : vector(loc,f)
   PARAMETER(S) : loc : flow field location;
                    f : node pointer.
 
        PURPOSE : returns the vector located at loc in f.

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

disp_vect_t vector(loc,f)
pos_t loc ;
qnode_ptr_t f ;

{ disp_vect_t u ;
  
  if ((loc.i != -1) && (loc.j != -1)) {
    u = (*f->flow_ptr)[f->res*loc.i + loc.j] ;
  }
  else {
    u.x = -100.0 ;
    u.y = 100.0 ;
  }
  return(u) ;
}

/* 
           NAME : regul1(f,n,sf,trsh)
   PARAMETER(S) : f : node pointer;
                  n : image area size;
                 sf : thresholding boolean value;
               trsh : threshold for Gaussian curvature.
 
        PURPOSE : operates procedure 1 (TR1 described in Uras et al.[88])
                  for regularization of flow field f.

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

regul1(f,n,sf,trsh)
qnode_ptr_t f ;
int n, sf ;
float trsh ;

{ pos_t min_cond(), sample[X], rep_loc ;
  disp_vect_t vector(), u ;
  float condition(), determinant(), cond_num, det ;
  int num, xf, yf, i, j, k ;

  xf = (int)(((f->sizx-f->ofst-NRADIUS-(f->ofst+NRADIUS))%n)/2.0 + 0.5) ;
  yf = (int)(((f->sizy-f->ofst-NRADIUS-(f->ofst+NRADIUS))%n)/2.0 + 0.5) ;
  KERNEL_X += xf ;
  KERNEL_Y += yf ;

  for (i = f->ofst+NRADIUS+yf ; i < f->sizy-f->ofst-NRADIUS-yf ; i += n) {
    for (j = f->ofst+NRADIUS+xf ; j < f->sizx-f->ofst-NRADIUS-xf ; j += n) {
      min_discr(f,i,j,n,sample,&num) ;
      rep_loc = min_cond(sample,f,n,num,trsh,sf) ;
      u = vector(rep_loc,f) ;
      cond_num = condition(rep_loc,f) ;
      det = determinant(rep_loc,f) ;
      propagate(u,det,cond_num,f,i,j,n) ;
    }
  }
}

/* 
           NAME : regul2(fl,ker)
   PARAMETER(S) :  fl : node pointer;
                  ker : 1D kernel for convolution of flow.
   
 
        PURPOSE : operates procedure 2 (TR2 described in Uras et al.[88])
                  on flow field fl.

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

regul2(fl,ker)
qnode_ptr_t fl ;
kernel_t ker ;

{ disp_field512_t *t ;
  disp_vect_t s ;
  int h, i, j, k ;

  h = ker.m/2 ;
  t = (disp_field512_t *)malloc(sizeof(disp_field512_t)) ;

  for (i = fl->ofst+NRADIUS ; i < fl->sizy-fl->ofst-NRADIUS ; i++) {
    for (j = fl->ofst+NRADIUS+h ; j < fl->sizx-fl->ofst-NRADIUS-h ; j++) {
      s.x = 0.0 ;
      s.y = 0.0 ;
      for (k = 0 ; k < ker.m ; k++) {
        s.x = s.x + (*fl->flow_ptr)[fl->res*i + j+k-h].x*ker.k[k] ;
        s.y = s.y + (*fl->flow_ptr)[fl->res*i + j+k-h].y*ker.k[k] ;
      }
      (*t)[fl->res*i + j].x = s.x/ker.f ;
      (*t)[fl->res*i + j].y = s.y/ker.f ;
    }
  }

     
  for (i = fl->ofst+NRADIUS+h ; i < fl->sizy-fl->ofst-NRADIUS-h ; i++) {
    for (j = fl->ofst+NRADIUS ; j < fl->sizx-fl->ofst-NRADIUS ; j++) {
      s.x = 0.0 ;
      s.y = 0.0 ;
      for (k = 0 ; k < ker.m ; k++) {
        s.x = s.x + (*t)[fl->res*(i+k-h) + j].x*ker.k[k] ;
        s.y = s.y + (*t)[fl->res*(i+k-h) + j].y*ker.k[k] ;
      }
      (*fl->flow_ptr)[fl->res*i + j].x = s.x/ker.f ;
      (*fl->flow_ptr)[fl->res*i + j].y = s.y/ker.f ;
    }
  }
  free((disp_field512_t *)t) ;
}

/*
           NAME : write_velocity(fn,flow_field)
   PARAMETER(S) : flow_field : the 2D array containing velocities 
                               or displacements.

        PURPOSE : Output field using Travis Burkitt's format.

         AUTHOR : Travis Burkitt, updated by Steven Beauchemin
             AT : University of Western Ontario
           DATE : May 7 1990
*/

write_velocity(fn,p)
qnode_ptr_t p ;
char *fn ;

{ extern int KERNEL_X, KERNEL_Y ;
  float x, y ; 
  int i, j, fdf, bytes ;

  if ((fdf=creat(fn,0600)) < 1) {
    error(6) ;
  }
  
  x = p->sizx ;
  y = p->sizy ;
  write(fdf,&x,4) ;
  write(fdf,&y,4) ;
  
  x = ((p->sizx-KERNEL_X-NRADIUS)-(KERNEL_X+NRADIUS)+SKIP-1)/SKIP ;
  y = ((p->sizy-KERNEL_Y-NRADIUS)-(KERNEL_Y+NRADIUS)+SKIP-1)/SKIP ;
  write(fdf,&x,4);
  write(fdf,&y,4);
  
  x = (KERNEL_X+NRADIUS+SKIP-1)/SKIP ;
  y = (KERNEL_Y+NRADIUS+SKIP-1)/SKIP ;
  write(fdf,&x,4) ;
  write(fdf,&y,4) ;
  bytes = 24 ;
  
  for(i = KERNEL_Y + NRADIUS ; i < p->sizy - KERNEL_Y - NRADIUS ; i++) {
    for(j = KERNEL_X + NRADIUS ; j < p->sizx - KERNEL_X - NRADIUS ; j++) {
      x = (*p->flow_ptr)[p->res*i + j].y ;
      y = -(*p->flow_ptr)[p->res*i + j].x ;
      write(fdf,&x,4) ;
      write(fdf,&y,4) ;
      bytes += 8 ;
    }
  }
  close(fdf) ;
}

/* 
           NAME : valid_option(argc,argv,in_path,out_path,sigma1,sigma2,
                  sigma3,histo,re,sf,trsh,i_fname,v_fname,c_fname,h_fname,
                  nbin_1,incr_1,nbin_2,incr_2,n_frame,binary,row,col) 

   PARAMETER(S) : argc : argument count;
                  argv : argument values;
               in_path : path name for input data;
              out_path : path name for output data;
                sigma1 : spatial sigma;
                sigma2 : temporal sigma;
                sigma3 : flow field sigma (TR2);
                 histo : error histogram option;
                    re : regularization option;
                    sf : filtering option;
                  trsh : threshold value for filtering;
               i_fname : input filename;
               v_fname : velocity filename;
               c_fname : correct velocities file name;
               h_fname : hisogram data file name ;
                nbin_1 : number of bins in histo 1;
                incr_1 : increment in histo 1 ;
                nbin_2 : number of bins in histo 2 ;
                incr_2 : increment in histo2;
               n_frame : number of required frames;
                binary : input files without header;
                   row : number of rows in input files;
                   col : number of cols in input files.
 
        PURPOSE : Validates line arguments.

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

valid_option(argc,argv,in_path,out_path,sigma1,sigma2,sigma3,histo,re,sf,trsh,
             i_fname,v_fname,c_fname,h_fname,nbin_1,incr_1,nbin_2,incr_2,
             n_frame,binary,row,col) 
int argc ;
char *argv[] ;
int *nbin_1, *nbin_2, *n_frame, *histo, *re, *sf, *binary, *row, *col ;
float *incr_1, *incr_2, *sigma1, *sigma2, *sigma3, *trsh ;
string in_path, out_path, i_fname, v_fname, c_fname, h_fname ;

{ int i ;
  string ext, t0, str_sg, str_tg, str_r2, str_trsh ;

  if (argc == 1) {
    usage() ;
  }
  if ((argc <= 23) && (argc >= 4)) {
    *binary = FALSE ;
    *re = TR1 ;
    *sf = FALSE ;
    *histo = FALSE ;
    *trsh = 0.0 ;
    *sigma1 = DEF_S1 ;
    strcpy(str_sg,"3.0") ;
    *sigma2 = DEF_S2 ;
    strcpy(str_tg,"1.5") ;
    *sigma3 = DEF_S3 ; 
    strcpy(str_r2,"3.0") ;
    *n_frame = 0 ;
    *nbin_1 = 0 ;
    *nbin_2 = 0 ;
    strcpy(i_fname,"Undefined") ;
    strcpy(v_fname,"Undefined") ;
    strcpy(c_fname,"Undefined") ;
    strcpy(h_fname,"Undefined") ;
  }
  else {
    error(1) ;
  }
  i = 4 ;
  while (i < argc) {
    if (strcmp("-R",argv[i]) == 0) {
      if (i + 1 < argc) {
        sscanf(argv[i+1],"%d",re) ;
        i += 2 ;
        if (*re == TR2) {
          if (i < argc) {
            sscanf(argv[i],"%f",sigma3) ;
            strcpy(str_r2,argv[i]) ;
            i += 1 ;
          }
          else {
            error(1) ;
          }
        }
      }
      else {
        error(1) ;
      }
    }
    else {
      if (strcmp("-F",argv[i]) == 0) {
        if (i + 1 < argc) {
          sscanf(argv[i+1],"%f",trsh) ;
          strcpy(str_trsh,argv[i+1]) ;
          *sf = TRUE ;
          i += 2 ;
        }
        else {
          error(1) ;
        }
      }
      else {
        if (strcmp("-C",argv[i]) == 0) {
          if (i + 5 < argc) {
            strcpy(c_fname,argv[i+1]) ;
            sscanf(argv[i+2],"%d",nbin_1) ;
            sscanf(argv[i+3],"%f",incr_1) ;
            sscanf(argv[i+4],"%d",nbin_2) ;
            sscanf(argv[i+5],"%f",incr_2) ;
            *histo = TRUE ;
            i += 6 ;
          }
          else {
            error(1) ;
          }
        } 
        else {
          if (strcmp("-SG", argv[i]) == 0) {
            if (i + 1 < argc) {
              sscanf(argv[i+1],"%f",sigma1) ;
              strcpy(str_sg,argv[i+1]) ;
              i += 2 ;
            }
            else {
              error(1) ;
            }
          }
          else {
            if (strcmp("-TG",argv[i]) == 0) {
              if (i + 1 < argc) {
                sscanf(argv[i+1],"%f",sigma2) ;
                strcpy(str_tg,argv[i+1]) ;
                i += 2 ;
              }
              else {
                error(1) ;
              } 
            }
            else {
              if (strcmp("-M",argv[i]) == 0) {
                if (i + 1 < argc) {
                  sscanf(argv[i+1],"%d",n_frame) ;
                  i += 2 ;
                }
                else {
                  error(1) ;
                }