1 #include "mrilib.h"
2
3 #ifdef USE_OMP
4 #include <omp.h>
5 #endif
6
7 #undef INMASK
8 #define INMASK(i) (mask == NULL || mask[i] != 0)
9
10 /*----------------------------------------------------------------------------*/
11 /*! Blur image, in place, confined to a mask, with blurring factors given
12 separately for each dimension xyz. A NULL blurring factor means don't
13 blur in that direction. Blurring factors can be thought of in 2 ways
14 - diffusion equation: fx = 0.5 * Delta_t * D_x / Delta_x**2
15 - FWHM blurring: fx = (Delta_FWHMx / Delta_x)**2 * 0.045084
16 - The fx (etc.) factors should be between 0 and 0.05 for stability;
17 for increasing the FWHM of the image, this means that the maximum
18 change in FWHM is about Delta_x in each step.
19 - Not all input fx,fy,fz factors should be NULL (or nothing will happen).
20 - Method: 3 point stencil (in each direction) conservative finite
21 difference approximation to du/dt = d/dx[ D_x(x,y,z) du/dx ] + ...
22 - Points not in the mask are not processed, and will be set to zero
23 in the output.
24 - The input mask can be NULL. If you really want speed, a special
25 version should be written for the mask=NULL case. Please send
26 pumpernickel bagels or dark chocolate covered cranberries.
27 - Author: Zhark, Emperor of the Galaxy!!! (Nov 2006)
28 ------------------------------------------------------------------------------*/
29
mri_blur3D_variable(MRI_IMAGE * im,byte * mask,MRI_IMAGE * fx,MRI_IMAGE * fy,MRI_IMAGE * fz)30 void mri_blur3D_variable( MRI_IMAGE *im , byte *mask ,
31 MRI_IMAGE *fx , MRI_IMAGE *fy , MRI_IMAGE *fz )
32 {
33 int nx,ny,nz,nxy,nxyz ;
34 float *iar , *fxar , *fyar , *fzar , *qar ;
35 int ijk , ii,jj,kk ;
36 float vcc , vsub , vout , vf ;
37
38 ENTRY("mri_blur3D_variable") ;
39
40 if( im == NULL ) EXRETURN ;
41 if( fx == NULL && fy == NULL && fz == NULL ) EXRETURN ;
42
43 nx = im->nx; ny = im->ny; nz = im->nz; nxy = nx*ny; nxyz = nxy*nz;
44
45 iar = MRI_FLOAT_PTR(im) ;
46 fxar = MRI_FLOAT_PTR(fx) ;
47 fyar = MRI_FLOAT_PTR(fy) ;
48 fzar = MRI_FLOAT_PTR(fz) ;
49 qar = (float *)calloc(sizeof(float),nxyz) ;
50
51 for( ijk=kk=0 ; kk < nz ; kk++ ){
52 for( jj=0 ; jj < ny ; jj++ ){
53 for( ii=0 ; ii < nx ; ii++,ijk++ ){
54 if( !INMASK(ijk) ) continue ;
55 vout = vcc = iar[ijk] ;
56 if( fxar != NULL ){ /* distribute (diffuse) in the x-direction */
57 vf = fxar[ijk] ;
58 if( ii-1 >= 0 && INMASK(ijk-1) ){
59 vsub = (fxar[ijk-1]+vf)*vcc; qar[ijk-1] += vsub; vout -= vsub;
60 }
61 if( ii+1 < nx && INMASK(ijk+1) ){
62 vsub = (fxar[ijk+1]+vf)*vcc; qar[ijk+1] += vsub; vout -= vsub;
63 }
64 }
65 if( fyar != NULL ){ /* distribute (diffuse) in the y-direction */
66 vf = fyar[ijk] ;
67 if( jj-1 >= 0 && INMASK(ijk-nx) ){
68 vsub = (fyar[ijk-nx]+vf)*vcc; qar[ijk-nx] += vsub; vout -= vsub;
69 }
70 if( jj+1 < ny && INMASK(ijk+nx) ){
71 vsub = (fyar[ijk+nx]+vf)*vcc; qar[ijk+nx] += vsub; vout -= vsub;
72 }
73 }
74 if( fzar != NULL ){ /* distribute (diffuse) in the z-direction */
75 vf = fzar[ijk] ;
76 if( kk-1 >= 0 && INMASK(ijk-nxy) ){
77 vsub = (fzar[ijk-nxy]+vf)*vcc; qar[ijk-nxy] += vsub; vout -= vsub;
78 }
79 if( kk+1 < nz && INMASK(ijk+nxy) ){
80 vsub = (fzar[ijk+nxy]+vf)*vcc; qar[ijk+nxy] += vsub; vout -= vsub;
81 }
82 }
83
84 qar[ijk] += vout ; /* whatever wasn't diffused away from this voxel */
85 }}}
86
87 AAmemcpy(iar,qar,sizeof(float)*nxyz) ;
88 free((void *)qar) ;
89 EXRETURN ;
90 }
91
92 /*----------------------------------------------------------------------------*/
93 /*! Similar to the above, but with fixed blurring factors in each direction.
94 But inside a mask. Again, 0 <= fx <= 0.05 for stability. [01 May 2009]
95 *//*--------------------------------------------------------------------------*/
96
mri_blur3D_inmask(MRI_IMAGE * im,byte * mask,float fx,float fy,float fz,int nrep)97 void mri_blur3D_inmask( MRI_IMAGE *im, byte *mask,
98 float fx,float fy,float fz, int nrep )
99 {
100 int nx,ny,nz,nxy,nxyz ;
101 float *iar , *qar ;
102 int ijk , ii,jj,kk , nn, nfloat_err =0 ;
103 register float vcc , vsub , vout , vx,vy,vz ;
104
105 ENTRY("mri_blur3D_inmask") ;
106
107 if( im == NULL || nrep <= 0 ) EXRETURN ;
108
109 nx = im->nx; ny = im->ny; nz = im->nz; nxy = nx*ny; nxyz = nxy*nz;
110
111 iar = MRI_FLOAT_PTR(im) ;
112 vx = 2.0f * fx ; if( nx < 2 ) vx = 0.0f ;
113 vy = 2.0f * fy ; if( ny < 2 ) vy = 0.0f ;
114 vz = 2.0f * fz ; if( nz < 2 ) vz = 0.0f ;
115 if( vx <= 0.0f && vy <= 0.0f && vz <= 0.0f ) EXRETURN ;
116
117 #pragma omp critical (MALLOC)
118 qar = (float *)calloc(sizeof(float),nxyz) ;
119
120 for( nn=0 ; nn < nrep ; nn++ ){
121 for( ijk=kk=0 ; kk < nz ; kk++ ){
122 for( jj=0 ; jj < ny ; jj++ ){
123 for( ii=0 ; ii < nx ; ii++,ijk++ ){
124 if( !INMASK(ijk) ) continue ;
125 vout = vcc = iar[ijk] ;
126 if( vx > 0.0f ){ /* distribute (diffuse) in the x-direction */
127 if( ii-1 >= 0 && INMASK(ijk-1) ){
128 vsub = vx*vcc; qar[ijk-1] += vsub; vout -= vsub;
129 }
130 if( ii+1 < nx && INMASK(ijk+1) ){
131 vsub = vx*vcc; qar[ijk+1] += vsub; vout -= vsub;
132 }
133 }
134 if( vy > 0.0f ){ /* distribute (diffuse) in the y-direction */
135 if( jj-1 >= 0 && INMASK(ijk-nx) ){
136 vsub = vy*vcc; qar[ijk-nx] += vsub; vout -= vsub;
137 }
138 if( jj+1 < ny && INMASK(ijk+nx) ){
139 vsub = vy*vcc; qar[ijk+nx] += vsub; vout -= vsub;
140 }
141 }
142 if( vz >= 0.0f ){ /* distribute (diffuse) in the z-direction */
143 if( kk-1 >= 0 && INMASK(ijk-nxy) ){
144 vsub = vz*vcc; qar[ijk-nxy] += vsub; vout -= vsub;
145 }
146 if( kk+1 < nz && INMASK(ijk+nxy) ){
147 vsub = vz*vcc; qar[ijk+nxy] += vsub; vout -= vsub;
148 }
149 }
150
151 qar[ijk] += vout ; /* whatever wasn't diffused away from this voxel */
152 }}}
153 AAmemcpy(iar,qar,sizeof(float)*nxyz) ;
154 if( nn != nrep-1 ){
155 AAmemset(qar,0,sizeof(float)*nxyz) ;
156 }
157 }
158
159 #pragma omp critical (MALLOC)
160 free((void *)qar) ;
161 EXRETURN ;
162 }
163
164 /*! This function can be a slightly faster (20%) version of mri_blur3D_inmask
165 under the following conditions:
166 1- The mask is big and has few or no holes in it, as a brain mask would
167 2- You have a volume with more than 2 voxels in each direction and you
168 are blurring in 3D
169 ZSS March 2011
170 */
mri_blur3D_inmask_speedy(MRI_IMAGE * im,byte * mask,float fx,float fy,float fz,int nrep)171 void mri_blur3D_inmask_speedy( MRI_IMAGE *im, byte *mask,
172 float fx,float fy,float fz, int nrep )
173 {
174 int nx,ny,nz,nxy,nxyz ;
175 float *iar , *qar ;
176 int ijk , ii,jj,kk , nn, ijkm, nfloat_err=0 ;
177 byte *skin = NULL;
178 register float vcc , vsub , vout , vx,vy,vz ;
179
180 ENTRY("mri_blur3D_inmask_speedy") ;
181
182 if( im == NULL || nrep <= 0 ) EXRETURN ;
183
184 nx = im->nx; ny = im->ny; nz = im->nz; nxy = nx*ny; nxyz = nxy*nz;
185
186 iar = MRI_FLOAT_PTR(im) ;
187 vx = 2.0f * fx ; if( nx < 2 ) vx = 0.0f ;
188 vy = 2.0f * fy ; if( ny < 2 ) vy = 0.0f ;
189 vz = 2.0f * fz ; if( nz < 2 ) vz = 0.0f ;
190 if( vx <= 0.0f || vy <= 0.0f || vz <= 0.0f ) {
191 ERROR_message("Cannot handle cases where v* <= 0.0");
192 EXRETURN ;
193 }
194
195 #pragma omp critical (MALLOC)
196 skin = (byte *)calloc(sizeof(byte), nxyz);
197 ijkm = nxyz-nxy-1;
198 for( ijk=kk=0 ; kk < nz ; kk++ ){
199 for( jj=0 ; jj < ny ; jj++ ){
200 for( ii=0 ; ii < nx ; ii++,ijk++ ){
201 if( !INMASK(ijk) ) continue ;
202 if ( ii == 0 || jj == 0 || kk == 0 ||
203 ii == nx-1 || jj == ny-1 || jj == nz-1 ||
204 ijk < nxy || ijk > ijkm) {
205 skin[ijk] = 1; /* in mask, on edge of volume,
206 or close to boundary slices */
207 } else if ( !INMASK(ijk-1) || !INMASK(ijk+1) ||
208 !INMASK(ijk-nx) || !INMASK(ijk+nx)||
209 !INMASK(ijk-nxy)|| !INMASK(ijk+nxy) ){
210 skin[ijk] = 1; /* on edge of mask */
211 }
212 }
213 }
214 }
215
216 #pragma omp critical (MALLOC)
217 qar = (float *)calloc(sizeof(float),nxyz) ;
218
219 for( nn=0 ; nn < nrep ; nn++ ){
220 for( ijk=kk=0 ; kk < nz ; kk++ ){
221 for( jj=0 ; jj < ny ; jj++ ){
222 for( ii=0 ; ii < nx ; ii++,ijk++ ){
223 if( !INMASK(ijk) ) continue ;
224 vout = vcc = iar[ijk] ;
225 if (!skin[ijk]) { /* Not skin, go fast */
226 { /* distribute (diffuse) in the x-direction */
227 {
228 vsub = vx*vcc; qar[ijk-1] += vsub; vout -= vsub;
229 }
230 {
231 vsub = vx*vcc; qar[ijk+1] += vsub; vout -= vsub;
232 }
233 }
234 { /* distribute (diffuse) in the y-direction */
235 {
236 vsub = vy*vcc; qar[ijk-nx] += vsub; vout -= vsub;
237 }
238 {
239 vsub = vy*vcc; qar[ijk+nx] += vsub; vout -= vsub;
240 }
241 }
242 { /* distribute (diffuse) in the z-direction */
243 {
244 vsub = vz*vcc; qar[ijk-nxy] += vsub; vout -= vsub;
245 }
246 {
247 vsub = vz*vcc; qar[ijk+nxy] += vsub; vout -= vsub;
248 }
249 }
250 } else { /* skin voxel, go slow */
251 { /* distribute (diffuse) in the x-direction */
252 if( ii-1 >= 0 && INMASK(ijk-1) ){
253 vsub = vx*vcc; qar[ijk-1] += vsub; vout -= vsub;
254 }
255 if( ii+1 < nx && INMASK(ijk+1) ){
256 vsub = vx*vcc; qar[ijk+1] += vsub; vout -= vsub;
257 }
258 }
259 { /* distribute (diffuse) in the y-direction */
260 if( jj-1 >= 0 && INMASK(ijk-nx) ){
261 vsub = vy*vcc; qar[ijk-nx] += vsub; vout -= vsub;
262 }
263 if( jj+1 < ny && INMASK(ijk+nx) ){
264 vsub = vy*vcc; qar[ijk+nx] += vsub; vout -= vsub;
265 }
266 }
267 { /* distribute (diffuse) in the z-direction */
268 if( kk-1 >= 0 && INMASK(ijk-nxy) ){
269 vsub = vz*vcc; qar[ijk-nxy] += vsub; vout -= vsub;
270 }
271 if( kk+1 < nz && INMASK(ijk+nxy) ){
272 vsub = vz*vcc; qar[ijk+nxy] += vsub; vout -= vsub;
273 }
274 }
275 }
276 qar[ijk] += vout ; /* whatever wasn't diffused away from this voxel */
277 }}}
278
279 AAmemcpy(iar,qar,sizeof(float)*nxyz) ;
280 if( nn != nrep-1 ){
281 AAmemset(qar,0,sizeof(float)*nxyz) ;
282 }
283 }
284
285 #pragma omp critical (MALLOC)
286 free((void *)qar) ;
287 #pragma omp critical (MALLOC)
288 free((void *)skin) ;
289 EXRETURN ;
290 }
291
292 /*
293 Nearest neighbor averaging, for faster but poorly titrated smoothing.
294 ZSS Nov 2011
295 */
mri_blur3D_inmask_NN(MRI_IMAGE * im,byte * mask,int nrep)296 void mri_blur3D_inmask_NN( MRI_IMAGE *im, byte *mask, int nrep )
297 {
298 int nx,ny,nz,nxy,nxyz ;
299 float *iar , *qar ;
300 int ijk , ii,jj,kk , nn, ijkm, nfloat_err=0 ;
301 byte *skin = NULL;
302 register float vout ;
303
304 ENTRY("mri_blur3D_inmask_NN") ;
305
306 if( im == NULL || nrep <= 0 ) EXRETURN ;
307
308 nx = im->nx; ny = im->ny; nz = im->nz; nxy = nx*ny; nxyz = nxy*nz;
309
310 iar = MRI_FLOAT_PTR(im) ;
311
312 #pragma omp critical (MALLOC)
313 skin = (byte *)calloc(sizeof(byte), nxyz);
314 ijkm = nxyz-nxy-1;
315 for( ijk=kk=0 ; kk < nz ; kk++ ){
316 for( jj=0 ; jj < ny ; jj++ ){
317 for( ii=0 ; ii < nx ; ii++,ijk++ ){
318 if( !INMASK(ijk) ) continue ;
319 if ( ii == 0 || jj == 0 || kk == 0 ||
320 ii == nx-1 || jj == ny-1 || jj == nz-1 ||
321 ijk < nxy || ijk > ijkm) {
322 skin[ijk] = 1; /* in mask, on edge of volume,
323 or close to boundary slices */
324 } else if ( !INMASK(ijk-1) || !INMASK(ijk+1) ||
325 !INMASK(ijk-nx) || !INMASK(ijk+nx)||
326 !INMASK(ijk-nxy)|| !INMASK(ijk+nxy) ){
327 skin[ijk] = 1; /* on edge of mask */
328 }
329 }
330 }
331 }
332
333 #pragma omp critical (MALLOC)
334 qar = (float *)calloc(sizeof(float),nxyz) ;
335
336 for( nn=0 ; nn < nrep ; nn++ ){
337 for( ijk=kk=0 ; kk < nz ; kk++ ){
338 for( jj=0 ; jj < ny ; jj++ ){
339 for( ii=0 ; ii < nx ; ii++,ijk++ ){
340 if( !INMASK(ijk) ) continue ;
341 vout = iar[ijk] ;
342 if (!skin[ijk]) { /* Not skin, go fast */
343 qar[ijk] = (iar[ijk] +
344 iar[ijk-1] + iar[ijk+1] +
345 iar[ijk-nx] + iar[ijk+nx] +
346 iar[ijk-nxy] + iar[ijk+nxy] ) / 7.0;
347 } else { /* skin voxel, go slow */
348 vout = 1.0;
349 qar[ijk] = iar[ijk];
350 if( ii-1 >= 0 && INMASK(ijk-1) ){
351 qar[ijk] += qar[ijk-1]; ++vout;
352 }
353 if( ii+1 < nx && INMASK(ijk+1) ){
354 qar[ijk] += qar[ijk+1]; ++vout;
355 }
356 if( jj-1 >= 0 && INMASK(ijk-nx) ){
357 qar[ijk] += qar[ijk-nx]; ++vout;
358 }
359 if( jj+1 < ny && INMASK(ijk+nx) ){
360 qar[ijk] += qar[ijk+nx]; ++vout;
361 }
362 if( kk-1 >= 0 && INMASK(ijk-nxy) ){
363 qar[ijk] += qar[ijk-nxy]; ++vout;
364 }
365 if( kk+1 < nz && INMASK(ijk+nxy) ){
366 qar[ijk] += qar[ijk+nxy]; ++vout;
367 }
368 qar[ijk] /= vout;
369 }
370 }}}
371
372 AAmemcpy(iar,qar,sizeof(float)*nxyz) ;
373 }
374
375 #pragma omp critical (MALLOC)
376 free((void *)qar) ;
377 #pragma omp critical (MALLOC)
378 free((void *)skin) ;
379 EXRETURN ;
380 }
381
382 /*----------------------------------------------------------------------------*/
383
384 #undef DMAX
385 #define DMAX 999999.9f
386
mri_blur3D_getfac(float fwhm,float dx,float dy,float dz,int * nrep,float * fx,float * fy,float * fz)387 void mri_blur3D_getfac( float fwhm , float dx, float dy, float dz ,
388 int *nrep , float *fx , float *fy , float *fz )
389 {
390 float dm = DMAX ;
391
392 if( fwhm <= 0.0f ) return ;
393
394 if( dx > 0.0f ) dm = dx ;
395 if( dy > 0.0f ) dm = MIN(dm,dy) ;
396 if( dz > 0.0f ) dm = MIN(dm,dz) ;
397 if( dm == DMAX ) return ;
398
399 *nrep = 2 + (int)(fwhm*fwhm/(dm*dm)) ;
400 if( dx > 0.0f ) *fx = 0.045084f*(fwhm*fwhm)/(*nrep*dx*dx) ;
401 else *fx = 0.0f ;
402 if( dy > 0.0f ) *fy = 0.045084f*(fwhm*fwhm)/(*nrep*dy*dy) ;
403 else *fy = 0.0f ;
404 if( dz > 0.0f ) *fz = 0.045084f*(fwhm*fwhm)/(*nrep*dz*dz) ;
405 else *fz = 0.0f ;
406
407 return ;
408 }
409
410 /*----------------------------------------------------------------------------*/
411
mri_blur3D_addfwhm(MRI_IMAGE * im,byte * mask,float fwhm)412 void mri_blur3D_addfwhm( MRI_IMAGE *im , byte *mask , float fwhm )
413 {
414 int nrep=-1 ;
415 float fx=-1.0f,fy=-1.0f,fz=-1.0f , dx,dy,dz ;
416
417 ENTRY("mri_blur3D_addfwhm") ;
418
419 if( im == NULL || fwhm <= 0.0f ) EXRETURN ;
420
421 dx = im->dx; if( dx == 0.0f ) dx = 1.0f; else if( dx < 0.0f ) dx = -dx;
422 dy = im->dy; if( dy == 0.0f ) dy = 1.0f; else if( dy < 0.0f ) dy = -dy;
423 dz = im->dz; if( dz == 0.0f ) dz = 1.0f; else if( dz < 0.0f ) dz = -dz;
424
425 mri_blur3D_getfac( fwhm , dx,dy,dz , &nrep,&fx,&fy,&fz ) ;
426
427 if( 1 || MRILIB_verb )
428 INFO_message("mri_blur3D: #iter=%d fx=%.5f fy=%.5f fz=%.5f",nrep,fx,fy,fz) ;
429 if( nrep < 0 || fx < 0.0f || fy < 0.0f || fz < 0.0f ) EXRETURN ;
430
431 mri_blur3D_inmask( im , mask , fx,fy,fz , nrep ) ;
432
433 EXRETURN ;
434 }
435
436 /*----------------------------------------------------------------------------*/
437 /* 24 Mar 2021 */
438
mri_blur3D_getfac3(float fwhmx,float fwhmy,float fwhmz,float dx,float dy,float dz,int * nrep,float * fx,float * fy,float * fz)439 void mri_blur3D_getfac3( float fwhmx, float fwhmy, float fwhmz ,
440 float dx , float dy , float dz ,
441 int *nrep , float *fx , float *fy , float *fz )
442 {
443 float fxx =0.0f , fyy =0.0f , fzz =0.0f ;
444 int nrepx=0 , nrepy=0 , nrepz=0 , nrepout=0 ;
445
446 if( fwhmx > 0.0f && dx > 0.0f ) nrepx = 2 + (int)(fwhmx*fwhmx/(dx*dx)) ;
447 if( fwhmy > 0.0f && dy > 0.0f ) nrepy = 2 + (int)(fwhmy*fwhmy/(dy*dy)) ;
448 if( fwhmz > 0.0f && dz > 0.0f ) nrepz = 2 + (int)(fwhmz*fwhmz/(dz*dz)) ;
449
450 if( nrepx > nrepout ) nrepout = nrepx ;
451 if( nrepy > nrepout ) nrepout = nrepy ;
452 if( nrepz > nrepout ) nrepout = nrepz ;
453
454 if( fwhmx > 0.0f && dx > 0.0f ) fxx = 0.045084f*(fwhmx*fwhmx)/(nrepout*dx*dx) ;
455 if( fwhmy > 0.0f && dy > 0.0f ) fyy = 0.045084f*(fwhmy*fwhmy)/(nrepout*dy*dy) ;
456 if( fwhmz > 0.0f && dz > 0.0f ) fzz = 0.045084f*(fwhmz*fwhmz)/(nrepout*dz*dz) ;
457
458 *fx = fxx ; *fy = fyy ; *fz = fzz ; *nrep = nrepout ;
459 return ;
460 }
461
462 /*----------------------------------------------------------------------------*/
463 /* 24 Mar 2021 */
464
mri_blur3D_addfwhm3(MRI_IMAGE * im,byte * mask,float fwhmx,float fwhmy,float fwhmz)465 void mri_blur3D_addfwhm3( MRI_IMAGE *im , byte *mask , float fwhmx,float fwhmy,float fwhmz )
466 {
467 int nrep=-1 ;
468 float fx=-1.0f,fy=-1.0f,fz=-1.0f , dx,dy,dz ;
469
470 ENTRY("mri_blur3D_addfwhm3") ;
471
472 if( im == NULL ) EXRETURN ;
473
474 #if 0
475 if( fwhmx == fwhmy && fwhmx == fwhmz ){
476 mri_blur3D_addfwhm( im , mask , fwhmx ) ;
477 EXRETURN ;
478 }
479 #endif
480
481 dx = im->dx; if( dx == 0.0f ) dx = 1.0f; else if( dx < 0.0f ) dx = -dx;
482 dy = im->dy; if( dy == 0.0f ) dy = 1.0f; else if( dy < 0.0f ) dy = -dy;
483 dz = im->dz; if( dz == 0.0f ) dz = 1.0f; else if( dz < 0.0f ) dz = -dz;
484
485 mri_blur3D_getfac3( fwhmx,fwhmy,fwhmz , dx,dy,dz , &nrep,&fx,&fy,&fz ) ;
486
487 if( 1 || MRILIB_verb )
488 INFO_message("mri_blur3D: #iter=%d fx=%.5f fy=%.5f fz=%.5f",nrep,fx,fy,fz) ;
489 if( nrep <= 0 || (fx <= 0.0f && fy <= 0.0f && fz <= 0.0f) ) EXRETURN ;
490
491 mri_blur3D_inmask( im , mask , fx,fy,fz , nrep ) ;
492
493 EXRETURN ;
494 }
495
496 /*----------------------------------------------------------------------------*/
497 /*! A version of mri_blur3D_addfwhm that can be about 20%faster under certain
498 conditions. See mri_blur3D_inmask_speedy for details.
499 ZSS March 2011
500 */
mri_blur3D_addfwhm_speedy(MRI_IMAGE * im,byte * mask,float fwhm)501 void mri_blur3D_addfwhm_speedy( MRI_IMAGE *im , byte *mask , float fwhm )
502 {
503 int nrep=-1 ;
504 float fx=-1.0f,fy=-1.0f,fz=-1.0f , dx,dy,dz ;
505
506 ENTRY("mri_blur3D_addfwhm_speedy") ;
507
508 if( im == NULL || fwhm <= 0.0f ) EXRETURN ;
509
510 dx = im->dx; if( dx == 0.0f ) dx = 1.0f; else if( dx < 0.0f ) dx = -dx;
511 dy = im->dy; if( dy == 0.0f ) dy = 1.0f; else if( dy < 0.0f ) dy = -dy;
512 dz = im->dz; if( dz == 0.0f ) dz = 1.0f; else if( dz < 0.0f ) dz = -dz;
513
514 mri_blur3D_getfac( fwhm , dx,dy,dz , &nrep,&fx,&fy,&fz ) ;
515 if( nrep < 0 || fx < 0.0f || fy < 0.0f || fz < 0.0f ) EXRETURN ;
516
517 if( MRILIB_verb )
518 INFO_message("mri_blur3D: #iter=%d fx=%.5f fy=%.5f fz=%.5f",nrep,fx,fy,fz) ;
519
520 if (fx > 0.0 && fy > 0.0 && fz > 0.0 &&
521 im->nx > 2 && im->ny > 2 && im->nz > 2) {
522 mri_blur3D_inmask_speedy( im , mask , fx,fy,fz , nrep ) ;
523 } else {
524 INFO_message("mri_blur3D_addfwhm_speedy:\n"
525 " Thin volume or 2D blurring, Going the slow route.");
526 mri_blur3D_inmask( im , mask , fx,fy,fz , nrep ) ;
527 }
528 EXRETURN ;
529 }
530
531
532 #ifndef MRILIB_MINI
533 /*----------------------------------------------------------------------------*/
534 /*! Blur a vectim, by converting each time point to a 3D image,
535 blurring that image, and then putting it back into the vectim.
536 *//*--------------------------------------------------------------------------*/
537
mri_blur3D_vectim(MRI_vectim * vim,float fwhm)538 void mri_blur3D_vectim( MRI_vectim *vim , float fwhm )
539 {
540 int nrep=-1 , nx,ny,nz , nvox , kk ;
541 float fx=-1.0f,fy=-1.0f,fz=-1.0f , dx,dy,dz ;
542 int *ivar ; byte *mmm ;
543
544 ENTRY("mri_blur3d_vectim") ;
545
546 if( vim == NULL || fwhm <= 0.0f ) EXRETURN ;
547
548 dx = vim->dx ; if( dx == 0.0f ) dx = 1.0f; else if( dx < 0.0f ) dx = -dx;
549 dy = vim->dy ; if( dy == 0.0f ) dy = 1.0f; else if( dy < 0.0f ) dy = -dy;
550 dz = vim->dz ; if( dz == 0.0f ) dz = 1.0f; else if( dz < 0.0f ) dz = -dz;
551
552 nx = vim->nx ; ny = vim->ny ; nz = vim->nz ; nvox = nx*ny*nz ;
553 if( nx < 1 || ny < 1 || nz < 1 ) EXRETURN ;
554
555 mri_blur3D_getfac( fwhm , dx,dy,dz , &nrep,&fx,&fy,&fz ) ;
556 if( nrep < 0 || fx < 0.0f || fy < 0.0f || fz < 0.0f ) EXRETURN ;
557
558 if( MRILIB_verb )
559 INFO_message("mri_blur3D: #iter=%d fx=%.5f fy=%.5f fz=%.5f",nrep,fx,fy,fz) ;
560
561 ivar = vim->ivec ;
562 mmm = (byte *)calloc(sizeof(byte),nvox) ;
563 for( kk=0 ; kk < vim->nvec ; kk++ ) mmm[ivar[kk]] = 1 ;
564
565 AFNI_OMP_START ;
566 #pragma omp parallel if( vim->nvals > 1 )
567 {
568 MRI_IMAGE *qim ;
569 float *qar , *var ;
570 int iv , jj ;
571
572 #pragma omp critical (BLUR3D_vectim)
573 { qim = mri_new_vol( nx,ny,nz , MRI_float ) ; qar = MRI_FLOAT_PTR(qim) ; }
574
575 #ifdef USE_OMP
576 if( MRILIB_verb ) ININFO_message(" start OpenMP thread %d",omp_get_thread_num());
577 #endif
578
579 #pragma omp for
580 for( iv=0 ; iv < vim->nvals ; iv++ ){
581 AAmemset( qar , 0 , sizeof(float)*nvox ) ;
582 for( jj=0 ; jj < vim->nvec ; jj++ ){
583 var = VECTIM_PTR(vim,jj) ; qar[ivar[jj]] = var[iv] ;
584 }
585 mri_blur3D_inmask( qim , mmm , fx,fy,fz , nrep ) ;
586 for( jj=0 ; jj < vim->nvec ; jj++ ){
587 var = VECTIM_PTR(vim,jj) ; var[iv] = qar[ivar[jj]] ;
588 }
589 }
590
591 #pragma omp critical (BLUR3D_vectim)
592 { mri_free(qim) ; }
593
594 } /* end OpenMP */
595 AFNI_OMP_END ;
596
597 free(mmm) ; EXRETURN ;
598 }
599 #endif /* MRILIB_MINI */
600