ccmean_shape.c

Computing 2D Skeleton

/* Program:  create average shapes                          */

/*                                                          */
/*                                                          */
/* Inputs:  radii of images and an axis  template           */

/* output:  average shapes: group1 gray,  grp II black      */
/*                                                          */
/*                                                          */
/* Date:     June, 2003                                     */
/*                                                          */
/*                                                          */





#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>

/* Schizophrenics vs Normals: ccSchzSchz vs ccSchzNorms */
/*
#define caseno 1
#define grp1 35
#define grp2 48
*/

/* Females vs Males - Adults:  ccSchzNorms + ccNorms  */
/*
#define caseno 2
#define grp1 31
#define grp2 37
*/

/* Females vs Males - Kids: ccKids */
/*
#define caseno 3
#define grp1 15
#define grp2 15
*/

/* Adults vs Kids -Males: ccSchzNorms + ccKids */
/*
#define caseno 4
#define grp1 27
#define grp2 15
*/

/* Adults vs Kids -Females: ccSchzNorms + ccKids */

#define caseno 5
#define grp1 21
#define grp2 15


 

#define ydim 364  
#define xdim 784 


#define num_axis_samples 72 /* main axis */
#define num_spl_samples 6  /* splenium */
#define num_hk_samples 34 /* genu hook */

#define num_axis_pts 134 /* num of pts on main axis template */
#define num_spl_pts 11 /* num of pts on splenium axis template */
#define num_hk_pts 68 /* num of pts on hook axis template */

#define axis_len 592.71 /* arclength of main axis template */
#define spl_len 42.28 /* arclength of splenium axis template */
#define hk_len 275.42 /* arclength of hk axis template */

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

/* file names */
FILE *infile,*outfile;
char inname[128],outname[128];

/* Program variables */
unsigned char outpic[ydim][xdim];
unsigned char ave_pic[2][ydim][xdim];

int Iarray[num_axis_pts],Jarray[num_axis_pts];
int num_template_pts,num_samples;
float template_len;
float arc_distance[xdim*ydim];
float rel_rho1[grp1][num_axis_samples];
float rel_rho2[grp2][num_axis_samples];
float rho[num_axis_pts];
float tmp[num_axis_samples];
float sample_ave[2][num_axis_samples];

