/*---------------------------------------------------------------------------*/
/*  vectorfield.c                                                            */
/*---------------------------------------------------------------------------*/

/*\
 *  This file is part of the Omnivis project, Harvard University 
 *  Copyright (C) December 2009  Oliver Knill and Jose Ramirez Herran 
 *  http://www.math.harvard.edu/~knill/3dscan2           (project page)
 *
 *  Usage:  gcc -o vectorfield -lm vectorfield.c;   
 *          ./vectorfield `ls *.pnm`
 *
 *  In a first run, it will produce auxiliary files in a directory tmp
 *  like curvature, gradient and smoothed pictures
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, write to the Free Software
 *  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
\*/

#include <stdio.h>                          // reading in standard libraries
#include <stdlib.h>                         // this compiles in any unix 
#include <string.h>                         // environment with gcc
#include <math.h>

#define PI 3.14159265358979323846           // pi
#define FILE_INDEX(x) "",(x)                // to generate file names 

#define MP 10000                            // maximal number of points to track and print
#define NF 300                              // maximal number of input files
#define S  4                                // how many times is the picture smoothed first
#define PS 4                                // patch inflation size for curvature
#define PS1 2                               // patch inflation size for correspondence  
#define RR 5                                // length of slope line which is drawn if drawn at all
#define M 10                                // max search distance for best nearby pattern 
#define N 5                                 // size of pattern patch to compare  
#define K 3                                 // search grid dimension, look at each K'th pixel
#define KK 1                                // grid size to average derivatives over
#define bb 20                               // boundary width to search 
#define merge 3                             // radius within which we merge points
double  threshold =  6000;                  // gradient      New: is now be determined
#define ethreshold 100                      // threshold for edges    New: will now be determined
#define scaledthresh 170                    // new threshold scaled to 255,  70 for florida, 150 for cubes
#define stroke  100                         // how long to draw red stream lines for Hamiltonian
double alpha=PI/2+0*PI/2;                   // PI/2 for Hamiltonian ,  PI for gradient
int  preprocessed;                          // =1 if pictures are in tmp, otherwise 0, meaning we have to pre process
int  smootheddone;                          // =1 if smoothed pictures are in tmp, otherwise 0

int  x_match[MP][NF];                       // x-coordinate of a matched point, [point number,camera number]
int  y_match[MP][NF];                       // y-coordinate of a matched point, [point number,camera number]
int  countpoints[NF];                       // number of points which qualify

FILE *out[NF];                              // output file
FILE *smoothout[NF];                        // output file for smooth pictures
FILE *outcurvature[NF];                     // output file for curvature pictures
FILE *outgradient[NF];                      // output file for gradient pictures
FILE *in[NF];                               // input file
FILE *insmoothed[NF];                       // smoothed input file 
FILE *incurvature[NF];                      // curvature input file 
FILE *ingradient[NF];                       // gradient input file 

char *infilename[NF];                       // input file name
char *outfilename[NF];                      // output file name
char *smoothfilename[NF];                   // output file name for smooth pictures
char *curvaturefilename[NF];                // output file name for curvature pictures
char *gradientfilename[NF];                 // output file name for gradient pictures
char *outputarray[NF];                      // pointers for output
char *outcurvaturearray[NF];                // pointers for curvature output
char *outgradientarray[NF];                 // pointers for gradient
char *inputarray[NF];                       // pointer for input
char *insmoothedarray[NF];                  // pointers for smoothed input
char *incurvaturearray[NF];                 // pointers for curvature input
char *ingradientarray[NF];                  // pointers for gradient input

int  x_size,y_size,c_size;                  // width, height and number of colors
int  midframe;                              // equal to 1 for NF=3,  equal to 1 for NF=4 

long gaussian[5][5];                        // Gaussian kernel
long laplace[5][5];                         // Laplace operator derivative
long hessianxx[5][5];                       // Hessian xx
long hessianxy[5][5];                       // Hessian xy
long hessianyy[5][5];                       // Hessian yy
long sobelx[5][5];                          // Sobel partial x derivative
long sobely[5][5];                          // Sobel partial y derivative
long rpatch[5][5];                          // Patch for smoothing 
long gpatch[5][5];                          // Patch for smoothing 
long bpatch[5][5];                          // Patch for smoothing 
long gaussiansum;                           // sum of Gaussian kernel entries

char *ptr;                                  // general pointer for output picture

long    rx,ry,gx,gy,bx,by;                  // red green blue value of two pictures
long    rthresh,rthreshx,rthreshy;
long    gthresh,gthreshx,gthreshy;
long    bthresh,bthreshx,bthreshy;          // threshold functions
double  thresh,edge,cross;                  // average over colors
double  maxthresh, maxedge,maxcross;        // maximal values 
int     imaxthresh, imaxedge,imaxcross;     // maximal values 
int     maxrthresh[NF]; 
int     maxgthresh[NF];
int     maxbthresh[NF];                     // maximal gradient values
int     maxredge[NF];
int     maxgedge[NF];
int     maxbedge[NF];                       // maximal curvature values
int     redge,gedge,bedge;                  // edge values
long    rcross,gcross,bcross;               // cross values 
double  rang,gang,bang;                     // red green blue angle
double  rcos,gcos,bcos;                     // cos(angles) of gradients
double  rsin,gsin,bsin;                     // sin(angles) of gradients

/*---------------------------------------------------------------------------*/
/*  Define first and second derivatives of the red,green and blue functions  */
/*---------------------------------------------------------------------------*/

