/* code for computing absolute distance function of a 2D shape */
/* the positon of the initial front and its position as it advances are
marked as points on grid lines */

/* input: image: piecewise constant /*

/* The boundary of the shape generically crosses gridlines
   between nodes */

/* output: absolute distance function  */


#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <fcntl.h>
#include <unistd.h>
#include <sioux.h>
#include <MacTypes.h>
#include <Fonts.h>
#include <math.h>

#define MAX(A,B) (A>B) ? (A) : (B)
#define MIN(A,B) (A<B) ? (A) : (B)
#define ABS(A) (A>=0) ? (A) : (-A)

/* imagesizes:
   ccSchzSchz:  376H x 796W; sample count: 35
   ccSchzNorms: 364H x 784W; sample count: 48
   ccNorms: 324H x 780W; sample count: 20 
   ccKids : 400H x 896W; sample count: 30
   ccDrugs: 324H x 708W; sample count: 14
 */

#define sample_count 14

#define ydim 324 /* image height - dimension for i */
#define xdim 708 /* image width - dimension for j */

#define MaxDist 10000

typedef struct element
{
	float dist;
	int row;
	int col;
} HeapNode;

FILE *infile,*outfile;
char inname[128],outname[128],dir_name[128];


/* row and col are indices Tmat, not in pic */

unsigned char pic[ydim][xdim],outpic[ydim][xdim];


int length, heaplength;
  
int IArray[ydim*xdim],JArray[ydim*xdim];
   /* IArray and JArray store the coordinates of 
   the bdry nodes the shape */

  unsigned char Cmat[ydim][xdim];
  int Bmat[ydim][xdim];
  float dist[ydim][xdim],Tmat[ydim+2][xdim+2];

/* Cmat marks 'far' points as 3, 'close' or 'trial points' as 2
and 'accepted pts' as 1.
far points = points not yet visited,
trial points: points next to the front and outside the front.
accepted pts = points inside the front.  */

/* Bmat stores backpointer from the grid points to their
index in the heap  */

/* Tmat stores d-values (=time when the front arrives at the node) */

  
  HeapNode H[ydim*xdim],TmpNode[1],Head[1];
  /* Note that H,TmpNode and Head are pointers to HeapNode */
/* array for storing the heap. Note that heap index = array index +1 */

  

void calc_dist();
void upheap(HeapNode *,int, HeapNode *);
void downheap(HeapNode *,int);
void locateboundary();
float upwinddist(int,int); 
void calc_dist();
void output_results();
void run();

