Skip to content
Snippets Groups Projects
Commit cba2a09b authored by Stephen Smith's avatar Stephen Smith
Browse files

*** empty log message ***

parent 021dfb4b
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 348 additions and 20 deletions
avwmaths.c
+ 348
− 20
View file @ cba2a09b
......@@ -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 */
......
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