Mex/C++/6D,blobs<->sinos: added C++11 datastructures and parallel functions

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

zUtil_Cxx/6D_ops/gt6DBlobsToMultipleSinos_c.cpp

+/*
+ * gt6DBlobsToSingleSino_c.cpp
+ *
+ *  Created on: Oct 30, 2014
+ *      Author: vigano
+ */
+
+#include "DctMatlabData.h"
+
+static const char * transforms_error_id = "C_FUN:gt6DBlobsSinogramsTransforms:wrong_argument";
+
+void mexFunction( int nlhs, mxArray * plhs[], int nrhs, const mxArray * prhs[] )
+{
+  if (nrhs < 3) {
+    mexErrMsgIdAndTxt(transforms_error_id,
+        "Not enough arguments!");
+    return;
+  }
+
+  const mxArray * const sinos_sizes = prhs[0];
+  const mxArray * const blobs_cells = prhs[1];
+  const mxArray * const coeffs = prhs[2];
+
+  if (!mxIsNumeric(sinos_sizes)) {
+    mexErrMsgIdAndTxt(transforms_error_id,
+        "The first argument should be an array (Sinogram sizes)");
+    return;
+  }
+  if (!mxIsCell(blobs_cells)) {
+    mexErrMsgIdAndTxt(transforms_error_id,
+        "The second argument should be a Cell array (Blobs)");
+    return;
+  }
+
+  const mxClassID datatype = dct::DctMatlabData::get_class_of_cell(blobs_cells);
+
+  if (datatype == mxUNKNOWN_CLASS)
+  {
+    mexErrMsgIdAndTxt(transforms_error_id,
+        "The blobs need to be a non empty Cell array of coherent floating point types");
+    return;
+  }
+
+  size_t num_threads = std::thread::hardware_concurrency();
+  if (nrhs >= 4)
+  {
+    num_threads = mxGetScalar(prhs[3]);
+  }
+
+  switch (datatype)
+  {
+    case mxDOUBLE_CLASS:
+    {
+      dct::DctProjectionTransform<double, dct::MatlabAllocator<double> > transformer(num_threads);
+
+      try {
+        std::vector<dct::DctProjectionTransformCoefficients<double> > && metad = dct::DctMatlabData::import_coefficients_list<double>(coeffs);
+        transformer.set_transforms(metad);
+
+        dct::DctProjectionVector<double, dct::MatlabAllocator<double> > && sinos = dct::DctMatlabData::create_collection_from_dims<double>(sinos_sizes);
+        transformer.set_sinos(sinos);
+
+        dct::DctProjectionVector<double, dct::MatlabAllocator<double> > && blobs = dct::DctMatlabData::import_collection<double>((mxArray *)blobs_cells);
+        transformer.set_blobs(blobs);
+
+        transformer.check_consistency();
+
+      } catch (const dct::BasicException & e) {
+        mexErrMsgIdAndTxt(transforms_error_id, e.what());
+        return;
+      }
+
+      transformer.produce_sinograms();
+
+      plhs[0] = dct::DctMatlabData::produce_collection(transformer.get_sinos());
+      break;
+    }
+    case mxSINGLE_CLASS:
+    {
+      dct::DctProjectionTransform<float, dct::MatlabAllocator<float> > transformer(num_threads);
+
+      try {
+        std::vector<dct::DctProjectionTransformCoefficients<float> > && metad = dct::DctMatlabData::import_coefficients_list<float>(coeffs);
+        transformer.set_transforms(metad);
+
+        dct::DctProjectionVector<float, dct::MatlabAllocator<float> > && sinos = dct::DctMatlabData::create_collection_from_dims<float>(sinos_sizes);
+        transformer.set_sinos(sinos);
+
+        dct::DctProjectionVector<float, dct::MatlabAllocator<float> > && blobs = dct::DctMatlabData::import_collection<float>((mxArray *)blobs_cells);
+        transformer.set_blobs(blobs);
+
+        transformer.check_consistency();
+
+      } catch (const dct::BasicException & e) {
+        mexErrMsgIdAndTxt(transforms_error_id, e.what());
+        return;
+      }
+
+      transformer.produce_sinograms();
+
+      plhs[0] = dct::DctMatlabData::produce_collection(transformer.get_sinos());
+      break;
+    }
+    default:
+    {
+      mexErrMsgIdAndTxt(transforms_error_id,
+          "The argument needs to be a Cell array of floating point numbers");
+      return;
+    }
+  }
+}
+
+

zUtil_Cxx/6D_ops/gt6DMultipleSinosToBlobs_c.cpp

+/*
+ * gt6DBlobsToSingleSino_c.cpp
+ *
+ *  Created on: Oct 30, 2014
+ *      Author: vigano
+ */
+
+#include "DctMatlabData.h"
+
+static const char * transforms_error_id = "C_FUN:gt6DBlobsSinogramsTransforms:wrong_argument";
+
+void mexFunction( int nlhs, mxArray * plhs[], int nrhs, const mxArray * prhs[] )
+{
+  if (nrhs < 3) {
+    mexErrMsgIdAndTxt(transforms_error_id,
+        "Not enough arguments!");
+    return;
+  }
+
+  const mxArray * const blobs_cells = prhs[0];
+  const mxArray * const sinos_cells = prhs[1];
+  const mxArray * const coeffs = prhs[2];
+
+  if (!mxIsCell(sinos_cells)) {
+    mexErrMsgIdAndTxt(transforms_error_id,
+        "The first argument should be a Cell array (Sinograms)");
+    return;
+  }
+  if (!mxIsCell(blobs_cells)) {
+    mexErrMsgIdAndTxt(transforms_error_id,
+        "The second argument should be a Cell array (Blobs)");
+    return;
+  }
+
+  const mxClassID datatype = dct::DctMatlabData::get_class_of_cell(blobs_cells);
+  const mxClassID datatype_sinos = dct::DctMatlabData::get_class_of_cell(sinos_cells);
+
+  if (datatype == mxUNKNOWN_CLASS)
+  {
+    mexErrMsgIdAndTxt(transforms_error_id,
+        "The blobs need to be a non empty Cell array of coherent floating point types");
+    return;
+  }
+  if (datatype != datatype_sinos)
+  {
+    mexErrMsgIdAndTxt(transforms_error_id,
+        "The blobs need to be a non empty Cell array of a coherent floating point types with the sinograms");
+    return;
+  }
+
+  size_t num_threads = std::thread::hardware_concurrency();
+  if (nrhs >= 4)
+  {
+    num_threads = mxGetScalar(prhs[3]);
+  }
+
+  plhs[0] = mxCreateSharedDataCopy(blobs_cells);
+
+  switch (datatype)
+  {
+    case mxDOUBLE_CLASS:
+    {
+      dct::DctProjectionTransform<double, dct::MatlabAllocator<double> > transformer(num_threads);
+
+      try {
+        auto && metad = dct::DctMatlabData::import_coefficients_list<double>(coeffs);
+        transformer.set_transforms(metad);
+
+        auto && sinos = dct::DctMatlabData::import_collection<double>((mxArray *)sinos_cells);
+        transformer.set_sinos(sinos);
+
+        auto && blobs = dct::DctMatlabData::import_collection<double>(plhs[0]);
+        transformer.set_blobs(blobs);
+
+        transformer.check_consistency();
+
+      } catch (const dct::BasicException & e) {
+        mexErrMsgIdAndTxt(transforms_error_id, e.what());
+        return;
+      }
+
+      transformer.produce_blobs();
+      break;
+    }
+    case mxSINGLE_CLASS:
+    {
+      dct::DctProjectionTransform<float, dct::MatlabAllocator<float> > transformer(num_threads);
+
+      try {
+        auto && metad = dct::DctMatlabData::import_coefficients_list<float>(coeffs);
+        transformer.set_transforms(metad);
+
+        auto && sinos = dct::DctMatlabData::import_collection<float>((mxArray *)sinos_cells);
+        transformer.set_sinos(sinos);
+
+        auto && blobs = dct::DctMatlabData::import_collection<float>(plhs[0]);
+        transformer.set_blobs(blobs);
+
+        transformer.check_consistency();
+
+      } catch (const dct::BasicException & e) {
+        mexErrMsgIdAndTxt(transforms_error_id, e.what());
+        return;
+      }
+
+      transformer.produce_blobs();
+      break;
+    }
+    default:
+    {
+      mexErrMsgIdAndTxt(transforms_error_id,
+          "The argument needs to be a Cell array of floating point numbers");
+      return;
+    }
+  }
+}
+
+

