TPCCLIB
Loading...
Searching...
No Matches
Functions
imgfilter.c File Reference

Gaussian IMG filters. More...

#include "libtpcimgp.h"

Go to the source code of this file.

Functions

int imgGaussianFIRFilter (IMG *img, float xsd, float ysd, float zsd, double tolerance, int verbose)
int imgGaussianAMFilter (IMG *img, float xsd, float ysd, float zsd, int passNr, double tolerance, int verbose)
int filterAMGaussian (double *a, int n, int passNr, double nu, double scale, int termNr)
int filterHSSBE (int n, int r)
double filterAMLB (double *a, int n, double nu, int termNr)
int imgGaussianEBoxFilter (IMG *img, float xsd, float ysd, float zsd, int passNr, int verbose)
int filterEbox (double *a, int n, double wob, double wib, int rib, int passNr)
int fMean1DFilter (float *data, const int n, const int s)
int imgMeanFilter (IMG *img, int xn, int yn, int zn, int tn, int verbose)

Detailed Description

Gaussian IMG filters.

Definition in file imgfilter.c.

Function Documentation

◆ filterAMGaussian()

int filterAMGaussian ( double * a,
int n,
int passNr,
double nu,
double scale,
int termNr )

1D Alvarez-Mazorra Gaussian filtering.

See also
imgGaussianAMFilter
Returns
If an error is encountered, function returns a non-zero value. Otherwise 0 is returned.
Parameters
aPointer to source data array, overwritten by the results.
nSample number.
passNrNumber of filtering passes; usually 4.
nuFilter coefficient.
scaleFilter scale.
termNrNumber of terms needed to approximate the sum with required accuracy for the left boundary.

Definition at line 314 of file imgfilter.c.

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}
filterAMLB
double filterAMLB(double *a, int n, double nu, int termNr)
Definition imgfilter.c:372

Referenced by imgGaussianAMFilter().

◆ filterAMLB()

double filterAMLB ( double * a,
int n,
double nu,
int termNr )

Left boundary for 1D Alvarez-Mazorra Gaussian filtering.

See also
filterAMGaussian, filterHSSBE
Returns
The sum approximating the first filtered sample value, or NaN in case of an error.
Parameters
aPointer to source data array; not modified.
nSample number.
nuFiltering parameter.
termNrNumber of terms.

Definition at line 372 of file imgfilter.c.

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}
filterHSSBE
int filterHSSBE(int n, int r)
Definition imgfilter.c:352

Referenced by filterAMGaussian().

◆ filterEbox()

int filterEbox ( double * a,
int n,
double wob,
double wib,
int rib,
int passNr )

1D extended box filtering .

See also
imgFast1DGaussianFilter
Returns
If an error is encountered, function returns a non-zero value. Otherwise 0 is returned.
Parameters
aPointer to source data array.
nSample number.
wobWeight of the outer box.
wibWeight of the inner box.
ribRadius of the inner box.
passNrNumber of filtering passes.

Definition at line 558 of file imgfilter.c.

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}

Referenced by imgGaussianEBoxFilter().

◆ filterHSSBE()

int filterHSSBE ( int n,
int r )
inline

Half-sample symmetric boundary extension.

See also
filterEbox, filterAMLB
Returns
Reflected sample in array[0..n-1].
Parameters
nSample number.
rRequested sample, possibly outside array[0..N-1].

Definition at line 352 of file imgfilter.c.

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}

Referenced by filterAMLB(), filterEbox(), and imgGaussianFIRFilter().

◆ fMean1DFilter()

int fMean1DFilter ( float * data,
const int n,
const int s )

In-place mean 1D-filtering.

See also
imgMeanFilter
Returns
If an error is encountered, function returns a non-zero value. Otherwise 0 is returned.
Parameters
dataPointer to array to be filtered.
nArray size.
sFilter size; should be odd number 3,5,7,... but not verified here.

Definition at line 595 of file imgfilter.c.

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}

Referenced by imgMeanFilter().

◆ imgGaussianAMFilter()

int imgGaussianAMFilter ( IMG * img,
float xsd,
float ysd,
float zsd,
int passNr,
double tolerance,
int verbose )

Alvarez-Mazorra approximate Gaussian image filtering in 1-3 dimensions.

References:

