diff options
Diffstat (limited to 'src/silx/math/fit/filters')
-rw-r--r-- | src/silx/math/fit/filters/include/filters.h | 45 | ||||
-rw-r--r-- | src/silx/math/fit/filters/src/smoothnd.c | 317 | ||||
-rw-r--r-- | src/silx/math/fit/filters/src/snip1d.c | 149 | ||||
-rw-r--r-- | src/silx/math/fit/filters/src/snip2d.c | 96 | ||||
-rw-r--r-- | src/silx/math/fit/filters/src/snip3d.c | 186 | ||||
-rw-r--r-- | src/silx/math/fit/filters/src/strip.c | 118 |
6 files changed, 911 insertions, 0 deletions
diff --git a/src/silx/math/fit/filters/include/filters.h b/src/silx/math/fit/filters/include/filters.h new file mode 100644 index 0000000..1ee9a95 --- /dev/null +++ b/src/silx/math/fit/filters/include/filters.h @@ -0,0 +1,45 @@ +/*########################################################################## +# Copyright (C) 2016 European Synchrotron Radiation Facility +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in +# all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +# THE SOFTWARE. +# +# ############################################################################*/ + +#ifndef FITFILTERS_H +#define FITFILTERS_H + +/* Background functions */ +void snip1d(double *data, int size, int width); +//void snip1d_multiple(double *data, int n_channels, int snip_width, int n_spectra); +void snip2d(double *data, int nrows, int ncolumns, int width); +void snip3d(double *data, int nx, int ny, int nz, int width); + +int strip(double* input, long len_input, double c, long niter, int deltai, + long* anchors, long len_anchors, double* output); + +/* Smoothing functions */ + +int SavitskyGolay(double* input, long len_input, int npoints, double* output); + +void smooth1d(double *data, int size); +void smooth2d(double *data, int size0, int size1); +void smooth3d(double *data, int size0, int size1, int size2); + + +#endif /* #define FITFILTERS_H */ diff --git a/src/silx/math/fit/filters/src/smoothnd.c b/src/silx/math/fit/filters/src/smoothnd.c new file mode 100644 index 0000000..cb96961 --- /dev/null +++ b/src/silx/math/fit/filters/src/smoothnd.c @@ -0,0 +1,317 @@ +#/*########################################################################## +# +# Copyright (c) 2004-2016 European Synchrotron Radiation Facility +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in +# all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +# THE SOFTWARE. +# +#############################################################################*/ +#include <stdlib.h> +#include <string.h> +#include <math.h> +#define MIN(x, y) (((x) < (y)) ? (x) : (y)) +#define MAX(x, y) (((x) > (y)) ? (x) : (y)) + +#define MAX_SAVITSKY_GOLAY_WIDTH 101 +#define MIN_SAVITSKY_GOLAY_WIDTH 3 + +/* Wrapped functions */ +void smooth1d(double *data, int size); +void smooth2d(double *data, int size0, int size1); +void smooth3d(double *data, int size0, int size1, int size2); +int SavitskyGolay(double* input, long len_input, int npoints, double* output); + +/* Internal functions */ +long index2d(long row_idx, long col_idx, long ncols); +long index3d(long x_idx, long y_idx, long z_idx, long ny, long nz); +void smooth1d_rows(double *data, long nrows, long ncols); +void smooth1d_cols(double *data, long nrows, long ncols); +void smooth1d_x(double *data, long nx, long ny, long nz); +void smooth1d_y(double *data, long nx, long ny, long nz); +void smooth1d_z(double *data, long nx, long ny, long nz); +void smooth2d_yzslice(double *data, long nx, long ny, long nz); +void smooth2d_xzslice(double *data, long nx, long ny, long nz); +void smooth2d_xyslice(double *data, long nx, long ny, long nz); + + +/* Simple smoothing of a 1D array */ +void smooth1d(double *data, int size) +{ + long i; + double prev_sample; + double next_sample; + + if (size < 3) + { + return; + } + prev_sample = data[0]; + for (i=0; i<(size-1); i++) + { + next_sample = 0.25 * (prev_sample + 2 * data[i] + data[i+1]); + prev_sample = data[i]; + data[i] = next_sample; + } + data[size-1] = 0.25 * prev_sample + 0.75 * data[size-1]; + return; +} + +/* Smoothing of a 2D array*/ +void smooth2d(double *data, int nrows, int ncols) +{ + /* smooth the first dimension (rows) */ + smooth1d_rows(data, nrows, ncols); + + /* smooth the 2nd dimension */ + smooth1d_cols(data, nrows, ncols); +} + +/* Smoothing of a 3DÂ array */ +void smooth3d(double *data, int nx, int ny, int nz) +{ + smooth2d_xyslice(data, nx, ny, nz); + smooth2d_xzslice(data, nx, ny, nz); + smooth2d_yzslice(data, nx, ny, nz); +} + +/* 1D Savitsky-Golay smoothing */ +int SavitskyGolay(double* input, long len_input, int npoints, double* output) +{ + + //double dpoints = 5.; + double coeff[MAX_SAVITSKY_GOLAY_WIDTH]; + int i, j, m; + double dhelp, den; + double *data; + + memcpy(output, input, len_input * sizeof(double)); + + if (!(npoints % 2)) npoints +=1; + + if((npoints < MIN_SAVITSKY_GOLAY_WIDTH) || (len_input < npoints) || \ + (npoints > MAX_SAVITSKY_GOLAY_WIDTH)) + { + /* do not smooth data */ + return 1; + } + + /* calculate the coefficients */ + m = (int) (npoints/2); + den = (double) ((2*m-1) * (2*m+1) * (2*m + 3)); + for (i=0; i<= m; i++){ + coeff[m+i] = (double) (3 * (3*m*m + 3*m - 1 - 5*i*i )); + coeff[m-i] = coeff[m+i]; + } + + /* simple smoothing at the beginning */ + for (j=0; j<=(int)(npoints/3); j++) + { + smooth1d(output, m); + } + + /* simple smoothing at the end */ + for (j=0; j<=(int)(npoints/3); j++) + { + smooth1d((output+len_input-m-1), m); + } + + /*one does not need the whole spectrum buffer, but code is clearer */ + data = (double *) malloc(len_input * sizeof(double)); + memcpy(data, output, len_input * sizeof(double)); + + /* the actual SG smoothing in the middle */ + for (i=m; i<(len_input-m); i++){ + dhelp = 0; + for (j=-m;j<=m;j++) { + dhelp += coeff[m+j] * (*(data+i+j)); + } + if(dhelp > 0.0){ + *(output+i) = dhelp / den; + } + } + free(data); + return (0); +} + +/*********************/ +/* Utility functions */ +/*********************/ + +long index2d(long row_idx, long col_idx, long ncols) +{ + return (row_idx*ncols+col_idx); +} + +/* Apply smooth 1d on all rows in a 2D array*/ +void smooth1d_rows(double *data, long nrows, long ncols) +{ + long row_idx; + + for (row_idx=0; row_idx < nrows; row_idx++) + { + smooth1d(&data[row_idx * ncols], ncols); + } +} + +/* Apply smooth 1d on all columns in a 2D array*/ +void smooth1d_cols(double *data, long nrows, long ncols) +{ + long row_idx, col_idx; + long this_idx2d, next_idx2d; + double prev_sample; + double next_sample; + + for (col_idx=0; col_idx < ncols; col_idx++) + { + prev_sample = data[index2d(0, col_idx, ncols)]; + for (row_idx=0; row_idx<(nrows-1); row_idx++) + { + this_idx2d = index2d(row_idx, col_idx, ncols); + next_idx2d = index2d(row_idx+1, col_idx, ncols); + + next_sample = 0.25 * (prev_sample + \ + 2 * data[this_idx2d] + \ + data[next_idx2d]); + prev_sample = data[this_idx2d]; + data[this_idx2d] = next_sample; + } + + this_idx2d = index2d(nrows-1, col_idx, ncols); + data[this_idx2d] = 0.25 * prev_sample + 0.75 * data[this_idx2d]; + } +} + +long index3d(long x_idx, long y_idx, long z_idx, long ny, long nz) +{ + return ((x_idx*ny + y_idx) * nz + z_idx); +} + +/* Apply smooth 1d along first dimension in a 3D array*/ +void smooth1d_x(double *data, long nx, long ny, long nz) +{ + long x_idx, y_idx, z_idx; + long this_idx3d, next_idx3d; + double prev_sample; + double next_sample; + + for (y_idx=0; y_idx < ny; y_idx++) + { + for (z_idx=0; z_idx < nz; z_idx++) + { + prev_sample = data[index3d(0, y_idx, z_idx, ny, nz)]; + for (x_idx=0; x_idx<(nx-1); x_idx++) + { + this_idx3d = index3d(x_idx, y_idx, z_idx, ny, nz); + next_idx3d = index3d(x_idx+1, y_idx, z_idx, ny, nz); + + next_sample = 0.25 * (prev_sample + \ + 2 * data[this_idx3d] + \ + data[next_idx3d]); + prev_sample = data[this_idx3d]; + data[this_idx3d] = next_sample; + } + + this_idx3d = index3d(nx-1, y_idx, z_idx, ny, nz); + data[this_idx3d] = 0.25 * prev_sample + 0.75 * data[this_idx3d]; + } + } +} + +/* Apply smooth 1d along second dimension in a 3D array*/ +void smooth1d_y(double *data, long nx, long ny, long nz) +{ + long x_idx, y_idx, z_idx; + long this_idx3d, next_idx3d; + double prev_sample; + double next_sample; + + for (x_idx=0; x_idx < nx; x_idx++) + { + for (z_idx=0; z_idx < nz; z_idx++) + { + prev_sample = data[index3d(x_idx, 0, z_idx, ny, nz)]; + for (y_idx=0; y_idx<(ny-1); y_idx++) + { + this_idx3d = index3d(x_idx, y_idx, z_idx, ny, nz); + next_idx3d = index3d(x_idx, y_idx+1, z_idx, ny, nz); + + next_sample = 0.25 * (prev_sample + \ + 2 * data[this_idx3d] + \ + data[next_idx3d]); + prev_sample = data[this_idx3d]; + data[this_idx3d] = next_sample; + } + + this_idx3d = index3d(x_idx, ny-1, z_idx, ny, nz); + data[this_idx3d] = 0.25 * prev_sample + 0.75 * data[this_idx3d]; + } + } +} + +/* Apply smooth 1d along third dimension in a 3D array*/ +void smooth1d_z(double *data, long nx, long ny, long nz) +{ + long x_idx, y_idx; + long idx3d_first_sample; + + for (x_idx=0; x_idx < nx; x_idx++) + { + for (y_idx=0; y_idx < ny; y_idx++) + { + idx3d_first_sample = index3d(x_idx, y_idx, 0, ny, nz); + /*We can use regular 1D smoothing function because z samples + are contiguous in memory*/ + smooth1d(&data[idx3d_first_sample], nz); + } + } +} + +/* 2D smoothing of a YZ slice in a 3D volume*/ +void smooth2d_yzslice(double *data, long nx, long ny, long nz) +{ + long x_idx; + long slice_size = ny * nz; + + /* a YZ slice is a "normal" 2D array of memory-contiguous data*/ + for (x_idx=0; x_idx < nx; x_idx++) + { + smooth2d(&data[x_idx*slice_size], ny, nz); + } +} + +/* 2D smoothing of a XZ slice in a 3D volume*/ +void smooth2d_xzslice(double *data, long nx, long ny, long nz) +{ + + /* smooth along the first dimension */ + smooth1d_x(data, nx, ny, nz); + + /* smooth along the third dimension */ + smooth1d_z(data, nx, ny, nz); +} + +/* 2D smoothing of a XY slice in a 3D volume*/ +void smooth2d_xyslice(double *data, long nx, long ny, long nz) +{ + /* smooth along the first dimension */ + smooth1d_x(data, nx, ny, nz); + + /* smooth along the second dimension */ + smooth1d_y(data, nx, ny, nz); +} + diff --git a/src/silx/math/fit/filters/src/snip1d.c b/src/silx/math/fit/filters/src/snip1d.c new file mode 100644 index 0000000..994a272 --- /dev/null +++ b/src/silx/math/fit/filters/src/snip1d.c @@ -0,0 +1,149 @@ +#/*########################################################################## +# Copyright (c) 2004-2016 European Synchrotron Radiation Facility +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in +# all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +# THE SOFTWARE. +# +#############################################################################*/ +/* + Implementation of the algorithm SNIP in 1D described in + Miroslav Morhac et al. Nucl. Instruments and Methods in Physics Research A401 (1997) 113-132. + + The original idea for 1D and the low-statistics-digital-filter (lsdf) come from + C.G. Ryan et al. Nucl. Instruments and Methods in Physics Research B34 (1988) 396-402. +*/ +#include <stdlib.h> +#include <string.h> +#include <math.h> + +#define MIN(x, y) (((x) < (y)) ? (x) : (y)) +#define MAX(x, y) (((x) > (y)) ? (x) : (y)) + +void lls(double *data, int size); +void lls_inv(double *data, int size); +void snip1d(double *data, int n_channels, int snip_width); +void snip1d_multiple(double *data, int n_channels, int snip_width, int n_spectra); +void lsdf(double *data, int size, int fwhm, double f, double A, double M, double ratio); + +void lls(double *data, int size) +{ + int i; + for (i=0; i< size; i++) + { + data[i] = log(log(sqrt(data[i]+1.0)+1.0)+1.0); + } +} + +void lls_inv(double *data, int size) +{ + int i; + double tmp; + for (i=0; i< size; i++) + { + /* slightly different than the published formula because + with the original formula: + + tmp = exp(exp(data[i]-1.0)-1.0); + data[i] = tmp * tmp - 1.0; + + one does not recover the original data */ + + tmp = exp(exp(data[i])-1.0)-1.0; + data[i] = tmp * tmp - 1.0; + } +} + +void lsdf(double *data, int size, int fwhm, double f, double A, double M, double ratio) +{ + int channel, i, j; + double L, R, S; + int width; + double dhelp; + + width = (int) (f * fwhm); + for (channel=width; channel<(size-width); channel++) + { + i = width; + while(i>0) + { + L=0; + R=0; + for(j=channel-i; j<channel; j++) + { + L += data[j]; + } + for(j=channel+1; j<channel+i; j++) + { + R += data[j]; + } + S = data[channel] + L + R; + if (S<M) + { + data[channel] = S /(2*i+1); + break; + } + dhelp = (R+1)/(L+1); + if ((dhelp < ratio) && (dhelp > (1/ratio))) + { + if (S<(A*sqrt(data[channel]))) + { + data[channel] = S /(2*i+1); + break; + } + } + i=i-1; + } + } +} + + +void snip1d(double *data, int n_channels, int snip_width) +{ + snip1d_multiple(data, n_channels, snip_width, 1); +} + +void snip1d_multiple(double *data, int n_channels, int snip_width, int n_spectra) +{ + int i; + int j; + int p; + int offset; + double *w; + + i = (int) (0.5 * snip_width); + /* lsdf(data, size, i, 1.5, 75., 10., 1.3); */ + + w = (double *) malloc(n_channels * sizeof(double)); + + for (j=0; j < n_spectra; j++) + { + offset = j * n_channels; + for (p = snip_width; p > 0; p--) + { + for (i=p; i<(n_channels - p); i++) + { + w[i] = MIN(data[i + offset], 0.5*(data[i + offset - p] + data[ i + offset + p])); + } + for (i=p; i<(n_channels - p); i++) + { + data[i+offset] = w[i]; + } + } + } + free(w); +} diff --git a/src/silx/math/fit/filters/src/snip2d.c b/src/silx/math/fit/filters/src/snip2d.c new file mode 100644 index 0000000..235759c --- /dev/null +++ b/src/silx/math/fit/filters/src/snip2d.c @@ -0,0 +1,96 @@ +#/*########################################################################## +# +# The PyMca X-Ray Fluorescence Toolkit +# +# Copyright (c) 2004-2014 European Synchrotron Radiation Facility +# +# This file is part of the PyMca X-ray Fluorescence Toolkit developed at +# the ESRF by the Software group. +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in +# all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +# THE SOFTWARE. +# +#############################################################################*/ +/* + Implementation of the algorithm SNIP in 2D described in + Miroslav Morhac et al. Nucl. Instruments and Methods in Physics Research A401 (1997) 113-132. +*/ +#include <stdlib.h> +#include <string.h> +#include <math.h> +#define MIN(x, y) (((x) < (y)) ? (x) : (y)) +#define MAX(x, y) (((x) > (y)) ? (x) : (y)) + +void lls(double *data, int size); +void lls_inv(double *data, int size); + +void snip2d(double *data, int nrows, int ncolumns, int width) +{ + int i, j; + int p; + int size; + double *w; + double P1, P2, P3, P4; + double S1, S2, S3, S4; + double dhelp; + int iminuspxncolumns; /* (i-p) * ncolumns */ + int ixncolumns; /* i * ncolumns */ + int ipluspxncolumns; /* (i+p) * ncolumns */ + + size = nrows * ncolumns; + w = (double *) malloc(size * sizeof(double)); + + for (p=width; p > 0; p--) + { + for (i=p; i<(nrows-p); i++) + { + iminuspxncolumns = (i-p) * ncolumns; + ixncolumns = i * ncolumns; + ipluspxncolumns = (i+p) * ncolumns; + for (j=p; j<(ncolumns-p); j++) + { + P4 = data[ iminuspxncolumns + (j-p)]; /* P4 = data[i-p][j-p] */ + S4 = data[ iminuspxncolumns + j]; /* S4 = data[i-p][j] */ + P2 = data[ iminuspxncolumns + (j+p)]; /* P2 = data[i-p][j+p] */ + S3 = data[ ixncolumns + (j-p)]; /* S3 = data[i][j-p] */ + S2 = data[ ixncolumns + (j+p)]; /* S2 = data[i][j+p] */ + P3 = data[ ipluspxncolumns + (j-p)]; /* P3 = data[i+p][j-p] */ + S1 = data[ ipluspxncolumns + j]; /* S1 = data[i+p][j] */ + P1 = data[ ipluspxncolumns + (j+p)]; /* P1 = data[i+p][j+p] */ + dhelp = 0.5*(P1+P3); + S1 = MAX(S1, dhelp) - dhelp; + dhelp = 0.5*(P1+P2); + S2 = MAX(S2, dhelp) - dhelp; + dhelp = 0.5*(P3+P4); + S3 = MAX(S3, dhelp) - dhelp; + dhelp = 0.5*(P2+P4); + S4 = MAX(S4, dhelp) - dhelp; + w[ixncolumns + j] = MIN(data[ixncolumns + j], 0.5 * (S1+S2+S3+S4) + 0.25 * (P1+P2+P3+P4)); + } + } + for (i=p; i<(nrows-p); i++) + { + ixncolumns = i * ncolumns; + for (j=p; j<(ncolumns-p); j++) + { + data[ixncolumns + j] = w[ixncolumns + j]; + } + } + } + free(w); +} diff --git a/src/silx/math/fit/filters/src/snip3d.c b/src/silx/math/fit/filters/src/snip3d.c new file mode 100644 index 0000000..cf48ee4 --- /dev/null +++ b/src/silx/math/fit/filters/src/snip3d.c @@ -0,0 +1,186 @@ +#/*########################################################################## +# +# The PyMca X-Ray Fluorescence Toolkit +# +# Copyright (c) 2004-2014 European Synchrotron Radiation Facility +# +# This file is part of the PyMca X-ray Fluorescence Toolkit developed at +# the ESRF by the Software group. +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in +# all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +# THE SOFTWARE. +# +#############################################################################*/ +/* + Implementation of the algorithm SNIP in 3D described in + Miroslav Morhac et al. Nucl. Instruments and Methods in Physics Research A401 (1997) 113-132. +*/ + +#include <stdlib.h> +#include <string.h> +#include <math.h> +#define MIN(x, y) (((x) < (y)) ? (x) : (y)) +#define MAX(x, y) (((x) > (y)) ? (x) : (y)) + +void lls(double *data, int size); +void lls_inv(double *data, int size); + +void snip3d(double *data, int nx, int ny, int nz, int width) +{ + int i, j, k; + int p; + int size; + double *w; + double P1, P2, P3, P4, P5, P6, P7, P8; + double R1, R2, R3, R4, R5, R6; + double S1, S2, S3, S4, S5, S6, S7, S8, S9, S10, S11, S12; + double dhelp; + long ioffset; + long iplus; + long imin; + long joffset; + long jplus; + long jmin; + + size = nx * ny * nz; + w = (double *) malloc(size * sizeof(double)); + + for (p=width; p > 0; p--) + { + for (i=p; i<(nx-p); i++) + { + ioffset = i * ny * nz; + iplus = (i + p) * ny * nz; + imin = (i - p) * ny * nz; + for (j=p; j<(ny-p); j++) + { + joffset = j * nz; + jplus = (j + p) * nz; + jmin = (j - p) * nz; + for (k=p; k<(nz-p); k++) + { + P1 = data[iplus + jplus + k-p]; /* P1 = data[i+p][j+p][k-p] */ + P2 = data[imin + jplus + k-p]; /* P2 = data[i-p][j+p][k-p] */ + P3 = data[iplus + jmin + k-p]; /* P3 = data[i+p][j-p][k-p] */ + P4 = data[imin + jmin + k-p]; /* P4 = data[i-p][j-p][k-p] */ + P5 = data[iplus + jplus + k+p]; /* P5 = data[i+p][j+p][k+p] */ + P6 = data[imin + jplus + k+p]; /* P6 = data[i-p][j+p][k+p] */ + P7 = data[imin + jmin + k+p]; /* P7 = data[i-p][j-p][k+p] */ + P8 = data[iplus + jmin + k+p]; /* P8 = data[i+p][j-p][k+p] */ + + S1 = data[iplus + joffset + k-p]; /* S1 = data[i+p][j][k-p] */ + S2 = data[ioffset + jmin + k-p]; /* S2 = data[i][j+p][k-p] */ + S3 = data[imin + joffset + k-p]; /* S3 = data[i-p][j][k-p] */ + S4 = data[ioffset + jmin + k-p]; /* S4 = data[i][j-p][k-p] */ + S5 = data[imin + joffset + k+p]; /* S5 = data[i-p][j][k+p] */ + S6 = data[ioffset + jplus + k+p]; /* S6 = data[i][j+p][k+p] */ + S7 = data[imin + joffset + k+p]; /* S7 = data[i-p][j][k+p] */ + S8 = data[ioffset + jmin + k+p]; /* S8 = data[i][j-p][k+p] */ + S9 = data[imin + jplus + k]; /* S9 = data[i-p][j+p][k] */ + S10 = data[imin + jmin + k]; /* S10 = data[i-p][j-p][k] */ + S11 = data[iplus + jmin + k]; /* S11 = data[i+p][j-p][k] */ + S12 = data[iplus + jplus + k]; /* S12 = data[i+p][j+p][k] */ + + R1 = data[ioffset + joffset + k-p]; /* R1 = data[i][j][k-p] */ + R2 = data[ioffset + joffset + k+p]; /* R2 = data[i][j][k+p] */ + R3 = data[imin + joffset + k]; /* R3 = data[i-p][j][k] */ + R4 = data[iplus + joffset + k]; /* R4 = data[i+p][j][k] */ + R5 = data[ioffset + jplus + k]; /* R5 = data[i][j+p][k] */ + R6 = data[ioffset + jmin + k]; /* R6 = data[i][j-p][k] */ + + dhelp = 0.5*(P1+P3); + S1 = MAX(S1, dhelp) - dhelp; + + dhelp = 0.5*(P1+P2); + S2 = MAX(S2, dhelp) - dhelp; + + dhelp = 0.5*(P2+P4); + S3 = MAX(S3, dhelp) - dhelp; + + dhelp = 0.5*(P3+P4); + S4 = MAX(S4, dhelp) - dhelp; + + dhelp = 0.5*(P5+P8); /* Different from paper (P5+P7) but according to drawing */ + S5 = MAX(S5, dhelp) - dhelp; + + dhelp = 0.5*(P5+P6); + S6 = MAX(S6, dhelp) - dhelp; + + dhelp = 0.5*(P6+P7); /* Different from paper (P6+P8) but according to drawing */ + S7 = MAX(S7, dhelp) - dhelp; + + dhelp = 0.5*(P7+P8); + S8 = MAX(S8, dhelp) - dhelp; + + dhelp = 0.5*(P2+P6); + S9 = MAX(S9, dhelp) - dhelp; + + dhelp = 0.5*(P4+P7); /* Different from paper (P4+P8) but according to drawing */ + S10 = MAX(S10, dhelp) - dhelp; + + dhelp = 0.5*(P3+P8); /* Different from paper (P1+P5) but according to drawing */ + S11 = MAX(S11, dhelp) - dhelp; + + dhelp = 0.5*(P1+P5); /* Different from paper (P3+P7) but according to drawing */ + S12 = MAX(S12, dhelp) - dhelp; + + /* The published formulae correspond to have: + P7 and P8 interchanged, and S11 and S12 interchanged + with respect to the published drawing */ + + dhelp = 0.5 * (S1+S2+S3+S4) + 0.25 * (P1+P2+P3+P4); + R1 = MAX(R1, dhelp) - dhelp; + + dhelp = 0.5 * (S5+S6+S7+S8) + 0.25 * (P5+P6+P7+P8); + R2 = MAX(R2, dhelp) - dhelp; + + dhelp = 0.5 * (S3+S7+S9+S10) + 0.25 * (P2+P4+P6+P7); /* Again same P7 and P8 change */ + R3 = MAX(R3, dhelp) - dhelp; + + dhelp = 0.5 * (S1+S5+S11+S12) + 0.25 * (P1+P3+P5+P8); /* Again same P7 and P8 change */ + R4 = MAX(R4, dhelp) - dhelp; + + dhelp = 0.5 * (S2+S6+S9+S12) + 0.25 * (P1+P2+P5+P6); /* Again same S11 and S12 change */ + R5 = MAX(R5, dhelp) - dhelp; + + dhelp = 0.5 * (S4+S8+S10+S11) + 0.25 * (P3+P4+P7+P8); /* Again same S11 and S12 change */ + R6 = MAX(R6, dhelp) - dhelp; + + dhelp = 0.5 * (R1 + R2 + R3 + R4 + R5 + R6) +\ + 0.25 * (S1 + S2 + S3 + S4 + S5 + S6) +\ + 0.25 * (S7 + S8 + S9 + S10 + S11 + S12) +\ + 0.125 * (P1 + P2 + P3 + P4 + P5 + P6 + P7 + P8); + w[ioffset + joffset + k] = MIN(data[ioffset + joffset + k], dhelp); + } + } + } + for (i=p; i<(nx-p); i++) + { + ioffset = i * ny * nz; + for (j=p; j<(ny-p); j++) + { + joffset = j * nz; + for (k=p; k<(nz-p); j++) + { + data[ioffset + joffset + k] = w[ioffset + joffset + k]; + } + } + } + } + free(w); +} diff --git a/src/silx/math/fit/filters/src/strip.c b/src/silx/math/fit/filters/src/strip.c new file mode 100644 index 0000000..dec0742 --- /dev/null +++ b/src/silx/math/fit/filters/src/strip.c @@ -0,0 +1,118 @@ +#/*########################################################################## +# Copyright (c) 2004-2016 European Synchrotron Radiation Facility +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in +# all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +# THE SOFTWARE. +# +#############################################################################*/ +/* + This file provides a background strip function, to isolate low frequency + background signal from a spectrum (and later substact it from the signal + to be left only with the peaks to be fitted). + + It is adapted from PyMca source file "SpecFitFuns.c". The main difference + with the original code is that this code does not handle the python + wrapping, which is done elsewhere using cython. + + Authors: V.A. Sole, P. Knobel + License: MIT + Last modified: 17/06/2016 +*/ + +#include <string.h> + +#include <stdio.h> + +/* strip(double* input, double c, long niter, double* output) + + The strip background is probably PyMca's most popular background model. + + In its simplest implementation it is just as an iterative procedure depending + on two parameters. These parameters are the strip background width w, and the + strip background number of iterations. At each iteration, if the contents of + channel i, y(i), is above the average of the contents of the channels at w + channels of distance, y(i-w) and y(i+w), y(i) is replaced by the average. + At the end of the process we are left with something that resembles a spectrum + in which the peaks have been "stripped". + + Parameters: + + - input: Input data array + - c: scaling factor applied to the average of y(i-w) and y(i+w) before + comparing to y(i) + - niter: number of iterations + - deltai: operator width (in number of channels) + - anchors: Array of anchors, indices of points that will not be + modified during the stripping procedure. + - output: output array + +*/ +int strip(double* input, long len_input, + double c, long niter, int deltai, + long* anchors, long len_anchors, + double* output) +{ + long iter_index, array_index, anchor_index, anchor; + int anchor_nearby_flag; + double t_mean; + + memcpy(output, input, len_input * sizeof(double)); + + if (deltai <=0) deltai = 1; + + if (len_input < (2*deltai+1)) return(-1); + + if (len_anchors > 0) { + for (iter_index = 0; iter_index < niter; iter_index++) { + for (array_index = deltai; array_index < len_input - deltai; array_index++) { + /* if index is within +- deltai of an anchor, don't do anything */ + anchor_nearby_flag = 0; + for (anchor_index=0; anchor_index<len_anchors; anchor_index++) + { + anchor = anchors[anchor_index]; + if (array_index > (anchor - deltai) && array_index < (anchor + deltai)) + { + anchor_nearby_flag = 1; + break; + } + } + /* skip this array_index index */ + if (anchor_nearby_flag) { + continue; + } + + t_mean = 0.5 * (input[array_index-deltai] + input[array_index+deltai]); + if (input[array_index] > (t_mean * c)) + output[array_index] = t_mean; + } + memcpy(input, output, len_input * sizeof(double)); + } + } + else { + for (iter_index = 0; iter_index < niter; iter_index++) { + for (array_index=deltai; array_index < len_input - deltai; array_index++) { + t_mean = 0.5 * (input[array_index-deltai] + input[array_index+deltai]); + + if (input[array_index] > (t_mean * c)) + output[array_index] = t_mean; + } + memcpy(input, output, len_input * sizeof(double)); + } + } + return(0); +} |