void definematrices() {
  gaussian[0][0]=1; gaussian[0][1]=3; gaussian[0][2]=4; gaussian[0][3]=3; gaussian[0][4]=1;  
  gaussian[1][0]=3; gaussian[1][1]=7; gaussian[1][2]=10;gaussian[1][3]=7; gaussian[1][4]=3;  
  gaussian[2][0]=4; gaussian[2][1]=10;gaussian[2][2]=14;gaussian[2][3]=10;gaussian[2][4]=4;  
  gaussian[3][0]=3; gaussian[3][1]=7; gaussian[3][2]=10;gaussian[3][3]=7; gaussian[3][4]=3;  
  gaussian[4][0]=1; gaussian[4][1]=3; gaussian[4][2]=4; gaussian[4][3]=3; gaussian[4][4]=1;
  gaussiansum=126; 

  sobely[0][0]=0; sobely[0][1]=0;  sobely[0][2]=0; sobely[0][3]=0; sobely[0][4]=0;
  sobely[1][0]=0; sobely[1][1]=-1; sobely[1][2]=0; sobely[1][3]=1; sobely[1][4]=0;
  sobely[2][0]=0; sobely[2][1]=-2; sobely[2][2]=0; sobely[2][3]=2; sobely[2][4]=0;
  sobely[3][0]=0; sobely[3][1]=-1; sobely[3][2]=0; sobely[3][3]=1; sobely[3][4]=0;
  sobely[4][0]=0; sobely[4][1]=0;  sobely[4][2]=0; sobely[4][3]=0; sobely[4][4]=0;

  sobelx[0][0]=0; sobelx[0][1]=0;  sobelx[0][2]=0;  sobelx[0][3]=0;  sobelx[0][4]=0;
  sobelx[1][0]=0; sobelx[1][1]=1;  sobelx[1][2]=2;  sobelx[1][3]=1;  sobelx[1][4]=0;
  sobelx[2][0]=0; sobelx[2][1]=0;  sobelx[2][2]=0;  sobelx[2][3]=0;  sobelx[2][4]=0;
  sobelx[3][0]=0; sobelx[3][1]=-1; sobelx[3][2]=-2; sobelx[3][3]=-1; sobelx[3][4]=0;
  sobelx[4][0]=0; sobelx[4][1]=0;  sobelx[4][2]=0;  sobelx[4][3]=0;  sobelx[4][4]=0;

  laplace[0][0]=0; laplace[0][1]=0;  laplace[0][2]=0;  laplace[0][3]=0;  laplace[0][4]=0;
  laplace[1][0]=0; laplace[1][1]=0;  laplace[1][2]=1;  laplace[1][3]=0;  laplace[1][4]=0;
  laplace[2][0]=0; laplace[2][1]=1;  laplace[2][2]=-4; laplace[2][3]=1;  laplace[2][4]=0;
  laplace[3][0]=0; laplace[3][1]=0;  laplace[3][2]=1;  laplace[3][3]=0;  laplace[3][4]=0;
  laplace[4][0]=0; laplace[4][1]=0;  laplace[4][2]=0;  laplace[4][3]=0;  laplace[4][4]=0;

  hessianyy[0][0]= 0; hessianyy[0][1]= 0;  hessianyy[0][2]= 0;  hessianyy[0][3]=0;  hessianyy[0][4]=0;
  hessianyy[1][0]= 0; hessianyy[1][1]= 1;  hessianyy[1][2]=-2;  hessianyy[1][3]=1;  hessianyy[1][4]=0;
  hessianyy[2][0]= 0; hessianyy[2][1]= 2;  hessianyy[2][2]=-4;  hessianyy[2][3]=2;  hessianyy[2][4]=0;
  hessianyy[3][0]= 0; hessianyy[3][1]= 1;  hessianyy[3][2]=-2;  hessianyy[3][3]=1;  hessianyy[3][4]=0;
  hessianyy[4][0]= 0; hessianyy[4][1]= 0;  hessianyy[4][2]= 0;  hessianyy[4][3]=0;  hessianyy[4][4]=0;

  hessianxx[0][0]=0; hessianxx[0][1]=0;  hessianxx[0][2]=0;  hessianxx[0][3]=0;  hessianxx[0][4]=0;
  hessianxx[1][0]=0; hessianxx[1][1]=1;  hessianxx[1][2]=2;  hessianxx[1][3]=1;  hessianxx[1][4]=0;
  hessianxx[2][0]=0; hessianxx[2][1]=-2; hessianxx[2][2]=-4; hessianxx[2][3]=-2; hessianxx[2][4]=0;
  hessianxx[3][0]=0; hessianxx[3][1]=1;  hessianxx[3][2]=2;  hessianxx[3][3]=1;  hessianxx[3][4]=0;
  hessianxx[4][0]=0; hessianxx[4][1]=0;  hessianxx[4][2]=0;  hessianxx[4][3]=0;  hessianxx[4][4]=0;

  hessianxy[0][0]=0; hessianxy[0][1]= 0; hessianxy[0][2]=0;  hessianxy[0][3]=0;  hessianxy[0][4]=0;
  hessianxy[1][0]=0; hessianxy[1][1]=-4; hessianxy[1][2]=0;  hessianxy[1][3]=4;  hessianxy[1][4]=0;
  hessianxy[2][0]=0; hessianxy[2][1]= 0; hessianxy[2][2]=0;  hessianxy[2][3]=0;  hessianxy[2][4]=0;
  hessianxy[3][0]=0; hessianxy[3][1]= 4; hessianxy[3][2]=0;  hessianxy[3][3]=-4; hessianxy[3][4]=0;
  hessianxy[4][0]=0; hessianxy[4][1]= 0; hessianxy[4][2]=0;  hessianxy[4][3]= 0; hessianxy[4][4]=0;
}

/*---------------------------------------------------------------------------*/
/*  Draw routines of points, lines and arrows                                */
/*---------------------------------------------------------------------------*/

void draw_point(int a,int xp,int yp,int r,int g,int b){
  ptr=outputarray[a]+3*((-yp)*x_size+xp); *ptr++=r; *ptr++=g; *ptr=b; 
}

void draw_line(int a,int xa,int ya, int xb,int yb,int r,int g,int b){
   int l; int L;
   int lx,ly; 
   if (xa<0) { xa=0; } if (xa> x_size) { xa= x_size; } if (xb< 1) { xb=0; } if (xb> x_size) { xb= x_size; } 
   if (ya>0) { ya=0; } if (ya<-y_size) { ya=-y_size; } if (yb>-1) { yb=0; } if (yb<-y_size) { yb=-y_size; } 
   L=abs(xb-xa)+abs(yb-ya); 
   if (L==0) { draw_point(a,xa,ya,r,g,b); } else {  
    for (l=0; l<=L; l++) { 
      lx=floor(1.0*l*(xb-xa)/L); 
      ly=floor(1.0*l*(yb-ya)/L);
      draw_point(a,xa+lx,ya+ly,r,g,b); 
    }
   }
}

void draw_arrow(int a,int xa,int ya, int xb,int yb,int r,int g,int b){
   int vx,vy,wx,wy,xc,yc,xd,yd,xe,ye;
     draw_line(a,xa,ya,xb,yb,r,g,b);                   // base of arrow
     vx=xb-xa; vy=yb-ya;  wx=-vy; wy=vx;
     xc=floor((xa+1*xb)/2); yc=floor((ya+1*yb)/2);
     xd=floor(xc+3*wx/4);   yd=floor(yc+3*wy/4);
     xe=floor(xc-3*wx/4);   ye=floor(yc-3*wy/4);
     draw_line(a,xe,ye,xd,yd,r,g,b);                   // head of arrow
     draw_line(a,xb,yb,xd,yd,r,g,b);                   // is a triangle
     draw_line(a,xb,yb,xe,ye,r,g,b);                   //
}


long  chartoint(char cen) { long a=cen; if (a<0) a+=256; return(a); }  //  convert char to integer
int xwrap(int x){int xx=(x+x_size) % x_size; return(xx); }             //  wrapping functions
int ywrap(int y){int yy=(y+y_size) % y_size; return(yy); }             //  wrapping function
int huetored(float h){float r=0.0; int u; float f,q; u=floor(5.999*h); f=5.999*h-u; q=1-f;
      if(u==0 || u==5) {r=1.0;} if(u==1){r=q;} if(u==4){r=f;} return(floor(255.0*r));}
