mitkCollectionStatistic.cpp

Go to the documentation of this file.

/*============================================================================

The Medical Imaging Interaction Toolkit (MITK)

Copyright (c) German Cancer Research Center (DKFZ)
All rights reserved.

Use of this source code is governed by a 3-clause BSD license that can be
found in the LICENSE file.

============================================================================*/

#include <mitkCollectionStatistic.h>
#include <mitkDataCollection.h>

// DataCollection Stuff
#include <mitkDataCollectionImageIterator.h>
//stl stuff
#include <sstream>

#include <itkShapedNeighborhoodIterator.h>
#include <itkSignedDanielssonDistanceMapImageFilter.h>


mitk::CollectionStatistic::CollectionStatistic()
  : m_GroundTruthValueToIndexMapper(nullptr)
  , m_TestValueToIndexMapper(nullptr)
{

}

mitk::CollectionStatistic::~CollectionStatistic()
{
  m_GroundTruthValueToIndexMapper = nullptr;
  m_TestValueToIndexMapper = nullptr;
}


void mitk::CollectionStatistic::SetCollection(mitk::DataCollection::Pointer collection)
{
  m_Collection = collection;
}

mitk::DataCollection::Pointer mitk::CollectionStatistic::GetCollection()
{
  return m_Collection;
}

void mitk::CollectionStatistic::SetClassCount (size_t count)
{
  m_ClassCount = count;
}

size_t mitk::CollectionStatistic::GetClassCount()
{
  return m_ClassCount;
}

void mitk::CollectionStatistic::SetGoldName(std::string name)
{
  m_GroundTruthName = name;
}

std::string mitk::CollectionStatistic::GetGoldName()
{
  return m_GroundTruthName;
}

void mitk::CollectionStatistic::SetTestName(std::string name)
{
  m_TestName = name;
}

std::string mitk::CollectionStatistic::GetTestName()
{
  return m_TestName;
}

void mitk::CollectionStatistic::SetGroundTruthValueToIndexMapper(const ValueToIndexMapper* mapper)
{
  m_GroundTruthValueToIndexMapper = mapper;
}

const mitk::ValueToIndexMapper* mitk::CollectionStatistic::GetGroundTruthValueToIndexMapper(void) const
{
  return m_GroundTruthValueToIndexMapper;
}

void mitk::CollectionStatistic::SetTestValueToIndexMapper(const ValueToIndexMapper* mapper)
{
  m_TestValueToIndexMapper = mapper;
}

const mitk::ValueToIndexMapper* mitk::CollectionStatistic::GetTestValueToIndexMapper(void) const
{
  return m_TestValueToIndexMapper;
}

int mitk::CollectionStatistic::IsInSameVirtualClass(unsigned char gold, unsigned char test)
{
  int resultClass = -1;
  for (std::size_t i = 0; i < m_ConnectionGold.size(); ++i)
  {
    if (m_ConnectionGold[i] == gold && m_ConnectionTest[i] == test)
    {
      resultClass = m_ConnectionClass[i];
      break;
    }
  }
  return resultClass;
}

