SoftThresholding: added fast soft-thresholding implementation

Can't be pushed upstream yet, because it doesn't already handle the pathological cases.

Signed-off-by: Nicola Vigano <nicola.vigano@esrf.fr>

git-svn-id: https://svn.code.sf.net/p/dct/code/trunk@675 4c865b51-4357-4376-afb4-474e03ccb993

zUtil_Cxx/CompressiveSensing/gtSoftThreshold.cpp

+/*
+ * gtSoftThreshold.cpp
+ *
+ *  Created on: Jul 17, 2012
+ *
+ * Nicola Vigano', 2012, ID11 @ ESRF vigano@esrf.eu
+ */
+
+#include "internal_array_defs.h"
+
+#include "mex.h"
+
+void mexFunction( int nlhs, mxArray * plhs[], int nrhs, const mxArray * prhs[] )
+{
+  if (!(nrhs == 3 || nrhs == 2)) {
+    mexErrMsgIdAndTxt("C_FUN:gtSoftThreshold:wrong_argument",
+        "Not enough arguments! You need to provide an array, a threshold and "
+        "optionally a flag that tells if making a copy is necessary (default: true)");
+  }
+  if (!(mxIsDouble(prhs[1]) || mxIsSingle(prhs[1]))) {
+    mexErrMsgIdAndTxt("C_FUN:gtSoftThreshold:wrong_argument",
+        "The second parameter should be a double containing the threshold");
+  }
+
+  mxLogical copy = true;
+  if (nrhs == 3) {
+    if (!mxIsLogicalScalar(prhs[2])) {
+      mexErrMsgIdAndTxt("C_FUN:gtSoftThreshold:wrong_argument",
+          "The third parameter should be a logical");
+    } else {
+      copy = mxIsLogicalScalarTrue(prhs[2]);
+    }
+  }
+
+  const double * const __restrict inData = mxGetPr(prhs[0]);
+  double * const __restrict outData = mxGetPr(prhs[1]);
+
+  soft_threshold<mwSize, double> thr(*mxGetPr(prhs[1]));
+
+  const mwSize numElems = mxGetNumberOfElements(prhs[0]);
+  if (copy) {
+    const mwSize numDims = mxGetNumberOfDimensions(prhs[0]);
+    const mwSize * dims = mxGetDimensions(prhs[0]);
+    plhs[0] = mxCreateNumericArray(numDims, dims, mxDOUBLE_CLASS, mxREAL);
+
+    thr(mxGetPr(plhs[0]), mxGetPr(prhs[0]), numElems);
+  } else {
+    plhs[0] = (mxArray *) prhs[0];
+
+    thr(mxGetPr(plhs[0]), numElems);
+  }
+}
+
+
+
+

zUtil_Cxx/include/internal_array_defs.h

+/*
+ * internal_array_defs.h
+ *
+ *  Created on: Jul 17, 2012
+ *
+ * Nicola Vigano', 2012, ID11 @ ESRF vigano@esrf.eu
+ */
+
+#ifndef INTERNAL_ARRAY_DEFS_H_
+#define INTERNAL_ARRAY_DEFS_H_
+
+#include <cmath>
+#include <omp.h>
+
+using namespace std;
+
+template<typename SizeType, typename Type>
+class soft_threshold {
+  const Type thr;
+public:
+  soft_threshold(const Type & _thr) : thr(_thr) { }
+
+  void operator()(Type * const __restrict array, const SizeType & numElems);
+  void operator()(Type * const __restrict outArray,
+      const Type * const __restrict inArray, const SizeType & numElems);
+};
+
+template<typename SizeType, typename Type>
+inline void
+soft_threshold<SizeType, Type>::operator()(Type * const __restrict array,
+    const SizeType & numElems)
+{
+#pragma omp prallel for
+  for(SizeType count = 0; count < numElems; count++)
+  {
+    const Type & oldVal = array[count];
+    const Type temp = oldVal - thr;
+    array[count] = copysign((temp + fabs(temp)) / 2, oldVal);
+  }
+}
+
+template<typename SizeType, typename Type>
+inline void
+soft_threshold<SizeType, Type>::operator()(Type * const __restrict outArray,
+    const Type * const __restrict inArray, const SizeType & numElems)
+{
+#pragma omp prallel for
+  for(SizeType count = 0; count < numElems; count++)
+  {
+    const Type & oldVal = inArray[count];
+    const Type temp = oldVal - thr;
+    outArray[count] = copysign((temp + fabs(temp)) / 2, oldVal);
+  }
+}
+
+
+#endif /* INTERNAL_ARRAY_DEFS_H_ */

zUtil_Imaging/gtSoftThreshold.m → zUtil_Imaging/gtSoftThreshold2.m

-function img = gtSoftThreshold(img, thr)
+function img = gtSoftThreshold2(img, thr)
 % GTSOFTTHRESHOLD Applies soft thresholding to the image
 %   img = gtSoftThreshold(img, thr)