mlGraphAnalyser.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_GRAPH_ANALYSER_H
#define ML_GRAPH_ANALYSER_H
 

 
// Include dll-specific settings.
#include "mlVesselGraphSystem.h"
#include "VesselGraphTypes.h"
 
// Include vessel graph class.
#include "VesselGraphDisableWarnings.h"
#include <map>
#include <memory.h>
#include "VesselGraphRestoreWarnings.h"
#include "mlGraphComponents.h"
 
#ifndef VESSELGRAPH_USE_MPIF
#include "AssocGraph.h"
#endif
 
 
 
// Be sure that all your modules are part of the ml. So collisions
// with names in system files or other libraries are minimized.
VESSELGRAPH_BEGIN_NAMESPACE
 
  // ------------------------------------------------------------------
  // ------------------------------------------------------------------
  typedef enum
  {
    DATA_NO         = 0,
    DATA_SCL,
    DATA_BARYCENTER,
    DATA_DISTANCE,
    DATA_FLAG,
    DATA_NODE,
    DATA_GET_FKT,
    DATA_SET_FKT,
    DATA_AV_LEN,
    OUTDATA_CLUSTER,
    OUTDATA_LENGTH,
    OUTDATA_VOLUME,
    OUTDATA_POS,
    OUTDATA_GLOBAL_MIN_DIST,
    OUTDATA_AV_MIN_DIST,
    OUTDATA_ROOT,
  } NData;
 
  typedef std::map<NData, void *> DataMap;
 
  //----------------------------------------------------------
  //----------------------------------------------------------
  typedef enum
  {
    PROC_NONE         = 0x000,
    PROC_LENGTH       = 0x001,
    PROC_VOLUME       = 0x002,
    PROC_WEIGHTED_POS = 0x004,
    PROC_AV_MIN_DIST  = 0x008,
    PROC_SCL          = 0x010,
    PROC_SKELETON_AREA= 0x020,
    PROC_ROOT         = 0x100,
    PROC_ST_VOLUME    = 0x202,
    PROC_ST_POS       = 0x204,
    PROC_ST_NODE_SUM  = 0x208
#ifndef VESSELGRAPH_USE_MPIF
    , PROC_CLUSTER      = 0x400
#endif
  } NProcessMode;
 
 
 

 
  // ------------------------------------------------------------------
  // ------------------------------------------------------------------
  struct CIdxFloat{
    int   _nIdx;
    float _fValue;
 
 
    CIdxFloat(const CIdxFloat & oIF) : _nIdx(oIF._nIdx), _fValue(oIF._fValue)  {};
    CIdxFloat(int nIdx = 0, float fValue = 0.0) : _nIdx(nIdx), _fValue(fValue) {};
    CIdxFloat& operator=(const CIdxFloat &oIF) {if (this != &oIF){ this->_nIdx = oIF._nIdx; this->_fValue=oIF._fValue;} return *this; }
 

    inline bool operator> (const CIdxFloat& oIF) const { if(_fValue != oIF._fValue){return _fValue > oIF._fValue;} else {return _nIdx > oIF._nIdx;} }
    inline bool operator< (const CIdxFloat& oIF) const { if(_fValue != oIF._fValue){return _fValue < oIF._fValue;} else {return _nIdx < oIF._nIdx;} }
    inline bool operator==(const CIdxFloat& oIF) const { return ((_nIdx == oIF._nIdx) && (_fValue == oIF._fValue)); }
  };

  struct CIdxIdx{
    int _n1;
    int _n2;
    CIdxIdx(const CIdxIdx & oII) : _n1(oII._n1), _n2(oII._n2)  {};
    CIdxIdx& operator=(const CIdxIdx &oII) {if (this != &oII){ this->_n1 = oII._n1; this->_n2=oII._n2;} return *this; }
    inline bool operator< (const CIdxIdx& oII) const { if(_n1 != oII._n1){return _n1 < oII._n1;} else {return _n2 < oII._n2;} }
    inline bool operator==(const CIdxIdx& oII) const { return ((_n1 == oII._n1) && (_n2 == oII._n2)); }
  };

 
 
 
  // ------------------------------------------------------------------
  // ------------------------------------------------------------------
  class VESSELGRAPH_EXPORT GraphAnalyser
  {
  public:
    //--------------------------------------------------------
    //--------------------------------------------------------
    explicit GraphAnalyser (Graph * pGraph) : _pGraph(pGraph) {};
 
    //----------------------------------------------------------
    //----------------------------------------------------------
    void measureLength();

    void measureVolume();

    void measureST_Volume(bool bIgnoreOrientation=false     
                         );

    float measureAvMinDistance();

    void measureSkeletonArea(int nAvLength=5);

    Vector3 measureBarycenter();

    void measureST_Barycenter(bool bIgnoreOrientation=false 
                             );
 
    void sumupST_Value(VesselNode::NodeGetFkt getFkt,       
                       VesselNode::NodeSetFkt setFkt,       
                       bool bIgnoreOrientation=false        
                       );

    void measureMaxPathValue(VesselNode* pRootNode, VesselEdge::EdgeGetFkt edgeGetFkt, VesselNode::NodeGetFkt nodeGetFkt, VesselNode::NodeSetFkt nodeSetFkt);

    void searchMaxSubTrees(std::vector<VesselEdge *>::iterator iBegin,
                           std::vector<VesselEdge *>::iterator iEnd,
                           int nNum,                                
                           std::vector<VesselEdge *> * pvMaxEdges,  
                           float * pfActValue,                      
                           float * pfActSum                         
                           );

    inline void searchMaxSubTrees(std::vector<VesselEdge *> vStartEdges,   
                           unsigned int nNum,                       
                           std::vector<VesselEdge *> * pvMaxEdges,  
                           float * pfActValue,                      
                           float * pfActSum                         
                           )
    {
      searchMaxSubTrees(vStartEdges.begin(), vStartEdges.end(), static_cast<int>(nNum), pvMaxEdges, pfActValue, pfActSum);
    }

    inline void searchMaxSubTrees(std::vector<VesselEdge *> vStartEdges,   
                           int nPos,                                
                           unsigned int nNum,                       
                           std::vector<VesselEdge *> * pvMaxEdges,  
                           float * pfActValue,                      
                           float * pfActSum                         
                           )
    {
      searchMaxSubTrees(vStartEdges.begin()+nPos, vStartEdges.begin()+nPos+1, static_cast<int>(nNum), pvMaxEdges, pfActValue, pfActSum);
    }
 
    // Starting at node pRoot, search for two subtrees (originating from a common node) with barycenters at a large distance
    // in separationDirection and large volumes. Push_back() them to vClassifiedEdges ordered by descending position in separationDirection.
    // If (pRoot==NULL) all nodes in the graph are considered.
    void searchSubTreeSeparation(VesselNode* pRoot, const Vector3& separationDirection, std::vector<VesselEdge*>& vClassifiedEdges);

    void searchST_MajorBifurcation(VesselEdge *pStartEdge,          
                                   float fMinorLimit,               
                                   std::vector<VesselEdge *> *pvEdges, 
                                   bool bClear = true               
                                   );

    void searchST_MajorBifurcation(VesselEdge *pStartEdge,          
                                   float fMinorLimitStart,          
                                   float fMinorLimitStep,           
                                   int   nMinorLimitNumber,         
                                   std::vector<VesselEdge *> *pvEdges, 
                                   bool bClear = true               
                                   );
 
 
 
    //----------------------------------------------------------
    //----------------------------------------------------------
    VesselNode* suggestRoot(float fSignificanceLevel = 0.7f, int nPreferedDirection = 0);

    const VesselNode * getParent(const VesselNode *pN1, const VesselNode *pN2) const;
    const VesselNode * getParent(const VesselEdge *pE1, const VesselEdge *pE2) const;
    VesselNode* getParent(VesselNode* pN1, VesselNode* pN2);
    VesselNode* getParent(VesselEdge *pE1, VesselEdge* pE2);
 
  private:
    //----------------------------------------------------------
    //----------------------------------------------------------
    Graph * _pGraph;

    DataMap _MData;
 
 
#ifndef VESSELGRAPH_USE_MPIF
public:
    inline void sumupST_DrainVolume(bool bIgnoreOrientation=false  
                                                            )
    { sumupST_Value (&VesselNode::getDrainVolume, &VesselNode::setDrainVolume, bIgnoreOrientation); }
 

    bool calcDistanceMatrix(MatrixTemplate<double>& dm, std::vector<int>& vId, int nMaxNodes);
 

    bool calcNodeClassifikator(CDoubleArray * pArAssoc, std::vector<int>* pvId, int nMaxNodes);
 

    bool calcNodeBaseFromPosition(const Vector3 & vPosition, const std::vector<int> & vId, std::vector<CIdxFloat> * pvResult);

    bool calcPositionFromNodeBase(std::vector<CIdxFloat> vBase, Vector3 * pvResult);

    double calcSpannedVolume(const std::vector<int> & vId);

    short sampleConnectedCluster();

    void tansferLabels(Skeleton::LABELTYPE NType, float fThresh = 0.1f);
    void tansferLabelsInverse(void);
    void scanSkeletonLabelMinMax( float *min, float *max );

    void collectLeafs(VesselNode *pRoot, std::vector<VesselNode*>& pLeafs);
 
private:
    //----------------------------------------------------------
    //----------------------------------------------------------
    void calcRootPath(int nParentStart, int nParentEnd, int * pFreeIdx, std::vector<PathInfo>::iterator iPathBase, std::vector<PathInfo>::iterator iPathEnd);
 
#endif
 
  };
 
 
 
 
 
  // ------------------------------------------------------------------
  // ------------------------------------------------------------------
  class VESSELGRAPH_EXPORT EdgeAnalyser
  {
  public:
 
    //--------------------------------------------------------
    //--------------------------------------------------------
    explicit EdgeAnalyser (long NMode, DataMap * pData = nullptr);
 
 
 
    //--------------------------------------------------------
    //--------------------------------------------------------
    // function \c operator() provide singular outer access point
    inline void operator() (VesselEdge *pEdge){(this->*_actMeasure)(pEdge);}
 
 
  private:
    //--------------------------------------------------------
    //--------------------------------------------------------
    void measureNothing(VesselEdge * /*pEdge*/) {}

    void measureLength      (VesselEdge *pEdge);

    void measureVolume      (VesselEdge *pEdge);

    void measurePosition    (VesselEdge *pEdge);

    void averageMinDistance (VesselEdge *pEdge);

    void measureST_Volume(VesselEdge *pEdge);

    void measureST_Barycenter(VesselEdge *pEdge);

    void measure_avSkeletonArea(VesselEdge *pEdge);
 

    void sumupST_NodeValue(VesselEdge *pEdge);

 
    //----------------------------------------------------------
    //----------------------------------------------------------
    void (EdgeAnalyser::*_actMeasure) (VesselEdge *pEdge);

    float     _fVoxVolume;
    DataMap * _pData;
    long _NMode;
  };
 
 
 
    typedef enum {
    CMP_LENGTH,
    CMP_VOLUME,
    CMP_STVOLUME,
    CMP_STCENTER_X,
    CMP_STCENTER_Y,
    CMP_STCENTER_Z,
    CMP_STCENTER_PLANE,
    CMP_LABEL
  } CMP_FEATURE;
 
    // ------------------------------------------------------------------
    // ------------------------------------------------------------------
  class compareEdges{
  public:
    explicit compareEdges (CMP_FEATURE NFeature, Vector3 planeNormal = Vector3(0)){
      _planeNormal = planeNormal;
      switch(NFeature){
        case (CMP_LENGTH):
          _actComperator = &compareEdges::compareLength;
          break;
        case (CMP_VOLUME):
          _actComperator = &compareEdges::compareVolume;
          break;
        case (CMP_STVOLUME):
          _actComperator = &compareEdges::compareSTVolume;
          break;
        case (CMP_STCENTER_X):
          _actComperator = &compareEdges::compareSTCenterX;
          break;
        case (CMP_STCENTER_Y):
          _actComperator = &compareEdges::compareSTCenterY;
          break;
        case (CMP_STCENTER_Z):
          _actComperator = &compareEdges::compareSTCenterZ;
          break;
        case (CMP_STCENTER_PLANE):
          _actComperator = &compareEdges::compareSTCenterPlane;
          break;
        case (CMP_LABEL):
        default:
          _actComperator = &compareEdges::compareLabel;
          break;
      }
    }
 
    int operator() (const VesselEdge *pE1, const VesselEdge *pE2) const
      { return (this->*_actComperator)(pE1, pE2); }
 
  private:
    int (compareEdges::*_actComperator) (const VesselEdge *pE1, const VesselEdge *pE2) const;
    Vector3 _planeNormal;
 
    int compareLabel         (const VesselEdge *pE1, const VesselEdge *pE2) const { return pE1->getLabel()             > pE2->getLabel(); }
    int compareLength        (const VesselEdge *pE1, const VesselEdge *pE2) const { return pE1->getLength()            > pE2->getLength(); }
    int compareVolume        (const VesselEdge *pE1, const VesselEdge *pE2) const { return pE1->getVolume()            > pE2->getVolume(); }
    int compareSTVolume      (const VesselEdge *pE1, const VesselEdge *pE2) const { return pE1->getSTVolume()          > pE2->getSTVolume(); }
    int compareSTCenterX     (const VesselEdge *pE1, const VesselEdge *pE2) const { return (pE1->getSTBarycenter())[0] > (pE2->getSTBarycenter())[0]; }
    int compareSTCenterY     (const VesselEdge *pE1, const VesselEdge *pE2) const { return (pE1->getSTBarycenter())[1] > (pE2->getSTBarycenter())[1]; }
    int compareSTCenterZ     (const VesselEdge *pE1, const VesselEdge *pE2) const { return (pE1->getSTBarycenter())[2] > (pE2->getSTBarycenter())[2]; }
    int compareSTCenterPlane (const VesselEdge *pE1, const VesselEdge *pE2) const { return (pE1->getSTBarycenter()-pE2->getSTBarycenter()).dot(_planeNormal) > 0; }
 
  };
 
 
  // ------------------------------------------------------------------
  // ------------------------------------------------------------------
  class VESSELGRAPH_EXPORT NodeAnalyser
  {
  public:
 
    //--------------------------------------------------------
    //--------------------------------------------------------
    explicit NodeAnalyser (long NMode, DataMap * pData = nullptr) :
    _actMeasure (&NodeAnalyser::doNothing),
    _pData(pData)
    {
      switch(NMode){
        case PROC_ROOT:
          if(_pData != nullptr){
            if((*pData)[OUTDATA_ROOT] != nullptr){
              *( reinterpret_cast<long *>((*pData)[OUTDATA_ROOT])) = 0;
              _fDist = 0.0;
              _actMeasure = &NodeAnalyser::suggestRoot;
            }
          } else {
            _actMeasure = &NodeAnalyser::doNothing;
          }
          break;
#ifndef VESSELGRAPH_USE_MPIF
        case PROC_CLUSTER:
          if(_pData != nullptr){
            if((*pData)[OUTDATA_CLUSTER] != nullptr){
              _actMeasure = &NodeAnalyser::setClusterLabel;
            }
          } else {
            _actMeasure = &NodeAnalyser::doNothing;
          }
          break;
#endif
        default:
          _actMeasure = &NodeAnalyser::doNothing;
      }
    }

 
 
    //--------------------------------------------------------
    //--------------------------------------------------------
    // function \c operator() provide singular outer access point
    inline void operator() (VesselNode *pNode){(this->*_actMeasure)(pNode);}
 
 
  private:
    //--------------------------------------------------------
    //--------------------------------------------------------
    void doNothing          (VesselNode * /*pNode*/) {}

    void suggestRoot        (VesselNode *pNode);
 
#ifndef VESSELGRAPH_USE_MPIF
    void setClusterLabel    (VesselNode *pNode);
#endif
 
    //----------------------------------------------------------
    //----------------------------------------------------------
    void (NodeAnalyser::*_actMeasure) (VesselNode *pNode);

    DataMap * _pData;
 
    float _fDist;
  };
 
  // ------------------------------------------------------------------
  // ------------------------------------------------------------------
  class measureSkeletonAgent
  {

    class Process
    {
      typedef void (measureSkeletonAgent::*ProcessFkt) (const Skeleton & oSkeleton, Process * pProcess);
    public:
        ProcessFkt      oProcess = &measureSkeletonAgent::doNothing; 
        void            *pResult = nullptr;         
    };

    friend class Process;

    typedef std::map<NProcessMode, Process>::iterator IProcess;
 
  public:
 
 
    //--------------------------------------------------------
    //--------------------------------------------------------
    measureSkeletonAgent(int nMode, DataMap * pData = nullptr) :
        _fLastLength(0.0),
        _pScale(nullptr),
        _pData(pData)
      {
        _MProcess.erase(_MProcess.begin(), _MProcess.end());
 
        if(_pData != nullptr){
          if(nMode & PROC_LENGTH)
          {
            Process oProcess;
            oProcess.pResult = (*_pData)[OUTDATA_LENGTH];
            if(nMode & PROC_SCL){
              _pScale = static_cast<Vector3 *>((*_pData)[DATA_SCL]);
              oProcess.oProcess = &measureSkeletonAgent::ProcessL1Scale;
            }
            else {
              oProcess.oProcess =&measureSkeletonAgent::ProcessL1;
            }
            _MProcess[PROC_LENGTH] = oProcess;
          }
 
          if(nMode & PROC_VOLUME)
          {
            Process oProcess;
            oProcess.pResult = (*_pData)[OUTDATA_VOLUME];
            oProcess.oProcess = &measureSkeletonAgent::ProcessV1;
            _MProcess[PROC_VOLUME] = oProcess;
          }
 
          if(nMode & PROC_WEIGHTED_POS)
          {
            Process oProcess;
            oProcess.pResult = (*_pData)[OUTDATA_POS];
            oProcess.oProcess = &measureSkeletonAgent::ProcessP1;
            _MProcess[PROC_WEIGHTED_POS] = oProcess;
          }
 
          if(nMode & PROC_AV_MIN_DIST)
          {
            Process oProcess;
            oProcess.pResult = (*_pData)[OUTDATA_AV_MIN_DIST];
            oProcess.oProcess = &measureSkeletonAgent::ProcessD1;
            _MProcess[PROC_AV_MIN_DIST] = oProcess;
          }
        }
      }

 
 
    //--------------------------------------------------------
    //--------------------------------------------------------
    inline void operator() (const Skeleton & oSkeleton)
    {
      for(IProcess iProcess = _MProcess.begin(); iProcess != _MProcess.end(); ++iProcess){
         (this->*((iProcess->second).oProcess))(oSkeleton, &(iProcess->second));
      }
    }
    inline void operator() (Skeleton* oSkeleton)
    {
      for(IProcess iProcess = _MProcess.begin(); iProcess != _MProcess.end(); ++iProcess){
        (this->*((iProcess->second).oProcess))(*oSkeleton, &(iProcess->second));
      }
    }

 
 
 
  private:
 
    //--------------------------------------------------------
    //--------------------------------------------------------
    void doNothing       (const Skeleton & /*oSkeleton*/, Process * /*pProcess*/) {}

    void ProcessL1       (const Skeleton & oSkeleton, Process * pProcess);
    void ProcessL1Scale  (const Skeleton & oSkeleton, Process * pProcess);
    void ProcessL2       (const Skeleton & oSkeleton, Process * pProcess);
    void ProcessL2Scale  (const Skeleton & oSkeleton, Process * pProcess);
    void ProcessL3       (const Skeleton & oSkeleton, Process * pProcess);
    void ProcessL3Scale  (const Skeleton & oSkeleton, Process * pProcess);

    void ProcessV1       (const Skeleton & oSkeleton, Process * pProcess);

    void ProcessP1       (const Skeleton & oSkeleton, Process * pProcess);

    void ProcessD1      (const Skeleton & oSkeleton, Process * pProcess);
 

 
    //----------------------------------------------------------
    //----------------------------------------------------------
    std::map<NProcessMode, Process> _MProcess;
 

    float _fLastLength;

    std::list<Vector3>  _LPoints;
    std::list<Vector3>::iterator IBack;
    std::list<Vector3>::iterator IMiddle;
    std::list<Vector3>::iterator IFront;
 
    // interpolated Skeleton position
    Vector3      _vMean;
    Vector3      * _pScale;

    DataMap   * _pData;
  };
 
VESSELGRAPH_END_NAMESPACE
 
#endif // __GRAPHANALYSER_H