void input_rho(int,char *);
void input_template(int,char *); 
void create_ave_shapes();
void calc_arc_distance();
void average_rho();
void assign_rho(int);
void ave_shapes(int);
void superpose();
void mark_sample_pts();
void output_results();
void run();


   
  void create_ave_shapes()
  {
     int  k;
     
     calc_arc_distance();
     average_rho(); 
     for (k=0;k<2;k++)
     {
       assign_rho(k);
       ave_shapes(k);
     }
  }




  void calc_arc_distance()
  {
     int m,p,q,p1,q1;
     float x;
     
     arc_distance[0]=0;
     for (m=1;m<num_template_pts;m++)
     {
        p = Iarray[m-1];
        q = Jarray[m-1];
        p1 = Iarray[m];
        q1 = Jarray[m];
        x = (p-p1)*(p-p1)+(q-q1)*(q-q1);
        if (x!=0) x = sqrt((float)x);
        arc_distance[m]=arc_distance[m-1]+x;
     }
  }
  
 
  
  
  void average_rho()
  {
     int m,k;
     float x,y,rho1,rho2;

     x = axis_len;
     
     for (m=0;m<num_samples;m++)
       sample_ave[0][m]=sample_ave[1][m]=0;
     
     for (m=0;m<num_samples;m++)
     {
       rho1=rho2=0;
       
       for (k=0;k<grp1;k++)
       {
          y = x*rel_rho1[k][m];
          rho1 += y;
       }
       
       for (k=0;k<grp2;k++)
       {
          y = x*rel_rho2[k][m];
          rho2 += y;
       }  
       sample_ave[0][m] = rho1/grp1;
       sample_ave[1][m] = rho2/grp2;
     }
  }
  
  
  
  
  void assign_rho(int k)
  {
     int m,k1,k2;
     float step;
     float v1,v2,x;
     
     x = axis_len; /* conversion to float */
     step = x/(num_axis_samples-1);
 
     for (m=0;m<num_template_pts;m++) rho[m]=0;
     rho[0] = sample_ave[k][0];
     for (m=1;m<num_template_pts-1;m++)
     {
       k1 = x = arc_distance[m]/step;
       k2 = k1+1;
       v1 = sample_ave[k][k1];
       v2 = sample_ave[k][k2];
       rho[m] = v2*(x-k1)+v1*(k2-x); /* interpolation */
     }
  }
  
  
  
  void ave_shapes(int k)
  {
     int i,j,m,p,q;
     float d;
     float dist[201][201];
  
     /* template */
     for (i=0;i<201;i++) for (j=0;j<201;j++)
     {
       d = (i-100)*(i-100)+(j-100)*(j-100);
       dist[i][j]=sqrt((double)d);
     }
     
     for (m=0;m<num_template_pts;m++)
     {    
       i = Iarray[m];
       j = Jarray[m];
       d = rho[m];
/* if (k==0) printf("(%d %d) %f\n",i,j,d); */
       if (d>=1)
       for (p=i-d;p<i+d+1;p++) for (q=j-d;q<j+d+1;q++)
       if (dist[p-i+100][q-j+100]<=d)
         ave_pic[k][p][q]=0;
     }
  }
  
  
  
  
  void superpose()
  {
    int i,j,k,clr[2];
    
	clr[0]=128; clr[1]=1;
	for (k=0;k<2;k++)
	{
	  for (i=0;i<ydim-1;i++) for (j=0;j<xdim-1;j++)
      {
        if (ave_pic[k][i][j]!=ave_pic[k][i+1][j])
        outpic[i][j] = outpic[i+1][j] = clr[k];
        if (ave_pic[k][i][j]!=ave_pic[k][i][j+1])
        outpic[i][j] = outpic[i][j+1] = clr[k];
      }
      for (i=0;i<ydim-1;i++)
      {
        j=xdim-1;
        if (ave_pic[k][i][j]!=ave_pic[k][i+1][j])
        outpic[i][j] = outpic[i+1][j] = clr[k];
      }
      for (j=0;j<xdim-1;j++)
      {
        i=ydim-1;
        if (ave_pic[k][i][j]!=ave_pic[k][i][j+1])
        outpic[i][j] = outpic[i][j+1] = clr[k];
      }
    }
  }
  
  
  
  
  void mark_sample_pts()
  {
     int m,i,j,i1,j1,i2,j2,k1,k2;
     int p,q;
     float d,d1,d2,len,step;
     float dist[201][201],x;
     
     /* template */
     for (i=0;i<201;i++) for (j=0;j<201;j++)
     {
       d = (i-100)*(i-100)+(j-100)*(j-100);
       dist[i][j]=sqrt((double)d);
     }  
     
     /* mark the main axis */
     x = len = axis_len; /* conversion to float */
     step = x/(num_axis_samples-1);
     
     k2 = 0;
     
     outpic[Iarray[0]][Jarray[0]] = 0;

     for (m=1;m<num_axis_samples-1;m++)
     if ((d=m*step)<=len)
     {
       while (arc_distance[k2]<d)
       {
         if (k2<num_axis_pts-1) k2++;
         else break;
       }
       k1 = k2-1;
       i1 = Iarray[k1]; j1 = Jarray[k1];
       i2 = Iarray[k2]; j2 = Jarray[k2];
       d1 = arc_distance[k1];
       d2 = arc_distance[k2];
       i = (i2*(d-d1)+i1*(d2-d))/(d2-d1);
       j = (j2*(d-d1)+j1*(d2-d))/(d2-d1);
       if (m%4==3)
       for (p=i-1;p<i+2;p++) for (q=j-1;q<j+2;q++) outpic[p][q]=0;
       else outpic[i][j]=0;
     }
     m = num_template_pts - 1;
     i = Iarray[m]; j = Jarray[m];
     for (p=i-1;p<i+2;p++) for (q=j-1;q<j+2;q++)
            outpic[p][q]=0;;     
  }
  
 
  

 void output_results()
 { 
     strcpy(outname,"::CorpusCallosum:Classification:Results:");
     if (caseno==1) strcat(outname,"ave_shapeNS");
     else if (caseno==2) strcat(outname,"ave_shapeFMadults");
     else if (caseno==3) strcat(outname,"ave_shapeFMkids");
     else if (caseno==4) strcat(outname,"ave_shapeAKmales"); 
     else if (caseno==5) strcat(outname,"ave_shapeAKfemales"); 
     
     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); 
 }
 
 
 
 void input_rho(num_pts,string)
 int num_pts;
 char *string;
 {
    int i,j;
    float x;
    
    if (caseno==1)
 /* Schizophrenics vs Normals: ccSchzSchz vs ccSchzNorms */
    {
      strcpy(inname, "::CorpusCallosum:ccSchzSchz:");
      strcat(inname, string);
      printf("Opening %s\n",inname);
      if ((infile = fopen(inname,"r+")) == NULL)
      {
	    printf("Cannot open %s\n",inname);
	    exit(1);
      }
      for (i=0;i<grp1;i++) for (j=0;j<num_pts;j++)
       fscanf(infile,"%f",rel_rho1[i]+j);
      fclose(infile);
   
      strcpy(inname, "::CorpusCallosum:ccSchzNorms:");
      strcat(inname, string);
      printf("Opening %s\n",inname);
      if ((infile = fopen(inname,"r+")) == NULL)
      {
	    printf("Cannot open %s\n",inname);
	    exit(1);
      }
      for (i=0;i<grp2;i++) for (j=0;j<num_pts;j++)
       fscanf(infile,"%f",rel_rho2[i]+j);
      fclose(infile);
    }
    
    else  if (caseno==2)
 /* Females vs Males - Adults:  ccSchzNorms + ccNorms  */
    {
      strcpy(inname, "::CorpusCallosum:ccSchzNorms:");
      strcat(inname, string);
      printf("Opening %s\n",inname);
      if ((infile = fopen(inname,"r+")) == NULL)
      {
	    printf("Cannot open %s\n",inname);
	    exit(1);
      }
      for (i=0;i<21;i++) for (j=0;j<num_pts;j++)
       fscanf(infile,"%f",rel_rho1[i]+j);
   
      for (i=0;i<27;i++) for (j=0;j<num_pts;j++)
       fscanf(infile,"%f",rel_rho2[i]+j);
      fclose(infile);
   
      strcpy(inname, "::CorpusCallosum:ccNorms:");
      strcat(inname, string);
      printf("Opening %s\n",inname);
      if ((infile = fopen(inname,"r+")) == NULL)
      {
	    printf("Cannot open %s\n",inname);
	    exit(1);
      }
      for (i=0;i<10;i++) for (j=0;j<num_pts;j++)
       fscanf(infile,"%f",rel_rho1[i+21]+j);
      for (i=0;i<10;i++) for (j=0;j<num_pts;j++)
       fscanf(infile,"%f",rel_rho2[i+27]+j);
      fclose(infile);
    }
    
/* Females vs Males - Kids:  ccSchzNorms + ccNorms  */    
    else if (caseno==3)
    {
      strcpy(inname, "::CorpusCallosum:ccKids:");
      strcat(inname, string);
      printf("Opening %s\n",inname);
      if ((infile = fopen(inname,"r+")) == NULL)
      {
	    printf("Cannot open %s\n",inname);
	    exit(1);
      }
   
      for (i=0;i<grp1;i++) for (j=0;j<num_pts;j++)
       fscanf(infile,"%f",rel_rho1[i]+j);

      for (i=0;i<grp2;i++) for (j=0;j<num_pts;j++)
       fscanf(infile,"%f",rel_rho2[i]+j);
      fclose(infile);
    }
    
    /* Adults vs Kids - Males:  ccSchzNorms + ccKids  */    
    else if (caseno==4)
    {
      strcpy(inname, "::CorpusCallosum:ccSchzNorms:");
      strcat(inname, string);
      printf("Opening %s\n",inname);
      if ((infile = fopen(inname,"r+")) == NULL)
      {
	    printf("Cannot open %s\n",inname);
	    exit(1);
      }
      /* skip first 21 records */
      for (i=0;i<21;i++) for (j=0;j<num_pts;j++)
        fscanf(infile,"%f",&x);
   
      for (i=0;i<27;i++) for (j=0;j<num_pts;j++)
        fscanf(infile,"%f",rel_rho1[i]+j);
      fclose(infile);
      
      strcpy(inname, "::CorpusCallosum:ccKids:");
      strcat(inname, string);
      printf("Opening %s\n",inname);
      if ((infile = fopen(inname,"r+")) == NULL)
      {
	    printf("Cannot open %s\n",inname);
	    exit(1);
      }
      /* skip first 15 records */
      for (i=0;i<15;i++) for (j=0;j<num_pts;j++)
        fscanf(infile,"%f",&x);
       
      for (i=0;i<15;i++) for (j=0;j<num_pts;j++)
        fscanf(infile,"%f",rel_rho2[i]+j);
      fclose(infile);
    }
    
        /* Adults vs Kids - Females:  ccSchzNorms + ccKids  */    
    else if (caseno==5)
    {
      strcpy(inname, "::CorpusCallosum:ccSchzNorms:");
      strcat(inname, string);
      printf("Opening %s\n",inname);
      if ((infile = fopen(inname,"r+")) == NULL)
      {
	    printf("Cannot open %s\n",inname);
	    exit(1);
      }
      for (i=0;i<grp1;i++) for (j=0;j<num_pts;j++)
        fscanf(infile,"%f",rel_rho1[i]+j);
      fclose(infile);
      
      strcpy(inname, "::CorpusCallosum:ccKids:");
      strcat(inname, string);
      printf("Opening %s\n",inname);
      if ((infile = fopen(inname,"r+")) == NULL)
      {
	    printf("Cannot open %s\n",inname);
	    exit(1);
      }
      for (i=0;i<grp2;i++) for (j=0;j<num_pts;j++)
        fscanf(infile,"%f",rel_rho2[i]+j);
      fclose(infile);
    }
 }    
 
 
 
 void input_template(num_pts,string)
 int num_pts;
 char *string;
 {
    int j;
    
    strcpy(inname, "::CorpusCallosum:");
    strcat(inname, string);
    printf("Opening %s\n",inname);
    if ((infile = fopen(inname,"r+")) == NULL)
    {
	  printf("Cannot open %s\n",inname);
	  exit(1);
    }   
    for (j=0;j<num_pts;j++)
    fscanf(infile,"%d %d",Iarray+j,Jarray+j);
    fclose(infile);  
    
 }
 
 

