uras.c

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

    }
  }
  return(d/Ix.f) ;
}

/* 
           NAME : dvy(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 v of flow p.

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

float dvy(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[j+h][i+h] ;
      }
    }
  }
  return(d/Ix.f) ;
}

/* 
           NAME : dx(p,q,x,y,z)
   PARAMETER(S) : p : image cube node;
                  q : parameters of flow (image derivatives);
              x,y,z : image location for computation.
 
        PURPOSE : computes 1st order partial derviative with respect
                  to x at x,y for image p->gauss_ptr[z].

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

float dx(p,q,x,y,z)
qnode_ptr_t p, q ;
int x, y, z ;

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

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

  if (overflow(x,y,p->sizy,p->sizx,p->ofst,h,1)) {
    (*q->param_ptr[z])[q->res*x + y].err = SA_OVERFLOW ;
  }
  else {
    (*q->param_ptr[z])[p->res*x + y].err = NO_ERROR ;
    for (i = -h ; i <= h ; i++) {
      d += (*p->gauss_ptr[z])[p->res*(x+i) + y]*C.k[i+h] ;
    }
    d /= C.f ;
  }
  return(d) ;
}

/* 
           NAME : dxt(p,x,y)
   PARAMETER(S) : p : image parameter node;
                x,y : image location for computation. 

        PURPOSE : computes partial derivative with respect to
                  x,t at x,y.

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

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

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

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

  if  ((*p->param_ptr[0])[p->res*x + y].err == NO_ERROR) {
    for (i = -h ; i <= h ; i++) {
      d += (*p->param_ptr[i+h])[p->res*x + y].fx*C.k[i+h] ;
    }
    d /= C.f ;
  }
  return(d) ;
}

/* 
           NAME : dxx(q,x,y,z)
   PARAMETER(S) : q : parameters of flow (image derivatives);
              x,y,z : image location for computation.
 
        PURPOSE : computes 2nd order partial derviative with respect
                  to x at x,y. 

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

float dxx(q,x,y,z)
qnode_ptr_t q ;
int x, y, z ;

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

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

  if (overflow(x,y,q->sizy,q->sizx,q->ofst,h,2)) {
    (*q->param_ptr[z])[q->res*x + y].err = SA_OVERFLOW ;
  }
  else {
    (*q->param_ptr[z])[q->res*x + y].err = NO_ERROR ;
    for (i = -h ; i <= h ; i++) {
      d += (*q->param_ptr[z])[q->res*(x+i) + y].fx*C.k[i+h] ;
    }
    d /= C.f ;
  }
  return(d) ;
}

/* 
           NAME : dxy(q,x,y,z)
   PARAMETER(S) : q : parameters of flow (image derivatives);
              x,y,z : image location for computation.
 
        PURPOSE : computes 2nd order partial derviative with respect
                  to x,y at x,y. 

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

float dxy(q,x,y,z)
qnode_ptr_t q ;
int x, y, z ;

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

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

  if (overflow(x,y,q->sizy,q->sizx,q->ofst,h,2)) {
    (*q->param_ptr[z])[q->res*x + y].err = SA_OVERFLOW ;
  }
  else {
    (*q->param_ptr[z])[q->res*x + y].err = NO_ERROR ;
    for (i = -h ; i <= h ; i++) {
      d += (*q->param_ptr[z])[q->res*(x+i) + y].fy*C.k[i+h] ;
    }
    d /= C.f ;
  }
  return(d) ;
}

/* 
           NAME : dy(p,q,x,y,z)
   PARAMETER(S) : p : image cube node;
                  q : parameters of flow (image derivatives);
              x,y,z : image location for computation.
 
        PURPOSE : computes 1st order partial derviative with respect
                  to y at x,y for image p->gauss_ptr[z].

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

float dy(p,q,x,y,z)
qnode_ptr_t p, q ;
int x, y, z ;

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

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

  if (overflow(x,y,p->sizy,p->sizx,p->ofst,h,1)) {
    (*q->param_ptr[z])[q->res*x + y].err = SA_OVERFLOW ;
  }
  else {
    (*q->param_ptr[z])[p->res*x + y].err = NO_ERROR ;
    for (i = -h ; i <= h ; i++) {
      d += (*p->gauss_ptr[z])[p->res*(x) + y + i]*C.k[i+h] ;
    }
    d /= C.f ;
  }
  return(d) ;
}

/* 
           NAME : dyt(p,x,y)
   PARAMETER(S) : p : image parameter node;
                x,y : image location for computation. 

        PURPOSE : computes partial derivative with respect to
                  y,t at x,y.

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

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

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

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

  if  ((*p->param_ptr[0])[p->res*x + y].err == NO_ERROR) {
    for (i = -h ; i <= h ; i++) {
      d += (*p->param_ptr[i+h])[p->res*x + y].fy*C.k[i+h] ;
    }
    d /= C.f ;
  }
  return(d) ;
}

/* 
           NAME : dyy(q,x,y,z)
   PARAMETER(S) : q : parameters of flow (image derivatives);
              x,y,z : image location for computation.
 
        PURPOSE : computes 2nd order partial derviative with respect
                  to y at x,y. 

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

float dyy(q,x,y,z)
qnode_ptr_t q ;
int x, y, z ;

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

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

  if (overflow(x,y,q->sizy,q->sizx,q->ofst,h,2)) {
    (*q->param_ptr[z])[q->res*x + y].err = SA_OVERFLOW ;
  }
  else {
    (*q->param_ptr[z])[q->res*x + y].err = NO_ERROR ;
    for (i = -h ; i <= h ; i++) {
      d += (*q->param_ptr[z])[q->res*x + y+i].fy*C.k[i+h] ;
    }
    d /= C.f ;
  }
  return(d) ;
}

/* 
           NAME : usage()
   PARAMETER(S) : None.
 
        PURPOSE : Prints the appropriate usage and 
                  stops the run.

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

usage()

{ 
  fprintf(stderr,"Usage: uras input-path output-path seq-name. [-SG n.n] [-TG n.n] [-M n] [-R 1 | -R 2 n.n] [-B rows cols] [-F n.n] [-C corr-vel-file nbins1 increment1 nbins2 increment2]\n") ;
  fprintf(stderr,"[-SG n.n]         : sigma value for spatial gaussian\n") ;
  fprintf(stderr,"                    default value is 3.0\n") ;
  fprintf(stderr,"[-TG n.n]         : sigma value for temporal gaussian\n") ; 
  fprintf(stderr,"                    default value is 1.5\n") ;
  fprintf(stderr,"[-M n]            : middle frame of image sequence\n") ;
  fprintf(stderr,"                    default value is 0\n") ;
  fprintf(stderr,"[-R 1 | -R 2 n.n] : n = 1 : applies TR1 regularization\n") ; 
  fprintf(stderr,"                    n = 2 : applies TR2 regularization\n") ;
  fprintf(stderr,"                    with sigma = n.n\n") ;
  fprintf(stderr,"[-B cols rows]    : for binary input files without header\n") ;
  fprintf(stderr,"[-F n.n]          : filters out unreliable estimates\n") ;
  fprintf(stderr,"                    using the Gaussian curv. (det(H))\n") ;
  fprintf(stderr,"[-C corr-vel-file nbins1 increment1 nbins2 increment2]\n") ;
  fprintf(stderr,"                  : error histogram\n") ;  
  fprintf(stderr,"                    corr-vel-file : file of correct velocities\n") ;
  fprintf(stderr,"                    nbins1, nbins2 : number of bins for histograms\n") ;
  fprintf(stderr,"                    increment1, increment2 : value of gap between bins\n") ;
  exit(-1) ;
}

/* 
           NAME : error(n)
   PARAMETER(S) : n : error number.
 
        PURPOSE : Prints the appropriate error message and 
                  stops the run.

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

error(n)
int n ;

{ 
  switch(n) {
    
     case 1 : fprintf(stderr,"error %d : wrong number of arguments\n",n) ;
              break ;
     case 2 : fprintf(stderr,"error %d : invalid option\n",n) ;
              break ;
     case 3 : fprintf(stderr,"error %d : maximum number of frames is %d\n",n,Z) ;
              break ;
     case 4 : fprintf(stderr,"error %d : read error\n",n) ;
              exit(-1) ;
              break ;
     case 5 : fprintf(stderr,"error %d : header read error\n",n) ;
              exit(-1) ;
              break ;
     case 6 : fprintf(stderr,"error %d : file open error\n",n) ;
              exit(-1) ;
              break ;
     case 7 : fprintf(stderr,"error %d : write error\n",n) ;
              exit(-1) ;
              break ;
     case 8 : fprintf(stderr,"error %d : header write error\n",n) ;
              exit(-1) ;
              break ;
     case 9 : fprintf(stderr,"error %d : file create error\n",n) ;
              exit(-1) ;
              break ;
     case 10: fprintf(stderr,"error %d : maximum image size is %d by %d\n",n,X,Y) ;
              break ;
     case 11: fprintf(stderr,"error %d : option -R n mandatory\n",n) ;
              break ;
     case 12: fprintf(stderr,"error %d : no filtering with -R 2\n",n) ;
              break ;
     case 13: fprintf(stderr,"error %d : no histogram with -R 2\n",n) ;
              break ;
     case 14: fprintf(stderr,"error %d : undefined error\n",n) ; 
              break ;
     case 15: fprintf(stderr,"error %d : maximum number of bins is %d\n",
              n,N_BINS) ;
              break ;
    default : fprintf(stderr,"error %d : undefined error\n",n) ;
              exit(-1) ;
              break ;
  }
  usage() ;
}

/* 
           NAME : frames(s)
   PARAMETER(S) : s : Length of temporal Gaussian.
 
        PURPOSE : returns the number of frames necessary for
                  the smoothing and the computation of
                  4-point central differences.

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

int frames(s)
int s ;

{ 
  if (s - 1 + SMSK > Z) {
    error(3) ;
  }
  else {
   return(s - 1  + SMSK) ;
  }
}

/* 
           NAME : start_number(middle,s)
   PARAMETER(S) : middle : number of middle frame;
                       s : length of temporal mask for derivatives;
 
        PURPOSE : computes the number of the start frame given
                  the middle frame and the length of the temporal
                  masks for derivatives. 
                 

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

int start_number(middle,s)
int middle, s ;

{ return(middle - (SMSK + s - 1)/2) ; }

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

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

free_cube(p,n)
qnode_ptr_t p ;

{ int i ;

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

/* 
           NAME : free_flow(p,n) ;
   PARAMETER(S) : p : pointer on an flow and parameter node.
 
        PURPOSE : deallocates a flow field to p->flow_ptr
                  and n arrays of flow parameters to 
                  p->param_ptr[i].

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

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

{ int i ;

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

/* 
           NAME : generate_gauss(ker,sig)
   PARAMETER(S) : ker : 1D  kernel;
                  sig : sigma.
 
        PURPOSE : Inits the kernel values with G(x,sig).

         AUTHOR : Steven Beauchemin
             AT : University of Western Ontario
           DATE : April 6 1990

*/

generate_gauss(ker,sig)
kernel_t *ker ;
float sig ;

{ int i, h ;

  h = (int)(sig*6.0 + 1.0) ;
  if (h%2 == 0) {
    h++ ;
  }
  (*ker).f = 0.0 ;
  if (sig != 0.0) {
    for (i = -(int)(h/2.0) ; i <= (int)(h/2.0) ; i++) {
      (*ker).k[i+(int)(h/2.0)] = exp(-pow((float)i,2.0)/(2.0*pow(sig,2.0)))/
      (sqrt(2.0*PI)*sig) ;
      (*ker).f += (*ker).k[i+(int)(h/2.0)] ;
    }
  }
  else {