diff --git a/Modules/PhotoacousticsLib/files.cmake b/Modules/PhotoacousticsLib/files.cmake index 9452595421..2e2d84c31d 100644 --- a/Modules/PhotoacousticsLib/files.cmake +++ b/Modules/PhotoacousticsLib/files.cmake @@ -1,58 +1,64 @@ SET(H_FILES include/mitkPAPropertyCalculator.h include/mitkPAVector.h include/mitkPATissueGeneratorParameters.h include/mitkPAInSilicoTissueVolume.h include/mitkPATissueGenerator.h include/mitkPAVesselTree.h include/mitkPAVessel.h include/mitkPAVesselMeanderStrategy.h include/mitkPANoiseGenerator.h include/mitkPAVolume.h include/mitkPAComposedVolume.h include/mitkPASlicedVolumeGenerator.h include/mitkPAProbe.h include/mitkPALightSource.h include/mitkPAIOUtil.h include/mitkPAMonteCarloThreadHandler.h include/mitkPASimulationBatchGenerator.h include/mitkPAFluenceYOffsetPair.h include/mitkPAVolumeManipulator.h include/mitkPAVesselProperties.h include/mitkPASimulationBatchGeneratorParameters.h include/mitkPAExceptions.h include/mitkPASpectralUnmixingFilterBase.h include/mitkPALinearSpectralUnmixingFilter.h include/mitkPASpectralUnmixingSO2.h + include/mitkPASpectralUnmixingFilterLagrange.h + include/mitkPASpectralUnmixingFilterSimplex.h + include/mitkPASpectralUnmixingFilterVigra.h ) set(CPP_FILES Domain/Vessel/mitkPAVesselTree.cpp Domain/Vessel/mitkPAVessel.cpp Domain/Vessel/mitkPAVesselMeanderStrategy.cpp Domain/Vessel/mitkPAVesselProperties.cpp Domain/Volume/mitkPAInSilicoTissueVolume.cpp Domain/Volume/mitkPAVolume.cpp Domain/Volume/mitkPAComposedVolume.cpp Domain/Volume/mitkPAFluenceYOffsetPair.cpp Generator/mitkPATissueGenerator.cpp Generator/mitkPANoiseGenerator.cpp Generator/mitkPASlicedVolumeGenerator.cpp Generator/mitkPASimulationBatchGenerator.cpp Generator/mitkPASimulationBatchGeneratorParameters.cpp IO/mitkPAIOUtil.cpp Utils/mitkPAPropertyCalculator.cpp Utils/mitkPAVector.cpp Utils/mitkPATissueGeneratorParameters.cpp Utils/mitkPAVolumeManipulator.cpp Utils/ProbeDesign/mitkPAProbe.cpp Utils/ProbeDesign/mitkPALightSource.cpp Utils/Thread/mitkPAMonteCarloThreadHandler.cpp SUFilter/mitkPASpectralUnmixingFilterBase.cpp SUFilter/mitkPALinearSpectralUnmixingFilter.cpp SUFilter/mitkPASpectralUnmixingSO2.cpp + SUFilter/mitkPASpectralUnmixingFilterSimplex.cpp + SUFilter/mitkPASpectralUnmixingFilterVigra.cpp + SUFilter/mitkPASpectralUnmixingFilterLagrange.cpp ) set(RESOURCE_FILES spectralLIB.dat ) diff --git a/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilterLagrange.h b/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilterLagrange.h new file mode 100644 index 0000000000..def0212178 --- /dev/null +++ b/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilterLagrange.h @@ -0,0 +1,41 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#ifndef MITKPHOTOACOUSTICSPECTRALUNMIXINGFILTERLAGRANGE_H +#define MITKPHOTOACOUSTICSPECTRALUNMIXINGFILTERLAGRANGE_H + +#include +#include + +namespace mitk { + namespace pa { + class MITKPHOTOACOUSTICSLIB_EXPORT SpectralUnmixingFilterLagrange : public SpectralUnmixingFilterBase + { + public: + + mitkClassMacro(SpectralUnmixingFilterLagrange, SpectralUnmixingFilterBase) + //itkFactorylessNewMacro(Self) + + protected: + SpectralUnmixingFilterLagrange(); + virtual ~SpectralUnmixingFilterLagrange(); + + private: + + }; + } +} +#endif // MITKPHOTOACOUSTICSPECTRALUNMIXINGFILTERLAGRANGE_H diff --git a/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilterSimplex.h b/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilterSimplex.h new file mode 100644 index 0000000000..47bb93f820 --- /dev/null +++ b/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilterSimplex.h @@ -0,0 +1,41 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#ifndef MITKPHOTOACOUSTICSPECTRALUNMIXINGFILTERSIMPLEX_H +#define MITKPHOTOACOUSTICSPECTRALUNMIXINGFILTERSIMPLEX_H + +#include +#include + +namespace mitk { + namespace pa { + class MITKPHOTOACOUSTICSLIB_EXPORT SpectralUnmixingFilterSimplex : public SpectralUnmixingFilterBase + { + public: + + mitkClassMacro(SpectralUnmixingFilterSimplex, SpectralUnmixingFilterBase) + //itkFactorylessNewMacro(Self) + + protected: + SpectralUnmixingFilterSimplex(); + virtual ~SpectralUnmixingFilterSimplex(); + + private: + + }; + } +} +#endif // MITKPHOTOACOUSTICSPECTRALUNMIXINGFILTERSIMPLEX_H diff --git a/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilterVigra.h b/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilterVigra.h new file mode 100644 index 0000000000..b25766cd86 --- /dev/null +++ b/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilterVigra.h @@ -0,0 +1,46 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#ifndef MITKPHOTOACOUSTICSPECTRALUNMIXINGFILTERVIGRA_H +#define MITKPHOTOACOUSTICSPECTRALUNMIXINGFILTERVIGRA_H + +#include +#include + +namespace mitk { + namespace pa { + class MITKPHOTOACOUSTICSLIB_EXPORT SpectralUnmixingFilterVigra : public SpectralUnmixingFilterBase + { + public: + + mitkClassMacro(SpectralUnmixingFilterVigra, SpectralUnmixingFilterBase) + //itkFactorylessNewMacro(Self) + + protected: + SpectralUnmixingFilterVigra(); + virtual ~SpectralUnmixingFilterVigra(); + + virtual Eigen::VectorXf NNLSLARS(Eigen::Matrix EndmemberMatrix, + Eigen::VectorXf inputVector); + virtual Eigen::VectorXf NNLSGoldfarb(Eigen::Matrix EndmemberMatrix, + Eigen::VectorXf inputVector); + + private: + + }; + } +} +#endif // MITKPHOTOACOUSTICSPECTRALUNMIXINGFILTERVIGRA_H diff --git a/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterLagrange.cpp b/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterLagrange.cpp new file mode 100644 index 0000000000..bfecfca51f --- /dev/null +++ b/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterLagrange.cpp @@ -0,0 +1,32 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#include "mitkPASpectralUnmixingFilterLagrange.h" + +// Includes for AddEnmemberMatrix +#include "mitkPAPropertyCalculator.h" +#include + + +mitk::pa::SpectralUnmixingFilterLagrange::SpectralUnmixingFilterLagrange() +{ + +} + +mitk::pa::SpectralUnmixingFilterLagrange::~SpectralUnmixingFilterLagrange() +{ + +} diff --git a/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterSimplex.cpp b/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterSimplex.cpp new file mode 100644 index 0000000000..552b11bd35 --- /dev/null +++ b/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterSimplex.cpp @@ -0,0 +1,32 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#include "mitkPASpectralUnmixingFilterSimplex.h" + +// Includes for AddEnmemberMatrix +#include "mitkPAPropertyCalculator.h" +#include + + +mitk::pa::SpectralUnmixingFilterSimplex::SpectralUnmixingFilterSimplex() +{ + +} + +mitk::pa::SpectralUnmixingFilterSimplex::~SpectralUnmixingFilterSimplex() +{ + +} diff --git a/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterVigra.cpp b/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterVigra.cpp new file mode 100644 index 0000000000..77f5b250a3 --- /dev/null +++ b/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterVigra.cpp @@ -0,0 +1,157 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#include "mitkPASpectralUnmixingFilterVigra.h" + +// Includes for AddEnmemberMatrix +#include "mitkPAPropertyCalculator.h" +#include + +// Testing algorithms +#include + +// ImageAccessor +#include +#include + +//vigra +/* +The VIGRA License +================= +(identical to the MIT X11 License) + +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 +#include +#include + +mitk::pa::SpectralUnmixingFilterVigra::SpectralUnmixingFilterVigra() +{ + +} + +mitk::pa::SpectralUnmixingFilterVigra::~SpectralUnmixingFilterVigra() +{ + +} + +Eigen::VectorXf mitk::pa::SpectralUnmixingFilterVigra::NNLSLARS( + Eigen::Matrix EndmemberMatrix, Eigen::VectorXf inputVector) +{ + using namespace vigra; + using namespace vigra::linalg; + + int numberOfWavelengths = EndmemberMatrix.rows(); + int numberOfChromophores = EndmemberMatrix.cols(); + + std::vector A_data; + std::vector B_data; + + for (int i = 0; i < numberOfWavelengths; ++i) + { + B_data.push_back(inputVector(i)); + for (int j = 0; j < numberOfChromophores; ++j) + { + A_data.push_back(EndmemberMatrix(i, j)); + } + } + + const double* Adat = A_data.data(); + const double* Bdat = B_data.data(); + + vigra::Matrix A(Shape2(numberOfWavelengths, numberOfChromophores), Adat); + vigra::Matrix b(Shape2(numberOfWavelengths, 1), Bdat); + vigra::Matrix x(Shape2(numberOfChromophores, 1)); + + // minimize (A*x-b)^2 s.t. x>=0 using least angle regression (LARS algorithm) + nonnegativeLeastSquares(A, b, x); + + Eigen::VectorXf result(numberOfChromophores); + for (int k = 0; k < numberOfChromophores; ++k) + { + float test = x(k, 0); + result[k] = test; + } + return result; +} + +Eigen::VectorXf mitk::pa::SpectralUnmixingFilterVigra::NNLSGoldfarb( + Eigen::Matrix EndmemberMatrix, Eigen::VectorXf inputVector) +{ + using namespace vigra; + using namespace vigra::linalg; + + int numberOfWavelengths = EndmemberMatrix.rows(); + int numberOfChromophores = EndmemberMatrix.cols(); + + std::vector A_data; + std::vector B_data; + + for (int i = 0; i < numberOfWavelengths; ++i) + { + B_data.push_back(inputVector(i)); + for (int j = 0; j < numberOfChromophores; ++j) + { + A_data.push_back(EndmemberMatrix(i, j)); + } + } + + const double* Adat = A_data.data(); + const double* Bdat = B_data.data(); + + vigra::Matrix A(Shape2(numberOfWavelengths, numberOfChromophores), Adat); + vigra::Matrix b(Shape2(numberOfWavelengths, 1), Bdat); + vigra::Matrix x(Shape2(numberOfChromophores, 1)); + + vigra::Matrix eye(identityMatrix(numberOfChromophores)), + zeros(Shape2(numberOfChromophores, 1)), + empty, + U = transpose(A)*A, + v = -transpose(A)*b; + x = 0; + + // minimize (A*x-b)^2 s.t. x>=0 using the Goldfarb-Idnani algorithm + quadraticProgramming(U, v, empty, empty, eye, zeros, x); + + Eigen::VectorXf result(numberOfChromophores); + for (int k = 0; k < numberOfChromophores; ++k) + { + float test = x(k, 0); + result[k] = test; + } + return result; +}