void run(void)
{
   int i,j;
   time_t systime;
   	   
   SIOUXSettings.asktosaveonclose = 0;
   /* SIOUXSettings.showstatusline = 1; */
	      systime = time(NULL);
	      puts(ctime(&systime));

   strcpy(inname,"::CorpusCallosum:template_axis_pic");
   if ((infile = fopen(inname,"r+")) == NULL)
   {
	 printf("Cannot open %s\n",inname);
	 exit(1);
   }
   fread(outpic,sizeof(unsigned char),ydim*xdim,infile);
   fclose(infile);
   
   for (i=0;i<ydim;i++) for (j=0;j<xdim;j++)
     ave_pic[0][i][j]=ave_pic[1][i][j]=255;
 
   num_template_pts = num_axis_pts;
   template_len = axis_len;
   input_template(num_template_pts,"main_axis_template");
   num_samples = num_axis_samples;
   input_rho(num_samples,"axis_rho");
   create_ave_shapes();
/*   mark_sample_pts(); */
   
   num_template_pts = num_spl_pts;
   template_len = spl_len;
   input_template(num_template_pts,"spl_axis_template");
   num_samples = num_spl_samples;
   input_rho(num_samples,"spl_rho");
   create_ave_shapes();
   
   num_template_pts = num_hk_pts;
   template_len = hk_len;
   input_template(num_template_pts,"hk_axis_template");
   num_samples = num_hk_samples;
   input_rho(num_samples,"hk_rho");
   create_ave_shapes();
   
   superpose(); 
   
   output_results(); 

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