int huetogreen(float h){float r=0.0; int u; float f,q; u=floor(5.999*h); f=5.999*h-u; q=1-f;
      if(u==1 || u==2) {r=1.0;} if(u==3){r=q;} if(u==0){r=f;} return(floor(255.0*r));}
int huetoblue(float h){float r=0.0; int u; float f,q; u=floor(5.999*h); f=5.999*h-u; q=1-f;
      if(u==3 || u==4) {r=1.0;} if(u==5){r=q;} if(u==2){r=f;} return(floor(255.0*r));}
double angle(int xx, int yy) {     // compute arg(x+i y), in C atan takes values in (-pi/2,pi/2)
double ang;
  if(yy==0) { if (xx>0) {ang=0.0;  } else {ang=PI;}     } else {
  if(xx==0) { if (yy>0) {ang=PI/2; } else {ang=3*PI/2;} } else {
  if(xx>0)  { ang=atan(1.0*yy/xx);}  else {ang=atan(1.0*yy/xx) + PI;}}}
  ang+=alpha;
return(ang); }

/*---------------------------------------------------------------------------*/
/*  Gradient and curvature computations  in color space                      */
/*---------------------------------------------------------------------------*/

void curvature1(int a, int iii, int jjj){                   
  long    rxx,gxx,bxx;                      // store variables
  long    ryy,gyy,byy;                      // store variables
  long    rxy,gxy,bxy;                      // store variables
  int     u,v;                              // summing over patch 
  int     ii,jj;                            // indices for patch 
  char    *inptr;                           // for gradient computation
  long    rthresh1,gthresh1,bthresh1;


    for (u=-2; u<=2; u++) {
      for (v=-2; v<=2; v++) {
        ii=xwrap(iii+PS*u);jj=-ywrap(-jjj+PS*v);
        inptr=inputarray[a]+3*((-jj)*x_size+ii);
        rpatch[u+2][v+2]=chartoint(*inptr++);gpatch[u+2][v+2]=chartoint(*inptr++);bpatch[u+2][v+2]=chartoint(*inptr);
      }
      }

      // first derivatives
      rx=0; ry=0; gx=0; gy=0; bx=0; by=0;
      for (u=1; u<=3; u+=2) {
      for (v=1; v<=3; v+=2) {
        rx += sobelx[u][v]*rpatch[u][v];   gx+=sobelx[u][v]*gpatch[u][v];   bx+=sobelx[u][v]*bpatch[u][v];
        ry += sobely[u][v]*rpatch[u][v];   gy+=sobely[u][v]*gpatch[u][v];   by+=sobely[u][v]*bpatch[u][v];
      }
      }

      // second derivatives
      rxx=0; ryy=0; rxy=0; gxx=0; gyy=0; gxy=0; bxx=0; byy=0; bxy=0;
      for (u=1; u<=3; u+=1) {
      for (v=1; v<=3; v+=1) {
        rxx+=hessianxx[u][v]*rpatch[u][v];  gxx+=hessianxx[u][v]*gpatch[u][v]; bxx+=hessianxx[u][v]*bpatch[u][v];
        ryy+=hessianyy[u][v]*rpatch[u][v];  gyy+=hessianyy[u][v]*gpatch[u][v]; byy+=hessianyy[u][v]*bpatch[u][v];
        rxy+=hessianxy[u][v]*rpatch[u][v];  gxy+=hessianxy[u][v]*gpatch[u][v]; bxy+=hessianxy[u][v]*bpatch[u][v];
      }
      }

      // unsigned curvature, where |grad(f)|^3 is replaced by |grad(f)|^2 to get larger integers
      // it just works so much better than the Kitchen-Rosenfeld curvature 
      rthresh=(rx*rx+ry*ry); gthresh=(gx*gx+gy*gy); bthresh=(bx*bx+by*by); 
      if (rthresh==0) {redge=0;} else {redge=floor(abs((rxx*ry*ry-2*rxy*rx*ry+ryy*rx*rx)*1.0)/rthresh);}
      if (gthresh==0) {gedge=0;} else {gedge=floor(abs((gxx*gy*gy-2*gxy*gx*gy+gyy*gx*gx)*1.0)/gthresh);}
      if (bthresh==0) {bedge=0;} else {bedge=floor(abs((bxx*by*by-2*bxy*bx*by+byy*bx*bx)*1.0)/bthresh);}
}

void curvature(int a, int iii, int jjj){
  int     u,v;
  int     KKK = (2*KK+1)*(2*KK+1);
  int     rrthresh,ggthresh,bbthresh,rredge,ggedge,bbedge;
  long    rrx,rry,ggx,ggy,bbx,bby;

  rrthresh=0;ggthresh=0;bbthresh=0;rredge=0;ggedge=0;bbedge=0;
  rrx=0; ggx=0; bbx=0; rry=0; ggy=0; bby=0;

  for (u=-KK; u<=KK; u++) {
  for (v=-KK; v<=KK; v++) {
      curvature1(a,iii+u,jjj+v);
      rrthresh+=rthresh; ggthresh+=gthresh; bbthresh+=bthresh;
      rredge  +=redge;   ggedge  +=gedge;   bbedge  +=bedge;
      rrx     +=rx;      ggx     +=gx;      bbx     +=bx;
      rry     +=ry;      ggy     +=gy;      bby     +=by;
  }
  }

  rthresh=floor(rrthresh/KKK); gthresh=floor(ggthresh/KKK); bthresh=floor(bbthresh/KKK);
  redge=floor(rredge/KKK);     gedge  =floor(ggedge/KKK);   bedge=floor(bbedge/KKK);
  rx=rrx/KKK; gx=ggx/KKK;      bx=bbx/KKK; ry=rry/KKK;      gy=ggy/KKK; by=bby/KKK;

  rang = angle(rx,ry);  if( rx==0 && ry==0) { rcos=0.0; rsin=0.0;} else { rcos = cos(rang); rsin = sin(rang); }
  gang = angle(rx,ry);  if( gx==0 && gy==0) { gcos=0.0; gsin=0.0;} else { gcos = cos(gang); gsin = sin(gang); }
  bang = angle(rx,ry);  if( bx==0 && by==0) { bcos=0.0; bsin=0.0;} else { bcos = cos(bang); bsin = sin(bang); }
}

/*---------------------------------------------------------------------------*/
/*  Side track fancy pictures, where red,green and blue pens along levels    */
/*---------------------------------------------------------------------------*/