zUtil_Cxx/include/DctMatlabData.h

+/*
+ * DctMatlabData.h
+ *
+ *  Created on: Jan 23, 2017
+ *      Author: vigano
+ */
+
+#ifndef ZUTIL_CXX_INCLUDE_DCTMATLABDATA_H_
+#define ZUTIL_CXX_INCLUDE_DCTMATLABDATA_H_
+
+#include <mex.h>
+
+extern "C"
+{
+  bool mxUnshareArray(mxArray *array_ptr, bool noDeepCopy);
+  mxArray *mxCreateSharedCopy(const mxArray *pr);
+  mxArray *mxCreateSharedDataCopy(const mxArray *pr);
+  mxArray *mxUnreference(const mxArray *pr);
+}
+
+#include "DctData.h"
+
+namespace dct {
+
+  template<typename Type>
+  class MatlabAllocator {
+  public:
+    Type * allocate(const size_t & num_elements){
+      return (Type *) mxMalloc(num_elements * sizeof(Type));
+    }
+    void deallocate(Type * p) { mxFree(p); }
+  };
+
+  template<typename Type>
+  class DctMatlabDatatype {
+  public:
+    static constexpr mxClassID get_matlab_type() noexcept;
+  };
+
+  template<>
+  inline constexpr mxClassID
+  DctMatlabDatatype<double>::get_matlab_type() noexcept { return mxDOUBLE_CLASS; }
+
+  template<>
+  inline constexpr mxClassID
+  DctMatlabDatatype<float>::get_matlab_type() noexcept { return mxSINGLE_CLASS; }
+
+  class DctMatlabData {
+  public:
+    static mxClassID get_class_of_cell(const mxArray * const cell);
+
+    template<typename Type>
+    static DctProjectionVector<Type, MatlabAllocator<Type> > import_collection(mxArray * const);
+
+    template<typename Type>
+    static DctProjectionVector<Type, MatlabAllocator<Type> > create_collection_from_dims(const mxArray * const);
+
+    template<typename Type>
+    static void import_coefficients(DctProjectionTransformCoefficients<Type> & out, const mxArray * const);
+    template<typename Type>
+    static std::vector<DctProjectionTransformCoefficients<Type> > import_coefficients_list(const mxArray * const);
+
+    static mxArray * vol_allocate(const size_t & num_dims, const size_t * const dims, const mxClassID & type);
+
+    template<typename Type>
+    static mxArray * produce_collection(const DctProjectionVector<Type, MatlabAllocator<Type> > &);
+  };
+
+} /* namespace dct */
+
+
+// Ugly trick to not deal with instantiation
+
+#include "exceptions.h"
+
+#include <emlrt.h>
+
+namespace dct {
+
+  mxClassID
+  DctMatlabData::get_class_of_cell(const mxArray * const cell_array)
+  {
+    const mwSize numElems = mxGetNumberOfElements(cell_array);
+    if (numElems > 0) {
+      const mxArray* const firstArray = mxGetCell(cell_array, 0);
+      const mxClassID first_class = mxGetClassID(firstArray);
+
+      if (first_class == mxSINGLE_CLASS || first_class == mxDOUBLE_CLASS) {
+        for (mwIndex count = 1; count < numElems; count++) {
+          const mxArray* const array = mxGetCell(cell_array, count);
+
+          if (first_class != mxGetClassID(array)) {
+            // Damn mixed types!
+            return mxUNKNOWN_CLASS;
+          }
+        }
+        // OK it is a float type!
+        return first_class;
+      } else {
+        // It is still an error!
+        return first_class;
+      }
+    } else {
+      return mxUNKNOWN_CLASS;
+    }
+  }
+
+  template<typename Type>
+  DctProjectionVector<Type, MatlabAllocator<Type> >
+  DctMatlabData::import_collection(mxArray * const mxdata)
+  {
+    DctProjectionVector<Type, MatlabAllocator<Type> > output;
+
+    bool is_cell = mxIsCell(mxdata);
+    const mxClassID datatype = is_cell
+        ? get_class_of_cell(mxdata) : mxGetClassID(mxdata);
+
+    if (DctMatlabDatatype<Type>::get_matlab_type() != datatype) {
+      throw WrongArgumentException("Input array should be of the same type as the object type");
+    }
+
+    if (is_cell) {
+      const mwSize tot_cells = mxGetNumberOfElements(mxdata);
+      output.reserve(tot_cells);
+      for (mwIndex num_cell = 0; num_cell < tot_cells; num_cell++) {
+        mxArray * cell = mxGetCell(mxdata, num_cell);
+
+        const mwSize * dims = mxGetDimensions(cell);
+        const mwSize num_dims = mxGetNumberOfDimensions(cell);
+        Type * const data = (Type *)mxGetData(cell);
+
+        output.emplace_back(dims, num_dims, 1, data);
+      }
+    } else {
+      const mwSize * dims = mxGetDimensions(mxdata);
+      const mwSize num_dims = mxGetNumberOfDimensions(mxdata);
+      Type * const data = (Type *)mxGetData(mxdata);
+
+      output.emplace_back(dims, num_dims, 1, data);
+    }
+
+    return output;
+  }
+
+  template<typename Type>
+  DctProjectionVector<Type, MatlabAllocator<Type> >
+  DctMatlabData::create_collection_from_dims(const mxArray * const mxdata)
+  {
+    DctProjectionVector<Type, MatlabAllocator<Type> > output;
+
+    const mwSize num_dims_dims = mxGetNumberOfDimensions(mxdata);
+    const mwSize * dims_dims = mxGetDimensions(mxdata);
+
+    const mwSize num_dims_sinos = dims_dims[1];
+
+    if (!(mxIsNumeric(mxdata) && num_dims_dims == 2 && num_dims_sinos == 3)) {
+      throw WrongArgumentException("Array dimensions should be in a Mx3 matrix, where M is the number of volumes");
+    }
+
+    const mxClassID datatype = mxGetClassID(mxdata);
+
+    if (mxDOUBLE_CLASS != datatype) {
+      throw WrongArgumentException("Array dimensions should be 'double'");
+    }
+//    if (mxDOUBLE_CLASS != datatype && mxSINGLE_CLASS != datatype) {
+//      throw WrongArgumentException("Array dimensions should either be 'single' or 'double'");
+//    }
+
+    const mwSize tot_projs = dims_dims[0];
+    output.reserve(tot_projs);
+
+    const double * const dims_sinos = (const double *) mxGetData(mxdata);
+
+    for (mwIndex num_proj = 0; num_proj < tot_projs; num_proj++) {
+      output.emplace_back(dims_sinos + num_proj, num_dims_sinos, tot_projs);
+    }
+
+    output.allocate();
+
+    return output;
+  }
+
+  template<typename Type>
+  void
+  DctMatlabData::import_coefficients(DctProjectionTransformCoefficients<Type> & output, const mxArray * const coeff)
+  {
+    if (!mxIsStruct(coeff)) {
+      throw WrongArgumentException("Invalid Coeffs: third argument should be a structure!");
+    }
+    const mxArray * const sino_offsets_a = mxGetField(coeff, 0, "sino_offsets");
+    const mxArray * const blob_offsets_a = mxGetField(coeff, 0, "blob_offsets");
+    const mxArray * const proj_offsets_a = mxGetField(coeff, 0, "proj_offsets");
+    const mxArray * const proj_coeffs_a = mxGetField(coeff, 0, "proj_coeffs");
+
+    if (!sino_offsets_a) {
+      throw WrongArgumentException("Invalid Coeffs: 'sino_offsets' is missing!");
+    }
+    if (!blob_offsets_a) {
+      throw WrongArgumentException("Invalid Coeffs: 'blob_offsets' is missing!");
+    }
+    if (!proj_offsets_a) {
+      throw WrongArgumentException("Invalid Coeffs: 'proj_offsets' is missing!");
+    }
+    if (!proj_coeffs_a) {
+      throw WrongArgumentException("Invalid Coeffs: 'proj_coeffs' is missing!");
+    }
+
+    const mwSize num_elems = mxGetNumberOfElements(sino_offsets_a);
+    if (num_elems != mxGetNumberOfElements(blob_offsets_a)) {
+      throw WrongArgumentException("Invalid Coeffs: number of 'blob_offsets' is different from 'sino_offsets'!");
+    }
+    if (num_elems != mxGetNumberOfElements(proj_offsets_a)) {
+      throw WrongArgumentException("Invalid Coeffs: number of 'proj_offsets' is different from 'sino_offsets'!");
+    }
+    if (num_elems != mxGetNumberOfElements(proj_coeffs_a)) {
+      throw WrongArgumentException("Invalid Coeffs: number of 'proj_coeffs' is different from 'sino_offsets'!");
+    }
+
+    if (mxDOUBLE_CLASS != mxGetClassID(sino_offsets_a)) {
+      throw WrongArgumentException("Invalid Coeffs: type of 'sino_offsets' should be 'double'!");
+    }
+    if (mxDOUBLE_CLASS != mxGetClassID(blob_offsets_a)) {
+      throw WrongArgumentException("Invalid Coeffs: type of 'blob_offsets' should be 'double'!");
+    }
+    if (mxDOUBLE_CLASS != mxGetClassID(proj_offsets_a)) {
+      throw WrongArgumentException("Invalid Coeffs: type of 'proj_offsets' should be 'double'!");
+    }
+    if (mxDOUBLE_CLASS != mxGetClassID(proj_coeffs_a)) {
+      throw WrongArgumentException("Invalid Coeffs: type of 'proj_coeffs' should be 'double'!");
+    }
+
+    output.resize(num_elems);
+
+    const double * const blobs_pos = (const double *) mxGetData(blob_offsets_a);
+    const double * const sinos_pos = (const double *) mxGetData(sino_offsets_a);
+    const double * const projs_pos = (const double *) mxGetData(proj_offsets_a);
+
+    // Converting from Matlab's indexing, starting from 1, to C's indexing starting from 0
+    OpUnitarySumScalar<double> op(-1);
+
+    std::transform(blobs_pos, blobs_pos + num_elems, output.blob_pos.begin(), op);
+    std::transform(sinos_pos, sinos_pos + num_elems, output.sino_pos.begin(), op);
+    std::transform(projs_pos, projs_pos + num_elems, output.proj_pos.begin(), op);
+
+    const double * const coeffs = (const double *) mxGetData(proj_coeffs_a);
+    std::copy(coeffs, coeffs + num_elems, output.coeff.begin());
+  }
+
+  template<typename Type>
+  std::vector<DctProjectionTransformCoefficients<Type> >
+  DctMatlabData::import_coefficients_list(const mxArray * const coeffs)
+  {
+    std::vector<DctProjectionTransformCoefficients<Type> > output;
+
+    const mwSize num_coeff_structs = mxGetNumberOfElements(coeffs);
+    output.resize(num_coeff_structs);
+
+    for (mwIndex num_struct = 0; num_struct < num_coeff_structs; num_struct++) {
+      const mxArray * const coeff_cell = mxGetCell(coeffs, num_struct);
+      DctMatlabData::import_coefficients<Type>(output[num_struct], coeff_cell);
+    }
+
+    return output;
+  }
+
+  mxArray *
+  DctMatlabData::vol_allocate(const size_t & num_dims, const size_t * const dims, const mxClassID & type)
+  {
+    #if defined(EMLRT_VERSION_INFO) && EMLRT_VERSION_INFO >= 0x2015a && !defined(MX_COMPAT_32)
+        return mxCreateUninitNumericArray(num_dims, (mwSize *)dims, type, mxREAL);
+    #else
+        const mwSize zero_dims[2] = {0, 0};
+        mxArray * const xOutput = mxCreateNumericArray(2, zero_dims, type, mxREAL);
+        mxSetDimensions(xOutput, dims, num_dims);
+        const mwSize num_elems = mxGetNumberOfElements(xOutput);
+        const mwSize elem_size = mxGetElementSize(xOutput);
+        mxSetData(xOutput, mxMalloc(elem_size * num_elems));
+        return xOutput;
+    #endif
+  }
+
+  template<typename Type>
+  mxArray *
+  DctMatlabData::produce_collection(const DctProjectionVector<Type, MatlabAllocator<Type> > & vec)
+  {
+    const mwSize zero_dims[2] = {0, 0};
+
+    const mxClassID datatype = DctMatlabDatatype<Type>::get_matlab_type();
+
+    const size_t tot_cells = vec.size();
+
+    if (tot_cells > 1) {
+      mxArray * const output = mxCreateCellMatrix(tot_cells, 1);
+
+      for (mwIndex num_cell = 0; num_cell < tot_cells; num_cell++) {
+        const size_t * const dims = vec[num_cell].get_dims();
+
+        mxArray * const cell = mxCreateNumericArray(2, zero_dims, datatype, mxREAL);
+        mxSetDimensions(cell, dims, 3);
+        mxSetData(cell, vec[num_cell].get_data());
+
+        mxSetCell(output, num_cell, cell);
+      }
+
+      return output;
+
+    } else {
+      const size_t * const dims = vec[0].get_dims();
+
+      mxArray * const output = mxCreateNumericArray(2, zero_dims, datatype, mxREAL);
+
+      mxSetDimensions(output, dims, 3);
+      mxSetData(output, vec[0].get_data());
+
+      return output;
+    }
+  }
+
+} /* namespace dct */
+
+#endif /* ZUTIL_CXX_INCLUDE_DCTMATLABDATA_H_ */

