Skip to content
Snippets Groups Projects
Commit f3bd8377 authored by Mark Jenkinson's avatar Mark Jenkinson
Browse files

Initial revision

parents
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 887 additions and 0 deletions
siena_diff.c 0 → 100644
+ 887
− 0
View file @ f3bd8377
/* {{{ Copyright etc. */
/* siena_diff - compute brain change using edge motion or segmentation
Stephen Smith, FMRIB Image Analysis Group
Copyright (C) 1999-2002 University of Oxford */
/* CCOPYRIGHT */
/* }}} */
/* {{{ defines, includes and typedefs */
#include "libss/libss.h"
#include "libss/libavw.h"
#include <unistd.h>
/* }}} */
/* {{{ usage */
void usage()
{
printf("\nUsage: siena_diff <input1_basename> <input2_basename> [options] [-s segmentation options]\n\n");
printf("[-d] debug - generate edge images and don't remove temporary images\n");
printf("[-2] don't segment grey+white separately (because there is poor grey-white contrast)\n\n");
printf("[-c <corr>] apply self-calibrating correction factor\n");
printf("[-e] erode joint mask a lot instead of dilating it slightly (ie find ventricle surface)\n");
printf("[-i] ignore flow in z (may be good if top of brain is missing)\n");
printf("[-m] apply <input1_basename>_talmask to brain edge points\n");
/* printf("[-t <n>] ignore top n slices (may be good if top of brain is missing)\n");*/
/* printf("[-b <n>] ignore bottom n slices (may be good if top of brain is missing)\n");*/
printf("[-s <options>] <options> to be passed to segmentation (type \"fast\" to get these)\n\n");
exit(1);
}
/* }}} */
/* {{{ main(argc, argv) */
#define CORRWIDTH 3
#define SEARCH 4
int main(argc, argv)
int argc;
char *argv [];
{
/* {{{ vars */
unsigned char *mask, *mask_dil;
char filename[10000], thestring[10000], segoptions[10000], fsldir[10000];
int x_size, y_size, z_size, size, x, y, z, i, count,
seg2=0, ignore_z=0, erode_mask=0, ignore_top_slices=0,
ignore_bottom_slices=0, debug=0, flow_output=1, edge_masking=0;
float *flow, *in1, *in2, tmpf, calib=1.0, *m1, *m2, *seg1, *arrA, *arrB, *arr1, *arr2;
double total, voxel_volume, voxel_area, ex, ey, ez, lowest;
image_struct im;
/* }}} */
/* {{{ process arguments */
if (argc<3)
usage();
sprintf(fsldir,"%s",getenv("FSLDIR"));
for (i = 3; i < argc; i++) {
if (!strcmp(argv[i], "-i"))
ignore_z=1;
else if (!strcmp(argv[i], "-e"))
erode_mask=1;
else if (!strcmp(argv[i], "-d"))
debug=1;
else if (!strcmp(argv[i], "-2"))
seg2=1;
else if (!strcmp(argv[i], "-c"))
/* {{{ apply self-calibrating factor */
{
i++;
if (argc<i+1)
{
printf("Error: no factor given following -c\n");
usage();
}
calib=atof(argv[i]);
}
/* }}} */
else if (!strcmp(argv[i], "-m"))
edge_masking=1;
else if (!strcmp(argv[i], "-t"))
/* {{{ ignore n slices at top */
{
i++;
if (argc<i+1)
{
printf("Error: no number of slices given following -t\n");
usage();
}
ignore_top_slices=atoi(argv[i]);
}
/* }}} */
else if (!strcmp(argv[i], "-b"))
/* {{{ ignore n slices at bottom */
{
i++;
if (argc<i+1)
{
printf("Error: no number of slices given following -b\n");
usage();
}
ignore_bottom_slices=atoi(argv[i]);
}
/* }}} */
else if (!strcmp(argv[i], "-s"))
/* {{{ segmentation options */
{
i++;
segoptions[0]=0;
while(i<argc)
{
strcat(segoptions,argv[i]);
strcat(segoptions," ");
i++;
}
}
/* }}} */
else
usage();
}
/* }}} */
/* {{{ transform images and masks */
sprintf(thestring,"%s/bin/flirt -o %s_halfwayto_%s -applyisoxfm 1 -paddingsize 0 -init %s_halfwayto_%s.mat -ref %s -in %s",
fsldir,argv[1],argv[2],argv[1],argv[2],argv[1],argv[1]);
printf("%s\n",thestring); system(thestring);
sprintf(thestring,"%s/bin/flirt -o %s_halfwayto_%s -applyisoxfm 1 -paddingsize 0 -init %s_halfwayto_%s.mat -ref %s -in %s",
fsldir,argv[2],argv[1],argv[2],argv[1],argv[1],argv[2]);
printf("%s\n",thestring); system(thestring);
sprintf(thestring,"%s/bin/flirt -o %s_halfwayto_%s_mask -applyisoxfm 1 -paddingsize 0 -init %s_halfwayto_%s.mat -ref %s -in %s_brain_mask",
fsldir,argv[1],argv[2],argv[1],argv[2],argv[1],argv[1]);
printf("%s\n",thestring); system(thestring);
sprintf(thestring,"%s/bin/flirt -o %s_halfwayto_%s_mask -applyisoxfm 1 -paddingsize 0 -init %s_halfwayto_%s.mat -ref %s -in %s_brain_mask",
fsldir,argv[2],argv[1],argv[2],argv[1],argv[1],argv[2]);
printf("%s\n",thestring); system(thestring);
if (edge_masking)
{
sprintf(thestring,"%s/bin/flirt -o %s_halfwayto_%s_valid_mask -applyisoxfm 1 -paddingsize 0 -init %s_halfwayto_%s.mat -ref %s -in %s_valid_mask_with_%s",
fsldir,argv[1],argv[2],argv[1],argv[2],argv[1],argv[1],argv[2]);
printf("%s\n",thestring); system(thestring);
}
/* }}} */
/* {{{ dilate masks, read transformed images and masks, and combine to jointly-masked transformed images */
printf("reading and combining transformed masks\n");
sprintf(filename,"%s_halfwayto_%s_mask",argv[1],argv[2]);
avw_read(filename,&im);
m1=im.i;
/* {{{ process header info */
x_size=im.x;
y_size=im.y;
z_size=im.z;
size=x_size*y_size*z_size;
voxel_volume = ABS( im.xv * im.yv * im.zv );
voxel_area = pow(voxel_volume,((double)0.6666667));
printf("final image dimensions = %d %d %d, voxel volume = %.3fmm^3, voxel area = %.3fmm^2\n",
x_size,y_size,z_size,voxel_volume,voxel_area);
/* }}} */
sprintf(filename,"%s_halfwayto_%s_mask",argv[2],argv[1]);
avw_read(filename,&im);
m2=im.i;
mask = (unsigned char*) malloc(sizeof(unsigned char)*x_size*y_size*z_size);
for(i=0;i<size;i++)
if (m1[i]+m2[i]>0.5) /* mask1 OR mask2 */
mask[i]=1;
else
mask[i]=0;
free(m1);
free(m2);
printf("dilating/eroding combined mask\n");
mask_dil = (unsigned char*) malloc(sizeof(unsigned char)*x_size*y_size*z_size);
memset(mask_dil,(unsigned char)0,sizeof(unsigned char)*x_size*y_size*z_size);
if (erode_mask)
/* {{{ erode mask LOTS */
/* note - this is a hack - depends on the scale of the images! */
{
/* {{{ mask -> mask_dil */
for (z=0; z<z_size; z++)
for (y=0; y<y_size; y++)
for (x=0; x<x_size; x++)
{
if ( (IA(mask,x,y,z)==1) &&
((x>0)&&(IA(mask,x-1,y,z)==1)) && ((x<x_size-1)&&(IA(mask,x+1,y,z)==1)) &&
((y>0)&&(IA(mask,x,y-1,z)==1)) && ((y<y_size-1)&&(IA(mask,x,y+1,z)==1)) &&
((z>0)&&(IA(mask,x,y,z-1)==1)) && ((z<z_size-1)&&(IA(mask,x,y,z+1)==1)) )
IA(mask_dil,x,y,z)=1;
else
IA(mask_dil,x,y,z)=0;
}
/* }}} */
/* {{{ mask_dil -> mask */
for (z=0; z<z_size; z++)
for (y=0; y<y_size; y++)
for (x=0; x<x_size; x++)
{
if ( (IA(mask_dil,x,y,z)==1) &&
((x>0)&&(IA(mask_dil,x-1,y,z)==1)) && ((x<x_size-1)&&(IA(mask_dil,x+1,y,z)==1)) &&
((y>0)&&(IA(mask_dil,x,y-1,z)==1)) && ((y<y_size-1)&&(IA(mask_dil,x,y+1,z)==1)) &&
((z>0)&&(IA(mask_dil,x,y,z-1)==1)) && ((z<z_size-1)&&(IA(mask_dil,x,y,z+1)==1)) )
IA(mask,x,y,z)=1;
else
IA(mask,x,y,z)=0;
}
/* }}} */
/* {{{ mask -> mask_dil */
for (z=0; z<z_size; z++)
for (y=0; y<y_size; y++)
for (x=0; x<x_size; x++)
{
if ( (IA(mask,x,y,z)==1) &&
((x>0)&&(IA(mask,x-1,y,z)==1)) && ((x<x_size-1)&&(IA(mask,x+1,y,z)==1)) &&
((y>0)&&(IA(mask,x,y-1,z)==1)) && ((y<y_size-1)&&(IA(mask,x,y+1,z)==1)) &&
((z>0)&&(IA(mask,x,y,z-1)==1)) && ((z<z_size-1)&&(IA(mask,x,y,z+1)==1)) )
IA(mask_dil,x,y,z)=1;
else
IA(mask_dil,x,y,z)=0;
}
/* }}} */
/* {{{ mask_dil -> mask */
for (z=0; z<z_size; z++)
for (y=0; y<y_size; y++)
for (x=0; x<x_size; x++)
{
if ( (IA(mask_dil,x,y,z)==1) &&
((x>0)&&(IA(mask_dil,x-1,y,z)==1)) && ((x<x_size-1)&&(IA(mask_dil,x+1,y,z)==1)) &&
((y>0)&&(IA(mask_dil,x,y-1,z)==1)) && ((y<y_size-1)&&(IA(mask_dil,x,y+1,z)==1)) &&
((z>0)&&(IA(mask_dil,x,y,z-1)==1)) && ((z<z_size-1)&&(IA(mask_dil,x,y,z+1)==1)) )
IA(mask,x,y,z)=1;
else
IA(mask,x,y,z)=0;
}
/* }}} */
/* {{{ mask -> mask_dil */
for (z=0; z<z_size; z++)
for (y=0; y<y_size; y++)
for (x=0; x<x_size; x++)
{
if ( (IA(mask,x,y,z)==1) &&
((x>0)&&(IA(mask,x-1,y,z)==1)) && ((x<x_size-1)&&(IA(mask,x+1,y,z)==1)) &&
((y>0)&&(IA(mask,x,y-1,z)==1)) && ((y<y_size-1)&&(IA(mask,x,y+1,z)==1)) &&
((z>0)&&(IA(mask,x,y,z-1)==1)) && ((z<z_size-1)&&(IA(mask,x,y,z+1)==1)) )
IA(mask_dil,x,y,z)=1;
else
IA(mask_dil,x,y,z)=0;
}
/* }}} */
/* {{{ mask_dil -> mask */
for (z=0; z<z_size; z++)
for (y=0; y<y_size; y++)
for (x=0; x<x_size; x++)
{
if ( (IA(mask_dil,x,y,z)==1) &&
((x>0)&&(IA(mask_dil,x-1,y,z)==1)) && ((x<x_size-1)&&(IA(mask_dil,x+1,y,z)==1)) &&
((y>0)&&(IA(mask_dil,x,y-1,z)==1)) && ((y<y_size-1)&&(IA(mask_dil,x,y+1,z)==1)) &&
((z>0)&&(IA(mask_dil,x,y,z-1)==1)) && ((z<z_size-1)&&(IA(mask_dil,x,y,z+1)==1)) )
IA(mask,x,y,z)=1;
else
IA(mask,x,y,z)=0;
}
/* }}} */
/* {{{ mask -> mask_dil */
for (z=0; z<z_size; z++)
for (y=0; y<y_size; y++)
for (x=0; x<x_size; x++)
{
if ( (IA(mask,x,y,z)==1) &&
((x>0)&&(IA(mask,x-1,y,z)==1)) && ((x<x_size-1)&&(IA(mask,x+1,y,z)==1)) &&
((y>0)&&(IA(mask,x,y-1,z)==1)) && ((y<y_size-1)&&(IA(mask,x,y+1,z)==1)) &&
((z>0)&&(IA(mask,x,y,z-1)==1)) && ((z<z_size-1)&&(IA(mask,x,y,z+1)==1)) )
IA(mask_dil,x,y,z)=1;
else
IA(mask_dil,x,y,z)=0;
}
/* }}} */
/* {{{ mask_dil -> mask */
for (z=0; z<z_size; z++)
for (y=0; y<y_size; y++)
for (x=0; x<x_size; x++)
{
if ( (IA(mask_dil,x,y,z)==1) &&
((x>0)&&(IA(mask_dil,x-1,y,z)==1)) && ((x<x_size-1)&&(IA(mask_dil,x+1,y,z)==1)) &&
((y>0)&&(IA(mask_dil,x,y-1,z)==1)) && ((y<y_size-1)&&(IA(mask_dil,x,y+1,z)==1)) &&
((z>0)&&(IA(mask_dil,x,y,z-1)==1)) && ((z<z_size-1)&&(IA(mask_dil,x,y,z+1)==1)) )
IA(mask,x,y,z)=1;
else
IA(mask,x,y,z)=0;
}
/* }}} */
/* {{{ mask -> mask_dil */
for (z=0; z<z_size; z++)
for (y=0; y<y_size; y++)
for (x=0; x<x_size; x++)
{
if ( (IA(mask,x,y,z)==1) &&
((x>0)&&(IA(mask,x-1,y,z)==1)) && ((x<x_size-1)&&(IA(mask,x+1,y,z)==1)) &&
((y>0)&&(IA(mask,x,y-1,z)==1)) && ((y<y_size-1)&&(IA(mask,x,y+1,z)==1)) &&
((z>0)&&(IA(mask,x,y,z-1)==1)) && ((z<z_size-1)&&(IA(mask,x,y,z+1)==1)) )
IA(mask_dil,x,y,z)=1;
else
IA(mask_dil,x,y,z)=0;
}
/* }}} */
/* {{{ mask_dil -> mask */
for (z=0; z<z_size; z++)
for (y=0; y<y_size; y++)
for (x=0; x<x_size; x++)
{
if ( (IA(mask_dil,x,y,z)==1) &&
((x>0)&&(IA(mask_dil,x-1,y,z)==1)) && ((x<x_size-1)&&(IA(mask_dil,x+1,y,z)==1)) &&
((y>0)&&(IA(mask_dil,x,y-1,z)==1)) && ((y<y_size-1)&&(IA(mask_dil,x,y+1,z)==1)) &&
((z>0)&&(IA(mask_dil,x,y,z-1)==1)) && ((z<z_size-1)&&(IA(mask_dil,x,y,z+1)==1)) )
IA(mask,x,y,z)=1;
else
IA(mask,x,y,z)=0;
}
/* }}} */
/* {{{ mask -> mask_dil */
for (z=0; z<z_size; z++)
for (y=0; y<y_size; y++)
for (x=0; x<x_size; x++)
{
if ( (IA(mask,x,y,z)==1) &&
((x>0)&&(IA(mask,x-1,y,z)==1)) && ((x<x_size-1)&&(IA(mask,x+1,y,z)==1)) &&
((y>0)&&(IA(mask,x,y-1,z)==1)) && ((y<y_size-1)&&(IA(mask,x,y+1,z)==1)) &&
((z>0)&&(IA(mask,x,y,z-1)==1)) && ((z<z_size-1)&&(IA(mask,x,y,z+1)==1)) )
IA(mask_dil,x,y,z)=1;
else
IA(mask_dil,x,y,z)=0;
}
/* }}} */
/* {{{ mask_dil -> mask */
for (z=0; z<z_size; z++)
for (y=0; y<y_size; y++)
for (x=0; x<x_size; x++)
{
if ( (IA(mask_dil,x,y,z)==1) &&
((x>0)&&(IA(mask_dil,x-1,y,z)==1)) && ((x<x_size-1)&&(IA(mask_dil,x+1,y,z)==1)) &&
((y>0)&&(IA(mask_dil,x,y-1,z)==1)) && ((y<y_size-1)&&(IA(mask_dil,x,y+1,z)==1)) &&
((z>0)&&(IA(mask_dil,x,y,z-1)==1)) && ((z<z_size-1)&&(IA(mask_dil,x,y,z+1)==1)) )
IA(mask,x,y,z)=1;
else
IA(mask,x,y,z)=0;
}
/* }}} */
/* {{{ mask -> mask_dil */
for (z=0; z<z_size; z++)
for (y=0; y<y_size; y++)
for (x=0; x<x_size; x++)
{
if ( (IA(mask,x,y,z)==1) &&
((x>0)&&(IA(mask,x-1,y,z)==1)) && ((x<x_size-1)&&(IA(mask,x+1,y,z)==1)) &&
((y>0)&&(IA(mask,x,y-1,z)==1)) && ((y<y_size-1)&&(IA(mask,x,y+1,z)==1)) &&
((z>0)&&(IA(mask,x,y,z-1)==1)) && ((z<z_size-1)&&(IA(mask,x,y,z+1)==1)) )
IA(mask_dil,x,y,z)=1;
else
IA(mask_dil,x,y,z)=0;
}
/* }}} */
/* {{{ mask_dil -> mask */
for (z=0; z<z_size; z++)
for (y=0; y<y_size; y++)
for (x=0; x<x_size; x++)
{
if ( (IA(mask_dil,x,y,z)==1) &&
((x>0)&&(IA(mask_dil,x-1,y,z)==1)) && ((x<x_size-1)&&(IA(mask_dil,x+1,y,z)==1)) &&
((y>0)&&(IA(mask_dil,x,y-1,z)==1)) && ((y<y_size-1)&&(IA(mask_dil,x,y+1,z)==1)) &&
((z>0)&&(IA(mask_dil,x,y,z-1)==1)) && ((z<z_size-1)&&(IA(mask_dil,x,y,z+1)==1)) )
IA(mask,x,y,z)=1;
else
IA(mask,x,y,z)=0;
}
/* }}} */
/* {{{ mask -> mask_dil */
for (z=0; z<z_size; z++)
for (y=0; y<y_size; y++)
for (x=0; x<x_size; x++)
{
if ( (IA(mask,x,y,z)==1) &&
((x>0)&&(IA(mask,x-1,y,z)==1)) && ((x<x_size-1)&&(IA(mask,x+1,y,z)==1)) &&
((y>0)&&(IA(mask,x,y-1,z)==1)) && ((y<y_size-1)&&(IA(mask,x,y+1,z)==1)) &&
((z>0)&&(IA(mask,x,y,z-1)==1)) && ((z<z_size-1)&&(IA(mask,x,y,z+1)==1)) )
IA(mask_dil,x,y,z)=1;
else
IA(mask_dil,x,y,z)=0;
}
/* }}} */
/* {{{ mask_dil -> mask */
for (z=0; z<z_size; z++)
for (y=0; y<y_size; y++)
for (x=0; x<x_size; x++)
{
if ( (IA(mask_dil,x,y,z)==1) &&
((x>0)&&(IA(mask_dil,x-1,y,z)==1)) && ((x<x_size-1)&&(IA(mask_dil,x+1,y,z)==1)) &&
((y>0)&&(IA(mask_dil,x,y-1,z)==1)) && ((y<y_size-1)&&(IA(mask_dil,x,y+1,z)==1)) &&
((z>0)&&(IA(mask_dil,x,y,z-1)==1)) && ((z<z_size-1)&&(IA(mask_dil,x,y,z+1)==1)) )
IA(mask,x,y,z)=1;
else
IA(mask,x,y,z)=0;
}
/* }}} */
/* {{{ mask -> mask_dil */
for (z=0; z<z_size; z++)
for (y=0; y<y_size; y++)
for (x=0; x<x_size; x++)
{
if ( (IA(mask,x,y,z)==1) &&
((x>0)&&(IA(mask,x-1,y,z)==1)) && ((x<x_size-1)&&(IA(mask,x+1,y,z)==1)) &&
((y>0)&&(IA(mask,x,y-1,z)==1)) && ((y<y_size-1)&&(IA(mask,x,y+1,z)==1)) &&
((z>0)&&(IA(mask,x,y,z-1)==1)) && ((z<z_size-1)&&(IA(mask,x,y,z+1)==1)) )
IA(mask_dil,x,y,z)=1;
else
IA(mask_dil,x,y,z)=0;
}
/* }}} */
}
/* }}} */
else
/* {{{ dilate mask by one */
{
for (z=0; z<z_size; z++)
for (y=0; y<y_size; y++)
for (x=0; x<x_size; x++)
{
if ( (IA(mask,x,y,z)==1) ||
((x>0)&&(IA(mask,x-1,y,z)==1)) || ((x<x_size-1)&&(IA(mask,x+1,y,z)==1)) ||
((y>0)&&(IA(mask,x,y-1,z)==1)) || ((y<y_size-1)&&(IA(mask,x,y+1,z)==1)) ||
((z>0)&&(IA(mask,x,y,z-1)==1)) || ((z<z_size-1)&&(IA(mask,x,y,z+1)==1)) )
IA(mask_dil,x,y,z)=1;
}
}
/* }}} */
free(mask);
printf("reading transformed images and applying mask\n");
sprintf(filename,"%s_halfwayto_%s",argv[1],argv[2]);
avw_read(filename,&im);
in1=im.i;
lowest=in1[0];
for(i=0;i<size;i++)
if (in1[i]<lowest) lowest=in1[i];
for(i=0;i<size;i++)
in1[i]=(in1[i]-lowest)*mask_dil[i];
free(mask_dil);
/* }}} */
/* {{{ mallocs, recentre 1D arrays */
if ((arrA = calloc(10*(CORRWIDTH+SEARCH),sizeof(float)))==NULL)
{
printf("aaargh\n");
exit(1);
}
if ((arrB = calloc(10*(CORRWIDTH+SEARCH),sizeof(float)))==NULL)
{
printf("aaargh\n");
exit(1);
}
if ((arr1 = calloc(10*(CORRWIDTH+SEARCH),sizeof(float)))==NULL)
{
printf("aaargh\n");
exit(1);
}
if ((arr2 = calloc(10*(CORRWIDTH+SEARCH),sizeof(float)))==NULL)
{
printf("aaargh\n");
exit(1);
}
arrA = arrA + 5*(CORRWIDTH+SEARCH);
arrB = arrB + 5*(CORRWIDTH+SEARCH);
arr1 = arr1 + 5*(CORRWIDTH+SEARCH);
arr2 = arr2 + 5*(CORRWIDTH+SEARCH);
/* }}} */
/* {{{ do segmentation on image 1 */
{
FILE *tmpfp;
sprintf(thestring,"%s_halfwayto_%s_brain_seg.hdr",argv[1],argv[2]);
/* if((tmpfp=fopen(thestring,"rb"))==NULL)*/ /* test for existing segmentation output */
if(1)
{
char segtype[100];
if (seg2) sprintf(segtype,"-c 2"); else segtype[0]=0;
printf("saving image 1 to disk prior to segmentation\n");
im.i=in1;
sprintf(thestring,"%s_halfwayto_%s_brain",argv[1],argv[2]);
avw_write(thestring,im);
free(in1);
sprintf(thestring,"%s/bin/fast %s %s %s_halfwayto_%s_brain > %s_halfwayto_%s_brain.vol 2>&1",
fsldir,segtype,segoptions,argv[1],argv[2],argv[1],argv[2]);
printf("%s\n",thestring);
system(thestring);
}
else
{
printf("using previously carried out segmentation\n");
free(in1);
}
}
/* }}} */
/* {{{ read segmentation output into edges1 and simplify; reread in1 and in2 */
printf("finding brain edges\n");
sprintf(filename,"%s_halfwayto_%s_brain_seg",argv[1],argv[2]);
avw_read(filename,&im);
seg1=im.i;
for(x=0;x<size;x++)
if (seg1[x]>1)
seg1[x]=1;
else
seg1[x]=0;
if (edge_masking)
{
sprintf(filename,"%s_halfwayto_%s_valid_mask",argv[1],argv[2]);
avw_read(filename,&im);
m1=im.i;
}
sprintf(filename,"%s_halfwayto_%s",argv[1],argv[2]);
avw_read(filename,&im);
in1=im.i;
sprintf(filename,"%s_halfwayto_%s",argv[2],argv[1]);
avw_read(filename,&im);
in2=im.i;
if ((flow = calloc(x_size*y_size*z_size,sizeof(float)))==NULL)
{
printf("aaargh\n");
exit(1);
}
/* }}} */
/* {{{ find segmentation-based edges in image 1 and flow */
printf("finding flow\n");
count=0;
total=0;
ignore_bottom_slices=MAX(1,ignore_bottom_slices);
ignore_top_slices=MAX(1,ignore_top_slices);
for (z=ignore_bottom_slices; z<z_size-ignore_top_slices; z++)
for (y=1; y<y_size-1; y++)
for (x=1; x<x_size-1; x++)
{
if ( (IA(seg1,x,y,z)>0.5) && /* not background or CSF */
( (IA(seg1,x+1,y,z)<0.5) || (IA(seg1,x-1,y,z)<0.5) ||
(IA(seg1,x,y+1,z)<0.5) || (IA(seg1,x,y-1,z)<0.5) ||
(IA(seg1,x,y,z+1)<0.5) || (IA(seg1,x,y,z-1)<0.5) ) &&
( ( ! edge_masking ) || ( IA(m1,x,y,z)>0 ) ) )
{
int pos, neg, r, rr, rrr, d, X, Y, Z;
float ss, maxss, segvalpos=0, segvalneg=0;
image_struct im1=im, im2=im;
/* {{{ find local gradient and derive unit normal */
ex = ( 10*(IA(in1,x+1,y,z)-IA(in1,x-1,y,z)) +
5*(IA(in1,x+1,y+1,z)+IA(in1,x+1,y-1,z)+IA(in1,x+1,y,z+1)+IA(in1,x+1,y,z-1)-
IA(in1,x-1,y+1,z)-IA(in1,x-1,y-1,z)-IA(in1,x-1,y,z+1)-IA(in1,x-1,y,z-1)) +
2*(IA(in1,x+1,y+1,z+1)+IA(in1,x+1,y-1,z+1)+IA(in1,x+1,y+1,z-1)+IA(in1,x+1,y-1,z-1)-
IA(in1,x-1,y+1,z+1)-IA(in1,x-1,y-1,z+1)-IA(in1,x-1,y+1,z-1)-IA(in1,x-1,y-1,z-1)) ) / 38;
ey = ( 10*(IA(in1,x,y+1,z)-IA(in1,x,y-1,z)) +
5*(IA(in1,x+1,y+1,z)+IA(in1,x-1,y+1,z)+IA(in1,x,y+1,z+1)+IA(in1,x,y+1,z-1)-
IA(in1,x+1,y-1,z)-IA(in1,x-1,y-1,z)-IA(in1,x,y-1,z+1)-IA(in1,x,y-1,z-1)) +
2*(IA(in1,x+1,y+1,z+1)+IA(in1,x-1,y+1,z+1)+IA(in1,x+1,y+1,z-1)+IA(in1,x-1,y+1,z-1)-
IA(in1,x+1,y-1,z+1)-IA(in1,x-1,y-1,z+1)-IA(in1,x+1,y-1,z-1)-IA(in1,x-1,y-1,z-1)) ) / 38;
ez = ( 10*(IA(in1,x,y,z+1)-IA(in1,x,y,z-1)) +
5*(IA(in1,x,y+1,z+1)+IA(in1,x,y-1,z+1)+IA(in1,x+1,y,z+1)+IA(in1,x-1,y,z+1)-
IA(in1,x,y+1,z-1)-IA(in1,x,y-1,z-1)-IA(in1,x+1,y,z-1)-IA(in1,x-1,y,z-1)) +
2*(IA(in1,x+1,y+1,z+1)+IA(in1,x+1,y-1,z+1)+IA(in1,x-1,y+1,z+1)+IA(in1,x-1,y-1,z+1)-
IA(in1,x+1,y+1,z-1)-IA(in1,x+1,y-1,z-1)-IA(in1,x-1,y+1,z-1)-IA(in1,x-1,y-1,z-1)) ) / 38;
tmpf = sqrt(ex*ex+ey*ey+ez*ez);
if (tmpf>0)
{
ex/=(double)tmpf;
ey/=(double)tmpf;
ez/=(double)tmpf;
}
/* }}} */
if ( (!ignore_z) ||
( (ABS(ez)<ABS(ex)) && (ABS(ez)<ABS(ey)) ) )
{
/* {{{ fill 1D arrays and differentiate TLI */
im1.i=in1;
im2.i=in2;
memset(arrA-5*(CORRWIDTH+SEARCH),0,sizeof(float)*10*(CORRWIDTH+SEARCH));
memset(arrB-5*(CORRWIDTH+SEARCH),0,sizeof(float)*10*(CORRWIDTH+SEARCH));
memset(arr1-5*(CORRWIDTH+SEARCH),0,sizeof(float)*10*(CORRWIDTH+SEARCH));
memset(arr2-5*(CORRWIDTH+SEARCH),0,sizeof(float)*10*(CORRWIDTH+SEARCH));
/*IA(flow,x,y,z) = 1;*/ /* DEBUG colour edge point */
/*if ((x==53)&&(y==61)&&(z==78)) {*/ /* DEBUG */
/* printf("normal=(%f %f %f) ",ex,ey,ez);*/ /* DEBUG */
arrA[0]=IA(in1,x,y,z); /* most odd - these lines give -32768 if full optimisation is turned on */
arrB[0]=IA(in2,x,y,z);
/*IA(flow,x,y,z) = 3;*/ /* DEBUG colour central point */
pos=0;
d=1; X=FTOI(x+d*ex); Y=FTOI(y+d*ey); Z=FTOI(z+d*ez);
if ( (X>0) && (X<x_size-1) && (Y>0) && (Y<y_size-1) && (Z>0) && (Z<z_size-1) )
{
arrA[1]=TLI(im1,x+d*ex,y+d*ey,z+d*ez);
arrB[1]=TLI(im2,x+d*ex,y+d*ey,z+d*ez);
pos=-1;
segvalpos = IA(seg1,X,Y,Z);
for(d=2;d<=CORRWIDTH+SEARCH+1;d++)
{
X=FTOI(x+d*ex); Y=FTOI(y+d*ey); Z=FTOI(z+d*ez);
if ( (X>0) && (X<x_size-1) && (Y>0) && (Y<y_size-1) && (Z>0) && (Z<z_size-1) )
{
if ( (pos<0) && (IA(seg1,X,Y,Z)!=segvalpos) )
pos=d-1;
arrA[d]=TLI(im1,x+d*ex,y+d*ey,z+d*ez);
arrB[d]=TLI(im2,x+d*ex,y+d*ey,z+d*ez);
}
else
break;
}
if ( (pos<0) || (pos>CORRWIDTH) )
pos=CORRWIDTH;
if (pos==d-1)
pos=d-2;
}
/* {{{ COMMENT DEBUG draw search space */
#ifdef FoldingComment
for(d=1;d<=SEARCH+pos;d++)
{
X=FTOI(x+d*ex); Y=FTOI(y+d*ey); Z=FTOI(z+d*ez);
if (d<=pos)
IA(flow,X,Y,Z) = 7;
else
IA(flow,X,Y,Z) = 5;
}
#endif
/* }}} */
neg=0;
d=-1; X=FTOI(x+d*ex); Y=FTOI(y+d*ey); Z=FTOI(z+d*ez);
if ( (X>0) && (X<x_size-1) && (Y>0) && (Y<y_size-1) && (Z>0) && (Z<z_size-1) )
{
arrA[-1]=TLI(im1,x+d*ex,y+d*ey,z+d*ez);
arrB[-1]=TLI(im2,x+d*ex,y+d*ey,z+d*ez);
neg=1;
segvalneg = IA(seg1,X,Y,Z);
for(d=-2;d>=-CORRWIDTH-SEARCH-1;d--)
{
X=FTOI(x+d*ex); Y=FTOI(y+d*ey); Z=FTOI(z+d*ez);
if ( (X>0) && (X<x_size-1) && (Y>0) && (Y<y_size-1) && (Z>0) && (Z<z_size-1) )
{
if ( (neg>0) && (IA(seg1,X,Y,Z)!=segvalneg) )
neg=d+1;
arrA[d]=TLI(im1,x+d*ex,y+d*ey,z+d*ez);
arrB[d]=TLI(im2,x+d*ex,y+d*ey,z+d*ez);
}
else
break;
}
if ( (neg>0) || (neg<-CORRWIDTH) )
neg=-CORRWIDTH;
if (neg==d+1)
neg=d+2;
}
/* {{{ COMMENT DEBUG draw search space */
#ifdef FoldingComment
for(d=-1;d>=-SEARCH+neg;d--)
{
X=FTOI(x+d*ex); Y=FTOI(y+d*ey); Z=FTOI(z+d*ez);
if (d>=neg)
IA(flow,X,Y,Z) = 7;
else
IA(flow,X,Y,Z) = 5;
}
#endif
/* }}} */
/*printf("<%d %d %d %d> ",neg,pos,(int)segvalneg,(int)segvalpos);*/ /* DEBUG*/
for(d=-SEARCH-CORRWIDTH-1;d<=SEARCH+CORRWIDTH+1;d++)
{
arr1[d]=exp(-0.5*pow((2.0*d-neg-pos)/(pos-neg),4.0)) * (arrA[d+1]-arrA[d-1]);
arr2[d]=exp(-0.5*pow((2.0*d-neg-pos)/(pos-neg),4.0)) * (arrB[d+1]-arrB[d-1]);
/*printf("[%d %d %d %d %d] ",d,(int)arrA[d],(int)arrB[d],(int)arr1[d],(int)arr2[d]);*/ /* DEBUG */
}
/* }}} */
/* {{{ find position of maximum correlation */
for(r=-SEARCH, maxss=0, rrr=0; r<=SEARCH; r++)
{
for(rr=neg, ss=0; rr<=pos; rr++)
ss+=arr1[rr]*arr2[rr+r];
arrA[r]=ss;
/* printf("[%d %.2f] ",r,ss);*/ /* DEBUG */
if ( (ss>maxss) && (r>-SEARCH) && (r<SEARCH) )
{
maxss=ss;
rrr=r;
}
}
/* now find this max to sub-voxel accuracy */
tmpf = arrA[rrr+1] + arrA[rrr-1] - 2*arrA[rrr];
if (tmpf!=0)
tmpf = 0.5 * (arrA[rrr-1]-arrA[rrr+1]) / tmpf;
if ( (tmpf<-0.5) || (tmpf>0.5) ) /* protect against sub-voxel fit not making sense */
tmpf=0;
else
tmpf+=rrr;
tmpf = (segvalneg-segvalpos)*tmpf; /* use segmentation info to get directionality */
/*printf(" tmpf=%f\n",tmpf);*/ /* DEBUG */
IA(flow,x,y,z) = tmpf; /* turn off if DEBUGging */
total += tmpf;
count ++;
/*}*/ /* DEBUG */
/* }}} */
}
}
}
/* }}} */
/* {{{ output flow image */
if (flow_output)
{
im.i=flow;
im.dt=DT_FLOAT;
sprintf(thestring,"%s_to_%s_flow",argv[1],argv[2]);
avw_write(thestring,im);
}
/* }}} */
/* {{{ final outputs */
printf("AREA %.4f mm^2\n", count*voxel_area );
printf("VOLC %.4f mm^3\n", total*voxel_volume );
printf("RATIO %.4f mm\n", (total*voxel_volume) / (count*voxel_area) ); /* mean perpendicular edge motion; l in the equations */
printf("PBVC %.4f %%\n", (calib*30*total*voxel_volume) / (count*voxel_area) );
/* }}} */
/* {{{ COMMENT remove resampled intermediate files */
#ifdef FoldingComment
if ( ! debug )
{
sprintf(thestring,"/bin/rm -f %s_halfwayto_%s.img %s_halfwayto_%s.hdr %s_halfwayto_%s.img %s_halfwayto_%s.hdr %s_halfwayto_%s_mask.* %s_halfwayto_%s_mask.* %s_halfwayto_%s_brain.img %s_halfwayto_%s_brain.hdr %s_halfwayto_%s_brain_seg.img %s_halfwayto_%s_brain_seg.hdr",
argv[1],argv[2],argv[1],argv[2],argv[2],argv[1],argv[2],argv[1],argv[1],argv[2],argv[2],argv[1],
argv[1],argv[2],argv[1],argv[2],argv[1],argv[2],argv[1],argv[2]);
system(thestring);
}
#endif
/* }}} */
return 0;
}
/* }}} */
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment