Doxygen Source Code Documentation
edt_blur.c File Reference
#include "mrilib.h"
Go to the source code of this file.
Defines | |
#define | GET_AS_BIG(name, type, dim) |
Functions | |
void | EDIT_blur_volume (int nx, int ny, int nz, float dx, float dy, float dz, int ftype, void *vfim, float sigma) |
void | EDIT_blur_volume_3d (int nx, int ny, int nz, float dx, float dy, float dz, int ftype, void *vfim, float sigmax, float sigmay, float sigmaz) |
Define Documentation
|
Value: do{ if( (dim) > name ## _size ){ \ if( name != NULL ) free(name) ; \ name = (type *) malloc( sizeof(type) * (dim) ) ; \ if( name == NULL ){ \ fprintf(stderr,"\n*** cannot malloc EDIT workspace!\n") ; \ EXIT(1) ; } \ name ## _size = (dim) ; } \ break ; } while(1) Definition at line 13 of file edt_blur.c. Referenced by EDIT_blur_volume_3d(), mri_filt_fft(), and shifter(). |
Function Documentation
|
Definition at line 24 of file edt_blur.c. References EDIT_blur_volume_3d(), and nz. Referenced by EDIT_one_dataset().
00027 { 00028 EDIT_blur_volume_3d (nx, ny, nz, dx, dy, dz, ftype, vfim, 00029 sigma, sigma, sigma); 00030 } |
|
Definition at line 40 of file edt_blur.c. References base, csfft_cox(), csfft_nextup_one35(), ENTRY, GET_AS_BIG, complex::i, nz, complex::r, and STATUS. Referenced by EDIT_blur_volume(), gaussian_filter(), mri_3dalign_setup(), mri_warp3d_align_one(), and mri_warp3D_align_setup().
00044 { 00045 int jj,kk , nxy , base,nby2 ; 00046 float dk , aa , k , fac ; 00047 register int ii , nup ; 00048 00049 static int cx_size = 0 ; /* workspaces */ 00050 static int gg_size = 0 ; 00051 static complex * cx = NULL ; 00052 static float * gg = NULL ; 00053 00054 byte * bfim = NULL ; 00055 short * sfim = NULL ; 00056 float * ffim = NULL ; 00057 complex * cfim = NULL ; 00058 00059 float fbot,ftop ; /* 10 Jan 2003 */ 00060 int nxyz ; 00061 00062 /*** initialize ***/ 00063 00064 ENTRY("EDIT_blur_volume_3d") ; 00065 00066 if( vfim == NULL || 00067 sigmax <= 0.0 || sigmay <= 0.0 || sigmaz <= 0.0) EXRETURN ; 00068 00069 switch( ftype ){ 00070 default: EXRETURN ; 00071 case MRI_short: sfim = (short *) vfim ; break ; 00072 case MRI_float: ffim = (float *) vfim ; break ; 00073 case MRI_byte: bfim = (byte *) vfim ; break ; 00074 case MRI_complex: cfim = (complex *) vfim ; break ; 00075 } 00076 nxy = nx * ny ; nxyz = nxy * nz ; 00077 00078 /*** 10 Jan 2003: find bot and top of data input */ 00079 00080 switch( ftype ){ 00081 default: 00082 fbot = ftop = 0.0 ; 00083 break ; 00084 00085 case MRI_short: 00086 fbot = ftop = sfim[0] ; 00087 for( ii=1 ; ii < nxyz ; ii++ ) 00088 if( sfim[ii] < fbot ) fbot = sfim[ii] ; 00089 else if( sfim[ii] > ftop ) ftop = sfim[ii] ; 00090 break ; 00091 00092 case MRI_float: 00093 fbot = ftop = ffim[0] ; 00094 for( ii=1 ; ii < nxyz ; ii++ ) 00095 if( ffim[ii] < fbot ) fbot = ffim[ii] ; 00096 else if( ffim[ii] > ftop ) ftop = ffim[ii] ; 00097 break ; 00098 00099 case MRI_byte: 00100 fbot = ftop = bfim[0] ; 00101 for( ii=1 ; ii < nxyz ; ii++ ) 00102 if( bfim[ii] < fbot ) fbot = bfim[ii] ; 00103 else if( bfim[ii] > ftop ) ftop = bfim[ii] ; 00104 break ; 00105 } 00106 00107 /*** do x-direction ***/ 00108 00109 STATUS("start x FFTs") ; 00110 00111 aa = sigmax * sigmax * 0.5 ; 00112 nup = nx + (int)(3.0 * sigmax / dx) ; /* min FFT length */ 00113 #if 0 00114 ii = 2 ; while( ii < nup ){ ii *= 2 ; } /* next power of 2 larger */ 00115 #else 00116 ii = csfft_nextup_one35(nup) ; 00117 #endif 00118 nup = ii ; nby2 = nup / 2 ; 00119 00120 GET_AS_BIG(cx,complex,nup) ; GET_AS_BIG(gg,float,nup) ; 00121 00122 dk = (2.0*PI) / (nup * dx) ; 00123 fac = 1.0 / nup ; 00124 gg[0] = fac ; 00125 for( ii=1 ; ii<=nby2 ; ii++ ){ k=ii*dk; gg[nup-ii]=gg[ii]=fac*exp(-aa*k*k); } 00126 00127 /** July 20: double up on FFTs **/ 00128 /** Feb 09: extend to other data types besides shorts; 00129 doubling up does not apply to complex data! **/ 00130 00131 switch( ftype ){ 00132 case MRI_short:{ 00133 register short * qfim ; 00134 for( kk=0 ; kk < nz ; kk++ ){ 00135 for( jj=0 ; jj < ny ; jj+=2 ){ 00136 base = jj*nx + kk*nxy ; 00137 qfim = sfim + base ; 00138 if( jj == ny-1 ) 00139 for( ii=0 ; ii<nx ; ii++){ cx[ii].r = qfim[ii] ; cx[ii].i = 0.0 ; } 00140 else 00141 for( ii=0 ; ii<nx ; ii++){ cx[ii].r = qfim[ii] ; cx[ii].i = qfim[ii+nx] ; } 00142 for( ii=nx; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; } 00143 csfft_cox( -1 , nup , cx ) ; 00144 for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; } 00145 csfft_cox( 1 , nup , cx ) ; 00146 if( jj == ny-1 ) 00147 for( ii=0 ; ii<nx ; ii++){ qfim[ii] = cx[ii].r ; } 00148 else 00149 for( ii=0 ; ii<nx ; ii++){ qfim[ii] = cx[ii].r ; qfim[ii+nx] = cx[ii].i ; } 00150 } 00151 } 00152 } 00153 break ; 00154 00155 case MRI_float:{ 00156 register float * qfim ; 00157 for( kk=0 ; kk < nz ; kk++ ){ 00158 for( jj=0 ; jj < ny ; jj+=2 ){ 00159 base = jj*nx + kk*nxy ; 00160 qfim = ffim + base ; 00161 if( jj == ny-1 ) 00162 for( ii=0 ; ii<nx ; ii++){ cx[ii].r = qfim[ii] ; cx[ii].i = 0.0 ; } 00163 else 00164 for( ii=0 ; ii<nx ; ii++){ cx[ii].r = qfim[ii] ; cx[ii].i = qfim[ii+nx] ; } 00165 for( ii=nx; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; } 00166 csfft_cox( -1 , nup , cx ) ; 00167 for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; } 00168 csfft_cox( 1 , nup , cx ) ; 00169 if( jj == ny-1 ) 00170 for( ii=0 ; ii<nx ; ii++){ qfim[ii] = cx[ii].r ; } 00171 else 00172 for( ii=0 ; ii<nx ; ii++){ qfim[ii] = cx[ii].r ; qfim[ii+nx] = cx[ii].i ; } 00173 } 00174 } 00175 } 00176 break ; 00177 00178 case MRI_byte:{ 00179 register byte * qfim ; 00180 for( kk=0 ; kk < nz ; kk++ ){ 00181 for( jj=0 ; jj < ny ; jj+=2 ){ 00182 base = jj*nx + kk*nxy ; 00183 qfim = bfim + base ; 00184 if( jj == ny-1 ) 00185 for( ii=0 ; ii<nx ; ii++){ cx[ii].r = qfim[ii] ; cx[ii].i = 0.0 ; } 00186 else 00187 for( ii=0 ; ii<nx ; ii++){ cx[ii].r = qfim[ii] ; cx[ii].i = qfim[ii+nx] ; } 00188 for( ii=nx; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; } 00189 csfft_cox( -1 , nup , cx ) ; 00190 for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; } 00191 csfft_cox( 1 , nup , cx ) ; 00192 if( jj == ny-1 ) 00193 for( ii=0 ; ii<nx ; ii++){ qfim[ii] = cx[ii].r ; } 00194 else 00195 for( ii=0 ; ii<nx ; ii++){ qfim[ii] = cx[ii].r ; qfim[ii+nx] = cx[ii].i ; } 00196 } 00197 } 00198 } 00199 break ; 00200 00201 case MRI_complex:{ 00202 register complex * qfim ; 00203 for( kk=0 ; kk < nz ; kk++ ){ 00204 for( jj=0 ; jj < ny ; jj++ ){ 00205 base = jj*nx + kk*nxy ; 00206 qfim = cfim + base ; 00207 for( ii=0 ; ii<nx ; ii++) { cx[ii] = qfim[ii] ; } 00208 for( ii=nx; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; } 00209 csfft_cox( -1 , nup , cx ) ; 00210 for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; } 00211 csfft_cox( 1 , nup , cx ) ; 00212 for( ii=0 ; ii<nx ; ii++){ qfim[ii] = cx[ii] ; } 00213 } 00214 } 00215 } 00216 break ; 00217 } 00218 00219 /*** do y-direction ***/ 00220 00221 STATUS("start y FFTs") ; 00222 00223 aa = sigmay * sigmay * 0.5 ; 00224 nup = ny + (int)(3.0 * sigmay / dy) ; /* min FFT length */ 00225 #if 0 00226 ii = 2 ; while( ii < nup ){ ii *= 2 ; } /* next power of 2 larger */ 00227 #else 00228 ii = csfft_nextup_one35(nup) ; 00229 #endif 00230 nup = ii ; nby2 = nup / 2 ; 00231 00232 GET_AS_BIG(cx,complex,nup) ; GET_AS_BIG(gg,float,nup) ; 00233 00234 dk = (2.0*PI) / (nup * dy) ; 00235 fac = 1.0 / nup ; 00236 gg[0] = fac ; 00237 for( ii=1 ; ii<=nby2 ; ii++ ){ k=ii*dk; gg[nup-ii]=gg[ii]=fac*exp(-aa*k*k); } 00238 00239 switch( ftype ){ 00240 case MRI_short:{ 00241 register short * qfim ; 00242 for( kk=0 ; kk < nz ; kk++ ){ 00243 for( jj=0 ; jj < nx ; jj+=2 ){ 00244 base = jj + kk*nxy ; 00245 qfim = sfim + base ; 00246 if( jj == nx-1 ) 00247 for( ii=0 ; ii<ny ; ii++){ cx[ii].r = qfim[ii*nx] ; cx[ii].i = 0.0 ; } 00248 else 00249 for( ii=0 ; ii<ny ; ii++){ cx[ii].r = qfim[ii*nx] ; cx[ii].i = qfim[ii*nx+1] ; } 00250 for( ii=ny; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; } 00251 csfft_cox( -1 , nup , cx ) ; 00252 for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; } 00253 csfft_cox( 1 , nup , cx ) ; 00254 if( jj == nx-1 ) 00255 for( ii=0 ; ii<ny ; ii++){ qfim[ii*nx] = cx[ii].r ; } 00256 else 00257 for( ii=0 ; ii<ny ; ii++){ qfim[ii*nx] = cx[ii].r ; qfim[ii*nx+1] = cx[ii].i ; } 00258 } 00259 } 00260 } 00261 break ; 00262 00263 case MRI_byte:{ 00264 register byte * qfim ; 00265 for( kk=0 ; kk < nz ; kk++ ){ 00266 for( jj=0 ; jj < nx ; jj+=2 ){ 00267 base = jj + kk*nxy ; 00268 qfim = bfim + base ; 00269 if( jj == nx-1 ) 00270 for( ii=0 ; ii<ny ; ii++){ cx[ii].r = qfim[ii*nx] ; cx[ii].i = 0.0 ; } 00271 else 00272 for( ii=0 ; ii<ny ; ii++){ cx[ii].r = qfim[ii*nx] ; cx[ii].i = qfim[ii*nx+1] ; } 00273 for( ii=ny; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; } 00274 csfft_cox( -1 , nup , cx ) ; 00275 for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; } 00276 csfft_cox( 1 , nup , cx ) ; 00277 if( jj == nx-1 ) 00278 for( ii=0 ; ii<ny ; ii++){ qfim[ii*nx] = cx[ii].r ; } 00279 else 00280 for( ii=0 ; ii<ny ; ii++){ qfim[ii*nx] = cx[ii].r ; qfim[ii*nx+1] = cx[ii].i ; } 00281 } 00282 } 00283 } 00284 break ; 00285 00286 case MRI_float:{ 00287 register float * qfim ; 00288 for( kk=0 ; kk < nz ; kk++ ){ 00289 for( jj=0 ; jj < nx ; jj+=2 ){ 00290 base = jj + kk*nxy ; 00291 qfim = ffim + base ; 00292 if( jj == nx-1 ) 00293 for( ii=0 ; ii<ny ; ii++){ cx[ii].r = qfim[ii*nx] ; cx[ii].i = 0.0 ; } 00294 else 00295 for( ii=0 ; ii<ny ; ii++){ cx[ii].r = qfim[ii*nx] ; cx[ii].i = qfim[ii*nx+1] ; } 00296 for( ii=ny; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; } 00297 csfft_cox( -1 , nup , cx ) ; 00298 for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; } 00299 csfft_cox( 1 , nup , cx ) ; 00300 if( jj == nx-1 ) 00301 for( ii=0 ; ii<ny ; ii++){ qfim[ii*nx] = cx[ii].r ; } 00302 else 00303 for( ii=0 ; ii<ny ; ii++){ qfim[ii*nx] = cx[ii].r ; qfim[ii*nx+1] = cx[ii].i ; } 00304 } 00305 } 00306 } 00307 break ; 00308 00309 case MRI_complex:{ 00310 register complex * qfim ; 00311 for( kk=0 ; kk < nz ; kk++ ){ 00312 for( jj=0 ; jj < nx ; jj++ ){ 00313 base = jj + kk*nxy ; 00314 qfim = cfim + base ; 00315 for( ii=0 ; ii<ny ; ii++){ cx[ii] = qfim[ii*nx] ; } 00316 for( ii=ny; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; } 00317 csfft_cox( -1 , nup , cx ) ; 00318 for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; } 00319 csfft_cox( 1 , nup , cx ) ; 00320 for( ii=0 ; ii<ny ; ii++){ qfim[ii*nx] = cx[ii] ; } 00321 } 00322 } 00323 } 00324 break ; 00325 } 00326 00327 /*** do z-direction ***/ 00328 00329 STATUS("start z FFTs") ; 00330 00331 aa = sigmaz * sigmaz * 0.5 ; 00332 nup = nz + (int)(3.0 * sigmaz / dz) ; /* min FFT length */ 00333 #if 0 00334 ii = 2 ; while( ii < nup ){ ii *= 2 ; } /* next power of 2 larger */ 00335 #else 00336 ii = csfft_nextup_one35(nup) ; 00337 #endif 00338 nup = ii ; nby2 = nup / 2 ; 00339 00340 GET_AS_BIG(cx,complex,nup) ; GET_AS_BIG(gg,float,nup) ; 00341 00342 dk = (2.0*PI) / (nup * dz) ; 00343 fac = 1.0 / nup ; 00344 gg[0] = fac ; 00345 for( ii=1 ; ii<=nby2 ; ii++ ){ k=ii*dk; gg[nup-ii]=gg[ii]=fac*exp(-aa*k*k); } 00346 00347 switch( ftype ){ 00348 case MRI_short:{ 00349 register short * qfim ; 00350 for( kk=0 ; kk < ny ; kk++ ){ 00351 for( jj=0 ; jj < nx ; jj+=2 ){ 00352 base = jj + kk*nx ; 00353 qfim = sfim + base ; 00354 if( jj == nx-1 ) 00355 for( ii=0 ; ii<nz ; ii++){ cx[ii].r = qfim[ii*nxy] ; cx[ii].i = 0.0 ; } 00356 else 00357 for( ii=0 ; ii<nz ; ii++){ cx[ii].r = qfim[ii*nxy] ; cx[ii].i = qfim[ii*nxy+1] ; } 00358 for( ii=nz; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; } 00359 csfft_cox( -1 , nup , cx ) ; 00360 for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; } 00361 csfft_cox( 1 , nup , cx ) ; 00362 if( jj == nx-1 ) 00363 for( ii=0 ; ii<nz ; ii++){ qfim[ii*nxy] = cx[ii].r ; } 00364 else 00365 for( ii=0 ; ii<nz ; ii++){ qfim[ii*nxy] = cx[ii].r ; qfim[ii*nxy+1] = cx[ii].i ; } 00366 } 00367 } 00368 } 00369 break ; 00370 00371 case MRI_float:{ 00372 register float * qfim ; 00373 for( kk=0 ; kk < ny ; kk++ ){ 00374 for( jj=0 ; jj < nx ; jj+=2 ){ 00375 base = jj + kk*nx ; 00376 qfim = ffim + base ; 00377 if( jj == nx-1 ) 00378 for( ii=0 ; ii<nz ; ii++){ cx[ii].r = qfim[ii*nxy] ; cx[ii].i = 0.0 ; } 00379 else 00380 for( ii=0 ; ii<nz ; ii++){ cx[ii].r = qfim[ii*nxy] ; cx[ii].i = qfim[ii*nxy+1] ; } 00381 for( ii=nz; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; } 00382 csfft_cox( -1 , nup , cx ) ; 00383 for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; } 00384 csfft_cox( 1 , nup , cx ) ; 00385 if( jj == nx-1 ) 00386 for( ii=0 ; ii<nz ; ii++){ qfim[ii*nxy] = cx[ii].r ; } 00387 else 00388 for( ii=0 ; ii<nz ; ii++){ qfim[ii*nxy] = cx[ii].r ; qfim[ii*nxy+1] = cx[ii].i ; } 00389 } 00390 } 00391 } 00392 break ; 00393 00394 case MRI_byte:{ 00395 register byte * qfim ; 00396 for( kk=0 ; kk < ny ; kk++ ){ 00397 for( jj=0 ; jj < nx ; jj+=2 ){ 00398 base = jj + kk*nx ; 00399 qfim = bfim + base ; 00400 if( jj == nx-1 ) 00401 for( ii=0 ; ii<nz ; ii++){ cx[ii].r = qfim[ii*nxy] ; cx[ii].i = 0.0 ; } 00402 else 00403 for( ii=0 ; ii<nz ; ii++){ cx[ii].r = qfim[ii*nxy] ; cx[ii].i = qfim[ii*nxy+1] ; } 00404 for( ii=nz; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; } 00405 csfft_cox( -1 , nup , cx ) ; 00406 for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; } 00407 csfft_cox( 1 , nup , cx ) ; 00408 if( jj == nx-1 ) 00409 for( ii=0 ; ii<nz ; ii++){ qfim[ii*nxy] = cx[ii].r ; } 00410 else 00411 for( ii=0 ; ii<nz ; ii++){ qfim[ii*nxy] = cx[ii].r ; qfim[ii*nxy+1] = cx[ii].i ; } 00412 } 00413 } 00414 } 00415 break ; 00416 00417 case MRI_complex:{ 00418 register complex * qfim ; 00419 for( kk=0 ; kk < ny ; kk++ ){ 00420 for( jj=0 ; jj < nx ; jj++ ){ 00421 base = jj + kk*nx ; 00422 qfim = cfim + base ; 00423 for( ii=0 ; ii<nz ; ii++){ cx[ii] = qfim[ii*nxy] ; } 00424 for( ii=nz; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; } 00425 csfft_cox( -1 , nup , cx ) ; 00426 for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; } 00427 csfft_cox( 1 , nup , cx ) ; 00428 for( ii=0 ; ii<nz ; ii++){ qfim[ii*nxy] = cx[ii] ; } 00429 } 00430 } 00431 } 00432 break ; 00433 } 00434 00435 /*** 10 Jan 2003: clip data to bot and top found above ***/ 00436 /*** to minimize Gibbs ringing artifacts ***/ 00437 00438 switch( ftype ){ 00439 00440 case MRI_short: 00441 for( ii=0 ; ii < nxyz ; ii++ ) 00442 if( sfim[ii] < fbot ) sfim[ii] = fbot ; 00443 else if( sfim[ii] > ftop ) sfim[ii] = ftop ; 00444 break ; 00445 00446 case MRI_float: 00447 for( ii=0 ; ii < nxyz ; ii++ ) 00448 if( ffim[ii] < fbot ) ffim[ii] = fbot ; 00449 else if( ffim[ii] > ftop ) ffim[ii] = ftop ; 00450 break ; 00451 00452 case MRI_byte: 00453 for( ii=0 ; ii < nxyz ; ii++ ) 00454 if( bfim[ii] < fbot ) bfim[ii] = fbot ; 00455 else if( bfim[ii] > ftop ) bfim[ii] = ftop ; 00456 break ; 00457 } 00458 00459 /*** done! ***/ 00460 00461 EXRETURN ; 00462 } |