summaryrefslogtreecommitdiff
path: root/silx/math/fit/filters
diff options
context:
space:
mode:
Diffstat (limited to 'silx/math/fit/filters')
-rw-r--r--silx/math/fit/filters/include/filters.h45
-rw-r--r--silx/math/fit/filters/src/smoothnd.c317
-rw-r--r--silx/math/fit/filters/src/snip1d.c149
-rw-r--r--silx/math/fit/filters/src/snip2d.c96
-rw-r--r--silx/math/fit/filters/src/snip3d.c186
-rw-r--r--silx/math/fit/filters/src/strip.c118
6 files changed, 0 insertions, 911 deletions
diff --git a/silx/math/fit/filters/include/filters.h b/silx/math/fit/filters/include/filters.h
deleted file mode 100644
index 1ee9a95..0000000
--- a/silx/math/fit/filters/include/filters.h
+++ /dev/null
@@ -1,45 +0,0 @@
-/*##########################################################################
-# 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/silx/math/fit/filters/src/smoothnd.c b/silx/math/fit/filters/src/smoothnd.c
deleted file mode 100644
index cb96961..0000000
--- a/silx/math/fit/filters/src/smoothnd.c
+++ /dev/null
@@ -1,317 +0,0 @@
-#/*##########################################################################
-#
-# 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/silx/math/fit/filters/src/snip1d.c b/silx/math/fit/filters/src/snip1d.c
deleted file mode 100644
index 994a272..0000000
--- a/silx/math/fit/filters/src/snip1d.c
+++ /dev/null
@@ -1,149 +0,0 @@
-#/*##########################################################################
-# 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/silx/math/fit/filters/src/snip2d.c b/silx/math/fit/filters/src/snip2d.c
deleted file mode 100644
index 235759c..0000000
--- a/silx/math/fit/filters/src/snip2d.c
+++ /dev/null
@@ -1,96 +0,0 @@
-#/*##########################################################################
-#
-# 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/silx/math/fit/filters/src/snip3d.c b/silx/math/fit/filters/src/snip3d.c
deleted file mode 100644
index cf48ee4..0000000
--- a/silx/math/fit/filters/src/snip3d.c
+++ /dev/null
@@ -1,186 +0,0 @@
-#/*##########################################################################
-#
-# 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/silx/math/fit/filters/src/strip.c b/silx/math/fit/filters/src/strip.c
deleted file mode 100644
index dec0742..0000000
--- a/silx/math/fit/filters/src/strip.c
+++ /dev/null
@@ -1,118 +0,0 @@
-#/*##########################################################################
-# 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);
-}
generated by cgit v1.2.3 (git 2.25.1) at 2024-05-29 00:35:37 +0000