void drawhamiltonian(a){
  int     i,j,k;                            // indices
  int     ii,jj,ii0,jj0,ii1,jj1;            // pixel positions
  double  rii,gii,bii,rjj,gjj,bjj;          // real positions

    for(j = -bb; ( (j > (-y_size+bb)) ); j-=K){
    for(i =  bb; ( (i < ( x_size-bb)) ); i+=K){

      curvature(a,i,j);
      rii=1.0*i; rjj=1.0*j; ii=i; jj=j; k=0;  // draw red flow lines
      while  ((rthresh>threshold) && (jj > (-y_size+bb)) && (jj<-bb) && (ii < (x_size-bb)) && (ii>bb) && k<stroke) {
        curvature(a,ii,jj); k++;
        if (rthresh>threshold) {
          rii+=rcos; rjj+=rsin; ii0=ii; jj0=jj; ii=(int) rii; jj=(int) rjj;
          ii1= (int) RR*rcos; jj1=(int) RR*rsin;
          draw_line(a,ii0,jj0,ii,jj,255,20,20);
        }
      }

      curvature(a,i,j);
      gii=1.0*i; gjj=1.0*j; ii=i; jj=j; k=0;  // draw green flow lines
      while  ((gthresh>threshold) && (jj > (-y_size+bb)) && (jj<-bb) && (ii < (x_size-bb)) && (ii>bb) && k<stroke) {
        curvature(a,ii,jj); k++;
        if (gthresh>threshold) {
          gii+=gcos; gjj+=gsin; ii0=ii; jj0=jj; ii=(int) gii; jj=(int) gjj;
          ii1= (int) RR*gcos; jj1=(int) RR*gsin;
          draw_line(a,ii0,jj0,ii,jj,20,255,20);
        }
      }

      curvature(a,i,j);
      bii=1.0*i; bjj=1.0*j; ii=i; jj=j; k=0;  // draw blue flow lines
      while  ((bthresh>threshold) && (jj > (-y_size+bb)) && (jj<-bb) && (ii < (x_size-bb)) && (ii>bb) && k<stroke) {
        curvature(a,ii,jj); k++;
        if (bthresh>threshold) {
          bii+=bcos; bjj+=bsin; ii0=ii; jj0=jj; ii=(int) bii; jj=(int) bjj;
          ii1= (int) RR*bcos; jj1=(int) RR*bsin;
          draw_line(a,ii0,jj0,ii,jj,20,20,255); 
        }
      }

   }                                          // end of i loop
   }                                          // end of j loop
}

/*---------------------------------------------------------------------------*/
/*  To gauge the curvature spectrum, we make a dry run                       */
/*---------------------------------------------------------------------------*/


void find_parameters(int a){                  // we have to know the maximal curvatures before
int     i,j;                                  // we can write the curvature files 
   fprintf(stderr,"Find parameters: ");

    maxthresh=0;     maxedge=0;
    maxrthresh[a]=1; maxgthresh[a]=1; maxbthresh[a]=1;
    maxredge[a]=1;   maxgedge[a]=1;   maxbedge[a]=1;

    for (j = -bb; ( (j > (-y_size+bb)) ); j-=K){                    // go through the entire picture 
    for (i =  bb; ( (i < ( x_size-bb)) ); i+=K){

      curvature1(a,i,j);

      if (rthresh>maxrthresh[a]) { maxrthresh[a]=rthresh; }         // maximal red green blue gradients
      if (gthresh>maxgthresh[a]) { maxgthresh[a]=gthresh; }
      if (bthresh>maxbthresh[a]) { maxbthresh[a]=bthresh; }

      if (redge>maxredge[a]) { maxredge[a]=redge; }                 // maximal red green blue curvatures
      if (gedge>maxgedge[a]) { maxgedge[a]=gedge; }
      if (bedge>maxbedge[a]) { maxbedge[a]=bedge; }
    }
    }

   fprintf(stderr,"Curvatures %i %i %i \n",maxredge[a],maxgedge[a],maxbedge[a]);
}

/*---------------------------------------------------------------------------*/
/*  At the beginning, we have to decide which points to take                 */
/*---------------------------------------------------------------------------*/

void initial_vectorfield(int a){ 
  int     i,j,k,l,m,p,q,u,v;                // indices
  int     ii,jj,kk,ll;                      // points in picture
  int     i0,j0,i2,j2;                      // points in previous and next picture
  long    laplacer,laplaceg,laplaceb;       // Laplace values
  double  diff[NF], mindiff[NF];            // best differences
  long    r0,r1,r2,g0,g1,g2,b0,b1,b2;       // red green blue value of two pictures 
  int     minp[NF], minq[NF];               // the best translation vector 
  int     mx,my;                            // for matching
  int     dist;                             // length of line in taximetric
  int     cp;                               // point index

    fprintf(stderr,"Compute initial vector field. ");

    cp=0; 
    for(j = -bb; ( (j > (-y_size+bb)) ); j-=K){               // go through the entire picture 
    for(i =  bb; ( (i < ( x_size-bb)) ); i+=K){

      ptr=incurvaturearray[a]+3*((-j)*x_size+i);
      redge=chartoint(*ptr++); gedge=chartoint(*ptr++); bedge=chartoint(*ptr);

      if ( redge>scaledthresh  || gedge>scaledthresh || bedge>scaledthresh )   {    //  interesting pts

          mindiff[a-1]=10000.0; mindiff[a+1]=10000.0; 
          mx=xwrap(i); my=-ywrap(-j);
          if( cp<MP) { 
             x_match[cp][a]=mx; y_match[cp][a]=my; cp++; countpoints[a]=cp;   
          }
          minp[a-1]=mx; minq[a-1]=my;
          minp[a+1]=mx; minq[a+1]=my;

          for (q =  M; q >= -M; q--) {                       // begin passing over displacement vector
          for (p = -M; p <=  M; p++) {

          mx=xwrap(i+q); my=-ywrap(-j-p); diff[a-1]=0.0;  diff[a+1]=0.0; 

          for (v =  N; v >= -N; v--) {                       // begin summing over patch 
          for (u = -N; u <=  N; u++) {
              ii=xwrap(i+PS1*u);  jj=-ywrap(-j-PS1*v);
              kk=xwrap(i+PS1*u+q);ll=-ywrap(-j-PS1*v-p);

              ptr=inputarray[a]+3*((-jj)*x_size+ii);
              r1=chartoint(*ptr++); g1=chartoint(*ptr++); b1=chartoint(*ptr);

              ptr=inputarray[a-1]+3*((-ll)*x_size+kk);
              r2=chartoint(*ptr++); g2=chartoint(*ptr++); b2=chartoint(*ptr);
              diff[a-1]+= abs(r1-r2)+abs(g1-g2)+abs(b1-b2);

              ptr=inputarray[a+1]+3*((-ll)*x_size+kk);
              r2=chartoint(*ptr++); g2=chartoint(*ptr++); b2=chartoint(*ptr);
              diff[a+1]+= abs(r1-r2)+abs(g1-g2)+abs(b1-b2);
          }
          }                                                   // end summing over patch 

          if (diff[a-1]<mindiff[a-1]) { minp[a-1]=mx; minq[a-1]=my; mindiff[a-1]=diff[a-1]; } 
          if (diff[a+1]<mindiff[a+1]) { minp[a+1]=mx; minq[a+1]=my; mindiff[a+1]=diff[a+1]; } 

          }                          
          }                                                   // end passing over displacement vector

       x_match[cp][a+1]=minp[a+1]; y_match[cp][a+1]=minq[a+1];
       dist= abs(minp[a-1]-minp[a+1]) +  abs(minq[a-1]-minq[a+1]); 
       if (mindiff[a-1]>1 && mindiff[a+1]>1 && dist<x_size/2 && dist<y_size/2) { 
         draw_line(1,minp[a-1],minq[a-1],minp[a+1],minq[a+1],255,  0,  0); 
       } 

     }

   }
   }

   fprintf(stderr,"Number of points: %i \n", countpoints[a]);
   fprintf(stderr,"Tracking: ");
}                                           //  end initiate_vectorfield

