#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 80typedef 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] ;      }