zUtil_Cxx/include/exceptions.h

+/*
+ * exceptions.h
+ *
+ *  Created on: Jan 22, 2017
+ *      Author: vigano
+ */
+
+#ifndef ZUTIL_CXX_INCLUDE_EXCEPTIONS_H_
+#define ZUTIL_CXX_INCLUDE_EXCEPTIONS_H_
+
+#include <stdexcept>
+#include <string>
+
+namespace dct {
+
+  class BasicException : public std::exception {
+    std::string message;
+  public:
+    BasicException() = default;
+    BasicException(const char * _mess) : message(_mess) { }
+    BasicException(const std::string& _mess) : message(_mess) { }
+    BasicException(const BasicException& _ex) : message(_ex.message) { }
+
+    virtual ~BasicException() throw() = default;
+
+    virtual const char * what() const throw() { return message.c_str(); }
+    const std::string& getMessage() const throw() { return message; }
+
+    virtual void setMessage(const std::string & _mess) { message = _mess; }
+    virtual void appendMessage(const std::string & _mess) { message += _mess; }
+    virtual void prefixMessage(const std::string & _mess) { message = _mess + message; }
+  };
+
+  class WrongArgumentException : public BasicException {
+  public:
+    WrongArgumentException() = default;
+    WrongArgumentException(const char * _mess) : BasicException(_mess) { }
+    WrongArgumentException(const std::string & _mess) : BasicException(_mess) { }
+  };
+
+}  // namespace dct
+
+#endif /* ZUTIL_CXX_INCLUDE_EXCEPTIONS_H_ */