/*---------------------------------------------------------------------------*/
/*  Now, we evolve the chosen points along the motion field                  */
/*---------------------------------------------------------------------------*/

void vectorfield(int a){ 
  int     i,j,k,l,m,p,q,u,v;                // indices
  int     ii,jj,kk,ll;                      // points in picture
  int     i0,j0,i2,j2;                      // points in previous and next picture
  long    laplacer,laplaceg,laplaceb;       // Laplace values
  double  thresh,edge;                      // average over colors
  double  diff[NF], mindiff[NF];            // best differences
  double  localthresh;                      // old threshold in neighborhood
  double  localedge;                        // old edge in neighborhood
  long    r0,r1,r2,g0,g1,g2,b0,b1,b2;       // red green blue value of two pictures 
  int     minp[NF], minq[NF];               // the best translations
  int     mx,my,mx0,my0;                    // for matching
  int     dist;                             // length of line in taximetric
  int     cp; 

  fprintf(stderr,".");
  for (cp=0; cp<countpoints[a]; cp++) {
    i = x_match[cp][a-1]; j=y_match[cp][a-1]; 

    ptr=incurvaturearray[a]+3*((-j)*x_size+i);
    redge=chartoint(*ptr++); gedge=chartoint(*ptr++); bedge=chartoint(*ptr);
    ptr=ingradientarray[a]+3*((-j)*x_size+i);
    rthresh=chartoint(*ptr++); gthresh=chartoint(*ptr++); bthresh=chartoint(*ptr);

    mx=xwrap(i); my=-ywrap(-j);  mx0=mx; my0=my; 
    minp[a]=mx; minq[a]=my; mindiff[a]=10000.0; 

          for (q =  M; q >= -M; q--) {                       // begin passing over displacement vector
          for (p = -M; p <=  M; p++) {

          mx=xwrap(i+q); my=-ywrap(-j-p); diff[a]=0.0; 

          for (v =  N; v >= -N; v--) {                       // begin summing over patch 
          for (u = -N; u <=  N; u++) {
              ii=xwrap(i+PS1*u);  jj=-ywrap(-j-PS1*v);
              kk=xwrap(i+PS1*u+q);ll=-ywrap(-j-PS1*v-p);
              ptr=inputarray[a-1]+3*((-jj)*x_size+ii);
              r1=chartoint(*ptr++); g1=chartoint(*ptr++); b1=chartoint(*ptr);
              ptr=inputarray[a]+3*((-ll)*x_size+kk);
              r2=chartoint(*ptr++); g2=chartoint(*ptr++); b2=chartoint(*ptr);
              diff[a]+= abs(r1-r2)+abs(g1-g2)+abs(b1-b2);
          }
          }                                                   // end summing over patch 

          if (diff[a]<mindiff[a]) { minp[a]=mx; minq[a]=my; mindiff[a]=diff[a]; } 

          }                          
          }                                                   // end passing over displacement vector

       x_match[cp][a]=minp[a]; y_match[cp][a]=minq[a];
       dist= abs(mx0-minp[a]) + abs(my0-minq[a]);

       if (mindiff[a]>1 && dist<x_size/2  && dist<y_size/2) {
          draw_line(a,mx0,my0,minp[a],minq[a],255,  0,  0);
          draw_line(1,mx0,my0,minp[a],minq[a],255,  0,  0);
       } 
    }
}

/*---------------------------------------------------------------------------*/
/*  Produce sequences of file names for various external files               */
/*---------------------------------------------------------------------------*/
  

char *make_filename_a(int a, char *ext){             // for output pictures
  char *basename = "a";
  char *filename;
  filename = (char *) malloc((strlen(basename)+strlen(ext)+5)*sizeof(char));
  sprintf(filename,"%s%s%i.%s",basename,FILE_INDEX(a),ext);
  return(filename);
}

char *make_filename_s(int a, char *ext){           // for smooth pictures
  char *basename = "tmp/s";
  char *filename;
  filename = (char *) malloc((strlen(basename)+strlen(ext)+9)*sizeof(char));
  sprintf(filename,"%s%s%i.%s",basename,FILE_INDEX(a),ext);
  return(filename);
}

char *make_filename_c(int a, char *ext){           // for curvature pictures
  char *basename = "tmp/c";
  char *filename;
  filename = (char *) malloc((strlen(basename)+strlen(ext)+9)*sizeof(char));
  sprintf(filename,"%s%s%i.%s",basename,FILE_INDEX(a),ext);
  return(filename);
}

char *make_filename_g(int a, char *ext){           // for curvature pictures
  char *basename = "tmp/g";
  char *filename;
  filename = (char *) malloc((strlen(basename)+strlen(ext)+9)*sizeof(char));
  sprintf(filename,"%s%s%i.%s",basename,FILE_INDEX(a),ext);
  return(filename);
}

void write_output(int a){                          // write the pnm outputfile
  int   i,j,k;
  char *ptr0;
  out[a]=fopen(outfilename[a],"wb");
  fprintf(out[a],"P6\n%i %i\n%i\n", x_size, y_size, c_size);  
    for(j = 0; j > -y_size; j--){
    for(i = 0; i <  x_size; i++){
      ptr0  = outputarray[a] + 3*((-j)*x_size+i);  
      fputc(*ptr0++,out[a]); fputc(*ptr0++,out[a]); fputc(*ptr0,out[a]);
    }
    }
  fclose(out[a]); 
}

/*---------------------------------------------------------------------------*/
/*  Write external files which are ppm pictures                              */
/*---------------------------------------------------------------------------*/

