|
Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
edt_filtervol.c File Reference#include "mrilib.h"
#include "parser.h"
Go to the source code of this file.
|
Defines |
#define | GOODVOX(ijk) (fmask==NULL || fmask[ijk]!=0) |
#define | BADVOX(ijk) (fmask!=NULL && fmask[ijk]==0) |
#define | II 8 |
#define | JJ 9 |
#define | KK 10 |
#define | RR 17 |
#define | XX 23 |
#define | YY 24 |
#define | ZZ 25 |
Functions |
void | EDIT_filter_volume (int nx, int ny, int nz, float dx, float dy, float dz, int fim_type, void *vfim, int filter_opt, float filter_rmm, byte *fmask, char *fexpr) |
void | EDIT_aver_fvol (int nx, int ny, int nz, float dx, float dy, float dz, float *fim, float rmm) |
Define Documentation
#define BADVOX |
( |
ijk |
|
) |
(fmask!=NULL && fmask[ijk]==0)
|
|
#define GOODVOX |
( |
ijk |
|
) |
(fmask==NULL || fmask[ijk]!=0)
|
|
Function Documentation
void EDIT_aver_fvol |
( |
int |
nx, |
|
|
int |
ny, |
|
|
int |
nz, |
|
|
float |
dx, |
|
|
float |
dy, |
|
|
float |
dz, |
|
|
float * |
fim, |
|
|
float |
rmm |
|
) |
|
|
|
Definition at line 391 of file edt_filtervol.c.
References abs, ENTRY, EXIT, fim, free, MCW_cluster::i, i, MCW_cluster::j, MCW_cluster::k, KILL_CLUSTER, malloc, MAX, MCW_build_mask(), MCW_cluster::num_pt, and nz.
Referenced by EDIT_filter_volume().
00393 {
00394 MCW_cluster * mask ;
00395 int i, j, k , ij , ii ;
00396 int jk,jkadd , nxadd,nyadd,nzadd , nxyz_add , mnum ;
00397 float * ffim ;
00398 int * madd ;
00399 float fac , sum ;
00400
00401 ENTRY("EDIT_aver_fvol") ;
00402
00403
00404
00405 mask = MCW_build_mask(nx,ny,nz, dx,dy,dz, rmm) ;
00406 if( mask == NULL || mask->num_pt < 2 ){
00407 fprintf(stderr,"Warning: EDIT_aver_volume has no effect.\n") ;
00408 EXRETURN ;
00409 }
00410 mnum = mask->num_pt ;
00411
00412
00413
00414 #if 1
00415 nxadd = nyadd = nzadd = 1 ;
00416 for( ii=0 ; ii < mnum ; ii++ ){
00417 i = abs((int)mask->i[ii]) ; nxadd = MAX(i,nxadd) ;
00418 j = abs((int)mask->j[ii]) ; nyadd = MAX(j,nyadd) ;
00419 k = abs((int)mask->k[ii]) ; nzadd = MAX(k,nzadd) ;
00420 }
00421 #else
00422 nxadd = (int)(rmm/dx) ;
00423 nyadd = (int)(rmm/dy) ;
00424 nzadd = (int)(rmm/dz) ;
00425 #endif
00426
00427 nxyz_add = (nx+2*nxadd) * (ny+2*nyadd) * (nz+2*nzadd) ;
00428
00429 ffim = (float *) malloc(sizeof(float) * nxyz_add) ;
00430 if(ffim == NULL){
00431 fprintf(stderr,"*** EDIT_aver_volume can't malloc workspace!\n") ;
00432 fprintf(stderr,"nx=%d ny=%d nz=%d nxadd=%d nyadd=%d nzadd=%d\n",
00433 nx,ny,nz , nxadd,nyadd,nzadd ) ;
00434 EXIT(1) ;
00435 }
00436 for( i=0 ; i < nxyz_add ; i++ ) ffim[i] = 0.0 ;
00437
00438 madd = (int *) malloc( sizeof(int) * (mnum+1) ) ;
00439 if( madd == NULL ){
00440 fprintf(stderr,"*** EDIT_aver_volume can't malloc workspace!\n") ;
00441 EXIT(1) ;
00442 }
00443 madd[0] = 0 ;
00444 for( ii=0 ; ii < mnum ; ii++ ){
00445 madd[ii+1] = mask->i[ii] +
00446 mask->j[ii] * (nx+2*nxadd) +
00447 mask->k[ii] * ((nx+2*nxadd)*(ny+2*nyadd)) ;
00448 }
00449 mnum++ ; fac = 1.0 / mnum ;
00450
00451 KILL_CLUSTER(mask) ;
00452
00453
00454
00455 for( k=0 ; k < nz ; k++ ){
00456 for( j=0 ; j < ny ; j++ ){
00457 jkadd = j * (nx+2*nxadd) + k * ((nx+2*nxadd)*(ny+2*nyadd)) ;
00458 jk = j * nx + k * (nx * ny) ;
00459 for( i=0 ; i < nx ; i++ ) ffim[i+jkadd] = fim[i+jk] ;
00460 }
00461 }
00462
00463
00464
00465 for( k=0 ; k < nz ; k++ ){
00466 for( j=0 ; j < ny ; j++ ){
00467 jkadd = j * (nx+2*nxadd) + k * ((nx+2*nxadd)*(ny+2*nyadd)) ;
00468 jk = j * nx + k * (nx * ny) ;
00469 for( i=0 ; i < nx ; i++ ){
00470 sum = 0.0 ; ij = i+jkadd ;
00471 for( ii=0 ; ii < mnum ; ii++ ) sum += ffim[ij+madd[ii]] ;
00472 fim[i+jk] = fac * sum ;
00473 }
00474 }
00475 }
00476
00477 free(ffim); free(madd);
00478 EXRETURN;
00479 }
|
void EDIT_filter_volume |
( |
int |
nx, |
|
|
int |
ny, |
|
|
int |
nz, |
|
|
float |
dx, |
|
|
float |
dy, |
|
|
float |
dz, |
|
|
int |
fim_type, |
|
|
void * |
vfim, |
|
|
int |
filter_opt, |
|
|
float |
filter_rmm, |
|
|
byte * |
fmask, |
|
|
char * |
fexpr |
|
) |
|
|
|
Definition at line 38 of file edt_filtervol.c.
References amax, EDIT_aver_fvol(), EDIT_coerce_autoscale(), EDIT_coerce_type(), ENTRY, EXIT, FCFLAG_AMAX, FCFLAG_AVER, FCFLAG_EXPR, FCFLAG_MAX, FCFLAG_MEAN, FCFLAG_NZMEAN, FCFLAG_ONE_STEP, FCFLAG_SMAX, FCFLAG_WINSOR, free, MCW_cluster::i, i, MCW_cluster::j, MCW_cluster::k, KILL_CLUSTER, malloc, MCW_build_mask(), MCW_cluster::num_pt, nz, PARSER_evaluate_one(), PARSER_generate_code(), qsort_float(), and THREE_TO_IJK.
Referenced by EDIT_one_dataset().
00041 {
00042 MCW_cluster * mask;
00043 int nxy, nxyz;
00044 int mnum;
00045 int i, j, k, ii, jj, kk;
00046 int ijkvox, ijkma, jma;
00047 float * ffim, * ffim_out;
00048
00049 float
00050 mag,
00051 sum,
00052 sumnz,
00053 mean,
00054 max,
00055 amax,
00056 smax;
00057 int
00058 npts, nznpts;
00059
00060 float wtsum ;
00061 float * wt=NULL ;
00062 PARSER_code * pcode ;
00063
00064 int nw , nnw , iw ;
00065 float * sw=NULL , vw ;
00066
00067 nxy = nx*ny; nxyz = nxy*nz;
00068
00069 #define GOODVOX(ijk) (fmask==NULL || fmask[ijk]!=0)
00070 #define BADVOX(ijk) (fmask!=NULL && fmask[ijk]==0)
00071
00072 ENTRY("EDIT_filter_volume") ;
00073
00074
00075
00076 if( fmask != NULL && filter_opt == FCFLAG_AVER ) filter_opt = FCFLAG_MEAN ;
00077
00078
00079
00080 if( filter_opt == FCFLAG_AVER ){
00081 if( fim_type != MRI_float ){
00082 ffim = (float *) malloc (sizeof(float) * nxyz);
00083 if( ffim == NULL ){
00084 fprintf(stderr,"EDIT_filter_volume: no workspace!\n") ;
00085 EXIT(1) ;
00086 }
00087 EDIT_coerce_type (nxyz, fim_type, vfim, MRI_float, ffim);
00088 } else {
00089 ffim = (float *) vfim ;
00090 }
00091 EDIT_aver_fvol(nx,ny,nz,dx,dy,dz,ffim,filter_rmm) ;
00092 if( ffim != vfim ){
00093 EDIT_coerce_autoscale(nxyz, MRI_float, ffim, fim_type, vfim);
00094 free(ffim) ;
00095 }
00096 EXRETURN ;
00097 }
00098
00099
00100
00101
00102
00103 mask = MCW_build_mask (nx, ny, nz, dx, dy, dz, filter_rmm);
00104 if (mask == NULL)
00105 {
00106 fprintf (stderr, "Warning: Filter option has no effect. \n");
00107 EXRETURN;
00108 }
00109 mnum = mask->num_pt;
00110
00111
00112
00113 if( filter_opt == FCFLAG_EXPR ){
00114 double atoz[26] ;
00115
00116 if( fexpr == NULL ){
00117 fprintf(stderr,"*** EDIT_filter_volume: no fexpr for FCFLAG_EXPR!\n");
00118 EXIT(1) ;
00119 }
00120
00121 pcode = PARSER_generate_code( fexpr ) ;
00122 if( pcode == NULL ){
00123 fprintf(stderr,"*** EDIT_filter_volume: illegal fexpr!\n"); EXIT(1);
00124 }
00125
00126 wt = (float *) malloc(sizeof(float)*(mnum+1)) ;
00127
00128 #define II 8
00129 #define JJ 9
00130 #define KK 10
00131 #define RR 17
00132 #define XX 23
00133 #define YY 24
00134 #define ZZ 25
00135
00136 for( ii=0 ; ii < 26 ; ii++ ) atoz[ii] = 0.0 ;
00137
00138 wt[0] = PARSER_evaluate_one( pcode , atoz ) ;
00139
00140 for (jma = 0; jma < mnum; jma++){
00141 atoz[II] = mask->i[jma] ;
00142 atoz[JJ] = mask->j[jma] ;
00143 atoz[KK] = mask->k[jma] ;
00144 atoz[XX] = atoz[II] * dx ;
00145 atoz[YY] = atoz[JJ] * dy ;
00146 atoz[ZZ] = atoz[KK] * dz ;
00147 atoz[RR] = sqrt(atoz[XX]*atoz[XX] + atoz[YY]*atoz[YY] + atoz[ZZ]*atoz[ZZ]) ;
00148 wt[jma+1] = PARSER_evaluate_one( pcode , atoz ) ;
00149 }
00150 free(pcode) ;
00151 }
00152
00153
00154
00155 if( filter_opt > FCFLAG_WINSOR &&
00156 filter_opt < FCFLAG_WINSOR+FCFLAG_ONE_STEP ){
00157
00158 static int first=1 ;
00159
00160 nw = filter_opt - FCFLAG_WINSOR ;
00161 nnw = mnum - nw ;
00162 filter_opt = FCFLAG_WINSOR ;
00163
00164 fprintf(stderr,"++ Winsor filter: N=%d nw=%d\n",mnum+1,nw) ;
00165 if( first || nnw < nw ){
00166 first = 0 ;
00167 if( nnw < nw ){
00168 fprintf(stderr,"** Illegal Winsor parameters - skipping!\n") ;
00169 EXRETURN ;
00170 }
00171 }
00172
00173 sw = (float *) malloc(sizeof(float)*(mnum+1)) ;
00174
00175 fmask = NULL ;
00176 }
00177
00178
00179
00180 ffim = (float *) malloc (sizeof(float) * nxyz);
00181 if (ffim == NULL)
00182 {
00183 fprintf (stderr, "\n Error: cannot allocate filter workspace! \n");
00184 EXIT(1);
00185 }
00186 ffim_out = (float *) malloc (sizeof(float) * nxyz);
00187 if (ffim_out == NULL)
00188 {
00189 fprintf (stderr, "\n Error: cannot allocate filter workspace! \n");
00190 EXIT(1);
00191 }
00192
00193
00194
00195 EDIT_coerce_type (nxyz, fim_type, vfim, MRI_float, ffim);
00196
00197
00198
00199 for (k = 0; k < nz; k++)
00200 {
00201 for (j = 0; j < ny; j++)
00202 {
00203 for (i = 0; i < nx; i++)
00204 {
00205
00206
00207 ijkvox = THREE_TO_IJK (i, j, k, nx, nxy);
00208 npts = nznpts = 0 ;
00209 sum = sumnz = max = amax = smax = wtsum = 0.0 ;
00210 if( GOODVOX(ijkvox) ){
00211 mag = ffim[ijkvox];
00212 switch (filter_opt)
00213 {
00214 case FCFLAG_MEAN:
00215 sum = mag; npts = 1; break;
00216 case FCFLAG_NZMEAN:
00217 if (mag != 0.0)
00218 {sumnz = mag; nznpts = 1;}
00219 break;
00220 case FCFLAG_MAX:
00221 max = mag; npts = 1 ; break;
00222 case FCFLAG_AMAX:
00223 amax = fabs(mag); npts = 1 ; break;
00224 case FCFLAG_SMAX:
00225 smax = mag; npts = 1 ; break;
00226 case FCFLAG_EXPR:
00227 sum = wt[0]*mag ; wtsum = wt[0] ; npts = 1 ; break ;
00228
00229 case FCFLAG_WINSOR:
00230 for( iw=0 ; iw <= mnum ; iw++ ) sw[iw] = mag ;
00231 vw = mag ; break ;
00232 }
00233 }
00234
00235
00236
00237 switch (filter_opt)
00238 {
00239 case FCFLAG_MEAN:
00240 for (jma = 0; jma < mnum; jma++)
00241 {
00242 ii = i + mask->i[jma];
00243 jj = j + mask->j[jma];
00244 kk = k + mask->k[jma];
00245 if (ii < 0 || jj < 0 || kk < 0 ||
00246 ii >= nx || jj >= ny || kk >= nz) continue;
00247 ijkma = THREE_TO_IJK (ii, jj, kk, nx, nxy);
00248 if( BADVOX(ijkma) ) continue ;
00249 mag = ffim[ijkma];
00250 sum += mag; npts++;
00251 }
00252 break;
00253 case FCFLAG_NZMEAN:
00254 for (jma = 0; jma < mnum; jma++)
00255 {
00256 ii = i + mask->i[jma];
00257 jj = j + mask->j[jma];
00258 kk = k + mask->k[jma];
00259 if (ii < 0 || jj < 0 || kk < 0 ||
00260 ii >= nx || jj >= ny || kk >= nz) continue;
00261 ijkma = THREE_TO_IJK (ii, jj, kk, nx, nxy);
00262 if( BADVOX(ijkma) ) continue ;
00263 mag = ffim[ijkma];
00264 if (mag != 0.0) {sumnz += mag; nznpts++;}
00265 }
00266 break;
00267 case FCFLAG_EXPR:
00268 for (jma = 0; jma < mnum; jma++)
00269 {
00270 ii = i + mask->i[jma];
00271 jj = j + mask->j[jma];
00272 kk = k + mask->k[jma];
00273 if (ii < 0 || jj < 0 || kk < 0 ||
00274 ii >= nx || jj >= ny || kk >= nz) continue;
00275 ijkma = THREE_TO_IJK (ii, jj, kk, nx, nxy);
00276 if( BADVOX(ijkma) ) continue ;
00277 mag = ffim[ijkma];
00278 sum += wt[jma+1]*mag; npts++; wtsum += wt[jma+1] ;
00279 }
00280 break ;
00281 case FCFLAG_MAX:
00282 for (jma = 0; jma < mnum; jma++)
00283 {
00284 ii = i + mask->i[jma];
00285 jj = j + mask->j[jma];
00286 kk = k + mask->k[jma];
00287 if (ii < 0 || jj < 0 || kk < 0 ||
00288 ii >= nx || jj >= ny || kk >= nz) continue;
00289 ijkma = THREE_TO_IJK (ii, jj, kk, nx, nxy);
00290 if( BADVOX(ijkma) ) continue ;
00291 mag = ffim[ijkma];
00292 if (npts == 0 || mag > max) max = mag;
00293 npts++ ;
00294 }
00295 break;
00296 case FCFLAG_AMAX:
00297 for (jma = 0; jma < mnum; jma++)
00298 {
00299 ii = i + mask->i[jma];
00300 jj = j + mask->j[jma];
00301 kk = k + mask->k[jma];
00302 if (ii < 0 || jj < 0 || kk < 0 ||
00303 ii >= nx || jj >= ny || kk >= nz) continue;
00304 ijkma = THREE_TO_IJK (ii, jj, kk, nx, nxy);
00305 if( BADVOX(ijkma) ) continue ;
00306 mag = ffim[ijkma];
00307 if (npts == 0 || fabs(mag) > amax) amax = fabs(mag);
00308 npts++ ;
00309 }
00310 break;
00311 case FCFLAG_SMAX:
00312 for (jma = 0; jma < mnum; jma++)
00313 {
00314 ii = i + mask->i[jma];
00315 jj = j + mask->j[jma];
00316 kk = k + mask->k[jma];
00317 if (ii < 0 || jj < 0 || kk < 0 ||
00318 ii >= nx || jj >= ny || kk >= nz) continue;
00319 ijkma = THREE_TO_IJK (ii, jj, kk, nx, nxy);
00320 if( BADVOX(ijkma) ) continue ;
00321 mag = ffim[ijkma];
00322 if (npts == 0 || fabs(mag) > fabs(smax)) smax = mag;
00323 npts++ ;
00324 }
00325 break;
00326 case FCFLAG_WINSOR:
00327 for (jma = 0; jma < mnum; jma++)
00328 {
00329 ii = i + mask->i[jma];
00330 jj = j + mask->j[jma];
00331 kk = k + mask->k[jma];
00332 if (ii < 0 || jj < 0 || kk < 0 ||
00333 ii >= nx || jj >= ny || kk >= nz) continue;
00334 ijkma = THREE_TO_IJK (ii, jj, kk, nx, nxy);
00335 sw[jma+1] = ffim[ijkma];
00336 }
00337 break;
00338
00339 default: break;
00340 }
00341
00342
00343
00344 switch (filter_opt)
00345 {
00346 case FCFLAG_MEAN:
00347 ffim_out[ijkvox] = (npts > 0) ? sum/npts : 0.0 ; break;
00348
00349 case FCFLAG_NZMEAN:
00350 ffim_out[ijkvox] = (nznpts > 0) ? sumnz/nznpts : 0.0 ; break ;
00351
00352 case FCFLAG_EXPR:
00353 ffim_out[ijkvox] = (wtsum != 0.0) ? sum/wtsum : 0.0 ; break ;
00354
00355 case FCFLAG_MAX: ffim_out[ijkvox] = max ; break;
00356 case FCFLAG_AMAX: ffim_out[ijkvox] = amax; break;
00357 case FCFLAG_SMAX: ffim_out[ijkvox] = smax; break;
00358
00359 case FCFLAG_WINSOR:
00360 qsort_float( mnum+1 , sw ) ;
00361 ffim_out[ijkvox] = (vw < sw[nw]) ? sw[nw]
00362 : (vw > sw[nnw]) ? sw[nnw]
00363 : vw ;
00364 break ;
00365
00366 default: break;
00367 }
00368 }
00369 }
00370 }
00371
00372
00373 EDIT_coerce_autoscale(nxyz, MRI_float, ffim_out, fim_type, vfim);
00374
00375
00376 KILL_CLUSTER (mask);
00377 free (ffim);
00378 free (ffim_out);
00379 if( wt != NULL ) free(wt) ;
00380 if( sw != NULL ) free(sw) ;
00381
00382 EXRETURN;
00383 }
|
|