bool mitk::CollectionStatistic::Update()
{
  if (m_GroundTruthValueToIndexMapper == nullptr)
  {
    MITK_ERROR << "m_GroundTruthValueToIndexMapper is nullptr";
    return false;
  }

  if (m_TestValueToIndexMapper == nullptr)
  {
    MITK_ERROR << "m_TestValueToIndexMapper is nullptr";
    return false;
  }

  DataCollectionImageIterator<unsigned char, 3> goldIter(m_Collection, m_GroundTruthName);
  DataCollectionImageIterator<unsigned char, 3> testIter(m_Collection, m_TestName);
  DataCollectionImageIterator<unsigned char, 3> maskIter(m_Collection, m_MaskName);

  int index = 0;

  while (!goldIter.IsAtEnd())
  {
    std::size_t imageIndex = goldIter.GetImageIndex();
    if (m_ImageClassStatistic.size() <= imageIndex)
    {
      MITK_INFO << "New Image: " << goldIter.GetFilePrefix();
      m_ImageNames.push_back(goldIter.GetFilePrefix());
      StatisticData statData;
      m_ImageStatistic.push_back(statData);
      DataVector data;
      for (std::size_t i = 0; i < m_ClassCount; ++i)
      {
        StatisticData stat;
        data.push_back(stat);
      }
      m_ImageClassStatistic.push_back(data);
    }

    if (maskIter.GetVoxel() <= 0)
    {
      ++goldIter;
      ++testIter;
      ++maskIter;
      continue;
    }

    ++index;
    unsigned char goldClass = m_GroundTruthValueToIndexMapper->operator()(goldIter.GetVoxel());
    unsigned char testClass = m_TestValueToIndexMapper->operator()(testIter.GetVoxel());
    if (goldClass == testClass) // True Positive
    {
      m_ImageStatistic[imageIndex].m_TruePositive += 1;
      for (std::size_t i = 0; i < m_ClassCount; ++i)
      {
        if (goldClass == i) // For the detected class it is a true positive
        {
          m_ImageClassStatistic[imageIndex][i].m_TruePositive += 1;
        } else // for all other classes than the detected it is a true negative
        {
          m_ImageClassStatistic[imageIndex][i].m_TrueNegative += 1;
        }
      }
    } else // No True Positive
    {
      m_ImageStatistic[imageIndex].m_FalseNegative += 1;
      m_ImageStatistic[imageIndex].m_FalsePositive += 1;
      for (std::size_t i = 0; i < m_ClassCount; ++i)
      {
        if (goldClass == i) // For the class in Goldstandard it is a false negative
        {
          m_ImageClassStatistic[imageIndex][i].m_FalseNegative += 1;
        } else if ( testClass == i) // For the test class it is a false positive
        {
          m_ImageClassStatistic[imageIndex][i].m_FalsePositive += 1;
        } else // For all other it is a true negative
        {
          m_ImageClassStatistic[imageIndex][i].m_TrueNegative += 1;
        }
      }
    }

    ++goldIter;
    ++testIter;
    ++maskIter;
  }
  MITK_INFO << "Evaluated " << index << " points";
  return true;
}

