uras.c

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

              }
              else {
                if (strcmp("-B",argv[i]) == 0) { 
                  if (i + 2 < argc) {
                    sscanf(argv[i+1],"%d",col) ;
                    sscanf(argv[i+2],"%d",row) ;
                    *binary = TRUE ;
                    i += 3 ;
                  }
                  else {
                    error(1) ;
                  }
                }
                else {
                  error(2) ;
                }  
              }
            }
          }
        }
      }
    }
  }
  if ((*sf) && (*re == TR2)) {
    error(12) ;
  }
  if ((*nbin_1 > N_BINS) || (*nbin_2 > N_BINS)) {
    error(15) ;
  }
  if ((*re == TR2) && (*histo)) {
    error(13) ;
  }
  if ((*re != TR1) && (*re != TR2) && (*re != NO_REG)) {
    error(2) ;
  }
  concat(argv[1],"/",in_path) ;
  concat(argv[2],"/",out_path) ;
  concat(in_path,argv[3],i_fname) ;
  strcpy(ext,"-sg") ;
  concat(ext,str_sg,ext) ;
  concat(ext,"-tg",ext) ;
  concat(ext,str_tg,ext) ;
  concat(ext,"-m",ext) ;
  itoa(*n_frame,t0) ;
  concat(ext,t0,ext) ;
  concat(ext,"-r",ext) ;
  if (*re == TR1) {
    concat(ext,"1",ext) ;
  }
  else {
    concat(ext,"2-",ext) ;
    concat(ext,str_r2,ext) ;
  }
  if (*sf) {
    concat(ext,"-f",ext) ;
    concat(ext,str_trsh,ext) ;
  }
  concat(out_path,argv[0],v_fname) ;
  concat(v_fname,".",v_fname) ;
  concat(v_fname,argv[3],v_fname) ;
  concat(v_fname,"F",v_fname) ;
  concat(v_fname,ext,v_fname) ;
  if (*histo) {
    concat(out_path,argv[0],h_fname) ;
    concat(h_fname,".",h_fname) ;
    concat(h_fname,argv[3],h_fname) ;
    concat(h_fname,"H",h_fname) ;
    concat(h_fname,ext,h_fname) ;
  }
  printf(" Generated velocity filename  : %s\n",v_fname) ;
  printf(" Generated histogram filename : %s\n",h_fname) ;
  printf(" Input stem name              : %s\n",i_fname) ;
  printf(" Correct velocity filename    : %s\n",c_fname) ;
  fflush(stdout) ;
}

/* 
           NAME : prod_histo(f,c_fn,h_fn,nbin_1,incr_1,nbin_2,incr_2,
                             ttl_err,ttl_std,dens)
   PARAMETER(S) :    f : flow field with parameters;
                  c_fn : correct velocities file name;
                  h_fn : histogram file name;
                nbin_1 : number of bins in histo 1;
                incr_1 : increment per bin in histo 1;
                nbin_2 : number of bins in histo 2;
                incr_2 : increment per bin in histo 2;
               ttl_err : error in flow field;
               ttl_std : standard deviation;
                  dens : density of flow field.
 
        PURPOSE : produces an error histogram
                  h1: error vs determinant;
                  h2: error vs condition number;
                  h3: h1 cumulated;
                  h4: h2 cumulated.

         AUTHOR : Steven Beauchemin
             AT : University of Western Ontario
           DATE : February 20 1992
*/

prod_histo(f,c_fn,h_fn,nbin_1,incr_1,nbin_2,incr_2,ttl_err,ttl_std,dens)
qnode_ptr_t f ;
string c_fn, h_fn ;
float incr_1, incr_2, *ttl_err, *ttl_std, *dens ;
int nbin_1, nbin_2 ;