See also
imgGaussianFIRFilter, imgSS
Returns
If an error is encountered, function returns a non-zero value. Otherwise 0 is returned.
Parameters
imgImage data to be processed; data is overwritten with filtered image. All image frames are processed.
xsdGaussian S.D. in pixels in x dimension; SD=FWHM/2.355. Enter 0 to not filter in x dimension.
ysdGaussian S.D. in pixels in y dimension; SD=FWHM/2.355. Enter 0 to not filter in y dimension.
zsdGaussian S.D. in pixels in z dimension; SD=FWHM/2.355. Enter 0 to not filter in z dimension.
passNrNumber of filtering passes, usually 4 or more; enter 0 to use the default steps (4).
toleranceTolerance; enter 0.0 to use the default (1.0E-06).
verboseVerbose level; if zero, then only warnings are printed into stderr.

Definition at line 173 of file imgfilter.c.

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}
filterAMGaussian
int filterAMGaussian(double *a, int n, int passNr, double nu, double scale, int termNr)
Definition imgfilter.c:314
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

◆ imgGaussianEBoxFilter()

int imgGaussianEBoxFilter ( IMG * img,
float xsd,
float ysd,
float zsd,
int passNr,
int verbose )

Gaussian extended box image filtering in 1-3 dimensions.

References:

See also
imgFast2DGaussianFilter, imgGaussianFilter, imgSS
Returns
If an error is encountered, function returns a non-zero value. Otherwise 0 is returned.
Parameters
imgImage data to be processed; data is overwritten with filtered image. All image frames are processed.
xsdGaussian S.D. in pixels in x dimension; SD=FWHM/2.355. Enter 0 to not filter in x dimension.
ysdGaussian S.D. in pixels in y dimension; SD=FWHM/2.355. Enter 0 to not filter in y dimension.
zsdGaussian S.D. in pixels in z dimension; SD=FWHM/2.355. Enter 0 to not filter in z dimension.
passNrNumber of filtering passes, usually 3-5; enter 0 to use the default steps (4).
verboseVerbose level; if zero, then only warnings are printed into stderr.

Definition at line 408 of file imgfilter.c.

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}
filterEbox
int filterEbox(double *a, int n, double wob, double wib, int rib, int passNr)
Definition imgfilter.c:558

◆ imgGaussianFIRFilter()

int imgGaussianFIRFilter ( IMG * img,
float xsd,
float ysd,
float zsd,
double tolerance,
int verbose )

Gaussian finite impulse response (FIR) isotropic or anisotropic image filtering in 1-3 dimensions.

Based on the code in: Getreuer P. A survey of Gaussian convolution algorithms. Image Processing On Line 2013;3:286–310. https://doi.org/10.5201/ipol.2013.87

See also
imgGaussianAMFilter, imgSS
Returns
If an error is encountered, function returns a non-zero value. Otherwise 0 is returned.
Parameters
imgImage data to be processed; data is overwritten with filtered image. All image frames are processed.
xsdGaussian S.D. in pixels in x dimension; SD=FWHM/2.355. Enter 0 to not filter in x dimension.
ysdGaussian S.D. in pixels in y dimension; SD=FWHM/2.355. Enter 0 to not filter in y dimension.
zsdGaussian S.D. in pixels in z dimension; SD=FWHM/2.355. Enter 0 to not filter in z dimension.
toleranceTolerance; enter 0.0 to use the default (1.0E-03).
verboseVerbose level; if zero, then only warnings are printed into stderr.

Definition at line 18 of file imgfilter.c.

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}
inverfc
double inverfc(double x)
Definition doubleutil.c:205

Referenced by imgPVCRRL(), and imgPVCRVC().

◆ imgMeanFilter()

int imgMeanFilter ( IMG * img,
int xn,
int yn,
int zn,
int tn,
int verbose )

Image mean filtering in 1-4 dimensions.

See also
imgSmoothOverFrames, imgGaussianAMFilter, imgSS, fMean1DFilter
Returns
If an error is encountered, function returns a non-zero value. Otherwise 0 is returned.
Parameters
imgImage data to be processed; data is overwritten with filtered image. All image frames are processed.
xnNumber of adjacent pixels in x dimension; Enter 0 to not filter in x dimension.
ynNumber of adjacent pixels in y dimension; Enter 0 to not filter in y dimension.
znNumber of adjacent pixels in z dimension; Enter 0 to not filter in z dimension.
tnNumber of adjacent pixels in t dimension; Enter 0 to not filter in t dimension.
verboseVerbose level; if zero, then only warnings are printed into stderr.

Definition at line 624 of file imgfilter.c.

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}
fMean1DFilter
int fMean1DFilter(float *data, const int n, const int s)
Definition imgfilter.c:595
Generated on for TPCCLIB by doxygen 1.15.0