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

*** empty log message ***

parent a6d29838
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 57 additions and 46 deletions
avwmaths.c
+ 57
− 46
View file @ 59b1e80e
...@@ -628,28 +628,21 @@ double tmpd; ...@@ -628,28 +628,21 @@ double tmpd;
{ {
int MAXZ=1; int MAXZ=1;
FDT *tmpim=calloc(size,sizeof(FDT)); FDT *tmpim=calloc(size,sizeof(FDT));
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));
if (im.z<3) if (im.z<3)
MAXZ=0; MAXZ=0;
for(t=0;t<im.t;t++) /* kill this - shouldn't do this over t */ /* {{{ estimate perp from CofG and curvature, and store this in X,Y,Z */
{
/* {{{ 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++) 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]; FDT theval = in[0][z*xysize+y*im.x+x];
if ( theval != 0 ) if ( theval != 0 )
{ {
...@@ -660,7 +653,7 @@ double tmpd; ...@@ -660,7 +653,7 @@ double tmpd;
for(yy=-1; yy<=1; yy++) for(yy=-1; yy<=1; yy++)
for(xx=-1; xx<=1; xx++) for(xx=-1; xx<=1; xx++)
{ {
float val = in[0][t*xyzsize+(z+zz)*xysize+(y+yy)*im.x+x+xx]; float val = in[0][(z+zz)*xysize+(y+yy)*im.x+x+xx];
Sum += val; Sum += val;
CofGx += xx * val; CofGy += yy * val; CofGz += zz * val; CofGx += xx * val; CofGy += yy * val; CofGz += zz * val;
} }
...@@ -687,8 +680,8 @@ double tmpd; ...@@ -687,8 +680,8 @@ double tmpd;
{ {
float weighting = pow( (float)(xx*xx+yy*yy+zz*zz) , -0.7 ); /* power is arbitrary: maybe test other functions here */ 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 cost = weighting * ( centreval
- (float)in[0][t*xyzsize+(z+zz)*xysize+(y+yy)*im.x+x+xx] - (float)in[0][(z+zz)*xysize+(y+yy)*im.x+x+xx]
- (float)in[0][t*xyzsize+(z-zz)*xysize+(y-yy)*im.x+x-xx] ); - (float)in[0][(z-zz)*xysize+(y-yy)*im.x+x-xx] );
if (cost>maxcost) if (cost>maxcost)
{ {
...@@ -747,7 +740,7 @@ double tmpd; ...@@ -747,7 +740,7 @@ double tmpd;
/* }}} */ /* }}} */
/* }}} */ /* }}} */
/* {{{ smooth X,Y,Z and store in XX,YY,ZZ */ /* {{{ 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++) for(z=MAXZ;z<im.z-MAXZ;z++) for(y=1;y<im.y-1;y++) for(x=1;x<im.x-1;x++)
{ {
...@@ -815,33 +808,35 @@ double tmpd; ...@@ -815,33 +808,35 @@ double tmpd;
/* }}} */ /* }}} */
free(X); free(Y); free(Z);
/* }}} */ /* }}} */
/* {{{ do non-max-suppression in the direction of perp */ /* {{{ 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++) 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]; FDT theval = in[0][z*xysize+y*im.x+x];
int xxx=XX[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 yyy=YY[z*xysize+y*im.x+x];
int zzz=ZZ[z*xysize+y*im.x+x]; int zzz=ZZ[z*xysize+y*im.x+x];
if ( ( (xxx!=0) || (yyy!=0) || (zzz!=0) ) && if ( ( (xxx!=0) || (yyy!=0) || (zzz!=0) ) &&
( theval >= in[0][t*xyzsize+(z+zzz)*xysize+(y+yyy)*im.x+x+xxx] ) && ( theval >= in[0][(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][(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][(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] ) ) ( theval > in[0][(z-2*zzz)*xysize+(y-2*yyy)*im.x+x-2*xxx] ) )
tmpim[t*xyzsize+z*xysize+y*im.x+x] = theval; tmpim[z*xysize+y*im.x+x] = theval;
} }
/* }}} */ /* }}} */
if (maxsearch==1) if (maxsearch==1)
/* {{{ do search in the direction of perp */ /* {{{ do search in the direction of perp */
#define SEARCHSIGMA 5 /* length in linear voxel dimensions */ #define SEARCHSIGMA 10 /* length in linear voxel dimensions */
#define MAXSEARCHLENGTH (3*SEARCHSIGMA) #define MAXSEARCHLENGTH (3*SEARCHSIGMA)
#define FLOWDEBUG /*#define FLOWDEBUG*/
{ {
int d, T, iters; int d, T, iters;
...@@ -853,7 +848,7 @@ double tmpd; ...@@ -853,7 +848,7 @@ double tmpd;
#endif #endif
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++) 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) if (tmpim[z*xysize+y*im.x+x] > origthresh)
{ {
int xxx=XX[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 yyy=YY[z*xysize+y*im.x+x];
...@@ -863,6 +858,8 @@ double tmpd; ...@@ -863,6 +858,8 @@ double tmpd;
short maxvalX=0, maxvalY=0, maxvalZ=0; short maxvalX=0, maxvalY=0, maxvalZ=0;
float exponentfactor = -0.5 * (xxx*xxx+yyy*yyy+zzz*zzz) / (float)(SEARCHSIGMA*SEARCHSIGMA); float exponentfactor = -0.5 * (xxx*xxx+yyy*yyy+zzz*zzz) / (float)(SEARCHSIGMA*SEARCHSIGMA);
/* need to put proper FOV bounds checks in both below searches !!! */
if (lowercingulum.i[z*xysize+y*im.x+x] == 0) if (lowercingulum.i[z*xysize+y*im.x+x] == 0)
/* {{{ search perp to sheet */ /* {{{ search perp to sheet */
...@@ -876,10 +873,10 @@ double tmpd; ...@@ -876,10 +873,10 @@ double tmpd;
int D=d; int D=d;
if (iters==1) 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) if (distancemap.i[(z+zzz*D)*xysize+(y+yyy*D)*im.x+x+xxx*D]>=distance)
{ {
float distanceweight = exp(d * d * exponentfactor); float distanceweight = exp(d * d * exponentfactor);
distance=distancemap.i[t*xyzsize+(z+zzz*D)*xysize+(y+yyy*D)*im.x+x+xxx*D]; distance=distancemap.i[(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) 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=thedata.i[T*xyzsize+(z+zzz*D)*xysize+(y+yyy*D)*im.x+x+xxx*D];
...@@ -900,24 +897,39 @@ double tmpd; ...@@ -900,24 +897,39 @@ double tmpd;
/* {{{ search all around tube */ /* {{{ search all around tube */
{ {
for(yyy=-5; yyy<=5; yyy++) for(xxx=-5; xxx<=5; xxx++) for(yyy=-MAXSEARCHLENGTH; yyy<=MAXSEARCHLENGTH; yyy++) for(xxx=-MAXSEARCHLENGTH; xxx<=MAXSEARCHLENGTH; xxx++)
{ {
float distanceweight = exp(-0.5 * (xxx*xxx+yyy*yyy) / (float)(SEARCHSIGMA*SEARCHSIGMA) ); float distanceweight = exp(-0.5 * (xxx*xxx+yyy*yyy) / (float)(SEARCHSIGMA*SEARCHSIGMA) );
float r=sqrt((float)(xxx*xxx+yyy*yyy)); /* now generate xxx2,yyy2 = next voxel outwards */ float r=sqrt((float)(xxx*xxx+yyy*yyy));
int xxx2=FTOI(xxx/r * (r+1));
int yyy2=FTOI(yyy/r * (r+1));
if ( ( (r==0) || ( distancemap.i[z*xysize+(y+yyy)*im.x+x+xxx] >= distancemap.i[z*xysize+(y+yyy2)*im.x+x+xxx2] ) ) && if (r>0)
( distanceweight * thedata.i[T*xyzsize+z*xysize+(y+yyy)*im.x+x+xxx] > maxval_weighted ) )
{ {
maxval=thedata.i[T*xyzsize+z*xysize+(y+yyy)*im.x+x+xxx]; float rr;
maxval_weighted=maxval*distanceweight; int allok=1;
maxvalX=xxx;
maxvalY=yyy; for(rr=1; rr<=r+0.1; rr++) /* search outwards from centre to current voxel - test that distancemap always increasing */
maxvalZ=0; {
int xxx1=FTOI(rr*xxx/r);
int yyy1=FTOI(rr*yyy/r);
int xxx2=FTOI((rr+1)*xxx/r);
int yyy2=FTOI((rr+1)*yyy/r);
if ( distancemap.i[z*xysize+(y+yyy1)*im.x+x+xxx1] > distancemap.i[z*xysize+(y+yyy2)*im.x+x+xxx2] )
allok=0;
}
if ( allok &&
( distanceweight * thedata.i[T*xyzsize+z*xysize+(y+yyy)*im.x+x+xxx] > maxval_weighted ) )
{
maxval=thedata.i[T*xyzsize+z*xysize+(y+yyy)*im.x+x+xxx];
maxval_weighted=maxval*distanceweight;
maxvalX=xxx;
maxvalY=yyy;
maxvalZ=0;
}
} }
} }
} }
/* }}} */ /* }}} */
...@@ -959,8 +971,7 @@ double tmpd; ...@@ -959,8 +971,7 @@ double tmpd;
/* }}} */ /* }}} */
free(X); free(Y); free(Z); free(XX); free(YY); free(ZZ); free(XX); free(YY); free(ZZ);
}
free(in[0]); free(in[0]);
in[0]=tmpim; in[0]=tmpim;
......
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