From cba2a09b42c867ccef2de3cbdd9c26d0af625fef Mon Sep 17 00:00:00 2001
From: Stephen Smith <steve@fmrib.ox.ac.uk>
Date: Sat, 4 Jun 2005 06:27:54 +0000
Subject: [PATCH] *** empty log message ***
---
avwmaths.c | 368 ++++++++++++++++++++++++++++++++++++++++++++++++++---
1 file changed, 348 insertions(+), 20 deletions(-)
diff --git a/avwmaths.c b/avwmaths.c
index 0fa4034..f3f36dd 100644
--- a/avwmaths.c
+++ b/avwmaths.c
@@ -113,7 +113,7 @@ double tmpd;
if (isupper((int)argv[i][1]))
/* {{{ group operations (dimensionality reduction) */
-#define CI FDT tmpval, *tmpdata=malloc(sizeof(FDT)*MAX(MAX(im.x,im.y),MAX(im.z,im.t))); double sum=0,sumsquares=0,min=im.dtmax,max=im.dtmin; int nmax=0, n
+#define CI FDT tmpval, *tmpdata=malloc(sizeof(FDT)*MAX(MAX(im.x,im.y),MAX(im.z,im.t))); double sum=0,sumsquares=0,min=im.dtmax,max=im.dtmin; int nmax=0, n=0
#define CU { if (argv[i][1] == 'T') n=t; \
else if (argv[i][1] == 'Z') n=z; \
@@ -364,8 +364,7 @@ double tmpd;
{
for(j=0;j<size;j++)
- if (in[0][j]<0)
- in[0][j]=-in[0][j];
+ in[0][j]=ABS(in[0][j]);
}
/* }}} */
@@ -601,7 +600,27 @@ double tmpd;
/* }}} */
else if (!strncmp(argv[i], "-nms", 4))
- /* {{{ nms */
+ /* {{{ non-max-suppression */
+
+{
+ int maxsearch=0;
+ float origthresh, *tmpdata;
+ image_struct distancemap, thedata;
+
+ if (!strncmp(argv[i], "-nmss", 5))
+ {
+ maxsearch=1;
+ i++;
+ origthresh=atof(argv[i]);
+ i++;
+ avw_read(argv[i],&distancemap);
+ i++;
+ avw_read(argv[i],&thedata);
+
+ tmpdata=calloc(thedata.t*size,sizeof(float));
+ }
+
+ /* {{{ nms 3x3x3 CofG etc with flow smoothing */
{
int MAXZ=1;
@@ -610,35 +629,344 @@ double tmpd;
if (im.z<3)
MAXZ=0;
- for(t=0;t<im.t;t++) for(z=MAXZ;z<im.z-MAXZ;z++) for(y=1;y<im.y-1;y++) for(x=1;x<im.x-1;x++)
+ for(t=0;t<im.t;t++) /* kill this - shouldn't do this over t */
+ {
+ /* {{{ setup temp images to store perp */
+
+ short *X=calloc(size,sizeof(short)),
+ *Y=calloc(size,sizeof(short)),
+ *Z=calloc(size,sizeof(short)),
+ *XX=calloc(size,sizeof(short)),
+ *YY=calloc(size,sizeof(short)),
+ *ZZ=calloc(size,sizeof(short));
+
+/* }}} */
+
+ /* {{{ estimate perp from CofG and curvature, and store this in X,Y,Z */
+
+ for(z=MAXZ;z<im.z-MAXZ;z++) for(y=1;y<im.y-1;y++) for(x=1;x<im.x-1;x++)
+ {
+ FDT theval = in[0][t*xyzsize+z*xysize+y*im.x+x];
+
+ if ( theval != 0 )
+ {
+ float CofGx=0, CofGy=0, CofGz=0, Sum=0, CofGl;
+ int xx, yy, zz, xxx=0, yyy=0, zzz=0;
+
+ for(zz=-MAXZ; zz<=MAXZ; zz++)
+ for(yy=-1; yy<=1; yy++)
+ for(xx=-1; xx<=1; xx++)
+ {
+ float val = in[0][t*xyzsize+(z+zz)*xysize+(y+yy)*im.x+x+xx];
+ Sum += val;
+ CofGx += xx * val; CofGy += yy * val; CofGz += zz * val;
+ }
+
+ CofGx /= Sum; CofGy /= Sum; CofGz /= Sum;
+ CofGl = sqrt(CofGx*CofGx+CofGy*CofGy+CofGz*CofGz);
+
+ if (CofGl > .1) /* is CofG far enough away from centre voxel? */
+ {
+ xxx = FTOI(CofGx/CofGl);
+ yyy = FTOI(CofGy/CofGl);
+ zzz = FTOI(CofGz/CofGl);
+ }
+ else
+ /* {{{ find direction of max curvature */
+
+ {
+ float maxcost=0, centreval=2.0 * (float)theval;
+
+ for(zz=0; zz<=MAXZ; zz++) /* note - starts at zero as we're only searching half the voxels */
+ for(yy=-1; yy<=1; yy++)
+ for(xx=-1; xx<=1; xx++)
+ if ( (zz==1) || (yy==1) || ((yy==0)&&(xx==1)) ) /* only search half the voxels */
+ {
+ float weighting = pow( (float)(xx*xx+yy*yy+zz*zz) , -0.7 ); /* power is arbitrary: maybe test other functions here */
+ float cost = weighting * ( centreval
+ - (float)in[0][t*xyzsize+(z+zz)*xysize+(y+yy)*im.x+x+xx]
+ - (float)in[0][t*xyzsize+(z-zz)*xysize+(y-yy)*im.x+x-xx] );
+
+ if (cost>maxcost)
+ {
+ maxcost=cost;
+ xxx=xx;
+ yyy=yy;
+ zzz=zz;
+ }
+ }
+ }
+
+/* }}} */
+
+ X[z*xysize+y*im.x+x]=xxx;
+ Y[z*xysize+y*im.x+x]=yyy;
+ Z[z*xysize+y*im.x+x]=zzz;
+ }
+
+ }
+
+ /* {{{ COMMENT save perp image */
+
+#ifdef FoldingComment
+
+{
+ char tmpname[10000];
+ image_struct tmpim=im;
+
+ tmpim.i=calloc(size*3,sizeof(float));
+ sprintf(tmpname,"%s_flow",argv[argc-1]);
+
+ tmpim.dt=DT_FLOAT;
+ tmpim.t=3;
+
+ for(z=0;z<im.z;z++) for(y=0;y<im.y;y++) for(x=0;x<im.x;x++)
+ {
+ float tmpX=X[z*xysize+y*im.x+x],
+ tmpY=Y[z*xysize+y*im.x+x],
+ tmpZ=Z[z*xysize+y*im.x+x],
+ tmpf=sqrt(tmpX*tmpX+tmpY*tmpY+tmpZ*tmpZ);
+
+ if (tmpf>0)
+ {
+ tmpim.i[ z*xysize+y*im.x+x]=tmpX/tmpf;
+ tmpim.i[ xyzsize+z*xysize+y*im.x+x]=tmpY/tmpf;
+ tmpim.i[2*xyzsize+z*xysize+y*im.x+x]=tmpZ/tmpf;
+ }
+ }
+
+ avw_write(tmpname,tmpim);
+ free(tmpim.i);
+}
+
+#endif
+
+/* }}} */
+
+/* }}} */
+ /* {{{ smooth X,Y,Z and store in XX,YY,ZZ */
+
+ for(z=MAXZ;z<im.z-MAXZ;z++) for(y=1;y<im.y-1;y++) for(x=1;x<im.x-1;x++)
+ {
+ int xx, yy, zz, *localsum=calloc(27,sizeof(int)), localmax=0, xxx, yyy, zzz;
+
+ for(zz=-MAXZ; zz<=MAXZ; zz++)
+ for(yy=-1; yy<=1; yy++)
+ for(xx=-1; xx<=1; xx++)
+ {
+ xxx = X[(z+zz)*xysize+(y+yy)*im.x+x+xx];
+ yyy = Y[(z+zz)*xysize+(y+yy)*im.x+x+xx];
+ zzz = Z[(z+zz)*xysize+(y+yy)*im.x+x+xx];
+ localsum[(1+zzz)*9+(1+yyy)*3+1+xxx]++;
+ localsum[(1-zzz)*9+(1-yyy)*3+1-xxx]++;
+ }
+
+ for(zz=-MAXZ; zz<=MAXZ; zz++)
+ for(yy=-1; yy<=1; yy++)
+ for(xx=-1; xx<=1; xx++)
+ {
+ if (localsum[(1+zz)*9+(1+yy)*3+1+xx]>localmax)
+ {
+ localmax=localsum[(1+zz)*9+(1+yy)*3+1+xx];
+ XX[z*xysize+y*im.x+x]=xx;
+ YY[z*xysize+y*im.x+x]=yy;
+ ZZ[z*xysize+y*im.x+x]=zz;
+ }
+ }
+ }
+
+ /* {{{ COMMENT save perp image */
+
+#ifdef FoldingComment
+
+{
+ char tmpname[10000];
+ image_struct tmpim=im;
+
+ tmpim.i=calloc(size*3,sizeof(float));
+ sprintf(tmpname,"%s_flowsmooth",argv[argc-1]);
+
+ tmpim.dt=DT_FLOAT;
+ tmpim.t=3;
+
+ for(z=0;z<im.z;z++) for(y=0;y<im.y;y++) for(x=0;x<im.x;x++)
{
- double cost, mincost=1e10;
- int xx, yy, zz, xxx=0, yyy=0, zzz=0;
+ float tmpX=XX[z*xysize+y*im.x+x],
+ tmpY=YY[z*xysize+y*im.x+x],
+ tmpZ=ZZ[z*xysize+y*im.x+x],
+ tmpf=sqrt(tmpX*tmpX+tmpY*tmpY+tmpZ*tmpZ);
+
+ if (tmpf>0)
+ {
+ tmpim.i[ z*xysize+y*im.x+x]=tmpX/tmpf;
+ tmpim.i[ xyzsize+z*xysize+y*im.x+x]=tmpY/tmpf;
+ tmpim.i[2*xyzsize+z*xysize+y*im.x+x]=tmpZ/tmpf;
+ }
+ }
+
+ avw_write(tmpname,tmpim);
+ free(tmpim.i);
+}
+
+#endif
+
+/* }}} */
+
+/* }}} */
+ /* {{{ do non-max-suppression in the direction of perp */
+
+ for(z=MAXZ;z<im.z-MAXZ;z++) for(y=1;y<im.y-1;y++) for(x=1;x<im.x-1;x++)
+ {
+ FDT theval = in[0][t*xyzsize+z*xysize+y*im.x+x];
+ int xxx=XX[z*xysize+y*im.x+x];
+ int yyy=YY[z*xysize+y*im.x+x];
+ int zzz=ZZ[z*xysize+y*im.x+x];
+
+ if ( ( (xxx!=0) || (yyy!=0) || (zzz!=0) ) &&
+ ( theval >= in[0][t*xyzsize+(z+zzz)*xysize+(y+yyy)*im.x+x+xxx] ) &&
+ ( theval > in[0][t*xyzsize+(z-zzz)*xysize+(y-yyy)*im.x+x-xxx] ) &&
+ ( theval >= in[0][t*xyzsize+(z+2*zzz)*xysize+(y+2*yyy)*im.x+x+2*xxx] ) &&
+ ( theval > in[0][t*xyzsize+(z-2*zzz)*xysize+(y-2*yyy)*im.x+x-2*xxx] ) )
+ tmpim[t*xyzsize+z*xysize+y*im.x+x] = theval;
+ }
+
+/* }}} */
+
+ if (maxsearch==1)
+ /* {{{ do search in the direction of perp */
+
+#define SEARCHSIGMA 5 /* length in linear voxel dimensions */
+#define MAXSEARCHLENGTH (3*SEARCHSIGMA)
+
+{
+ int d, T, iters;
+
+ for(T=0;T<thedata.t;T++) for(z=MAXZ;z<im.z-MAXZ;z++) for(y=1;y<im.y-1;y++) for(x=1;x<im.x-1;x++)
+ if (tmpim[t*xyzsize+z*xysize+y*im.x+x] > origthresh)
+ {
+ int xxx=XX[z*xysize+y*im.x+x];
+ int yyy=YY[z*xysize+y*im.x+x];
+ int zzz=ZZ[z*xysize+y*im.x+x];
+ float maxval=thedata.i[T*xyzsize+z*xysize+y*im.x+x];
+ float maxval_weighted=maxval;
+ int maxvalX=x, maxvalY=y, maxvalZ=z;
+ float exponentfactor = -0.5 * (xxx*xxx+yyy*yyy+zzz*zzz) / (float)(SEARCHSIGMA*SEARCHSIGMA);
+
+ for(iters=0;iters<2;iters++)
+ {
+ float distance=0;
+
+ for(d=1;d<MAXSEARCHLENGTH;d++)
+ {
+ int D=d;
+ if (iters==1) D=-d;
+
+ if (distancemap.i[t*xyzsize+(z+zzz*D)*xysize+(y+yyy*D)*im.x+x+xxx*D]>=distance)
+ {
+ float distanceweight = exp(d * d * exponentfactor);
+ distance=distancemap.i[t*xyzsize+(z+zzz*D)*xysize+(y+yyy*D)*im.x+x+xxx*D];
+ if (distanceweight * thedata.i[T*xyzsize+(z+zzz*D)*xysize+(y+yyy*D)*im.x+x+xxx*D]>maxval_weighted)
+ {
+ maxval=thedata.i[T*xyzsize+(z+zzz*D)*xysize+(y+yyy*D)*im.x+x+xxx*D];
+ maxval_weighted=maxval*distanceweight;
+ maxvalX=x+xxx*D;
+ maxvalY=y+yyy*D;
+ maxvalZ=z+zzz*D;
+ }
+ }
+ else
+ d=MAXSEARCHLENGTH;
+ }
+ }
+
+ tmpdata[T*xyzsize+z*xysize+y*im.x+x]=maxval; /* output maxsearch data */
- for(zz=0; zz<=MAXZ; zz++)
- for(yy=-1; yy<=1; yy++)
- for(xx=-1; xx<=1; xx++)
- if ( (zz==1) || (yy==1) || ((yy==0)&&(xx==1)) )
+ /*tmpdata[T*xyzsize+z*xysize+y*im.x+x]=maxvalX-x;*/ /* output X component of search vector */
+ /*tmpdata[T*xyzsize+z*xysize+y*im.x+x]=maxvalY-y;*/ /* output Y component of search vector */
+ /*tmpdata[T*xyzsize+z*xysize+y*im.x+x]=maxvalZ-z;*/ /* output Z component of search vector */
+ }
+
+ thedata.i=tmpdata;
+ avw_write(argv[argc-1],thedata);
+
+ exit(0);
+}
+
+/* }}} */
+ /* {{{ COMMENT do search in the direction of perp */
+
+#ifdef FoldingComment
+
+/*#define MAXSEARCHLENGTH 20*/
+
+#define SEARCHSIGMA 5 /* length in linear voxel dimensions */
+#define MAXSEARCHLENGTH (3*SEARCHSIGMA)
+
+{
+ int d, T, iters;
+
+ for(T=0;T<thedata.t;T++) for(z=MAXZ;z<im.z-MAXZ;z++) for(y=1;y<im.y-1;y++) for(x=1;x<im.x-1;x++)
+ if (tmpim[t*xyzsize+z*xysize+y*im.x+x] > origthresh)
+ {
+ int xxx=XX[z*xysize+y*im.x+x];
+ int yyy=YY[z*xysize+y*im.x+x];
+ int zzz=ZZ[z*xysize+y*im.x+x];
+ float maxval=thedata.i[T*xyzsize+z*xysize+y*im.x+x];
+ int maxvalX=x, maxvalY=y, maxvalZ=z;
+
+ for(iters=0;iters<2;iters++)
+ {
+ float distance=0;
+
+ for(d=1;d<MAXSEARCHLENGTH;d++)
{
- cost = in[0][t*xyzsize+(z+zz)*xysize+(y+yy)*im.x+x+xx] + in[0][t*xyzsize+(z-zz)*xysize+(y-yy)*im.x+x-xx];
- if (cost<mincost)
+ int D=d;
+ if (iters==1) D=-d;
+
+ if (distancemap.i[t*xyzsize+(z+zzz*D)*xysize+(y+yyy*D)*im.x+x+xxx*D]>=distance)
{
- mincost=cost;
- xxx=xx;
- yyy=yy;
- zzz=zz;
+ distance=distancemap.i[t*xyzsize+(z+zzz*D)*xysize+(y+yyy*D)*im.x+x+xxx*D];
+ if (thedata.i[T*xyzsize+(z+zzz*D)*xysize+(y+yyy*D)*im.x+x+xxx*D]>maxval)
+ {
+ maxval=thedata.i[T*xyzsize+(z+zzz*D)*xysize+(y+yyy*D)*im.x+x+xxx*D];
+ maxvalX=x+xxx*D;
+ maxvalY=y+yyy*D;
+ maxvalZ=z+zzz*D;
+ }
}
+ else
+ d=MAXSEARCHLENGTH;
}
+ }
+
+ tmpdata[T*xyzsize+z*xysize+y*im.x+x]=maxval; /* output maxsearch data */
+
+ /*tmpdata[T*xyzsize+z*xysize+y*im.x+x]=maxvalX-x;*/ /* output X component of search vector */
+ /*tmpdata[T*xyzsize+z*xysize+y*im.x+x]=maxvalY-y;*/ /* output Y component of search vector */
+ /*tmpdata[T*xyzsize+z*xysize+y*im.x+x]=maxvalZ-z;*/ /* output Z component of search vector */
+ /*tmpdata[T*xyzsize+maxvalZ*xysize+maxvalY*im.x+maxvalX]++;*/ /* output number of times skel has looked in this original space voxel */
+ }
+
+ thedata.i=tmpdata;
+ avw_write(argv[argc-1],thedata);
+
+ exit(0);
+}
- if ( ( in[0][t*xyzsize+z*xysize+y*im.x+x] >= in[0][t*xyzsize+(z+zzz)*xysize+(y+yyy)*im.x+x+xxx] ) &&
- ( in[0][t*xyzsize+z*xysize+y*im.x+x] > in[0][t*xyzsize+(z-zzz)*xysize+(y-yyy)*im.x+x-xxx] ) )
- tmpim[t*xyzsize+z*xysize+y*im.x+x] = in[0][t*xyzsize+z*xysize+y*im.x+x];
+#endif
+
+/* }}} */
+
+ free(X); free(Y); free(Z); free(XX); free(YY); free(ZZ);
}
free(in[0]);
in[0]=tmpim;
}
+/* }}} */
+}
+
/* }}} */
else if (!strncmp(argv[i], "-nanm", 5))
/* {{{ NaN mask */
--
GitLab