void mitk::CollectionStatistic::Print(std::ostream& out, std::ostream& sout, bool withHeader, std::string label)
{
  assert(m_ImageClassStatistic.size() > 0);
  assert(m_ImageClassStatistic[0].size() == m_ClassCount);

  if (withHeader)
  {
    sout << "Label;ImageName;";
    for (std::size_t i = 0; i < m_ClassCount; ++i)
    {
      sout << "DICE-Class-"<< i << ";";
      sout << "Jaccard-Class-"<< i << ";";
      sout << "Sensitivity-Class-"<< i << ";";
      sout << "Specificity-Class-"<< i << ";";
      sout << "TP-Class-"<< i << ";";
      sout << "TN-Class-"<< i << ";";
      sout << "FP-Class-"<< i << ";";
      sout << "FN-Class-"<< i << ";";
    }
    sout << "DICE-MEAN"<< ";";
    sout << "Jaccard-MEAN"<< ";";
    sout << "Sensitivity-MEAN"<< ";";
    sout << "Specificity-MEAN"<< ";";
    sout << "TP-MEAN"<< ";";
    sout << "TN-MEAN"<< ";";
    sout << "FP-MEAN"<< ";";
    sout << "FN-MEAN"<< ";";
    sout << "DICE-WMEAN"<< ";";
    sout << "Jaccard-WMEAN"<< ";";
    sout << "Sensitivity-WMEAN"<< ";";
    sout << "Specificity-WMEAN"<< ";";
    sout << "TP-WMEAN"<< ";";
    sout << "TN-WMEAN"<< ";";
    sout << "FP-WMEAN"<< ";";
    sout << "FN-WMEAN"<< ";";
    sout << "COMPLETE-TRUE/FALSE"<< ";";
    sout << "COMPLETE-TRUES"<< ";";
    sout << "COMPLETE_FALSE"<< ";";
    sout << std::endl;
  }
  out << std::setprecision(5);


  MITK_INFO << "m_ImageClassStatistic.size(): " << m_ImageClassStatistic.size();

  for (std::size_t i = 0; i < m_ImageClassStatistic.size(); ++i)
  {
    sout << label << ";"<< m_ImageNames[i]<<";";
    StatisticData meanStat;
    StatisticData wMeanStat;
    double pointsSum = 0;

    out << "======================================================== Image " << std::setw(3) << i << " ========================================================" << std::endl;
    out << " Image ID : " << m_ImageNames[i] <<std::endl;

    out << "|--------|-------------|-------------|-------------|-------------|-------------|-------------|-------------|--------------|" << std::endl;
    out << "| Class  |    DICE     |   Jaccard   | Sensitivity | Specificity |     TP      |     TN      |     FP      |      FN      |" << std::endl;
    out << "|--------|-------------|-------------|-------------|-------------|-------------|-------------|-------------|--------------|" << std::endl;

    for (std::size_t j =0; j < m_ImageClassStatistic[i].size(); ++j)
    {
      StatisticData& stat = m_ImageClassStatistic[i][j];
      stat.m_DICE  = std::max(0.0,(2.0 * stat.m_TruePositive) / (2.0 * stat.m_TruePositive + stat.m_FalseNegative + stat.m_FalsePositive));
      stat.m_Jaccard = std::max(0.0,(1.0 * stat.m_TruePositive) / (1.0 * stat.m_TruePositive + stat.m_FalseNegative + stat.m_FalsePositive));
      stat.m_Sensitivity = std::max(0.0,(1.0 * stat.m_TruePositive) / (stat.m_TruePositive + stat.m_FalseNegative));
      stat.m_Specificity = std::max(0.0,(1.0 * stat.m_TrueNegative)  / ( stat.m_FalsePositive + stat.m_TrueNegative));

      meanStat.m_DICE += std::max(stat.m_DICE,0.0);
      meanStat.m_Jaccard += std::max(stat.m_Jaccard,0.0);
      meanStat.m_Sensitivity += std::max(stat.m_Sensitivity,0.0);
      meanStat.m_Specificity += std::max(stat.m_Specificity,0.0);
      meanStat.m_TruePositive += stat.m_TruePositive;
      meanStat.m_TrueNegative += stat.m_TrueNegative;
      meanStat.m_FalsePositive += stat.m_FalsePositive;
      meanStat.m_FalseNegative += stat.m_FalseNegative;

      double points = stat.m_TruePositive + stat.m_FalseNegative;
      pointsSum += points;

      wMeanStat.m_DICE += std::max(stat.m_DICE,0.0) * points;
      wMeanStat.m_Jaccard += std::max(stat.m_Jaccard,0.0) * points;
      wMeanStat.m_Sensitivity += std::max(stat.m_Sensitivity,0.0) * points;
      wMeanStat.m_Specificity += std::max(stat.m_Specificity,0.0) * points;
      wMeanStat.m_TruePositive += stat.m_TruePositive;
      wMeanStat.m_TrueNegative += stat.m_TrueNegative;
      wMeanStat.m_FalsePositive += stat.m_FalsePositive;
      wMeanStat.m_FalseNegative += stat.m_FalseNegative;

      out<< "|" << std::setw(7) << j << " | ";
      out << std::setw(11) << stat.m_DICE << " | ";
      out << std::setw(11) << stat.m_Jaccard << " | ";
      out << std::setw(11) << stat.m_Sensitivity  << " | ";
      out << std::setw(11) << stat.m_Specificity << " | ";
      out << std::setw(11) << stat.m_TruePositive << " | ";
      out << std::setw(11) << stat.m_TrueNegative<< " | ";
      out << std::setw(11) << stat.m_FalsePositive<< " | ";
      out << std::setw(11) << stat.m_FalseNegative << "  |"<< std::endl;
      sout << stat.m_DICE << ";";
      sout << stat.m_Jaccard << ";";
      sout << stat.m_Sensitivity  << ";";
      sout << stat.m_Specificity << ";";
      sout << stat.m_TruePositive << ";";
      sout << stat.m_TrueNegative<< ";";
      sout << stat.m_FalsePositive<< ";";
      sout << stat.m_FalseNegative << ";";
      sout << stat.m_RMSD << ";";
    }

    meanStat.m_DICE /=  m_ImageClassStatistic[i].size();
    meanStat.m_Jaccard /=  m_ImageClassStatistic[i].size();
    meanStat.m_Sensitivity /=  m_ImageClassStatistic[i].size();
    meanStat.m_Specificity /=  m_ImageClassStatistic[i].size();
    meanStat.m_TruePositive /=  m_ImageClassStatistic[i].size();
    meanStat.m_TrueNegative /=  m_ImageClassStatistic[i].size();
    meanStat.m_FalsePositive /= m_ImageClassStatistic[i].size();
    meanStat.m_FalseNegative /=  m_ImageClassStatistic[i].size();

    out << "|--------|-------------|-------------|-------------|-------------|-------------|-------------|-------------|--------------|" << std::endl;
    out << std::setw(7) << "| Mean "<< "  | ";
    out << std::setw(11) << meanStat.m_DICE << " | ";
    out << std::setw(11) << meanStat.m_Jaccard << " | ";
    out << std::setw(11) << meanStat.m_Sensitivity  << " | ";
    out << std::setw(11) << meanStat.m_Specificity << " | ";
    out << std::setw(11) << meanStat.m_TruePositive << " | ";
    out << std::setw(11) << meanStat.m_TrueNegative<< " | ";
    out << std::setw(11) << meanStat.m_FalsePositive<< " | ";
    out << std::setw(11) << meanStat.m_FalseNegative << "  |"<< std::endl;
    sout << meanStat.m_DICE << ";";
    sout << meanStat.m_Jaccard << ";";
    sout << meanStat.m_Sensitivity  << ";";
    sout << meanStat.m_Specificity << ";";
    sout << meanStat.m_TruePositive << ";";
    sout << meanStat.m_TrueNegative<< ";";
    sout << meanStat.m_FalsePositive<< ";";
    sout << meanStat.m_FalseNegative << ";";

    wMeanStat.m_DICE /=  pointsSum;
    wMeanStat.m_Jaccard /=  pointsSum;
    wMeanStat.m_Sensitivity /=  pointsSum;
    wMeanStat.m_Specificity /=  pointsSum;
    wMeanStat.m_TruePositive /=  pointsSum;
    wMeanStat.m_TrueNegative /= pointsSum;
    wMeanStat.m_FalsePositive /= pointsSum;
    wMeanStat.m_FalseNegative /=  pointsSum;

    out << std::setw(7) << "| W-Mean"<< " | ";
    out << std::setw(11) << wMeanStat.m_DICE << " | ";
    out << std::setw(11) << wMeanStat.m_Jaccard << " | ";
    out << std::setw(11) << wMeanStat.m_Sensitivity  << " | ";
    out << std::setw(11) << wMeanStat.m_Specificity << " | ";
    out << std::setw(11) << wMeanStat.m_TruePositive << " | ";
    out << std::setw(11) << wMeanStat.m_TrueNegative<< " | ";
    out << std::setw(11) << wMeanStat.m_FalsePositive<< " | ";
    out << std::setw(11) << wMeanStat.m_FalseNegative << "  |"<< std::endl;
    sout << wMeanStat.m_DICE << ";";
    sout << wMeanStat.m_Jaccard << ";";
    sout << wMeanStat.m_Sensitivity  << ";";
    sout << wMeanStat.m_Specificity << ";";
    sout << wMeanStat.m_TruePositive << ";";
    sout << wMeanStat.m_TrueNegative<< ";";
    sout << wMeanStat.m_FalsePositive<< ";";
    sout << wMeanStat.m_FalseNegative << ";";

    m_ImageStatistic[i].m_Sensitivity = std::max(0.0,(1.0 * m_ImageStatistic[i].m_TruePositive) / (m_ImageStatistic[i].m_TruePositive + m_ImageStatistic[i].m_FalseNegative));
    m_ImageStatistic[i].m_Specificity = std::max(0.0,(1.0 * m_ImageStatistic[i].m_TrueNegative)  / ( m_ImageStatistic[i].m_FalsePositive + m_ImageStatistic[i].m_TrueNegative));

    out << std::setw(7) << "| Compl."<< " | ";
    out << std::setw(11) << "     x     " << " | ";
    out << std::setw(11) << "     x     " << " | ";
    out << std::setw(11) << m_ImageStatistic[i].m_Sensitivity  << " | ";
    out << std::setw(11) << "     x     " << " | ";
    out << std::setw(11) << m_ImageStatistic[i].m_TruePositive << " | ";
    out << std::setw(11) << m_ImageStatistic[i].m_TrueNegative<< " | ";
    out << std::setw(11) << m_ImageStatistic[i].m_FalsePositive<< " | ";
    out << std::setw(11) << m_ImageStatistic[i].m_FalseNegative << "  |"<< std::endl;
    out << "|--------|-------------|-------------|-------------|-------------|-------------|-------------|-------------|--------------|" << std::endl;

    sout << m_ImageStatistic[i].m_Sensitivity  << ";";
    sout << m_ImageStatistic[i].m_TruePositive << ";";
    sout << m_ImageStatistic[i].m_FalsePositive<< std::endl;
    out << std::endl;
  }
}

