mlClusterAlgorithm.h

Go to the documentation of this file.

/*************************************************************************************
**
** Copyright 2007, MeVis Medical Solutions AG
**
** The user may use this file in accordance with the license agreement provided with
** the Software or, alternatively, in accordance with the terms contained in a
** written agreement between the user and MeVis Medical Solutions AG.
**
** For further information use the contact form at https://www.mevislab.de/contact
**
**************************************************************************************/
 
#ifndef ML_CLUSTER_ALGORITHM_H
#define ML_CLUSTER_ALGORITHM_H
 
#include "MLClusterAlgorithmSystem.h"
 
#include "mlClusters.h"
#include "mlClusterInfo.h"
#include "mlClusterRefVolume.h"
#include "mlClusterRefCollection.h"
 
#include <mlModuleIncludes.h>
#include <mlProcessAllPagesHandler.h>
#include <mlHost.h>
 
ML_START_NAMESPACE
 
struct ComputeClusterParameters
{
  ComputeClusterParameters()
  {
    clusterMode = CLUSTER_MODE_IdenticalIntensities;
    similarityToleranceIntensities = 0;
    similarityToleranceVectors = 0.98;
    neighborhoodRelation = NBH_3D_26_XYZ;
    contentImage = nullptr;
    maskImage = nullptr;
    useBackgroundValue = true;
    backgroundValue = 0;
    useImageValueAsUserData = false;
    progressCB = nullptr;
    progressCBUserData = nullptr;
  }
 
  CLUSTER_MODE_TYPE clusterMode;
  MLdouble similarityToleranceIntensities;
  MLdouble similarityToleranceVectors;
  NBH_TYPE neighborhoodRelation;
  PagedImage* contentImage;
  PagedImage* maskImage;
  bool useBackgroundValue;
  MLdouble backgroundValue;
  bool useImageValueAsUserData;
  MLRequestProgressCB* progressCB;
  void* progressCBUserData;
};
 
 
template <typename CLUSTERVOXELTYPE, typename DerivedAlgorithm> class ClusterAlgorithm;
 
template <typename CLUSTERVOXELTYPE, typename DerivedAlgorithm>
class ClusterHandler : public OrderedProcessAllPagesHandler
{
  ML_DISALLOW_COPY_AND_ASSIGN(ClusterHandler)
 
public:
  ClusterHandler(ClusterAlgorithm<CLUSTERVOXELTYPE, DerivedAlgorithm>* algorithm)
  {
    _algorithm = algorithm;
    _error = ML_RESULT_OK;
    _useMask = false;
    _dataType = -1;
  }
 
  void setParameters(bool useMask, MLDataType type)
  {
    _useMask = useMask;
    _dataType = type;
  }
 
  void calculateOutputImageProperties(PagedImage* image) override
  {
    OrderedProcessAllPagesHandler::calculateOutputImageProperties(image);
    image->setDataType(_dataType);
    if (_useMask)
    {
      image->setInputSubImageDataType(1, MLint8Type);
    }
    image->setInputSubImagesAreReadOnly(true);
    ImageVector imageExt = image->getImageExtent();
    imageExt.c = 1;
    imageExt.t = 1;
    imageExt.u = 1;
    image->setImageExtent(imageExt);
    ImageVector pageExt = imageExt;
    pageExt.z = 1;
    image->setPageExtent(pageExt);
    _inImageExtent = imageExt;
  }
 
  SubImageBox calculateInputSubImageBox(int inputIndex, const SubImageBox& outputSubImageBox) override
  {
    if ((inputIndex == 0) || (inputIndex == 1 && _useMask))
    {
      return outputSubImageBox;
    }
    else
    {
      return SubImageBox();
    }
  }
 
  MLErrorCode processTiles(SubImage* images) override
  {
    TSubImage<MLint8> maskImage;
    if (_useMask)
    {
      maskImage = TSubImage<MLint8>(images[1]);
    }
    _algorithm->storeNextSlice(images[0]);
    _error = _algorithm->processSlice(images[0].getOrigin().z, _inImageExtent, maskImage);
    return _error;
  }
 
  MLErrorCode getErrorCode() const
  {
    return _error;
  }
 
private:
  ClusterAlgorithm<CLUSTERVOXELTYPE, DerivedAlgorithm>* _algorithm;
 
  ImageVector _inImageExtent;
  MLDataType _dataType;
  bool _useMask;
  MLErrorCode _error;
};
 