void calc_dist()
{
  int i,j;
  int Iind, Jind,NIind,NJind,pos;
  /*Initialize */
  
    for(i=0;i<ydim;i++) for(j=0;j<xdim;j++) Cmat[i][j]=3;

    for(i=0;i<ydim+2;i++) for(j=0;j<xdim+2;j++) Tmat[i][j] = MaxDist;
    
    /*find the points neighboring the bdry and
    calculate their distances, store them in Tmat */
    /* bdry points are the grid nodes next to the initial front
    which generically crosses the grid lines between grid nodes */
  
    length = 0;

    locateboundary();
    
    /* Begin advancing the front */
    
    /* first, tag bdry points as accepted points */
    for(i=0;i<length;i++)
    {
      Iind = IArray[i];
      Jind = JArray[i];
      Cmat[Iind-1][Jind-1] = 1; /* point part of the front */
    }

    /* Start of the heap:  neighbors of bdry pts are initial trial points */
      heaplength = 0; /* number of elements on the heap */
      for(i=0;i<length;i++)
      {
      Iind = IArray[i];
      Jind = JArray[i];
      /*West Neighbor*/
      NIind = Iind;
      NJind = Jind-1;
      if(NJind > 0)
	if (Cmat[NIind-1][NJind-1] == 3)
	/*  its d-value is unassigned */
	  {
	    TmpNode->dist = upwinddist(NIind,NJind);
	    TmpNode->row = NIind;
	    TmpNode->col = NJind;
	    Cmat[NIind-1][NJind-1] = 2; /* a new trial point */
	    upheap(H,heaplength,TmpNode);
	    heaplength++;
	  }

      /*North Neighbor*/
      NIind = Iind-1;
      NJind = Jind;
      if(NIind > 0)
	if (Cmat[NIind-1][NJind-1] == 3)
	  {
	    TmpNode->dist = upwinddist(NIind,NJind);
	    TmpNode->row = NIind;
	    TmpNode->col = NJind;
	    Cmat[NIind-1][NJind-1] = 2;
	    upheap(H,heaplength,TmpNode);
	    heaplength++;
	  }


      /*East Neighbor*/
      NIind = Iind;
      NJind = Jind+1;
      if(NJind <= xdim)
      if (Cmat[NIind-1][NJind-1] == 3)
	  {
	    TmpNode->dist = upwinddist(NIind,NJind);
	    TmpNode->row = NIind;
	    TmpNode->col = NJind;
	    Cmat[NIind-1][NJind-1] = 2;
	    upheap(H,heaplength,TmpNode);
	    heaplength++;
	  }

      
      /*South Neighbor*/
      NIind = Iind+1;
      NJind = Jind;
      if(NIind <= ydim)
      if (Cmat[NIind-1][NJind-1] == 3)
	  {
	    TmpNode->dist = upwinddist(NIind,NJind);
	    TmpNode->row = NIind;
	    TmpNode->col = NJind;
	    Cmat[NIind-1][NJind-1] = 2;
	    upheap(H,heaplength,TmpNode);
	    heaplength++;
	  }

    }

   /*march the front forward */
   while(heaplength>0)
   {
     *Head = *H;
     downheap(H,heaplength);
     heaplength--;
     Iind = Head->row;
     Jind = Head->col;
     Cmat[Iind-1][Jind-1] = 1;
     Tmat[Iind][Jind] = Head->dist;

     /*West Neighbor*/
     NIind = Iind;
     NJind = Jind-1;
     if((NJind>0) && (Cmat[NIind-1][NJind-1]>1))
     {
	   TmpNode->dist = upwinddist(NIind,NJind);
	   TmpNode->row = NIind;
	   TmpNode->col = NJind;
	   /*Update the heap*/
	 if (Cmat[NIind-1][NJind-1]==2) /*already a trial point*/
	   {
	     pos = Bmat[NIind-1][NJind-1]; /* heap index of the point */
	     upheap(H,pos-1,TmpNode); /* pos-1 is the array index of the point */
	   }
	 else
	   {
	     Cmat[NIind-1][NJind-1] = 2;
	     upheap(H,heaplength,TmpNode);
	     heaplength++;
	   }  
     }
     /*North Neighbor*/
     NIind = Iind-1;
     NJind = Jind;
     if((NIind>0) && (Cmat[NIind-1][NJind-1]>1))
     {
	   TmpNode->dist = upwinddist(NIind,NJind);
	   TmpNode->row = NIind;
	   TmpNode->col = NJind;
	   
	 if (Cmat[NIind-1][NJind-1]==2) 
	   {
	     pos = Bmat[NIind-1][NJind-1];
	     upheap(H,pos-1,TmpNode);
	   }
	 else
	   {
	     Cmat[NIind-1][NJind-1] = 2;
	     upheap(H,heaplength,TmpNode);
	     heaplength++;
	   }
     }
     
     /*East Neighbor*/
     NIind = Iind;
     NJind = Jind+1;
     if((NJind<=xdim) && (Cmat[NIind-1][NJind-1]>1))
     {
	   TmpNode->dist = upwinddist(NIind,NJind);
	   TmpNode->row = NIind;
	   TmpNode->col = NJind;
	   
	 if (Cmat[NIind-1][NJind-1]==2) 
	   {
	     pos = Bmat[NIind-1][NJind-1];
	     upheap(H,pos-1,TmpNode);
	   }
	 else
	   {
	     Cmat[NIind-1][NJind-1] = 2;
	     upheap(H,heaplength,TmpNode);
	     heaplength++;
	   }
     }
     
     /*South Neighbor*/
     NIind = Iind+1;
     NJind = Jind;
     if((NIind<=ydim) && (Cmat[NIind-1][NJind-1]>1))
     {
	   TmpNode->dist = upwinddist(NIind,NJind);
	   TmpNode->row = NIind;
	   TmpNode->col = NJind;
	   
	 if (Cmat[NIind-1][NJind-1]==2) 
	   {
	     pos = Bmat[NIind-1][NJind-1];
	     upheap(H,pos-1,TmpNode);
	   }
	 else
	   {
	     Cmat[NIind-1][NJind-1] = 2;
	     upheap(H,heaplength,TmpNode);
	     heaplength++;
	   }
     }
   }

   for(i=0;i<ydim;i++) for(j=0;j<xdim;j++)
   {
      if (Tmat[i+1][j+1]<MaxDist)
      dist[i][j] = Tmat[i+1][j+1];
      else printf("(%d %d) unassigned\n",i,j);
   }
}



void locateboundary() 
{
  /* Locate the pts neighboring the actual bdry.
  Calculate their distance from the bdry. */
 
  int NIind, NJind,i,j;
  int  VP,VW,VN,VE,VS,SW,SN,SE,SS;
  float d;
  
  for(i=2;i<=ydim-1;i++) for(j=2;j<=xdim-1;j++)
  /* indices (i,j) belong to array Tmat, not pic */
  {
	VP = pic[i-1][j-1];
	/*West Neighbor*/
    NIind = i;
	NJind = j-1;
	VW = pic[NIind-1][NJind-1];
	/*North Neighbor*/
	NIind = i-1;
	NJind = j;
	VN = pic[NIind-1][NJind-1];
	/*East Neighbor*/
	NIind = i;
	NJind = j+1;
	VE = pic[NIind-1][NJind-1];
	/*South Neighbor*/
	NIind = i+1;
	NJind = j;
	VS = pic[NIind-1][NJind-1];

	SW = SN = SE = SS = 0;
    if(VW!=VP) SW = 1;
    if(VN!=VP) SN = 1;
	if(VE!=VP) SE = 1;
	if(VS!=VP) SS = 1;
	
	if (SW+SE+SN+SS==0) d = 0;
	else
	{
	  IArray[length] = i;
	  JArray[length] = j;
	  length++;
	  if ((SW+SE==2)||(SN+SS==2)||(SW+SE+SN+SS==1))
		d = 0.5;
	  else d = 0.35355;
	  Tmat[i][j] = d;
	}
  }
}

float upwinddist(int Iind,int Jind)
{
  /*Calculate the value at (Iind,Jind) using the upwind scheme*/
  float v1,v2,d;
  double w1,w2;

  v1 = MIN(Tmat[Iind][Jind-1],Tmat[Iind][Jind+1]);
  v2 = MIN(Tmat[Iind-1][Jind],Tmat[Iind+1][Jind]);
  w1 = MAX(v1,v2);
  w2 = MIN(v1,v2);
  if(w1-w2>=1) d = w2+1;
  else d = (w1+w2+sqrt(2-(w1-w2)*(w1-w2)))/2;

  return d;			
}



