TPCCLIB
Loading...
Searching...
No Matches
imgfilter.c
Go to the documentation of this file.
1
4#include "libtpcimgp.h"
5/*****************************************************************************/
6
7/*****************************************************************************/
18int imgGaussianFIRFilter(
21 IMG *img,
23 float xsd,
25 float ysd,
27 float zsd,
29 double tolerance,
31 int verbose
32) {
33 if(verbose>0) printf("%s(*img, %g, %g, %g, %g, %d)\n", __func__, xsd, ysd, zsd, tolerance, verbose);
34
35 if(img==NULL || img->dimz<1 || img->dimy<1 || img->dimx<1 || img->dimt<1) return 1;
36 int dimx=img->dimx, dimy=img->dimy, dimz=img->dimz, dimt=img->dimt;
37 if(dimx==1) xsd=0.0;
38 if(dimy==1) ysd=0.0;
39 if(dimz==1) zsd=0.0;
40 if(xsd>0.0 && (dimx<4 || !(img->sizex>0.0))) return(2);
41 if(ysd>0.0 && (dimy<4 || !(img->sizey>0.0))) return(3);
42 if(zsd>0.0 && (dimz<4 || !(img->sizez>0.0))) return(4);
43 if(!(xsd>0.0) && !(ysd>0.0) && !(zsd>0.0)) {
44 if(verbose>1) printf(" no filtering done.\n");
45 return(0);
46 }
47 if(!(tolerance>0.0)) tolerance=1.0E-03;
48 if(tolerance>=1.0) return(5);
49
50 double sqrt_pi2=sqrt(M_PI)/M_SQRT2; // sqrt(pi)/sqrt(2)
51
52 /*
53 * Filtering in x dimension
54 */
55 if(xsd>0.0) {
56 if(verbose>1) {printf(" filtering in x dimension\n"); fflush(stdout);}
57 int r=ceil(M_SQRT2*xsd*inverfc(0.5*tolerance)); // Radius of the filter
58 double gtrunc[r+1]; // Truncated Gaussian filter for FIR convolution
59 double sdsqrt2=M_SQRT2*xsd;
60 gtrunc[0]=1.0;
61// for(int i=1; i<=r; i++) gtrunc[i]=exp(-0.5*((double)i/xsd)*((double)i/xsd));
62 for(int i=1; i<=r; i++)
63 gtrunc[i]=sqrt_pi2*xsd*(erf(((double)i+0.5)/sdsqrt2) - erf(((double)i-0.5)/sdsqrt2));
64 double sum=gtrunc[0]; for(int i=1; i<=r; i++) sum+=2.0*gtrunc[i];
65 for(int i=0; i<=r; i++) gtrunc[i]/=sum;
66 if(verbose>2) printf(" r=%d \n", r);
67 if(verbose>3) {
68 printf(" gtrunc[]= %g", gtrunc[0]); for(int i=1; i<=r; i++) printf(",%g", gtrunc[i]);
69 printf("\n");
70 }
71
72 for(int t=0; t<dimt; t++) { // Image time frames are processed as separate images
73 if(verbose>5) {printf(" frame %d\n", 1+t); fflush(stdout);}
74 for(int z=0; z<dimz; z++) for(int y=0; y<dimy; y++) {
75 double a[dimx];
76 for(int x=0; x<dimx; x++) a[x]=img->m[z][y][x][t];
77 for(int x=0; x<dimx; x++) {
78 double v=gtrunc[0]*a[x];
79 for(int m=1; m<=r; m++)
80 v+=gtrunc[m]*(a[filterHSSBE(dimx, x-m)] + a[filterHSSBE(dimx, x+m)]);
81 img->m[z][y][x][t]=v;
82 }
83 }
84 }
85 }
86
87
88 /*
89 * Filtering in y dimension
90 */
91 if(ysd>0.0) {
92 if(verbose>1) {printf(" filtering in y dimension\n"); fflush(stdout);}
93 int r=ceil(M_SQRT2*ysd*inverfc(0.5*tolerance)); // Radius of the filter
94 double gtrunc[r+1]; // Truncated Gaussian filter for FIR convolution
95 double sdsqrt2=M_SQRT2*ysd;
96 gtrunc[0]=1.0;
97// for(int i=1; i<=r; i++) gtrunc[i]=exp(-0.5*((double)i/ysd)*((double)i/ysd));
98 for(int i=1; i<=r; i++)
99 gtrunc[i]=sqrt_pi2*ysd*(erf(((double)i+0.5)/sdsqrt2) - erf(((double)i-0.5)/sdsqrt2));
100 double sum=gtrunc[0]; for(int i=1; i<=r; i++) sum+=2.0*gtrunc[i];
101 for(int i=0; i<=r; i++) gtrunc[i]/=sum;
102 if(verbose>2) printf(" r=%d \n", r);
103 if(verbose>3) {
104 printf(" gtrunc[]= %g", gtrunc[0]); for(int i=1; i<=r; i++) printf(",%g", gtrunc[i]);
105 printf("\n");
106 }
107
108 for(int t=0; t<dimt; t++) { // Image time frames are processed as separate images
109 if(verbose>5) {printf(" frame %d\n", 1+t); fflush(stdout);}
110 for(int z=0; z<dimz; z++) for(int x=0; x<dimx; x++) {
111 double a[dimy];
112 for(int y=0; y<dimy; y++) a[y]=img->m[z][y][x][t];
113 for(int y=0; y<dimy; y++) {
114 double v=gtrunc[0]*a[y];
115 for(int m=1; m<=r; m++)
116 v+=gtrunc[m]*(a[filterHSSBE(dimy, y-m)] + a[filterHSSBE(dimy, y+m)]);
117 img->m[z][y][x][t]=v;
118 }
119 }
120 }
121 }
122
123
124 /*
125 * Filtering in z dimension
126 */
127 if(zsd>0.0) {
128 if(verbose>1) {printf(" filtering in z dimension\n"); fflush(stdout);}
129 int r=ceil(M_SQRT2*zsd*inverfc(0.5*tolerance)); // Radius of the filter
130 double gtrunc[r+1]; // Truncated Gaussian filter for FIR convolution
131 double sdsqrt2=M_SQRT2*zsd;
132 gtrunc[0]=1.0;
133// for(int i=1; i<=r; i++) gtrunc[i]=exp(-0.5*((double)i/zsd)*((double)i/zsd));
134 for(int i=1; i<=r; i++)
135 gtrunc[i]=sqrt_pi2*zsd*(erf(((double)i+0.5)/sdsqrt2) - erf(((double)i-0.5)/sdsqrt2));
136 double sum=gtrunc[0]; for(int i=1; i<=r; i++) sum+=2.0*gtrunc[i];
137 for(int i=0; i<=r; i++) gtrunc[i]/=sum;
138 if(verbose>2) printf(" r=%d \n", r);
139 if(verbose>3) {
140 printf(" gtrunc[]= %g", gtrunc[0]); for(int i=1; i<=r; i++) printf(",%g", gtrunc[i]);
141 printf("\n");
142 }
143
144 for(int t=0; t<dimt; t++) { // Image time frames are processed as separate images
145 if(verbose>5) {printf(" frame %d\n", 1+t); fflush(stdout);}
146 for(int y=0; y<dimy; y++) for(int x=0; x<dimx; x++) {
147 double a[dimz];
148 for(int z=0; z<dimz; z++) a[z]=img->m[z][y][x][t];
149 for(int z=0; z<dimz; z++) {
150 double v=gtrunc[0]*a[z];
151 for(int m=1; m<=r; m++)
152 v+=gtrunc[m]*(a[filterHSSBE(dimz, z-m)] + a[filterHSSBE(dimz, z+m)]);
153 img->m[z][y][x][t]=v;
154 }
155 }
156 }
157 }
158 return(0);
159}
160/*****************************************************************************/
161
162/*****************************************************************************/
173int imgGaussianAMFilter(
176 IMG *img,
178 float xsd,
180 float ysd,
182 float zsd,
184 int passNr,
186 double tolerance,
188 int verbose
189) {
190 if(verbose>0) printf("%s(*img, %g, %g, %g, %d, %d)\n", __func__, xsd, ysd, zsd, passNr, verbose);
191
192 if(img==NULL || img->dimz<1 || img->dimy<1 || img->dimx<1 || img->dimt<1) return 1;
193 int dimx=img->dimx, dimy=img->dimy, dimz=img->dimz, dimt=img->dimt;
194 if(dimx==1) xsd=0.0;
195 if(dimy==1) ysd=0.0;
196 if(dimz==1) zsd=0.0;
197 if(xsd>0.0 && (dimx<4 || !(img->sizex>0.0))) return(2);
198 if(ysd>0.0 && (dimy<4 || !(img->sizey>0.0))) return(3);
199 if(zsd>0.0 && (dimz<4 || !(img->sizez>0.0))) return(4);
200 if(!(xsd>0.0) && !(ysd>0.0) && !(zsd>0.0)) {
201 if(verbose>1) printf(" no filtering done.\n");
202 fflush(stdout); return(0);
203 }
204 if(passNr<=0) passNr=4;
205 if(!(tolerance>0.0)) tolerance=1.0E-06;
206
207
208 /*
209 * Filtering in x dimension
210 */
211 if(xsd>0.0) {
212 if(verbose>1) {printf(" filtering in x dimension\n"); fflush(stdout);}
213 double q=xsd;
214 if(0) { // Adjust SD for improved accuracy at low step numbers
215 double f=0.7818+(double)passNr;
216 q*= 1.0 + (0.3165*(double)passNr + 0.5695)/(f*f);
217 }
218 double lambda=(q*q)/(2.0*(double)passNr);
219 double nu=(1.0 + 2.0*lambda - sqrt(1.0 + 4.0*lambda))/(2.0*lambda);
220 double scale=pow(nu/lambda, passNr);
221 int termNr=ceil(log((1.0-nu)*tolerance)/log(nu));
222 if(verbose>2) printf(" nu=%g scale=%g termNr=%d\n", nu, scale, termNr);
223
224 for(int t=0; t<dimt; t++) { // Image time frames are processed as separate images
225 if(verbose>5) {printf(" frame %d\n", 1+t); fflush(stdout);}
226 double sum1=0.0, sum2=0.0;
227 for(int z=0; z<dimz; z++) for(int y=0; y<dimy; y++) {
228 double d[dimx];
229 for(int x=0; x<dimx; x++) sum1+=d[x]=img->m[z][y][x][t];
230 filterAMGaussian(d, dimx, passNr, nu, scale, termNr);
231 for(int x=0; x<dimx; x++) sum2+=img->m[z][y][x][t]=d[x];
232 }
233 double f=sum1/sum2;
234 if(verbose>1 && fabs(1.0-f)>1.0E-10) {
235 printf(" frame %d: sum of voxels differs pre/post filtering %g/%g=%g\n", 1+t, sum1, sum2, f);
236 }
237 }
238 }
239
240 /*
241 * Filtering in y dimension
242 */
243 if(ysd>0.0) {
244 if(verbose>1) {printf(" filtering in y dimension\n"); fflush(stdout);}
245 double q=ysd;
246 if(0) { // Adjust SD for improved accuracy at low step numbers
247 double f=0.7818+(double)passNr;
248 q*= 1.0 + (0.3165*(double)passNr + 0.5695)/(f*f);
249 }
250 double lambda=(q*q)/(2.0*(double)passNr);
251 double nu=(1.0 + 2.0*lambda - sqrt(1.0 + 4.0*lambda))/(2.0*lambda);
252 double scale=pow(nu/lambda, passNr);
253 int termNr=ceil(log((1.0-nu)*tolerance)/log(nu));
254 if(verbose>2) printf(" nu=%g scale=%g termNr=%d\n", nu, scale, termNr);
255
256 for(int t=0; t<dimt; t++) { // Image time frames are processed as separate images
257 if(verbose>5) {printf(" frame %d\n", 1+t); fflush(stdout);}
258 double sum1=0.0, sum2=0.0;
259 for(int z=0; z<dimz; z++) for(int x=0; x<dimx; x++) {
260 double d[dimy];
261 for(int y=0; y<dimy; y++) sum1+=d[y]=img->m[z][y][x][t];
262 filterAMGaussian(d, dimy, passNr, nu, scale, termNr);
263 for(int y=0; y<dimy; y++) sum2+=img->m[z][y][x][t]=d[y];
264 }
265 double f=sum1/sum2;
266 if(verbose>1 && fabs(1.0-f)>1.0E-10) {
267 printf(" frame %d: sum of voxels differs pre/post filtering %g/%g=%g\n", 1+t, sum1, sum2, f);
268 }
269 }
270 }
271
272 /*
273 * Filtering in z dimension
274 */
275 if(zsd>0.0) {
276 if(verbose>1) {printf(" filtering in z dimension\n"); fflush(stdout);}
277 double q=xsd;
278 if(0) { // Adjust SD for improved accuracy at low step numbers
279 double f=0.7818+(double)passNr;
280 q*= 1.0 + (0.3165*(double)passNr + 0.5695)/(f*f);
281 }
282 double lambda=(q*q)/(2.0*(double)passNr);
283 double nu=(1.0 + 2.0*lambda - sqrt(1.0 + 4.0*lambda))/(2.0*lambda);
284 double scale=pow(nu/lambda, passNr);
285 int termNr=ceil(log((1.0-nu)*tolerance)/log(nu));
286 if(verbose>2) printf(" nu=%g scale=%g termNr=%d\n", nu, scale, termNr);
287
288 for(int t=0; t<dimt; t++) { // Image time frames are processed as separate images
289 if(verbose>5) {printf(" frame %d\n", 1+t); fflush(stdout);}
290 double sum1=0.0, sum2=0.0;
291 for(int y=0; y<dimy; y++) for(int x=0; x<dimy; x++) {
292 double d[dimz];
293 for(int z=0; z<dimz; z++) sum1+=d[z]=img->m[z][y][x][t];
294 filterAMGaussian(d, dimz, passNr, nu, scale, termNr);
295 for(int z=0; z<dimz; z++) sum2+=img->m[z][y][x][t]=d[z];
296 }
297 double f=sum1/sum2;
298 if(verbose>1 && fabs(1.0-f)>1.0E-10) {
299 printf(" frame %d: sum of voxels differs pre/post filtering %g/%g=%g\n", 1+t, sum1, sum2, f);
300 }
301 }
302 }
303
304
305 return(0);
306}
307/*****************************************************************************/
308
309/*****************************************************************************/
314int filterAMGaussian(
316 double *a,
318 int n,
320 int passNr,
322 double nu,
324 double scale,
326 int termNr
327) {
328 if(a==NULL || n<3) return(1);
329 if(!(nu>0.0)) return(2);
330 if(passNr<1) passNr=4;
331 if(termNr<6) termNr=6;
332
333 /* Scale the data. */
334 for(int i=0; i<n; i++) a[i]*=scale;
335
336 /* Loop through the passes of filtering */
337 for(int pass=0; pass<passNr; pass++) {
338 a[0]=filterAMLB(a, n, nu, termNr);
339 for(int i=1; i<n; i++) a[i]+=nu*a[i-1]; // causal filter
340 a[n-1]/=(1.0-nu);
341 for(int i=n-1; i>0; i--) a[i-1]+=nu*a[i]; // anti-causal filter
342 }
343 return(0);
344}
345/*****************************************************************************/
346
347/*****************************************************************************/
352inline int filterHSSBE(
354 int n,
356 int r
357) {
358 while(1) {
359 if(r<0) r=-1-r;
360 else if(r>=n) r=2*n-1-r;
361 else break;
362 }
363 return(r);
364}
365/*****************************************************************************/
366
367/*****************************************************************************/
372double filterAMLB(
374 double *a,
376 int n,
378 double nu,
380 int termNr
381) {
382 if(a==NULL || n<1) return(nan(""));
383 double h=1.0, sum=a[0];
384 for(int m=1; m<termNr; m++) {
385 h*=nu;
386 sum+=h*a[filterHSSBE(n, -m)];
387 }
388 return(sum);
389}
390/*****************************************************************************/
391
392/*****************************************************************************/
408int imgGaussianEBoxFilter(
411 IMG *img,
413 float xsd,
415 float ysd,
417 float zsd,
419 int passNr,
421 int verbose
422) {
423 if(verbose>0) printf("%s(*img, %g, %g, %g, %d, %d)\n", __func__, xsd, ysd, zsd, passNr, verbose);
424
425 if(img==NULL || img->dimz<1 || img->dimy<1 || img->dimx<1 || img->dimt<1) return 1;
426 int dimx=img->dimx, dimy=img->dimy, dimz=img->dimz, dimt=img->dimt;
427 if(dimx==1) xsd=0.0;
428 if(dimy==1) ysd=0.0;
429 if(dimz==1) zsd=0.0;
430 if(xsd>0.0 && (dimx<4 || !(img->sizex>0.0))) return(2);
431 if(ysd>0.0 && (dimy<4 || !(img->sizey>0.0))) return(3);
432 if(zsd>0.0 && (dimz<4 || !(img->sizez>0.0))) return(4);
433 if(!(xsd>0.0) && !(ysd>0.0) && !(zsd>0.0)) {
434 if(verbose>1) printf(" no filtering done.\n");
435 return(0);
436 }
437 if(passNr<=0) passNr=4;
438
439
440 /*
441 * Filtering in x dimension
442 */
443 if(xsd>0.0) {
444 if(verbose>1) {printf(" filtering in x dimension\n"); fflush(stdout);}
445 int rib;
446 double wob, wib, s2p=xsd*xsd/(double)passNr;
447 /* Calculate inner box radius (by Wells) */
448 rib=0.5*sqrt(1.0 + 12.0*s2p) - 0.5;
449 if(verbose>3) printf(" precise_r := %f\n", 0.5*sqrt(1.0+12.0*s2p));
450 /* Alpha */
451 double a=(2*rib+1)*(rib*(rib+1) - 3.0*s2p) / (6.0*(s2p - (rib+1)*(rib+1)));
452 /* Outer box weight */
453 wob = a/(2.0*(a+rib)+1.0);
454 /* Inner box weight */
455 wib = 1.0-wob;
456 if(verbose>2) {printf(" rib=%d wob=%g wib=%g a=%g\n", rib, wob, wib, a); fflush(stdout);}
457
458 for(int t=0; t<dimt; t++) { // Image time frames are processed as separate images
459 if(verbose>5) {printf(" frame %d\n", 1+t); fflush(stdout);}
460 double sum1=0.0, sum2=0.0;
461 for(int z=0; z<dimz; z++) for(int y=0; y<dimy; y++) {
462 double d[dimx];
463 for(int x=0; x<dimx; x++) sum1+=d[x]=img->m[z][y][x][t];
464 filterEbox(d, dimx, wob, wib, rib, passNr);
465 for(int x=0; x<dimx; x++) sum2+=img->m[z][y][x][t]=d[x];
466 }
467 double f=sum1/sum2;
468 if(fabs(1.0-f)>1.0E-10) {
469 if(verbose>1)
470 printf(" frame %d: sum of voxels differs pre/post filtering %g/%g=%g\n", 1+t, sum1, sum2, f);
471 for(int z=0; z<dimz; z++) for(int y=0; y<dimy; y++) for(int x=0; x<dimx; x++)
472 img->m[z][y][x][t]*=f;
473 }
474 }
475 }
476
477 /*
478 * Filtering in y dimension
479 */
480 if(ysd>0.0) {
481 if(verbose>1) {printf(" filtering in y dimension\n"); fflush(stdout);}
482 int rib;
483 double wob, wib, s2p=ysd*ysd/(double)passNr;
484 /* Calculate inner box radius */
485 rib=0.5*sqrt(1.0 + 12.0*s2p) - 0.5;
486 /* Alpha */
487 double a=(2*rib+1)*(rib*(rib+1) - 3.0*s2p) / (6.0*(s2p - (rib+1)*(rib+1)));
488 /* Outer box weight */
489 wob = a/(2.0*(a+rib)+1);
490 /* Inner box weight */
491 wib = 1.0-wob;
492 if(verbose>2) {printf(" rib=%d wob=%g wib=%g a=%g\n", rib, wob, wib, a); fflush(stdout);}
493
494 for(int t=0; t<dimt; t++) { // Image time frames are processed as separate images
495 if(verbose>5) {printf(" frame %d\n", 1+t); fflush(stdout);}
496 double sum1=0.0, sum2=0.0;
497 for(int z=0; z<dimz; z++) for(int x=0; x<dimx; x++) {
498 double d[dimy];
499 for(int y=0; y<dimy; y++) sum1+=d[y]=img->m[z][y][x][t];
500 filterEbox(d, dimy, wob, wib, rib, passNr);
501 for(int y=0; y<dimy; y++) sum2+=img->m[z][y][x][t]=d[y];
502 }
503 double f=sum1/sum2;
504 if(fabs(1.0-f)>1.0E-10) {
505 if(verbose>1)
506 printf(" frame %d: sum of voxels differs pre/post filtering %g/%g=%g\n", 1+t, sum1, sum2, f);
507 for(int z=0; z<dimz; z++) for(int y=0; y<dimy; y++) for(int x=0; x<dimx; x++)
508 img->m[z][y][x][t]*=f;
509 }
510 }
511 }
512
513 /*
514 * Filtering in z dimension
515 */
516 if(zsd>0.0) {
517 if(verbose>1) {printf(" filtering in z dimension\n"); fflush(stdout);}
518 int rib;
519 double wob, wib, s2p=zsd*zsd/(double)passNr;
520 /* Calculate inner box radius */
521 rib=0.5*sqrt(1.0 + 12.0*s2p) - 0.5;
522 /* Alpha */
523 double a=(2*rib+1)*(rib*(rib+1) - 3.0*s2p) / (6.0*(s2p - (rib+1)*(rib+1)));
524 /* Outer box weight */
525 wob = a/(2.0*(a+rib)+1);
526 /* Inner box weight */
527 wib = 1.0-wob;
528 if(verbose>2) {printf(" rib=%d wob=%g wib=%g a=%g\n", rib, wob, wib, a); fflush(stdout);}
529
530 for(int t=0; t<dimt; t++) { // Image time frames are processed as separate images
531 if(verbose>5) {printf(" frame %d\n", 1+t); fflush(stdout);}
532 double sum1=0.0, sum2=0.0;
533 for(int y=0; y<dimy; y++) for(int x=0; x<dimx; x++) {
534 double d[dimz];
535 for(int z=0; z<dimz; z++) sum1+=d[z]=img->m[z][y][x][t];
536 filterEbox(d, dimz, wob, wib, rib, passNr);
537 for(int z=0; z<dimz; z++) sum2+=img->m[z][y][x][t]=d[z];
538 }
539 double f=sum1/sum2;
540 if(fabs(1.0-f)>1.0E-10) {
541 if(verbose>1)
542 printf(" frame %d: sum of voxels differs pre/post filtering %g/%g=%g\n", 1+t, sum1, sum2, f);
543 for(int z=0; z<dimz; z++) for(int y=0; y<dimy; y++) for(int x=0; x<dimx; x++)
544 img->m[z][y][x][t]*=f;
545 }
546 }
547 }
548
549 return(0);
550}
551/*****************************************************************************/
552
553/*****************************************************************************/
558int filterEbox(
560 double *a,
562 int n,
564 double wob,
566 double wib,
568 int rib,
570 int passNr
571) {
572 if(a==NULL || n<3 || wob<=0.0 || wib<=0.0 || rib<0) return(1);
573 double buf[n];
574 for(int pass=0; pass<passNr; pass++) {
575 double accum=0.0;
576 for(int i=-rib; i<=rib; i++) accum+=a[filterHSSBE(n, i)];
577 buf[0]=wob*(a[filterHSSBE(n, rib+1)] + a[filterHSSBE(n, -rib-1)]) + (wob+wib)*accum;
578 accum=buf[0];
579 for(int i=1; i<n; i++) {
580 accum += wob*(a[filterHSSBE(n, i+rib+1)] - a[filterHSSBE(n, i-rib-2)])
581 + wib*(a[filterHSSBE(n, i+rib)] - a[filterHSSBE(n, i-rib-1)]);
582 buf[i]=accum;
583 }
584 for(int i=0; i<n; i++) a[i]=buf[i];
585 }
586 return(0);
587}
588/*****************************************************************************/
589
590/*****************************************************************************/
595int fMean1DFilter(
597 float *data,
599 const int n,
601 const int s
602) {
603 if(data==NULL || n<1 || s<1) return(1);
604 float buf[n];
605 for(int i=0; i<n; i++) buf[i]=data[i];
606 for(int i=0; i<n; i++) {
607 data[i]=0.0;
608 int a=i-(s-1)/2; if(a<0) a=0;
609 int b=i+(s-1)/2; if(b>=n) b=n-1;
610 int m=1+b-a;
611 for(int j=a; j<=b; j++) data[i]+=buf[j];
612 data[i]/=(float)m;
613 }
614 return(0);
615}
616/*****************************************************************************/
617
618/*****************************************************************************/
624int imgMeanFilter(
627 IMG *img,
629 int xn,
631 int yn,
633 int zn,
635 int tn,
637 int verbose
638) {
639 if(verbose>0) printf("%s(*img, %d, %d, %d, %d)\n", __func__, xn, yn, zn, verbose);
640
641 if(img==NULL || img->dimz<1 || img->dimy<1 || img->dimx<1 || img->dimt<1) return(1);
642 int dimx=img->dimx, dimy=img->dimy, dimz=img->dimz, dimt=img->dimt;
643 if(dimx==1) xn=0;
644 if(dimy==1) yn=0;
645 if(dimz==1) zn=0;
646 if(dimt==1) tn=0;
647 if(!(xn>0) && !(yn>0) && !(zn>0) && !(tn>0)) {
648 if(verbose>1) printf(" no filtering done.\n");
649 return(0);
650 }
651 if((xn>0 && xn<3) || (yn>0 && yn<3) || (zn>0 && zn<3) || (tn>0 && tn<3)) {
652 if(verbose>1) printf(" invalid pixel number.\n");
653 return(1);
654 }
655 if((xn>0 && xn%2==0) || (yn>0 && yn%2==0) || (zn>0 && zn%2==0) || (tn>0 && tn%2==0)) {
656 if(verbose>1) printf(" invalid pixel number.\n");
657 return(1);
658 }
659 if((xn>0 && dimx<=xn) || (yn>0 && dimy<=yn) || (zn>0 && dimz<=zn) || (tn>0 && dimt<=tn)) {
660 if(verbose>1) printf(" invalid pixel number.\n");
661 return(1);
662 }
663
664
665 /*
666 * Filtering in x dimension
667 */
668 if(xn>=3) {
669 if(verbose>1) {printf(" filtering in x dimension\n"); fflush(stdout);}
670 float buf[dimx];
671 for(int z=0; z<dimz; z++)
672 for(int y=0; y<dimy; y++)
673 for(int t=0; t<dimt; t++) {
674 for(int x=0; x<dimx; x++) buf[x]=img->m[z][y][x][t];
675 if(!fMean1DFilter(buf, dimx, xn))
676 for(int x=0; x<dimx; x++) img->m[z][y][x][t]=buf[x];
677 }
678 }
679
680 /*
681 * Filtering in y dimension
682 */
683 if(yn>=3) {
684 if(verbose>1) {printf(" filtering in y dimension\n"); fflush(stdout);}
685 float buf[dimy];
686 for(int z=0; z<dimz; z++)
687 for(int x=0; x<dimx; x++)
688 for(int t=0; t<dimt; t++) {
689 for(int y=0; y<dimy; y++) buf[y]=img->m[z][y][x][t];
690 if(!fMean1DFilter(buf, dimy, yn))
691 for(int y=0; y<dimy; y++) img->m[z][y][x][t]=buf[y];
692 }
693 }
694
695 /*
696 * Filtering in z dimension
697 */
698 if(zn>=3) {
699 if(verbose>1) {printf(" filtering in z dimension\n"); fflush(stdout);}
700 float buf[dimz];
701 for(int y=0; y<dimy; y++)
702 for(int x=0; x<dimx; x++)
703 for(int t=0; t<dimt; t++) {
704 for(int z=0; z<dimz; z++) buf[z]=img->m[z][y][x][t];
705 if(!fMean1DFilter(buf, dimz, zn))
706 for(int z=0; z<dimz; z++) img->m[z][y][x][t]=buf[z];
707 }
708 }
709
710 /*
711 * Filtering in t dimension
712 */
713 if(tn>=3) {
714 if(verbose>1) {printf(" filtering in t dimension\n"); fflush(stdout);}
715 float buf[dimt];
716 for(int z=0; z<dimz; z++)
717 for(int y=0; y<dimy; y++)
718 for(int x=0; x<dimx; x++) {
719 for(int t=0; t<dimt; t++) buf[t]=img->m[z][y][x][t];
720 if(!fMean1DFilter(buf, dimt, tn))
721 for(int t=0; t<dimt; t++) img->m[z][y][x][t]=buf[t];
722 }
723 }
724
725
726 return(0);
727}
728/*****************************************************************************/
729
730/*****************************************************************************/
inverfc
double inverfc(double x)
Definition doubleutil.c:205
filterHSSBE
int filterHSSBE(int n, int r)
Definition imgfilter.c:352
fMean1DFilter
int fMean1DFilter(float *data, const int n, const int s)
Definition imgfilter.c:595
imgGaussianEBoxFilter
int imgGaussianEBoxFilter(IMG *img, float xsd, float ysd, float zsd, int passNr, int verbose)
Definition imgfilter.c:408
imgMeanFilter
int imgMeanFilter(IMG *img, int xn, int yn, int zn, int tn, int verbose)
Definition imgfilter.c:624
imgGaussianFIRFilter
int imgGaussianFIRFilter(IMG *img, float xsd, float ysd, float zsd, double tolerance, int verbose)
Definition imgfilter.c:18
filterAMGaussian
int filterAMGaussian(double *a, int n, int passNr, double nu, double scale, int termNr)
Definition imgfilter.c:314
filterAMLB
double filterAMLB(double *a, int n, double nu, int termNr)
Definition imgfilter.c:372
filterEbox
int filterEbox(double *a, int n, double wob, double wib, int rib, int passNr)
Definition imgfilter.c:558
imgGaussianAMFilter
int imgGaussianAMFilter(IMG *img, float xsd, float ysd, float zsd, int passNr, double tolerance, int verbose)
Definition imgfilter.c:173
libtpcimgp.h
Header file for libtpcimgp.
IMG::sizex
float sizex
Definition libtpcimgio.h:1714
IMG::dimx
unsigned short int dimx
Definition libtpcimgio.h:1780
IMG::m
float **** m
Definition libtpcimgio.h:1799
IMG::dimt
unsigned short int dimt
Definition libtpcimgio.h:1778
IMG::sizey
float sizey
Definition libtpcimgio.h:1716
IMG::dimz
unsigned short int dimz
Definition libtpcimgio.h:1784
IMG::dimy
unsigned short int dimy
Definition libtpcimgio.h:1782
IMG::sizez
float sizez
Definition libtpcimgio.h:1718
Generated on Fri Sep 19 2025 13:32:05 for TPCCLIB by doxygen 1.12.0