//-----------------------------------------------------------------------------------------------
//-----------------------------------------------------------------------------------------------
class ClusterAlgorithmBase
{
public:
 
  ClusterAlgorithmBase(const ComputeClusterParameters& parameters, Clusters* clusters);
 
 
  virtual ~ClusterAlgorithmBase();
 
 
  bool hasResults() const { return _hasResults; }
 
 
protected:
  enum { _NR_OF_SUBIMAGES = 2 } NumInImagesEnum;
 
  void setInSliceVoxelType(MLDataType dt);
 
  void freeInSliceData();
 
  void storeNextSlice(SubImage& slice);
 
  MLErrorCode loadNextSlice(MLint sliceNr, const ImageVector& inImgExt);
 
  MLErrorCode loadMaskSlice(TSubImage<MLint8>& slice, MLint sliceNr, const ImageVector& inImgExt) const;
 
  MLuint mergeClusterReference(MLint x, MLint y, MLint z, MLuint currClusterRef);
 
  void setResultFlag(bool result) { _hasResults = result; }
 
  Clusters* _clusters;
 
  SubImage _currentSlices[_NR_OF_SUBIMAGES];
 
  ComputeClusterParameters _parameters;
 
private:
  bool _hasResults;
 
  ClusterAlgorithmBase& operator=(const ClusterAlgorithmBase&);
 
 
  template <typename CLUSTERVOXELTYPE, typename DerivedAlgorithm>
  friend class ClusterHandler;
};
 
 
//-----------------------------------------------------------------------------------------------
//-----------------------------------------------------------------------------------------------
template <typename CLUSTERVOXELTYPE, typename DerivedAlgorithm>
class ClusterAlgorithm : public ClusterAlgorithmBase
{
public:
 
  typedef CLUSTERVOXELTYPE ClusterVoxelType;
 
  ClusterAlgorithm(const ComputeClusterParameters& parameters, Clusters* clusters) : ClusterAlgorithmBase(parameters, clusters)
  {
  }
 
  MLErrorCode runClustering()
  {
    return static_cast<DerivedAlgorithm*>(this)->calcClustersForSlices();
  }
 
protected:
 
 
  inline MLuint updateClusterRefForNeighbor(void* neighborVoxelPtr, MLint x, MLint y, MLint z,
                                            CLUSTERVOXELTYPE currentValue,
                                            MLuint currentClusterRef)
  {
    // Get predecessor voxel value.
    CLUSTERVOXELTYPE neighborVoxel = DerivedAlgorithm::getVoxel(neighborVoxelPtr);
    const bool isInToleranceFlag = static_cast<DerivedAlgorithm*>(this)->isInTolerance(neighborVoxel, currentValue);
    if (isInToleranceFlag)
    {
      currentClusterRef = mergeClusterReference(x, y, z, currentClusterRef);
    }
    return currentClusterRef;
  }
 
  inline static CLUSTERVOXELTYPE convertBackgroundValue(MLdouble backgroundValue)
  {
    return static_cast<CLUSTERVOXELTYPE>(backgroundValue);
  }
 
