MeVisLab/Standard/Sources/ML/MLVesselGraph/mlGraphAnalyser.h

Go to the documentation of this file.

// **InsertLicense** code
//----------------------------------------------------------------------------------

//----------------------------------------------------------------------------------
#ifndef __GraphAnalyser_H
#define __GraphAnalyser_H


// Include dll-specific settings.
#include "mlVesselGraphSystem.h"

// Include ml Macros for implementation support.
#include "mlDebug.h"

// Include vessel graph class.
#include "mlSystemWarningsDisable.h"
#include <map>
#include <memory.h>
#include "mlSystemWarningsRestore.h"
#include "mlGraph.h"
#include "mlGraphComponents.h"
#include "mlFloatingPointVector.h"
#include "mlMatrixTemplate.h"
#include "AssocGraph.h"



// Be sure that all your modules are part of the ml. So collisions
// with names in system files or other libraries are minimized.
ML_START_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,
    PROC_CLUSTER      = 0x400
  } 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 

                             );

    inline void sumupST_DrainVolume(bool bIgnoreOrientation=false  

                                                            )
    { sumupST_Value (&VesselNode::getDrainVolume, &VesselNode::setDrainVolume, bIgnoreOrientation); }

    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(), (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, (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               
                                   );


    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);



    //----------------------------------------------------------

    //----------------------------------------------------------
    VesselNode* suggestRoot(float fSignificanceLevel = 0.7, int nPreferedDirection = 0);


    short sampleConnectedCluster();

    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);

    void tansferLabels(Skeleton::LABELTYPE NType, float fThresh = 0.1);
    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);
    


    //----------------------------------------------------------

    //----------------------------------------------------------
    Graph * _pGraph;

    DataMap _MData;

  };





  // ------------------------------------------------------------------
  // ------------------------------------------------------------------
  class VESSELGRAPH_EXPORT EdgeAnalyser
  {
  public:    
    
    //--------------------------------------------------------

    //--------------------------------------------------------
    explicit EdgeAnalyser (long NMode, DataMap * pData = NULL);



    //--------------------------------------------------------

    //--------------------------------------------------------
    // 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 = NULL) : 
    _actMeasure (&NodeAnalyser::doNothing),
    _pData(pData)
    { 
      switch(NMode){
        case PROC_ROOT:
          if(_pData != NULL){
            if((*pData)[OUTDATA_ROOT] != NULL){
              *((long *) (*pData)[OUTDATA_ROOT]) = 0;
              _fDist = 0.0;
              _actMeasure = &NodeAnalyser::suggestRoot;
            }
          } else {
            _actMeasure = &NodeAnalyser::doNothing;
          }
          break;
        case PROC_CLUSTER:
          if(_pData != NULL){
            if((*pData)[OUTDATA_CLUSTER] != NULL){
              _actMeasure = &NodeAnalyser::setClusterLabel;
            }
          } else {
            _actMeasure = &NodeAnalyser::doNothing;
          }
          break;
        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);

    void setClusterLabel    (VesselNode *pNode);

    //----------------------------------------------------------

    //----------------------------------------------------------
    void (NodeAnalyser::*_actMeasure) (VesselNode *pNode);

    DataMap * _pData;

    float _fDist;
  };

  // ------------------------------------------------------------------
  // ------------------------------------------------------------------
  class measureSkeletonAgent
  {

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

        ProcessFkt      oProcess;         
        void            *pResult;         
    };

    friend class Process;

    typedef std::map<NProcessMode, Process>::iterator IProcess;

  public:


    //--------------------------------------------------------

    //--------------------------------------------------------
    measureSkeletonAgent(int nMode, DataMap * pData = NULL) : 
        _fLastLength(0.0),
        _pScale(NULL),
        _pData(pData)
      { 
        _MProcess.erase(_MProcess.begin(), _MProcess.end());

        if(_pData != NULL){
          if(nMode & PROC_LENGTH)
          {
            Process oProcess;
            oProcess.pResult = (*_pData)[OUTDATA_LENGTH];
            if(nMode & PROC_SCL){
              _pScale = (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;
  };

  ML_END_NAMESPACE

#endif // __GRAPHANALYSER_H