{ qnode_ptr_t create_node(), q ;
  disp_vect_t u ;
  histo_t histo1[N_BINS], histo2[N_BINS], c_histo1[N_BINS], c_histo2[N_BINS] ;
  float max_det, max_cond, psi_error(), actual_x, actual_y, size_x, size_y, 
        offset_x, offset_y, det, cond, density, err, avg_err, t, x, y ;
  int fdf, nbytes, i, j, k, index1, index2, ttl_freq, abs_freq ;
  FILE *fdp ;
  extern int KERNEL_X, KERNEL_Y ;

  q = create_node(0,f->res,f->sizx,f->sizy,f->sizz,0) ;
  q->flow_ptr = (disp_field512_t *)malloc(sizeof(disp_field512_t)) ;
  if ((fdf = open(c_fn,O_RDONLY)) <= 0) {
    error(6) ;
  }

  nbytes = 0 ;
  nbytes += read(fdf,&actual_x,sizeof(float)) ;
  nbytes += read(fdf,&actual_y,sizeof(float)) ;
  nbytes += read(fdf,&size_x,sizeof(float)) ;
  nbytes += read(fdf,&size_y,sizeof(float)) ;
  nbytes += read(fdf,&offset_x,sizeof(float)) ;
  nbytes += read(fdf,&offset_y,sizeof(float)) ;
  
  for (i = (int)offset_y ; i < (int)actual_y ; i++) {
    for (j = (int)offset_x ; j < (int)actual_x ; j++) {
      nbytes += read(fdf,&y,sizeof(float)) ;
      nbytes += read(fdf,&x,sizeof(float)) ;
      (*q->flow_ptr)[q->res*i + j].y = y ;
      (*q->flow_ptr)[q->res*i + j].x = -x ;
    }
  }
  close(fdf) ;

  max_det = incr_1*(float)nbin_1 ;
  max_cond = incr_2*(float)nbin_2 ;
  for (i = 0 ; i < N_BINS ; i++) {
    histo1[i].avg = 0.0 ;
    histo1[i].std = 0.0 ;
    histo1[i].freq = 0 ;
    histo2[i].avg = 0.0 ;
    histo2[i].std = 0.0 ;
    histo2[i].freq = 0 ;
    c_histo1[i].avg = 0.0 ;
    c_histo1[i].std = 0.0 ;
    c_histo1[i].freq = 0 ;
    c_histo2[i].avg = 0.0 ;
    c_histo2[i].std = 0.0 ;
    c_histo2[i].freq = 0 ;
    abs_freq = 0 ;
    ttl_freq = 0 ;
    avg_err = 0.0 ;
    *ttl_err = 0.0 ;
    *ttl_std = 0.0 ;
  }
  for (i = KERNEL_Y + NRADIUS ; i < f->sizy - KERNEL_Y - NRADIUS ; i++) {
    for (j = KERNEL_X + NRADIUS ; j < f->sizx -  KERNEL_X - NRADIUS ; j++) {
      det = (*f->param_ptr[f->sizz/2])[f->res*i + j].gauss ;
      cond = (*f->param_ptr[f->sizz/2])[f->res*i + j].cond ;
      index1 = (int)((det/max_det)*(float)nbin_1) ;
      index2 = (int)((cond/max_cond)*(float)nbin_2) ;
      u = (*f->flow_ptr)[f->res*i+j] ;
      abs_freq++ ;
      if (u.x != -100.0 && u.y != 100.0) {
        err = psi_error((*f->flow_ptr)[f->res*i+j],(*q->flow_ptr)[q->res*i+j]) ;
        avg_err += err ;
        ttl_freq++ ;
        if (index1 >= nbin_1) {
          index1 = nbin_1 - 1 ;
        }
        histo1[index1].avg += err ;
        histo1[index1].freq++ ;
        if (index2 >= nbin_2) {
          index2 = nbin_2 - 1 ;
        }
        histo2[index2].avg += err ;
        histo2[index2].freq++ ;
        for (k = 0 ; k <= index1 ; k++) {
          c_histo1[k].avg += err ;
          c_histo1[k].freq++ ;
        }
        for (k = index2 ; k < nbin_2 ; k++) {
          c_histo2[k].avg += err ;
          c_histo2[k].freq++ ;
        }
      } 
    }
  }
  if (abs_freq != 0) {
    *dens = (float)ttl_freq/(float)abs_freq ;
  }
  else {
    *dens = 0.0 ;
  }
  *ttl_err = avg_err/(float)ttl_freq ;
  avg_err = *ttl_err ;
  for (i = 0 ; i < nbin_1 ; i++) {
    if (histo1[i].freq != 0) {
      histo1[i].avg /= (float)histo1[i].freq ;
    }
    if (c_histo1[i].freq != 0) {
      c_histo1[i].avg /= (float)c_histo1[i].freq ;
    } 
  }
  for (i = 0 ; i < nbin_2 ; i++) {
    if (histo2[i].freq != 0) {
      histo2[i].avg /= (float)histo2[i].freq ;
    }
    if (c_histo2[i].freq != 0) {
      c_histo2[i].avg /= (float)c_histo2[i].freq ;
    }
  }
  for (i = KERNEL_Y + NRADIUS ; i < f->sizy - KERNEL_Y - NRADIUS ; i++) {
    for (j = KERNEL_X + NRADIUS ; j < f->sizx - KERNEL_X - NRADIUS ; j++) {
      det = (*f->param_ptr[f->sizz/2])[f->res*i + j].gauss ;
      cond = (*f->param_ptr[f->sizz/2])[f->res*i + j].cond ;
      index1 = (int)((det/max_det)*(float)nbin_1) ;
      index2 = (int)((cond/max_cond)*(float)nbin_2) ;
      u = (*f->flow_ptr)[f->res*i+j] ;
      if (u.x != -100.0 && u.y != 100.0) {
        err = psi_error((*f->flow_ptr)[f->res*i+j],(*q->flow_ptr)[q->res*i+j]) ;
        *ttl_std += pow(avg_err - err,2.0) ;
        if (index1 >= nbin_1) {
          index1 = nbin_1 - 1 ;
        }
        histo1[index1].std += pow(histo1[index1].avg - err,2.0) ;
        if (index2 >= nbin_2) {
          index2 = nbin_2 - 1 ;
        }
        histo2[index2].std += pow(histo2[index2].avg - err,2.0) ;
        for (k = 0 ; k <= index1 ; k++) {
          c_histo1[k].std += pow(c_histo1[k].avg - err,2.0) ;
        }
        for (k = index2 ; k < nbin_2 ; k++) {
          c_histo2[k].std += pow(c_histo2[k].avg - err,2.0) ;
        }
      }
    }
  }
  *ttl_std = sqrt(*ttl_std/(float)ttl_freq) ;
  for (i = 0 ; i < nbin_1 ; i++) {
    if (histo1[i].freq != 0) {
      histo1[i].std = sqrt(histo1[i].std/(float)histo1[i].freq) ;
    }
    if (c_histo1[i].freq != 0) {
      c_histo1[i].std = sqrt(c_histo1[i].std/(float)c_histo1[i].freq) ;
    } 
  }
  for (i = 0 ; i < nbin_2 ; i++) {
    if (histo2[i].freq != 0) { 
      histo2[i].std = sqrt(histo2[i].std/(float)histo2[i].freq) ;
    }
    if (c_histo2[i].freq != 0) { 
      c_histo2[i].std = sqrt(c_histo2[i].std/(float)c_histo2[i].freq) ;
    }
  }
  if ((fdp = fopen(h_fn,"w")) == NULL) {
    error(6) ;
  }

  density = c_histo1[0].freq ;
  fprintf(fdp,"%2d\n",N_HISTO) ;
  fprintf(fdp,"%3d\n",nbin_1) ;
  fprintf(fdp,"%5.7f\n",max_det - incr_1/2.0) ;
  for (i = 0 ; i < nbin_1 ; i++) {
    t = max_det/(float)nbin_1 ;
    fprintf(fdp,"%5.7f %10.7f %10.7f %3.7f\n", t*(float)(i)+t/2.0, 
    histo1[i].avg, histo1[i].std, (float)histo1[i].freq/density) ;
  }
  fprintf(fdp,"\n\n\n") ;

  fprintf(fdp,"%3d\n",nbin_1) ;
  fprintf(fdp,"%5.7f\n",max_det - incr_1/2.0) ;
  for (i = 0 ; i < nbin_1 ; i++) {
    t = max_det/(float)nbin_1 ; 
    fprintf(fdp,"%5.7f %10.7f %10.7f %3.7f\n", t*(float)(i)+t/2.0, 
    c_histo1[i].avg, c_histo1[i].std, (float)c_histo1[i].freq/density) ;
  }
  fprintf(fdp,"\n\n\n") ;

  fprintf(fdp,"%3d\n",nbin_2) ;
  fprintf(fdp,"%5.7f\n",max_cond - incr_2/2.0) ;
  for (i = 0 ; i < nbin_2 ; i++) {
    t = max_cond/(float)nbin_2 ;
    fprintf(fdp,"%5.7f %10.7f %10.7f %3.7f\n", t*(float)(i)+t/2.0, 
    histo2[i].avg, histo2[i].std, (float)histo2[i].freq/density) ;
  }
  fprintf(fdp,"\n\n\n") ;

  fprintf(fdp,"%3d\n",nbin_2) ;
  fprintf(fdp,"%5.7f\n",max_cond - incr_2/2.0) ;
  for (i = 0 ; i < nbin_2 ; i++) {
    t = max_cond/(float)nbin_2 ; 
    fprintf(fdp,"%5.7f %10.7f %10.7f %3.7f\n", t*(float)(i)+t/2.0, 
    c_histo2[i].avg, c_histo2[i].std, (float)c_histo2[i].freq/density) ;
  }
  fclose(fdp) ;
}