  MLErrorCode calcClustersForSlices()
  {
    // Allocate memory with subimages of the requested voxel type.
    setInSliceVoxelType(_parameters.contentImage->getDataType());
 
    ImageVector inputImageExtent = _parameters.contentImage->getImageExtent();
 
    TSubImage<MLint8> inMaskSlice;
 
    // Initialize data structures for clustering.
    MLErrorCode errorCode = _clusters->init(inputImageExtent.x, inputImageExtent.y, inputImageExtent.z);
 
    if (ML_RESULT_OK == errorCode)
    {
      if ((_parameters.neighborhoodRelation == NBH_2D_4_XY) || (_parameters.neighborhoodRelation == NBH_2D_8_XY))
      {
        _clusters->setUse2DNeighborhood();
      }
 
      const bool useMask = _parameters.maskImage != nullptr;
 
      if (!Host::getDefaultHost().getUseClassicHost())
      {
        ClusterHandler<CLUSTERVOXELTYPE, DerivedAlgorithm> handler(this);
        handler.setParameters(useMask, _currentSlices[0].getDataType());
 
        std::vector<PagedImage*> images;
        images.push_back(_parameters.contentImage);
        if (useMask)
        {
          images.push_back(_parameters.maskImage);
        }
 
        errorCode = Host::getDefaultHost().processAllPagesWithInputImages(images, handler, SubImageBox(), _parameters.progressCB, _parameters.progressCBUserData);
        if (errorCode != ML_RESULT_OK)
        {
          errorCode = handler.getErrorCode();
        }
      }
      else
      {
        // Main loop, proceed slice-wise only on success.
        for (MLint zi = 0; (ML_RESULT_OK == errorCode) && (zi < inputImageExtent.z); zi++)
        {
          errorCode = loadNextSlice(zi, inputImageExtent);
          if (ML_RESULT_OK == errorCode && useMask)
          {
            errorCode = loadMaskSlice(inMaskSlice, zi, inputImageExtent);
          }
 
          if (ML_RESULT_OK == errorCode)
          {
            errorCode = processSlice(zi, inputImageExtent, inMaskSlice);
          }
        }
      }
 
      if (ML_RESULT_OK == errorCode)
      {
        ClusterUserDataParameters userDataParameters;
        userDataParameters.useImageValueAsUserData = _parameters.useImageValueAsUserData;
        userDataParameters.voxelSizeInMilliliters = _clusters->getVoxelVolumeInMilliliters();
        userDataParameters.backgroundValue = _parameters.backgroundValue;
 
        _clusters->calcClusterSizes(userDataParameters);
      }
    }
 
    setResultFlag(ML_RESULT_OK == errorCode);
    freeInSliceData();
 
    return errorCode;
  }
 
public:
 
