/***************************************************************************** Major portions of this software are copyrighted by the Medical College of Wisconsin, 1994-2000, and are released under the Gnu General Public License, Version 2. See the file README.Copyright for details. ******************************************************************************/ #include "mrilib.h" /**************************************************************************/ /***** prototypes for FIR blurring functions at the end of this file ******/ void fir_blurx( int m, float *wt,int nx, int ny, int nz, float *f ) ; void fir_blury( int m, float *wt,int nx, int ny, int nz, float *f ) ; void fir_blurz( int m, float *wt,int nx, int ny, int nz, float *f ) ; static void fir_gaussian_load( int m, float dx, float *ff ) ; /* 03 Apr 2007 */ #undef FIR_MAX #define FIR_MAX 35 /** max length of FIR filter to use instead of FFTs **/ #undef DO_INFO #define DO_INFO 0 /** for debugging and timing **/ /*! If set to 0, EDIT_blur_volume() will not use the fir_blur? functions. */ static int allow_fir = 1 ; /*! Controls whether EDIT_blur_volume() will use the fir_blur? functions. */ void EDIT_blur_allow_fir( int i ){ allow_fir = i; } /* mirror blur array at ends, all the way */ /* (as macro for speed) */ /* pad for full length m [9 Feb 2022 rickr] */ /* RRR - probably 'r' array */ /* MMM - padding length (on each side) */ /* NDIM - unnpadded length, NX, NY or NZ */ #ifdef BMIRROR #undef BMIRROR #endif #define BMIRROR(RRR,MMM,NDIM) \ do { \ for( mm=1; mm <= MMM; mm++) { \ RRR[MMM-mm] = RRR[MMM+mm]; \ RRR[MMM+NDIM-1+mm] = RRR[MMM+NDIM-mm]; \ } \ } while(0) /*************************************************************************/ /*! Routine to blur a 3D volume with a Gaussian, using FFTs. If the blurring parameter (sigma) is small enough, will use real-space FIR filter instead. This function actually just calls EDIT_blur_volume_3d() to do the real work. **************************************************************************/ void EDIT_blur_volume( int nx, int ny, int nz, float dx, float dy, float dz, int ftype , void * vfim , float sigma ) { EDIT_blur_volume_3d (nx, ny, nz, dx, dy, dz, ftype, vfim, sigma, sigma, sigma); } /**************************************************************************/ /*! The following slightly modified version of EDIT_blur_volume allows independent specification of Gaussian filter widths along the three perpendicular axes. - BDW - 21 Feb 1997 - RWC - 04 Feb 2005: use fir_blur? function if sigma? is small - also see EDIT_blur_allow_fir() and FIR_blur_volume_3d() --------------------------------------------------------------------------*/ 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 ) { int jj,kk , nxy , base,nby2 ; float dk , aa , k , fac ; register int ii , nup ; complex *cx = NULL ; float *gg = NULL ; int ncg=0 ; byte *bfim = NULL ; /* pointers to data array vfim */ short *sfim = NULL ; float *ffim = NULL ; complex *cfim = NULL ; float fbot=0.0,ftop=0.0 ; /* 10 Jan 2003: for clipping results */ int nxyz ; /* number of voxels */ int fir_m , fir_num=0 ; /* 03 Oct 2005: for fir_blur? filtering */ float fir_wt[FIR_MAX+1] ; int all_fir=0 ; /* 06 Oct 2005 */ int as_float = 0 ; /* process as float (to apply FIR) 15 Jun 2009 [rickr] */ /* - convert short/byte to float, blur, convert back */ int dtype = ftype ; /* processing datatype (ftype, if no conversion) */ double sfac = AFNI_numenv("AFNI_BLUR_FIRFAC") ; /* <2 -> <=0 to allow for truncated blurs [8 May 2019 rickr] */ /* */ /* This can be abused for a "fast ANATICOR", for example. */ /* Since sigma = 0.4246609 * fwhm, consider using: */ /* sfac = 1/(2*.0.4246609) = 1.17741 */ /* That number of sigmas should match the half width at half max, */ /* which should terminate the blur just after a half height. */ /* */ /* Or use 2*FWHM and sfac = 1.17741/2 = 0.588705 to make it more flat,*/ /* with a min contribution of ~0.84, rather than 0.5, yet limiting */ /* the output to the same HWHM radius (e.g. FWHM=80mm with sfac=0.589 */ /* results in a fairly flat blur out to a radius of ~20 mm). */ /* */ /* By default, sfac = 2.5 leads to FWHM ~= sfac*sigma. Since FIR */ /* length is a radius, FIR length ends twice as far out as HWHM. */ /* So if e^(-x^2) = 0.5, e^(-(2x)^2) = (e^(-x^2))^4 = 0.5^4 = 0.0625, */ /* which is the fractional limit of the FIR. */ if( sfac <= 0.0 ) sfac = 2.5 ; /***---------- initialize ----------***/ ENTRY("EDIT_blur_volume_3d") ; if( vfim == NULL ) EXRETURN ; /* no data? */ if( sigmax <= 0.0 && sigmay <= 0.0 && sigmaz <= 0.0 ) EXRETURN ; if( dx <= 0.0 ) dx = 1.0 ; /* 03 Oct 2005: regularize grid sizes */ if( dy <= 0.0 ) dy = dx ; if( dz <= 0.0 ) dz = dx ; switch( ftype ){ /* cast pointer to correct type */ default: EXRETURN ; case MRI_short: sfim = (short *) vfim ; break ; case MRI_float: ffim = (float *) vfim ; break ; case MRI_byte: bfim = (byte *) vfim ; break ; case MRI_complex: cfim = (complex *) vfim ; break ; } nxy = nx * ny ; nxyz = nxy * nz ; /*** 10 Jan 2003: find bot and top of data input */ if( allow_fir ){ ii = (int) ceil( sfac * sigmax / dx ) ; jj = (int) ceil( sfac * sigmay / dy ) ; kk = (int) ceil( sfac * sigmaz / dz ) ; if( ii <= FIR_MAX && jj <= FIR_MAX && kk <= FIR_MAX ) all_fir = 1 ; { /* print some details, once per execution [14 Feb 2025 rickr] */ static int bcomment = 0; if( ! bcomment ) { INFO_message("edt_blur: sigx %g, dx %g, fwhm %g, sfac %g, firx %d\n", sigmax, dx, sigmax/0.4246609, sfac, ii); bcomment = 1; } } /* try to process most data using Finite Impule Response, not with Fourier interpolation (leads to ringing problems both inside and outside a mask) - set AFNI_BLUR_INTS_AS_OLD=YES to use old method 15 Jun 2009 [rickr] */ if( (ftype == MRI_short || ftype == MRI_byte) && ! AFNI_yesenv("AFNI_BLUR_INTS_AS_OLD") ) { ffim = (float *)malloc(nxyz * sizeof(float)); as_float = 1; /* flag, redundant but clear */ dtype = MRI_float; /* process data as float */ } } if( ftype != MRI_float ) all_fir = 0 ; /* 17 Nov 2005: oopsie */ if( !all_fir ){ switch( ftype ){ default: fbot = ftop = 0.0 ; /* for complex */ break ; case MRI_short: fbot = ftop = sfim[0] ; for( ii=1 ; ii < nxyz ; ii++ ) { if( sfim[ii] < fbot ) fbot = sfim[ii] ; else if( sfim[ii] > ftop ) ftop = sfim[ii] ; if( as_float ) ffim[ii] = (float)sfim[ii] ; /* copy data as int */ } break ; case MRI_float: fbot = ftop = ffim[0] ; for( ii=1 ; ii < nxyz ; ii++ ) if( ffim[ii] < fbot ) fbot = ffim[ii] ; else if( ffim[ii] > ftop ) ftop = ffim[ii] ; break ; case MRI_byte: fbot = ftop = bfim[0] ; for( ii=1 ; ii < nxyz ; ii++ ) { if( bfim[ii] < fbot ) fbot = bfim[ii] ; else if( bfim[ii] > ftop ) ftop = bfim[ii] ; if( as_float ) ffim[ii] = (float)bfim[ii] ; /* copy data as int */ } break ; } } /*** do x-direction ***/ /** 03 Oct 2005: perhaps do the x-blur in real-space? **/ if( nx < 2 || sigmax <= 0.0 ){ STATUS("skipping x blur") ; fir_num++ ; goto DO_Y_BLUR ; } fir_m = (int) ceil( sfac * sigmax / dx ) ; if( allow_fir && dtype == MRI_float && fir_m <= FIR_MAX ){ STATUS("start x FIR") ; if( fir_m < 1 ) fir_m = 1 ; fir_gaussian_load( fir_m , dx/sigmax , fir_wt ) ; #if 0 fac = fir_wt[0] = 1.0f ; for( ii=1 ; ii <= fir_m ; ii++ ){ fir_wt[ii] = exp(-0.5*(ii*dx)*(ii*dx)/(sigmax*sigmax)) ; fac += 2.0f * fir_wt[ii] ; } fac = 1.0f / fac ; for( ii=0 ; ii <= fir_m ; ii++ ) fir_wt[ii] *= fac ; #endif fir_blurx( fir_m , fir_wt , nx,ny,nz , ffim ) ; fir_num++ ; goto DO_Y_BLUR ; } STATUS("start x FFTs") ; aa = sigmax * sigmax * 0.5 ; nup = nx + (int)(3.0 * sigmax / dx) ; /* min FFT length */ nup = csfft_nextup_one35(nup) ; nby2 = nup / 2 ; if(DO_INFO) INFO_message("x FFT(%d)",nup) ; if( nup > ncg ){ cx = (complex *)realloc(cx,sizeof(complex)*nup) ; gg = (float *)realloc(gg,sizeof(float )*nup) ; ncg = nup ; } dk = (2.0*PI) / (nup * dx) ; fac = 1.0 / nup ; gg[0] = fac ; for( ii=1 ; ii<=nby2 ; ii++ ){ k=ii*dk; gg[nup-ii]=gg[ii]=fac*exp(-aa*k*k); } /** July 20: double up on FFTs **/ /** Feb 09: extend to other data types besides shorts; doubling up does not apply to complex data! **/ switch( dtype ){ case MRI_short:{ register short *qfim ; for( kk=0 ; kk < nz ; kk++ ){ for( jj=0 ; jj < ny ; jj+=2 ){ base = jj*nx + kk*nxy ; qfim = sfim + base ; if( jj == ny-1 ) for( ii=0 ; ii ncg ){ cx = (complex *)realloc(cx,sizeof(complex)*nup) ; gg = (float *)realloc(gg,sizeof(float )*nup) ; ncg = nup ; } dk = (2.0*PI) / (nup * dy) ; fac = 1.0 / nup ; gg[0] = fac ; for( ii=1 ; ii<=nby2 ; ii++ ){ k=ii*dk; gg[nup-ii]=gg[ii]=fac*exp(-aa*k*k); } switch( dtype ){ case MRI_short:{ register short *qfim ; for( kk=0 ; kk < nz ; kk++ ){ for( jj=0 ; jj < nx ; jj+=2 ){ base = jj + kk*nxy ; qfim = sfim + base ; if( jj == nx-1 ) for( ii=0 ; ii ncg ){ cx = (complex *)realloc(cx,sizeof(complex)*nup) ; gg = (float *)realloc(gg,sizeof(float )*nup) ; ncg = nup ; } dk = (2.0*PI) / (nup * dz) ; fac = 1.0 / nup ; gg[0] = fac ; for( ii=1 ; ii<=nby2 ; ii++ ){ k=ii*dk; gg[nup-ii]=gg[ii]=fac*exp(-aa*k*k); } switch( dtype ){ case MRI_short:{ register short *qfim ; for( kk=0 ; kk < ny ; kk++ ){ for( jj=0 ; jj < nx ; jj+=2 ){ base = jj + kk*nx ; qfim = sfim + base ; if( jj == nx-1 ) for( ii=0 ; ii= 0.0f ) sfim[ii] = (short)(ffim[ii]+0.5) ; else sfim[ii] = (short)(ffim[ii]-0.5) ; } if( sfim[ii] < fbot ) sfim[ii] = fbot ; else if( sfim[ii] > ftop ) sfim[ii] = ftop ; } break ; case MRI_float: for( ii=0 ; ii < nxyz ; ii++ ) if( ffim[ii] < fbot ) ffim[ii] = fbot ; else if( ffim[ii] > ftop ) ffim[ii] = ftop ; break ; case MRI_byte: for( ii=0 ; ii < nxyz ; ii++ ) { if( as_float ) bfim[ii] = (byte)(ffim[ii]+0.5) ; /* return to byte */ if( bfim[ii] < fbot ) bfim[ii] = fbot ; else if( bfim[ii] > ftop ) bfim[ii] = ftop ; } break ; } } /*** done! ***/ if( cx != NULL ) free(cx) ; if( gg != NULL ) free(gg) ; if( as_float && ffim ) free(ffim) ; EXRETURN ; } #if 0 /*-------------------------------------------------------------------*/ /*! Partition of unity function: n=+infinity sum partu(x-n) = 1 identically n=-infinity and support of partu(x) is interval [-1,1] *//*-----------------------------------------------------------------*/ static INLINE float partu( float x ) { register float ax ; ax = fabsf(x) ; if( ax >= 1.0f ) return 0.0f ; return (x*x*(20.0f*ax+10.0f)+4.0f*ax+1.0f) ; } #endif /*-------------------------------------------------------------------*/ /*! Function to blur a 3D volume in-place with a symmetric FIR filter along the x-direction. - m = stencil size (+/- m points in each direction; m >= 1) - wt = array of weights - nx,ny,nz = dimensions of 3D array - f = 3D array f_out[i,j,k] = wt[0] * f_in[i,j,k] + wt[1] *(f_in[i+1,j,k]+f_in[i-1,j,k)) + wt[2] *(f_in[i+2,j,k]+f_in[i-2,j,k)) + ... + wt[m] *(f_in[i+m,j,k]+f_in[i-m,j,k)) Similar routines for blurring along the y- and z-directions are fir_blury() and fir_blurz(). -- RWCox - 03 Oct 2005 - trying for some speedup for Daniel Glen ---------------------------------------------------------------------*/ void fir_blurx( int m, float *wt,int nx, int ny, int nz, float *f ) { int ii,jj,kk,qq , nxy=nx*ny , off ; int mm ; /* for padding ends */ float *r , wt0=0.0,wt1=0.0,wt2=0.0,wt3=0.0,wt4=0.0,wt5=0.0,wt6=0.0,wt7=0.0 , sum , *ff ; ENTRY("fir_blurx") ; if(PRINT_TRACING){char str[256];sprintf(str,"m=%d",m);STATUS(str);} if( m < 1 || wt == NULL || nx < (m+1) || f == NULL ) EXRETURN ; if( ny <= 0 || nz <= 0 ) EXRETURN ; switch(m){ /** assign weights to variables not arrays **/ case 7: wt7 = wt[7]; /* let cases fall through to next case to assign weights */ case 6: wt6 = wt[6]; case 5: wt5 = wt[5]; case 4: wt4 = wt[4]; case 3: wt3 = wt[3]; case 2: wt2 = wt[2]; case 1: wt1 = wt[1]; case 0: wt0 = wt[0]; default: break ; } /* 1 row workspace, with m-long buffers at each end (so that the i-th element of the row is in r[i+m]) */ r = (float *)calloc(sizeof(float),(nx+2*m)) ; switch( m ){ default: for( kk=0 ; kk < nz ; kk++ ){ for( jj=0 ; jj < ny ; jj++ ){ off = jj*nx + kk*nxy ; ff = f+off ; /* ff = ptr to this row */ memcpy( r+m , ff , sizeof(float)*nx ) ; /* copy row into workspace */ /* replace: r[m-1] = r[m+1] ; r[nx+m] = r[nx+m-2] ; */ BMIRROR(r,m,nx); /* fully mirror at ends [9 Feb 2022 rickr]*/ for( ii=0 ; ii < nx ; ii++ ){ /* filter at ii-th location */ sum = wt[0]*r[ii+m] ; for( qq=1 ; qq <= m ; qq++ ) sum += wt[qq] * ( r[ii+m-qq] + r[ii+m+qq] ) ; ff[ii] = sum ; /* save result back in input array */ } }} break ; /** for the cases below, the innermost loop (qq, above) is completely unrolled for speedup **/ #undef M #define M 7 case 7: for( kk=0 ; kk < nz ; kk++ ){ for( jj=0 ; jj < ny ; jj++ ){ off = jj*nx + kk*nxy ; ff = f+off ; memcpy( r+m , ff , sizeof(float)*nx ) ; BMIRROR(r,M,nx); /* mirror at ends, all the way */ for( ii=0 ; ii < nx ; ii++ ) ff[ii] = wt7*(r[ii ]+r[ii+14]) +wt6*(r[ii+1]+r[ii+13]) +wt5*(r[ii+2]+r[ii+12]) +wt4*(r[ii+3]+r[ii+11]) +wt3*(r[ii+4]+r[ii+10]) +wt2*(r[ii+5]+r[ii+ 9]) +wt1*(r[ii+6]+r[ii+ 8])+wt0*r[ii+7] ; }} break ; #undef M #define M 6 case 6: for( kk=0 ; kk < nz ; kk++ ){ for( jj=0 ; jj < ny ; jj++ ){ off = jj*nx + kk*nxy ; ff = f+off ; memcpy( r+m , ff , sizeof(float)*nx ) ; BMIRROR(r,M,nx); /* mirror at ends, all the way */ for( ii=0 ; ii < nx ; ii++ ) ff[ii] = wt6*(r[ii ]+r[ii+12]) +wt5*(r[ii+1]+r[ii+11]) +wt4*(r[ii+2]+r[ii+10]) +wt3*(r[ii+3]+r[ii+ 9]) +wt2*(r[ii+4]+r[ii+ 8]) +wt1*(r[ii+5]+r[ii+ 7])+wt0*r[ii+6] ; }} break ; #undef M #define M 5 case 5: for( kk=0 ; kk < nz ; kk++ ){ for( jj=0 ; jj < ny ; jj++ ){ off = jj*nx + kk*nxy ; ff = f+off ; memcpy( r+m , ff , sizeof(float)*nx ) ; BMIRROR(r,M,nx); /* mirror at ends, all the way */ for( ii=0 ; ii < nx ; ii++ ) ff[ii] = wt5*(r[ii ]+r[ii+10]) +wt4*(r[ii+1]+r[ii+ 9]) +wt3*(r[ii+2]+r[ii+ 8]) +wt2*(r[ii+3]+r[ii+ 7]) +wt1*(r[ii+4]+r[ii+ 6])+wt0*r[ii+5] ; }} break ; #undef M #define M 4 case 4: for( kk=0 ; kk < nz ; kk++ ){ for( jj=0 ; jj < ny ; jj++ ){ off = jj*nx + kk*nxy ; ff = f+off ; memcpy( r+m , ff , sizeof(float)*nx ) ; BMIRROR(r,M,nx); /* mirror at ends, all the way */ for( ii=0 ; ii < nx ; ii++ ) ff[ii] = wt4*(r[ii ]+r[ii+ 8]) +wt3*(r[ii+1]+r[ii+ 7]) +wt2*(r[ii+2]+r[ii+ 6]) +wt1*(r[ii+3]+r[ii+ 5])+wt0*r[ii+4] ; }} break ; #undef M #define M 3 case 3: for( kk=0 ; kk < nz ; kk++ ){ for( jj=0 ; jj < ny ; jj++ ){ off = jj*nx + kk*nxy ; ff = f+off ; memcpy( r+m , ff , sizeof(float)*nx ) ; BMIRROR(r,M,nx); /* mirror at ends, all the way */ for( ii=0 ; ii < nx ; ii++ ) ff[ii] = wt3*(r[ii ]+r[ii+ 6]) +wt2*(r[ii+1]+r[ii+ 5]) +wt1*(r[ii+2]+r[ii+ 4])+wt0*r[ii+3] ; }} break ; #undef M #define M 2 case 2: for( kk=0 ; kk < nz ; kk++ ){ for( jj=0 ; jj < ny ; jj++ ){ off = jj*nx + kk*nxy ; ff = f+off ; memcpy( r+m , ff , sizeof(float)*nx ) ; BMIRROR(r,M,nx); /* mirror at ends, all the way */ for( ii=0 ; ii < nx ; ii++ ) ff[ii] = wt2*(r[ii ]+r[ii+ 4]) +wt1*(r[ii+1]+r[ii+ 3])+wt0*r[ii+2] ; }} break ; #undef M #define M 1 case 1: for( kk=0 ; kk < nz ; kk++ ){ for( jj=0 ; jj < ny ; jj++ ){ off = jj*nx + kk*nxy ; ff = f+off ; memcpy( r+m , ff , sizeof(float)*nx ) ; BMIRROR(r,M,nx); /* mirror at ends, all the way */ for( ii=0 ; ii < nx ; ii++ ) ff[ii] = wt1*(r[ii]+r[ii+2])+wt0*r[ii+1] ; }} break ; } /* end of switch on m */ free((void *)r) ; EXRETURN ; } /*-------------------------------------------------------------------*/ #undef D #define D nx /* stride along y-axis */ /*! Similar to fir_blurx(), but along the y-axis. For further comments, see fir_blurx() source code. */ void fir_blury( int m, float *wt,int nx, int ny, int nz, float *f ) { int ii,jj,kk,qq , nxy=nx*ny , off ; int mm ; /* for padding ends */ float *r, wt0=0.0,wt1=0.0,wt2=0.0,wt3=0.0,wt4=0.0,wt5=0.0,wt6=0.0,wt7=0.0 , sum , *ff ; float *rr, *ss; int ny2m = ny+2*m; ENTRY("fir_blury") ; if(PRINT_TRACING){char str[256];sprintf(str,"m=%d",m);STATUS(str);} if( m < 1 || wt == NULL || ny < (m+1) || f == NULL ) EXRETURN ; if( nx <= 0 || nz <= 0 ) EXRETURN ; switch(m){ /**assign weights to variables not arrays **/ case 7: wt7 = wt[7]; /* let cases fall through to next case to assign weights */ case 6: wt6 = wt[6]; case 5: wt5 = wt[5]; case 4: wt4 = wt[4]; case 3: wt3 = wt[3]; case 2: wt2 = wt[2]; case 1: wt1 = wt[1]; case 0: wt0 = wt[0]; default: break ; } if( nx < 512 ) goto SMALLIMAGE; /* In this function, for each value of kk (z index), we extract a 2D (y,x) slice, with m-long buffers on each side in the y-direction. The purpose of this is to get multiple lines of y-direction data into the CPU cache, to speed up processing (a lot). For the x-axis, this was unneeded, since the x-rows are contiguous in memory. For data at 256x256 this 2D extract/process/insert trick was nugatory. However, for 512x512 data this trick becomes important. The same method is used for the z-axis in fir_blurz(). */ /* macro to access the input data 2D slice: (i,j) = (x,y) indexes */ #undef RR #define RR(i,j) rr[(j)+m+(i)*ny2m] /*** 0 <= i <= nx-1 ; -m <= m <= ny-1+m ***/ /* macro to access the output data 2D slice */ #undef SS #define SS(i,k) ss[(k)+(i)*ny] rr = (float *)calloc(sizeof(float),ny2m*nx) ; /* ny2m = ny+2*m */ ss = (float *)malloc(sizeof(float)*ny *nx) ; for( kk=0 ; kk < nz ; kk++ ){ /* loop in z-direction (an xy slice at a time) */ off = kk*nxy ; ff = f+off ; /* ff = ptr to start of this 2D slice */ /* load data into 2D (y+2m,x) slice from 3D (x,y,z) array; inner loop is over ii so as to access in the most contiguous way */ for( jj=0 ; jj < ny ; jj++ ){ for( ii=0 ; ii < nx ; ii++ ) RR(ii,jj) = ff[ii+D*jj] ; /* D = nx here */ } for( ii=0 ; ii < nx ; ii++ ){ RR(ii,-1) = RR(ii,1) ; RR(ii,ny) = RR(ii,ny-2) ; /* edge reflection - */ /* only 1 point reflected*/ } /* filter data in RR along y-direction, put into 2D SS array */ switch(m){ /** for small m, unroll the inner loop for speed **/ default: for( ii=0 ; ii < nx ; ii++ ){ for( jj=0 ; jj < ny ; jj++ ){ sum = wt[0]*RR(ii,jj) ; for( qq=1 ; qq <= m ; qq++ ) sum += wt[qq] * ( RR(ii,jj+qq) + RR(ii,jj-qq) ) ; SS(ii,jj) = sum ; }} break ; case 7: for( ii=0 ; ii < nx ; ii++ ){ for( jj=0 ; jj < ny ; jj++ ){ SS(ii,jj) = wt7 * ( RR(ii,jj+7) + RR(ii,jj-7) ) + wt6 * ( RR(ii,jj+6) + RR(ii,jj-6) ) + wt5 * ( RR(ii,jj+5) + RR(ii,jj-5) ) + wt4 * ( RR(ii,jj+4) + RR(ii,jj-4) ) + wt3 * ( RR(ii,jj+3) + RR(ii,jj-3) ) + wt2 * ( RR(ii,jj+2) + RR(ii,jj-2) ) + wt1 * ( RR(ii,jj+1) + RR(ii,jj-1) ) + wt0 * RR(ii,jj) ; }} break ; case 6: for( ii=0 ; ii < nx ; ii++ ){ for( jj=0 ; jj < ny ; jj++ ){ SS(ii,jj) = wt6 * ( RR(ii,jj+6) + RR(ii,jj-6) ) + wt5 * ( RR(ii,jj+5) + RR(ii,jj-5) ) + wt4 * ( RR(ii,jj+4) + RR(ii,jj-4) ) + wt3 * ( RR(ii,jj+3) + RR(ii,jj-3) ) + wt2 * ( RR(ii,jj+2) + RR(ii,jj-2) ) + wt1 * ( RR(ii,jj+1) + RR(ii,jj-1) ) + wt0 * RR(ii,jj) ; }} break ; case 5: for( ii=0 ; ii < nx ; ii++ ){ for( jj=0 ; jj < ny ; jj++ ){ SS(ii,jj) = wt5 * ( RR(ii,jj+5) + RR(ii,jj-5) ) + wt4 * ( RR(ii,jj+4) + RR(ii,jj-4) ) + wt3 * ( RR(ii,jj+3) + RR(ii,jj-3) ) + wt2 * ( RR(ii,jj+2) + RR(ii,jj-2) ) + wt1 * ( RR(ii,jj+1) + RR(ii,jj-1) ) + wt0 * RR(ii,jj) ; }} break ; case 4: for( ii=0 ; ii < nx ; ii++ ){ for( jj=0 ; jj < ny ; jj++ ){ SS(ii,jj) = wt4 * ( RR(ii,jj+4) + RR(ii,jj-4) ) + wt3 * ( RR(ii,jj+3) + RR(ii,jj-3) ) + wt2 * ( RR(ii,jj+2) + RR(ii,jj-2) ) + wt1 * ( RR(ii,jj+1) + RR(ii,jj-1) ) + wt0 * RR(ii,jj) ; }} break ; case 3: for( ii=0 ; ii < nx ; ii++ ){ for( jj=0 ; jj < ny ; jj++ ){ SS(ii,jj) = wt3 * ( RR(ii,jj+3) + RR(ii,jj-3) ) + wt2 * ( RR(ii,jj+2) + RR(ii,jj-2) ) + wt1 * ( RR(ii,jj+1) + RR(ii,jj-1) ) + wt0 * RR(ii,jj) ; }} break ; case 2: for( ii=0 ; ii < nx ; ii++ ){ for( jj=0 ; jj < ny ; jj++ ){ SS(ii,jj) = wt2 * ( RR(ii,jj+2) + RR(ii,jj-2) ) + wt1 * ( RR(ii,jj+1) + RR(ii,jj-1) ) + wt0 * RR(ii,jj) ; }} break ; case 1: for( ii=0 ; ii < nx ; ii++ ){ for( jj=0 ; jj < ny ; jj++ ){ SS(ii,jj) = wt1 * ( RR(ii,jj+1) + RR(ii,jj-1) ) + wt0 * RR(ii,jj) ; }} break ; } /* end of special cases of m */ /* put SS array back into input 3D array; again, inner loop over ii for most contiguous access to f[] array */ for( jj=0 ; jj < ny ; jj++ ){ for( ii=0 ; ii < nx ; ii++ ) ff[ii+D*jj] = SS(ii,jj) ; } } /* end of loop over y-direction (zz) */ /*** finito, cara mia mine, oh, oh, oh, each time we part, my heart wants to die...***/ free((void *)ss) ; free((void *)rr) ; EXRETURN ; /* for small images (nx<512), use slice as is and don't use reslicing trick*/ SMALLIMAGE: r = (float *)calloc(sizeof(float),(ny+2*m)) ; switch( m ){ default: for( kk=0 ; kk < nz ; kk++ ){ for( ii=0 ; ii < nx ; ii++ ){ off = ii + kk*nxy ; ff = f+off ; for( jj=0 ; jj < ny ; jj++ ) r[jj+m] = ff[D*jj] ; BMIRROR(r,m,ny); /* mirror at ends, all the way */ for( jj=0 ; jj < ny ; jj++ ){ sum = wt[0]*r[jj+m] ; for( qq=1 ; qq <= m ; qq++ ) sum += wt[qq] * ( r[jj+m-qq] + r[jj+m+qq] ) ; ff[D*jj] = sum ; } }} break ; #undef M #define M 7 case 7: for( kk=0 ; kk < nz ; kk++ ){ for( ii=0 ; ii < nx ; ii++ ){ off = ii + kk*nxy ; ff = f+off ; for( jj=0 ; jj < ny ; jj++ ) r[jj+M] = ff[D*jj] ; BMIRROR(r,M,ny); /* mirror at ends, all the way */ for( jj=0 ; jj < ny ; jj++ ) ff[D*jj] = wt7*(r[jj ]+r[jj+14]) +wt6*(r[jj+1]+r[jj+13]) +wt5*(r[jj+2]+r[jj+12]) +wt4*(r[jj+3]+r[jj+11]) +wt3*(r[jj+4]+r[jj+10]) +wt2*(r[jj+5]+r[jj+ 9]) +wt1*(r[jj+6]+r[jj+ 8])+wt0*r[jj+7] ; }} break ; #undef M #define M 6 case 6: for( kk=0 ; kk < nz ; kk++ ){ for( ii=0 ; ii < nx ; ii++ ){ off = ii + kk*nxy ; ff = f+off ; for( jj=0 ; jj < ny ; jj++ ) r[jj+M] = ff[D*jj] ; BMIRROR(r,M,ny); /* mirror at ends, all the way */ for( jj=0 ; jj < ny ; jj++ ) ff[D*jj] = wt6*(r[jj ]+r[jj+12]) +wt5*(r[jj+1]+r[jj+11]) +wt4*(r[jj+2]+r[jj+10]) +wt3*(r[jj+3]+r[jj+ 9]) +wt2*(r[jj+4]+r[jj+ 8]) +wt1*(r[jj+5]+r[jj+ 7])+wt0*r[jj+6] ; }} break ; #undef M #define M 5 case 5: for( kk=0 ; kk < nz ; kk++ ){ for( ii=0 ; ii < nx ; ii++ ){ off = ii + kk*nxy ; ff = f+off ; for( jj=0 ; jj < ny ; jj++ ) r[jj+M] = ff[D*jj] ; BMIRROR(r,M,ny); /* mirror at ends, all the way */ for( jj=0 ; jj < ny ; jj++ ) ff[D*jj] = wt5*(r[jj ]+r[jj+10]) +wt4*(r[jj+1]+r[jj+ 9]) +wt3*(r[jj+2]+r[jj+ 8]) +wt2*(r[jj+3]+r[jj+ 7]) +wt1*(r[jj+4]+r[jj+ 6])+wt0*r[jj+5] ; }} break ; #undef M #define M 4 case 4: for( kk=0 ; kk < nz ; kk++ ){ for( ii=0 ; ii < nx ; ii++ ){ off = ii + kk*nxy ; ff = f+off ; for( jj=0 ; jj < ny ; jj++ ) r[jj+M] = ff[D*jj] ; BMIRROR(r,M,ny); /* mirror at ends, all the way */ for( jj=0 ; jj < ny ; jj++ ) ff[D*jj] = wt4*(r[jj ]+r[jj+ 8]) +wt3*(r[jj+1]+r[jj+ 7]) +wt2*(r[jj+2]+r[jj+ 6]) +wt1*(r[jj+3]+r[jj+ 5])+wt0*r[jj+4] ; }} break ; #undef M #define M 3 case 3: for( kk=0 ; kk < nz ; kk++ ){ for( ii=0 ; ii < nx ; ii++ ){ off = ii + kk*nxy ; ff = f+off ; for( jj=0 ; jj < ny ; jj++ ) r[jj+M] = ff[D*jj] ; BMIRROR(r,M,ny); /* mirror at ends, all the way */ for( jj=0 ; jj < ny ; jj++ ) ff[D*jj] = wt3*(r[jj ]+r[jj+ 6]) +wt2*(r[jj+1]+r[jj+ 5]) +wt1*(r[jj+2]+r[jj+ 4])+wt0*r[jj+3] ; }} break ; #undef M #define M 2 case 2: for( kk=0 ; kk < nz ; kk++ ){ for( ii=0 ; ii < nx ; ii++ ){ off = ii + kk*nxy ; ff = f+off ; for( jj=0 ; jj < ny ; jj++ ) r[jj+M] = ff[D*jj] ; BMIRROR(r,M,ny); /* mirror at ends, all the way */ for( jj=0 ; jj < ny ; jj++ ) ff[D*jj] = wt2*(r[jj ]+r[jj+ 4]) +wt1*(r[jj+1]+r[jj+ 3])+wt0*r[jj+2] ; }} break ; #undef M #define M 1 case 1: for( kk=0 ; kk < nz ; kk++ ){ for( ii=0 ; ii < nx ; ii++ ){ off = ii + kk*nxy ; ff = f+off ; for( jj=0 ; jj < ny ; jj++ ) r[jj+M] = ff[D*jj] ; BMIRROR(r,M,ny); /* mirror at ends, all the way */ for( jj=0 ; jj < ny ; jj++ ) ff[D*jj] = wt1*(r[jj]+r[jj+2])+wt0*r[jj+1] ; }} break ; } /* end of switch on m */ free((void *)r) ; EXRETURN ; } /*-------------------------------------------------------------------*/ #undef D #define D nxy /* stride along z-axis */ /*! Similar to fir_blurx(), but along z-axis. */ void fir_blurz( int m, float *wt,int nx, int ny, int nz, float *f ) { int ii,jj,kk,qq , nxy=nx*ny , off ; int mm ; /* for padding ends */ float *rr,*ss , wt0=0.0,wt1=0.0,wt2=0.0,wt3=0.0,wt4=0.0,wt5=0.0,wt6=0.0,wt7=0.0 , sum , *ff ; int nz2m = nz+2*m ; ENTRY("fir_blurz") ; if(PRINT_TRACING){char str[256];sprintf(str,"m=%d",m);STATUS(str);} if( m < 1 || wt == NULL || nz < (m+1) || f == NULL ) EXRETURN ; if( nxy <= 0 ) EXRETURN ; /* In this function, for each value of jj (y index), we extract a 2D (z,x) slice, with m-long buffers on each side in the z-direction. The purpose of this is to get multiple lines of z-direction data into the CPU cache, to speed up processing (a lot). For the x-axis, this was unneeded, since the x-rows are contiguous in memory. For the y-axis, this trick might help, but only if a single (x,y) plane doesn't fit into cache. For nx=ny=256, 1 plane is 256 KB, so I decided that this 2D extract/process/insert trick was nugatory. */ switch(m){ /**assign weights to variables not arrays **/ case 7: wt7 = wt[7]; /* let cases fall through to next case to assign weights */ case 6: wt6 = wt[6]; case 5: wt5 = wt[5]; case 4: wt4 = wt[4]; case 3: wt3 = wt[3]; case 2: wt2 = wt[2]; case 1: wt1 = wt[1]; case 0: wt0 = wt[0]; default: break ; } /* macro to access the input data 2D slice: (i,k) = (x,z) indexes */ #undef RR #define RR(i,k) rr[(k)+m+(i)*nz2m] /*** 0 <= i <= nx-1 ; -m <= k <= nz-1+m ***/ /* macro to access the output data 2D slice */ #undef SS #define SS(i,k) ss[(k)+(i)*nz] /* mirror into blur pads, like BMIRROR, but for RR access [9 Feb 2022 rickr] */ #ifdef RMIRROR #undef RMIRROR #endif #define RMIRROR(MMM,NDIM) \ do { \ for( mm=1; mm <= MMM; mm++) { \ RR(ii,-mm) = RR(ii,mm); \ RR(ii,NDIM-1+mm) = RR(ii,NDIM-1-mm); \ } \ } while(0) rr = (float *)calloc(sizeof(float),nz2m*nx) ; /* nz2m = nz+2*m */ ss = (float *)malloc(sizeof(float)*nz *nx) ; for( jj=0 ; jj < ny ; jj++ ){ /* loop in y-direction */ off = jj*nx ; ff = f+off ; /* ff = ptr to start of this 2D slice */ /* load data into 2D (z,x) slice from 3D (x,y,z) array; inner loop is over ii so as to access in the most contiguous way */ for( kk=0 ; kk < nz ; kk++ ){ for( ii=0 ; ii < nx ; ii++ ) RR(ii,kk) = ff[ii+D*kk] ; } for( ii=0 ; ii < nx ; ii++ ){ /* replace: RR(ii,-1) = RR(ii,1) ; RR(ii,nz) = RR(ii,nz-2) ; */ RMIRROR(m,nz); /* edge reflection in z [9 Feb 2022 rickr] */ } /* filter data in RR along z-direction, put into 2D SS array */ switch(m){ /** for small m, unroll the inner loop for speed **/ default: for( ii=0 ; ii < nx ; ii++ ){ for( kk=0 ; kk < nz ; kk++ ){ sum = wt[0]*RR(ii,kk) ; for( qq=1 ; qq <= m ; qq++ ) sum += wt[qq] * ( RR(ii,kk+qq) + RR(ii,kk-qq) ) ; SS(ii,kk) = sum ; }} break ; case 7: for( ii=0 ; ii < nx ; ii++ ){ for( kk=0 ; kk < nz ; kk++ ){ SS(ii,kk) = wt7 * ( RR(ii,kk+7) + RR(ii,kk-7) ) + wt6 * ( RR(ii,kk+6) + RR(ii,kk-6) ) + wt5 * ( RR(ii,kk+5) + RR(ii,kk-5) ) + wt4 * ( RR(ii,kk+4) + RR(ii,kk-4) ) + wt3 * ( RR(ii,kk+3) + RR(ii,kk-3) ) + wt2 * ( RR(ii,kk+2) + RR(ii,kk-2) ) + wt1 * ( RR(ii,kk+1) + RR(ii,kk-1) ) + wt0 * RR(ii,kk) ; }} break ; case 6: for( ii=0 ; ii < nx ; ii++ ){ for( kk=0 ; kk < nz ; kk++ ){ SS(ii,kk) = wt6 * ( RR(ii,kk+6) + RR(ii,kk-6) ) + wt5 * ( RR(ii,kk+5) + RR(ii,kk-5) ) + wt4 * ( RR(ii,kk+4) + RR(ii,kk-4) ) + wt3 * ( RR(ii,kk+3) + RR(ii,kk-3) ) + wt2 * ( RR(ii,kk+2) + RR(ii,kk-2) ) + wt1 * ( RR(ii,kk+1) + RR(ii,kk-1) ) + wt0 * RR(ii,kk) ; }} break ; case 5: for( ii=0 ; ii < nx ; ii++ ){ for( kk=0 ; kk < nz ; kk++ ){ SS(ii,kk) = wt5 * ( RR(ii,kk+5) + RR(ii,kk-5) ) + wt4 * ( RR(ii,kk+4) + RR(ii,kk-4) ) + wt3 * ( RR(ii,kk+3) + RR(ii,kk-3) ) + wt2 * ( RR(ii,kk+2) + RR(ii,kk-2) ) + wt1 * ( RR(ii,kk+1) + RR(ii,kk-1) ) + wt0 * RR(ii,kk) ; }} break ; case 4: for( ii=0 ; ii < nx ; ii++ ){ for( kk=0 ; kk < nz ; kk++ ){ SS(ii,kk) = wt4 * ( RR(ii,kk+4) + RR(ii,kk-4) ) + wt3 * ( RR(ii,kk+3) + RR(ii,kk-3) ) + wt2 * ( RR(ii,kk+2) + RR(ii,kk-2) ) + wt1 * ( RR(ii,kk+1) + RR(ii,kk-1) ) + wt0 * RR(ii,kk) ; }} break ; case 3: for( ii=0 ; ii < nx ; ii++ ){ for( kk=0 ; kk < nz ; kk++ ){ SS(ii,kk) = wt3 * ( RR(ii,kk+3) + RR(ii,kk-3) ) + wt2 * ( RR(ii,kk+2) + RR(ii,kk-2) ) + wt1 * ( RR(ii,kk+1) + RR(ii,kk-1) ) + wt0 * RR(ii,kk) ; }} break ; case 2: for( ii=0 ; ii < nx ; ii++ ){ for( kk=0 ; kk < nz ; kk++ ){ SS(ii,kk) = wt2 * ( RR(ii,kk+2) + RR(ii,kk-2) ) + wt1 * ( RR(ii,kk+1) + RR(ii,kk-1) ) + wt0 * RR(ii,kk) ; }} break ; case 1: for( ii=0 ; ii < nx ; ii++ ){ for( kk=0 ; kk < nz ; kk++ ){ SS(ii,kk) = wt1 * ( RR(ii,kk+1) + RR(ii,kk-1) ) + wt0 * RR(ii,kk) ; }} break ; } /* end of special cases of m */ /* put SS array back into input 3D array; again, inner loop over ii for most contiguous access to f[] array */ for( kk=0 ; kk < nz ; kk++ ){ for( ii=0 ; ii < nx ; ii++ ) ff[ii+D*kk] = SS(ii,kk) ; } } /* end of loop over y-direction (jj) */ /*** finito, cara mia ***/ free((void *)ss) ; free((void *)rr) ; EXRETURN ; } /*-------------------------------------------------------------------*/ /*! Like EDIT_blur_volume(), using FIR only (no FFTs) and only for float arrays. ---------------------------------------------------------------------*/ void FIR_blur_volume( int nx, int ny, int nz, float dx, float dy, float dz, float *ffim , float sigma ) { if( ffim != NULL && sigma > 0.0 ) FIR_blur_volume_3d(nx,ny,nz, dx,dy,dz, ffim, sigma,sigma,sigma) ; } /*-------------------------------------------------------------------*/ /* Gaussian blur in real space, of a float volume; like EDIT_blur_volume_3d(), but using FIR only. - nx,ny,nz = dimensions of array - dx,dy,dz = grid step sizes - ffim = array - sigmax = stdev for blur along x; if 0, no blurring in x; etc. ---------------------------------------------------------------------*/ void FIR_blur_volume_3d( int nx, int ny, int nz, float dx, float dy, float dz, float *ffim , float sigmax, float sigmay, float sigmaz ) { int fir_m , ii ; float *fir_wt=NULL , fac ; double sfac = AFNI_numenv("AFNI_BLUR_FIRFAC") ; if( sfac < 2.0 ) sfac = 2.5 ; ENTRY("FIR_blur_volume_3d") ; if( ffim == NULL ) EXRETURN ; if( sigmax <= 0.0 && sigmay <= 0.0 && sigmaz <= 0.0 ) EXRETURN ; if( dx <= 0.0 ) dx = 1.0 ; if( dy <= 0.0 ) dy = dx ; if( dz <= 0.0 ) dz = dx ; /*-- blur along x --*/ if( sigmax > 0.0 && nx > 1 ){ fir_m = (int) ceil( sfac * sigmax / dx ) ; /* about the 5% level */ if( fir_m < 1 ) fir_m = 1 ; if( fir_m > nx/2 ) fir_m = nx/2 ; fir_wt = (float *)malloc(sizeof(float)*(fir_m+1)) ; fir_gaussian_load( fir_m , dx/sigmax , fir_wt ) ; fir_blurx( fir_m , fir_wt , nx,ny,nz , ffim ) ; } /*-- blur along y --*/ if( sigmay > 0.0 && ny > 1 ){ fir_m = (int) ceil( sfac * sigmay / dy ) ; if( fir_m < 1 ) fir_m = 1 ; if( fir_m > ny/2 ) fir_m = ny/2 ; fir_wt = (float *)realloc(fir_wt,sizeof(float)*(fir_m+1)) ; fir_gaussian_load( fir_m , dy/sigmay , fir_wt ) ; fir_blury( fir_m , fir_wt , nx,ny,nz , ffim ) ; } /*-- blur along z --*/ if( sigmaz > 0.0 && nz > 1 ){ fir_m = (int) ceil( sfac * sigmaz / dz ) ; if( fir_m < 1 ) fir_m = 1 ; if( fir_m > nz/2 ) fir_m = nz/2 ; fir_wt = (float *)realloc(fir_wt,sizeof(float)*(fir_m+1)) ; fir_gaussian_load( fir_m , dz/sigmaz , fir_wt ) ; fir_blurz( fir_m , fir_wt , nx,ny,nz , ffim ) ; } if( fir_wt != NULL ) free(fir_wt) ; EXRETURN ; } /*-----------------------------------------------------------------------------*/ /*! Blurring, applied to a 2D RGB image. */ MRI_IMAGE * mri_rgb_blur2D( float sig , MRI_IMAGE *im ) { MRI_IMARR *imar ; MRI_IMAGE *rim , *gim , *bim , *newim ; ENTRY("mri_rgb_blur2D") ; if( im == NULL || im->kind != MRI_rgb || sig <= 0.0f ) RETURN(NULL) ; imar = mri_rgb_to_3float(im) ; rim = IMARR_SUBIM(imar,0) ; gim = IMARR_SUBIM(imar,1) ; bim = IMARR_SUBIM(imar,2) ; FIR_blur_volume_3d( rim->nx,rim->ny,1 , 1.0f,1.0f,1.0f , MRI_FLOAT_PTR(rim) , sig,sig,0.0f ) ; FIR_blur_volume_3d( gim->nx,gim->ny,1 , 1.0f,1.0f,1.0f , MRI_FLOAT_PTR(gim) , sig,sig,0.0f ) ; FIR_blur_volume_3d( bim->nx,bim->ny,1 , 1.0f,1.0f,1.0f , MRI_FLOAT_PTR(bim) , sig,sig,0.0f ) ; newim = mri_3to_rgb(rim,gim,bim) ; MRI_COPY_AUX(newim,im) ; DESTROY_IMARR(imar) ; RETURN(newim) ; } /*-----------------------------------------------------------------------------*/ /*! Blurring, applied to a 2D byte image. */ MRI_IMAGE * mri_byte_blur2D( float sig , MRI_IMAGE *im ) { MRI_IMAGE *rim, *newim ; ENTRY("mri_byte_blur2D") ; if( im == NULL || im->kind != MRI_byte || sig <= 0.0f ) RETURN(NULL) ; rim = mri_to_mri( MRI_float, im ) ; FIR_blur_volume_3d( rim->nx,rim->ny,1 , 1.0f,1.0f,1.0f , MRI_FLOAT_PTR(rim) , sig,sig,0.0f ) ; newim = mri_to_mri( MRI_byte , rim ); MRI_COPY_AUX(newim,im) ; mri_free(rim); RETURN(newim) ; } /*-----------------------------------------------------------------------------*/ /*! Blurring, applied to a 2D float image. */ MRI_IMAGE * mri_float_blur2D( float sig , MRI_IMAGE *im ) { MRI_IMAGE *newim ; ENTRY("mri_float_blur2D") ; if( im == NULL || im->kind != MRI_float || sig <= 0.0f ) RETURN(NULL) ; newim = mri_to_float(im) ; FIR_blur_volume_3d( newim->nx,newim->ny,1 , 1.0f,1.0f,1.0f , MRI_FLOAT_PTR(newim) , sig,sig,0.0f ) ; RETURN(newim) ; } /*-----------------------------------------------------------------------------*/ /*! Blurring, applied to a 3D float image (sig is in voxels, not mm). */ MRI_IMAGE * mri_float_blur3D( float sig , MRI_IMAGE *im ) { MRI_IMAGE *newim ; ENTRY("mri_float_blur3D") ; if( im == NULL || im->kind != MRI_float || sig <= 0.0f ) RETURN(NULL) ; newim = mri_to_float(im) ; FIR_blur_volume_3d( newim->nx,newim->ny,newim->nz , 1.0f,1.0f,1.0f , MRI_FLOAT_PTR(newim) , sig,sig,sig ) ; RETURN(newim) ; } /*-----------------------------------------------------------------------------*/ #undef NG #define NG 11 static void fir_gaussian_load( int m, float dx, float *ff ) { int ii ; float fac ; double xx ; if( m < 0 || dx <= 0.0f || ff == NULL ) return ; if(DO_INFO) INFO_message("fir_gaussian_load(%d,%g)",m,dx) ; if( AFNI_yesenv("AFNI_BLUR_FIROLD") ){ /** the olde waye **/ fac = ff[0] = 1.0f ; for( ii=1 ; ii <= m ; ii++ ){ xx = ii*dx ; ff[ii] = exp(-0.5*xx*xx ) ; /* center of each cell */ fac += 2.0f * ff[ii] ; } } else { /** the newe bettere waye **/ double sum , dxj , ee ; int jj ; fac = 0.0f ; dxj = dx/(2.0*NG) ; /* subdivision of each cell */ for( ii=0 ; ii <= m ; ii++ ){ sum = 0.0 ; for( jj=-NG ; jj <= NG ; jj++ ){ /* sum over subdivisions */ xx = ii*dx + jj*dxj ; ee = exp( -0.5*xx*xx ) ; if( jj==-NG || jj==NG ) sum += 0.5*ee ; else sum += ee ; } ff[ii] = sum ; if( ii==0 ) fac += sum ; else fac += 2.0*sum ; } } fac = 1.0f / fac ; for( ii=0 ; ii <= m ; ii++ ) ff[ii] *= fac ; if(DO_INFO) { fprintf(stderr,"-- fir weights:\n"); for( ii=0 ; ii <= m ; ii++ ) fprintf(stderr," %g\n", ff[ii]); } return ; } /**************************************************************************/ /**************************************************************************/ /**************************************************************************/ #if 0 /*------------------------------------------------------------------------*/ /* Below is a test program that can be used to evaluate the relative speed of FFT and FIR blurring. If compiled into 'tblur', then a command line would be 'tblur 240 13.0' meaning blur a 240x240x240 volume with sigma=13.0. The output is a line indicating the input parameters and the CPU time ratio (FFT time)/(FIR time). Compilation should be something like so: make tblur.o cc -o tblur tblur.o -L. -L/usr/X11R6/lib -lmri -lXt -lm --------------------------------------------------------------------------*/ #include "mrilib.h" int main( int argc , char *argv[] ) { int nx , nxxx , ii , jj ; float sigma , *gar , *far ; double c1,c2,c3 , cg,cf ; EDIT_blur_allow_fir(0) ; /* turn off auto-FIR in EDIT_blur_volume() */ /* get command line params */ if( argc < 3 ) exit(0) ; nx = (int)strtod( argv[1] , NULL) ; if( nx < 32 ) nx = 128 ; sigma = (float)strtod( argv[2] , NULL ) ; if( sigma <= 0.0 ) exit(0) ; /* allocate volumes */ nxxx = nx*nx*nx ; gar = (float *)malloc( sizeof(float)*nxxx ) ; far = (float *)malloc( sizeof(float)*nxxx ) ; /* load volume for FFT and then blur 5 times */ for( ii=0 ; ii < nxxx ; ii++ ) gar[ii] = cos(0.37382*ii) ; c1 = COX_cpu_time() ; for( jj=0 ; jj < 5 ; jj++ ) EDIT_blur_volume( nx,nx,nx , 1.0,1.0,1.0 , MRI_float,gar , sigma ) ; c2 = COX_cpu_time() ; cg = c2-c1 ; /** printf("nx=%d sigma=%.2f Gaussian cpu=%.3f\n",nx,sigma,cg) ; **/ /* load volume for FIR and then blur 5 times */ for( ii=0 ; ii < nxxx ; ii++ ) far[ii] = cos(0.37382*ii) ; c1 = COX_cpu_time() ; for( jj=0 ; jj < 5 ; jj++ ) FIR_blur_volume( nx,nx,nx , 1.0,1.0,1.0 , far , sigma ) ; c2 = COX_cpu_time() ; cf = c2-c1 ; /** printf("nx=%d sigma=%.2f FIR_blur cpu=%.3f\n",nx,sigma,cf) ; **/ /* output CPU time ratio */ printf("ratio: %d %.2f %.2f\n",nx,sigma,cg/cf) ; /* compute mean difference between the 2 approaches */ /** c1 = 0.0 ; for( ii=0 ; ii < nxxx ; ii++ ) c1 += fabs(far[ii]-gar[ii]) ; c1 = c1 / nxxx ; printf("mean abs diff = %g\n",c1) ; **/ exit(0) ; } #endif