/* 
           NAME : load_frames(fn,start,n_frame,x,y,cub,bin)
   PARAMETER(S) : fn : filename if the image sequence;
               start : start number of frames;
             n_frame : number of frames needed;
                 x,y : image size; 
                 cub : pointer on an image cube;
                 bin : binary input files.
 
        PURPOSE : Loads the image sequence from frame start to
                  frame start + n_frame into image cube cub1.

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

load_frames(fn,start,n_frame,x,y,cub,bin)
string fn ;
int start, n_frame, x, y, bin ;
qnode_ptr_t cub ;

{ unsigned char hd[H] ;
  string s1, s2 ;
  int i ;

  for (i = 0 ; i < n_frame ; i++) {
    itoa(start,s1) ;
    concat(fn,s1,s2) ;
    if (!bin) {
      pgetrast(s2,hd,cub->gauss_ptr[i],x,y,X) ;
    }
    else {
      Bpgetrast(s2,cub->gauss_ptr[i],x,y,X) ;
    }
    start++ ;
  }
}

/* 
           NAME : compute_deriv(c,f)
   PARAMETER(S) :     c : pointer on an image cube node;
                      f : pointer on a flow field and
                          its parameters.
 
        PURPOSE : Computes image derivatives needed for
                  the second-order method of Uras et al [88].

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

compute_deriv(c,f)
qnode_ptr_t c, f ;

{ float dxx(), dyy(), dxy(), dx(), dy(), dt(), dxt(), dyt() ;
  int i, j, k ;
param_t P ;
  for (i = f->ofst ; i < f->sizy - f->ofst ; i++) {
    for (j = f->ofst ; j < f->sizx - f->ofst ; j++) {
      for (k = 0 ; k < f->sizz ; k++) {
        (*f->param_ptr[k])[f->res*i + j].fx = dx(c,f,i,j,k) ;
        (*f->param_ptr[k])[f->res*i + j].fy = dy(c,f,i,j,k) ;
      }
    }
  }
  k = c->sizz/2.0 ; 
  for (i = f->ofst ; i < f->sizy - f->ofst ; i++) {
    for (j = f->ofst ; j < f->sizx - f->ofst ; j++) { 
      if ((*f->param_ptr[k])[f->res*i + j].err == NO_ERROR) { 
        (*f->param_ptr[k])[f->res*i + j].ft = dt(c,f,i,j) ;
        (*f->param_ptr[k])[f->res*i + j].fxt = dxt(f,i,j) ;
        (*f->param_ptr[k])[f->res*i + j].fyt = dyt(f,i,j) ;
      }
    }
  } 
  for (i = f->ofst ; i < f->sizy - f->ofst ; i++) {
    for (j = f->ofst ; j < f->sizx - f->ofst ; j++) { 
      if ((*f->param_ptr[k])[f->res*i + j].err == NO_ERROR) { 
        (*f->param_ptr[k])[f->res*i + j].fxx = dxx(f,i,j,k) ;
        (*f->param_ptr[k])[f->res*i + j].fyy = dyy(f,i,j,k) ;
        (*f->param_ptr[k])[f->res*i + j].fxy = dxy(f,i,j,k) ;
      }
    }
  }
}

/* 
           NAME : compute_discr(fl)
   PARAMETER(S) : flow field node.
 
        PURPOSE : computes ||Mt*gradient(U)||/||gradient(It)||
                  and gaussian curvature for all non-error points.

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

compute_discr(fl)
qnode_ptr_t fl ;

{ param_t P ;
  float dux(), duy(), dvx(), dvy(), ux, uy, vx, vy, MtU, It, cmax, cmin ;
  int i, j ;

  for (i = fl->ofst ; i < fl->sizy - fl->ofst ; i++) {
    for (j = fl->ofst ; j < fl->sizx - fl->ofst ; j++) {
      P = (*fl->param_ptr[fl->sizz/2])[fl->res*i + j] ;
      if (P.err == NO_ERROR) {
        ux = dux(fl,i,j) ;
        uy = duy(fl,i,j) ;
        vx = dvx(fl,i,j) ;
        vy = dvy(fl,i,j) ;
        MtU = pow(P.fx*ux + P.fy*vx,2.0) + pow(P.fx*uy + P.fy*vy,2.0) ;
        It = pow(P.fxt,2.0) + pow(P.fyt,2.0) ;
        P.discr = sqrt(MtU/It) ;
        P.gauss = fabs(P.fxx*P.fyy - pow(P.fxy,2.0)) ;
        cmax = 0.5*((P.fxx+P.fyy)+sqrt(pow(P.fxx-P.fyy,2.0)+4.0*P.fxy*P.fxy)) ;
        cmin = 0.5*((P.fxx+P.fyy)-sqrt(pow(P.fxx-P.fyy,2.0)+4.0*P.fxy*P.fxy)) ;
        if (fabs(cmin) < fabs(cmax)) {
          if (cmin != 0.0) {
            P.cond = fabs(cmax)/fabs(cmin) ;
          }
          else {
            P.cond = MAX_COND ;
          }
        }
        else {
          if (cmax != 0.0) {
            P.cond = fabs(cmin)/fabs(cmax) ;
          }
          else {
            P.cond = MAX_COND ;
          }
        }
        (*fl->param_ptr[fl->sizz/2])[fl->res*i + j] = P ;
      }
    }
  }
}

/* 
           NAME : compute_flow(fl)
   PARAMETER(S) : fl : pointer on a flow field and its parameters.
 
        PURPOSE : Computes image displacements with
                  the second-order method of Nagel[87].

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

compute_flow(fl)
qnode_ptr_t fl ;

{ disp_vect_t u, v ;
  param_t P ;
  float mag() ;
  int i, j ;

  for (i = fl->ofst ; i < fl->sizy - fl->ofst ; i++) {
    for (j = fl->ofst ; j < fl->sizx - fl->ofst ; j++) { 
      P = (*fl->param_ptr[fl->sizz/2])[fl->res*i + j] ;
      if (P.err == NO_ERROR) {      
        if (P.fxx*P.fyy - pow(P.fxy,2.0) != 0.0) {
          u.x = (P.fyt*P.fxy - P.fxt*P.fyy)/(P.fxx*P.fyy - pow(P.fxy,2.0)) ;
          u.y = (P.fxt*P.fxy - P.fyt*P.fxx)/(P.fxx*P.fyy - pow(P.fxy,2.0)) ;
          if (mag(u) > MAXFLOW) { 
            v.x = (u.x/mag(u))*MAXFLOW ;
            v.y = (u.y/mag(u))*MAXFLOW ;
            u = v ;
          }
        }
        else {  
          u.x = 0.0 ;
          u.y = 0.0 ;
          (*fl->param_ptr[fl->sizz/2])[fl->res*i + j].err = NO_FLOW ;
        }
      }
      else {
        u.x = 0.0 ;
        u.y = 0.0 ; 
      }
      (*fl->flow_ptr)[fl->res*i + j] = u ;
    }
  }
}

/* 
           NAME : uras.c
   PARAMETER(S) : line parameters.

        PURPOSE : Computes optic displacement
                  between two image frames.

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

main(argc,argv)
int argc ;
char *argv[] ;

{ qnode_ptr_t create_node(), cub1h, cub1q, cub2h, cub2q, flh