  MLErrorCode processSlice(MLint zi, const ImageVector& inputImageExtent, TSubImage<MLint8>& inMaskSlice)
  {
    MLErrorCode errorCode = ML_RESULT_OK;
 
    const MLint xres = inputImageExtent.x;
    const MLint yres = inputImageExtent.y;
 
    const ClusterVoxelType backgroundValue = DerivedAlgorithm::convertBackgroundValue(_parameters.backgroundValue);
 
    MLTypeData* sliceInternalPtr = reinterpret_cast<MLTypeData*>(_currentSlices[1].getData());
    MLTypeData* previousSliceInternalPtr = reinterpret_cast<MLTypeData*>(_currentSlices[0].getData());
    const MLDataType dataType = _currentSlices[1].getDataType();
    size_t dataTypeSize = MLSizeOf(dataType);
 
    MLint8* maskPtr = inMaskSlice.getData();
 
    size_t xoffset = dataTypeSize;
    size_t yoffset = dataTypeSize * xres;
 
    const bool useMask = _parameters.maskImage != nullptr;
    const bool useBackgroundValue = _parameters.useBackgroundValue;
 
    // iterate over slice voxels
    for (MLint yi = 0; yi < yres; yi++)
    {
      for (MLint xi = 0; xi < xres; xi++)
      {
        MLTypeData* slicePtr = sliceInternalPtr;
        MLTypeData* previousSlicePtr = previousSliceInternalPtr;
        // advance the internal pointers, we can't advance slicePtr/previousSlicePtr,
        // because they are accessed below and we don't want to place the advancing into each
        // continue statement.
        sliceInternalPtr += dataTypeSize;
        previousSliceInternalPtr += dataTypeSize;
 
        // get current voxel value from untyped image.
        CLUSTERVOXELTYPE currentValue = DerivedAlgorithm::getVoxel(slicePtr);
 
        // get corresponding mask value
        if (useMask)
        {
          MLint8 maskValue = *maskPtr++;
          if (maskValue == 0)
          {
            continue;
          }
        }
 
        // if background skip voxel
        if (useBackgroundValue && (currentValue == backgroundValue))
        {
          continue;
        }
 
        // reset cluster reference
        MLuint currentClusterRef = 0;
 
        // left and top neighbor are always evaluated
        if (xi > 0)
        {
          currentClusterRef = updateClusterRefForNeighbor(slicePtr - xoffset, xi - 1, yi, zi, currentValue, currentClusterRef);
        }
        if (yi > 0)
        {
          currentClusterRef = updateClusterRefForNeighbor(slicePtr - yoffset, xi, yi - 1, zi, currentValue, currentClusterRef);
        }
 
        switch (_parameters.neighborhoodRelation)
        {
        case NBH_2D_8_XY:
          if ((xi > 0) && (yi > 0))
          {
            currentClusterRef = updateClusterRefForNeighbor(slicePtr - xoffset - yoffset, xi - 1, yi - 1, zi, currentValue, currentClusterRef);
          }
          if ((xi<xres - 1) && (yi>0))
          {
            currentClusterRef = updateClusterRefForNeighbor(slicePtr + xoffset - yoffset, xi + 1, yi - 1, zi, currentValue, currentClusterRef);
          }
          break;
          // The missing breaks in the following 3d cases have been left out intentionally to be able to deal with
          // the different 3d neighborhoods without doubling code lines (so complete 26-3d-neighborhood calculation is
          // done by case NBH_3D_26_XYZ +  NBH_3D_18_XYZ + NBH_3D_6_XYZ)
        case NBH_3D_26_XYZ:
          if ((xi > 0) && (yi > 0) && (zi > 0))
          {
            currentClusterRef = updateClusterRefForNeighbor(previousSlicePtr - xoffset - yoffset, xi - 1, yi - 1, zi - 1, currentValue, currentClusterRef);
          }
          if ((xi<xres - 1) && (yi>0) && (zi > 0))
          {
            currentClusterRef = updateClusterRefForNeighbor(previousSlicePtr + xoffset - yoffset, xi + 1, yi - 1, zi - 1, currentValue, currentClusterRef);
          }
          if ((xi > 0) && (yi<yres - 1) && (zi>0))
          {
            currentClusterRef = updateClusterRefForNeighbor(previousSlicePtr - xoffset + yoffset, xi - 1, yi + 1, zi - 1, currentValue, currentClusterRef);
          }
          if ((xi<xres - 1) && (yi<yres - 1) && (zi>0))
          {
            currentClusterRef = updateClusterRefForNeighbor(previousSlicePtr + xoffset + yoffset, xi + 1, yi + 1, zi - 1, currentValue, currentClusterRef);
          }
          [[fallthrough]];
        case NBH_3D_18_XYZ:
          if ((xi>0) && (yi > 0))
          {
            currentClusterRef = updateClusterRefForNeighbor(slicePtr - xoffset - yoffset, xi - 1, yi - 1, zi, currentValue, currentClusterRef);
          }
          if ((xi<xres - 1) && (yi>0))
          {
            currentClusterRef = updateClusterRefForNeighbor(slicePtr + xoffset - yoffset, xi + 1, yi - 1, zi, currentValue, currentClusterRef);
          }
          if ((zi > 0) && (xi > 0))
          {
            currentClusterRef = updateClusterRefForNeighbor(previousSlicePtr - xoffset, xi - 1, yi, zi - 1, currentValue, currentClusterRef);
          }
          if ((zi > 0) && (xi<xres - 1))
          {
            currentClusterRef = updateClusterRefForNeighbor(previousSlicePtr + xoffset, xi + 1, yi, zi - 1, currentValue, currentClusterRef);
          }
          if ((zi>0) && (yi > 0))
          {
            currentClusterRef = updateClusterRefForNeighbor(previousSlicePtr - yoffset, xi, yi - 1, zi - 1, currentValue, currentClusterRef);
          }
          if ((zi > 0) && (yi<yres - 1))
          {
            currentClusterRef = updateClusterRefForNeighbor(previousSlicePtr + yoffset, xi, yi + 1, zi - 1, currentValue, currentClusterRef);
          }
          [[fallthrough]];
        case NBH_3D_6_XYZ:
          if (zi>0)
          {
            currentClusterRef = updateClusterRefForNeighbor(previousSlicePtr, xi, yi, zi - 1, currentValue, currentClusterRef);
          }
          break;
        case NBH_2D_4_XY:
          // nothing to do, because 2 neighbors are already evaluated before the switch
          break;
        case NBH_COUNT:
          errorCode = ML_PROGRAMMING_ERROR;
          break;
        default:
          errorCode = ML_PROGRAMMING_ERROR;
          break;
        }
 
        // if no matching neighbor found create new clusterRef
        if (currentClusterRef == 0)
        {
          MLdouble imageValue = 0;
 
          if (_parameters.useImageValueAsUserData)
          {
            imageValue = DerivedAlgorithm::getVoxelAsDouble(slicePtr);
          }
 
          currentClusterRef = _clusters->getNewReference(imageValue);
        }
        _clusters->setClusterRef(xi, yi, zi, currentClusterRef);
      }
    }
    return errorCode;
  }
 
 
private:
 
  ClusterAlgorithm& operator=(const ClusterAlgorithm&);
};
 
ML_END_NAMESPACE
 
#endif