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 }
|
