diff --git a/Makefile b/Makefile
index 8b59d6fc4fabe9e0d5a89a284885e807a4bb5f6e..a2e50c3aeaa6d4d2b74bc439d84da140721d1eda 100644
--- a/Makefile
+++ b/Makefile
@@ -19,10 +19,8 @@ SCRIPTS = siena siena_flirt siena_cal sienax siena_flow2std
all: ${XFILES} docscripts
# DON'T REMOVE THE -O0!!!! siena_diff doesn't work with optimisation
-siena_diff: siena_diff.c
- $(CC) $(CFLAGS) -O0 -DFDT="float" -o siena_diff siena_diff.c $(LDFLAGS) $(LIBS)
-siena_diffCC: siena_diff.cc
- ${CXX} $(LDFLAGS) -o $@ ${SN_OBJS} ${LIBCC}
+siena_diff: siena_diff.cc
+ ${CXX} ${CXXFLAGS} ${LDFLAGS} -O0 -o siena_diff siena_diff.cc ${LIBCC}
surface_norm: ${SN_OBJS}
${CXX} ${LDFLAGS} -o $@ ${SN_OBJS} ${LIBCC}
diff --git a/siena b/siena
index f4e9b495a0efa017a3a79d5998ffc505fd13dd45..9669ae24abcdbdb456f4c2e5b88a4c4d4dc35c3b 100755
--- a/siena
+++ b/siena
@@ -141,6 +141,8 @@ if [ $dostd = 1 ] ; then
if [ $stdmask = 1 ] ; then
${FSLDIR}/bin/flirt -in ${FSLDIR}/etc/standard/avg152T1_brain_mask_dil2 -ref $A -out ${A}_stdmask -applyxfm -init ${A}_to_std_inv.mat
${FSLDIR}/bin/flirt -in ${FSLDIR}/etc/standard/avg152T1_brain_mask_dil2 -ref $B -out ${B}_stdmask -applyxfm -init ${B}_to_std_inv.mat
+ ${FSLDIR}/bin/avwmaths ${A}_stdmask -thr 0.5 -bin ${A}_stdmask
+ ${FSLDIR}/bin/avwmaths ${B}_stdmask -thr 0.5 -bin ${B}_stdmask
${FSLDIR}/bin/avwmaths ${A}_brain_mask -mas ${A}_stdmask ${A}_brain_mask
${FSLDIR}/bin/avwmaths ${B}_brain_mask -mas ${B}_stdmask ${B}_brain_mask
fi
@@ -149,6 +151,8 @@ if [ $dostd = 1 ] ; then
${FSLDIR}/bin/avwmaths ${FSLDIR}/etc/standard/avg152T1_brain_mask -mul 0 -add 1 $stdroi ${A}_and_${B}_stdmask
${FSLDIR}/bin/flirt -in ${A}_and_${B}_stdmask -ref $A -out ${A}_stdmask -applyxfm -init ${A}_to_std_inv.mat
${FSLDIR}/bin/flirt -in ${A}_and_${B}_stdmask -ref $B -out ${B}_stdmask -applyxfm -init ${B}_to_std_inv.mat
+ ${FSLDIR}/bin/avwmaths ${A}_stdmask -thr 0.5 -bin ${A}_stdmask
+ ${FSLDIR}/bin/avwmaths ${B}_stdmask -thr 0.5 -bin ${B}_stdmask
${FSLDIR}/bin/avwmaths ${A}_valid_mask_with_$B -mul ${A}_stdmask ${A}_valid_mask_with_$B
${FSLDIR}/bin/avwmaths ${B}_valid_mask_with_$A -mul ${B}_stdmask ${B}_valid_mask_with_$A
fi
diff --git a/siena_diff.c b/siena_diff.c
deleted file mode 100644
index 3176c42ebf82cec7abef7fe1ae2fc705cbc309bc..0000000000000000000000000000000000000000
--- a/siena_diff.c
+++ /dev/null
@@ -1,871 +0,0 @@
-/* {{{ 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>_stdmask 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)
- /*if (seg1[x]>2)*/
- 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) );
-
-/* }}} */
-
- return 0;
-}
-
-/* }}} */
diff --git a/siena_diff.cc b/siena_diff.cc
new file mode 100644
index 0000000000000000000000000000000000000000..d37f9db63089cedb6a75f4b3d491109e99278059
--- /dev/null
+++ b/siena_diff.cc
@@ -0,0 +1,492 @@
+// {{{ Copyright etc.
+
+/* siena_diff - compute brain change using edge motion or segmentation
+
+ Stephen Smith, FMRIB Image Analysis Group
+
+ Copyright (C) 1999-2006 University of Oxford */
+
+/* CCOPYRIGHT */
+
+// }}}
+// {{{ includes and options
+
+#define _GNU_SOURCE 1
+#define POSIX_SOURCE 1
+
+#include "newimage/newimageall.h"
+#include "miscmaths/miscmaths.h"
+
+using namespace MISCMATHS;
+using namespace NEWIMAGE;
+
+// }}}
+// {{{ usage
+
+void usage()
+{
+ cout << "\nUsage: siena_diff <input1_basename> <input2_basename> [options] [-s segmentation options]\n\n" <<
+ "[-d] debug - generate edge images and don't remove temporary images\n" <<
+ "[-2] don't segment grey+white separately (because there is poor grey-white contrast)\n" <<
+ "[-c <corr>] apply self-calibrating correction factor\n" <<
+ // "[-e] erode joint mask a lot instead of dilating it slightly (ie find ventricle surface)\n" <<
+ "[-i] ignore flow in z (may be good if top of brain is missing)\n" <<
+ "[-m] apply <input1_basename>_stdmask to brain edge points\n" <<
+ "[-s <options>] <options> to be passed to segmentation (type \"fast\" to get these)\n\n" << endl;
+ exit(1);
+}
+
+// }}}
+// {{{ main(argc, argv)
+
+#define CORRWIDTH 3
+#define SEARCH 4
+#define CS (10*(CORRWIDTH+SEARCH))
+
+int main(int argc,char *argv[])
+{
+ // {{{ vars
+
+char thestring[10000], segoptions[10000], fsldir[10000];
+int x_size, y_size, z_size, size, x, y, z, i, count,
+ seg2=0, ignore_z=0, ignore_top_slices=0, //erode_mask=0,
+ ignore_bottom_slices=0, debug=0, flow_output=1, edge_masking=0;
+float tmpf, calib=1.0, ex, ey, ez;
+ColumnVector arrA(2*CS+1), arrB(2*CS+1), arr1(2*CS+1), arr2(2*CS+1);
+double total, voxel_volume, voxel_area;
+
+// }}}
+
+ // {{{ process arguments
+
+if (argc<3)
+ usage();
+
+string argv1(argv[1]), argv2(argv[2]);
+
+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
+
+cout << "reading and combining transformed masks" << endl;
+volume<float> mask;
+read_volume(mask,argv1+"_halfwayto_"+argv2+"_mask");
+
+// setup header sizes etc.
+x_size=mask.xsize();
+y_size=mask.ysize();
+z_size=mask.zsize();
+size=x_size*y_size*z_size;
+voxel_volume = abs( mask.xdim() * mask.ydim() * mask.zdim() );
+voxel_area = pow(voxel_volume,((double)0.6666667));
+cout << "final image dimensions = " << x_size << " " << y_size << " " << z_size << ", voxel volume = " << voxel_volume << "mm^3, voxel area = " << voxel_area << "mm^2" << endl;
+
+// read mask 2 and combine with mask 1
+volume<float> mask2;
+read_volume(mask2,argv2+"_halfwayto_"+argv1+"_mask");
+mask=mask+mask2;
+mask2.destroy();
+mask.binarise(0.5);
+
+cout << "dilating/eroding combined mask" << endl;
+// if (erode_mask)
+// {
+// volume<float>kernel=spherical_kernel(17,mask.xdim(),mask.ydim(),mask.zdim());
+// mask=morphfilter(mask,kernel,"erodeS");
+// }
+// else
+ {
+ volume<float>kernel=box_kernel(3,3,3);
+ mask=morphfilter(mask,kernel,"dilate");
+ }
+
+cout << "reading transformed images and applying mask" << endl;
+volume<float> in1;
+read_volume(in1,argv1+"_halfwayto_"+argv2);
+in1 = (in1-in1.min()) * mask;
+mask.destroy();
+
+// }}}
+ // {{{ 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) // always done unless the above uncommented and used instead of this test
+ {
+ char segtype[100];
+ if (seg2) sprintf(segtype,"-c 2"); else segtype[0]=0;
+ cout << "saving image 1 to disk prior to segmentation" << endl;
+ save_volume(in1,argv1+"_halfwayto_"+argv2+"_brain");
+ in1.destroy();
+ 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]);
+ cout << thestring << endl;
+ // system(thestring);
+ }
+ else
+ {
+ cout << "using previously carried out segmentation" << endl;
+ in1.destroy();
+ }
+
+// }}}
+ // {{{ read segmentation output into edges1 and simplify; reread in1 and in2
+
+printf("finding brain edges\n");
+
+volume<float> seg1;
+read_volume(seg1,argv1+"_halfwayto_"+argv2+"_brain_seg");
+seg1.binarise(1.5);
+
+volume<float> m1;
+if (edge_masking)
+ read_volume(m1,argv1+"_halfwayto_"+argv2+"_valid_mask");
+
+read_volume(in1,argv1+"_halfwayto_"+argv2);
+in1.setinterpolationmethod(trilinear);
+
+volume<float> in2;
+read_volume(in2,argv2+"_halfwayto_"+argv1);
+in2.setinterpolationmethod(trilinear);
+
+// }}}
+ // {{{ find segmentation-based edges in image 1 and flow
+
+printf("finding flow\n");
+
+volume<float> flow=in1;
+flow=0;
+
+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 ( (seg1(x,y,z)>0.5) && /* not background or CSF */
+ ( (seg1(x+1,y,z)<0.5) || (seg1(x-1,y,z)<0.5) ||
+ (seg1(x,y+1,z)<0.5) || (seg1(x,y-1,z)<0.5) ||
+ (seg1(x,y,z+1)<0.5) || (seg1(x,y,z-1)<0.5) ) &&
+ ( ( ! edge_masking ) || ( m1(x,y,z)>0 ) ) )
+ {
+ int pos, neg, r, rr, rrr, d, X, Y, Z;
+ float ss, maxss, segvalpos=0, segvalneg=0;
+
+ // {{{ find local gradient and derive unit normal
+
+ ex = ( 10*(in1(x+1,y,z)-in1(x-1,y,z)) +
+ 5*(in1(x+1,y+1,z)+in1(x+1,y-1,z)+in1(x+1,y,z+1)+in1(x+1,y,z-1)-
+ in1(x-1,y+1,z)-in1(x-1,y-1,z)-in1(x-1,y,z+1)-in1(x-1,y,z-1)) +
+ 2*(in1(x+1,y+1,z+1)+in1(x+1,y-1,z+1)+in1(x+1,y+1,z-1)+in1(x+1,y-1,z-1)-
+ in1(x-1,y+1,z+1)-in1(x-1,y-1,z+1)-in1(x-1,y+1,z-1)-in1(x-1,y-1,z-1)) ) / 38;
+ ey = ( 10*(in1(x,y+1,z)-in1(x,y-1,z)) +
+ 5*(in1(x+1,y+1,z)+in1(x-1,y+1,z)+in1(x,y+1,z+1)+in1(x,y+1,z-1)-
+ in1(x+1,y-1,z)-in1(x-1,y-1,z)-in1(x,y-1,z+1)-in1(x,y-1,z-1)) +
+ 2*(in1(x+1,y+1,z+1)+in1(x-1,y+1,z+1)+in1(x+1,y+1,z-1)+in1(x-1,y+1,z-1)-
+ in1(x+1,y-1,z+1)-in1(x-1,y-1,z+1)-in1(x+1,y-1,z-1)-in1(x-1,y-1,z-1)) ) / 38;
+ ez = ( 10*(in1(x,y,z+1)-in1(x,y,z-1)) +
+ 5*(in1(x,y+1,z+1)+in1(x,y-1,z+1)+in1(x+1,y,z+1)+in1(x-1,y,z+1)-
+ in1(x,y+1,z-1)-in1(x,y-1,z-1)-in1(x+1,y,z-1)-in1(x-1,y,z-1)) +
+ 2*(in1(x+1,y+1,z+1)+in1(x+1,y-1,z+1)+in1(x-1,y+1,z+1)+in1(x-1,y-1,z+1)-
+ in1(x+1,y+1,z-1)-in1(x+1,y-1,z-1)-in1(x-1,y+1,z-1)-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
+
+arrA=0; arrB=0; arr1=0; arr2=0;
+
+/*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(CS)=in1(x,y,z);
+arrB(CS)=in2(x,y,z);
+
+/*flow(x,y,z) = 3;*/ /* DEBUG colour central point */
+
+ pos=0;
+ d=1; X=round(x+d*ex); Y=round(y+d*ey); Z=round(z+d*ez);
+ if ( (X>0) && (X<x_size-1) && (Y>0) && (Y<y_size-1) && (Z>0) && (Z<z_size-1) )
+ {
+ arrA(CS+1)=in1.interpolate(x+d*ex,y+d*ey,z+d*ez);
+ arrB(CS+1)=in2.interpolate(x+d*ex,y+d*ey,z+d*ez);
+ pos=-1;
+ segvalpos = seg1(X,Y,Z);
+ for(d=2;d<=CORRWIDTH+SEARCH+1;d++)
+ {
+ X=round(x+d*ex); Y=round(y+d*ey); Z=round(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) && (seg1(X,Y,Z)!=segvalpos) )
+ pos=d-1;
+ arrA(CS+d)=in1.interpolate(x+d*ex,y+d*ey,z+d*ez);
+ arrB(CS+d)=in2.interpolate(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=round(x+d*ex); Y=round(y+d*ey); Z=round(z+d*ez);
+ if (d<=pos)
+ flow(X,Y,Z) = 7;
+ else
+ flow(X,Y,Z) = 5;
+}
+
+#endif
+
+// }}}
+
+ neg=0;
+ d=-1; X=round(x+d*ex); Y=round(y+d*ey); Z=round(z+d*ez);
+ if ( (X>0) && (X<x_size-1) && (Y>0) && (Y<y_size-1) && (Z>0) && (Z<z_size-1) )
+ {
+ arrA(CS-1)=in1.interpolate(x+d*ex,y+d*ey,z+d*ez);
+ arrB(CS-1)=in2.interpolate(x+d*ex,y+d*ey,z+d*ez);
+ neg=1;
+ segvalneg = seg1(X,Y,Z);
+ for(d=-2;d>=-CORRWIDTH-SEARCH-1;d--)
+ {
+ X=round(x+d*ex); Y=round(y+d*ey); Z=round(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) && (seg1(X,Y,Z)!=segvalneg) )
+ neg=d+1;
+ arrA(CS+d)=in1.interpolate(x+d*ex,y+d*ey,z+d*ez);
+ arrB(CS+d)=in2.interpolate(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=round(x+d*ex); Y=round(y+d*ey); Z=round(z+d*ez);
+ if (d>=neg)
+ flow(X,Y,Z) = 7;
+ else
+ 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++)
+ {
+ float denom = max(1,pos-neg);
+ arr1(CS+d)=exp(-0.5*pow((2.0*d-neg-pos)/denom,4.0)) * (arrA(CS+d+1)-arrA(CS+d-1));
+ arr2(CS+d)=exp(-0.5*pow((2.0*d-neg-pos)/denom,4.0)) * (arrB(CS+d+1)-arrB(CS+d-1));
+ }
+
+// }}}
+ // {{{ 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(CS+rr)*arr2(CS+rr+r);
+
+ arrA(CS+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(CS+rrr+1) + arrA(CS+rrr-1) - 2*arrA(CS+rrr);
+if (tmpf!=0)
+ tmpf = 0.5 * (arrA(CS+rrr-1)-arrA(CS+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 */
+
+flow(x,y,z) = tmpf; /* turn off if DEBUGging */
+total += tmpf;
+count ++;
+
+/*}*/ /* DEBUG */
+
+// }}}
+ }
+ }
+ }
+
+// }}}
+ // {{{ final outputs
+
+if (flow_output)
+ save_volume(flow,argv1+"_to_"+argv2+"_flow");
+
+cout << "AREA " << count*voxel_area << " mm^2" << endl;
+cout << "VOLC " << total*voxel_volume << " mm^3" << endl;
+cout << "RATIO " << (total*voxel_volume) / (count*voxel_area) << " mm" << endl; /* mean perpendicular edge motion; l in the equations */
+cout << "PBVC " << (calib*30*total*voxel_volume) / (count*voxel_area) << " %%" << endl;
+
+// }}}
+ return 0;
+}
+
+// }}}
+
diff --git a/sienax b/sienax
index 01e737f32e0c7c4a6ac1800d26780f2a25d184f4..f5a05f4febec00a36053effc6aef7ec14c974856 100755
--- a/sienax
+++ b/sienax
@@ -123,12 +123,15 @@ if [ "$stdroi" != "" ] ; then
MASK=${I}_stdmaskroi
fi
${FSLDIR}/bin/flirt -in $MASK -ref ${I}_brain -out ${I}_stdmask -applyxfm -init ${I}2std_inv.mat
+${FSLDIR}/bin/avwmaths ${I}_stdmask -thr 0.5 -bin ${I}_stdmask
${FSLDIR}/bin/avwmaths ${I}_brain -mask ${I}_stdmask ${I}_stdmaskbrain
if [ $regional = 1 ] ; then
${FSLDIR}/bin/flirt -in ${FSLDIR}/etc/standard/avg152T1_strucseg_periph -ref ${I}_brain -out ${I}_stdmask_segperiph -applyxfm -init ${I}2std_inv.mat
+ ${FSLDIR}/bin/avwmaths ${I}_stdmask_segperiph -thr 0.5 -bin ${I}_stdmask_segperiph
${FSLDIR}/bin/avwmaths ${FSLDIR}/etc/standard/avg152T1_strucseg -thr 4.5 -bin ${I}_tmpmask
${FSLDIR}/bin/flirt -in ${I}_tmpmask -ref ${I}_brain -out ${I}_stdmask_segvent -applyxfm -init ${I}2std_inv.mat
+ ${FSLDIR}/bin/avwmaths ${I}_stdmask_segvent -thr 0.5 -bin ${I}_stdmask_segvent
/bin/rm ${I}_tmpmask*
fi