void upheap(HeapNode *init, int length,HeapNode *newnode)
{
    /* init points to the root node. length = heap array index of the hole.
    heap index numbering starts at 1. (i.e. the numbering of tree nodes.)
    The index of the heap array = heap index - 1 */
    
    /* For a heap structure, two objects are needed to describe it:
    the array of HeapNodes and the length of the array. 
    length indicates where the upheap should begin, which means that at 
	length+1 the heap array has hole and the upheap should begin there.

    The change of the number of nodes in the heap is maintained separately.
    We should not depend on either upheap or downheap to update it. So that
    the two function can be used flexibly.
	*/
	int flag,index,parentindex,Iind,Jind;

	index = length+1; /* heap index of the hole */

    flag = 1;

	do{
		parentindex =  (int)floor(index/2);
		if(parentindex > 0)
		{
			if(newnode->dist<(init+parentindex-1)->dist)
			/* move the parent down */
			{
				*(init+index-1) = *(init+parentindex-1);
				Iind = (init+parentindex-1)->row;
				Jind = (init+parentindex-1)->col;
				Bmat[Iind-1][Jind-1] = index; /* update backpointer */
				index = parentindex; /* new location of the hole */
			}
			else /* insert the newnode */
			{
				*(init+index-1) = *newnode;
				Iind = newnode->row;
				Jind = newnode->col;
				Bmat[Iind-1][Jind-1] = index;
				flag = 0;
			}
		}
		else /* arrived at the root node */
		{
			*(init+index-1) = *newnode;
			Iind = newnode->row;
			Jind = newnode->col;
			Bmat[Iind-1][Jind-1] = index;
			flag = 0;
		}

	} while(flag == 1);

}

void downheap(HeapNode *init,int length)
/* downheap is called after removing the root node. */
{
	int  index,Iind,Jind;
     
	index=1; /* heap index of the hole */

	while(length>=2*index) /* not a leaf node */
	{
		if(length>=(2*index+1)) /* both children present */
		{
			if((init+2*index-1)->dist <= (init+2*index)->dist)
			/* left child smaller than right child, promote left */
			{
				*(init+index-1) = *(init+2*index-1);
				Iind = (init+2*index-1)->row;
				Jind = (init+2*index-1)->col;
				Bmat[Iind-1][Jind-1] = index; /* update backpointer */
				index = 2*index; /* new location of the hole */
			}
			else /* promote right child */
			{
				*(init+index-1) = *(init+2*index);
				Iind = (init+2*index)->row;
				Jind = (init+2*index)->col;
				Bmat[Iind-1][Jind-1] = index;
				index = 2*index+1;
			}
		}
		else /* only one child present, always on the left */
		{
				*(init+index-1) = *(init+2*index-1);
				Iind = (init+2*index-1)->row;
				Jind = (init+2*index-1)->col;
				Bmat[Iind-1][Jind-1] = index;
				index = 2*index;
		}
	}

	if (index < length)
	/* hole at a leaf which is not the last leaf.
	   insert the last leaf node here and then migrate upward  */
	   upheap(init,index-1,init+length-1);
}


 void output_results()
 {
    /*
    int i,j;
    float x,mx=0,mn=MaxDist+1,scale;

    for (i=0;i<ydim;i++) for (j=0;j<xdim;j++)
    {
       x = dist[i][j];
       mx = MAX(mx,x);
       mn = MIN(mn,x);
    }
    x = mx-mn;
    if (x>1) scale=255/x;
    else scale = 1;
    for (i=0;i<ydim;i++) for (j=0;j<xdim;j++)
      outpic[i][j]=(int) ((dist[i][j]-mn)*scale);

     printf("Writing %s\n",outname);
     if ((outfile = fopen(outname,"w+")) == NULL)
     {
	   printf("Cannot create %s\n",outname);
	   exit(1);
     }
     fwrite(outpic,sizeof(unsigned char),xdim*ydim,outfile);
     fclose(outfile); 
     */
    
     
     printf("Writing %s\n",outname);
     if ((outfile = fopen(outname,"w+")) == NULL)
     {
	   printf("Cannot create %s\n",outname);
	   exit(1);
     }
     fwrite(dist,sizeof(float),xdim*ydim,outfile); 
     fclose(outfile); 
 } 
 



void run(void)
{
   int m,n;
   char str[100][128];
   char sstr[100][128];
   time_t systime;
   	   
   SIOUXSettings.asktosaveonclose = 0;
   systime = time(NULL);
   puts(ctime(&systime));
   
  /* FILE NAMES */
 /* Strings str are input files; sstr are output files */
 
       
   /* SCHIZOPHRENICS */
  /*
   strcpy(dir_name,"::CorpusCallosum:ccSchzSchz:");

   strcpy(str[0],"::CorpusCallosum:ccSchzSchz:1030_2_128.buchar");
   strcpy(str[1],"::CorpusCallosum:ccSchzSchz:1033_2_128.buchar");
   strcpy(str[2],"::CorpusCallosum:ccSchzSchz:1035_2_128.buchar");
   strcpy(str[3],"::CorpusCallosum:ccSchzSchz:1036_2_128.buchar");
   strcpy(str[4],"::CorpusCallosum:ccSchzSchz:1037_2_128.buchar");
   strcpy(str[5],"::CorpusCallosum:ccSchzSchz:1038_4_128.buchar");
   strcpy(str[6],"::CorpusCallosum:ccSchzSchz:468_6_128.buchar");
   strcpy(str[7],"::CorpusCallosum:ccSchzSchz:474_4_128.buchar");
   strcpy(str[8],"::CorpusCallosum:ccSchzSchz:476_6_128.buchar");
   strcpy(str[9],"::CorpusCallosum:ccSchzSchz:479_6_128.buchar");
   strcpy(str[10],"::CorpusCallosum:ccSchzSchz:491_4_128.buchar");
   strcpy(str[11],"::CorpusCallosum:ccSchzSchz:510_4_128.buchar");
   strcpy(str[12],"::CorpusCallosum:ccSchzSchz:511_4_128.buchar");
   strcpy(str[13],"::CorpusCallosum:ccSchzSchz:519_4_128.buchar");
   strcpy(str[14],"::CorpusCallosum:ccSchzSchz:578_4_128.buchar");
   strcpy(str[15],"::CorpusCallosum:ccSchzSchz:579_4_128.buchar");
   strcpy(str[16],"::CorpusCallosum:ccSchzSchz:580_4_128.buchar");
   strcpy(str[17],"::CorpusCallosum:ccSchzSchz:597_4_128.buchar");
   strcpy(str[18],"::CorpusCallosum:ccSchzSchz:614_4_128.buchar");
   strcpy(str[19],"::CorpusCallosum:ccSchzSchz:615_4_128.buchar");
   strcpy(str[20],"::CorpusCallosum:ccSchzSchz:631_4_128.buchar");
   strcpy(str[21],"::CorpusCallosum:ccSchzSchz:635_5_128.buchar");
   strcpy(str[22],"::CorpusCallosum:ccSchzSchz:647_4_128.buchar");
   strcpy(str[23],"::CorpusCallosum:ccSchzSchz:672_4_128.buchar");
   strcpy(str[24],"::CorpusCallosum:ccSchzSchz:677_4_128.buchar");
   strcpy(str[25],"::CorpusCallosum:ccSchzSchz:686_4_128.buchar");
   strcpy(str[26],"::CorpusCallosum:ccSchzSchz:689_4_128.buchar");
   strcpy(str[27],"::CorpusCallosum:ccSchzSchz:690_4_128.buchar");
   strcpy(str[28],"::CorpusCallosum:ccSchzSchz:772_4_128.buchar");
   strcpy(str[29],"::CorpusCallosum:ccSchzSchz:789_4_128.buchar");
   strcpy(str[30],"::CorpusCallosum:ccSchzSchz:790_4_128.buchar");
   strcpy(str[31],"::CorpusCallosum:ccSchzSchz:791_4_128.buchar");
   strcpy(str[32],"::CorpusCallosum:ccSchzSchz:792_4_128.buchar");
   strcpy(str[33],"::CorpusCallosum:ccSchzSchz:793_4_128.buchar");
   strcpy(str[34],"::CorpusCallosum:ccSchzSchz:853_4_128.buchar");
   
   strcpy(sstr[0],"::CorpusCallosum:ccSchzSchz:1030");
   strcpy(sstr[1],"::CorpusCallosum:ccSchzSchz:1033");
   strcpy(sstr[2],"::CorpusCallosum:ccSchzSchz:1035");
   strcpy(sstr[3],"::CorpusCallosum:ccSchzSchz:1036");
   strcpy(sstr[4],"::CorpusCallosum:ccSchzSchz:1037");
   strcpy(sstr[5],"::CorpusCallosum:ccSchzSchz:1038");
   strcpy(sstr[6],"::CorpusCallosum:ccSchzSchz:468");
   strcpy(sstr[7],"::CorpusCallosum:ccSchzSchz:474");
   strcpy(sstr[8],"::CorpusCallosum:ccSchzSchz:476");
   strcpy(sstr[9],"::CorpusCallosum:ccSchzSchz:479");
   strcpy(sstr[10],"::CorpusCallosum:ccSchzSchz:491");
   strcpy(sstr[11],"::CorpusCallosum:ccSchzSchz:510");
   strcpy(sstr[12],"::CorpusCallosum:ccSchzSchz:511");
   strcpy(sstr[13],"::CorpusCallosum:ccSchzSchz:519");
   strcpy(sstr[14],"::CorpusCallosum:ccSchzSchz:578");
   strcpy(sstr[15],"::CorpusCallosum:ccSchzSchz:579");
   strcpy(sstr[16],"::CorpusCallosum:ccSchzSchz:580");
   strcpy(sstr[17],"::CorpusCallosum:ccSchzSchz:597");
   strcpy(sstr[18],"::CorpusCallosum:ccSchzSchz:614");
   strcpy(sstr[19],"::CorpusCallosum:ccSchzSchz:615");
   strcpy(sstr[20],"::CorpusCallosum:ccSchzSchz:631");
   strcpy(sstr[21],"::CorpusCallosum:ccSchzSchz:635");
   strcpy(sstr[22],"::CorpusCallosum:ccSchzSchz:647");
   strcpy(sstr[23],"::CorpusCallosum:ccSchzSchz:672");
   strcpy(sstr[24],"::CorpusCallosum:ccSchzSchz:677");
   strcpy(sstr[25],"::CorpusCallosum:ccSchzSchz:686");
   strcpy(sstr[26],"::CorpusCallosum:ccSchzSchz:689");
   strcpy(sstr[27],"::CorpusCallosum:ccSchzSchz:690");
   strcpy(sstr[28],"::CorpusCallosum:ccSchzSchz:772");
   strcpy(sstr[29],"::CorpusCallosum:ccSchzSchz:789");
   strcpy(sstr[30],"::CorpusCallosum:ccSchzSchz:790");
   strcpy(sstr[31],"::CorpusCallosum:ccSchzSchz:791");
   strcpy(sstr[32],"::CorpusCallosum:ccSchzSchz:792");
   strcpy(sstr[33],"::CorpusCallosum:ccSchzSchz:793");
   strcpy(sstr[34],"::CorpusCallosum:ccSchzSchz:853");
  */
   
 
 /* NORMALS: BASE LINE FOR SCHIZOPHRENICS */
 /*  
   strcpy(dir_name,"::CorpusCallosum:ccSchzNorms:");

   strcpy(str[0],"::CorpusCallosum:ccSchzNorms:1019_2_128.buchar");
   strcpy(str[1],"::CorpusCallosum:ccSchzNorms:1039_4_128.buchar");
   strcpy(str[2],"::CorpusCallosum:ccSchzNorms:1042_2_128.buchar");
   strcpy(str[3],"::CorpusCallosum:ccSchzNorms:1044_2_128.buchar");
   strcpy(str[4],"::CorpusCallosum:ccSchzNorms:1045_2_128.buchar");
   strcpy(str[5],"::CorpusCallosum:ccSchzNorms:1046_2_128.buchar");
   strcpy(str[6],"::CorpusCallosum:ccSchzNorms:1050_3_128.buchar");
   strcpy(str[7],"::CorpusCallosum:ccSchzNorms:1069_2_128.buchar");
   strcpy(str[8],"::CorpusCallosum:ccSchzNorms:1077_2_128.buchar");
   strcpy(str[9],"::CorpusCallosum:ccSchzNorms:1110_2_128.buchar");
   strcpy(str[10],"::CorpusCallosum:ccSchzNorms:1112_2_128.buchar");
   strcpy(str[11],"::CorpusCallosum:ccSchzNorms:1120_2_128.buchar");
   strcpy(str[12],"::CorpusCallosum:ccSchzNorms:1134_2_128.buchar");
   strcpy(str[13],"::CorpusCallosum:ccSchzNorms:1150_2_128.buchar");
   strcpy(str[14],"::CorpusCallosum:ccSchzNorms:1151_2_128.buchar");
   strcpy(str[15],"::CorpusCallosum:ccSchzNorms:1152_2_128.buchar");
   strcpy(str[16],"::CorpusCallosum:ccSchzNorms:1153_2_128.buchar");
   strcpy(str[17],"::CorpusCallosum:ccSchzNorms:1154_2_128.buchar");
   strcpy(str[18],"::CorpusCallosum:ccSchzNorms:1155_2_128.buchar");
   strcpy(str[19],"::CorpusCallosum:ccSchzNorms:1156_2_128.buchar");
   strcpy(str[20],"::CorpusCallosum:ccSchzNorms:1157_2_128.buchar");
   strcpy(str[21],"::CorpusCallosum:ccSchzNorms:1158_2_128.buchar");
   strcpy(str[22],"::CorpusCallosum:ccSchzNorms:1159_2_128.buchar");
   strcpy(str[23],"::CorpusCallosum:ccSchzNorms:1182_3_128.buchar");
   strcpy(str[24],"::CorpusCallosum:ccSchzNorms:521_4_130.buchar");
   strcpy(str[25],"::CorpusCallosum:ccSchzNorms:522_4_128.buchar");
   strcpy(str[26],"::CorpusCallosum:ccSchzNorms:538_4_128.buchar");
   strcpy(str[27],"::CorpusCallosum:ccSchzNorms:546_4_128.buchar");
   strcpy(str[28],"::CorpusCallosum:ccSchzNorms:594_5_128.buchar");
   strcpy(str[29],"::CorpusCallosum:ccSchzNorms:595_4_128.buchar");
   strcpy(str[30],"::CorpusCallosum:ccSchzNorms:632_4_128.buchar");
   strcpy(str[31],"::CorpusCallosum:ccSchzNorms:634_4_128.buchar");
   strcpy(str[32],"::CorpusCallosum:ccSchzNorms:638_4_128.buchar");
   strcpy(str[33],"::CorpusCallosum:ccSchzNorms:649_4_128.buchar");
   strcpy(str[34],"::CorpusCallosum:ccSchzNorms:652_4_128.buchar");
   strcpy(str[35],"::CorpusCallosum:ccSchzNorms:653_4_128.buchar");
   strcpy(str[36],"::CorpusCallosum:ccSchzNorms:654_4_128.buchar");
   strcpy(str[37],"::CorpusCallosum:ccSchzNorms:673_4_128.buchar");
   strcpy(str[38],"::CorpusCallosum:ccSchzNorms:674_4_128.buchar");
   strcpy(str[39],"::CorpusCallosum:ccSchzNorms:681_4_128.buchar");
   strcpy(str[40],"::CorpusCallosum:ccSchzNorms:688_4_128.buchar");
   strcpy(str[41],"::CorpusCallosum:ccSchzNorms:794_4_128.buchar");
   strcpy(str[42],"::CorpusCallosum:ccSchzNorms:847_7_128.buchar");
   strcpy(str[43],"::CorpusCallosum:ccSchzNorms:851_6_128.buchar");
   strcpy(str[44],"::CorpusCallosum:ccSchzNorms:860_4_128.buchar");
   strcpy(str[45],"::CorpusCallosum:ccSchzNorms:871_4_128.buchar");
   strcpy(str[46],"::CorpusCallosum:ccSchzNorms:874_4_128.buchar");
   strcpy(str[47],"::CorpusCallosum:ccSchzNorms:880_4_128.buchar");

   strcpy(sstr[0],"::CorpusCallosum:ccSchzNorms:1019");
   strcpy(sstr[1],"::CorpusCallosum:ccSchzNorms:1039");
   strcpy(sstr[2],"::CorpusCallosum:ccSchzNorms:1042");
   strcpy(sstr[3],"::CorpusCallosum:ccSchzNorms:1044");
   strcpy(sstr[4],"::CorpusCallosum:ccSchzNorms:1045");
   strcpy(sstr[5],"::CorpusCallosum:ccSchzNorms:1046");
   strcpy(sstr[6],"::CorpusCallosum:ccSchzNorms:1050");
   strcpy(sstr[7],"::CorpusCallosum:ccSchzNorms:1069");
   strcpy(sstr[8],"::CorpusCallosum:ccSchzNorms:1077");
   strcpy(sstr[9],"::CorpusCallosum:ccSchzNorms:1110");
   strcpy(sstr[10],"::CorpusCallosum:ccSchzNorms:1112");
   strcpy(sstr[11],"::CorpusCallosum:ccSchzNorms:1120");
   strcpy(sstr[12],"::CorpusCallosum:ccSchzNorms:1134");
   strcpy(sstr[13],"::CorpusCallosum:ccSchzNorms:1150");
   strcpy(sstr[14],"::CorpusCallosum:ccSchzNorms:1151");
   strcpy(sstr[15],"::CorpusCallosum:ccSchzNorms:1152");
   strcpy(sstr[16],"::CorpusCallosum:ccSchzNorms:1153");
   strcpy(sstr[17],"::CorpusCallosum:ccSchzNorms:1154");
   strcpy(sstr[18],"::CorpusCallosum:ccSchzNorms:1155");
   strcpy(sstr[19],"::CorpusCallosum:ccSchzNorms:1156");
   strcpy(sstr[20],"::CorpusCallosum:ccSchzNorms:1157");
   strcpy(sstr[21],"::CorpusCallosum:ccSchzNorms:1158");
   strcpy(sstr[22],"::CorpusCallosum:ccSchzNorms:1159");
   strcpy(sstr[23],"::CorpusCallosum:ccSchzNorms:1182");
   strcpy(sstr[24],"::CorpusCallosum:ccSchzNorms:521");
   strcpy(sstr[25],"::CorpusCallosum:ccSchzNorms:522");
   strcpy(sstr[26],"::CorpusCallosum:ccSchzNorms:538");
   strcpy(sstr[27],"::CorpusCallosum:ccSchzNorms:546");
   strcpy(sstr[28],"::CorpusCallosum:ccSchzNorms:594");
   strcpy(sstr[29],"::CorpusCallosum:ccSchzNorms:595");
   strcpy(sstr[30],"::CorpusCallosum:ccSchzNorms:632");
   strcpy(sstr[31],"::CorpusCallosum:ccSchzNorms:634");
   strcpy(sstr[32],"::CorpusCallosum:ccSchzNorms:638");
   strcpy(sstr[33],"::CorpusCallosum:ccSchzNorms:649");
   strcpy(sstr[34],"::CorpusCallosum:ccSchzNorms:652");
   strcpy(sstr[35],"::CorpusCallosum:ccSchzNorms:653");
   strcpy(sstr[36],"::CorpusCallosum:ccSchzNorms:654");
   strcpy(sstr[37],"::CorpusCallosum:ccSchzNorms:673");
   strcpy(sstr[38],"::CorpusCallosum:ccSchzNorms:674");
   strcpy(sstr[39],"::CorpusCallosum:ccSchzNorms:681");
   strcpy(sstr[40],"::CorpusCallosum:ccSchzNorms:688");
   strcpy(sstr[41],"::CorpusCallosum:ccSchzNorms:794");
   strcpy(sstr[42],"::CorpusCallosum:ccSchzNorms:847");
   strcpy(sstr[43],"::CorpusCallosum:ccSchzNorms:851");
   strcpy(sstr[44],"::CorpusCallosum:ccSchzNorms:860");
   strcpy(sstr[45],"::CorpusCallosum:ccSchzNorms:871");
   strcpy(sstr[46],"::CorpusCallosum:ccSchzNorms:874");
   strcpy(sstr[47],"::CorpusCallosum:ccSchzNorms:880");
  */
  
 
 
 /* NORMALS */
 /* The first ten are females, the last ten males.*/
 
 /*
   strcpy(dir_name,"::CorpusCallosum:ccNorms:");

   strcpy(str[0],"::CorpusCallosum:ccNorms:100_23_127.buchar");
   strcpy(str[1],"::CorpusCallosum:ccNorms:16_3_134.buchar");
   strcpy(str[2],"::CorpusCallosum:ccNorms:17_3_132.buchar");
   strcpy(str[3],"::CorpusCallosum:ccNorms:191_3_133.buchar");
   strcpy(str[4],"::CorpusCallosum:ccNorms:202_3_131.buchar");
   strcpy(str[5],"::CorpusCallosum:ccNorms:205_3_120.buchar");
   strcpy(str[6],"::CorpusCallosum:ccNorms:1_24_130.buchar");
   strcpy(str[7],"::CorpusCallosum:ccNorms:2_4_131.buchar");
   strcpy(str[8],"::CorpusCallosum:ccNorms:4_8_128.buchar");
   strcpy(str[9],"::CorpusCallosum:ccNorms:5_8_131.buchar");

   strcpy(str[10],"::CorpusCallosum:ccNorms:110_3_127.buchar");
   strcpy(str[11],"::CorpusCallosum:ccNorms:111_2_122.buchar");
   strcpy(str[12],"::CorpusCallosum:ccNorms:112_2_121.buchar");
   strcpy(str[13],"::CorpusCallosum:ccNorms:11_3_131.buchar");
   strcpy(str[14],"::CorpusCallosum:ccNorms:12_3_137.buchar");
   strcpy(str[15],"::CorpusCallosum:ccNorms:13_3_127.buchar");
   strcpy(str[16],"::CorpusCallosum:ccNorms:15_3_137.buchar");
   strcpy(str[17],"::CorpusCallosum:ccNorms:6_10_131.buchar");
   strcpy(str[18],"::CorpusCallosum:ccNorms:7_8_132.buchar");
   strcpy(str[19],"::CorpusCallosum:ccNorms:8_4_133.buchar");

   strcpy(sstr[0],"::CorpusCallosum:ccNorms:100");
   strcpy(sstr[1],"::CorpusCallosum:ccNorms:16");
   strcpy(sstr[2],"::CorpusCallosum:ccNorms:17");
   strcpy(sstr[3],"::CorpusCallosum:ccNorms:191");
   strcpy(sstr[4],"::CorpusCallosum:ccNorms:202");
   strcpy(sstr[5],"::CorpusCallosum:ccNorms:205");
   strcpy(sstr[6],"::CorpusCallosum:ccNorms:1");
   strcpy(sstr[7],"::CorpusCallosum:ccNorms:2");
   strcpy(sstr[8],"::CorpusCallosum:ccNorms:4");
   strcpy(sstr[9],"::CorpusCallosum:ccNorms:5");

   strcpy(sstr[10],"::CorpusCallosum:ccNorms:110");
   strcpy(sstr[11],"::CorpusCallosum:ccNorms:111");
   strcpy(sstr[12],"::CorpusCallosum:ccNorms:112");
   strcpy(sstr[13],"::CorpusCallosum:ccNorms:11");
   strcpy(sstr[14],"::CorpusCallosum:ccNorms:12");
   strcpy(sstr[15],"::CorpusCallosum:ccNorms:13");
   strcpy(sstr[16],"::CorpusCallosum:ccNorms:15");
   strcpy(sstr[17],"::CorpusCallosum:ccNorms:6");
   strcpy(sstr[18],"::CorpusCallosum:ccNorms:7");
   strcpy(sstr[19],"::CorpusCallosum:ccNorms:8");
*/


/* KIDS */
/* The first 15 are females, the last 15 males. */

/*  strcpy(dir_name,"::CorpusCallosum:ccKids:"); 
  
   strcpy(str[0],"::CorpusCallosum:ccKids:132_5_123.buchar");
   strcpy(str[1],"::CorpusCallosum:ccKids:158_4_129.buchar");
   strcpy(str[2],"::CorpusCallosum:ccKids:161_5_126.buchar");
   strcpy(str[3],"::CorpusCallosum:ccKids:169_4_127.buchar");
   strcpy(str[4],"::CorpusCallosum:ccKids:267_6_135.buchar");
   strcpy(str[5],"::CorpusCallosum:ccKids:277_5_90.buchar");
   strcpy(str[6],"::CorpusCallosum:ccKids:278_5_124.buchar");
   strcpy(str[7],"::CorpusCallosum:ccKids:291_5_129.buchar");
   strcpy(str[8],"::CorpusCallosum:ccKids:292_6_127.buchar");
   strcpy(str[9],"::CorpusCallosum:ccKids:349_4_103.buchar");
   strcpy(str[10],"::CorpusCallosum:ccKids:376_4_128.buchar");
   strcpy(str[11],"::CorpusCallosum:ccKids:386_4_124.buchar");
   strcpy(str[12],"::CorpusCallosum:ccKids:392_4_131.buchar");
   strcpy(str[13],"::CorpusCallosum:ccKids:401_4_124.buchar");
   strcpy(str[14],"::CorpusCallosum:ccKids:9_12_132.buchar");

   strcpy(str[15],"::CorpusCallosum:ccKids:139_4_125.buchar");
   strcpy(str[16],"::CorpusCallosum:ccKids:140_4_125.buchar");
   strcpy(str[17],"::CorpusCallosum:ccKids:146_4_128.buchar");
   strcpy(str[18],"::CorpusCallosum:ccKids:159_4_130.buchar");
   strcpy(str[19],"::CorpusCallosum:ccKids:160_4_124.buchar");
   strcpy(str[20],"::CorpusCallosum:ccKids:295_4_129.buchar");
   strcpy(str[21],"::CorpusCallosum:ccKids:308_6_129.buchar");
   strcpy(str[22],"::CorpusCallosum:ccKids:309_7_127.buchar");
   strcpy(str[23],"::CorpusCallosum:ccKids:310_4_119.buchar");
   strcpy(str[24],"::CorpusCallosum:ccKids:312_4_135.buchar");
   strcpy(str[25],"::CorpusCallosum:ccKids:323_4_127.buchar");
   strcpy(str[26],"::CorpusCallosum:ccKids:330_5_125.buchar");
   strcpy(str[27],"::CorpusCallosum:ccKids:385_4_121.buchar");
   strcpy(str[28],"::CorpusCallosum:ccKids:387_4_130.buchar");
   strcpy(str[29],"::CorpusCallosum:ccKids:415_4_129.buchar");

   strcpy(sstr[0],"::CorpusCallosum:ccKids:132");
   strcpy(sstr[1],"::CorpusCallosum:ccKids:158");
   strcpy(sstr[2],"::CorpusCallosum:ccKids:161");
   strcpy(sstr[3],"::CorpusCallosum:ccKids:169");
   strcpy(sstr[4],"::CorpusCallosum:ccKids:267");
   strcpy(sstr[5],"::CorpusCallosum:ccKids:277");
   strcpy(sstr[6],"::CorpusCallosum:ccKids:278");
   strcpy(sstr[7],"::CorpusCallosum:ccKids:291");
   strcpy(sstr[8],"::CorpusCallosum:ccKids:292");
   strcpy(sstr[9],"::CorpusCallosum:ccKids:349");
   strcpy(sstr[10],"::CorpusCallosum:ccKids:376");
   strcpy(sstr[11],"::CorpusCallosum:ccKids:386");
   strcpy(sstr[12],"::CorpusCallosum:ccKids:392");
   strcpy(sstr[13],"::CorpusCallosum:ccKids:401");
   strcpy(sstr[14],"::CorpusCallosum:ccKids:9");


   strcpy(sstr[15],"::CorpusCallosum:ccKids:139");
   strcpy(sstr[16],"::CorpusCallosum:ccKids:140");
   strcpy(sstr[17],"::CorpusCallosum:ccKids:146");
   strcpy(sstr[18],"::CorpusCallosum:ccKids:159");
   strcpy(sstr[19],"::CorpusCallosum:ccKids:160");
   strcpy(sstr[20],"::CorpusCallosum:ccKids:295");
   strcpy(sstr[21],"::CorpusCallosum:ccKids:308");
   strcpy(sstr[22],"::CorpusCallosum:ccKids:309");
   strcpy(sstr[23],"::CorpusCallosum:ccKids:310");
   strcpy(sstr[24],"::CorpusCallosum:ccKids:312");
   strcpy(sstr[25],"::CorpusCallosum:ccKids:323");
   strcpy(sstr[26],"::CorpusCallosum:ccKids:330");
   strcpy(sstr[27],"::CorpusCallosum:ccKids:385");
   strcpy(sstr[28],"::CorpusCallosum:ccKids:387");
   strcpy(sstr[29],"::CorpusCallosum:ccKids:415");
*/

/* DRUGS */
/* first 7 belong to group I; the next 6 belong to
group II and the last belongs to group III  */

   strcpy(dir_name,"::CorpusCallosum:ccDrugs:");
   
   strcpy(str[0],"::CorpusCallosum:ccDrugs:2560_2");
   strcpy(str[1],"::CorpusCallosum:ccDrugs:2563_2");
   strcpy(str[2],"::CorpusCallosum:ccDrugs:2564_2");
   strcpy(str[3],"::CorpusCallosum:ccDrugs:2569_2");
   strcpy(str[4],"::CorpusCallosum:ccDrugs:2570_2");
   strcpy(str[5],"::CorpusCallosum:ccDrugs:2573_2");
   strcpy(str[6],"::CorpusCallosum:ccDrugs:2574_2");
   
   strcpy(str[7],"::CorpusCallosum:ccDrugs:2561_2");
   strcpy(str[8],"::CorpusCallosum:ccDrugs:2562_2");
   strcpy(str[9],"::CorpusCallosum:ccDrugs:2565_2");
   strcpy(str[10],"::CorpusCallosum:ccDrugs:2566_2");
   strcpy(str[11],"::CorpusCallosum:ccDrugs:2567_2");
   strcpy(str[12],"::CorpusCallosum:ccDrugs:2568_2");
   strcpy(str[13],"::CorpusCallosum:ccDrugs:2572_2");
   
   strcpy(sstr[0],"::CorpusCallosum:ccDrugs:2560");
   strcpy(sstr[1],"::CorpusCallosum:ccDrugs:2563");
   strcpy(sstr[2],"::CorpusCallosum:ccDrugs:2564");
   strcpy(sstr[3],"::CorpusCallosum:ccDrugs:2569");
   strcpy(sstr[4],"::CorpusCallosum:ccDrugs:2570");
   strcpy(sstr[5],"::CorpusCallosum:ccDrugs:2573");
   strcpy(sstr[6],"::CorpusCallosum:ccDrugs:2574");
   
   strcpy(sstr[7],"::CorpusCallosum:ccDrugs:2561");
   strcpy(sstr[8],"::CorpusCallosum:ccDrugs:2562");
   strcpy(sstr[9],"::CorpusCallosum:ccDrugs:2565");
   strcpy(sstr[10],"::CorpusCallosum:ccDrugs:2566");
   strcpy(sstr[11],"::CorpusCallosum:ccDrugs:2567");
   strcpy(sstr[12],"::CorpusCallosum:ccDrugs:2568");
   strcpy(sstr[13],"::CorpusCallosum:ccDrugs:2572");
  
   
   for (m=0;m<sample_count;m++) 
   {
      printf("m = %d\n",m);
	  strcpy(inname,str[m]);
	  strcpy(outname, sstr[m]);
	  strcat(outname, "dist"); 
   
   /* READ IN THE IMAGE */
	
	  printf("Opening %s\n",inname);
      if ((infile = fopen(inname,"r+")) == NULL)
      {
	    printf("Cannot open %s\n",inname);
	    exit(1);
      }   

	  n = fread(pic,sizeof(unsigned char),ydim*xdim,infile);
	  printf("%d bytes in image read\n",n);
	  fclose(infile);
	  
      calc_dist(); 
      output_results();
   }

   printf("\n\n done \n");
   systime = time(NULL);
   puts(ctime(&systime));
}