Medical Imaging Interaction Toolkit  2018.4.99-1bab67a2
Medical Imaging Interaction Toolkit
mitkPASpectralUnmixingFilterBase.cpp
Go to the documentation of this file.
1 /*============================================================================
2 
3 The Medical Imaging Interaction Toolkit (MITK)
4 
5 Copyright (c) German Cancer Research Center (DKFZ)
6 All rights reserved.
7 
8 Use of this source code is governed by a 3-clause BSD license that can be
9 found in the LICENSE file.
10 
11 ============================================================================*/
12 
14 
15 // Includes for AddEnmemberMatrix
17 #include <eigen3/Eigen/Dense>
18 
19 // ImageAccessor
20 #include <mitkImageReadAccessor.h>
21 #include <mitkImageWriteAccessor.h>
22 
24 {
25  m_PropertyCalculatorEigen = mitk::pa::PropertyCalculator::New();
26 }
27 
29 {
30 
31 }
32 
34 {
35  if (outputs == 0)
36  mitkThrow() << "ADD OUTPUTS HAS TO BE LARGER THEN ZERO!";
37  this->SetNumberOfIndexedOutputs(outputs);
38  for (unsigned int i = 0; i<GetNumberOfIndexedOutputs(); i++)
39  this->SetNthOutput(i, mitk::Image::New());
40 }
41 
43 {
44  m_Wavelength.push_back(wavelength);
45 }
46 
48 {
49  m_Chromophore.push_back(chromophore);
50 }
51 
53 {
55 }
56 
58 {
59  m_RelativeError = relativeError;
60 }
61 
63 {
64  m_RelativeErrorSettings.push_back(value);
65 }
66 
67 void mitk::pa::SpectralUnmixingFilterBase::GenerateData()
68 {
69  MITK_INFO(m_Verbose) << "GENERATING DATA..";
70 
71  mitk::Image::Pointer input = GetInput(0);
72  CheckPreConditions(input);
73 
74  unsigned int xDim = input->GetDimensions()[0];
75  unsigned int yDim = input->GetDimensions()[1];
76  unsigned int numberOfInputImages = input->GetDimensions()[2];
77 
78  MITK_INFO(m_Verbose) << "x dimension: " << xDim;
79  MITK_INFO(m_Verbose) << "y dimension: " << yDim;
80  MITK_INFO(m_Verbose) << "z dimension: " << numberOfInputImages;
81 
82  // Copy input image into array
83  mitk::ImageReadAccessor readAccess(input);
84  const float* inputDataArray = ((const float*)readAccess.GetData());
85 
86  unsigned int sequenceSize = m_Wavelength.size();
87  unsigned int totalNumberOfSequences = numberOfInputImages / sequenceSize;
88  MITK_INFO(m_Verbose) << "TotalNumberOfSequences: " << totalNumberOfSequences;
89 
90  InitializeOutputs(totalNumberOfSequences);
91 
92  auto endmemberMatrix = CalculateEndmemberMatrix(m_Chromophore, m_Wavelength);
93 
94  // test to see pixel values @ txt file
95  myfile.open("SimplexNormalisation.txt");
96 
97  unsigned int outputCounter = GetNumberOfIndexedOutputs();
98  std::vector<float*> writteBufferVector;
99  for (unsigned int i = 0; i < outputCounter; ++i)
100  {
101  auto output = GetOutput(i);
102  mitk::ImageWriteAccessor writeOutput(output);
103  float* writeBuffer = (float *)writeOutput.GetData();
104  writteBufferVector.push_back(writeBuffer);
105  }
106 
107  if (m_RelativeError == true)
108  {
109  // -1 because rel error is output[IndexedOutputs() - 1] and loop over chromophore outputs has to end at [IndexedOutputs() - 2]
110  outputCounter -= 1;
111  }
112 
113  for (unsigned int sequenceCounter = 0; sequenceCounter < totalNumberOfSequences; ++sequenceCounter)
114  {
115  MITK_INFO(m_Verbose) << "SequenceCounter: " << sequenceCounter;
116  //loop over every pixel in XY-plane
117  for (unsigned int x = 0; x < xDim; x++)
118  {
119  for (unsigned int y = 0; y < yDim; y++)
120  {
121  Eigen::VectorXf inputVector(sequenceSize);
122  for (unsigned int z = 0; z < sequenceSize; z++)
123  {
129  unsigned int pixelNumber = (xDim*yDim*(z+sequenceCounter*sequenceSize)) + x * yDim + y;
130  auto pixel = inputDataArray[pixelNumber];
131 
132  inputVector[z] = pixel;
133  }
134  Eigen::VectorXf resultVector = SpectralUnmixingAlgorithm(endmemberMatrix, inputVector);
135 
136  if (m_RelativeError == true)
137  {
138  float relativeError = CalculateRelativeError(endmemberMatrix, inputVector, resultVector);
139  writteBufferVector[outputCounter][(xDim*yDim * sequenceCounter) + x * yDim + y] = relativeError;
140  }
141 
142  for (unsigned int outputIdx = 0; outputIdx < outputCounter; ++outputIdx)
143  {
144  writteBufferVector[outputIdx][(xDim*yDim * sequenceCounter) + x * yDim + y] = resultVector[outputIdx];
145  }
146  }
147  }
148  }
149  MITK_INFO(m_Verbose) << "GENERATING DATA...[DONE]";
150  myfile.close();
151 }
152 
153 void mitk::pa::SpectralUnmixingFilterBase::CheckPreConditions(mitk::Image::Pointer input)
154 {
155  MITK_INFO(m_Verbose) << "CHECK PRECONDITIONS ...";
156 
157  if (m_Chromophore.size() == 0 || m_Wavelength.size() == 0)
158  mitkThrow() << "NO WAVELENGHTS/CHROMOPHORES SELECTED!";
159 
160  if (m_Wavelength.size() < input->GetDimensions()[2])
161  MITK_WARN(m_Verbose) << "NUMBER OF WAVELENGTHS < NUMBER OF INPUT IMAGES";
162 
163  if (m_Chromophore.size() > m_Wavelength.size())
164  mitkThrow() << "ADD MORE WAVELENGTHS OR REMOVE ENDMEMBERS!";
165 
166  if (input->GetPixelType() != mitk::MakeScalarPixelType<float>())
167  mitkThrow() << "PIXELTYPE ERROR! FLOAT REQUIRED";
168 
169  if ((m_Chromophore.size()+ m_RelativeError )!= GetNumberOfIndexedOutputs() || input->GetDimensions()[2] < GetNumberOfIndexedOutputs())
170  mitkThrow() << "INDEX ERROR! NUMBER OF OUTPUTS DOESN'T FIT TO OTHER SETTIGNS!";
171 
172  MITK_INFO(m_Verbose) << "...[DONE]";
173 }
174 
175 void mitk::pa::SpectralUnmixingFilterBase::InitializeOutputs(unsigned int totalNumberOfSequences)
176 {
177  MITK_INFO(m_Verbose) << "Initialize Outputs ...";
178 
179  unsigned int numberOfInputs = GetNumberOfIndexedInputs();
180  unsigned int numberOfOutputs = GetNumberOfIndexedOutputs();
181  MITK_INFO(m_Verbose) << "Inputs: " << numberOfInputs << " Outputs: " << numberOfOutputs;
182 
183  mitk::PixelType pixelType = mitk::MakeScalarPixelType<float>();
184  const int NUMBER_OF_SPATIAL_DIMENSIONS = 3;
185  auto* dimensions = new unsigned int[NUMBER_OF_SPATIAL_DIMENSIONS];
186  for (unsigned int dimIdx = 0; dimIdx < 2; dimIdx++)
187  dimensions[dimIdx] = GetInput()->GetDimensions()[dimIdx];
188  dimensions[2] = totalNumberOfSequences;
189 
190  for (unsigned int outputIdx = 0; outputIdx < numberOfOutputs; outputIdx++)
191  GetOutput(outputIdx)->Initialize(pixelType, NUMBER_OF_SPATIAL_DIMENSIONS, dimensions);
192  MITK_INFO(m_Verbose) << "...[DONE]";
193 }
194 
195 Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic> mitk::pa::SpectralUnmixingFilterBase::CalculateEndmemberMatrix(
196  std::vector<mitk::pa::PropertyCalculator::ChromophoreType> m_Chromophore, std::vector<int> m_Wavelength)
197 {
198  unsigned int numberOfChromophores = m_Chromophore.size(); //columns
199  unsigned int numberOfWavelengths = m_Wavelength.size(); //rows
200  Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic> endmemberMatrixEigen(numberOfWavelengths, numberOfChromophores);
201 
202  for (unsigned int j = 0; j < numberOfChromophores; ++j)
203  {
204  for (unsigned int i = 0; i < numberOfWavelengths; ++i)
205  endmemberMatrixEigen(i, j) = PropertyElement(m_Chromophore[j], m_Wavelength[i]);
206  }
207  MITK_INFO(m_Verbose) << "GENERATING ENMEMBERMATRIX [DONE]";
208  return endmemberMatrixEigen;
209 }
210 
211 float mitk::pa::SpectralUnmixingFilterBase::PropertyElement(mitk::pa::PropertyCalculator::ChromophoreType chromophore, int wavelength)
212 {
213  if (chromophore == mitk::pa::PropertyCalculator::ChromophoreType::ONEENDMEMBER)
214  return 1;
215  else
216  {
217  float value = m_PropertyCalculatorEigen->GetAbsorptionForWavelength(chromophore, wavelength);
218  if (value == 0)
219  mitkThrow() << "WAVELENGTH " << wavelength << "nm NOT SUPPORTED!";
220  else
221  return value;
222  }
223 }
224 
225 float mitk::pa::SpectralUnmixingFilterBase::CalculateRelativeError(Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic> endmemberMatrix,
226  Eigen::VectorXf inputVector, Eigen::VectorXf resultVector)
227 {
228  float relativeError = (endmemberMatrix*resultVector - inputVector).norm() / inputVector.norm();
229  for (int i = 0; i < 2; ++i)
230  {
231  if (resultVector[i] < m_RelativeErrorSettings[i])
232  return 0;
233  }
234  return relativeError;
235 }
#define MITK_INFO
Definition: mitkLogMacros.h:18
void AddWavelength(int wavelength)
AddWavelength takes integers and writes them at the end of the m_Wavelength vector. The first call of the method then corresponds to the first input image and so on.
unsigned int * GetDimensions() const
Get the sizes of all dimensions as an integer-array.
Definition: mitkImage.cpp:1309
void * GetData()
Gives full data access.
virtual void AddOutputs(unsigned int outputs)
AddOutputs takes an integer and sets indexed outputs.
virtual void AddRelativeErrorSettings(int value)
AddRelativeErrorSettings takes integers and writes them at the end of the m_RelativeErrorSettings vec...
std::vector< mitk::pa::PropertyCalculator::ChromophoreType > m_Chromophore
#define MITK_WARN
Definition: mitkLogMacros.h:19
void AddChromophore(mitk::pa::PropertyCalculator::ChromophoreType chromophore)
AddChromophore takes mitk::pa::PropertyCalculator::ChromophoreType and writes them at the end of the ...
#define mitkThrow()
bool verbose(false)
static Pointer New()
InputImageType * GetInput(void)
SpectralUnmixingFilterBase()
Constructor creats proptery calculater smart pointer new()
OutputType * GetOutput()
Get the output data of this image source object.
ImageWriteAccessor class to get locked write-access for a particular image part.
virtual Eigen::VectorXf SpectralUnmixingAlgorithm(Eigen::Matrix< float, Eigen::Dynamic, Eigen::Dynamic > endmemberMatrix, Eigen::VectorXf inputVector)=0
The subclasses will override the mehtod to calculate the spectral unmixing result vector...
ImageReadAccessor class to get locked read access for a particular image part.
const void * GetData() const
Gives const access to the data.
Class for defining the data type of pixels.
Definition: mitkPixelType.h:51