void write_smoothed(int a){                        // write the pnm smoothed file
  int   i,j,k;
  char *ptr0;
  smoothout[a]=fopen(smoothfilename[a],"wb");
  fprintf(smoothout[a],"P6\n%i %i\n%i\n", x_size, y_size, c_size);
    for(j = 0; j > -y_size; j--){
    for(i = 0; i <  x_size; i++){
      ptr0  = outputarray[a] + 3*((-j)*x_size+i);
      fputc(*ptr0++,smoothout[a]); 
      fputc(*ptr0++,smoothout[a]); 
      fputc(*ptr0,smoothout[a]);
    }
    }
  fclose(smoothout[a]);
}

void write_curvature(int a){                       // write the curvature file
  int   i,j,k;
  char *ptr0;
  outcurvature[a]=fopen(curvaturefilename[a],"wb");
  fprintf(outcurvature[a],"P6\n%i %i\n%i\n", x_size, y_size, c_size);
    for(j = 0; j > -y_size; j--){
    for(i = 0; i <  x_size; i++){
      ptr0  = outcurvaturearray[a] + 3*((-j)*x_size+i);
      fputc(*ptr0++,outcurvature[a]); 
      fputc(*ptr0++,outcurvature[a]); 
      fputc(*ptr0,outcurvature[a]);
    }
    }
  fclose(outcurvature[a]);
}



void write_gradient(int a){                        // write the gradients
  int   i,j,k;
  char *ptr0;
  outgradient[a]=fopen(gradientfilename[a],"wb");
  fprintf(outgradient[a],"P6\n%i %i\n%i\n", x_size, y_size, c_size);
    for(j = 0; j > -y_size; j--){
    for(i = 0; i <  x_size; i++){
      ptr0  = outgradientarray[a] + 3*((-j)*x_size+i);
      fputc(*ptr0++,outgradient[a]); 
      fputc(*ptr0++,outgradient[a]); 
      fputc(*ptr0,outgradient[a]);
    }
    }
  fclose(outgradient[a]);
}

/*---------------------------------------------------------------------------*/
/*  Smoothing and export of smoothed pictures onto files                    */
/*---------------------------------------------------------------------------*/

void smoothedpic(int a){                          // write smooth picture into inital input bin
char *inptr,*outptr;
int i,j,ii,jj,u,v,k;
long r2,g2,b2;
int r3,g3,b3;
  for(j = 0; j > -y_size; j--){
  for(i = 0; i <  x_size; i++){
      for (u=-2; u<=2; u++) {
      for (v=-2; v<=2; v++) {
        ii=xwrap(i+u);jj=-ywrap(-j+v);
        inptr=outputarray[a]+3*((-jj)*x_size+ii);
        rpatch[u+2][v+2]=chartoint(*inptr++);gpatch[u+2][v+2]=chartoint(*inptr++);bpatch[u+2][v+2]=chartoint(*inptr);
      }
      }
      r2=0; g2=0; b2=0;
      for (u=0; u<=4; u++) {
      for (v=0; v<=4; v++) {
        r2 += gaussian[u][v]*rpatch[u][v]; g2+=gaussian[u][v]*gpatch[u][v]; b2+=gaussian[u][v]*bpatch[u][v];
      }
      }
      r3 = floor(r2/gaussiansum); g3 = floor(g2/gaussiansum); b3 = floor(b2/gaussiansum);
      outptr    = inputarray[a] + 3*((-j)*x_size+i );
      *outptr++ = r3; *outptr++ = g3; *outptr++ = b3;
  }
  }
}

/*---------------------------------------------------------------------------*/
/*  Curvature computation and writing out to files                           */
/*---------------------------------------------------------------------------*/

void compute_curvature(int a, int b){         // write curvature picture into allocated curvature memory 
char *inptr,*curvatureptr,*gradientptr;
int i,j,ii,jj,u,v,k;
int  r2,g2,b2;
int  r3,g3,b3;
int  newrmaxthresh,newgmaxthresh,newbmaxthresh;     // adapt the thresholds
int  newrmaxedge,newgmaxedge,newbmaxedge; 

fprintf(stderr,"Compute curvatures and gradients:");
for (k=a; k<=b; k++) {
  fprintf(stderr," %i ",k);
  newrmaxthresh=1; newgmaxthresh=1; newbmaxthresh=1; 
  newrmaxedge=1;   newgmaxedge=1;   newbmaxedge=1; 
  for(j = 0; j > -y_size; j--){
  for(i = 0; i <  x_size; i++){
        curvatureptr    = outcurvaturearray[k] + 3*((-j)*x_size+i );
        gradientptr     = outgradientarray[k]  + 3*((-j)*x_size+i );
        if ( (j > (-y_size+bb)) && (j<-bb) && (i < (x_size-bb)) && (i>bb) ){
          curvature1(k,i,j); 
          r2 = floor((255.0*rthresh)/maxrthresh[k]); g2 = floor((255.0*gthresh)/maxgthresh[k]); b2 = floor((255.0*bthresh)/maxbthresh[k]);
          r3 = floor((255.0*redge)/maxredge[k]);     g3 = floor((255.0*gedge)/maxgedge[k]);     b3 = floor((255.0*bedge)/maxbedge[k]);
          *gradientptr++  = r2;  *gradientptr++  = g2; *gradientptr++ = b2;
          *curvatureptr++ = r3;  *curvatureptr++ = g3; *curvatureptr++ = b3;
        } else  {
          *curvatureptr++ = 0; *curvatureptr++ = 0; *curvatureptr++ = 0;
          *gradientptr++  = 0; *gradientptr++  = 0; *gradientptr++  = 0;
        }
       if (rthresh>newrmaxthresh) { newrmaxthresh=rthresh; } 
       if (gthresh>newgmaxthresh) { newgmaxthresh=gthresh; } 
       if (bthresh>newbmaxthresh) { newbmaxthresh=bthresh; } 
       if (redge>newrmaxedge) { newrmaxedge=redge; } 
       if (gedge>newgmaxedge) { newgmaxedge=gedge; } 
       if (bedge>newbmaxedge) { newbmaxedge=bedge; } 
  }
  }
  if (k<b) {
    maxrthresh[k+1]=newrmaxthresh; maxgthresh[k+1]=newgmaxthresh; maxbthresh[k+1]=newbmaxthresh; 
    maxredge[k+1]=newrmaxedge; maxgedge[k+1]=newgmaxedge; maxbedge[k+1]=newbmaxedge; 
  }

  write_curvature(k);
  write_gradient(k);
}
fprintf(stderr,"\n");
}

/*---------------------------------------------------------------------------*/
/*  Shuffling around some memory  for example to iterate smoothing           */
/*---------------------------------------------------------------------------*/

void copypic(int a){                  // copy the picture from input to output
char *inptr,*outptr; 
int i,j,k;
  for(j = 0; j > -y_size; j--){
  for(i = 0; i <  x_size; i++){
      inptr     = inputarray[a]  + 3*((-j)*x_size+i );
      outptr    = outputarray[a] + 3*((-j)*x_size+i );
      *outptr++ = *inptr++; *outptr++ = *inptr++; *outptr++ = *inptr++; 
  }
  }
}