zUtil_Cxx/include/parallelization.h

+/*
+ * parallelization.h
+ *
+ *  Created on: Jan 25, 2017
+ *      Author: vigano
+ */
+
+#ifndef ZUTIL_CXX_INCLUDE_PARALLELIZATION_H_
+#define ZUTIL_CXX_INCLUDE_PARALLELIZATION_H_
+
+
+#include <queue>
+#include <memory>
+#include <mutex>
+#include <condition_variable>
+
+#include <thread>
+#include <functional>
+#include <atomic>
+
+namespace dct {
+
+  template<typename Type>
+    class ThreadSafeQueue {
+    private:
+      mutable std::mutex mut;
+      std::queue<Type> data_queue;
+      std::condition_variable data_cond;
+
+    public:
+      ThreadSafeQueue() { }
+      ThreadSafeQueue(ThreadSafeQueue const & other) {
+        std::lock_guard<std::mutex> lk(other.mut);
+        data_queue = other.data_queue;
+      }
+
+      void push(const Type & new_value) {
+        std::lock_guard<std::mutex> lk(mut);
+        data_queue.emplace(new_value);
+        data_cond.notify_one();
+      }
+
+      void push(const std::vector<Type> & new_values) {
+        std::lock_guard<std::mutex> lk(mut);
+        for (const auto & item : new_values) {
+          data_queue.emplace(item);
+        }
+        data_cond.notify_one();
+      }
+
+      bool try_pop(Type & value) {
+        std::lock_guard<std::mutex> lk(mut);
+        if (data_queue.empty()) {
+          return false;
+        } else {
+          value = data_queue.front();
+          data_queue.pop();
+          return true;
+        }
+      }
+
+      bool empty() const {
+        std::lock_guard<std::mutex> lk(mut);
+        return data_queue.empty();
+      }
+    };
+
+  class JoinThreads {
+    std::vector<std::thread> & threads;
+  public:
+    explicit
+    JoinThreads(std::vector<std::thread> & threads_)
+    : threads(threads_) { }
+
+    ~JoinThreads() {
+      for (size_t t = 0; t < threads.size(); t++) {
+        if (threads[t].joinable()) {
+          threads[t].join();
+        }
+      }
+    }
+  };
+
+  class ThreadPool {
+    std::atomic_bool done;
+
+    std::vector<std::thread> threads;
+    JoinThreads joiner;
+
+    ThreadSafeQueue< std::function<void()> > work_queue;
+
+    void
+    worker_thread()
+    {
+      while (!done) {
+        std::function < void() > task;
+        if (work_queue.try_pop(task)) {
+          task();
+        } else {
+          std::this_thread::yield();
+        }
+      }
+    }
+
+  public:
+    ThreadPool(const size_t & tot_threads = std::thread::hardware_concurrency())
+    : done(false), joiner(threads)
+    {
+      try {
+        for (size_t t = 0; t < tot_threads; t++) {
+          threads.push_back(std::thread(&ThreadPool::worker_thread, this));
+        }
+      } catch (...) {
+        done = true;
+        throw;
+      }
+    }
+
+    ~ThreadPool() {
+      done = true;
+    }
+
+    template<typename FunctionType>
+    void
+    submit_multiple(const std::vector<FunctionType> & fs) {
+      work_queue.push(fs);
+    }
+
+    template<typename FunctionType>
+    void
+    submit(FunctionType & f) {
+      work_queue.push(std::function < void() > (f));
+    }
+
+    bool work_queue_is_empty() const { return work_queue.empty(); }
+
+    const size_t get_num_threads() const noexcept { return threads.size(); }
+  };
+
+}  // namespace dct
+
+
+
+#endif /* ZUTIL_CXX_INCLUDE_PARALLELIZATION_H_ */

zUtil_Python/dct_utils_platform.py

@@ -215,7 +215,7 @@ class UtilsLinux(UtilsInterface):
 
     def resetMatlabMexopts(self, matlabVersion, matlabPath = ""):
         if (self.isMatlabOlderThan(matlabVersion, "2014a")):
-            searchPatterns = [('CXXOPTIMFLAGS="-O -DNDEBUG"', 'CXXOPTIMFLAGS="-O -DNDEBUG -fopenmp"'), \
+            searchPatterns = [('CXXOPTIMFLAGS="-O -DNDEBUG"', 'CXXOPTIMFLAGS="-O -DNDEBUG -fopenmp -std=c++11"'), \
                               ('LDOPTIMFLAGS="-O"', 'LDOPTIMFLAGS="-O -fopenmp"')]
         else:
             searchPatterns = [("COPTIMFLAGS='-O -DNDEBUG'", "COPTIMFLAGS='-O -DNDEBUG -fopenmp'"), \