std::vector<mitk::StatisticData> mitk::CollectionStatistic::GetStatisticData(unsigned char c) const
{
  std::vector<StatisticData> statistics;

  for (size_t i = 0; i < m_ImageClassStatistic.size(); i++)
  {
    statistics.push_back(m_ImageClassStatistic[i][c]);
  }

  return statistics;
}

void mitk::CollectionStatistic::ComputeRMSD()
{
  assert(m_ClassCount == 2);
  assert(m_GroundTruthValueToIndexMapper != nullptr);
  assert(m_TestValueToIndexMapper != nullptr);

  DataCollectionImageIterator<unsigned char, 3> groundTruthIter(m_Collection, m_GroundTruthName);
  DataCollectionImageIterator<unsigned char, 3> testIter(m_Collection, m_TestName);
  DataCollectionImageIterator<unsigned char, 3> maskIter(m_Collection, m_MaskName);

  typedef itk::Image<unsigned char, 3> LabelImage;
  typedef itk::Image<double, 3> ImageType;
  typedef itk::SignedDanielssonDistanceMapImageFilter<LabelImage, ImageType, ImageType> DistanceMapFilterType;
  typedef itk::ConstantBoundaryCondition<LabelImage> BoundaryConditionType;
  typedef itk::ConstShapedNeighborhoodIterator<LabelImage, BoundaryConditionType> ConstNeighborhoodIteratorType;

  // Build up 6-neighborhood. Diagonal voxel are maybe not considered for distance map computation.
  // So 6-neighborhood avoids inconsistencies.
  ConstNeighborhoodIteratorType::OffsetType offset0 = {{ 0,  0, -1}};
  ConstNeighborhoodIteratorType::OffsetType offset1 = {{ 0,  0,  1}};
  ConstNeighborhoodIteratorType::OffsetType offset2 = {{ 0, -1,  0}};
  ConstNeighborhoodIteratorType::OffsetType offset3 = {{ 0,  1,  0}};
  ConstNeighborhoodIteratorType::OffsetType offset4 = {{-1,  0,  0}};
  ConstNeighborhoodIteratorType::OffsetType offset5 = {{ 1,  0,  0}};

  const int outsideVal = 17;
  itk::NeighborhoodIterator<LabelImage>::RadiusType radius;
  radius.Fill(1);

  BoundaryConditionType bc;
  // if a neighborhood voxel is outside of the image region a default value is returned
  bc.SetConstant(outsideVal);

  ConstNeighborhoodIteratorType neighborhoodIter;

  ImageType::Pointer distanceImage;
  std::vector<mitk::StatisticData>* currentImageStatistics = nullptr;

  unsigned int distanceBorderSamples = 0;
  double totalBorderRMSDistance = 0;

  int previousImageIndex = -1;

  while (!testIter.IsAtEnd())
  {
    int currentImageIndex = testIter.GetImageIndex();

    // prepare data for next image
    if (previousImageIndex != currentImageIndex)
    {
      previousImageIndex = currentImageIndex;

      currentImageStatistics = &(m_ImageClassStatistic.at(currentImageIndex));

      distanceBorderSamples = 0;
      totalBorderRMSDistance = 0;

      // generate distance map from gold standard image
      DistanceMapFilterType::Pointer distanceMapFilter = DistanceMapFilterType::New();
      distanceMapFilter->SetInput(groundTruthIter.GetImageIterator().GetImage());
      distanceMapFilter->SetUseImageSpacing(true);
      distanceMapFilter->Update();
      distanceImage = distanceMapFilter->GetOutput();

      neighborhoodIter = ConstNeighborhoodIteratorType(radius, testIter.GetImageIterator().GetImage(), testIter.GetImageIterator().GetImage()->GetRequestedRegion());
      neighborhoodIter.OverrideBoundaryCondition(&bc);

      // configure 6-neighborhood
      neighborhoodIter.ActivateOffset(offset0);
      neighborhoodIter.ActivateOffset(offset1);
      neighborhoodIter.ActivateOffset(offset2);
      neighborhoodIter.ActivateOffset(offset3);
      neighborhoodIter.ActivateOffset(offset4);
      neighborhoodIter.ActivateOffset(offset5);
    }

    unsigned char testClass = m_TestValueToIndexMapper->operator()(testIter.GetVoxel());

    if ( maskIter.GetVoxel() > 0 && testClass != 0)
    {
      // check if it is a border voxel
      neighborhoodIter.SetLocation(testIter.GetImageIterator().GetIndex());
      bool border = false;

      ConstNeighborhoodIteratorType::ConstIterator iter;
      for (iter = neighborhoodIter.Begin(); !iter.IsAtEnd(); iter++)
      {
        if (iter.Get() == outsideVal)
        {
          continue;
        }

        if (m_TestValueToIndexMapper->operator()(iter.Get()) != 1 )
        {
          border = true;
          break;
        }
      }

      if (border)
      {
        double currentDistance = distanceImage->GetPixel(testIter.GetImageIterator().GetIndex());
        totalBorderRMSDistance += currentDistance * currentDistance;
        ++distanceBorderSamples;

        // update immediately, so the final value is set after the iterator of the last image has reached the end
        double rmsd = std::sqrt(totalBorderRMSDistance / (double) distanceBorderSamples);
        currentImageStatistics->at(1).m_RMSD = rmsd;
      }
    }
    ++groundTruthIter;
    ++testIter;
    ++maskIter;
  }
}