void copyallpic(int a, int b){
int k;
for (k=a; k<=b; k++) { copypic(k); }
}

void iterated_smoothedpic(int a, int b){
int k,l; 
fprintf(stderr,"Produce smoothed pictures:  ");
for (l=a; l<=b; l++) { 
   fprintf(stderr," %i ",l);
   for (k=0; k<S; k++) {                       // smooth S times, switching forth and back from input to output
        smoothedpic(l);
        copypic(l);
    }
   write_smoothed(l);
}
fprintf(stderr,"\n");
}


/*---------------------------------------------------------------------------*/
/*  End in various files                                                     */
/*---------------------------------------------------------------------------*/


void readpics(int a, int b){
char buffer[1025];
char *ptr;
int i,k;
  fprintf(stderr,"Read in original pictures:");
  for (k=a; k<=b; k++) {
    fprintf(stderr,".");
    fgets(buffer,1025,in[k]);
    if(strncmp(buffer,"P6",2)){ fprintf(stderr,"Error (need PPM raw)\n"); exit(1); }
    do fgets(buffer,1025,in[k]); while(*buffer == '#');
    x_size = atoi(strtok(buffer," ")); y_size = atoi(strtok(NULL," "));
    fgets(buffer,1025,in[k]); c_size = atoi(buffer);
    if((inputarray[k]= (char *) malloc(3*x_size*y_size*sizeof(char)))==NULL){ fprintf(stderr,"Err\n"); exit(1);}
    i = 0; ptr = inputarray[k];
    while(!feof(in[k]) && i<3*x_size*y_size){ *ptr++ = fgetc(in[k]); i++;} fclose(in[k]);
    if(i<x_size*y_size){ fprintf(stderr,"Premature end of image\n"); exit(1); }
  }
  fprintf(stderr,"\n");
}

void readsmoothedpics(int a, int b){
char buffer[1025];
char *ptr;
int i,k;
  if (preprocessed==1) {                       // only if preprocessing has happened 
  fprintf(stderr,"Read in smoothed pictures:");
  for (k=a; k<=b; k++) {
    fprintf(stderr,".");
    fgets(buffer,1025,insmoothed[k]);
    if(strncmp(buffer,"P6",2)){ fprintf(stderr,"Error (need PPM raw)\n"); exit(1); }
    do fgets(buffer,1025,insmoothed[k]); while(*buffer == '#');
    x_size = atoi(strtok(buffer," ")); y_size = atoi(strtok(NULL," "));
    fgets(buffer,1025,insmoothed[k]); c_size = atoi(buffer);
    if((insmoothedarray[k]= (char *) malloc(3*x_size*y_size*sizeof(char)))==NULL){ fprintf(stderr,"Err\n"); exit(1);}
    i = 0; ptr = insmoothedarray[k];
    while(!feof(insmoothed[k]) && i<3*x_size*y_size){ *ptr++ = fgetc(insmoothed[k]); i++;} fclose(insmoothed[k]);
    if(i<x_size*y_size){ fprintf(stderr,"Premature end of image\n"); exit(1); }
  }
  }
  fprintf(stderr,"\n");
}

void readcurvaturepics(int a, int b){
char buffer[1025];
char *ptr;
int i,k;
  if (preprocessed==1) {                       // only if preprocessing has happened 
  fprintf(stderr,"Read in curvature pictures:");
  for (k=a; k<=b; k++) {
    fprintf(stderr,".");
    fgets(buffer,1025,incurvature[k]);
    if(strncmp(buffer,"P6",2)){ fprintf(stderr,"Error (need PPM raw)\n"); exit(1); }
    do fgets(buffer,1025,incurvature[k]); while(*buffer == '#');
    x_size = atoi(strtok(buffer," ")); y_size = atoi(strtok(NULL," "));
    fgets(buffer,1025,incurvature[k]); c_size = atoi(buffer);
    if((incurvaturearray[k]= (char *) malloc(3*x_size*y_size*sizeof(char)))==NULL){ fprintf(stderr,"Err\n"); exit(1);}
    i = 0; ptr = incurvaturearray[k];
    while(!feof(incurvature[k]) && i<3*x_size*y_size){ *ptr++ = fgetc(incurvature[k]); i++;} fclose(incurvature[k]);
    if(i<x_size*y_size){ fprintf(stderr,"Premature end of image\n"); exit(1); }
  }
  }
  fprintf(stderr,"\n");
}

void readgradientpics(int a, int b){
char buffer[1025];
char *ptr;
int i,k;
  if (preprocessed==1) {                       // only if preprocessing has happened
  fprintf(stderr,"Read in gradient pictures:");
  for (k=a; k<=b; k++) {
    fprintf(stderr,".");
    fgets(buffer,1025,ingradient[k]);
    if(strncmp(buffer,"P6",2)){ fprintf(stderr,"Error (need PPM raw)\n"); exit(1); }
    do fgets(buffer,1025,ingradient[k]); while(*buffer == '#');
    x_size = atoi(strtok(buffer," ")); y_size = atoi(strtok(NULL," "));
    fgets(buffer,1025,ingradient[k]); c_size = atoi(buffer);
    if((ingradientarray[k]= (char *) malloc(3*x_size*y_size*sizeof(char)))==NULL){ fprintf(stderr,"Err\n"); exit(1);}
    i = 0; ptr = ingradientarray[k];
    while(!feof(ingradient[k]) && i<3*x_size*y_size){ *ptr++ = fgetc(ingradient[k]); i++;} fclose(ingradient[k]);
    if(i<x_size*y_size){ fprintf(stderr,"Premature end of image\n"); exit(1); }
  }
  }
  fprintf(stderr,"\n");
}

/*---------------------------------------------------------------------------*/
/*  Memory allocation for various picture arrays (movies)                    */
/*---------------------------------------------------------------------------*/


void memoryalloc(int a, int b){
int k; 
fprintf(stderr,"Allocate memory:\n");
  for (k=a; k<=b; k++) {
    if((outputarray[k]       = (char *) malloc(3*x_size*y_size*sizeof(char)))==NULL){fprintf(stderr,"Err\n");exit(1);}
    if((outcurvaturearray[k] = (char *) malloc(3*x_size*y_size*sizeof(char)))==NULL){fprintf(stderr,"Err\n");exit(1);}
    if((outgradientarray[k]  = (char *) malloc(3*x_size*y_size*sizeof(char)))==NULL){fprintf(stderr,"Err\n");exit(1);}
  }
}

/*---------------------------------------------------------------------------*/
/*  Merging points is nifty but makes export of paths a bit tricky           */
/*---------------------------------------------------------------------------*/


void mergepoints(int a){                    // merge points with similar position and velocity 
                                            // and where curvature values have not deteriorated a lot 
int cp, newcp, cpl;                         // countpoint index
int take;                                   // take the point (1) or not (0)
int  x1_match[MP], y1_match[MP];            // point position
int  v1_match[MP], w1_match[MP];            // point velocity 
double  thresh,edge;                        // average over colors


for (cp=0; cp<countpoints[a]; cp++)  {      // copy points over and compute velocities
   x1_match[cp]=x_match[cp][a];   
   y1_match[cp]=y_match[cp][a]; 
   v1_match[cp]=x_match[cp][a-1]; 
   w1_match[cp]=y_match[cp][a-1]; 
}

newcp=0; 
for (cp=0; cp<countpoints[a]; cp++) {
    take=1; cpl=0; 
    curvature(a,x1_match[cp],y1_match[cp]);
    thresh=(rthresh+gthresh+bthresh)/3; edge=(redge+gedge+bedge)/3;
    if  ( thresh<threshold/3  || edge < ethreshold/3 )  { take =0; }
    while (cpl<newcp && take==1) { 
       if (
           ( abs(x1_match[cpl]-x_match[cp][a])<merge) &&
           ( abs(y1_match[cpl]-y_match[cp][a])<merge) &&
           ( abs(v1_match[cpl]-x_match[cp][a-1])<merge) &&
           ( abs(w1_match[cpl]-y_match[cp][a-1])<merge) 
         ) { take=0; } cpl++; 
    }
    if (take==1) { 
      x_match[newcp][a]=x1_match[cp]; y_match[newcp][a]=y1_match[cp];   newcp++; 
    }
}
countpoints[a]=newcp; 
fprintf(stderr,"New number of points for %i frame: %i\n", a,countpoints[a]);
}

/*---------------------------------------------------------------------------*/
/*  Write the tracked points onto a file                                     */
/*---------------------------------------------------------------------------*/

void write_points(int aa1, int bb1){
FILE *res;                                  // outputfile
int  k,cp;                                  // index for counting points
char *match_file="points.dat";              // name of output file
  res=fopen(match_file,"w");
  fprintf(res,"0 0 0 0 %i %i \n",x_size,y_size);    // print first the size of the picture
  for (k=aa1; k<=bb1; k++) {
    for (cp=1; cp<countpoints[k]; cp++) { // camera number k, first frame aa1, last frame bb1, point number k, xcoord, ycoord
      fprintf(res,"%i %i %i %i %i %i \n",k,aa1,bb1,cp,x_match[cp][k],-y_match[cp][k]);
    }
  }
  fclose(res);
}

/*---------------------------------------------------------------------------*/
/*  Various initializations like file names or checks that pictures exist    */
/*---------------------------------------------------------------------------*/

void initializations(int a, int b) {
int k; 
if (b>NF) {fprintf(stderr,"Use less pictures\n\n"); }
if (fopen("tmp/c1000.pnm","rb")==NULL) { preprocessed=0; } else { preprocessed=1; }
if (fopen("tmp/s1000.pnm","rb")==NULL) { smootheddone=0; } else { smootheddone=1; }

for (k=a; k<=b; k++) {
    outfilename[k]      =make_filename_a(1000+k,"pnm");
    smoothfilename[k]   =make_filename_s(1000+k,"pnm");
    curvaturefilename[k]=make_filename_c(1000+k,"pnm");
    gradientfilename[k] =make_filename_g(1000+k,"pnm");

    if((in[k] = fopen(infilename[k],"rb")) == NULL){              // check  for input files
        fprintf(stderr,"Error %s\n", infilename[k]); exit(1);
    }

    if (preprocessed==1) {                                        // check for curvature files
      if((incurvature[k] = fopen(curvaturefilename[k],"rb")) == NULL){
         fprintf(stderr,"Error %s\n", curvaturefilename[k]); exit(1);
      }

                                                                 //  smoothed files 
      if((insmoothed[k] = fopen(smoothfilename[k],"rb")) == NULL){
         fprintf(stderr,"Error %s\n", smoothfilename[k]); exit(1);
      }

                                                                  //  ... and gradient files
      if((ingradient[k] = fopen(gradientfilename[k],"rb")) == NULL){
         fprintf(stderr,"Error %s\n", gradientfilename[k]); exit(1);
      }

    }
}
}

/*---------------------------------------------------------------------------*/
/*  Do preprosessing                                                         */
/*---------------------------------------------------------------------------*/

void smoothpictures(int a, int b) {
int k; 

if (smootheddone==0) {                          // preprocessing is necessary
    iterated_smoothedpic(a,b);                  // computed smooth pictures
} else { fprintf(stderr,"No smoothing necessary\n");  }
}


void preprocesscurvature(int a, int b) {
int k; 
if (preprocessed==0) {                          // preprocessing is necessary
    find_parameters(a);                         // determine parameters from first picture
    for (k=a+1; k<=b; k++) { 
         maxrthresh[k]= maxrthresh[a]; maxgthresh[k]= maxgthresh[a]; maxbthresh[k]= maxbthresh[a]; 
         maxredge[k]  = maxredge[a];   maxgedge[k]  = maxgedge[a];   maxbedge[k]  = maxbedge[a];
    }
    compute_curvature(a,b);                     // compute the curvatures 
  } else { fprintf(stderr,"No curvature computations necessary\n"); }
}

void processvectorfields(int a, int b) {        // the main tracking procedure
int k; 
if (preprocessed==1) { 
  initial_vectorfield(a+1);                     // start vector field and point choice
  for (k=a+2; k<b; k++) {
     vectorfield(k);
     countpoints[k]  =countpoints[k-1]; 
     countpoints[k+1]=countpoints[k-1]; 
     // mergepoints(k);
     write_output(k);
  }
  fprintf(stderr,"\n");
}
}

/*---------------------------------------------------------------------------*/
/*  The main file                                                            */
/*---------------------------------------------------------------------------*/

int main(int argc, char *argv[]){
int a,b,k; 
  a=0; b=argc-2; 
  if(b<=NF) { for (k=a; k<=b; k++) { infilename[k]=argv[k+1]; } }
  initializations(a,b);                       // make filenames 
  readpics(a,b);                              // read in the pictures 
  readsmoothedpics(a,b);                      // read in the smoothed pictures if exist
  readcurvaturepics(a,b);                     // read in the curvature pictures if exist
  readgradientpics(a,b);                      // read in the gradient pictures if existent
  memoryalloc(a,b);                           // memory allocation 
  definematrices();                           // define various matrices
  copyallpic(a,b);                            // copy pictures over 
  smoothpictures(a,b);                        // smooth pictures if necessary
  preprocesscurvature(a,b);                   // compute curvature values if necessary
  processvectorfields(a,b);                   // track the points along the vector fields
  write_points(a+3,b-1);                      // write points to text file 
  write_output(1);                            // write output file with trajectories
  return 0;                                   
}