## =============
OPTION(BUILD_PLUGINS "Build plugins" ON)
+OPTION(BUILD_TempAirwaysAppli "Build temporary appli" OFF)
## =========================
## == Include directories ==
cmake
lib
plugins
- # appli
+ appli
)
## eof - $RCSfile$
-IF(BUILD_EXAMPLES)
+#IF(BUILD_EXAMPLES)
+#SUBDIRS(
+#examples
+## fpaLab
+#)
+#ENDIF(BUILD_EXAMPLES)
+
+IF(BUILD_TempAirwaysAppli)
SUBDIRS(
- examples
- ## fpaLab
+ TempAirwaysAppli
)
-ENDIF(BUILD_EXAMPLES)
+ENDIF(BUILD_TempAirwaysAppli)
## eof - $RCSfile$
--- /dev/null
+## =============================
+## = Set names and directories =
+## =============================
+
+SET(lib_NAME AirwaysLib)
+
+INCLUDE_DIRECTORIES(
+ ${PROJECT_SOURCE_DIR}/appli/TempAirwaysAppli/MathLib
+ ${PROJECT_SOURCE_DIR}/appli/TempAirwaysAppli/AirwaysLib
+ ${PROJECT_BINARY_DIR}/appli/TempAirwaysAppli/MathLib
+ ${PROJECT_BINARY_DIR}/appli/TempAirwaysAppli/AirwaysLib
+ )
+
+## ===============
+## = Source code =
+## ===============
+
+FILE(GLOB lib_HEADERS_H "*.h")
+FILE(GLOB lib_HEADERS_HPP "*.hpp")
+FILE(GLOB lib_HEADERS_HXX "*.hxx")
+FILE(GLOB lib_SOURCES_C "*.c")
+FILE(GLOB lib_SOURCES_CPP "*.cpp")
+FILE(GLOB lib_SOURCES_CXX "*.cxx")
+
+## =====================
+## = Compilation rules =
+## =====================
+
+ADD_LIBRARY(
+ ${lib_NAME}
+ SHARED
+ ${lib_SOURCES_C}
+ ${lib_SOURCES_CPP}
+ ${lib_SOURCES_CXX}
+ )
+GENERATE_EXPORT_HEADER(
+ ${lib_NAME}
+ BASE_NAME ${lib_NAME}
+ EXPORT_MACRO_NAME ${lib_NAME}_EXPORT
+ EXPORT_FILE_NAME ${CMAKE_CURRENT_BINARY_DIR}/${lib_NAME}_Export.h
+ STATIC_DEFINE ${lib_NAME}_BUILT_AS_STATIC
+ )
+TARGET_LINK_LIBRARIES(${lib_NAME} MathLib ${ITK_LIBRARIES} ${VTK_LIBRARIES})
+
+## eof - $RCSfile$
--- /dev/null
+/*
+ * airwaysEdge.cxx
+ *
+ * Created on: May 12, 2014
+ * Author: caceres
+ */
+
+#include "airwaysEdge.h"
+#include <iostream>
+
+namespace airways
+{
+
+Edge::Edge() : m_id(0), m_source(NULL), m_target(NULL), m_mark(false), m_angle(0), m_length(
+ 0), m_eDistance(0), m_aRadius(0), m_minRadius(0), m_maxRadius(0)
+{
+
+}
+
+Edge::Edge(Edge* edge) : m_id(0), m_source(NULL), m_target(NULL), m_mark(false), m_angle(0), m_length(
+ 0), m_eDistance(0), m_aRadius(0), m_minRadius(0), m_maxRadius(0)
+{
+ this->m_id = edge->m_id;
+ this->m_source = edge->m_source;
+ this->m_target = edge->m_target;
+ this->m_angle = edge->m_angle;
+ this->m_length = edge->m_length;
+ this->m_eDistance = edge->m_eDistance;
+ this->m_vec_pair_posVox_rad = edge->m_vec_pair_posVox_rad;
+ this->m_aRadius = edge->m_aRadius;
+ this->m_minRadius = edge->m_minRadius;
+ this->m_maxRadius = edge->m_maxRadius;
+}
+
+Edge::~Edge()
+{
+}
+
+int Edge::GetAngle() const
+{
+ return this->m_angle;
+}
+
+double Edge::GetARadius() const
+{
+ return this->m_aRadius;
+}
+
+double Edge::GetEDistance() const
+{
+ return this->m_eDistance;
+}
+
+double Edge::GetLength() const
+{
+ return this->m_length;
+}
+
+double Edge::GetMaxRadius() const
+{
+ return this->m_maxRadius;
+}
+
+double Edge::GetMinRadius() const
+{
+ return this->m_minRadius;
+}
+
+Node* Edge::GetSource() const
+{
+ return this->m_source;
+}
+
+Node* Edge::GetTarget() const
+{
+ return this->m_target;
+}
+
+const vec_pair_posVox_rad& Edge::GetEdgeInfo() const
+{
+ return m_vec_pair_posVox_rad;
+}
+
+bool Edge::IsMarked() const
+{
+ return this->m_mark;
+}
+
+bool Edge::IsPointInfluencedByEdge(float point_x, float point_y, float point_z) const
+{
+ // Variables
+ float minDistance = -1.0;
+ bool first = true;
+ float distanceTemp = 0.0;
+ bool influenced = false;
+ Vec3 vector(point_x, point_y, point_z);
+
+ // Iterate over all voxel in the edge
+ vec_pair_posVox_rad::const_iterator it = this->m_vec_pair_posVox_rad.begin();
+ for (; it != this->m_vec_pair_posVox_rad.end() && !influenced; ++it)
+ {
+ // Get the vector from the compared source to each voxel
+ Vec3 voxel_actual = ((*it).first);
+ //Vec3 source_actual = this->m_source->GetCoords( );
+ //Vec3 vector_actual(voxel_actual - source_actual);
+
+ // Get the difference vector
+ Vec3 vector_diff = vector - voxel_actual;
+
+ // Get the distance
+ distanceTemp = vector_diff.Norm();
+
+ // Get the minimum
+ if(first)
+ {
+ minDistance = distanceTemp;
+ first = false;
+ if(minDistance < ((*it).second*1.8))
+ influenced = true;
+ }
+ else if(distanceTemp < minDistance)
+ {
+ minDistance = distanceTemp;
+ if(minDistance < ((*it).second)*1.8)
+ influenced = true;
+ }
+ }
+
+ return influenced;
+}
+
+void Edge::SetAngle(const double& angle)
+{
+ this->m_angle = angle;
+}
+
+void Edge::SetARadius(const double& aRadius)
+{
+ this->m_aRadius = aRadius;
+}
+
+void Edge::SetEDistance(const double& eDistance)
+{
+ this->m_eDistance = eDistance;
+}
+
+void Edge::SetLength(const double& length)
+{
+ this->m_length = length;
+}
+
+void Edge::SetMaxRadius(const unsigned int& maxRadius)
+{
+ this->m_maxRadius = maxRadius;
+}
+
+void Edge::SetMinRadius(const unsigned int& minRadius)
+{
+ this->m_minRadius = minRadius;
+}
+
+void Edge::SetSource(Node* source)
+{
+ this->m_source = source;
+}
+
+void Edge::SetTarget(Node* target)
+{
+ this->m_target = target;
+}
+
+void Edge::AddSkeletonPairInfo(const pair_posVox_rad& skPairInfo)
+{
+ this->m_vec_pair_posVox_rad.push_back(skPairInfo);
+}
+
+void Edge::SetSkeletonPairVector(const vec_pair_posVox_rad& skPInfoVector)
+{
+ this->m_vec_pair_posVox_rad = skPInfoVector;
+ UpdateEdgeInfo();
+}
+
+void Edge::UpdateEdgeInfo()
+{
+ if (this->m_vec_pair_posVox_rad.empty())
+ {
+ std::cout << " ------------------------ This edge has no information" << std::endl;
+ std::cout << "Source:" << this->m_source->GetCoords() << std::endl;
+ std::cout << "Target:" << this->m_target->GetCoords() << std::endl;
+ return;
+ }
+ double min = std::numeric_limits<unsigned int>::max();
+ double max = std::numeric_limits<unsigned int>::min();
+ double average = 0;
+ vec_pair_posVox_rad::iterator it = this->m_vec_pair_posVox_rad.begin();
+ for (; it != this->m_vec_pair_posVox_rad.end(); ++it)
+ {
+ pair_posVox_rad info = *it;
+ average += info.second;
+ if (info.second < min)
+ min = info.second;
+ else if (info.second > max)
+ max = info.second;
+ }
+ this->m_aRadius = average / this->m_vec_pair_posVox_rad.size();
+ this->m_minRadius = min;
+ this->m_maxRadius = max;
+}
+
+double Edge::CompareWith(Edge* edge)
+{
+ //std::cout << "Comparing ... " << std::endl;
+ // To compared the edges, the distance in each component
+ // is calculated and a distance function between edges is calculated.
+ // The components are:
+ // - Each component of the quaternion (son-parent quaterion). (4 components)
+ // - Edge distance
+ // - Edge average radius
+
+ // Difference in each quaternion component
+ // p
+ // |
+ // |
+ // |
+ // s
+ // /
+ // /
+ // /
+ // g
+ /*const Edge* edgeActualParent = this->m_source->GetEdge( );
+ Vec3 sourceActualParent = edgeActualParent->m_source->GetCoords( );
+ Vec3 targetActualParent = edgeActualParent->m_target->GetCoords( );
+ Vec3 sourceActual = this->m_source->GetCoords( );
+ Vec3 targetActual = this->m_target->GetCoords( );
+
+ Vec3 sp_actual( sourceActualParent - targetActualParent );
+ Vec3 sg_actual( targetActual - sourceActual );
+
+ const Edge* edgeComparedParent = edge.m_source->GetEdge( );
+ Vec3 sourceComparedParent = edgeComparedParent->m_source->GetCoords( );
+ Vec3 targetComparedParent = edgeComparedParent->m_target->GetCoords( );
+ Vec3 sourceCompared = edge.m_source->GetCoords( );
+ Vec3 targetCompared = edge.m_target->GetCoords( );
+
+ Vec3 sp_compared( sourceComparedParent - targetComparedParent );
+ Vec3 sg_compared( targetCompared - sourceCompared );
+
+ Quaternion q_actual(sp_actual, sg_actual);
+ Quaternion q_compared(sp_compared, sg_compared);
+
+ float dif_r = q_actual.getR() - q_compared.getR();
+ float dif_i = q_actual.getI() - q_compared.getI();
+ float dif_j = q_actual.getJ() - q_compared.getJ();
+ float dif_k = q_actual.getK() - q_compared.getK();
+ float dif_distance = this->m_length - edge.m_length;
+ float dif_radius = this->m_aRadius - edge.m_aRadius;
+ */
+ double distanceEdge = GetDistanceToEdge(edge);
+
+ //std::cout << "Comparing ... OK" << std::endl;
+
+ return distanceEdge;
+}
+
+float Edge::GetDistanceToEdge(Edge* edge)
+{
+ // a. The parent edges must be aligned before calculate the distance.
+ // b. Parent end-points are aligned also.
+ // Steps a and b represent a rotation and a translation respectively.
+ // c. For each point in the actual Edge (this) the closest point in the
+ // compared edge is found. This distance is calculated and then averaged
+ // for all the points. Must be done in both ways?
+ // We can take the maximum of both distances max(D(a,b),D(b,a))
+
+ // Variables
+ float averageDistance = 0.0;
+ float numPoints = 0;
+
+ // Get the parent Edges
+ /*const Edge* edgeActualParent = this->m_source->GetEdge( );
+ Vec3 sourceActualParent = edgeActualParent->m_source->GetCoords( );
+ Vec3 targetActualParent = edgeActualParent->m_target->GetCoords( );
+ Vec3 sourceActual = this->m_source->GetCoords( );
+ Vec3 targetActual = this->m_target->GetCoords( );
+
+ Vec3 sp_actual( sourceActualParent - targetActualParent );
+ Vec3 sg_actual( targetActual - sourceActual );
+
+ const Edge* edgeComparedParent = edge.m_source->GetEdge( );
+ Vec3 sourceComparedParent = edgeComparedParent->m_source->GetCoords( );
+ Vec3 targetComparedParent = edgeComparedParent->m_target->GetCoords( );
+ Vec3 sourceCompared = edge.m_source->GetCoords( );
+ Vec3 targetCompared = edge.m_target->GetCoords( );
+
+ Vec3 sp_compared( sourceComparedParent - targetComparedParent );
+ Vec3 sg_compared( targetCompared - sourceCompared );
+
+ //Get the quaternion between parents
+ Quaternion q_parent(sp_actual, sp_compared);
+
+ // Get the translation vector to translate the compared parent to the actual parent
+ const Vec3 transParents( targetActualParent - targetComparedParent);
+ */
+
+ // Iterate over the voxels in the compared edge
+ for(vec_pair_posVox_rad::iterator it = edge->m_vec_pair_posVox_rad.begin(); it != edge->m_vec_pair_posVox_rad.end(); ++it)
+ {
+ // Get the vector from the compared source to each voxel in the compared edge
+ Vec3 voxel_actual = (*it).first;
+ /*Vec3 vector_compared(voxel_actual - sourceCompared);
+
+ // Rotate the vector using the parental quaternion
+ Vec3 vector_compared_rotated = q_parent.rotateVector(vector_compared);
+ */
+
+ // Find the distance to the closest point in the actual edge
+ //averageDistance += getDistanceToClosestPoint(sourceActual+vector_compared_rotated);
+ averageDistance += this->GetSmallestDistanceToPoint(voxel_actual);
+
+ ++numPoints;
+ }
+
+ // Get the average value
+ if(numPoints > 0)
+ averageDistance = averageDistance / numPoints;
+
+ return averageDistance;
+}
+
+float Edge::GetDistanceToTranslatedEdge(Edge* edge, Vec3 vector_translation)
+{
+ // a. The parent edges must be aligned before calculate the distance.
+ // b. Parent end-points are aligned also.
+ // Steps a and b represent a rotation and a translation respectively.
+ // c. For each point in the actual Edge (this) the closest point in the
+ // compared edge is found. This distance is calculated and then averaged
+ // for all the points. Must be done in both ways?
+ // We can take the maximum of both distances max(D(a,b),D(b,a))
+
+ // Variables
+ float averageDistance = 0.0;
+ float numPoints = 0;
+
+ // Iterate over the voxels in the compared edge
+ for(vec_pair_posVox_rad::iterator it = edge->m_vec_pair_posVox_rad.begin(); it != edge->m_vec_pair_posVox_rad.end(); ++it)
+ {
+ // Get the vector from the compared source to each voxel in the compared edge
+ Vec3 voxel_actual = (*it).first;
+
+ // Translate the voxel
+ Vec3 voxel_translated = voxel_actual + vector_translation;
+
+ // Find the distance to the closest point in the actual edge
+ averageDistance += this->GetSmallestDistanceToPoint(voxel_translated);
+
+ ++numPoints;
+ }
+
+ // Get the average value
+ if(numPoints > 0)
+ averageDistance = averageDistance / numPoints;
+
+ return averageDistance;
+}
+
+float Edge::GetDistanceWeigthedToTranslatedEdge(Edge* edge, Vec3 vector_translation)
+{
+ // a. The parent edges must be aligned before calculate the distance.
+ // b. Parent end-points are aligned also.
+ // Steps a and b represent a rotation and a translation respectively.
+ // c. For each point in the actual Edge (this) the closest point in the
+ // compared edge is found. This distance is calculated and then averaged
+ // for all the points. Must be done in both ways?
+ // We can take the maximum of both distances max(D(a,b),D(b,a))
+
+ // Variables
+ double alfa = 4; // Weight for duplicated closest points
+ float distance_average = 0.0;
+ float distance_actual = 0.0;
+ float numPoints = 0;
+ std::map<Vec3, vec_pair_posVox_rad > map_vector_pair_vector_distance; // Map to save the correspondences between the given edge voxels and this voxels, and the distance between them.
+
+ // Iterate over the voxels in the compared edge
+ for(vec_pair_posVox_rad::iterator it = edge->m_vec_pair_posVox_rad.begin(); it != edge->m_vec_pair_posVox_rad.end(); ++it)
+ {
+ // Get the vector from the compared source to each voxel in the compared edge
+ Vec3 voxel_actual = (*it).first;
+
+ // Translate the voxel
+ Vec3 voxel_translated = voxel_actual + vector_translation;
+
+ // Find the distance to the closest point in the actual edge
+ Vec3 point_closest = this->GetClosestPoint(voxel_translated, distance_actual);
+ distance_average += distance_actual;
+
+ // Add the link
+ pair_posVox_rad pair_vector_distance(voxel_actual, distance_actual);
+ map_vector_pair_vector_distance[point_closest].push_back(pair_vector_distance);
+
+ ++numPoints;
+ }
+
+ // Add the weight to shared closest points
+ float distance_weighted = 0.0;
+ for(std::map<Vec3, vec_pair_posVox_rad>::iterator it = map_vector_pair_vector_distance.begin(); it != map_vector_pair_vector_distance.end(); ++it)
+ {
+ vec_pair_posVox_rad vector_pair_vector_distance = (*it).second;
+ if(vector_pair_vector_distance.size() > 1)
+ {
+ for(vec_pair_posVox_rad::iterator it_vector = vector_pair_vector_distance.begin(); it_vector != vector_pair_vector_distance.end(); ++it_vector)
+ distance_weighted += (*it_vector).second * alfa;
+ }
+ else
+ {
+ for(vec_pair_posVox_rad::iterator it_vector = vector_pair_vector_distance.begin(); it_vector != vector_pair_vector_distance.end(); ++it_vector)
+ distance_weighted += (*it_vector).second;
+ }
+ }
+
+ // Get the average value
+ if(numPoints > 0)
+ {
+ distance_average = distance_average / numPoints;
+ distance_weighted = distance_weighted / numPoints;
+ }
+
+ //std::cout << "[DistanceW;DistAvg][" << distance_weighted << ";" << distance_average << "];;";
+ //return distance_average;
+ return distance_weighted;
+}
+
+float Edge::GetDistanceToPoint(Vec3 point)
+{
+ // Variables
+ bool first = true;
+ float distance = 0.0;
+
+ // Iterate over all voxel in the edge
+ vec_pair_posVox_rad::iterator it = this->m_vec_pair_posVox_rad.begin();
+ for (; it != this->m_vec_pair_posVox_rad.end(); ++it)
+ {
+ // Get the vector from the compared source to each voxel
+ Vec3 voxel_actual = ((*it).first);
+
+ // Get the difference vector
+ Vec3 vector_diff = point - voxel_actual;
+
+ // Get the distance
+ distance += vector_diff.Norm();
+ }
+
+ return distance;
+}
+
+
+float Edge::GetSmallestDistanceToPoint(Vec3 vector)
+{
+ // Variables
+ float minDistance = -1.0;
+ bool first = true;
+ float distanceTemp = 0.0;
+
+ // Iterate over all voxel in the edge
+ vec_pair_posVox_rad::iterator it = this->m_vec_pair_posVox_rad.begin();
+ for (; it != this->m_vec_pair_posVox_rad.end(); ++it)
+ {
+ // Get the vector from the compared source to each voxel
+ Vec3 voxel_actual = ((*it).first);
+ //Vec3 source_actual = this->m_source->GetCoords( );
+ //Vec3 vector_actual(voxel_actual - source_actual);
+
+ // Get the difference vector
+ Vec3 vector_diff = vector - voxel_actual;
+
+ // Get the distance
+ distanceTemp = vector_diff.Norm();
+
+ // Get the minimum
+ if(first)
+ {
+ minDistance = distanceTemp;
+ first = false;
+ }
+ else if(distanceTemp < minDistance)
+ minDistance = distanceTemp;
+ }
+
+ return minDistance;
+}
+
+Vec3 Edge::GetClosestPoint(Vec3 vector, float& distance)
+{
+ // Variables
+ Vec3 point_closest;
+ bool first = true;
+ float distanceTemp = 0.0;
+
+ // Iterate over all voxel in the edge
+ vec_pair_posVox_rad::iterator it = this->m_vec_pair_posVox_rad.begin();
+ for (; it != this->m_vec_pair_posVox_rad.end(); ++it)
+ {
+ // Get the vector from the compared source to each voxel
+ Vec3 voxel_actual = ((*it).first);
+ //Vec3 source_actual = this->m_source->GetCoords( );
+ //Vec3 vector_actual(voxel_actual - source_actual);
+
+ // Get the difference vector
+ Vec3 vector_diff = vector - voxel_actual;
+
+ // Get the distance
+ distanceTemp = vector_diff.Norm();
+
+ // Get the minimum
+ if(first)
+ {
+ point_closest = voxel_actual;
+ distance = distanceTemp;
+ first = false;
+ }
+ else if(distanceTemp < distance)
+ {
+ point_closest = voxel_actual;
+ distance = distanceTemp;
+ }
+ }
+
+ return point_closest;
+}
+
+Edge* Edge::ConcatenateToSuperior(Edge* superior)
+{
+ if(superior == NULL)
+ return new Edge(this);
+
+ // Check concatenation condition in the nodes
+ if(superior->GetTarget()->GetCoords() != this->GetSource()->GetCoords())
+ {
+ std::cout << "Edge::ConcatenateToSuperior - superior.target != actualEdge.source. Actual idNode:" << this->m_id << std::endl;
+ std::cout << "[superior.target;actualEdge.source]: [" << superior->GetTarget()->GetCoords() << ";" << this->GetSource()->GetCoords() << "]" << std::endl;
+ return NULL;
+ }
+
+ // Check concatenation condition in the edge information
+ vec_pair_posVox_rad vectorInfo_thisEdge = this->GetEdgeInfo(); // Actual edge information
+ vec_pair_posVox_rad vectorInfo_superiorEdge = superior->GetEdgeInfo(); // Superior edge information
+
+ // Check that last superior edge position is equal to first initial actual edge position
+ if(vectorInfo_superiorEdge[vectorInfo_superiorEdge.size()-1].first != vectorInfo_thisEdge[0].first)
+ {
+ std::cout << "Edge::ConcatenateToSuperior - superior.info[end] != actualEdge.info[begin]. Actual idNode:" << this->m_id << std::endl;
+ return NULL;
+ }
+
+ // Change the superior edge target, the superior length, and the euclidian distance
+ superior->SetTarget(this->GetTarget());
+ superior->SetLength(superior->GetLength()+this->GetLength()-1);
+ superior->SetEDistance(superior->GetEDistance()+this->GetEDistance());
+
+ // Add the position information
+ vec_pair_posVox_rad::iterator it_actualInfo = vectorInfo_thisEdge.begin();
+ ++it_actualInfo; // Skip the first position because is it shared between both edges
+ for(; it_actualInfo != vectorInfo_thisEdge.end(); ++it_actualInfo)
+ vectorInfo_superiorEdge.push_back((*it_actualInfo));
+
+ // Set the new information
+ superior->SetSkeletonPairVector(vectorInfo_superiorEdge);
+
+ return superior;
+
+}
+
+Edge& Edge::operator=(Edge& edge)
+{
+ this->m_mark = edge.m_mark;
+ this->m_angle = edge.m_angle;
+ this->m_length = edge.m_length;
+ this->m_eDistance = edge.m_eDistance;
+ this->m_aRadius = edge.m_aRadius;
+ this->m_minRadius = edge.m_minRadius;
+ this->m_maxRadius = edge.m_maxRadius;
+ this->m_source = edge.m_source;
+ this->m_target = edge.m_target;
+ this->m_vec_pair_posVox_rad = edge.m_vec_pair_posVox_rad;
+ return *this;
+}
+
+const Edge& Edge::operator=(const Edge& edge)
+{
+ this->m_mark = edge.m_mark;
+ this->m_angle = edge.m_angle;
+ this->m_length = edge.m_length;
+ this->m_eDistance = edge.m_eDistance;
+ this->m_aRadius = edge.m_aRadius;
+ this->m_minRadius = edge.m_minRadius;
+ this->m_maxRadius = edge.m_maxRadius;
+ this->m_source = edge.m_source;
+ this->m_target = edge.m_target;
+ this->m_vec_pair_posVox_rad = edge.m_vec_pair_posVox_rad;
+ return *this;
+}
+
+bool Edge::operator==(const Edge& edge)
+ {
+ bool cmp1 = true;
+ bool cmp2 = true;
+ bool cmp3 = true;
+ bool cmp4 = true;
+ bool cmp5 = true;
+ bool cmp6 = true;
+ /*if (this->m_angle != edge.m_angle)
+ {
+ double max =
+ this->m_angle >= edge.m_angle ? this->m_angle : edge.m_angle;
+ double tol = 0.3 * max;
+ double diff = this->m_angle - edge.m_angle;
+ if (!((-tol <= diff) && (diff <= tol)))
+ cmp1 = false;
+ }*/
+ if (this->m_length != edge.m_length)
+ {
+ double max =
+ this->m_length >= edge.m_length ? this->m_length : edge.m_length;
+ double tol = 0.3 * max;
+ double diff = abs(this->m_length - edge.m_length);
+ if (diff > tol)
+ cmp2 = false;
+ }
+ if (this->m_eDistance != edge.m_eDistance)
+ {
+ double max =
+ this->m_eDistance >= edge.m_eDistance ?
+ this->m_eDistance : edge.m_eDistance;
+ double tol = 0.3 * max;
+ double diff = abs(this->m_eDistance - edge.m_eDistance);
+ if (diff > tol)
+ cmp3 = false;
+ }
+ /*if (this->m_aRadius != edge.m_aRadius)
+ {
+ double max =
+ this->m_aRadius >= edge.m_aRadius ?
+ this->m_aRadius : edge.m_aRadius;
+ double tol = 0.3 * max;
+ double diff = abs(this->m_aRadius - edge.m_aRadius);
+ if (diff > tol)
+ cmp4 = false;
+ }
+ if (this->m_minRadius != edge.m_minRadius)
+ {
+ double max =
+ this->m_minRadius >= edge.m_minRadius ?
+ this->m_minRadius : edge.m_minRadius;
+ double tol = 0.3 * max;
+ double diff = abs(this->m_minRadius - edge.m_minRadius);
+ if (diff > tol)
+ cmp5 = false;
+ }
+ if (this->m_maxRadius != edge.m_maxRadius)
+ {
+ double max =
+ this->m_maxRadius >= edge.m_maxRadius ?
+ this->m_maxRadius : edge.m_maxRadius;
+ double tol = 0.3 * max;
+ double diff = abs(this->m_maxRadius - edge.m_maxRadius);
+ if (diff > tol)
+ cmp6 = false;
+ }*/
+ return (cmp1 && cmp2 && cmp3 && cmp4 && cmp5 && cmp6);
+ }
+
+bool Edge::operator!=(const Edge& edge)
+ {
+ return !(*this == edge);
+ }
+
+bool Edge::operator<(const Edge& edge)
+{
+ bool cmp1 = true;
+ bool cmp2 = true;
+ if (this->m_length != edge.m_length)
+ {
+ double max =
+ this->m_length >= edge.m_length ? this->m_length : edge.m_length;
+ double tol = 0.3 * max;
+ double diff = abs(this->m_length - edge.m_length);
+ if (diff > tol)
+ cmp2 = false;
+ }
+ if (this->m_eDistance != edge.m_eDistance)
+ {
+ double max =
+ this->m_eDistance >= edge.m_eDistance ?
+ this->m_eDistance : edge.m_eDistance;
+ double tol = 0.3 * max;
+ double diff = abs(this->m_eDistance - edge.m_eDistance);
+ if (diff > tol)
+ cmp2 = false;
+ }
+ if (cmp1 && cmp2)
+ return false;
+ if ((this->m_length < edge.m_length)
+ || (this->m_eDistance < edge.m_eDistance))
+ return true;
+ return false;
+
+}
+
+bool Edge::operator>(const Edge& edge)
+{
+ bool cmp1 = true;
+ bool cmp2 = true;
+ if (this->m_length != edge.m_length)
+ {
+ double max =
+ this->m_length >= edge.m_length ? this->m_length : edge.m_length;
+ double tol = 0.3 * max;
+ double diff = abs(this->m_length - edge.m_length);
+ if (diff > tol)
+ cmp2 = false;
+ }
+ if (this->m_eDistance != edge.m_eDistance)
+ {
+ double max =
+ this->m_eDistance >= edge.m_eDistance ?
+ this->m_eDistance : edge.m_eDistance;
+ double tol = 0.3 * max;
+ double diff = abs(this->m_eDistance - edge.m_eDistance);
+ if (diff > tol)
+ cmp2 = false;
+ }
+ if (cmp1 && cmp2)
+ return false;
+ if ((this->m_length > edge.m_length)
+ || (this->m_eDistance > edge.m_eDistance))
+ return true;
+ return false;
+}
+
+} /* namespace airways */
--- /dev/null
+/*
+ * airwaysEdge.h
+ *
+ * Created on: May 12, 2014
+ * Author: caceres
+ */
+
+#ifndef AIRWAYSEDGE_H_
+#define AIRWAYSEDGE_H_
+
+#include <cstddef>
+#include <cstdlib>
+#include <vector>
+#include <utility>
+#include <limits>
+
+#include "airwaysNode.h"
+#include "Quaternion.h"
+#include <appli/TempAirwaysAppli/AirwaysLib/AirwaysLib_Export.h>
+
+/**
+ * @brief Airways project namespace
+ */
+namespace airways
+{
+
+/*
+ * Pair of <Vec3,double> which means the position of the skeleton point
+ * between two nodes (including branch-points) and it radius (obtained from
+ * the Danielsson's distance map algorithm)
+ */
+typedef std::pair<Vec3, double> pair_posVox_rad;
+
+/*
+ * A vector stocking all the points between two nodes
+ */
+typedef std::vector<pair_posVox_rad> vec_pair_posVox_rad;
+
+class Node;
+
+/*
+ * This class represents an edge in a tree. Particularly, this edge represents the relation from child to father.
+ * This means that the edge source corresponds to the father and the target to the child. This implies that the root
+ * node has a NULL edge.
+ */
+class AirwaysLib_EXPORT Edge
+{
+public:
+ /*!
+ * Default constructor
+ */
+ Edge();
+
+ /*!
+ * Alternative constructor (copy)
+ * @param edge
+ */
+ Edge(Edge* edge);
+
+ /*!
+ * Destructor
+ */
+ virtual ~Edge();
+
+
+ //Getters
+
+ /*!
+ * This method returns the angle between two vectors using its parent as
+ * reference
+ * @return
+ */
+ int GetAngle() const;
+
+ /*!
+ * This method returns the Average radius between two branch-points
+ * @return
+ */
+ double GetARadius() const;
+
+ /*!
+ * This method returns the euclidean distance between two branch-points
+ * @return
+ */
+ double GetEDistance() const;
+
+ /*!
+ * This method returns the geodesic distance between two branch-points
+ * @return
+ */
+ double GetLength() const;
+
+ /*!
+ * This method returns the maximum radius between two branch-points
+ * @return
+ */
+ double GetMaxRadius() const;
+
+ /*!
+ * This method returns the minimum radius between two branch-points
+ * @return
+ */
+ double GetMinRadius() const;
+
+ /*!
+ * This method returns the source node of an edge
+ * @return
+ */
+ Node* GetSource() const;
+
+ /*!
+ * This method returns the target node of an edge
+ * @return
+ */
+ Node* GetTarget() const;
+
+ /*!
+ * This method returns the SKPairInfoVector of an edge
+ * @return
+ */
+ const vec_pair_posVox_rad& GetEdgeInfo() const;
+
+ /*!
+ * This method returns true if an edge is marked
+ * @return
+ */
+ bool IsMarked() const;
+
+ /**
+ * Method that returns if a voxel is influenced by this edge, i.e. if the distance
+ * to the closest voxel of the edge is smaller than the radius on the closest voxel.
+ * @param x,y,z are the voxel positions in coordinates x, y, and z respectively.
+ * @return true is the voxel is influenced, false otherwise
+ */
+ bool IsPointInfluencedByEdge(float point_x, float point_y, float point_z) const;
+
+ //Setters and Modifiers
+
+ /*!
+ * This method sets the angle between two vectors using its parent
+ * as reference
+ * @param angle
+ */
+ void SetAngle(const double& angle);
+
+ /*!
+ * This method sets the average radius between two nodes
+ * @param aRadius
+ */
+ void SetARadius(const double& aRadius);
+
+
+ /*!
+ * This method sets the Euclidean distance between two nodes
+ * @param eDistance
+ */
+ void SetEDistance(const double& eDistance);
+
+ /*!
+ * This method sets the geodesical length between two nodes
+ * @param length
+ */
+ void SetLength(const double& length);
+
+ /*!
+ * This method sets the maximum radius between two nodes
+ * @param maxRadius
+ */
+ void SetMaxRadius(const unsigned int& maxRadius);
+
+ /*!
+ * This method sets the minimum radius between two nodes
+ * @param minRadius
+ */
+ void SetMinRadius(const unsigned int& minRadius);
+
+ /*!
+ * This method sets the source node of an edge
+ * @param source
+ */
+ void SetSource(Node* source);
+
+ /*!
+ * This method sets the target node of an edge
+ * @param target
+ */
+ void SetTarget(Node* target);
+
+ /*!
+ * This method adds a pair <Coordinate, Radius> into skPInfoVector
+ * @param skPairInfo
+ */
+ void AddSkeletonPairInfo(const pair_posVox_rad& skPairInfo);
+
+ /*!
+ * This method sets the skPInfoVector
+ * @param skPInfoVector
+ */
+ void SetSkeletonPairVector(const vec_pair_posVox_rad& skPInfoVector);
+
+ /*!
+ * This method updates all the attributes of an edge using the information
+ * provided from m_skInfoPairVector attribute.
+ */
+ void UpdateEdgeInfo();
+
+ /**
+ * Method that compares two edges.
+ * @pre The actual and the edge given by parameter both have a parent edge
+ * @param edge is the edge to be compared with
+ * @return the evalation value
+ */
+ double CompareWith(Edge* edge);
+
+ /**
+ * Method that obtains the distance between two edges
+ * @param edge is the edge to be compared with
+ * @return the average distance from all points in the given edge
+ * to this edge. The distance from each in point in the edge to this
+ * edge is defined as the minimum distance to all the point in this
+ * edge.
+ */
+ float GetDistanceToEdge(Edge* edge);
+
+ float GetDistanceToTranslatedEdge(Edge* edge, Vec3 vector_translation);
+
+ /**
+ * Method that returns the distance point to point to the edge given by parameter.
+ * Additionally, the points that share the same closest point receive more weight in the distance.
+ * @param edge is the edge to be compared with
+ * @param vector_translation is the translation vector to translate the given edge
+ * @return the weigthed distance to the given edge
+ */
+ float GetDistanceWeigthedToTranslatedEdge(Edge* edge, Vec3 vector_translation);
+
+ /**
+ * Method that calculates the distance from edge to one point in the space
+ * @point is the point in the space
+ * @return the distance from edge to one point in the space
+ */
+ float GetDistanceToPoint(Vec3 point);
+
+ /**
+ * Method that returns the smallest distance of the edge voxels to the given voxel (represented as a vector)
+ * @param vector is the point to be compared
+ * @return the smallest distance of the edge voxels to the given voxel
+ */
+ float GetSmallestDistanceToPoint(Vec3 vector);
+
+ /**
+ * Method that returns the closest point in the edge to the given point (represented as a vector) and the distance between them (by parameter)
+ * @param vector is the point to be compared
+ * @param distance is the returned distance to the closes point
+ * @return the closest vector and the minimum distance of the point to the points in the edge by parameter.
+ */
+ Vec3 GetClosestPoint(Vec3 vector, float& distance);
+
+ /**
+ * Method that concatenates the actual edge to a superior one given by parameter
+ * @param superior is the superio edge. The target of the superior edge must be the source of the actual edge
+ * in terms of voxel position.
+ * @return the composed edge to the superior edge
+ */
+ Edge* ConcatenateToSuperior(Edge* superior);
+
+ //Operator Overload
+
+ Edge& operator=(Edge& edge);
+
+ const Edge& operator=(const Edge& edge);
+
+ bool operator==(const Edge& edge);
+
+ bool operator!=(const Edge& edge);
+
+ bool operator<(const Edge& edge);
+
+ bool operator>(const Edge& edge);
+
+protected:
+
+ unsigned int m_id;
+
+ /**
+ * Source node representing the father
+ */
+ Node* m_source;
+
+ /**
+ * Target node representing the child
+ */
+ Node* m_target;
+
+ bool m_mark; //!Mark attribute
+
+ double m_angle; //!Angle between two vectors having its parent as reference
+
+ double m_length; //!Geodesic distance
+
+ double m_eDistance; //!Euclidean distance
+
+ double m_aRadius; //!Average radius between two nodes
+
+ unsigned int m_minRadius; //!Minimum radius between two nodes
+
+ unsigned int m_maxRadius; //!Maximum radius between two nodes
+
+ /*!
+ * The m_vec_pair_posVox_rad is a vector containing all centerline points,
+ * in real coordinates of the image, and radius information of an edge.
+ * i.e., for each point of the skeleton between two nodes there is a pair
+ * <position real coordinates, radius>
+ */
+ vec_pair_posVox_rad m_vec_pair_posVox_rad;
+
+ friend class AirwaysTree; //!Friend class AirwaysTree
+
+ friend class Node; //!Friend class Node
+};
+
+} /* namespace airways */
+
+#endif /* AIRWAYSEDGE_H_ */
--- /dev/null
+/*
+ * airwaysNode.cxx
+ *
+ * Created on: May 12, 2014
+ * Author: Diego Cáceres
+ *
+ * Modified by: Alfredo Morales Pinzón
+ */
+
+#include <iostream>
+
+#include "airwaysNode.h"
+
+namespace airways
+{
+
+Node::Node() : m_id(0), m_children(), m_edge(NULL), m_coords(0, 0, 0), m_level(-1), m_mark(false)
+{
+}
+
+Node::Node(const Vec3& coord) : m_id(0), m_children(), m_edge(NULL), m_level(-1), m_mark(false)
+{
+ this->m_coords = coord;
+ this->m_level = 0;
+}
+
+Node::Node(Node* node) : m_id(0), m_children(), m_mark(false)
+{
+ this->m_coords = node->m_coords;
+ this->m_level = node->m_level;
+ this->m_edge = new Edge(node->m_edge);
+}
+
+Node::~Node()
+{
+}
+
+//Getters
+
+const std::vector<Node*>& Node::GetChildren() const
+{
+ return this->m_children;
+}
+
+const unsigned int Node::GetId() const
+{
+ return this->m_id;
+}
+
+const Vec3& Node::GetCoords() const
+{
+ return this->m_coords;
+}
+
+Edge* Node::GetEdge() const
+{
+ return this->m_edge;
+}
+
+Node* Node::GetFather() const
+{
+ return this->father;
+}
+
+unsigned int Node::GetLevel() const
+{
+ return this->m_level;
+}
+
+int Node::GetDepthById(unsigned int idNode) const
+{
+ if( this->m_id == idNode)
+ return 0;
+
+ vec_nodes::const_iterator it = this->m_children.begin();
+ for( ; it != this->m_children.end(); ++it)
+ {
+ if((*it)->HasNodeId(idNode))
+ return 1 + (*it)->GetDepthById(idNode);
+ }
+
+ std::cout << "GetDepthById - The idNode: " << idNode << "does not exist" << std::endl;
+
+ return -1;
+}
+
+unsigned int Node::GetWeight() const
+{
+ if(this->IsLeaf())
+ return 1;
+
+ unsigned int weight = 0;
+
+ vec_nodes::const_iterator it = this->m_children.begin();
+ for( ; it != this->m_children.end(); ++it)
+ weight += (*it)->GetWeight();
+
+ return 1 + weight;
+}
+
+unsigned int Node::GetHeight() const
+{
+ unsigned int heightChildren = 0;
+
+ vec_nodes::const_iterator it=this->m_children.begin();
+ for( ; it != this->m_children.end(); ++it)
+ {
+ unsigned int actualHeight = (*it)->GetHeight();
+ if(actualHeight > heightChildren)
+ heightChildren = actualHeight;
+ }
+
+ // Return my heigh plus my children height
+ return 1 + heightChildren;
+}
+
+unsigned int Node::GetNumberOfChildren() const
+{
+ return this->m_children.size();
+}
+
+unsigned int Node::GetNumberOfNonMarkedChildren( unsigned int depth ) const
+{
+ if(depth > 0)
+ {
+ unsigned int nonMarked = 0;
+
+ for(vec_nodes::const_iterator it = this->m_children.begin(); it != this->m_children.end(); ++it)
+ {
+ if( !(*it)->m_mark )
+ ++nonMarked;
+ nonMarked += (*it)->GetNumberOfNonMarkedChildren( depth - 1 );
+ }
+ return nonMarked;
+ }
+ return 0;
+}
+
+unsigned int Node::GetNumberOfLeafs() const
+{
+ if(this->IsLeaf())
+ return 1;
+
+ unsigned int leafs = 0;
+ vec_nodes::const_iterator it = this->m_children.begin();
+ for( ; it != this->m_children.end(); ++it)
+ leafs += (*it)->GetNumberOfLeafs();
+
+ return leafs;
+}
+
+Node* Node::GetNodeById(unsigned int nId)
+{
+ // Base case
+ if( m_id == nId)
+ return this;
+
+ Node* node_searched = NULL;
+
+ // Search in the children
+ vec_nodes::iterator it_children = m_children.begin();
+ for( ; it_children != m_children.end() && node_searched == NULL; ++it_children)
+ node_searched = (*it_children)->GetNodeById(nId);
+
+ return node_searched;
+}
+
+Node* Node::GetClosestBranchNodeToPoint(float point_x, float point_y, float point_z, float &distance)
+{
+ // Get the actual distance
+ distance = GetBranchDistanceToPoint(point_x, point_y, point_z);
+
+ // Make as closest "this" node
+ Node* node_closest = this;
+
+ // Iterate over the children
+ vec_nodes::iterator it_children = m_children.begin();
+ for( ; it_children != m_children.end(); ++it_children)
+ {
+ // Get actual child distance
+ float distance_actual = 0;
+ Node* bestNodeChild = (*it_children)->GetClosestBranchNodeToPoint(point_x, point_y, point_z, distance_actual);
+
+ // If the child distance is smaller then switch
+ if(distance_actual < distance)
+ {
+ node_closest = bestNodeChild;
+ distance = distance_actual;
+ }
+ }
+
+ return node_closest;
+}
+
+Node* Node::GetClosestNodeToPoint(Vec3 position, double& distance)
+{
+ // Calculate actual distance to point
+ Vec3 vector_diff = position - this->m_coords;
+ distance = vector_diff.Norm();
+
+ if(this->IsLeaf())
+ return this;
+
+ // Variables
+ double actualDistance = -1;
+ Node* node_Best = this;
+ Node* node_Actual = NULL;
+
+ // Check the children and get the best
+ for(vec_nodes::const_iterator it = this->m_children.begin(); it != this->m_children.end(); ++it)
+ {
+ node_Actual = (*it)->GetClosestNodeToPoint(position, actualDistance);
+
+ if(actualDistance < distance)
+ {
+ node_Best = node_Actual;
+ distance = actualDistance;
+ }
+ }
+
+ if(actualDistance == -1)
+ std::cout << "GetClosestNodeToPoint NOT FOUND!" << std::endl;
+
+ return node_Best;
+}
+
+float Node::GetBranchDistanceToPoint(int voxel_x, int voxel_y, int voxel_z)
+{
+ Vec3 voxel(voxel_x,voxel_y, voxel_z);
+ if(this->m_edge != NULL)
+ return this->m_edge->GetSmallestDistanceToPoint(voxel);
+ else
+ {
+ Vec3 vector_diff = this->m_coords - voxel;
+ return vector_diff.Norm();
+ }
+}
+
+bool Node::HasMinimumLevels(int levels)
+{
+ if(levels == 1)
+ return true;
+
+ bool minimum = false;
+
+ vec_nodes::iterator it=this->m_children.begin();
+ for( ; it != this->m_children.end() && !minimum; ++it)
+ if((*it)->HasMinimumLevels(levels-1))
+ minimum = true;
+
+ return minimum;
+}
+
+bool Node::IsLeaf() const
+{
+ return this->m_children.empty();
+}
+
+bool Node::IsMarked() const
+{
+ return this->m_mark;
+}
+
+bool Node::HasMarkedChildren() const
+{
+ if(this->IsLeaf())
+ return false;
+
+ // Check the children
+ vec_nodes::const_iterator it = this->m_children.begin();
+ for( ; it != this->m_children.end(); ++it)
+ if( (*it)->IsMarked() || (*it)->HasMarkedChildren())
+ return true;
+
+ return false;
+}
+
+bool Node::HasNode(Node* node_searched) const
+{
+ if(this->m_id == node_searched->m_id)
+ return true;
+
+ bool found = false;
+
+ vec_nodes::const_iterator it = this->m_children.begin();
+ for( ; it != this->m_children.end() && !found; ++it)
+ if( (*it)->HasNode(node_searched) )
+ found = true;
+
+ return found;
+}
+
+bool Node::HasNodeId(unsigned int id_node_searched) const
+{
+ // If this is the node then return true
+ if(this->m_id == id_node_searched)
+ return true;
+
+ // Search in the children
+ bool found = false;
+ vec_nodes::const_iterator it = this->m_children.begin();
+ for( ; it != this->m_children.end() && !found; ++it)
+ if( (*it)->HasNodeId(id_node_searched) )
+ found = true;
+
+ return found;
+}
+
+bool Node::IsPointInfluencedByNode(float point_x, float point_y, float point_z) const
+{
+ if(this->m_edge)
+ return m_edge->IsPointInfluencedByEdge(point_x, point_y, point_z);
+
+ return false;
+}
+
+void Node::GetBrancheNodesInfluencingPoint(float point_x, float point_y, float point_z, vec_nodes& vec_nodes_influence)
+{
+ if(this->IsPointInfluencedByNode(point_x, point_y, point_z))
+ vec_nodes_influence.push_back(this);
+
+ vec_nodes::const_iterator it = this->m_children.begin();
+ for( ; it != this->m_children.end(); ++it)
+ (*it)->GetBrancheNodesInfluencingPoint(point_x, point_y, point_z, vec_nodes_influence);
+}
+
+void Node::GetChildrenUptoDepth(unsigned int Q, vec_nodes& fathers)
+{
+ if(Q > 0)
+ {
+ vec_nodes::const_iterator it = this->m_children.begin();
+ for( ; it != this->m_children.end(); ++it)
+ {
+ fathers.push_back((*it));
+ (*it)->GetChildrenUptoDepth(Q-1, fathers);
+ }
+ /*// Add this node
+ fathers.push_back(this);
+
+ // If the depth is greater that 1 then add the children of this node
+ if(Q > 1)
+ {
+ vec_nodes::const_iterator it = this->m_children.begin();
+ for( ; it != this->m_children.end(); ++it)
+ (*it)->GetChildrenUptoDepth(Q-1, fathers);
+ }*/
+ }
+}
+
+void Node::CompareSubTreesOrkiszMorales(Node* nodeB, unsigned int Q, unsigned int F, vec_nodes& commonA, vec_nodes& commonB, vec_nodes& nonCommonA, vec_nodes& nonCommonB, map_id_vecnodes& map_A_to_B, map_id_vecnodes& map_B_to_A, std::vector< std::pair< pair_nodes, pair_nodes > >& vector_pair_edges)
+{
+ //std::cout << "---------------------------------------------------------" << std::endl;
+ //std::cout << "Comparing [A,B]: " << this->m_id << "," << nodeB->m_id << std::endl;
+
+ // Variables
+ double distanceAtoB = 0;
+ double distanceBtoA = 0;
+ Node* node_bestA = NULL;
+ Node* node_bestB = NULL;
+
+ // Algorithm
+ // The algorithm is recursive it is implemented the node class. The idea is to overlap the initial nodes of the subtrees
+ // and find the best branch correspondence for the first level of the trees.
+ // Hypothesis: It is supposed that the roots of the trees are common, so the algorithm does not run in the root nodes
+ // 1. Translate one of the subtree so that the initial nodes overlap.
+ // 2. While one of the trees has non-marked child
+ // 2.1 Compare each child to all the branches is the other tree.
+ // 2.2 Select the pair, child and branch in the other tree, that has the lowest distance function.
+ // 2.3 Add the pair and the distance to the final result.
+ // 2.4 Mark the nodes
+ // 2.5 Find and mark the omitted nodes in the similar branch. If the selected branch has more that 2 nodes it means that intermediary nodes do not have correspondence.
+ // 2.6 Run the process (First step) using the found pair of nodes.
+
+ //1. Translate one of the subtree so that the initial nodes overlap.
+ // The translation must be applied later
+
+ Vec3 vector_trans_B_to_A = this->m_coords - nodeB->m_coords;
+
+ // ------------------------
+ // ------------------------
+ // Making the matching faster
+ std::multimap<double, pair_nodes> map_dist_pairNodes_A_to_B = this->CalculateDistanceToAllBranches_FatherAndFamily(nodeB, Q, F);
+ std::multimap<double, pair_nodes> map_dist_pairNodes_B_to_A = nodeB->CalculateDistanceToAllBranches_FatherAndFamily(this, Q, F);
+
+ // ------------------------
+ // ------------------------
+
+ while( this->GetNumberOfNonMarkedChildren( Q ) > 0 && nodeB->GetNumberOfNonMarkedChildren( Q ) > 0 )
+ {
+ // 2.1 Compare each child to all the branches is the other tree.
+
+ // 2.1.1 Get best translated branch from A to B
+ //std::pair<pair_nodes,double> pair_A_to_B = this->GetClosestBranch_FatherAndFamily(nodeB);
+ // 2.1.2 Get best translated branch from B to A
+ //std::pair<pair_nodes,double> pair_B_to_A = nodeB->GetClosestBranch_FatherAndFamily(this);
+
+ // ------------------------
+ // ------------------------
+ // Making the matching faster
+ // Get the next best pair
+ std::pair<double,pair_nodes> pair_A_to_B = GetPairWithClosestDistance_FirstNodeNonMarked(map_dist_pairNodes_A_to_B);
+ std::pair<double,pair_nodes> pair_B_to_A = GetPairWithClosestDistance_FirstNodeNonMarked(map_dist_pairNodes_B_to_A);
+
+ if( ! pair_A_to_B.second.first || ! pair_A_to_B.second.second )
+ std::cout << "There is no closest-non-marked branch!" << std::endl;
+
+ distanceAtoB = pair_A_to_B.first;
+ distanceBtoA = pair_B_to_A.first;
+
+ //std::map<double, pair_nodes> map_dist_pairNodes_A_to_B = this->CalculateDistanceToAllBranches_FatherAndFamily(nodeB);
+ //std::map<double, pair_nodes> map_dist_pairNodes_B_to_A = nodeB->CalculateDistanceToAllBranches_FatherAndFamily(this);
+ // ------------------------
+ // ------------------------
+
+
+ // Assign the distances
+ //distanceAtoB = pair_A_to_B.second;
+ //distanceBtoA = pair_B_to_A.second;
+
+ // 2.2 Select the pair, child and branch in the other tree, that has the lowest distance function.
+ if(distanceAtoB < distanceBtoA)
+ {
+ //node_bestA = pair_A_to_B.first.first;
+ //node_bestB = pair_A_to_B.first.second;
+ node_bestA = pair_A_to_B.second.first;
+ node_bestB = pair_A_to_B.second.second;
+
+ //std::cout << "Best Pair [A',B]: " << node_bestA->m_id << "," << node_bestB->m_id << std::endl;
+
+ // 2.2.1 Check that there are no topological problems before adding the pair
+ // Search all the nodes in the tree A linked to the found best node in B
+ vec_nodes nodesA_linkedToBestB = map_B_to_A[node_bestB->m_id];
+
+ // Look for topological problems
+ bool topo_problem = false;
+ for(vec_nodes::iterator it = nodesA_linkedToBestB.begin(); it != nodesA_linkedToBestB.end() && !topo_problem; ++it)
+ {
+ if( ! ( (*it)->HasNode( node_bestA ) || node_bestA->HasNode( (*it) ) ) )
+ {
+ topo_problem = true;
+ //std::cout << "Topological problem adding node:" << node_bestA->m_id << std::endl;
+ }
+ }
+
+ if(! topo_problem)
+ {
+ // 2.3 Add the pair and the distance to the final result.
+ commonA.push_back(node_bestA);
+ commonB.push_back(node_bestB);
+ map_A_to_B[node_bestA->m_id].push_back(node_bestB);
+ map_B_to_A[node_bestB->m_id].push_back(node_bestA);
+
+ // Add the edge link
+ pair_nodes pair_edgeA(this, node_bestA);
+ pair_nodes pair_edgeB(nodeB, node_bestB);
+ std::pair< pair_nodes, pair_nodes > pair_edge(pair_edgeA, pair_edgeB);
+ vector_pair_edges.push_back(pair_edge);
+
+
+ // 2.4 Mark the nodes
+ node_bestA->Mark();
+ node_bestB->Mark();
+
+ // 2.5 Find and mark the omitted nodes in the similar branch. If the selected branch has more that 2 nodes it means that intermediary nodes do not have correspondence.
+ Node* node_nextA = this->GetNextNodeInPathToNode(node_bestA);
+ //bool markedA = node_nextA->MarkNodesUntilNode(node_bestA);
+ //if(!markedA)
+ // std::cout << "NodeA not found for marking [A]:" << node_bestA->m_id << std::endl;
+
+ Node* node_nextB = nodeB->GetNextNodeInPathToNode(node_bestB);
+ //bool markedB = node_nextB->MarkNodesUntilNode(node_bestB);
+ //if(!markedB)
+ // std::cout << "NodeB not found for marking [B]:" << node_bestB->m_id << std::endl;
+
+ // 2.6 Run the process (First step) using the found pair of nodes.
+ node_bestA->CompareSubTreesOrkiszMorales(node_bestB, Q, F, commonA, commonB, nonCommonA, nonCommonB, map_A_to_B, map_B_to_A, vector_pair_edges);
+ }
+ else
+ {
+ // 2.4 Mark the wrong node
+ node_bestA->Mark();
+ nonCommonA.push_back(node_bestA);
+ }
+ }
+ else
+ {
+ //node_bestA = pair_B_to_A.first.second;
+ //node_bestB = pair_B_to_A.first.first;
+ node_bestA = pair_B_to_A.second.second;
+ node_bestB = pair_B_to_A.second.first;
+
+ //std::cout << "Best Pair [A,B']: " << node_bestA->m_id << "," << node_bestB->m_id << std::endl;
+
+ // 2.29 Check that there are no topological problems before adding the pair
+ // 2.291 Search all the nodes in the A tree linked to the found best node in B
+ vec_nodes nodesB_linkedToBestA = map_A_to_B[node_bestA->m_id];
+
+ // 2.292 Look for topological problems
+ bool topo_problem = false;
+ for(vec_nodes::iterator it = nodesB_linkedToBestA.begin(); it != nodesB_linkedToBestA.end() && !topo_problem; ++it)
+ {
+ if( ! ( (*it)->HasNode(node_bestB) || node_bestB->HasNode( (*it) ) ) )
+ {
+ topo_problem = true;
+ //std::cout << "Topological problem adding node:" << node_bestB->m_id << std::endl;
+ }
+ }
+
+ if(! topo_problem)
+ {
+ // 2.3 Add the pair and the distance to the final result.
+ commonA.push_back(node_bestA);
+ commonB.push_back(node_bestB);
+ map_A_to_B[node_bestA->m_id].push_back(node_bestB);
+ map_B_to_A[node_bestB->m_id].push_back(node_bestA);
+
+ // Add the edge link
+ pair_nodes pair_edgeA(this, node_bestA);
+ pair_nodes pair_edgeB(nodeB, node_bestB);
+ std::pair< pair_nodes, pair_nodes > pair_edge(pair_edgeA, pair_edgeB);
+ vector_pair_edges.push_back(pair_edge);
+
+ // 2.4 Mark the nodes
+ node_bestA->Mark();
+ node_bestB->Mark();
+
+ // 2.5 Find and mark the omitted nodes in the similar branch. If the selected branch has more that 2 nodes it means that intermediary nodes do not have correspondence.
+ Node* node_nextA = this->GetNextNodeInPathToNode(node_bestA);
+ //bool markedA = node_nextA->MarkNodesUntilNode(node_bestA);
+ //if(!markedA)
+ // std::cout << "NodeA not found for marking [A]:" << node_bestA->m_id << std::endl;
+
+ Node* node_nextB = nodeB->GetNextNodeInPathToNode(node_bestB);
+ //bool markedB = node_nextB->MarkNodesUntilNode(node_bestB);
+ //if(!markedB)
+ // std::cout << "NodeB not found for marking [B]:" << node_bestB->m_id << std::endl;
+
+ // 2.6 Run the process (First step) using the found pair of nodes.
+ //node_bestB->CompareSubTreesOrkiszMorales(node_bestA, commonB, commonA, nonCommonB, nonCommonA, map_B_to_A, map_A_to_B, vector_pair_edges);
+ node_bestA->CompareSubTreesOrkiszMorales(node_bestB, Q, F, commonA, commonB, nonCommonA, nonCommonB, map_A_to_B, map_B_to_A, vector_pair_edges);
+ }
+ else
+ {
+ // 2.4 Mark the wrong node
+ node_bestB->Mark();
+ nonCommonB.push_back(node_bestB);
+ }
+ }
+
+ }
+ //std::cout << "OUT Comparing [A,B]: " << this->m_id << "," << nodeB->m_id << std::endl;
+ //std::cout << "---------------------------------------------------------" << std::endl;
+}
+
+std::pair< std::pair<Node*, Node*>, double> Node::GetBestBranches_ActualTranslatedOneLevelBranches_To_AllPossibleBranches(Node* nodeB)
+{
+ // Variables
+ double distance = -1;
+ double distance_bestFinal = -1;
+ Node* node_A_bestFinal = NULL;
+ Node* node_B_bestFinal = NULL;
+ bool first = true;
+
+ // Get the translation vector for root ovelapping
+ Vec3 vector_trans_A_to_B = nodeB->m_coords - this->m_coords;
+
+ // Iterate over my children
+ for(vec_nodes::iterator it_A = this->m_children.begin(); it_A != this->m_children.end(); ++it_A)
+ {
+ if( !(*it_A)->IsMarked() )
+ {
+ // Get the position of the actual child of this node
+ Vec3 position_actual = (*it_A)->m_coords;
+
+ // Translate the position
+ Vec3 position_translated = position_actual + vector_trans_A_to_B;
+
+ // Find the "closest" node in the subtree nodeB using a distance function
+ // a. Distance to ending point
+ //Node* bestNode_actual = nodeB->GetClosestRelativeToPoint(position_translated, distance);
+ // b. Distance point to point
+ Node* bestNode_actual = nodeB->GetClosestRelativeToBranch( (*it_A)->m_edge, vector_trans_A_to_B, distance);
+
+ if(first)
+ {
+ first = false;
+ node_A_bestFinal = (*it_A);
+ node_B_bestFinal = bestNode_actual;
+ distance_bestFinal = distance;
+ }
+ else if(distance < distance_bestFinal)
+ {
+ node_A_bestFinal = (*it_A);
+ node_B_bestFinal = bestNode_actual;
+ distance_bestFinal = distance;
+ }
+ }
+ }
+
+ if(first)
+ std::cout << "GetBestTranslatedBranchFromOtherTree NOT FOUND!" << std::endl;
+
+ pair_nodes pair_best(node_A_bestFinal, node_B_bestFinal);
+ std::pair<pair_nodes,double> best_pair_score(pair_best, distance_bestFinal);
+ return best_pair_score;
+}
+
+std::multimap< double, std::pair<Node*, Node*> > Node::CalculateDistanceToAllBranches_FatherAndFamily(Node* node_gradfather, unsigned int Q, unsigned int F)
+{
+ // "this" is a grandfather node
+ // Variables
+ std::multimap< double, pair_nodes > map_dist_pairNodeNode;
+ double distance_actual = -1;
+
+ // Get the translation vector for root overlapping
+ Vec3 vector_trans_A_to_B_grandfathers = node_gradfather->m_coords - this->m_coords;
+
+ // Get all the possible "fathers" based on Q
+ vec_nodes fathers;
+ this->GetChildrenUptoDepth(Q, fathers);
+ //fathers = this->m_children;
+
+ // Iterate over my children which actually are fathers
+
+ //vec_nodes::iterator it_father = this->m_children.begin();
+ vec_nodes::iterator it_father = fathers.begin();
+ //for( ; it_father != this->m_children.end(); ++it_father)
+ for( ; it_father != fathers.end(); ++it_father)
+ {
+ if( ! (*it_father)->IsMarked() )
+ {
+
+ Node* node_actualBest = node_gradfather->GetClosestRelativeFatherAndFamily( this, (*it_father), F, vector_trans_A_to_B_grandfathers, distance_actual);
+
+ pair_nodes pair_actual((*it_father), node_actualBest);
+
+ std::pair<double,pair_nodes> actual_dist_pair(distance_actual,pair_actual);
+ map_dist_pairNodeNode.insert(actual_dist_pair);
+ }
+ }
+
+ return map_dist_pairNodeNode;
+}
+
+
+/*std::pair< std::pair<Node*, Node*>, double> Node::GetClosestBranch_FatherAndFamily(Node* node_gradfather)
+{
+ // "this" is a grandfather node
+ // Variables
+ double distance_actual = -1;
+ double distance_bestFinal = -1;
+ Node* node_A_bestFinal = NULL;
+ Node* node_B_bestFinal = NULL;
+ bool first = true;
+
+ // Get the translation vector for root overlapping
+ Vec3 vector_trans_A_to_B_grandfathers = node_gradfather->m_coords - this->m_coords;
+
+ // Iterate over my children which actually are fathers
+ vec_nodes::iterator it_father = this->m_children.begin();
+ for( ; it_father != this->m_children.end(); ++it_father)
+ {
+ if( !(*it_father)->IsMarked() )
+ {
+ // Get the position of the actual child of this node
+ //Vec3 position_actualFather = (*it_father)->m_coords;
+
+ // Find the "closest" node in the subtree nodeB using a distance function
+ // a. Distance to ending point
+ // Translate the position
+ //Vec3 position_translated = position_actualFather + vector_trans_A_to_B_grandfathers;
+ //Node* bestNode_actual = nodeB->GetClosestRelativeToPoint(position_translated, distance);
+ // b. Distance point to point
+ Node* node_actualBest = node_gradfather->GetClosestRelativeFatherAndFamily( (*it_father), vector_trans_A_to_B_grandfathers, distance_actual);
+
+ if(first)
+ {
+ first = false;
+ node_A_bestFinal = (*it_father);
+ node_B_bestFinal = node_actualBest;
+ distance_bestFinal = distance_actual;
+ }
+ else if(distance_actual < distance_bestFinal)
+ {
+ node_A_bestFinal = (*it_father);
+ node_B_bestFinal = node_actualBest;
+ distance_bestFinal = distance_actual;
+ }
+ }
+ }
+
+ if(first)
+ std::cout << "GetBestTranslatedBranchFromOtherTree NOT FOUND!" << std::endl;
+
+ pair_nodes pair_best(node_A_bestFinal, node_B_bestFinal);
+ std::pair<pair_nodes,double> best_pair_score(pair_best, distance_bestFinal);
+ return best_pair_score;
+}
+*/
+
+Node* Node::GetClosestRelativeToPoint(Vec3 position, double& distance)
+{
+ if(this->IsLeaf())
+ {
+ distance = -1;
+ return NULL;
+ }
+
+ distance = -1;
+ bool first = true;
+ double actualDistance = -1;
+ Node* node_Best = NULL;
+ Node* node_Actual = NULL;
+
+ for(vec_nodes::const_iterator it = this->m_children.begin(); it != this->m_children.end(); ++it)
+ {
+ node_Actual = (*it)->GetClosestNodeToPoint(position, actualDistance);
+
+ if(first)
+ {
+ first = false;
+ node_Best = node_Actual;
+ distance = actualDistance;
+ }
+ else if(actualDistance < distance)
+ {
+ node_Best = node_Actual;
+ distance = actualDistance;
+ }
+ }
+
+ if(first)
+ std::cout << "GetClosestRelativeToPoint NOT FOUND!" << std::endl;
+
+ if(!node_Best)
+ std::cout << "GetClosestRelativeToPoint NULL!" << std::endl;
+
+ return node_Best;
+}
+
+Node* Node::GetClosestRelativeToBranch(Edge* branch, Vec3 vector_translation, double& distance)
+{
+ if(this->IsLeaf())
+ {
+ distance = -1;
+ //std::cout << "GetClosestRelativeToBranch -- This node is a leaf" << std::endl;
+ return NULL;
+ }
+
+ distance = -1;
+ bool first = true;
+ double actualDistance = -1;
+ Node* node_Best = NULL;
+ Node* node_Actual = NULL;
+
+ for(vec_nodes::const_iterator it = this->m_children.begin(); it != this->m_children.end(); ++it)
+ {
+ node_Actual = (*it)->GetClosestBranchToTranslatedBranch(branch,vector_translation, actualDistance, NULL);
+
+ if(first)
+ {
+ first = false;
+ node_Best = node_Actual;
+ distance = actualDistance;
+ }
+ else if(actualDistance < distance)
+ {
+ node_Best = node_Actual;
+ distance = actualDistance;
+ }
+ }
+
+ if(first)
+ std::cout << "GetClosestRelativeToPoint NOT FOUND!" << std::endl;
+
+ if(!node_Best)
+ std::cout << "GetClosestRelativeToPoint NULL!" << std::endl;
+
+ return node_Best;
+}
+
+Node* Node::GetClosestRelativeFatherAndFamily(Node* node_grandfather, Node* node_father, unsigned int F, Vec3 vector_translation, double& distance)
+{
+ // "this" is a grandfather node
+
+ distance = -1;
+ if(this->IsLeaf())
+ return NULL;
+
+ bool first = true;
+ double actualDistance = -1;
+ Node* node_Best = NULL;
+ Node* node_Actual = NULL;
+
+ // Iterate over the children of this granfather which means
+ // iterate over the fathers
+ vec_nodes::const_iterator it_father = this->m_children.begin();
+
+ for( ; it_father != this->m_children.end(); ++it_father)
+ {
+ node_Actual = (*it_father)->GetClosestFatherAndFamilyToTranslatedFatherAndFamily(node_grandfather, node_father, F, vector_translation, actualDistance, NULL);
+
+ if(first)
+ {
+ first = false;
+ node_Best = node_Actual;
+ distance = actualDistance;
+ }
+ else if(actualDistance < distance)
+ {
+ node_Best = node_Actual;
+ distance = actualDistance;
+ }
+ }
+
+ if(first)
+ std::cout << "GetClosestRelativeToPoint NOT FOUND!" << std::endl;
+
+ if(!node_Best)
+ std::cout << "GetClosestRelativeToPoint NULL!" << std::endl;
+
+ return node_Best;
+}
+
+Node* Node::GetClosestBranchToTranslatedBranch(Edge* branch, Vec3 vector_translation, double& distance, Edge* edgeSuperior)
+{
+ Node* node_best = this;
+ Node* node_Actual = NULL;
+ Edge* edge_composed = NULL;
+
+ // Joint the superior edge with this edge
+ edge_composed = this->ConcatenateEdgeToSuperiorEdge(edgeSuperior);
+
+ distance = edge_composed->GetDistanceToTranslatedEdge(branch, vector_translation);
+ distance += branch->GetDistanceToTranslatedEdge(edge_composed, -vector_translation);
+ //distance = edge_composed->GetDistanceWeigthedToTranslatedEdge(branch, vector_translation);
+ //distance += branch->GetDistanceWeigthedToTranslatedEdge(edge_composed, -vector_translation);
+ //std::cout << "DW[thisId,toBranchTargetId,Dist]:[" << this->m_id << "," << branch->GetTarget()->GetId() << "," << distance << "]";
+
+ // Lance the recursion and the the minumum distance and node
+ double distance_actual = -1;
+
+ // Check the children and get the best
+ for(vec_nodes::iterator it = this->m_children.begin(); it != this->m_children.end(); ++it)
+ {
+ Edge* edge_forChild = new Edge(edge_composed);
+ node_Actual = (*it)->GetClosestBranchToTranslatedBranch(branch, vector_translation, distance_actual, edge_forChild);
+
+ if(distance_actual < distance)
+ {
+ node_best = node_Actual;
+ distance = distance_actual;
+ }
+ delete(edge_forChild);
+ }
+
+ return node_best;
+}
+
+Node* Node::GetClosestFatherAndFamilyToTranslatedFatherAndFamily(Node* node_grandfather, Node* node_father, unsigned int F, Vec3 vector_translation, double& distance, Edge* edgeSuperior)
+{
+ // "this" is a father node_father
+ Node* node_best = this; // Best matching node is "this" in the beginnig.
+ Node* node_Actual = NULL;
+ Edge* edge_composed = NULL;
+
+ // Joint the superior edge with this edge
+ edge_composed = this->ConcatenateEdgeToSuperiorEdge(edgeSuperior);
+
+ //distance = edge_composed->GetDistanceToTranslatedEdge(branch, vector_translation);
+ //distance += branch->GetDistanceToTranslatedEdge(edge_composed, -vector_translation);
+
+ //Edge* branch_father = node_father->m_edge;
+ Edge* branch_father = node_father->GetEdgeToAncestor(node_grandfather, NULL);
+ if( branch_father == NULL )
+ std::cout << "No branch father" << std::endl;
+
+
+ distance = edge_composed->GetDistanceWeigthedToTranslatedEdge(branch_father, vector_translation);
+ distance += branch_father->GetDistanceWeigthedToTranslatedEdge(edge_composed, -vector_translation);
+
+ // Find the family distance from the node_father family to this family
+ double distance_family_A_to_B = 0;
+ // Add the distances to the family, up to one generation by the moment
+ Vec3 vector_translation_child_A_to_B = this->m_coords - node_father->m_coords; // Vector that align the fathers
+
+ // Get the family up to depth F
+ vec_nodes family;
+ node_father->GetChildrenUptoDepth(F, family);
+ //vec_nodes family = node_father->m_children;
+
+ vec_nodes::iterator it_child_other = family.begin();
+ for(; it_child_other != family.end(); ++it_child_other)
+ {
+ // Create the "child" edge
+ Edge* childEdge = (*it_child_other)->GetEdgeToAncestor(node_father, NULL);
+ //Edge* childEdge = (*it_child_other)->m_edge;
+
+ double distance_child = 0;
+ if(this->IsLeaf())
+ {
+ // The distance of the each point edge to the root point because the fathers are aligned
+ distance_child = childEdge->GetDistanceToPoint(this->m_coords);
+ }
+ else
+ {
+ this->GetClosestRelativeToBranch(childEdge,vector_translation_child_A_to_B, distance_child);
+ }
+ distance_family_A_to_B += distance_child;
+ }
+
+ distance += distance_family_A_to_B;
+
+ // Find the family distance from this family to node_father family
+ double distance_family_B_to_A = 0;
+
+ // Add the distances to the family, up to one generation by the moment
+ Vec3 vector_translation_child_B_to_A = node_father->m_coords - this->m_coords; // Vector that align the fathers
+
+ // Get the family up to depth F
+ vec_nodes my_family;
+ this->GetChildrenUptoDepth(F, my_family);
+ //vec_nodes my_family = this->m_children;
+
+ vec_nodes::iterator it_child_my = my_family.begin();
+ for(; it_child_my != my_family.end(); ++it_child_my)
+ {
+ // Create the "child" edge
+ Edge* my_childEdge = (*it_child_my)->GetEdgeToAncestor(this, NULL);
+ //Edge* my_childEdge = (*it_child_my)->m_edge;
+
+ double distance_child = 0;
+ if(node_father->IsLeaf())
+ {
+ // The distance of the each point edge to the root point because the fathers are aligned
+ distance_child = my_childEdge->GetDistanceToPoint(node_father->m_coords);
+ }
+ else
+ {
+ node_father->GetClosestRelativeToBranch(my_childEdge, vector_translation_child_B_to_A, distance_child);
+ }
+ distance_family_B_to_A += distance_child;
+ }
+
+ distance += distance_family_B_to_A;
+
+ //std::cout << "DW[thisId,toBranchTargetId,Dist]:[" << this->m_id << "," << branch->GetTarget()->GetId() << "," << distance << "]";
+
+ // Lance the recursion and the minimum distance and node_father
+ double distance_actual = -1;
+
+ // Check the children and get the best
+ vec_nodes::iterator it = this->m_children.begin();
+ for( ; it != this->m_children.end(); ++it)
+ {
+ Edge* edge_forChild = new Edge(edge_composed);
+ node_Actual = (*it)->GetClosestFatherAndFamilyToTranslatedFatherAndFamily(node_grandfather, node_father, F, vector_translation, distance_actual, edge_forChild);
+
+ if(distance_actual < distance)
+ {
+ node_best = node_Actual;
+ distance = distance_actual;
+ }
+ delete(edge_forChild);
+ }
+
+ return node_best;
+}
+
+const Node* Node::GetClosestCommonAscendingNode() const
+{
+ Node* node_closestCommonAscendingNode = NULL;
+
+ if(this->m_edge)
+ {
+ if(this->m_edge->GetSource()->IsMarked() || this->m_edge->GetSource()->HasMarkedChildren())
+ node_closestCommonAscendingNode = this->m_edge->GetSource();
+ else
+ return this->m_edge->GetSource()->GetClosestCommonAscendingNode();
+ }
+
+ return node_closestCommonAscendingNode;
+}
+
+bool Node::MarkNodesUntilNode(Node* nodeB)
+{
+ this->Mark();
+
+ if(this->m_id == nodeB->m_id)
+ return true;
+
+ bool found = false;
+ bool actualFound = false;
+ for(vec_nodes::const_iterator it = this->m_children.begin(); it != this->m_children.end(); ++it)
+ {
+ actualFound = (*it)->MarkNodesUntilNode(nodeB);
+
+ if(actualFound)
+ found = true;
+ }
+
+ return found;
+}
+
+Node* Node::GetNextNodeInPathToNode(Node* node_searched)
+{
+ Node* nextNode = NULL;
+ for(vec_nodes::const_iterator it = this->m_children.begin(); it != this->m_children.end() && !nextNode; ++it)
+ {
+ if( (*it)->HasNode(node_searched) )
+ nextNode = (*it);
+ }
+
+ return nextNode;
+}
+
+void Node::GetPathToNode(Node* node_end, vec_nodes& vector_pathNodes)
+{
+ if(!this->HasNode(node_end))
+ return;
+
+ vector_pathNodes.push_back(this);
+ for(vec_nodes::const_iterator it = this->m_children.begin(); it != this->m_children.end(); ++it)
+ (*it)->GetPathToNode(node_end, vector_pathNodes);
+}
+
+void Node::MarkPathFromNodeToNode(Node* node_begin, Node* node_end)
+{
+ if(this->m_id == node_begin->m_id)
+ {
+ if(this->HasNode(node_end))
+ this->MarkPathNodesToNode(node_end);
+ }
+ else
+ for(vec_nodes::const_iterator it = this->m_children.begin(); it != this->m_children.end(); ++it)
+ (*it)->MarkPathFromNodeToNode(node_begin, node_end);
+}
+
+void Node::MarkPathNodesToNode(Node* node_end)
+{
+ if(this->HasNode(node_end))
+ this->Mark();
+
+ for(vec_nodes::const_iterator it = this->m_children.begin(); it != this->m_children.end(); ++it)
+ (*it)->MarkPathNodesToNode(node_end);
+}
+
+Node* Node::GetDetachedRootCopy()
+{
+ Node* node_copy = new Node(this);
+
+ // Change properties for root node
+ Node* root_detached = new Node();
+ root_detached->m_coords = node_copy->m_edge->GetSource()->GetCoords();
+ root_detached->m_edge = NULL;
+ vec_nodes vector_new;
+ vector_new.push_back(node_copy);
+ root_detached->m_children = vector_new;
+
+ // Change properties for the copy node
+ node_copy->m_edge->SetSource(root_detached);
+
+ return root_detached;
+}
+
+bool Node::CompareNodes(Node* nodeTreeB)
+{
+ if (this == NULL && nodeTreeB == NULL)
+ return true;
+ else if (this == NULL || nodeTreeB == NULL)
+ return false;
+ vec_nodes::iterator it = this->m_children.begin();
+ for (; it != this->m_children.end(); ++it)
+ {
+ if ((*it)->m_mark)
+ continue;
+ bool suc = false;
+ vec_nodes::iterator it2 = nodeTreeB->m_children.begin();
+ for (; it2 != nodeTreeB->m_children.end(); ++it2)
+ suc = (*it)->CompareNodes(*it2);
+ if (!suc)
+ return false;
+ } //for
+ if (*this != *nodeTreeB)
+ return false;
+ this->m_mark = true;
+ nodeTreeB->m_mark = true;
+ return true;
+}
+
+Edge* Node::GetComposedEdge(unsigned int idNode)
+{
+ // If this is the node, then a copy edge (to the parent) is returned
+ if(this->m_id == idNode)
+ {
+ //std::cout << "This is the searched idNode: " << this->m_id << std::endl;
+ Edge* nEdge = new Edge(this->m_edge);
+ return nEdge;
+ }
+
+
+ Edge* composedEdge = NULL;
+
+ for(std::vector<Node*>::iterator it = this->m_children.begin(); it != this->m_children.end(); ++it)
+ {
+ // If one of my children has the composed node, then my edge must be added
+ // and the composed edge is returned
+ Edge* childEdge = (*it)->GetComposedEdge(idNode);
+ if(childEdge)
+ {
+ //std::cout << "My id: " << this->m_id << ", Looking for: " << idNode << std::endl;
+ return AddChildEdgeInformation(childEdge);
+ }
+ }
+
+ return composedEdge;
+}
+
+Edge* Node::AddChildEdgeInformation(Edge* nEdge)
+{
+ // If this is the root then the edge to be compose is null, so the same edge must be returned
+ if(this->m_edge == NULL)
+ {
+ //std::cout << "This is the root node" << std::endl;
+ return nEdge;
+ }
+
+ // New information vector
+ vec_pair_posVox_rad vector_newInformation;
+ vec_pair_posVox_rad actualEdgeInformation = this->m_edge->GetEdgeInfo();
+ vec_pair_posVox_rad composedInformation = nEdge->GetEdgeInfo();
+
+ // Check that last and initial information, from this and the other edge respectively, are equal.
+ if(actualEdgeInformation[actualEdgeInformation.size()-1].first != composedInformation[0].first)
+ {
+
+ std::cout << "actualEdgeInformation[0;end]: [" << actualEdgeInformation[0].first << ";" << actualEdgeInformation[actualEdgeInformation.size()-1].first << "], Points: " << actualEdgeInformation.size() << std::endl;
+ std::cout << "ActualEdge[source;target]: [" << this->m_edge->GetSource()->GetCoords() << "];[" << this->m_edge->GetTarget()->GetCoords() << "]" << std::endl;
+ std::cout << "-- -- -- " << std::endl;
+ std::cout << "composedInformation[0;end]: [" << composedInformation[0].first << ";" << composedInformation[composedInformation.size()-1].first << "], Points: " << composedInformation.size() << std::endl;
+ std::cout << "AddChildEdge[source;target]: [" << nEdge->GetSource()->GetCoords() << "];[" << nEdge->GetTarget()->GetCoords() << "]" << std::endl;
+ std::cout << "End and initial information are not equal! : " << actualEdgeInformation[actualEdgeInformation.size()-1].first << " , " << composedInformation[0].first << std::endl;
+
+
+ std::cout << "Print composed edge information" << std::endl;
+ std::cout << "Composed [source, target] " << nEdge->GetSource()->GetCoords() << "];[" << nEdge->GetTarget()->GetCoords() << "]" << std::endl;
+ std::cout << "Composed information using iterator:" << std::endl;
+
+ for(vec_pair_posVox_rad::const_iterator it_composed = nEdge->GetEdgeInfo().begin(); it_composed != nEdge->GetEdgeInfo().end(); ++it_composed)
+ std::cout << "[" << (*it_composed).first << "], ";
+ }
+
+ // Change the source
+ nEdge->SetSource(this->m_edge->GetSource());
+
+ // Change properties
+ nEdge->SetLength(nEdge->GetLength()+this->m_edge->GetLength());
+ nEdge->SetEDistance(nEdge->GetEDistance()+this->m_edge->GetEDistance());
+
+ // Create the new information vector and calculate the new radii attributes
+ for(vec_pair_posVox_rad::iterator it_actual = actualEdgeInformation.begin(); it_actual != actualEdgeInformation.end(); ++it_actual)
+ vector_newInformation.push_back((*it_actual));
+
+ bool connectionVoxel = true;
+ for(vec_pair_posVox_rad::iterator it_composed = composedInformation.begin(); it_composed != composedInformation.end(); ++it_composed)
+ {
+ if(connectionVoxel)
+ connectionVoxel = false;
+ else
+ vector_newInformation.push_back((*it_composed));
+ }
+
+ // Set the new information
+ nEdge->SetSkeletonPairVector(vector_newInformation);
+
+ return nEdge;
+}
+
+Edge* Node::GetEdgeToAncestor(Node* node_ancestor, Edge* edge)
+{
+ Edge* edge_composed = NULL;
+
+ // Compose the edge
+ if(edge == NULL)
+ edge_composed = new Edge(this->m_edge);
+ else
+ edge_composed = edge->ConcatenateToSuperior(new Edge(this->m_edge));
+
+ // Check if the father is the ancestor
+ if(this->m_edge->GetSource()->GetId() == node_ancestor->GetId())
+ return edge_composed;
+
+ return this->m_edge->GetSource()->GetEdgeToAncestor(node_ancestor, edge_composed);
+}
+
+Edge* Node::ConcatenateEdgeToSuperiorEdge(Edge* superiorEdge)
+{
+ if(superiorEdge == NULL)
+ return new Edge(this->m_edge);
+
+ // Check concatenation condition in the nodes
+ if(superiorEdge->GetTarget()->GetCoords() != this->m_edge->GetSource()->GetCoords())
+ {
+ std::cout << "Node::ConcatenateEdgeToSuperiorEdge - superiorEdge.target != actualEdge.source. Actual idNode:" << this->m_id << std::endl;
+ std::cout << "[superiorEdge.target;actualEdge.source]: [" << superiorEdge->GetTarget()->GetCoords() << ";" << this->m_edge->GetSource()->GetCoords() << "]" << std::endl;
+ return NULL;
+ }
+
+ // Check concatenation condition in the edge information
+ vec_pair_posVox_rad vectorInfo_thisEdge = this->m_edge->GetEdgeInfo(); // Actual edge information
+ vec_pair_posVox_rad vectorInfo_superiorEdge = superiorEdge->GetEdgeInfo(); // Superior edge information
+
+ // Check that last superior edge position is equal to first initial actual edge position
+ if(vectorInfo_superiorEdge[vectorInfo_superiorEdge.size()-1].first != vectorInfo_thisEdge[0].first)
+ {
+ std::cout << "Node::ConcatenateEdgeToSuperiorEdge - superiorEdge.info[end] != actualEdge.info[begin]. Actual idNode:" << this->m_id << std::endl;
+ return NULL;
+ }
+
+ // Change the superior edge target, the superior length, and the euclidian distance
+ superiorEdge->SetTarget(this->m_edge->GetTarget());
+ superiorEdge->SetLength(superiorEdge->GetLength()+this->m_edge->GetLength()-1);
+ superiorEdge->SetEDistance(superiorEdge->GetEDistance()+this->m_edge->GetEDistance());
+
+ // Add the position information
+ vec_pair_posVox_rad::iterator it_actualInfo = vectorInfo_thisEdge.begin();
+ ++it_actualInfo; // Skip the first position because is it shared between both edges
+ for(; it_actualInfo != vectorInfo_thisEdge.end(); ++it_actualInfo)
+ vectorInfo_superiorEdge.push_back((*it_actualInfo));
+
+ // Set the new information
+ superiorEdge->SetSkeletonPairVector(vectorInfo_superiorEdge);
+
+ return superiorEdge;
+}
+
+bool Node::FindSimilarEdgeByAccumulation(Node* nodeB, Node* nodeBAcc)
+{
+ bool firstTime = false;
+ if (nodeBAcc == NULL)
+ {
+ nodeBAcc = new Node();
+ nodeBAcc->m_edge = new Edge();
+ firstTime = true;
+ } //if
+ nodeBAcc->m_coords = nodeB->m_coords;
+ nodeBAcc->m_level = this->m_level;
+ Edge* edgeBAcc = nodeBAcc->m_edge;
+ Edge* edgeB = nodeB->m_edge;
+ edgeBAcc->m_angle = edgeB->m_angle;
+ edgeBAcc->m_eDistance = edgeB->m_eDistance;
+ if (!firstTime)
+ {
+ edgeBAcc->m_length += edgeB->m_length;
+ //edgeBAcc->m_aRadius = (edgeB->m_aRadius + edgeBAcc->m_aRadius) / 2.0;
+ edgeBAcc->m_maxRadius = edgeB->m_maxRadius > edgeBAcc->m_maxRadius ? edgeB->m_maxRadius : edgeBAcc->m_maxRadius;
+ edgeBAcc->m_minRadius = edgeB->m_minRadius < edgeBAcc->m_minRadius ? edgeB->m_minRadius : edgeBAcc->m_minRadius;
+ } //fi
+ else
+ {
+ edgeBAcc->m_length = edgeB->m_length;
+ edgeBAcc->m_aRadius = edgeB->m_aRadius;
+ edgeBAcc->m_maxRadius = edgeB->m_maxRadius;
+ edgeBAcc->m_minRadius = edgeB->m_minRadius;
+ } //else
+
+ if (*this == *nodeBAcc)
+ return true;
+ Node* aux = NULL;
+ bool suc = false;
+ vec_nodes::iterator it = nodeB->m_children.begin();
+ for (; it != nodeB->m_children.end(); ++it)
+ {
+ aux = *it;
+ if (this->FindSimilarEdgeByAccumulation(aux, nodeBAcc))
+ {
+ suc = true;
+ break;
+ } //if
+ } //for
+ if (suc)
+ {
+ nodeB = aux;
+ return true;
+ }//if
+ return false;
+}
+
+bool Node::FindMarked(Node* result)
+{
+ vec_nodes::iterator it = this->m_children.begin();
+ for (; it != this->m_children.end(); ++it)
+ if ((*it)->m_mark)
+ {
+ result = *it;
+ return true;
+ } //if
+ else if ((*it)->FindMarked(result))
+ return true;
+ result = NULL;
+ return false;
+}
+
+bool Node::MarkPath(Node* nodeB)
+{
+ if (this == NULL || nodeB == NULL)
+ return false;
+ //comparing addresses
+ if (this == nodeB)
+ {
+ nodeB->m_mark = true;
+ return true;
+ } //if
+ vec_nodes::iterator it = this->m_children.begin();
+ for (; it != this->m_children.end(); ++it)
+ if ((*it)->MarkPath(nodeB))
+ {
+ (*it)->m_mark = true;
+ return true;
+ } //if
+ return false;
+
+}
+
+unsigned int Node::GetCost(Node* nodeB)
+{
+ if (this == NULL || nodeB == NULL)
+ return 0;
+ unsigned int maxCost = 0;
+ if (*this == *nodeB)
+ {
+ vec_nodes::iterator it = this->m_children.begin();
+ for (; it != this->m_children.end(); ++it)
+ {
+ vec_nodes::iterator it2 = nodeB->m_children.begin();
+ for (; it2 != nodeB->m_children.end(); ++it2)
+ {
+ unsigned int currentCost = (*it)->GetCost(*it2);
+ if (currentCost > maxCost)
+ maxCost = currentCost;
+ } //for
+ } //for
+ maxCost += 1;
+ } //if
+ return maxCost;
+}
+
+
+//Setters and Modifiers
+
+void Node::SetId(unsigned int nM_id)
+{
+ m_id = nM_id;
+}
+
+
+void Node::AddChild(Node* node)
+{
+ this->m_children.push_back(node);
+}
+
+void Node::SetChildren(std::vector<Node*>& children)
+{
+ this->m_children = children;
+}
+
+void Node::SetFather(Node* nFather)
+{
+ father = nFather;
+}
+
+void Node::SetEdge(Edge* edge)
+{
+ this->m_edge = edge;
+}
+
+void Node::SetCoords(const Vec3& coords)
+{
+ this->m_coords = coords;
+}
+
+void Node::SetLevel(const int& level)
+{
+ this->m_level = level;
+}
+
+bool Node::DeleteNodeById(unsigned int nId)
+{
+ if(this->IsLeaf())
+ return false;
+
+ bool deleted = false;
+
+ vec_nodes::iterator it = this->m_children.begin();
+ while(it != this->m_children.end() && !deleted)
+ {
+ if((*it)->m_id == nId)
+ {
+ it = m_children.erase(it);
+ deleted = true;
+ }
+ else
+ ++it;
+ }
+
+ if(!deleted)
+ {
+ for(vec_nodes::const_iterator it = this->m_children.begin(); it != this->m_children.end() && !deleted; ++it)
+ {
+ deleted = (*it)->DeleteNodeById(nId);
+ }
+ }
+
+ return deleted;
+}
+
+void Node::MergeNodes(Node* nodeB)
+{
+ if (nodeB == NULL)
+ return;
+ this->m_coords = nodeB->m_coords;
+ Edge* edgeA = this->m_edge;
+ Edge* edgeB = nodeB->m_edge;
+ edgeA->m_angle = edgeB->m_angle;
+ edgeA->m_eDistance = edgeB->m_eDistance;
+ edgeA->m_length += edgeB->m_length;
+ edgeA->m_aRadius = (edgeB->m_aRadius + edgeA->m_aRadius) / 2.0;
+ edgeA->m_maxRadius = edgeB->m_maxRadius > edgeA->m_maxRadius ? edgeB->m_maxRadius : edgeA->m_maxRadius;
+ edgeA->m_minRadius = edgeB->m_minRadius < edgeA->m_minRadius ? edgeB->m_minRadius : edgeA->m_minRadius;
+ edgeA->m_vec_pair_posVox_rad.insert(edgeA->m_vec_pair_posVox_rad.end(), edgeB->m_vec_pair_posVox_rad.begin(), edgeB->m_vec_pair_posVox_rad.end());
+}
+
+void Node::Mark()
+{
+ this->m_mark = true;
+}
+
+void Node::UnMark()
+{
+ this->m_mark = false;
+}
+
+void Node::UpdateLevels(unsigned int level)
+{
+ this->m_level = level;
+ level++;
+ vec_nodes::iterator it = this->m_children.begin();
+ for (; it != this->m_children.end(); ++it)
+ (*it)->UpdateLevels(level);
+}
+
+void Node::printUpToLevel(int level)
+{
+ // Print myself
+ cout << "level" << level << " " << m_coords << " || " << endl;
+ if(level == 1)
+ return;
+
+ // Print my children
+ for (vec_nodes::iterator it = this->m_children.begin(); it != this->m_children.end(); ++it)
+ {
+ (*it)->printUpToLevel(level-1);
+ }
+}
+
+void Node::printNodeAndChildrenIds( )
+{
+ std::cout << "------------------------------------------" << std::endl;
+ std::cout << "Id: " << this->m_id << std::endl;
+
+ // Print children
+ for(vec_nodes::const_iterator it = this->m_children.begin(); it != this->m_children.end(); ++it)
+ std::cout << (*it)->m_id << ";";
+
+ // Recursion
+ for(vec_nodes::const_iterator it = this->m_children.begin(); it != this->m_children.end(); ++it)
+ (*it)->printNodeAndChildrenIds( );
+}
+
+void Node::printChildrenIds( )
+{
+ std::cout << "------------------------------------------" << std::endl;
+ std::cout << "Id: " << this->m_id << std::endl;
+
+ // Print children
+ for(vec_nodes::const_iterator it = this->m_children.begin(); it != this->m_children.end(); ++it)
+ std::cout << (*it)->m_id << ";";
+}
+
+void Node::createQuaternionsUpToLevel(int level, int requiredLevel)
+{
+ //cout << "Level:" << requiredLevel-level <<endl;
+
+ if(level >= 2 && this->HasMinimumLevels(2))
+ {
+ // 1. Take the vector to each son (SP - The vector that points from the son to the parent [P-S]).
+ // 2. For each son take each vector from son to grandson (SG - the vector that points from the son to the grandson [G-S])
+ // 3. For each combination SP-SG calculate the quaternion
+ for(vec_nodes::iterator it = this->m_children.begin(); it != this->m_children.end(); ++it)
+ {
+ // Create vector from son to parent (SP = [P - S])
+ Vec3 sp(this->m_coords - (*it)->m_coords);
+
+ // Get the mean radius for the parent-son branch
+ double radMeanSP=(*it)->m_edge->GetARadius();
+
+ // Get distance son to parent
+ double distSonParent = (*it)->m_edge->GetEDistance();
+
+ // Get the vectors from son to each grandson (SG = [G - S] )
+ for(vec_nodes::iterator itps = (*it)->m_children.begin(); itps != (*it)->m_children.end(); ++itps)
+ {
+ // Build the sg vector and get the quaternion
+ Vec3 sg((*itps)->m_coords - (*it)->m_coords);
+ Quaternion q(sp, sg);
+
+ // Get distance son to grandson
+ double distSonGrandson = (*itps)->m_edge->GetEDistance();
+
+ // Get the mean radius for the son-grandson branch
+ double radMeanSG=(*itps)->m_edge->GetARadius();
+
+ //cout << "Vectors:" << ps << sg <<endl;
+ //cout << "Q:" << q << ", Distance=" << (*itps)->m_edge->GetEDistance() << endl << endl;
+ // Printing
+ // Level | PositionSon | SPvector | SGvector | Quaternion | distanceSonParent | distanceSonGrandSon | meanRadiusSP | meanRadiusSG
+ cout << requiredLevel-level << " " << (*it)->m_coords << " " << sp << " " << sg << " " << q << " " << distSonParent << " " << distSonGrandson << " " << radMeanSP << " " << radMeanSG << endl;
+ }
+ }
+
+ for(vec_nodes::iterator it = this->m_children.begin(); it != this->m_children.end(); ++it)
+ {
+ (*it)->createQuaternionsUpToLevel(level-1, requiredLevel);
+ }
+ }
+}
+
+std::pair<double, std::pair<Node*, Node*> > Node::GetPairWithClosestDistance_FirstNodeNonMarked(std::multimap<double, pair_nodes> map_dist_pairNodes)
+{
+ // Variables
+ bool pairFound = false;
+ std::pair<Node*, Node*> pairNodes_null (NULL, NULL);
+ std::pair<double, std::pair<Node*, Node*> > pair_closest_dist_pairNodes (-1, pairNodes_null);
+
+ std::multimap<double, pair_nodes>::iterator it = map_dist_pairNodes.begin();
+
+ for( ; it != map_dist_pairNodes.end() && !pairFound; ++it)
+ {
+ if( ! (*it).second.first->IsMarked() )
+ {
+ pair_closest_dist_pairNodes.first=(*it).first;
+ pair_closest_dist_pairNodes.second=(*it).second;
+ pairFound = true;
+ }
+ }
+
+ return pair_closest_dist_pairNodes;
+}
+
+void Node::getStasticsBifurcations(int* p)
+{
+ int numBif = this->m_children.size();
+ if(numBif>0 && numBif<=6)
+ p[numBif-1]=p[numBif-1]+1;
+
+ for(vec_nodes::iterator it = this->m_children.begin(); it != this->m_children.end(); ++it)
+ {
+ (*it)->getStasticsBifurcations(p);
+ }
+}
+
+//Operator Overloads
+
+Node &Node::operator=(Node & node)
+{
+ this->m_children = node.m_children;
+ this->m_edge = node.m_edge;
+ this->m_coords = node.m_coords;
+ this->m_level = node.m_level;
+ this->m_mark = node.m_mark;
+ return *this;
+}
+
+const Node&Node::operator=(const Node& node)
+{
+ this->m_children = node.m_children;
+ this->m_edge = node.m_edge;
+ this->m_coords = node.m_coords;
+ this->m_level = node.m_level;
+ this->m_mark = node.m_mark;
+ return *this;
+}
+
+bool Node::operator ==(const Node& node)
+{
+ bool comp = false;
+ if (this->m_edge == NULL && this->m_edge != NULL)
+ return false;
+ else if (this->m_edge != NULL && this->m_edge == NULL)
+ return false;
+ else if (this->m_edge == NULL && this->m_edge == NULL)
+ comp = true;
+ else
+ comp = ((*this->m_edge) == (*node.m_edge));
+ bool ret = (/*(this->m_coords == node.m_coords)
+ && (this->m_children.size() == node.m_children.size())
+ && */(this->m_level == node.m_level) && comp);
+ return ret;
+}
+
+bool Node::operator !=(const Node& node)
+{
+ bool comp = false;
+ if (this->m_edge == NULL && this->m_edge != NULL)
+ return false;
+ else if (this->m_edge != NULL && this->m_edge == NULL)
+ return false;
+ else if (this->m_edge == NULL && this->m_edge == NULL)
+ comp = true;
+ else
+ comp = ((*this->m_edge) == (*node.m_edge));
+ return (/*(this->m_coords != node.m_coords)
+ && (this->m_children.size() != node.m_children.size())
+ && */(this->m_level != node.m_level) && !comp);
+}
+
+bool Node::operator <(const Node& node)
+{
+ bool comp = false;
+ if (this->m_edge == NULL && this->m_edge != NULL)
+ return false;
+ else if (this->m_edge != NULL && this->m_edge == NULL)
+ return false;
+ else if (this->m_edge == NULL && this->m_edge == NULL)
+ comp = false;
+ else
+ comp = ((*this->m_edge) < (*node.m_edge));
+ return comp;
+}
+
+bool Node::operator >(const Node& node)
+{
+ bool comp = false;
+ if (this->m_edge == NULL && this->m_edge != NULL)
+ return false;
+ else if (this->m_edge != NULL && this->m_edge == NULL)
+ return false;
+ else if (this->m_edge == NULL && this->m_edge == NULL)
+ comp = false;
+ else
+ comp = ((*this->m_edge) > (*node.m_edge));
+ return comp;
+}
+
+} /* namespace airways */
--- /dev/null
+/*
+ * airwaysNode.h
+ *
+ * Created on: May 12, 2014
+ * Author: caceres
+ *
+ * Modified by: Alfredo Morales Pinzón
+ */
+
+#ifndef AIRWAYSNODE_H_
+#define AIRWAYSNODE_H_
+
+#include <vector>
+#include <map>
+
+#include "../MathLib/vec3.h"
+#include "airwaysEdge.h"
+#include "Quaternion.h"
+#include <appli/TempAirwaysAppli/AirwaysLib/AirwaysLib_Export.h>
+
+// Namespaces
+using namespace std;
+
+namespace airways
+{
+
+class AirwaysTree;
+class Edge;
+
+/*
+ * This class represents a node in a tree.
+ */
+class AirwaysLib_EXPORT Node
+{
+
+public:
+
+ // Pair of nodes
+ typedef std::pair<Node*, Node*> pair_nodes;
+
+ // Vector of nodes
+ typedef std::vector<Node*> vec_nodes;
+
+ // Map to link an id node to a list of nodes
+ typedef std::map< unsigned int, vec_nodes > map_id_vecnodes;
+
+ /*!
+ * This is the default constructor
+ */
+ Node();
+
+ /*!
+ * This is the alternative constructor
+ * @param coord
+ */
+ Node(const Vec3& coord);
+
+ /*!
+ * This is the alternative constructor (copy)
+ * @param node
+ */
+ Node(Node* node);
+
+ /*!
+ * This is the destructor
+ */
+ virtual ~Node();
+
+ //Getters
+
+ /*!
+ * This method returns the vector of the node children (const)
+ * @return
+ */
+ const std::vector<Node*>& GetChildren() const;
+
+ /**
+ * Method that returns the node id
+ * @return the id of the node
+ */
+ const unsigned int GetId() const;
+
+ /*!
+ * This method returns the spatial position of a node
+ * @return
+ */
+ const Vec3& GetCoords() const;
+
+ /*
+ * Method that returns the edge between the current node and its parent
+ * @return the edge
+ */
+ Edge* GetEdge() const;
+
+ /*
+ * Method that returns the father
+ * @return the father
+ */
+ Node* GetFather() const;
+
+ /*
+ * This method returns the level of the node
+ * @return the level of the node
+ */
+ unsigned int GetLevel() const;
+
+ /**
+ * Method that returns the depth of the node whose id is given by parameter
+ * @pre the node exist in the sub-tree
+ * @param idNode is the search node id
+ * @return the depth of the node up to the actual node
+ */
+ int GetDepthById(unsigned int idNode) const;
+
+ /**
+ * Method that returns the weight of the tree
+ * @return the weight of the tree which is the number of nodes in the tree.
+ */
+ unsigned int GetWeight() const;
+
+ /**
+ * Method that returns the height of the tree
+ * @return the height of the tree
+ */
+ unsigned int GetHeight() const;
+
+ /*!
+ * This method returns the size of the children vector
+ * @return the number of children
+ */
+ unsigned int GetNumberOfChildren() const;
+
+ /**
+ * Method that returns the number of children non-marked
+ * @param depth is the depth up to which the non-marked child are taken into account
+ * @return the number of children non-marked
+ */
+ unsigned int GetNumberOfNonMarkedChildren( unsigned int depth ) const;
+
+ /**
+ * Method that returns the number of leafs in the node and its children
+ * @return the number of leafs in the in the node and its children
+ */
+ unsigned int GetNumberOfLeafs( ) const;
+
+ /**
+ * This method returns the node that has the requested id
+ * @param nId is the requested id
+ * @return is the node that has the requested id, NULL if the node is not found.
+ */
+ Node* GetNodeById(unsigned int nId);
+
+ /**
+ * Method that returns the closest node to the givel voxel
+ * The closes node is the one whose path (skeleton) to the father has the closest point to the given point.
+ * @param voxel_x x componente of the voxel
+ * @param voxel_y y componente of the voxel
+ * @param voxel_z z componente of the voxel
+ * @param distance return the distance to the closest node
+ * @return the closest node to the given voxel
+ */
+ Node* GetClosestBranchNodeToPoint(float point_x, float point_y, float point_z, float &distance);
+
+ /**
+ * Method that returns the closest node to the position given by parameter
+ * @param position is the position to be compares with
+ * @param distance is the minimum distance found
+ * return the closest node to the position given by parameter
+ */
+ Node* GetClosestNodeToPoint(Vec3 position, double& distance);
+
+ /**
+ * Method that return the closest distance of the node skeleton
+ * @param voxel_x x componente of the voxel
+ * @param voxel_y y componente of the voxel
+ * @param voxel_z z componente of the voxel
+ * @return the closest distance of the node skeleton
+ */
+ float GetBranchDistanceToPoint(int voxel_x, int voxel_y, int voxel_z);
+
+ /**
+ * Method that returns if the tree has a minimum of levels
+ * @return true if the tree has the minimum of levels, false otherwise
+ */
+ bool HasMinimumLevels(int levels);
+
+
+ /*!
+ * This method returns true if the node is a leaf
+ * @return
+ */
+ bool IsLeaf() const;
+
+ /*!
+ * This method returns true if the node is marked
+ * @return
+ */
+ bool IsMarked() const;
+
+ /**
+ * Method that tells if a node has marked children
+ * @return true if the node has marked children, false otherwise
+ */
+ bool HasMarkedChildren() const;
+
+ /**
+ * Method that returns if a node is in the descending set of nodes
+ * @return true is the given node is in the descending set of nodes, false otherwise
+ */
+ bool HasNode(Node* node_searched) const;
+
+ /**
+ * Method that return if a node is in the tree
+ * @param id_node_searched id of the searched node
+ * @return true is the node is in the tree, false otherwise
+ */
+ bool HasNodeId(unsigned int id_node_searched) const;
+
+ /**
+ * Method that returns if a voxel is influenced by a node, i.e. if the distance
+ * to the closest voxel of the edge is smaller than the radius on the closest voxel.
+ * @param x,y,z are the voxel positions in coordinates x, y, and z respectively.
+ * @return true is the voxel is influenced, false otherwise.
+ */
+ bool IsPointInfluencedByNode(float point_x, float point_y, float point_z) const;
+
+ /**
+ * Method that returns the list of nodes that influence the given point
+ * @param point_x, point_y, point_z are the position of the point in real coordinates x, y, and z respectively.
+ * @param vec_nodes_influence vector to save the nodes
+ */
+ void GetBrancheNodesInfluencingPoint(float point_x, float point_y, float point_z, vec_nodes& vec_nodes_influence);
+
+ /**
+ * Method that adds to a vector all the children with a maximum depth Q from actual node.
+ * If Q = 0, no nodes are added.
+ * @param Q the depth to find the children
+ * @param fathers vector to return the children
+ */
+ void GetChildrenUptoDepth(unsigned int Q, vec_nodes& fathers);
+
+ /**
+ * This method compares the subtrees that start in this node and one given by parameter.
+ * @param nodeB is the stating node of the subtree to be compared with
+ * @param Q is the depth to select "fathers" nodes
+ * @param F is the depth to select "family" nodes
+ * @param commonA vector to save the common nodes for this subtree
+ * @param commonA vector to save the common nodes for the subtree given by parameter
+ * @param nonCommonA vector to save the non-common nodes for this subtree
+ * @param nonCommonB vector to save the non-common nodes for the subtree given by parameter
+ */
+ void CompareSubTreesOrkiszMorales(Node* nodeB, unsigned int Q, unsigned int F, vec_nodes& commonA, vec_nodes& commonB, vec_nodes& nonCommonA, vec_nodes& nonCommonB, map_id_vecnodes& map_A_to_B, map_id_vecnodes& map_B_to_A, std::vector< std::pair< pair_nodes, pair_nodes > >& vector_pair_edges);
+
+ /**
+ * Method that returns the pair of branches (one node for each tree) that has the minimum distance after root overlapping.
+ * Valid branches for the actual node are the one level branches those the branches up to this child (height = 2).
+ * Root overlapping means the translation of branch so that root overlap.
+ * @pre: The actual branch is not marked.
+ * @param nodeB is the subtree root node to be compared with
+ * @return the pair of branches (one node for each tree) that has the minimum distance after root overlapping
+ */
+ std::pair< std::pair<Node*, Node*>, double> GetBestBranches_ActualTranslatedOneLevelBranches_To_AllPossibleBranches(Node* nodeB);
+
+ /**
+ * Method that returns all the distances to all the non-marked branches, i.e., the end node of the branch is non-marked,
+ * inside a multimap (key, value) structure where the value elementent corresponds to the compared branches.
+ * @param node_gradfather is the grandfather node to be compared
+ * @param Q is the depth to select "fathers"
+ * @param F is the depth to select "family" nodes
+ * @return multimap where the key is the distance between branches and the value is the pair of compared branches end-nodes.
+ */
+ std::multimap< double, std::pair<Node*, Node*> > CalculateDistanceToAllBranches_FatherAndFamily(Node* node_gradfather, unsigned int Q, unsigned int F);
+
+ /**
+ * Method is executed by grandfather nodes. It returns the closest pair of nodes,
+ * one must be my child and the other one of the relatives given node.
+ * The closest pair is the one that has the minimal distance between the children and the families.
+ * The distance is the sum of fathers distances and family distances.
+ * @param nodeB is the node to be compared with
+ * @return the closest pair of nodes, one must be my child and the other on of the relatives given node.
+ */
+ //std::pair< std::pair<Node*, Node*>, double> GetClosestBranch_FatherAndFamily(Node* nodeB);
+
+ /**
+ * Method that returns the closest relative node, all the nodes but this, whose position is the closest to the point given by parameter
+ * @param position is the point to be compared with
+ * @param distance is the minimum distance found in the relatives
+ * @return the closest relative node, all the nodes but this, whose position is the closest to the point given by parameter
+ */
+ Node* GetClosestRelativeToPoint(Vec3 position, double& distance);
+
+
+ Node* GetClosestRelativeToBranch(Edge* branch, Vec3 vector_translation, double& distance);
+
+ /**
+ * Method that returns the closest node to the given by node.
+ * Closest involve the distance between fathers and family.
+ * @param node_grandfather is the father of the "node_father" which can be different from the direct father.
+ * this is the case when the father distance is the algorithm, given by Q, is greater than 1.
+ * @param node_father is the node to be compared with
+ * @param F is the depth to select "family" nodes
+ * @param vector_translation is the translation vector to translate the "father" branch
+ * @param distance is the distance to be returned by parameter
+ * @return the closest node to the given by node.
+ */
+ Node* GetClosestRelativeFatherAndFamily(Node* node_grandfather, Node* node_father, unsigned int F, Vec3 vector_translation, double& distance);
+
+ Node* GetClosestBranchToTranslatedBranch(Edge* branch, Vec3 vector_translation, double& distance, Edge* edgeSuperior);
+
+ /**
+ * This method returns the closest father and family, in this subtree, that is the closest to the subtree given by parameter
+ * The returned node represents the father.
+ * @param node_grandfather is the father of the "node_father" which can be different from the direct father.
+ * this is the case when the father distance is the algorithm, given by Q, is greater than 1.
+ * @param node_father is the node father to be compared with
+ * @param F is the depth to select "family" nodes
+ * @param vector_translation is the translation vector to translation the father given by parameter
+ * @param distance is the closest distances to be returned by parameter
+ * @param edgeSuperior is a superior edge where the actual edge must be attached
+ * @returns the closest father and family, in this subtree, that is the closest to the subtree given by parameter
+ */
+ Node* GetClosestFatherAndFamilyToTranslatedFatherAndFamily(Node* node_grandfather, Node* node_father, unsigned int F, Vec3 vector_translation, double& distance, Edge* edgeSuperior);
+
+ /**
+ * Method that return the closest ascending node that is marked as common.
+ * A common node may common if:
+ * (1) it is marked or (2) one of their descendants is marked
+ * return the closest ascending common node. NULL is there is no such node.
+ */
+ const Node* GetClosestCommonAscendingNode() const;
+
+ /**
+ * Method that marks all nodes between this node and the one given by parameter. If the searched node is not found then all nodes are marked.
+ * @param nodeB is the searched node
+ * @return true is the node was found, false otherwise
+ */
+ bool MarkNodesUntilNode(Node* nodeB);
+
+ Node* GetNextNodeInPathToNode(Node* node_searched);
+
+ /**
+ * Method that returns the path from the actual node to the searched node
+ * @param node_end is the searched node
+ * @param vector_pathNodes vector_pathNodesis the vector to return the path
+ */
+ void GetPathToNode(Node* node_end, vec_nodes& vector_pathNodes);
+
+ void MarkPathFromNodeToNode(Node* node_begin, Node* node_end);
+
+ void MarkPathNodesToNode(Node* node_end);
+
+ /**
+ * Method that returns a root with a copy of this node as single child
+ * @return a root with a copy of this node as single child
+ */
+ Node* GetDetachedRootCopy();
+
+ /*!
+ * This method return true if two nodes are equals
+ * @param nodeTreeB
+ * @return
+ */
+ bool CompareNodes(Node* nodeTreeB);
+
+ /**
+ * Method that returns a new composed edge from this node to the searched node.
+ * @param idNode is the searched id node
+ * @return the edge to the searched node. Null if the node is not found.
+ */
+ Edge* GetComposedEdge(unsigned int idNode);
+
+ /**
+ * Method that composes the actual edge with the edge given by parameter
+ * @param nEdge edge to be composed with
+ * @return final composed edge
+ */
+ Edge* AddChildEdgeInformation(Edge* nEdge);
+
+ /**
+ * Method that returns a "composed edge" up to the ancestor
+ * @pre the ancesto node must be a real ancestor
+ * @param node_ancestor is the ancestor node to with the edge must be composed
+ * @param edge is the actual composed edge. NULL if the fist iteration
+ * @return the composed edge up to the ancesto node
+ */
+ Edge* GetEdgeToAncestor(Node* node_ancestor, Edge* edge);
+
+ /**
+ * Method that concatenates the actual edge to a superior edge given by parameter.
+ * As the given edge is superior then the actual edge must be concatenate at the end of the superior edge,
+ * this means that superiorEdge.target must equal to the actualEdge.source.
+ * @param superiorEdge is the superior edge where the actual edge must be added
+ * @return - a copy of the actual edge if the superior edge is NULL
+ * - if superiorEdge.target == actualEdge.source the the superior edge with the actual edge concatenated
+ * - NULL otherwise (superiorEdge.target != actualEdge.source)
+ */
+ Edge* ConcatenateEdgeToSuperiorEdge(Edge* superiorEdge);
+
+ /*!
+ * This method returns true if the accumulation of a set of nodes are equals
+ * (in tolerance) to the compared node (this). This method tries to find
+ * edges in the case of edge losses.
+ * When this method returns true, NodeB pointer is replaced to the node (target)
+ * found by accumulation (See the report to understand what this means).
+ * @param nodeB - Current node to compare
+ * @param nodeBAcc - Accumulator Node (Optional)
+ * @return
+ */
+ bool FindSimilarEdgeByAccumulation(Node* nodeB, Node* nodeBAcc = NULL);
+
+ /*!
+ * This method returns true if there is a marked node along the hierarchy
+ * @param result
+ * @return
+ */
+ bool FindMarked(Node* result);
+
+ /*!
+ * This method marks the path between two nodes and returns true if a path
+ * has found.
+ * @param nodeB
+ * @return
+ */
+ bool MarkPath(Node* nodeB);
+
+ /*!
+ * This method returns the cost between two nodes.
+ * @param nodeB
+ * @return
+ */
+ unsigned int GetCost(Node* nodeB);
+
+ //Setters and Modifiers
+
+ /*
+ * Method that sets the id of the node
+ * @param nM_id is the new id
+ */
+ void SetId(unsigned int nM_id);
+
+ /*!
+ * This method adds a node in the children vector
+ * @param node is the node to be added
+ */
+ void AddChild(Node* node);
+
+ /*!
+ * This method sets the children vector
+ * @param children
+ */
+ void SetChildren(std::vector<Node*>& children);
+
+ /**
+ * Method that sets the father
+ * @param nFather is the father
+ */
+ void SetFather(Node* nFather);
+
+ /*!
+ * This method sets an edge between the current node and its parent
+ * @param edge
+ */
+ void SetEdge(Edge* edge);
+
+ /*!
+ * This method sets the coordinates of the node
+ * @param coords
+ */
+ void SetCoords(const Vec3& coords);
+
+ /**
+ * Method that sets the level of the node
+ * @param level is the level of the node
+ */
+ void SetLevel(const int& level);
+
+ /**
+ * Method that deleted a node and all his children, from the tree
+ * @return true if the node was deleted (found), false otherwise
+ */
+ bool DeleteNodeById(unsigned int nId);
+
+ /*!
+ * This method merges the nodeB into this
+ * @param nodeB
+ */
+ void MergeNodes(Node* nodeB);
+
+ /*!
+ * This method marks a node
+ */
+ void Mark();
+
+ /*!
+ * This method unmarks a node
+ */
+ void UnMark();
+
+ /*!
+ * This method update all levels along the hierarchy
+ * @param level
+ */
+ void UpdateLevels(unsigned int level = 0);
+
+
+ /**
+ * Method that prints the tree up to a given level
+ * @param level where the print must finish
+ */
+ void printUpToLevel(int level);
+
+ void printNodeAndChildrenIds( );
+
+ void printChildrenIds( );
+
+ /**
+ * Method that creates the quaternions of the airways up the maxLevel
+ * @param level is the level where the extraction must end
+ * @param requiredLevel is the required level until the extraction must be done
+ * @pre level >= 2
+ */
+ void createQuaternionsUpToLevel(int level, int requiredLevel);
+
+ /**
+ * Method that returns the closest pair node, based on their distance, having the first node of the pair as non-marked.
+ * The pair is found on the given a multimap of distances and pair of nodes. It is supposed that the multimap sorts
+ * the elements by its key. In this case the key corresponds to the distance.
+ * @param map_dist_pairNodes is the map containing the distances (keys) and the pair of nodes (value)
+ * @return the closest pair nodes, based on their distance, given a multimap of distances and pair of nodes. If
+ * no pair was found with the pair node nonmaked, a null pair is returned.
+ * Actually the method returns the pair and its distances in a pair <distance, pair>
+ */
+ std::pair<double, std::pair<Node*, Node*> > GetPairWithClosestDistance_FirstNodeNonMarked(std::multimap<double, pair_nodes> map_dist_pairNodes);
+
+ /**
+ * Add the statistics of number of bifurcations in the tree to the array given by parameter
+ * @param p array to save the statistics
+ */
+ void getStasticsBifurcations(int* p);
+
+ //Operator Overloads
+
+ Node& operator=(Node& node);
+
+ const Node& operator=(const Node& node);
+
+ bool operator ==(const Node& node);
+
+ bool operator !=(const Node& node);
+
+ bool operator <(const Node& node);
+
+ bool operator >(const Node& node);
+
+protected:
+
+ /**
+ * Node id
+ */
+ unsigned int m_id;
+
+ /**
+ * Children nodes vector
+ */
+ vec_nodes m_children;
+
+ /**
+ * Father node
+ */
+ Node* father;
+
+ /**
+ * Edge between the current node and its parent.
+ * The source of the edge corresponds to the parent en the target to the child,
+ * in this case to "this" node.
+ */
+ Edge* m_edge;
+
+ /**
+ * Spatial 3D position of the node in real coordinates of the image
+ */
+ Vec3 m_coords;
+
+ /**
+ * Node level
+ */
+ int m_level;
+
+ /**
+ * Node mark
+ */
+ bool m_mark;
+};
+
+} /* namespace airways */
+
+#endif /* AIRWAYSNODE_H_ */
--- /dev/null
+/*
+ * airwaysTree.txx
+ *
+ * Created on: May 3, 2014
+ * Author: caceres@creatis.insa-lyon.fr
+ * Modifications by: Alfredo Morales Pinzón
+ */
+
+#ifndef _AIRWAYS_TREE_TXX_
+#define _AIRWAYS_TREE_TXX_
+
+#include "airwaysTree.h"
+
+namespace airways
+{
+
+AirwaysTree::AirwaysTree() : m_root(NULL)
+{
+
+}
+
+//Ceron
+AirwaysTree::AirwaysTree(TInputImage::Pointer img, TInputImage::Pointer img_skele, Node* root, bool shouldUpdate){
+ this->m_root = root;
+ this->m_root->SetLevel(0);
+ this->m_img = img;
+ this->m_skeleton = img_skele;
+ //TODO: Should this calculate levels and populate edges and nodes vectors?
+ FillTree(root);
+ if(shouldUpdate) this->UpdateEdges();
+}
+
+void AirwaysTree::FillTree(Node* node){
+ //std::cout << "Calls fill tree for " << node->GetCoords() << std::endl;
+ this->m_nodes.push_back(node);
+ node->SetId(this->m_nodes.size());
+ if(node->GetEdge() != NULL){
+ this->m_edges.push_back(node->GetEdge());
+ node->GetEdge()->UpdateEdgeInfo();
+ }
+ int level = node->GetLevel();
+ //std::cout << "Node is at level " << level << std::endl;
+ vec_nodes::const_iterator it = node->GetChildren().begin();
+ for (; it != node->GetChildren().end(); ++it){
+ (*it)->SetLevel(level+1);
+ this->FillTree(*it);
+ }
+}
+
+AirwaysTree::AirwaysTree(TInputImage::Pointer img, TInputImage::Pointer img_sk, Node* root)
+{
+ this->m_root = root;
+ this->m_root->SetLevel(0);
+ this->m_skeleton = img_sk;
+ this->m_img = img;
+ this->UpdateEdges();
+}
+
+AirwaysTree::AirwaysTree(TInputImage::Pointer image_airwaySegmentation, TInputImage::Pointer image_skeletonDM, EdgeVector& vector_edges, Edge* edge_trachea)
+{
+ std::cout << "AirwaysTree::AirwaysTree(TInputImage::Pointer img, TInputImage::Pointer img_sk, EdgeVector& tEdgeVector, Edge* trachea)" << std::endl;
+ std::cout << "Edges to add:" << vector_edges.size() << std::endl;
+
+ if(edge_trachea == NULL)
+ std::cout << "The trachea edge is NULL" << std::endl;
+
+ // Set the attribtes
+ this->m_skeleton = image_skeletonDM;
+ this->m_img = image_airwaySegmentation;
+ this->m_root = edge_trachea->m_source;
+ this->m_root->SetLevel(0);
+
+ this->m_nodes.push_back(edge_trachea->m_source); // Add the first node, the root to the nodes of the tree
+ edge_trachea->m_source->SetId(this->m_nodes.size());// Set the root id
+
+ // Insert the trachea edge
+ this->Insert(edge_trachea);
+
+ // Insert all edges
+ unsigned int iniSize = vector_edges.size();
+
+ while (!vector_edges.empty())
+ {
+ bool oneEdgeAdded = false;
+ EdgeVector::iterator it = vector_edges.begin();
+ while (it != vector_edges.end())
+ {
+ if (this->Insert(*it))
+ {
+ std::cout << "Edges added: [" << (*it)->GetSource()->GetCoords() << "]->[" << (*it)->GetTarget()->GetCoords() << "], missing: " << vector_edges.size()-1 << std::endl;
+ oneEdgeAdded = true;
+ it = vector_edges.erase(it);
+ }
+ else
+ ++it;
+ } //while
+
+ // If there are edges to be added but they could not be added
+ if(!oneEdgeAdded)
+ {
+ std::cout << "Edges not added:" << vector_edges.size() << std::endl;
+ for (EdgeVector::iterator it = vector_edges.begin(); it != vector_edges.end(); ++it)
+ std::cout << "Edge missing: [" << (*it)->GetSource()->GetCoords() << "]->[" << (*it)->GetTarget()->GetCoords() << "]" << std::endl;
+ break;
+ }
+ // AMP: I did not find any explanation for this "exception"
+ /*count++;
+ if (count == iniSize)
+ {
+ std::cout << "handled exception: arwaysTree.txx, AirwaysTree" << std::endl;
+ break;
+ } //if
+ */
+ } //while
+ this->UpdateEdges();
+
+ std::cout << "AirwaysTree::AirwaysTree(TInputImage::Pointer img, TInputImage::Pointer img_sk, EdgeVector& tEdgeVector, Edge* trachea) ... OK" << std::endl;
+}
+
+AirwaysTree::~AirwaysTree()
+{
+
+}
+
+Node* AirwaysTree::GetRoot() const
+{
+ return this->m_root;
+}
+
+std::string AirwaysTree::GetImageName() const
+{
+ return (this->m_img_name);
+}
+
+std::string AirwaysTree::GetResultPath() const
+{
+ return (this->m_result_path);
+}
+
+TInputImage::Pointer AirwaysTree::GetSegmentedImage() const
+{
+ return this->m_img;
+}
+
+TInputImage::Pointer AirwaysTree::GetSkeletonImage() const
+{
+ return this->m_skeleton;
+}
+
+const Node* AirwaysTree::GetNode(const Vec3 nodeCoords) const
+{
+ return this->GetNode(nodeCoords);
+}
+
+Node* AirwaysTree::GetNodeById(unsigned int nId)
+{
+ return this->m_root->GetNodeById(nId);
+}
+
+int AirwaysTree::GetDepthById(unsigned int nId)
+{
+ return this->m_root->GetDepthById(nId);
+}
+
+unsigned int AirwaysTree::CountMarkedNodes(Node* current) const
+{
+ if (current == NULL)
+ current = this->m_root;
+ Node* aux = NULL;
+ bool found = current->FindMarked(aux);
+ if (!found && !current->IsMarked())
+ return 0;
+ if (aux != NULL && found)
+ current = aux;
+ unsigned int count = 1;
+ vec_nodes::const_iterator it = current->GetChildren().begin();
+ for (; it != current->GetChildren().end(); ++it)
+ count += this->CountMarkedNodes(*it);
+ return count;
+}
+
+unsigned int AirwaysTree::GetWeight( ) const
+{
+ if (this->m_root == NULL)
+ return 0;
+
+ return this->m_root->GetWeight();
+}
+
+unsigned int AirwaysTree::GetHeight(Node* current) const
+{
+ if (current == NULL)
+ current = this->m_root;
+ if (current->IsLeaf())
+ return 1;
+ unsigned int max = 0;
+ vec_nodes::const_iterator it = current->GetChildren().begin();
+ for (; it != current->GetChildren().end(); ++it)
+ {
+ unsigned int cMax = this->GetHeight(*it);
+ if (cMax > max)
+ max = cMax;
+ } //rof
+ if (current == this->m_root)
+ return max;
+ return ++max;
+}
+
+unsigned int AirwaysTree::GetLevels(Node* current, unsigned int cLevel) const
+{
+ if (current == NULL)
+ current = this->m_root;
+ if (current->IsLeaf())
+ return current->GetLevel();
+ unsigned int maxLevel = cLevel;
+ vec_nodes::const_iterator it = current->GetChildren().begin();
+ for (; it != current->GetChildren().end(); ++it)
+ {
+ unsigned int level = this->GetLevels(*it, maxLevel);
+ if (level > maxLevel)
+ maxLevel = level;
+ } //rof
+ return maxLevel;
+}
+
+unsigned int AirwaysTree::GetNumberOfLeafs( ) const
+{
+ if(!this->m_root)
+ return 0;
+
+ return this->m_root->GetNumberOfLeafs( );
+}
+
+const EdgeVector& AirwaysTree::GetEdges() const
+{
+ return this->m_edges;
+}
+
+const vec_nodes& AirwaysTree::GetNodes() const
+{
+ return this->m_nodes;
+}
+
+const Node* AirwaysTree::GetClosestBrachNodeToIndex(float point_x, float point_y, float point_z) const
+{
+ float distance = 0;
+ Node* closestNode = this->m_root->GetClosestBranchNodeToPoint(point_x, point_y, point_z, distance);
+ return closestNode;
+}
+
+void AirwaysTree::GetBrancheNodesInfluencingPoint(float point_x, float point_y, float point_z, vec_nodes& vec_nodes_influence)
+{
+ if(this->m_root)
+ {
+ std::cout << "AirwaysTree::GetBrancheNodesInfluencingPoint" << std::endl;
+ this->m_root->GetBrancheNodesInfluencingPoint(point_x, point_y, point_z, vec_nodes_influence);
+ }
+}
+
+AirwaysTree& AirwaysTree::GetTreeFromMarkedNodes(Node* currentNode, AirwaysTree* tree)
+{
+ if (tree == NULL)
+ {
+ if (!this->m_root->IsMarked())
+ return *(new AirwaysTree());
+ currentNode = this->m_root;
+ tree = new AirwaysTree();
+ Node* current = new Node(currentNode->GetCoords());
+ tree->SetRoot(current);
+ vec_nodes::const_iterator it = currentNode->GetChildren().begin();
+ for (; it != currentNode->GetChildren().end(); ++it)
+ this->GetTreeFromMarkedNodes(*it, tree);
+ this->ImageReconstruction(tree->m_img);
+ //this->ImageReconstruction(tree->m_skeleton);
+ return *tree;
+ } //fi
+ if (currentNode->IsMarked())
+ {
+ Edge* edge = new Edge();
+ *edge = *currentNode->GetEdge();
+ Node* target = new Node(currentNode->GetCoords());
+ edge->m_target = target;
+ tree->Insert(edge);
+ vec_nodes::const_iterator it = currentNode->GetChildren().begin();
+ for (; it != currentNode->GetChildren().end(); ++it)
+ this->GetTreeFromMarkedNodes(*it, tree);
+ } //fi
+ return *(new AirwaysTree());
+}
+
+AirwaysTree& AirwaysTree::GetSubTreeFromNode(Node* node, AirwaysTree* tree)
+{
+ if (tree == NULL)
+ {
+ tree = new AirwaysTree();
+ Node* current = new Node(node->GetCoords());
+ tree->SetRoot(current);
+ vec_nodes::const_iterator it = node->GetChildren().begin();
+ for (; it != node->GetChildren().end(); ++it)
+ this->GetSubTreeFromNode(*it, tree);
+ this->ImageReconstructionFromNode(tree->m_img, node);
+ //this->ImageReconstructionFromNode(tree->m_skeleton, node, true);
+ return *tree;
+ } //fi
+ Edge* edge = new Edge();
+ *edge = *node->GetEdge();
+ Node* target = new Node(node->GetCoords());
+ edge->m_target = target;
+ tree->Insert(edge);
+ vec_nodes::const_iterator it = node->GetChildren().begin();
+ for (; it != node->GetChildren().end(); ++it)
+ this->GetSubTreeFromNode(*it, tree);
+
+ return *(new AirwaysTree());
+}
+
+AirwaysTree& AirwaysTree::GetSubTreeByLevels(const int& level, Node* node,
+ AirwaysTree* tree)
+{
+ if (tree == NULL)
+ {
+ node = this->m_root;
+ tree = new AirwaysTree();
+ Node* current = new Node(node->GetCoords());
+ tree->SetRoot(current);
+ if (node->GetLevel() <= level)
+ {
+ node->Mark();
+ vec_nodes::const_iterator it = node->GetChildren().begin();
+ for (; it != node->GetChildren().end(); ++it)
+ this->GetSubTreeByLevels(level, *it, tree);
+ } //fi
+ if (node->GetLevel() == level)
+ {
+ DuplicatorType::Pointer duplicator = DuplicatorType::New();
+ duplicator->SetInputImage(this->m_img);
+ duplicator->Update();
+ tree->m_img = duplicator->GetOutput();
+ duplicator->SetInputImage(this->m_skeleton);
+ duplicator->Update();
+ tree->m_skeleton = duplicator->GetOutput();
+ } //if
+ else
+ {
+ this->ImageReconstruction(tree->m_img);
+ //this->ImageReconstruction(tree->m_skeleton, true);
+ } //else
+ this->UnMarkAll();
+ return *tree;
+ } //fi
+ if (node->GetLevel() <= level)
+ {
+ node->Mark();
+ Edge* edge = new Edge();
+ *edge = *node->GetEdge();
+ Node* target = new Node(node->GetCoords());
+ edge->m_target = target;
+ tree->Insert(edge);
+ vec_nodes::const_iterator it = node->GetChildren().begin();
+ for (; it != node->GetChildren().end(); ++it)
+ this->GetSubTreeByLevels(level, *it, tree);
+ if (node->GetChildren().size() == 1)
+ {
+ Node* child = *(node->GetChildren().begin());
+ node->MergeNodes(child);
+ const Edge* edge = (child)->GetEdge();
+ //removing from nodeVector and edgeVector
+ EdgeVector::iterator it2 = this->m_edges.begin();
+ for (; it2 != this->m_edges.end(); ++it2)
+ if (edge == *it2)
+ {
+ this->m_edges.erase(it2);
+ break;
+ } //if
+ vec_nodes::iterator it3 = this->m_nodes.begin();
+ for (; it3 != this->m_nodes.end(); ++it3)
+ if (child == *it3)
+ {
+ this->m_nodes.erase(it3);
+ break;
+ } //if
+ } //if
+ } //if
+ return *(new AirwaysTree());
+}
+
+AirwaysTree& AirwaysTree::GetSubTreeByLength(const double& length, Node* node,
+ AirwaysTree* tree, double cLenght)
+{
+ if (tree == NULL)
+ {
+ node = this->m_root;
+ tree = new AirwaysTree();
+ Node* current = new Node(node->GetCoords());
+ tree->SetRoot(current);
+ node->Mark();
+ vec_nodes::const_iterator it = node->GetChildren().begin();
+ for (; it != node->GetChildren().end(); ++it)
+ this->GetSubTreeByLength(length, *it, tree);
+ this->ImageReconstruction(tree->m_img);
+ //this->ImageReconstruction(tree->m_skeleton, true);
+ this->UnMarkAll();
+ return *tree;
+ } //fi
+ cLenght += node->GetEdge()->m_length;
+ if (cLenght <= length)
+ {
+ node->Mark();
+ Edge* edge = new Edge();
+ *edge = *node->GetEdge();
+ Node* target = new Node(node->GetCoords());
+ edge->m_target = target;
+ tree->Insert(edge);
+ vec_nodes::const_iterator it = node->GetChildren().begin();
+ for (; it != node->GetChildren().end(); ++it)
+ this->GetSubTreeByLength(length, *it, tree, cLenght);
+ if (node->GetChildren().size() == 1)
+ {
+ Node* child = *(node->GetChildren().begin());
+ node->MergeNodes(child);
+ const Edge* edge = (child)->GetEdge();
+ //removing from nodeVector and edgeVector
+ EdgeVector::iterator it2 = this->m_edges.begin();
+ for (; it2 != this->m_edges.end(); ++it2)
+ if (edge == *it2)
+ {
+ this->m_edges.erase(it2);
+ break;
+ } //if
+ vec_nodes::iterator it3 = this->m_nodes.begin();
+ for (; it3 != this->m_nodes.end(); ++it3)
+ if (child == *it3)
+ {
+ this->m_nodes.erase(it3);
+ break;
+ } //if
+ } //if
+ } //if
+ return *(new AirwaysTree());
+}
+
+AirwaysTree& AirwaysTree::GetSubTreeByRadius(const double& radius, Node* node,
+ AirwaysTree* tree, double cRadius)
+{
+ if (tree == NULL)
+ {
+ node = this->m_root;
+ tree = new AirwaysTree();
+ Node* current = new Node(node->GetCoords());
+ tree->SetRoot(current);
+ node->Mark();
+ vec_nodes::const_iterator it = node->GetChildren().begin();
+ for (; it != node->GetChildren().end(); ++it)
+ this->GetSubTreeByRadius(radius, *it, tree);
+ this->ImageReconstruction(tree->m_img);
+ //this->ImageReconstruction(tree->m_skeleton, true);
+ this->UnMarkAll();
+ return *tree;
+ } //fi
+ cRadius += node->GetEdge()->m_aRadius;
+ if (cRadius >= radius)
+ {
+ node->Mark();
+ Edge* edge = new Edge();
+ *edge = *node->GetEdge();
+ Node* target = new Node(node->GetCoords());
+ edge->m_target = target;
+ tree->Insert(edge);
+ vec_nodes::const_iterator it = node->GetChildren().begin();
+ for (; it != node->GetChildren().end(); ++it)
+ this->GetSubTreeByRadius(radius, *it, tree, cRadius);
+ if (node->GetChildren().size() == 1)
+ {
+ Node* child = *(node->GetChildren().begin());
+ node->MergeNodes(child);
+ const Edge* edge = (child)->GetEdge();
+ //removing from nodeVector and edgeVector
+ EdgeVector::iterator it2 = this->m_edges.begin();
+ for (; it2 != this->m_edges.end(); ++it2)
+ if (edge == *it2)
+ {
+ this->m_edges.erase(it2);
+ break;
+ } //if
+ vec_nodes::iterator it3 = this->m_nodes.begin();
+ for (; it3 != this->m_nodes.end(); ++it3)
+ if (child == *it3)
+ {
+ this->m_nodes.erase(it3);
+ break;
+ } //if
+ } //if
+ } //if
+ return *(new AirwaysTree());
+}
+
+void AirwaysTree::MarkLevels(const int& level, Node* currentNode)
+{
+ if (currentNode == NULL)
+ currentNode = this->m_root;
+ if (currentNode->GetLevel() <= level)
+ {
+ currentNode->Mark();
+ vec_nodes::const_iterator it = currentNode->GetChildren().begin();
+ for (; it != currentNode->GetChildren().end(); it++)
+ this->MarkLevels(level, *it);
+ } //if
+}
+
+bool AirwaysTree::IsEmpty() const
+{
+ return this->m_root == NULL;
+}
+
+bool AirwaysTree::IsNodeInsidePath(unsigned int idNode_starting, unsigned int idNode_ending, unsigned int idNode_searched)
+{
+ // Check if the starting point exist
+ Node* node_starting = this->m_root->GetNodeById(idNode_starting);
+ if(!node_starting)
+ return false;
+
+ // Check if the searched node exist in the subtree of the starting node
+ Node* node_searched = node_starting->GetNodeById(idNode_searched);
+ if(!node_searched)
+ return false;
+
+ // Check if the ending node exist in the subtree of the searched node
+ if(!node_searched->HasNodeId(idNode_ending))
+ return false;
+
+ // If all conditions were true then the node is in the path
+ return true;
+}
+
+void AirwaysTree::SetRoot(Node* root)
+{
+ this->m_root = root;
+ this->m_root->SetLevel(0);
+ this->m_nodes.push_back(this->m_root);
+}
+
+void AirwaysTree::SetImageName(const std::string& img_name)
+{
+ this->m_img_name = img_name;
+}
+
+void AirwaysTree::SetResultPath(const std::string& folderPath)
+{
+ this->m_result_path = folderPath;
+}
+
+bool AirwaysTree::Insert(Edge* edge)
+{
+ if (this->m_root == NULL)
+ return false;
+ // Search if the source or target coordinates exist in the tree
+ Node* source = this->GetNode(edge->m_source->GetCoords());
+ Node* target = this->GetNode(edge->m_target->GetCoords());
+
+ // If the source and the target coordinates do not exist then the
+ // edge can not be added, also if both exist in the tree.
+ if ( ( source == NULL && target == NULL ) || ( source != NULL && target != NULL) )
+ return false;
+ else if (source == NULL && target != NULL) // If target was found
+ {
+ // Change the direction of the edge
+ edge->m_target = edge->m_source;
+ edge->m_source = target;
+ }
+ else
+ edge->m_source = source;
+
+ // Set the target and source nodes for the Edge
+ edge->m_target->SetEdge(edge);
+ edge->m_target->SetLevel(edge->m_source->GetLevel() + 1);
+ edge->m_source->AddChild(edge->m_target);
+
+
+ // Add the new child node
+ this->m_nodes.push_back(edge->m_target);
+ this->m_edges.push_back(edge);
+
+ // Set the id node
+ edge->m_target->SetId(this->m_nodes.size());
+ return true;
+}
+
+void AirwaysTree::UnMarkAll(Node* currentNode)
+{
+ if (currentNode == NULL)
+ currentNode = this->m_root;
+ currentNode->UnMark();
+
+ vec_nodes::const_iterator it = currentNode->GetChildren().begin();
+ for (; it != currentNode->GetChildren().end(); ++it)
+ this->UnMarkAll(*it);
+}
+
+void AirwaysTree::UpdateLevels(unsigned int levels)
+{
+ this->m_root->UpdateLevels(levels);
+}
+
+void AirwaysTree::CompareTreesOrkiszMorales(AirwaysTree& treeB, unsigned int Q, unsigned int F, vec_nodes& commonA, vec_nodes& commonB, vec_nodes& nonCommonA, vec_nodes& nonCommonB, std::map< unsigned int, std::vector<Node*> >& map_A_to_B, std::map< unsigned int, std::vector<Node*> >& map_B_to_A, std::vector< std::pair< std::pair<Node*, Node*>, std::pair<Node*, Node*> > >& vector_pair_edges_A_to_B)
+{
+ // Algorithm
+ // The algorithm is recursive it is implemented the node class. The idea is to overlap the initial nodes of the subtrees
+ // and find the best branch correspondence for the first level of the trees.
+ // Hypothesis: It is supposed that the roots of the trees are common, so the algorithm does not run in the root nodes
+ // 1. Translate one of the subtree so that the initial nodes overlap.
+ // 2. While one of the trees has non-marked child
+ // 2.1 Compare each child to all the branches is the other tree.
+ // 2.2 Select the pair, child and branch in the other tree, that has the lowest distance function.
+ // 2.3 Add the pair and the distance to the final result.
+ // 2.4 Mark the nodes
+ // 2.5 Find and mark the omitted nodes in the similar branch. If the selected branch has more that 2 nodes it means that intemediary nodes do not have correpondece.
+ // 2.6 Run the process (First step) using the found pair of nodes.
+
+ if( !this->IsEmpty() && !treeB.IsEmpty() )
+ {
+
+ // ----------------------------
+ // ---- Pre-Matching steps ----
+
+ // Get the root nodes
+ Node* root_A = this->m_root;
+ Node* root_B = treeB.m_root;
+
+ // Mark the roots
+ root_A->Mark();
+ root_B->Mark();
+
+ // Add the root match
+ map_A_to_B[root_A->GetId()].push_back(root_B);
+ map_B_to_A[root_B->GetId()].push_back(root_A);
+
+ // Add the root nodes match
+ pair_nodes pair_edgeA(root_A, root_A);
+ pair_nodes pair_edgeB(root_B, root_B);
+ std::pair< pair_nodes, pair_nodes > pairRoorNodes_edge(pair_edgeA, pair_edgeB);
+ vector_pair_edges_A_to_B.push_back(pairRoorNodes_edge);
+
+ // ----------------------------
+ // ----- Run the comparison----
+
+ root_A->CompareSubTreesOrkiszMorales(root_B, Q, F, commonA, commonB, nonCommonA, nonCommonB, map_A_to_B, map_B_to_A, vector_pair_edges_A_to_B);
+
+ }
+ else
+ std::cout << "Trees are empties." << std::endl;
+}
+
+void AirwaysTree::MarkPathFromNodeToNode(Node* node_begin, Node* node_end)
+{
+ if(this->m_root)
+ this->m_root->MarkPathFromNodeToNode(node_begin, node_end);
+}
+
+void AirwaysTree::CompareTrees(AirwaysTree& treeB)
+{
+ // Compare all the "branches" (paths) in treeA with all branches in treeB
+ Vector_Pair_PairNodes_Score vectorPairNodesScore = this->CompareAllBranches(treeB);
+
+ std::cout << "Vector pairs size:" << vectorPairNodesScore.size() << std::endl;
+
+ // Algorithm to find the comment branches (paths)
+ // Until one of the trees is empty
+ // 1. Find the pair, with at least one leaf in the pair, with the lowest score
+ // 2. Compare the score of the father node of the leaf, or one of the leafs, with the correspondent pair leaf
+ // If the score is lower that the score of the actual pair
+ // - The compared leaf is not common. Mark the compared leaf, delete the compared leave and save the leaf as not common.
+ // If the score is greater
+ // - The leafs are common then mark both leafs as commons. Mark both leaves, delete them and save them as common leafs.
+
+ // Variables
+ //AirwaysTree* tree_copyA = this->GetCopy();
+ //AirwaysTree* tree_copyB = treeB.GetCopy();
+
+ //NodeVector nonCommonNodes_A;
+ //NodeVector commonNodes_A;
+ //NodeVector nonCommonNodes_B;
+ //NodeVector commonNodes_B;
+
+ std::map<unsigned int, Node*> nonCommonNodes_A;
+ std::map<unsigned int, Node*> commonNodes_A;
+ std::map<unsigned int, Node*> nonCommonNodes_B;
+ std::map<unsigned int, Node*> commonNodes_B;
+
+ // Until one of the trees is not empty
+ while( this->GetWeight() > 1 && treeB.GetWeight() > 1 )
+ {
+ // Print the weight of the trees
+ std::cout << "Weight [A,B]: [" << this->GetWeight() << "," << treeB.GetWeight() << "]" << std::endl;
+ std::cout << "Leafs [A,B]: [" << this->GetNumberOfLeafs() << "," << treeB.GetNumberOfLeafs() << "]" << std::endl;
+ std::cout << "Vector pairs size:" << vectorPairNodesScore.size() << std::endl;
+
+
+ // 1. Find the pair, with at least one leaf in the pair, with the lowest score
+ Pair_PairNodes_Score pair_nonMarked = this->GetBestLeaf(vectorPairNodesScore);
+ Node* node_actualA = pair_nonMarked.first.first;
+ Node* node_actualB = pair_nonMarked.first.second;
+ double score_actual = pair_nonMarked.second;
+
+ // 2. Take the leaf (L) and his father (F), the other node is called (O). Compare the score of the father node (F) with the correspondent pair leaf (O).
+ if(node_actualA->IsLeaf())
+ {
+ //const Edge* edgeActual = node_actualA->GetEdge();
+ //Node* targetActual = edgeActual->GetTarget();
+ //const unsigned int id_father = targetActual->GetId();
+ //double score_fatherA_to_B = this->GetPairNodes_Score(vectorPairNodesScore, 0, node_actualB->GetId());
+ double score_fatherA_to_B = this->GetPairNodes_Score(vectorPairNodesScore, node_actualA->GetEdge()->GetTarget()->GetId(), node_actualB->GetId());
+ //double score_fatherA_to_B = pair_fatherA_to_B.second;
+ // If the score is lower than the score of the actual pair
+ if(score_fatherA_to_B < score_actual)
+ {
+ //- The compared leaf is not common. Mark the compared leaf, delete the compared leave and save the leaf as not common.
+ node_actualA->Mark();
+ std::map<unsigned int, Node*>::iterator it;
+ it=nonCommonNodes_A.find(node_actualA->GetId());
+ if(it==nonCommonNodes_A.end())
+ nonCommonNodes_A[node_actualA->GetId()] = node_actualA;
+ //nonCommonNodes_A.push_back(node_actualA);
+ }
+ // If the score is greater
+ else
+ {
+ // - The leafs are common then mark both leafs as commons. Mark both leaves, delete them and save them as common leafs.
+ node_actualA->Mark();
+ node_actualB->Mark();
+
+ std::map<unsigned int, Node*>::iterator it;
+ it=commonNodes_B.find(node_actualB->GetId());
+ if(it==commonNodes_B.end())
+ commonNodes_B[node_actualB->GetId()] = node_actualB;
+
+ it=commonNodes_A.find(node_actualA->GetId());
+ if(it==commonNodes_A.end())
+ {
+ std::cout << "CommonNode A:" << node_actualA->GetId() << std::endl;
+ commonNodes_A[node_actualA->GetId()] = node_actualA;
+ }
+
+ //commonNodes_B.push_back(node_actualB);
+ //commonNodes_A.push_back(node_actualA);
+ }
+ this->DeleteNodeById(node_actualA->GetId());
+ this->DeleteEntriesByIdNode(vectorPairNodesScore, 0, node_actualA->GetId());
+ }
+ else
+ {
+ double score_fatherB_to_A = this->GetPairNodes_Score(vectorPairNodesScore, node_actualA->GetId(), node_actualB->GetEdge()->GetTarget()->GetId());
+ //double score_fatherB_to_A = pair_fatherB_to_A.second;
+ // If the score is lower that the score of the actual pair
+ if(score_fatherB_to_A < score_actual)
+ {
+ //- The compared leaf is not common. Mark the compared leaf, delete the compared leave and save the leaf as not common.
+ node_actualB->Mark();
+ std::map<unsigned int, Node*>::iterator it;
+ it=nonCommonNodes_B.find(node_actualB->GetId());
+ if(it==nonCommonNodes_B.end())
+ nonCommonNodes_B[node_actualB->GetId()] = node_actualB;
+ //nonCommonNodes_B.push_back(node_actualB);
+ }
+ // If the score is greater
+ else
+ {
+ // - The leafs are common then mark both leafs as commons. Mark both leaves, delete them and save them as common leafs.
+ node_actualA->Mark();
+ node_actualB->Mark();
+ //commonNodes_B.push_back(node_actualB);
+ //commonNodes_A.push_back(node_actualA);
+
+ std::map<unsigned int, Node*>::iterator it;
+ it=commonNodes_B.find(node_actualB->GetId());
+ if(it==commonNodes_B.end())
+ commonNodes_B[node_actualB->GetId()] = node_actualB;
+
+ it=commonNodes_A.find(node_actualA->GetId());
+ if(it==commonNodes_A.end())
+ commonNodes_A[node_actualA->GetId()] = node_actualA;
+ }
+ treeB.DeleteNodeById(node_actualB->GetId());
+ this->DeleteEntriesByIdNode(vectorPairNodesScore, 1, node_actualB->GetId());
+ }
+ }
+
+ // Print the results
+ std::cout << "Non-CommonNodes [A,B]" << nonCommonNodes_A.size() << ", " << nonCommonNodes_B.size() << std::endl;
+ std::cout << "CommonNodes [A,B]" << commonNodes_A.size() << ", " << commonNodes_B.size() << std::endl;
+
+ // Get the branches with the lowest difference
+ /*bool first = true;
+ double lowestValue = 0;
+ Pair_Node_Node bestPairNodes;
+ int iteration = 1;
+ for(Vector_Pair_PairNodes_Score::iterator it = vectorPairNodesScore.begin(); it != vectorPairNodesScore.end(); ++it )
+ {
+ if(first)
+ {
+ lowestValue = (*it).second;
+ bestPairNodes = (*it).first;
+ first = false;
+ }
+ else if( (*it).second < lowestValue )
+ {
+ lowestValue = (*it).second;
+ bestPairNodes = (*it).first;
+ }
+ }
+
+ // Print the best pair of nodes
+ if(vectorPairNodesScore.size() > 0)
+ std::cout << "Best pair: Node A:" << (bestPairNodes).first->GetId() << ", Node B:" << (bestPairNodes).second->GetId() << std::endl;
+ */
+
+ /*// AMP - This code was commented after a meeting with Leonardo Florez
+ // We make the hypothesis that both roots are equal and
+ // that the root has only one child.
+ // The comparison of two branches implies the alignment
+ // of the mother branches and then the comparison of the
+ // daughters branches using intrinsic characteristics of
+ // each branch and a distance function between both.
+
+ if (this->IsEmpty() || treeB.IsEmpty())
+ {
+ std::cout << "One or both airways are empty." << std::endl;
+ return;
+ }
+
+ Edge actualRootEdge = this->m_root->GetEdge();
+ Edge comparedRootEdge = treeB.m_root->GetEdge();
+
+ // By the hypothesis, both root edges are similar. So we have to compared
+ // the edged of the root edges.
+ actualRootEdge.compareChildEdges(comparedRootEdge);
+ */
+}
+
+Pair_PairNodes_Score AirwaysTree::GetBestLeaf(Vector_Pair_PairNodes_Score vectorPairNodesScore)
+{
+ // Variables
+ Pair_PairNodes_Score pairNodesScore_final;
+ bool first = true;
+
+ for(Vector_Pair_PairNodes_Score::iterator it = vectorPairNodesScore.begin(); it != vectorPairNodesScore.end(); ++it)
+ {
+ // Check that one of the node is a leaf
+ pair_nodes pair_actual = (*it).first;
+ Node* a = pair_actual.first;
+ Node* b = pair_actual.second;
+ if( a->IsLeaf() || b->IsLeaf() )
+ {
+ if(first)
+ {
+ first = false;
+ pairNodesScore_final = (*it);
+ }
+ else if((*it).second < pairNodesScore_final.second)
+ {
+ pairNodesScore_final = (*it);
+ }
+ }
+ }
+
+ if(first)
+ std::cout << "Best leaf not found" << std::endl;
+
+ return pairNodesScore_final;
+}
+
+double AirwaysTree::GetPairNodes_Score(Vector_Pair_PairNodes_Score vectorPairNodesScore, const unsigned int idA, const unsigned int idB)
+{
+ bool found = false;
+ double score = -1;
+
+ for(Vector_Pair_PairNodes_Score::iterator it = vectorPairNodesScore.begin(); it != vectorPairNodesScore.end(); ++it)
+ {
+ // Check that one of the node is a leaf
+ pair_nodes pair_actual = (*it).first;
+ Node* a = pair_actual.first;
+ Node* b = pair_actual.second;
+ if( a->GetId() == idA && b->GetId() == idB )
+ {
+ found = true;
+ score = (*it).second;
+ }
+ }
+
+ if(!found)
+ std::cout << "Pair of Ids not found: " << idA << ", " << idB << std::endl;
+
+ return score;
+}
+
+bool AirwaysTree::DeleteNodeById(unsigned int nId)
+{
+ bool deleted = false;
+ if(m_root->GetId() == nId)
+ {
+ delete(m_root);
+ m_root = NULL;
+ return true;
+ }
+
+ return this->m_root->DeleteNodeById(nId);
+}
+
+unsigned int AirwaysTree::DeleteEntriesByIdNode(Vector_Pair_PairNodes_Score& vectorPairNodesScore, unsigned int component, unsigned int nId)
+{
+ unsigned int pairs_deleted = 0;
+
+ Vector_Pair_PairNodes_Score::iterator it = vectorPairNodesScore.begin();
+ while(it != vectorPairNodesScore.end())
+ {
+ // Check that one of the node is a leaf
+ pair_nodes pair_actual = (*it).first;
+ unsigned int id_actual = -1;
+ if(component == 0)
+ id_actual = pair_actual.first->GetId();
+ else
+ id_actual = pair_actual.second->GetId();
+
+ if( id_actual == nId)
+ {
+ it = vectorPairNodesScore.erase(it);
+ ++pairs_deleted;
+ }
+ else
+ ++it;
+ }
+
+ std::cout << "Pairs deleted: " << pairs_deleted << std::endl;
+
+ return pairs_deleted;
+}
+
+AirwaysTree* AirwaysTree::GetCopy()
+{
+ return NULL;
+}
+
+Vector_Pair_PairNodes_Score AirwaysTree::CompareAllBranches(AirwaysTree& treeB)
+{
+ std::cout << "Compare all branches ..." << std::endl;
+ // Variables
+ int nodesA = this->m_nodes.size();
+ int nodesB = treeB.m_nodes.size();
+
+ Vector_Pair_PairNodes_Score vectorPairNodesScore;
+
+ // For all the actual branches ( nodes different from the root)
+ for(unsigned int idA = 2; idA <= nodesA; ++idA)
+ {
+ // Get the actual branch in A
+ Edge* branchA = this->GetComposedEdge(idA);
+
+ //Get the final actual node in A
+ Node* node_actualA = this->GetNodeById(idA);
+
+ std::cout << "IdA: " << idA << std::endl;
+
+ // For all the branches in B tree
+ for(unsigned int idB = 2; idB <= nodesB; ++idB)
+ {
+ // Get the actual branch in B
+ Edge* branchB = treeB.GetComposedEdge(idB);
+
+ //Get the final actual node in B
+ Node* node_actualB = treeB.GetNodeById(idB);
+
+ // Get the score from A to B and from B to A
+ double scoreA2B = branchA->CompareWith(branchB);
+ double scoreB2A = branchB->CompareWith(branchA);
+
+ // Get the final score
+ // Options:
+ // A. Maximum score
+ //double score = scoreB2A > scoreA2B ? scoreB2A : scoreA2B ;
+
+ // B. Sum of scores
+ double score = scoreB2A + scoreA2B;
+
+ // Build the pair of nodes
+ pair_nodes actualPairNodes(node_actualA, node_actualB);
+
+ Pair_PairNodes_Score actualPairNodesScore(actualPairNodes, score);
+ vectorPairNodesScore.push_back(actualPairNodesScore);
+
+ //std::cout << "IdA: " << idA << ", IdB: " << idB << ", score: " << score << ", ... OK" << std::endl;
+ delete(branchB);
+ }
+ delete(branchA);
+ }
+
+ std::cout << "Compare all branches ... OK" << std::endl;
+ return vectorPairNodesScore;
+}
+
+Edge* AirwaysTree::GetComposedEdge(unsigned int idNode)
+{
+ return this->m_root->GetComposedEdge(idNode);
+}
+
+AirwaysTree* AirwaysTree::GetSingleDetachedTreeNodeById(unsigned int nId)
+{
+ AirwaysTree* tree_detached = NULL;
+ Node* node_searched = this->GetNodeById(nId);
+ if(node_searched)
+ {
+ Node* node_root_detached = node_searched->GetDetachedRootCopy();
+ tree_detached = new AirwaysTree();
+ tree_detached->m_root = node_root_detached;
+ }
+ else
+ std::cout << "Node not found" << std::endl;
+ return tree_detached;
+}
+
+void AirwaysTree::SubIsomorphism(AirwaysTree& treeB)
+{
+ if (this->IsEmpty() || treeB.IsEmpty())
+ return;
+ //Marking the root node
+ this->m_root->Mark();
+ treeB.m_root->Mark();
+
+ //iterating children - first step
+ vec_nodes::const_iterator it = this->m_root->GetChildren().begin();
+ for (; it != this->m_root->GetChildren().end(); ++it)
+ {
+ unsigned int totalCost = 0;
+ Node* currentB = NULL;
+ Node* currentA = *it;
+ if (currentA->IsMarked())
+ continue;
+ vec_nodes::const_iterator it2 = treeB.m_root->GetChildren().begin();
+ for (; it2 != treeB.m_root->GetChildren().end(); ++it2)
+ {
+ if ((*it2)->IsMarked())
+ continue;
+ unsigned int currentCost = (*it)->GetCost(*it2);
+ if (totalCost < currentCost)
+ {
+ currentB = *it2;
+ totalCost = currentCost;
+ } //if
+ } //for
+ if (currentB != NULL)
+ if (*currentA == *currentB)
+ this->GetSubTreeByDetachingNode(currentA).SubIsomorphism(
+ treeB.GetSubTreeByDetachingNode(currentB));
+ } //for
+ //end of first step
+ //second step
+ it = this->m_root->GetChildren().begin();
+ for (; it != this->m_root->GetChildren().end(); ++it)
+ {
+ if ((*it)->IsMarked())
+ continue;
+ vec_nodes::const_iterator it2 = treeB.m_root->GetChildren().begin();
+ for (; it2 != treeB.m_root->GetChildren().end(); ++it2)
+ {
+ Node* currentB = *it2;
+ Node* currentA = *it;
+ if (*currentB > *currentA)
+ continue;
+ if (currentA->FindSimilarEdgeByAccumulation(currentB))
+ {
+ AirwaysTree subTreeA = this->GetSubTreeByDetachingNode(
+ currentA);
+ AirwaysTree subTreeB = treeB.GetSubTreeByDetachingNode(
+ currentB);
+ subTreeB.UpdateLevels((*it2)->GetLevel());
+ subTreeA.SubIsomorphism(subTreeB);
+ treeB.m_root->UpdateLevels(treeB.m_root->GetLevel());
+ currentA->Mark();
+ if (!(*it2)->MarkPath(currentB))
+ std::cout << "Error finding path" << std::endl;
+ break;
+ } //if
+ } //for
+ } //for
+ vec_nodes::const_iterator it2 = treeB.m_root->GetChildren().begin();
+ for (; it2 != treeB.m_root->GetChildren().end(); ++it2)
+ {
+ if ((*it2)->IsMarked())
+ continue;
+ it = this->m_root->GetChildren().begin();
+ for (; it != this->m_root->GetChildren().end(); ++it)
+ {
+ Node* currentB = *it2;
+ Node* currentA = *it;
+ if (*currentA > *currentB)
+ continue;
+ if (currentB->FindSimilarEdgeByAccumulation(currentA))
+ {
+ AirwaysTree subTreeA = this->GetSubTreeByDetachingNode(
+ currentA);
+ AirwaysTree subTreeB = treeB.GetSubTreeByDetachingNode(
+ currentB);
+ subTreeA.UpdateLevels((*it)->GetLevel());
+ subTreeA.SubIsomorphism(subTreeB);
+ this->m_root->UpdateLevels(this->m_root->GetLevel());
+ currentB->Mark();
+ if (!(*it)->MarkPath(currentA))
+ std::cout << "Error finding path" << std::endl;
+ break;
+ } //if
+ } //for
+ }
+ //End of second step
+}
+
+bool AirwaysTree::Isomorphism(AirwaysTree& tree)
+{
+ if ((this->GetWeight() != tree.GetWeight())
+ || (this->GetHeight() != tree.GetHeight()))
+ return false;
+ return this->m_root->CompareNodes(tree.m_root);
+}
+
+void AirwaysTree::ImageReconstructionFromNode(TInputImage::Pointer& result, Node* node, bool skeleton)
+{
+ if (node == NULL)
+ return;
+ //Finding the center of the sphere
+ const Edge* edge = node->GetEdge();
+ if (edge == NULL)
+ return;
+ //making a copy of the current image
+ TInputImage::Pointer copy;
+ DuplicatorType::Pointer duplicator = DuplicatorType::New();
+ if (skeleton)
+ duplicator->SetInputImage(this->m_skeleton);
+ else
+ duplicator->SetInputImage(this->m_img);
+ duplicator->Update();
+ copy = duplicator->GetOutput();
+ //First remove the node brothers -- the following code in a way does not
+ // make sense but I'll just comment it if it is necessary
+ /*Node* parent = edge->m_source;
+ NodeVector::iterator pIt = parent->m_children.begin();
+ for (; pIt != parent->m_children.end(); ++pIt)
+ if ((*pIt) != node)
+ this->RemoveBranchFromImage(copy, node);
+ //end of brother removal*/
+ pair_posVox_rad center;
+ vec_pair_posVox_rad::const_iterator it = edge->m_vec_pair_posVox_rad.begin();
+ for (; it != edge->m_vec_pair_posVox_rad.end(); ++it)
+ if (edge->m_source->GetCoords() == (*it).first)
+ center = *it;
+ //converting to index
+ TInputImage::PointType point;
+ point[0] = center.first[0];
+ point[1] = center.first[1];
+ point[2] = center.first[2];
+ TInputImage::IndexType indCenter;
+ if (skeleton)
+ this->m_skeleton->TransformPhysicalPointToIndex(point, indCenter);
+ else
+ this->m_img->TransformPhysicalPointToIndex(point, indCenter);
+ //creating itk sphere
+ SphereType::Pointer sphere = SphereType::New();
+ sphere->SetRadius(center.second * 4);
+ SphereType::InputType sphereCenter;
+ sphereCenter[0] = indCenter[0];
+ sphereCenter[1] = indCenter[1];
+ sphereCenter[2] = indCenter[2];
+ sphere->SetCenter(sphereCenter);
+
+ //testing code
+ // std::cout<<"center = "<< indCenter << std::endl;
+ //std::cout << "Radio = " << center.second << std::endl;
+ TInputImage::SizeType regionSize;
+ regionSize[0] = 8 * center.second;
+ regionSize[1] = 8 * center.second;
+ regionSize[2] = 8 * center.second;
+ //std::cout << "RegionSize = " << regionSize << std::endl;
+
+ TInputImage::IndexType regionIndex;
+ regionIndex[0] = indCenter[0] - center.second * 4;
+ regionIndex[1] = indCenter[1] - center.second * 4;
+ regionIndex[2] = indCenter[2] - center.second * 4;
+ // std::cout << "regionIndex = " << regionIndex << std::endl;
+
+ TInputImage::RegionType region;
+ region.SetSize(regionSize);
+ region.SetIndex(regionIndex);
+ RegionIterator it2(copy, region);
+ while (!it2.IsAtEnd())
+ {
+ TInputImage::IndexType ind = it2.GetIndex();
+ SphereType::InputType point;
+ point[0] = ind[0];
+ point[1] = ind[1];
+ point[2] = ind[2];
+ SphereType::OutputType output = sphere->Evaluate(point);
+ if (output)
+ it2.Set(0);
+ ++it2;
+ }
+
+ /*
+ * First region growing - Segmenting the wanted branch to be removed
+ * Removing the branch from the image
+ */
+ if (!ComputeSimpleRegionGrowing(copy, this->m_root->GetCoords()))
+ return;
+ SubtractImageFilterType::Pointer subtractFilter =
+ SubtractImageFilterType::New();
+ if (skeleton)
+ subtractFilter->SetInput1(this->m_skeleton);
+ else
+ subtractFilter->SetInput1(this->m_img);
+ subtractFilter->SetInput2(copy);
+ subtractFilter->Update();
+ TInputImage::Pointer diff = subtractFilter->GetOutput();
+ /*
+ * Second region growing - Cleaning the resulting image
+ */
+ Node* leaf = this->GetLeaf(node);
+ if (ComputeSimpleRegionGrowing(diff, leaf->GetCoords()))
+ result = diff;
+ else
+ std::cout << "problemas" << std::endl;
+ //cin.ignore().get(); //Pause Command for Linux Terminal
+
+}
+
+void AirwaysTree::ImageReconstruction(TInputImage::Pointer& result, bool skeleton, Node* node)
+{
+ if (node == NULL)
+ {
+ DuplicatorType::Pointer duplicator = DuplicatorType::New();
+ if (skeleton)
+ duplicator->SetInputImage(this->m_skeleton);
+ else
+ duplicator->SetInputImage(this->m_img);
+ duplicator->Update();
+ result = duplicator->GetOutput();
+ node = this->m_root;
+ } //if
+ if (!node->IsMarked())
+ {
+ this->RemoveBranchFromImage(result, node);
+ return;
+ }
+ vec_nodes::const_iterator it = node->GetChildren().begin();
+ for (; it != node->GetChildren().end(); ++it)
+ this->ImageReconstruction(result, skeleton, *it);
+
+}
+
+void AirwaysTree::ImageReconstructionBySpheres(TInputImage::Pointer& result, Node* node, bool useMarks)
+{
+ if (node == NULL)
+ {
+ node = this->m_root;
+ DuplicatorType::Pointer duplicator = DuplicatorType::New();
+ duplicator->SetInputImage(m_img);
+ duplicator->Update();
+ result = duplicator->GetOutput();
+ //Removing white pixels -- cleaning image
+ TInputImage::RegionType lRegion = result->GetLargestPossibleRegion();
+ RegionIterator it(result, lRegion);
+ while (!it.IsAtEnd())
+ {
+ it.Set(0);
+ ++it;
+ } //while
+ } //if
+ if (node->IsMarked() || !useMarks)
+ {
+ this->CreateSpheres(result, node);
+ vec_nodes::const_iterator it = node->GetChildren().begin();
+ for (; it != node->GetChildren().end(); ++it)
+ this->ImageReconstructionBySpheres(result, *it);
+ } //if
+}
+
+Node* AirwaysTree::GetNode(const Vec3 nodeCoords, Node* current)
+{
+ if (current == NULL)
+ current = this->m_root;
+ if (current->GetCoords() == nodeCoords)
+ return current;
+ vec_nodes::const_iterator it = current->GetChildren().begin();
+ for (; it != current->GetChildren().end(); ++it)
+ {
+ Node* child = GetNode(nodeCoords, *it);
+ if (child != NULL)
+ return child;
+ } //rof
+ return NULL;
+}
+
+Node* AirwaysTree::GetLeaf(Node* current)
+{
+ if (current == NULL)
+ current = this->m_root;
+ if (current->IsLeaf())
+ return current;
+ vec_nodes::const_iterator it = current->GetChildren().begin();
+ Node* node = NULL;
+ unsigned int maxHeight = 0;
+ for (; it != current->GetChildren().end(); ++it)
+ {
+ unsigned int cHeight = this->GetHeight(*it);
+ if (maxHeight <= cHeight)
+ {
+ maxHeight = cHeight;
+ node = *it;
+ } //if
+ } //for
+ if (node != NULL)
+ return this->GetLeaf(node);
+ return NULL;
+}
+
+void AirwaysTree::UpdateEdges(Node* node)
+{
+ if (node == NULL)
+ node = this->m_root;
+
+ if (node->GetEdge() != NULL)
+ {
+ Edge* edge = node->GetEdge();
+ //edge->m_skInfoPairVector = this->GetSkeletonInfoFromEdge(edge->m_source->GetCoords(), edge->m_target->GetCoords());
+ edge->SetSkeletonPairVector( this->GetSkeletonInfoFromEdge(edge->m_source->GetCoords(), edge->m_target->GetCoords()) );
+ node->GetEdge()->UpdateEdgeInfo();
+ }
+ vec_nodes::const_iterator it = node->GetChildren().begin();
+ for (; it != node->GetChildren().end(); ++it)
+ this->UpdateEdges(*it);
+}
+
+AirwaysTree& AirwaysTree::GetSubTreeByDetachingNode(Node* node)
+{
+ AirwaysTree* tree = new AirwaysTree();
+ tree->m_root = node;
+ return *tree;
+}
+
+DiGraph_EdgePair AirwaysTree::AddBoostEdge(const pair_posVox_rad& source, const pair_posVox_rad& target, DiGraph& graph)
+{
+ bool vFound[2] = { false, false };
+ DiGraph_VertexDescriptor vDescriptor[2];
+
+ //Remember that DiGraph_VertexIteratorPair, first = begin, second = end
+ DiGraph_VertexIndexMap index = get(boost::vertex_index, graph);
+
+ //Search the vertex in the graph
+ for (DiGraph_VertexIteratorPair vPair = boost::vertices(graph); vPair.first != vPair.second && (!vFound[0] || !vFound[1]); ++vPair.first)
+ {
+ // Get the current information [vector, intensity]
+ pair_posVox_rad current = graph[index[*vPair.first]];
+
+ // Check if current is equal to the source or target vertex
+ if (current.first == source.first)
+ {
+ vDescriptor[0] = index[*vPair.first];
+ vFound[0] = true;
+ } //if
+ else if (current.first == target.first)
+ {
+ vDescriptor[1] = index[*vPair.first];
+ vFound[1] = true;
+ } //if
+ } //for
+
+ // If one of the vertices is not in the graph, adds the new vertex
+
+ // Add the vertices if they are not in the graph
+ if (!vFound[0])
+ vDescriptor[0] = boost::add_vertex(source, graph);
+ if (!vFound[1])
+ vDescriptor[1] = boost::add_vertex(target, graph);
+
+ // Check if the edge exist in the graph
+ // boost::edge returns the pair<edge_descriptor, bool>, with bool=true is the edge exists
+ DiGraph_EdgePair edgeC1 = boost::edge(vDescriptor[0], vDescriptor[1],graph);
+ DiGraph_EdgePair edgeC2 = boost::edge(vDescriptor[1], vDescriptor[0],graph);
+
+ // If the edge does not exist, then add it
+ if (!edgeC1.second && !edgeC2.second)
+ {
+ DiGraph_EdgePair nEdge = boost::add_edge(vDescriptor[0], vDescriptor[1], graph);
+ return nEdge;
+ } //if
+
+ // If an edge has not been added then -> second = false
+ // Remove the added vertices because the edge already exists.
+ if (!vFound[0])
+ boost::remove_vertex(vDescriptor[0], graph);
+
+ if (!vFound[1])
+ boost::remove_vertex(vDescriptor[1], graph);
+
+ edgeC1.second = false;
+ edgeC1.second = false;
+
+ return edgeC1;
+}
+
+bool AirwaysTree::FindNeighborhood(DiGraph& graph, Vec3Queue& queue, Vec3List& used, const Vec3& end, ConstNeighborhoodIterator& iterator)
+{
+ if (queue.empty())
+ return false;
+
+ // Get the next pixel in the queue
+ Vec3 current = queue.front();
+
+ // Put the pixel in the used list
+ used.push_back(current);
+
+ // Delete next pixel from the queu
+ queue.pop();
+
+ //AMP//iterator.GoToBegin();
+ //AMP//while (!iterator.IsAtEnd())
+ //AMP//{
+
+ // AMP commented
+ /*unsigned int centerIndex = iterator.GetCenterNeighborhoodIndex();
+ TInputImage::IndexType index = iterator.GetIndex(centerIndex);
+ SKPairInfo src(Vec3(index[0], index[1], index[2]),
+ iterator.GetPixel(iterator.GetCenterNeighborhoodIndex()));
+ if (src.first != current)
+ {
+ ++iterator;
+ continue;
+ } //if
+ */
+ //PMA commented
+
+ // AMP
+
+ // Set the current index to the actual pixel
+ TInputImage::IndexType indexActual;
+ indexActual[0] = current[0];
+ indexActual[1] = current[1];
+ indexActual[2] = current[2];
+
+ // Place the iterator in the actual index
+ iterator.SetLocation(indexActual);
+
+ // Create the pair [vector, intensity] for the actual index
+ pair_posVox_rad src(current, iterator.GetPixel(iterator.GetCenterNeighborhoodIndex()));
+
+ // PMA
+
+ // Iterate over the neighbors
+ for (unsigned int i = 0; i < iterator.Size(); i++)
+ {
+ //AMP//if (iterator.GetPixel(i) != 0 && i != centerIndex)
+
+ // If the intensity is not zero and it's the center
+ if (iterator.GetPixel(i) != 0 && i != iterator.GetCenterNeighborhoodIndex())
+ {
+ // Get the index of the iterator
+ TInputImage::IndexType ind = iterator.GetIndex(i);
+
+ // Create the target pair [vector, intensity] for the target
+ pair_posVox_rad target(Vec3(ind[0], ind[1], ind[2]), iterator.GetPixel(i));
+
+ // Add the edge to the graph
+ bool added = AddBoostEdge(src, target, graph).second;
+
+ // If it could not be added then it is because it already exist, then go to next neighboor
+ if (!added)
+ continue;
+
+ // The edge was added. Stop the adding process if the target vertex was reached.
+ if (target.first == end)
+ return true;
+
+ // If it the edge was added and it is not the target then add it to the queue, if and only if
+ // it was not already used.
+ if (std::find(used.begin(), used.end(), target.first)== used.end())
+ queue.push(target.first);
+ }//if
+ }//for
+ //AMP//break;
+ //AMP//} //while
+ return false;
+}
+
+DiGraph AirwaysTree::GetGraphShortestPath(const DiGraph& graph, const Vec3& begin, const Vec3& end)
+{
+ // Variables
+ DiGraph ret; // Returning graph with the shortest path
+ DiGraph_EdgeIterator ei, eo; // Initial and final iterators
+
+ // Iterate over the edges of the graph
+ boost::tie(ei, eo) = boost::edges(graph);
+ for (; ei != eo; ++ei)
+ {
+ // Get the source of the actual edge
+ pair_posVox_rad src = graph[boost::source(*ei, graph)];
+
+ // If the source corresponds to the begin pixel
+ if (src.first == begin)
+ {
+ // Get the target of the actual edge
+ pair_posVox_rad target = graph[boost::target(*ei, graph)];
+
+ // If the actual target corresponds to the end pixel
+ if (target.first == end)
+ {
+ // Add the edge to the final path
+ AddBoostEdge(src, target, ret);
+
+ //AMP
+ std::cout << "Add final target [src]-[ target]:[" << src.first << "]-[" << target.first << "], numVertices:" << boost::num_vertices(ret) << std::endl;
+ //PMA
+
+ // Return the path
+ return ret;
+ }
+
+ // If the actual target does not correspond to the end pixel, then add it to the returning path
+ // and run the shortest path with the target of the actual edge.
+ ret = GetGraphShortestPath(graph, target.first, end);
+
+ // If the path is empty then the iteration continues
+ if (boost::num_vertices(ret) == 0)
+ continue;
+
+ // if not return ret + new edge
+ // If the path was found then add the actual edge and return the path
+ AddBoostEdge(src, target, ret);
+ //AMP
+ //std::cout << "Add found edge [src]-[ target]:[" << src.first << "]-[" << target.first << "], numVertices:" << boost::num_vertices(ret) << std::endl;
+ //PMA
+ return ret;
+
+ }//fi
+ }
+ return ret;
+}
+
+vec_pair_posVox_rad AirwaysTree::GetGraphShortestPath_AMP(const DiGraph& graph, const Vec3& begin, const Vec3& end)
+{
+ //std::cout << "SKPairInfoVector AirwaysTree::GetGraphShortestPath_AMP(const DiGraph& graph, const Vec3& begin, const Vec3& end)" << std::endl;
+ // Variables
+ vec_pair_posVox_rad path;
+ DiGraph ret; // Returning graph with the shortest path
+ DiGraph_EdgeIterator ei, eo; // Initial and final iterators
+
+ // Iterate over the edges of the graph
+ boost::tie(ei, eo) = boost::edges(graph);
+ for (; ei != eo; ++ei)
+ {
+ // Get the source of the actual edge
+ pair_posVox_rad src = graph[boost::source(*ei, graph)];
+
+ // If the source corresponds to the begin pixel
+ if (src.first == begin)
+ {
+ // Get the target of the actual edge
+ pair_posVox_rad target = graph[boost::target(*ei, graph)];
+
+ // If the actual target corresponds to the end pixel
+ if (target.first == end)
+ {
+ // Add the edge to the final path
+ AddBoostEdge(src, target, ret);
+ path.push_back(target);
+ path.push_back(src);
+
+ //AMP
+ //std::cout << "Add final target [src]-[ target]:[" << src.first << "]-[" << target.first << "], numVertices:" << boost::num_vertices(ret) << ", path_size:" << path.size() << std::endl;
+ //PMA
+
+ // Return the path
+ //std::cout << "SKPairInfoVector AirwaysTree::GetGraphShortestPath_AMP(const DiGraph& graph, const Vec3& begin, const Vec3& end) ... FIRST ... OK" << std::endl;
+ return path;
+ }
+
+ //std::cout << "Not yet in the target voxel" << std::endl;
+ // If the actual target does not correspond to the end pixel, then add it to the returning path
+ // and run the shortest path with the target of the actual edge.
+ path = GetGraphShortestPath_AMP(graph, target.first, end);
+ //std::cout << "Path found, path size:" << path.size() << std::endl;
+
+ // If the path is empty then the iteration continues
+ if (path.size() == 0)
+ continue;
+
+ // if not return ret + new edge
+ // If the path was found then add the actual edge and return the path
+ AddBoostEdge(src, target, ret);
+ path.push_back(src);
+ //AMP
+ //std::cout << "Add found edge [src]-[ target]:[" << src.first << "]-[" << target.first << "], numVertices:" << boost::num_vertices(ret) << ", path_size:" << path.size() << std::endl;
+ //PMA
+ //std::cout << "SKPairInfoVector AirwaysTree::GetGraphShortestPath_AMP(const DiGraph& graph, const Vec3& begin, const Vec3& end) ... SECOND ... OK" << std::endl;
+ return path;
+
+ }//fi
+ }
+
+ //std::cout << "SKPairInfoVector AirwaysTree::GetGraphShortestPath_AMP(const DiGraph& graph, const Vec3& begin, const Vec3& end) ... END ... OK" << std::endl;
+ return path;
+}
+
+vec_pair_posVox_rad AirwaysTree::GetSkeletonInfoFromEdge(const Vec3& source, const Vec3& target)
+{
+ //converting to pixel
+ TInputImage::PointType pointSource;
+ pointSource[0] = source[0];
+ pointSource[1] = source[1];
+ pointSource[2] = source[2];
+
+ TInputImage::PointType pointTarget;
+ pointTarget[0] = target[0];
+ pointTarget[1] = target[1];
+ pointTarget[2] = target[2];
+
+ TInputImage::IndexType indexSource;
+ TInputImage::IndexType indexTarget;
+
+ this->m_skeleton->TransformPhysicalPointToIndex(pointSource, indexSource);
+ this->m_skeleton->TransformPhysicalPointToIndex(pointTarget, indexTarget);
+ //end of conversion
+
+ // AMP // Print the initial and final voxels to verify their positions
+ //std::cout << "Begin: [" << indexSource[0] << "," << indexSource[1] << "," << indexSource[2] << "] -- " << this->m_skeleton->GetPixel(indexSource) << std::endl;
+ //std::cout << "End: [" << indexTarget[0] << "," << indexTarget[1] << "," << indexTarget[2] << "] -- " << this->m_skeleton->GetPixel(indexTarget) << std::endl;
+ //std::cout << "BeginVoxel: [" << pointSource[0] << "," << pointSource[1] << "," << pointSource[2] << "]"<< std::endl;
+ //std::cout << "EndVoxel: [" << pointTarget[0] << "," << pointTarget[1] << "," << pointTarget[2] << "]" << std::endl;
+ // PMA //
+
+ Vec3 sourcePxl(indexSource[0], indexSource[1], indexSource[2]);
+ Vec3 targetPxl(indexTarget[0], indexTarget[1], indexTarget[2]);
+
+ // AMP // Commented
+ /*
+ // Compute the distance in each direction and add "21" in each one
+ float dstX = beginPxl.EculideanDistance(
+ Vec3(endPxl[0], beginPxl[1], beginPxl[2])) + 21;
+ float dstY = beginPxl.EculideanDistance(
+ Vec3(beginPxl[0], endPxl[1], beginPxl[2])) + 21;
+ float dstZ = beginPxl.EculideanDistance(
+ Vec3(beginPxl[0], beginPxl[1], endPxl[2])) + 21;
+
+ //start point
+ TInputImage::IndexType index;
+ index[0] = beginPxl[0] < endPxl[0] ? beginPxl[0] : endPxl[0];
+ index[1] = beginPxl[1] < endPxl[1] ? beginPxl[1] : endPxl[1];
+ index[2] = beginPxl[2] < endPxl[2] ? beginPxl[2] : endPxl[2];
+
+ index[0] -= index[0] <= 10 ? 0 : 10;
+ index[1] -= index[1] <= 10 ? 0 : 10;
+ index[2] -= index[2] <= 10 ? 0 : 10;
+
+
+ //size of region
+ TInputImage::SizeType size;
+ size[0] = ceil(dstX);
+ size[1] = ceil(dstY);
+ size[2] = ceil(dstZ);
+
+ TInputImage::RegionType region;
+ region.SetSize(size);
+ region.SetIndex(index);
+
+ //ConstNeighborhoodIterator iterator(radius, this->m_skeleton, region);
+ // PMA // Commented
+ */
+
+ TInputImage::SizeType radius;
+ radius[0] = 1;
+ radius[1] = 1;
+ radius[2] = 1;
+ ConstNeighborhoodIterator iterator(radius, this->m_skeleton, this->m_skeleton->GetLargestPossibleRegion());
+
+ DiGraph graph;
+ Vec3Queue queue;
+ Vec3List used;
+
+ // Set the first pixel as starting pixel
+ queue.push(sourcePxl);
+
+ //FindNeighborhood(graph, queue, used, endPxl, iterator);
+
+ // Add to the graph all the connected voxel until the ending pixel is found
+ bool endFound = false;
+ while (!queue.empty() && !endFound)
+ endFound = FindNeighborhood(graph, queue, used, targetPxl, iterator);
+
+ if(!endFound)
+ std::cout << "***************** END NOT FOUND!!!!" << std::endl;
+ else
+ std::cout << "***************** END FOUND!!!! - Graph Vertices:" << boost::num_vertices(graph) << ", Graph Edges:" << boost::num_edges(graph) << std::endl;
+
+ // Get the shortest path between the begin and end pixels
+ //std::cout << "GetGraphShortestPath [from]-[to]:[" << sourcePxl << "]-[" << targetPxl << "]" << std::endl;
+ //DiGraph result = GetGraphShortestPath(graph, sourcePxl, targetPxl);
+ vec_pair_posVox_rad path = GetGraphShortestPath_AMP(graph, sourcePxl, targetPxl);
+ //std::cout << "GetGraphShortestPath ... OK , numVertices:" << boost::num_vertices(result) << std::endl;
+
+ //AMP
+ bool beginPixelFound = false;
+ bool endPixelFound = false;
+ vec_pair_posVox_rad vector;
+
+ for(vec_pair_posVox_rad::reverse_iterator it_path = path.rbegin(); it_path != path.rend(); ++it_path)
+ {
+ TInputImage::IndexType ind;
+ ind[0] = (*it_path).first[0];
+ ind[1] = (*it_path).first[1];
+ ind[2] = (*it_path).first[2];
+ //std::cout << "["<< (*it_path).first << "], ";
+
+ // To avoid precision problems
+ if (indexSource == ind)
+ {
+ beginPixelFound = true;
+ (*it_path).first = source;
+ }
+ else if (indexTarget == ind)
+ {
+ endPixelFound = true;
+ (*it_path).first = target;
+ }
+ else
+ {
+ TInputImage::PointType pnt;
+ this->m_skeleton->TransformIndexToPhysicalPoint(ind, pnt);
+ (*it_path).first = Vec3(pnt[0], pnt[1], pnt[2]);
+ }
+ vector.push_back((*it_path));
+ }
+ //vector = path;
+ //PMA
+
+ // Variables to iterate over the shortest path and save the results
+ /*// AMP
+ DiGraph_VertexIterator i, e;
+ SKPairInfoVector vector;
+
+ // Iterate over the shortest path and add each vertex in real coordinates
+ boost::tie(i, e) = boost::vertices(result);
+ std::cout << "Inserting ... [begin,end]: [" << begin << "], [" << end << "]" << std::endl;
+ for (; i != e; ++i)
+ //for (; e != i; --e)
+ {
+ // Get the actual vertex
+ SKPairInfo vertex = result[*i];
+ //SKPairInfo vertex = result[*e];
+ std::cout << "["<< vertex.first << "], ";
+
+ // Get the index (voxel position) of the actual vertex
+ TInputImage::IndexType ind;
+ ind[0] = vertex.first[0];
+ ind[1] = vertex.first[1];
+ ind[2] = vertex.first[2];
+
+ // To avoid precision problems
+ if (indexBegin == ind)
+ {
+ beginPixelFound = true;
+ vertex.first = begin;
+ }
+ else if (indexEnd == ind)
+ {
+ endPixelFound = true;
+ vertex.first = end;
+ }
+ else
+ {
+ TInputImage::PointType pnt;
+ this->m_skeleton->TransformIndexToPhysicalPoint(ind, pnt);
+ vertex.first = Vec3(pnt[0], pnt[1], pnt[2]);
+ }
+ vector.push_back(vertex);
+ }
+ std::cout << "Inserting ... OK" << std::endl;
+
+
+ *///AMP
+
+
+ if(!beginPixelFound)
+ std::cout << "Begin not found: "<< source << std::endl;
+ if(!endPixelFound)
+ std::cout << "End not found:" << target << std::endl;
+
+ return vector;
+}
+
+bool AirwaysTree::ComputeSimpleRegionGrowing(TInputImage::Pointer& img,
+ const Vec3 seedPtn)
+{
+ TInputImage::PointType targetPtn;
+ targetPtn[0] = seedPtn[0];
+ targetPtn[1] = seedPtn[1];
+ targetPtn[2] = seedPtn[2];
+ TInputImage::IndexType seed;
+ img->TransformPhysicalPointToIndex(targetPtn, seed);
+ if (img->GetPixel(seed) == 0)
+ {
+ //std::cout << "seed is 0" << std::endl;
+ return false;
+ }
+ ConnectedFilterType::Pointer neighborhoodConnected =
+ ConnectedFilterType::New();
+ neighborhoodConnected->SetLower(1);
+ neighborhoodConnected->SetUpper(1);
+ neighborhoodConnected->SetReplaceValue(1);
+ neighborhoodConnected->SetSeed(seed);
+ neighborhoodConnected->SetInput(img);
+ neighborhoodConnected->Update();
+ img = neighborhoodConnected->GetOutput();
+ return true;
+}
+
+void AirwaysTree::RemoveBranchFromImage(TInputImage::Pointer& img, Node* node)
+{
+ //typedef itk::ImageFileWriter<TInputImage> WriterType;
+ //Finding the center of the sphere
+ const Edge* edge = node->GetEdge();
+ if (edge == NULL)
+ return;
+ pair_posVox_rad center;
+ vec_pair_posVox_rad::const_iterator it = edge->m_vec_pair_posVox_rad.begin();
+ for (; it != edge->m_vec_pair_posVox_rad.end(); ++it)
+ if (edge->m_source->GetCoords() == (*it).first)
+ center = *it;
+ //converting to index
+ TInputImage::PointType point;
+ point[0] = center.first[0];
+ point[1] = center.first[1];
+ point[2] = center.first[2];
+ TInputImage::IndexType indCenter;
+ img->TransformPhysicalPointToIndex(point, indCenter);
+ double radius = center.second;
+ //creating itk sphere
+ SphereType::Pointer sphere = SphereType::New();
+ SphereType::InputType sphereCenter;
+ sphereCenter[0] = indCenter[0];
+ sphereCenter[1] = indCenter[1];
+ sphereCenter[2] = indCenter[2];
+ sphere->SetCenter(sphereCenter);
+ unsigned int objCount = 0;
+ unsigned int obj3Count = 0; // to assure 3 Connected obj
+ //Number of connected components, children_size + parent
+ unsigned int nCConnex = edge->m_source->GetNumberOfChildren() + 1;
+ TInputImage::Pointer copy;
+ //std::cout << "entró" << std::endl;
+ while (objCount < nCConnex || obj3Count <= nCConnex)
+ {
+ if (objCount >= 3)
+ obj3Count++;
+ radius += 1.0;
+ sphere->SetRadius(radius);
+ //making a copy of the current image
+ DuplicatorType::Pointer duplicator = DuplicatorType::New();
+ duplicator->SetInputImage(img);
+ duplicator->Update();
+ copy = duplicator->GetOutput();
+
+ //testing code
+ // std::cout<<"center = "<< indCenter << std::endl;
+ //std::cout << "Radio = " << center.second << std::endl;
+ TInputImage::SizeType regionSize;
+ regionSize[0] = 4 * radius;
+ regionSize[1] = 4 * radius;
+ regionSize[2] = 4 * radius;
+ //std::cout << "RegionSize = " << regionSize << std::endl;
+ TInputImage::SizeType size = img->GetLargestPossibleRegion().GetSize();
+ TInputImage::IndexType regionIndex;
+ regionIndex[0] = indCenter[0] - (radius * 2);
+ if (regionIndex[0] < 0)
+ regionIndex[0] = indCenter[0];
+ else if (regionIndex[0] > (unsigned int) size[0])
+ regionIndex[0] = size[0];
+ regionIndex[1] = indCenter[1] - (radius * 2);
+ if (regionIndex[1] < 0)
+ regionIndex[1] = indCenter[1];
+ else if (regionIndex[1] > (unsigned int) size[1])
+ regionIndex[1] = size[1];
+ regionIndex[2] = indCenter[2] - (radius * 2);
+ if (regionIndex[2] < 0)
+ regionIndex[2] = indCenter[2];
+ else if (regionIndex[2] > (unsigned int) size[2])
+ regionIndex[2] = size[2];
+ int diff[3];
+ diff[0] = size[0] - (regionSize[0] + regionIndex[0]);
+ diff[1] = size[1] - (regionSize[1] + regionIndex[1]);
+ diff[2] = size[2] - (regionSize[2] + regionIndex[2]);
+ //std::cout << "(1) Region Index = " << regionIndex << std::endl;
+ //std::cout << "(1) Region size = " << regionSize << std::endl;
+ if (diff[0] < 0)
+ regionSize[0] = size[0] - regionIndex[0];
+ if (diff[1] < 0)
+ regionSize[1] = size[1] - regionIndex[1];
+ if (diff[2] < 0)
+ regionSize[2] = size[2] - regionIndex[2];
+ // std::cout << "regionIndex = " << regionIndex << std::endl;
+ //std::cout << "(2) Region Index = " << regionIndex << std::endl;
+ //std::cout << "(2) Region size = " << regionSize << std::endl;
+ TInputImage::RegionType region;
+ region.SetSize(regionSize);
+ region.SetIndex(regionIndex);
+ RegionIterator it2(copy, region);
+ while (!it2.IsAtEnd())
+ {
+ TInputImage::IndexType ind = it2.GetIndex();
+ SphereType::InputType point;
+ point[0] = ind[0];
+ point[1] = ind[1];
+ point[2] = ind[2];
+ SphereType::OutputType output = sphere->Evaluate(point);
+ if (output)
+ it2.Set(0);
+ ++it2;
+ } //while
+ if (node->IsLeaf())
+ {
+ TInputImage::PointType ptn;
+ ptn[0] = node->GetCoords()[0];
+ ptn[1] = node->GetCoords()[1];
+ ptn[2] = node->GetCoords()[2];
+ TInputImage::IndexType tgt;
+ img->TransformPhysicalPointToIndex(ptn, tgt);
+ if (img->GetPixel(tgt) == 0)
+ {
+ std::cout << "It entered here!!!!" << std::endl;
+ return;
+ }
+ } //if
+ CastFilterType::Pointer castFilter = CastFilterType::New();
+ castFilter->SetInput(copy);
+ castFilter->Update();
+ ConnectedComponentImageFilterType::Pointer connFilter =
+ ConnectedComponentImageFilterType::New();
+ connFilter->SetFullyConnected(true);
+ connFilter->SetInput(castFilter->GetOutput());
+ connFilter->Update();
+ objCount = connFilter->GetObjectCount();
+ //test
+ /*WriterType::Pointer writer = WriterType::New();
+ writer->SetInput(copy);
+ writer->SetFileName("output_antes.mhd");
+ writer->Update();*/
+ } //while
+ //std::cout << "salio" << std::endl;
+ /*
+ * First region growing - Segmenting the wanted branch to be removed
+ * Removing the branch from the image
+ */
+ Node* leaf = this->GetLeaf(node);
+ if (!ComputeSimpleRegionGrowing(copy, leaf->GetCoords()))
+ return;
+ SubtractImageFilterType::Pointer subtractFilter =
+ SubtractImageFilterType::New();
+ subtractFilter->SetInput1(img);
+ subtractFilter->SetInput2(copy);
+ subtractFilter->Update();
+ TInputImage::Pointer diff = subtractFilter->GetOutput();
+ /*
+ * Second region growing - Cleaning the resulting image
+ */
+ if (ComputeSimpleRegionGrowing(diff, this->m_root->GetCoords()))
+ img = diff;
+ /*else
+ std::cout << "problemas" << std::endl;*/
+ //cin.ignore().get(); //Pause Command for Linux Terminal
+}
+
+//Metodo obsoleto, ya no nos interesa la reconstruccion por esferas
+void AirwaysTree::CreateSpheres(TInputImage::Pointer& img, Node* node)
+{
+ if (node->GetEdge() == NULL)
+ return;
+ const Edge* edge = node->GetEdge();
+ vec_pair_posVox_rad::const_iterator it = edge->m_vec_pair_posVox_rad.begin();
+ for (; it != edge->m_vec_pair_posVox_rad.end(); ++it)
+ {
+ //converting to index
+ TInputImage::PointType point;
+ point[0] = (*it).first[0];
+ point[1] = (*it).first[1];
+ point[2] = (*it).first[2];
+ TInputImage::IndexType indexPoint;
+ img->TransformPhysicalPointToIndex(point, indexPoint);
+ //creating itk sphere
+ SphereType::Pointer sphere = SphereType::New();
+ sphere->SetRadius((*it).second);
+ SphereType::InputType sphereCenter;
+ sphereCenter[0] = indexPoint[0];
+ sphereCenter[1] = indexPoint[1];
+ sphereCenter[2] = indexPoint[2];
+ sphere->SetCenter(sphereCenter);
+ //creating region
+ TInputImage::SizeType regionSize;
+ regionSize[0] = 4 * (*it).second;
+ regionSize[1] = 4 * (*it).second;
+ regionSize[2] = 4 * (*it).second;
+ TInputImage::IndexType regionIndex;
+ regionIndex[0] = indexPoint[0] - 2 * ceil((*it).second);
+ regionIndex[1] = indexPoint[1] - 2 * ceil((*it).second);
+ regionIndex[2] = indexPoint[2] - 2 * ceil((*it).second);
+ TInputImage::RegionType region;
+ region.SetSize(regionSize);
+ region.SetIndex(regionIndex);
+ RegionIterator it2(img, region);
+ while (!it2.IsAtEnd())
+ {
+ TInputImage::IndexType ind = it2.GetIndex();
+ SphereType::InputType point;
+ point[0] = ind[0];
+ point[1] = ind[1];
+ point[2] = ind[2];
+ SphereType::OutputType output = sphere->Evaluate(point);
+ if (output)
+ it2.Set(1);
+ ++it2;
+ } //while
+ } //for
+}
+
+void AirwaysTree::printUpToLevel(int level)
+{
+ cout << "Printing up to level: " << level <<endl;
+ this->m_root->printUpToLevel(level);
+ cout << "Printing done" <<endl;
+}
+
+void AirwaysTree::printNodeAndChildrenIds()
+{
+ cout << "Printing nodes and children: " <<endl;
+ this->m_root->printNodeAndChildrenIds( );
+ cout << "Printing nodes and children DONE" <<endl;
+}
+
+void AirwaysTree::createQuaternionsUpToLevel(int level)
+{
+ cout << "Creating quaternions up to level: " << level <<endl;
+ if(level >= 2)
+ {
+ // Print header
+ cout << "Level" << " " << "SonPositionX" << " " << "SonPositionY" << " " << "SonPositionZ" << " "
+ "ParentVectorX" << " " << "ParentVectorY" << " " << "ParentVectorZ" << " " <<
+ "ChildVectorX" << " " << "ChildVectorY" << " " << "ChildVectorZ"<< " " <<
+ "Qr" << " " << "Qx" << " " << "Qy" " " << "Qz" " " << "Qtheta" << " " <<
+ "Wx" << " " << "Wy" << " " << "Wz" << " " <<
+ "EuclDistSonParent" << " " <<"EuclDistSonGrandson" << " " <<
+ "RadMeanSonParent" << " " << "RadMeanSonGrandson" <<endl;
+
+ this->m_root->createQuaternionsUpToLevel(level, level);
+ }
+ else
+ cout << "Quaternions must be created from level 2." <<endl;
+ cout << "Creating quaternions done" <<endl;
+}
+
+void AirwaysTree::getStatisticsBifurcations()
+{
+ cout << "Statistics bifurcations ... " <<endl;
+
+ int p[6] = {0,0,0,0,0,0};
+
+ this->m_root->getStasticsBifurcations(p);
+
+ cout << "Bifurcations:" << endl;
+ for(int i = 0; i < 6; i++)
+ {
+ cout << i+1 << "->" << p[i] << endl;
+ }
+
+ cout << "Statistics bifurcations done " <<endl;
+}
+
+void AirwaysTree::saveAirwayToImage(std::string filename, int dims[3], double spc[3], double nOrigin[3])
+{
+ TInputImage::Pointer imagePointer;
+ imagePointer = TInputImage::New();
+
+ TInputImage::IndexType start;
+ start[0] = 0;
+ start[1] = 0;
+ start[2] = 0;
+
+ TInputImage::SizeType size;
+ size[0] = dims[0];
+ size[1] = dims[1];
+ size[2] = dims[2];
+
+ TInputImage::RegionType region;
+ region.SetSize( size );
+ region.SetIndex( start );
+
+ TInputImage::SpacingType spacing;
+ spacing[0] = spc[0];
+ spacing[1] = spc[1];
+ spacing[2] = spc[2];
+
+ TInputImage::PointType origin;
+ origin[0] = nOrigin[0];
+ origin[1] = nOrigin[1];
+ origin[2] = nOrigin[2];
+
+ imagePointer->SetRegions( region );
+ imagePointer->SetSpacing( spacing );
+ imagePointer->SetOrigin( origin );
+ imagePointer->Allocate();
+
+ //this->m_root->saveToImage(_imagePointer);
+ vec_nodes nodes_children_root = this->m_root->GetChildren();
+ vec_nodes::const_iterator it_children = nodes_children_root.begin();
+ for( ; it_children != nodes_children_root.end(); ++it_children)
+ {
+ Vec3 coords_father = (*it_children)->GetEdge()->GetSource()->GetCoords();
+ Vec3 coords_child = (*it_children)->GetCoords();
+ drawLineInImage(coords_father, coords_child, imagePointer);
+ }
+
+ cout << "Writing ..." << endl;
+ typedef itk::ImageFileWriter<TInputImage> WriterType;
+ WriterType::Pointer writer = WriterType::New();
+ writer->SetInput(imagePointer);
+ writer->SetFileName(filename);
+ writer->Update();
+ cout << "Writing ... done" << endl;
+}
+
+void AirwaysTree::drawLineInImage(Vec3 initVoxel, Vec3 endVoxel, TInputImage::Pointer imagePointer)
+{
+ TInputImage::RegionType region = imagePointer->GetLargestPossibleRegion();
+
+ TInputImage::SpacingType spacing = imagePointer->GetSpacing();
+
+ TInputImage::PointType origin = imagePointer->GetOrigin();
+
+ TInputImage::IndexType init;
+ init[0]=(initVoxel[0]-origin[0])/spacing[0];
+ init[1]=(initVoxel[1]-origin[1])/spacing[1];
+ init[2]=(initVoxel[2]-origin[2])/spacing[2];
+
+ TInputImage::IndexType end;
+ end[0]=(endVoxel[0]-origin[0])/spacing[0];
+ end[1]=(endVoxel[1]-origin[1])/spacing[1];
+ end[2]=(endVoxel[2]-origin[2])/spacing[2];
+
+ TPixel pixel;
+ pixel = 1;
+
+ itk::LineIterator<TInputImage> it(imagePointer, init, end);
+ it.GoToBegin();
+ while (!it.IsAtEnd())
+ {
+ it.Set(pixel);
+ ++it;
+ }
+}
+
+//To store the graph in a vtk
+void AirwaysTree::saveAirwayAsVTK(const std::string filename, bool common)
+{
+ vtkSmartPointer<vtkUnsignedCharArray> colors = vtkSmartPointer<vtkUnsignedCharArray>::New();
+ colors->SetNumberOfComponents(3);
+ colors->SetName("Colors");
+
+ //pointer of lines and points
+ vtkSmartPointer<vtkCellArray> lines = vtkSmartPointer<vtkCellArray>::New();
+ vtkSmartPointer<vtkPoints> pts = vtkSmartPointer<vtkPoints>::New();
+
+ Vec3 root = this->m_root->GetCoords();
+
+ //UndirectedGraph
+
+ srand(time(NULL));
+ unsigned int id = 1;
+
+ CalculateVTKLinesFromEdges(this->m_root, 0, id, pts, lines, colors, common);
+
+ vtkSmartPointer<vtkPolyData> linesPolyData = vtkSmartPointer<vtkPolyData>::New();
+ linesPolyData->SetPoints(pts);
+ linesPolyData->SetLines(lines);
+
+ linesPolyData->GetCellData()->SetScalars(colors);
+
+ // Write the file
+ vtkSmartPointer<vtkPolyDataWriter> writer = vtkSmartPointer<vtkPolyDataWriter>::New();
+ writer->SetFileName(filename.c_str());
+#if VTK_MAJOR_VERSION <= 5
+ writer->SetInput(linesPolyData);
+#else
+ writer->SetInputData(linesPolyData);
+#endif
+ writer->Write();
+}
+
+//To store the graph in a vtk
+void AirwaysTree::CalculateVTKLinesFromEdges(const Node* node, const unsigned int& parentId,
+ unsigned int& cId, vtkSmartPointer<vtkPoints>& pts,
+ vtkSmartPointer<vtkCellArray>& lines,
+ vtkSmartPointer<vtkUnsignedCharArray>& colors, bool common)
+{
+ if (node == NULL)
+ return;
+
+ // Colors
+ unsigned char red[3] = { 255, 0, 0 };
+ unsigned char green[3] = { 0, 255, 0 };
+ unsigned char blue[3] = { 0, 0, 255 };
+ unsigned char yellow[3] = { 255, 255, 0 };
+
+ const vec_nodes children = node->GetChildren();
+
+ for(vec_nodes::const_iterator it = children.begin(); it != children.end(); ++it)
+ {
+ if (!(*it)->IsMarked() && common)
+ continue;
+
+ const Edge* edge = (*it)->GetEdge();
+ //pts->InsertNextPoint(edge->GetSource()->GetCoords().GetVec3());
+ //pts->InsertNextPoint(edge->GetTarget()->GetCoords().GetVec3());
+
+ pts->InsertNextPoint(node->GetCoords().GetVec3());
+ pts->InsertNextPoint((*it)->GetCoords().GetVec3());
+
+ // Set color to be used
+ int numColor = node->GetLevel() % 4;
+ if(numColor == 0)
+ colors->InsertNextTupleValue(green);
+ else if(numColor == 1)
+ colors->InsertNextTupleValue(red);
+ else if(numColor == 2)
+ colors->InsertNextTupleValue(blue);
+ else
+ colors->InsertNextTupleValue(yellow);
+
+ vtkSmartPointer<vtkLine> line = vtkSmartPointer<vtkLine>::New();
+ line->GetPointIds()->SetId(0, parentId);
+ line->GetPointIds()->SetId(1, cId++);
+ lines->InsertNextCell(line);
+ unsigned int newParent = cId++;
+ CalculateVTKLinesFromEdges(*it, newParent, cId, pts, lines, colors, common);
+ } //for
+}
+
+
+} /* namespace airways */
+
+
+#endif
+
--- /dev/null
+/*
+ * airwaysTree.h
+ *
+ * Created on: May 3, 2014
+ * Author: caceres@creatis.insa-lyon.fr
+ */
+
+#ifndef _AIRWAYS_TREE_H_
+#define _AIRWAYS_TREE_H_
+
+#include "airwaysTreeTypeDefinition.h"
+
+// VTK
+#include <vtkVersion.h>
+#include <vtkSmartPointer.h>
+#include <vtkCellArray.h>
+#include <vtkCellData.h>
+#include <vtkLine.h>
+#include <vtkPolyDataMapper.h>
+#include <vtkRenderWindow.h>
+#include <vtkRenderer.h>
+#include <vtkRenderWindowInteractor.h>
+#include <vtkPolyDataWriter.h>
+#include <vtkSphereSource.h>
+
+// ITK
+#include "itkImage.h"
+#include "itkImageRegionIterator.h"
+#include "itkLineIterator.h"
+#include <itkImageFileWriter.h>
+#include <appli/TempAirwaysAppli/AirwaysLib/AirwaysLib_Export.h>
+
+namespace airways
+{
+
+// Types
+
+// Vector of nodes
+typedef std::vector<Node*> vec_nodes;
+
+// Pair of nodes
+typedef std::pair<Node*, Node*> pair_nodes;
+
+//Pair of a pair of nodes and a score
+typedef std::pair<pair_nodes,double> Pair_PairNodes_Score;
+
+/**
+ *
+ */
+//typedef std::pair<pair_nodes,double> Pair_PairNodes_Score;
+
+/**
+ * Vector of pairs of a pair of nodes and a score
+ */
+typedef std::vector<Pair_PairNodes_Score> Vector_Pair_PairNodes_Score;
+
+/*!
+ * A vector stocking all the points between two nodes
+ */
+typedef std::vector<pair_posVox_rad> vec_pair_posVox_rad;
+
+class AirwaysLib_EXPORT AirwaysTree
+{
+public:
+
+ AirwaysTree();
+
+ AirwaysTree(TInputImage::Pointer img, TInputImage::Pointer img_sk, Node* root);
+
+ AirwaysTree(TInputImage::Pointer image_airwaySegmentation, TInputImage::Pointer image_skeletonDM, EdgeVector& vector_edges, Edge* edge_trachea);
+
+
+ //CERON
+ // New constructor that doesn't update path info between nodes.
+ AirwaysTree(TInputImage::Pointer image, TInputImage::Pointer image_skeleton, Node* root, bool shouldUpdate);
+ virtual ~AirwaysTree();
+
+ /**
+ * Getters
+ */
+
+ Node* GetRoot() const;
+
+ std::string GetImageName() const;
+
+ std::string GetResultPath() const;
+
+ TInputImage::Pointer GetSegmentedImage() const;
+
+ TInputImage::Pointer GetSkeletonImage() const;
+
+ const Node* GetNode(const Vec3 nodeCoords) const;
+
+ /**
+ * Method that return the node with the specified id
+ * @param nId is id of the requested node
+ * @return the requested id, NULL if not found
+ */
+ Node* GetNodeById(unsigned int nId);
+
+ /**
+ * Method that returns the depth of the node whose id is given by parameter
+ * @pre the node exist in the sub-tree
+ * @param idNode is the search node id
+ * @return the depth of the node up to the actual node
+ */
+ int GetDepthById(unsigned int nId);
+
+ unsigned int CountMarkedNodes(Node* current = NULL) const;
+
+ /**
+ * Method that returns the weight of the tree
+ * @return the weight of the tree which is the number of nodes in the tree.
+ */
+ unsigned int GetWeight( ) const;
+
+ unsigned int GetHeight(Node* current = NULL) const;
+
+ unsigned int GetLevels(Node* current = NULL, unsigned int cLevel = 0) const;
+
+ /**
+ * Method that returns the number of leafs in the tree
+ * @return the number of leafs in the tree
+ */
+ unsigned int GetNumberOfLeafs() const;
+
+ const EdgeVector& GetEdges() const;
+
+ const vec_nodes& GetNodes() const;
+
+ /**
+ * Method that returns the closest node of a given voxel.
+ * The closes node is the one whose path to the father has the closest point to the given point
+ * in real coordinates of the image.
+ * @param voxel_x x componente of the voxel
+ * @param voxel_y y componente of the voxel
+ * @param voxel_z z componente of the voxel
+ * @return the closest node to the given voxel
+ */
+ const Node* GetClosestBrachNodeToIndex(float point_x, float point_y, float point_z) const;
+
+ /**
+ * Method that returns the list of nodes that influence the given point
+ * @param point_x, point_y, point_z are the position of the point in real coordinates x, y, and z respectively.
+ * @param vec_nodes_influence vector to save the nodes
+ */
+ void GetBrancheNodesInfluencingPoint(float point_x, float point_y, float point_z, vec_nodes& vec_nodes_influence);
+
+ AirwaysTree& GetTreeFromMarkedNodes(Node* currentNode = NULL, AirwaysTree* tree = NULL);
+
+ AirwaysTree& GetSubTreeFromNode(Node* node, AirwaysTree* tree = NULL);
+
+ AirwaysTree& GetSubTreeByLevels(const int& level, Node* node = NULL, AirwaysTree* tree = NULL);
+
+ AirwaysTree& GetSubTreeByLength(const double& length, Node* node = NULL, AirwaysTree* tree = NULL, double cLenght = 0.0);
+
+ AirwaysTree& GetSubTreeByRadius(const double& radius, Node* node = NULL, AirwaysTree* tree = NULL, double cRadius = 0.0);
+
+ void MarkLevels(const int& level, Node* currentNode = NULL);
+
+ bool IsEmpty() const;
+
+ /**
+ * Method that search a node inside a path. The path is defined by the starting and ending nodes
+ * @param idNode_starting is the id of the starting node
+ * @param idNode_ending is the id of the ending node
+ * @param idNode_searched is the of the searched node
+ * @return true if the searched node is inside the path, false otherwise
+ */
+ bool IsNodeInsidePath(unsigned int idNode_starting, unsigned int idNode_ending, unsigned int idNode_searched);
+
+ void SetRoot(Node* root);
+
+ void SetImageName(const std::string& img_name);
+
+ void SetResultPath(const std::string& folderPath);
+
+ bool Insert(Edge* edge);
+
+ void UnMarkAll(Node* currentNode = NULL);
+
+ /**
+ * Method that compares two tree starting from the root node.
+ * @param treeB is the tree to be compared with
+ * @param Q is the depth to select "fathers" nodes
+ * @param F is the depth to select "family" nodes
+ * @param commonA vector to save the common nodes for this subtree
+ * @param commonA vector to save the common nodes for the subtree given by parameter
+ * @param nonCommonA vector to save the non-common nodes for this subtree
+ * @param nonCommonB vector to save the non-common nodes for the subtree given by parameter
+ */
+ void CompareTreesOrkiszMorales(AirwaysTree& treeB, unsigned int Q, unsigned int F, vec_nodes& commonA, vec_nodes& commonB, vec_nodes& nonCommonA, vec_nodes& nonCommonB, std::map< unsigned int, std::vector<Node*> >& map_A_to_B, std::map< unsigned int, std::vector<Node*> >& map_B_to_A, std::vector< std::pair< std::pair<Node*, Node*>, std::pair<Node*, Node*> > >& vector_pair_edges_A_to_B);
+
+ void MarkPathFromNodeToNode(Node* node_begin, Node* node_end);
+
+ /**
+ * Method that compares the actual tree with another given by parameter
+ * @param treeB the tree to be compared with
+ */
+ void CompareTrees(AirwaysTree& treeB);
+
+ /**
+ * Method that finds the pair, at least one is a leaf, with the lowest score
+ * @param vectorPairNodesScore is the vector of pairs
+ * @return the pair, at least one is a leaf, with the lowest score
+ */
+ Pair_PairNodes_Score GetBestLeaf(Vector_Pair_PairNodes_Score vectorPairNodesScore);
+
+ /**
+ * Method that returns the pair of pair nodes and score for the given ids.
+ * @param vectorPairNodesScore vector of pair nodes and score
+ * @param idA is the first id if the pair nodes
+ * @param idB is the second id if the pair nodes
+ * @return the score for the given pair.
+ */
+ double GetPairNodes_Score(Vector_Pair_PairNodes_Score vectorPairNodesScore, const unsigned int idA, const unsigned int idB);
+
+ /**
+ * Method that returns a copy of the tree
+ * @return the copy of this tree
+ */
+ AirwaysTree* GetCopy();
+
+ /**
+ * Method that deleted a node, all his children, from the tree
+ * @return true if the node was deleted (found), false otherwise
+ */
+ bool DeleteNodeById(unsigned int nId);
+
+ /**
+ * Method that removes all the pairs that contains the a given id is a given component position.
+ * @param vectorPairNodesScore the vector of pair nodes and score
+ * @param component is the component position to be analyzed
+ * @param nId is the searched id
+ * @return the number of deleted pairs
+ */
+ unsigned int DeleteEntriesByIdNode(Vector_Pair_PairNodes_Score& vectorPairNodesScore, unsigned int component, unsigned int nId);
+
+ /**
+ * Method that compares all the actual branches to all branches in the tree given by parameter
+ * @param treeB in the tree to be compared with
+ */
+ Vector_Pair_PairNodes_Score CompareAllBranches(AirwaysTree& treeB);
+
+ /**
+ * Method that return the composed edge from the root node to the node whose id is given by parameter
+ * @param idNode is the id of the node where the complex edge ends.
+ */
+ Edge* GetComposedEdge(unsigned int idNode);
+
+ /**
+ * Method that detaches a node and return a tree with a root and just the detaches node
+ * @param nId is the id node to be searched
+ * @return a tree with a root and just the detaches node
+ */
+ AirwaysTree* GetSingleDetachedTreeNodeById(unsigned int nId);
+
+ void UpdateLevels(unsigned int levels = 0);
+
+ void SubIsomorphism(AirwaysTree& treeB);
+
+ bool Isomorphism(AirwaysTree& tree);
+
+ void ImageReconstructionFromNode(TInputImage::Pointer& result, Node* node, bool skeleton = false);
+
+ void ImageReconstruction(TInputImage::Pointer& result, bool skeleton = false, Node* node = NULL);
+
+ void ImageReconstructionBySpheres(TInputImage::Pointer& result, Node* node = NULL, bool useMarks = false);
+
+
+ /**
+ * Method that prints the airway up to a given level
+ * @level is the level where the printing stops
+ */
+ void printUpToLevel(int level);
+
+ void printNodeAndChildrenIds();
+
+ /**
+ * Method that creates the quaternions of the airways up to a given level
+ * @level is the level where the creation stops.
+ * @pre level >= 2
+ */
+ void createQuaternionsUpToLevel(int level);
+
+ /**
+ * Method that prints the statistics of number of bifurcations in each node
+ */
+ void getStatisticsBifurcations();
+
+ void saveAirwayAsVTK(const std::string filename, bool common = false);
+
+ /**
+ * Method that save the airways in a image
+ * @param filename is the name of the image to be saved
+ * @param dims is the image dimensions
+ * @param spacing is the image spacing
+ * @param origin is the image origin
+ */
+ void saveAirwayToImage(std::string filename, int dims[3], double spc[3], double origin[3]);
+
+ /**
+ * Method that draws a line in an image
+ * @initVoxel is the initial voxel
+ * @endVoxel is the final voxel
+ * @imagePointer is the image pointer
+ */
+ void drawLineInImage(Vec3 initVoxel, Vec3 endVoxel, TInputImage::Pointer imagePointer);
+
+ void drawLine(int* pixelPointInit, int* pixelPointEnd, TInputImage::Pointer _imagePointer);
+
+
+protected:
+
+ Node* GetNode(const Vec3 nodeCoords, Node* current = NULL);
+
+ Node* GetLeaf(Node* current = NULL);
+
+ void UpdateEdges(Node* node = NULL);
+
+ AirwaysTree& GetSubTreeByDetachingNode(Node* node);
+
+ /**
+ * Method that adds an edge to the graph given by parameter
+ * @param source is the source vertex of the edge
+ * @param target is the target vertex of the edge
+ * @param graph is the graph to add the new edge
+ */
+ DiGraph_EdgePair AddBoostEdge(const pair_posVox_rad& source, const pair_posVox_rad& target, DiGraph& graph);
+
+ bool FindNeighborhood(DiGraph& graph, Vec3Queue& queue, Vec3List& used, const Vec3& end, ConstNeighborhoodIterator& iterator);
+
+ /**
+ * Method that return the shortest* path between two nodes, begin and end. *As we are dealing with arborescent trees
+ * then there is just one path between two nodes!
+ * The path goes from the end to the begin pixel!
+ * @param graph containing all the nodes and edges
+ * @param begin is the starting node
+ * @param end is the ending node
+ * @return the shortest path between the starting and ending nodes or vertices.
+ */
+ DiGraph GetGraphShortestPath(const DiGraph& graph, const Vec3& begin, const Vec3& end);
+
+ vec_pair_posVox_rad GetGraphShortestPath_AMP(const DiGraph& graph, const Vec3& begin, const Vec3& end);
+
+ vec_pair_posVox_rad GetSkeletonInfoFromEdge(const Vec3& begin, const Vec3& end);
+
+ bool ComputeSimpleRegionGrowing(TInputImage::Pointer& img, const Vec3 seedPtn);
+
+ void RemoveBranchFromImage(TInputImage::Pointer& img, Node* node);
+
+ void CreateSpheres(TInputImage::Pointer& img, Node* node);
+
+ void CalculateVTKLinesFromEdges(const Node* node, const unsigned int& parentId,
+ unsigned int& cId, vtkSmartPointer<vtkPoints>& pts,
+ vtkSmartPointer<vtkCellArray>& lines,
+ vtkSmartPointer<vtkUnsignedCharArray>& colors, bool common);
+ void FillTree(Node* node);
+
+protected:
+
+ std::string m_img_name;
+ std::string m_result_path;
+ TInputImage::Pointer m_skeleton;
+ TInputImage::Pointer m_img;
+ Node* m_root;
+ vec_nodes m_nodes;
+ EdgeVector m_edges;
+};
+} /* namespace airways */
+
+#endif /* AIRWAYSTREE_H_ */
--- /dev/null
+/*
+ * airwaysTreeTypeDefinition.h
+ *
+ * Created on: May 19, 2014
+ * Author: caceres
+ */
+
+#ifndef AIRWAYSTREETYPEDEFINITION_H_
+#define AIRWAYSTREETYPEDEFINITION_H_
+
+//std
+#include <vector>
+#include <queue>
+#include <cstddef>
+#include <utility>
+//boost
+#include <boost/graph/adjacency_list.hpp>
+//itk
+#include "itkImage.h"
+#include "itkConnectedComponentImageFilter.h"
+#include "itkConstNeighborhoodIterator.h"
+#include "itkImageDuplicator.h"
+#include "itkSphereSpatialFunction.h"
+#include "itkConnectedThresholdImageFilter.h"
+#include "itkSubtractImageFilter.h"
+#include "itkCastImageFilter.h"
+//airways
+#include "../MathLib/vec3.h"
+#include "airwaysNode.h"
+#include "airwaysEdge.h"
+
+namespace airways
+{
+ typedef std::pair<Vec3, double> pair_posVox_rad;
+ typedef std::vector<pair_posVox_rad> vec_pair_posVox_rad;
+
+ /**
+ * Pair of id nodes
+ */
+
+
+ /**
+ * Pair of a pair of nodes and a score
+ */
+ //typedef std::pair<pair_node_node,double> Pair_PairNodes_Score;
+
+ //----------------------------------------------------------------------------
+ /*
+ * Typedef declaration
+ */
+ /**
+ * @typedef Defines a graph using the Boost DiGraph library
+ */
+ typedef boost::adjacency_list<boost::vecS, boost::vecS, boost::directedS,
+ pair_posVox_rad> DiGraph;
+ //end typedef Graph
+ //--------------------------------------------------------------------------
+ /**
+ * @typedef Defines boost::vertex_descriptor of DiGraph
+ */
+ typedef boost::graph_traits<DiGraph>::vertex_descriptor DiGraph_VertexDescriptor;
+ //--------------------------------------------------------------------------
+ ///DiGraph edge_descriptor
+ typedef boost::graph_traits<DiGraph>::edge_descriptor DiGraph_EdgeDescriptor;
+ //--------------------------------------------------------------------------
+ ///DiGraph edge pair
+ typedef std::pair<DiGraph_EdgeDescriptor, bool> DiGraph_EdgePair;
+ //--------------------------------------------------------------------------
+ ///DiGraph Vertex IndexMap
+ typedef boost::property_map<DiGraph, boost::vertex_index_t>::type DiGraph_VertexIndexMap;
+ //--------------------------------------------------------------------------
+ ///DiGraph VertexIterator
+ typedef boost::graph_traits<DiGraph>::vertex_iterator DiGraph_VertexIterator;
+ //--------------------------------------------------------------------------
+ ///DiGraph Vertex iterator pair
+ typedef std::pair<DiGraph_VertexIterator, DiGraph_VertexIterator> DiGraph_VertexIteratorPair;
+ //--------------------------------------------------------------------------
+ ///DiGraph Edge iterator
+ typedef boost::graph_traits<DiGraph>::edge_iterator DiGraph_EdgeIterator;
+ //--------------------------------------------------------------------------
+ ///DiGraph In Edge iterator
+ typedef boost::graph_traits<DiGraph>::in_edge_iterator DiGraph_InEdgeIterator;
+ //--------------------------------------------------------------------------
+ ///DiGraph Out Edge iterator
+ typedef boost::graph_traits<DiGraph>::out_edge_iterator DiGraph_OutEdgeIterator;
+ //--------------------------------------------------------------------------
+ ///DiGraph Adjacency iterator
+ typedef boost::graph_traits<DiGraph>::adjacency_iterator DiGraph_AdjacencyIterator;
+ //--------------------------------------------------------------------------
+ ///DiGraph Vertex descriptor vector
+ typedef std::vector<DiGraph_VertexDescriptor> DiGraph_VertexVector;
+ //--------------------------------------------------------------------------
+ ///DiGraph Vertex descriptor queue
+ typedef std::queue<DiGraph_VertexDescriptor> DiGraph_VertexQueue;
+ //--------------------------------------------------------------------------
+ ///DiGraph Edge descriptor vector
+ typedef std::vector<DiGraph_EdgeDescriptor> DiGraph_EdgeVector;
+ //----------------------------------------------------------------------------
+ const unsigned int Dim = 3;
+ typedef double TPixel;
+ typedef unsigned char TPixel2;
+ typedef itk::Image<TPixel, Dim> TInputImage;
+ typedef itk::Image<TPixel2, Dim> TInputImage2;
+ typedef itk::Image<unsigned short, Dim> TInputImage3;
+ //----------------------------------------------------------------------------
+ typedef itk::ImageDuplicator<TInputImage> DuplicatorType;
+ //--------------------------------------------------------------------------
+ typedef itk::ConstNeighborhoodIterator<TInputImage> ConstNeighborhoodIterator;
+ typedef itk::ImageRegionIterator<TInputImage> RegionIterator;
+ //----------------------------------------------------------------------------
+ typedef itk::SphereSpatialFunction<> SphereType;
+ //----------------------------------------------------------------------------
+ typedef itk::ConnectedThresholdImageFilter<TInputImage, TInputImage> ConnectedFilterType;
+ //----------------------------------------------------------------------------
+ typedef itk::SubtractImageFilter<TInputImage, TInputImage> SubtractImageFilterType;
+ //----------------------------------------------------------------------------
+ typedef itk::CastImageFilter<TInputImage, TInputImage2> CastFilterType;
+ typedef itk::ConnectedComponentImageFilter<TInputImage2, TInputImage3> ConnectedComponentImageFilterType;
+ //----------------------------------------------------------------------------
+
+ //----------------------------------------------------------------------------
+ typedef std::vector<Edge*> EdgeVector;
+ //----------------------------------------------------------------------------
+ typedef std::queue<Vec3> Vec3Queue;
+ typedef std::vector<Vec3> Vec3List;
+//
+/*
+ * End of typedef declaration
+ */
+}
+
+#endif /* AIRWAYSTREETYPEDEFINITION_H_ */
--- /dev/null
+
+## Find Boost
+FIND_PACKAGE(Boost REQUIRED system filesystem program_options)
+INCLUDE_DIRECTORIES(${Boost_INCLUDE_DIRS})
+
+SUBDIRS(
+ MathLib
+ AirwaysLib
+ )
+INCLUDE_DIRECTORIES(
+ MathLib
+ AirwaysLib
+ )
+ADD_EXECUTABLE(TempAirwaysAppli TempAirwaysAppli.cxx)
+TARGET_LINK_LIBRARIES(
+ TempAirwaysAppli
+ ${Boost_LIBRARIES} MathLib AirwaysLib
+ cpPlugins
+ ${fpa_LIBRARIES}
+ )
+
+## eof - $RCSfile$
--- /dev/null
+## =============================
+## = Set names and directories =
+## =============================
+
+SET(lib_NAME MathLib)
+
+## ===============
+## = Source code =
+## ===============
+
+FILE(GLOB lib_HEADERS_H "*.h")
+FILE(GLOB lib_HEADERS_HPP "*.hpp")
+FILE(GLOB lib_HEADERS_HXX "*.hxx")
+FILE(GLOB lib_SOURCES_C "*.c")
+FILE(GLOB lib_SOURCES_CPP "*.cpp")
+FILE(GLOB lib_SOURCES_CXX "*.cxx")
+
+## =====================
+## = Compilation rules =
+## =====================
+
+ADD_LIBRARY(
+ ${lib_NAME}
+ SHARED
+ ${lib_SOURCES_C}
+ ${lib_SOURCES_CPP}
+ ${lib_SOURCES_CXX}
+ )
+GENERATE_EXPORT_HEADER(
+ ${lib_NAME}
+ BASE_NAME ${lib_NAME}
+ EXPORT_MACRO_NAME ${lib_NAME}_EXPORT
+ EXPORT_FILE_NAME ${CMAKE_CURRENT_BINARY_DIR}/${lib_NAME}_Export.h
+ STATIC_DEFINE ${lib_NAME}_BUILT_AS_STATIC
+ )
+##TARGET_LINK_LIBRARIES(${lib_NAME})
+
+## eof - $RCSfile$
--- /dev/null
+//----------------------------------------------------------------------------
+// Class definition include
+//----------------------------------------------------------------------------
+
+#include "Quaternion.h"
+
+//----------------------------------------------------------------------------
+// Class implementation
+//----------------------------------------------------------------------------
+//----------------------------------------------------------------------------
+// Builder
+//----------------------------------------------------------------------------
+
+Quaternion::Quaternion( )
+{
+ theta = 0.0;
+ w = Vec3(0,0,0);
+ r = cos(theta/2);
+ double sinHalfTheta = sin(theta/2);
+ im = w*sinHalfTheta;
+}
+
+Quaternion::Quaternion(Vec3 v1, Vec3 v2)
+{
+ // Normalize the vectors
+ v1.Normalize();
+ v2.Normalize();
+
+ // Dot product
+ float cosTheta = v1.Dot(v2);
+ theta = acos(cosTheta);
+
+ // Cross product
+ w = v1.Cross(v2);
+
+ // Normalization
+ w.Normalize();
+
+ // Build the quaternion
+ r = cos(theta/2);
+ double sinHalfTheta = sin(theta/2);
+ im = w*sinHalfTheta;
+}
+
+Quaternion::Quaternion(float nTheta, Vec3 nW)
+{
+ theta = nTheta;
+ w = nW;
+ w.Normalize();
+
+ // Build the quaternion
+ r = cos(theta/2);
+ double sinHalfTheta = sin(theta/2);
+ im = w*sinHalfTheta;
+}
+
+Quaternion::Quaternion(float nR, float nI, float nJ, float nK)
+{
+ r = nR;
+ im = Vec3(nI, nJ, nK);
+
+ theta = 2*acos(r);
+ double sinHalfTheta = sin(theta/2);
+ w = im/sinHalfTheta;
+}
+
+//----------------------------------------------------------------------------
+// Destructor
+//----------------------------------------------------------------------------
+
+Quaternion::~Quaternion( )
+{
+
+}
+//----------------------------------------------------------------------------
+// Methods
+//----------------------------------------------------------------------------
+
+float Quaternion::getR()
+{
+ return r;
+}
+
+float Quaternion::getI()
+{
+ return w[0];
+}
+
+float Quaternion::getJ()
+{
+ return w[1];
+}
+
+float Quaternion::getK()
+{
+ return w[2];
+}
+
+Vec3 Quaternion::getW()
+{
+ return w;
+}
+
+Quaternion Quaternion::conjugate( )
+{
+ return Quaternion(theta, -w);
+}
+
+Vec3 Quaternion::rotateVector(Vec3 vector)
+{
+ // Quaterion representing the vector
+ vector.Normalize();
+ Quaternion qv(0,vector[0],vector[1],vector[2]);
+ Quaternion qthis(r,im[0],im[1],im[2]);
+ Quaternion qthis_conj = qthis.conjugate();
+
+
+ cout << "Qv: [" << qv << "]" << endl;
+ cout << "Qthis: [" << qthis << "]" << endl;
+ cout << "QthisConj: [" << qthis_conj << "]" << endl;
+
+ // Final quaternion
+ Quaternion Qtemp = qthis*qv;
+ cout << "Qtemp: [" << Qtemp << "]" << endl;
+
+ Quaternion qf = Qtemp*qthis_conj;
+
+ cout << "Qf: [" << qf << "]" << endl;
+
+ return qf.im;
+ //return Vec3(0,0,0);
+}
+
+//----------------------------------------------------------------------------
+// Operators
+//----------------------------------------------------------------------------
+
+const float& Quaternion::operator[](unsigned int i) const
+{
+ return (i == 0 ? r : (im)[i-1]);
+}
+
+Quaternion Quaternion::operator*(Quaternion q)
+{
+ float a1 = this->r;
+ float b1 = this->w[0];
+ float c1 = this->w[1];
+ float d1 = this->w[2];
+
+ float a2 = q.r;
+ float b2 = q.w[0];
+ float c2 = q.w[1];
+ float d2 = q.w[2];
+
+ float r_n = a1*a2 - b1*b2 - c1*c2 - d1*d2;
+ float i_n = a1*b2 + a2*b1 + c1*d2 - c2*d1;
+ float j_n = a1*c2 - b1*d2 + a2*c1 + b2*d1;
+ float k_n = a1*d2 + b1*c2 - b2*c1 + a2*d1;
+
+ return Quaternion(r_n, i_n, j_n, k_n);
+}
+
+std::ostream& operator<<(std::ostream& os, const Quaternion& q)
+{
+ os << q[0] << " " << q[1] << " " << q[2] << " " << q[3];
+ os << " " <<q.theta*180/M_PI << " " << q.w;
+ return os;
+}
--- /dev/null
+/**
+ * Class that represents a Quaternion
+ * Author: Alfredo Morales Pinzón
+ * Date: October 7, 2014
+ */
+
+#ifndef __Quaternion_h__
+#define __Quaternion_h__
+
+//-----------------
+// C++
+//-----------------
+#include <math.h>
+
+//-----------------
+// OWN
+//-----------------
+#include <appli/TempAirwaysAppli/MathLib/MathLib_Export.h>
+#include "vec3.h"
+
+// Namespaces
+using namespace std;
+
+//------------------------------------------------------------------------------
+// Class that represents a Quaterion with all its operations
+//------------------------------------------------------------------------------
+class MathLib_EXPORT Quaternion
+{
+public:
+
+ //------------------------------------------------------------
+ //Builder
+ //------------------------------------------------------------
+
+ /*
+ * Class builder
+ */
+ Quaternion( );
+
+ /*
+ * Class builder using two vectors
+ * @param v1 is the first vector
+ * @param v2 is the second vector
+ */
+ Quaternion(Vec3 v1, Vec3 v2);
+
+ /*
+ * Class builder using theta and w vector
+ * @param nTheta is the theta angle
+ * @param nW is the w vector
+ */
+ Quaternion(float nTheta, Vec3 nW);
+
+ /*
+ * Class builder using r and im
+ * @param nR in the real part of the quaternion
+ * @param nI, nJ and nK are the imaginary parts of the quaternion
+ */
+ Quaternion(float nR, float nI, float nJ, float nK);
+
+ //------------------------------------------------------------
+ //Destructor
+ //------------------------------------------------------------
+
+ ~Quaternion();
+
+ //------------------------------------------------------------
+ //Public methods
+ //------------------------------------------------------------
+
+ /**
+ * Method that returns the real (r) component
+ * @return the r component
+ */
+ float getR();
+
+ /**
+ * Method that returns the i component
+ * @return the i component
+ */
+ float getI();
+
+ /**
+ * Method that returns the j component
+ * @return the j component
+ */
+ float getJ();
+
+ /**
+ * Method that returns the k component
+ * @return the k component
+ */
+ float getK();
+
+ /**
+ * Method that returns the w vector
+ * @return the w vector
+ */
+ Vec3 getW();
+
+ /**
+ * Method that returns the conjugate of the quaternion
+ * @return the conjugate of the quaternion
+ */
+ Quaternion conjugate( );
+
+ /**
+ * Method that returns a vector transformed by the quaternion
+ * @param vector to be rotated
+ * @return the transformed vector
+ */
+ Vec3 rotateVector(Vec3 vector);
+
+ //------------------------------------------------------------
+ //Operators
+ //------------------------------------------------------------
+
+ /*
+ * Definition of operator *
+ * @param is the quaternion to be multiplied with
+ * @return the multiplication quaternion
+ */
+ Quaternion operator*(Quaternion q);
+
+ /**
+ * Definition of operator []
+ * @param Position of the quanternion to be returned. Position 0 corresponds
+ * to the real component and 1, 2, and 3 to the complex components.
+ * @param The value on the specified position is returned.
+ *
+ */
+ const float& operator[](unsigned int i) const;
+
+ /**
+ * Operator to print
+ */
+ friend std::ostream& operator<<(std::ostream& os, const Quaternion& coord);
+
+
+private:
+
+ //------------------------------------------------------------
+ //Private methods
+ //------------------------------------------------------------
+
+
+ //------------------------------------------------------------
+ //Atributes
+ //------------------------------------------------------------
+
+ /**
+ * Real part of the quaternion
+ */
+ float r;
+
+ /**
+ * Angle between vectors
+ */
+ float theta;
+
+ /**
+ * Rotation vector
+ */
+ Vec3 w;
+
+ /**
+ * Vector the quaternion
+ */
+ Vec3 im;
+
+};
+
+//------------------------------------------------------------------------------
+#endif
--- /dev/null
+#include "vec3.h"
+
+
+Vec3::Vec3() : x(0.0), y(0.0), z(0.0)
+{
+}
+
+Vec3::Vec3(const float& x, const float& y, const float& z)
+{
+ this->x = x;
+ this->y = y;
+ this->z = z;
+}
+
+const float* Vec3::GetVec3() const
+{
+ float* ret = new float[3];
+ ret[0] = this->x;
+ ret[1] = this->y;
+ ret[2] = this->z;
+ return ret;
+}
+
+float* Vec3::GetVec3()
+{
+ float* ret = new float[3];
+ ret[0] = this->x;
+ ret[1] = this->y;
+ ret[2] = this->z;
+ return ret;
+}
+
+Vec3& Vec3::operator=(const Vec3& vec)
+{
+ this->x = vec.x;
+ this->y = vec.y;
+ this->z = vec.z;
+ return *this;
+}
+
+float& Vec3::operator[](unsigned int i)
+{
+ return (i == 0 ? x : i == 1 ? y : z);
+}
+
+const float& Vec3::operator[](unsigned int i) const
+{
+ return (i == 0 ? this->x : i == 1 ? this->y : this->z);
+}
+
+Vec3 Vec3::operator+() const
+{
+ return *this;
+}
+
+Vec3 Vec3::operator-() const
+{
+ return Vec3(-this->x, -this->y, -this->z);
+}
+
+Vec3& Vec3::operator+=(const Vec3& b)
+ {
+ this->x += b.x;
+ this->y += b.y;
+ this->z += b.z;
+ return *this;
+ }
+
+Vec3& Vec3::operator-=(const Vec3& b)
+ {
+ this->x -= b.x;
+ this->y -= b.y;
+ this->z -= b.z;
+ return *this;
+ }
+
+Vec3& Vec3::operator*=(const float& k)
+ {
+ this->x *= k;
+ this->y *= k;
+ this->z *= k;
+ return *this;
+ }
+
+Vec3& Vec3::operator/=(const float& k)
+ {
+ this->x /= k;
+ this->y /= k;
+ this->z /= k;
+ return *this;
+ }
+
+Vec3 Vec3::operator+(const Vec3& b)
+{
+ return Vec3(this->x + b.x, this->y + b.y, this->z + b.z);
+}
+
+Vec3 Vec3::operator-(const Vec3& b)
+{
+ return Vec3(this->x - b.x, this->y - b.y, this->z - b.z);
+}
+
+Vec3 Vec3::operator*(const float& k)
+{
+ return Vec3(this->x * k, this->y * k, this->z * k);
+}
+
+Vec3 Vec3::operator/(const float& k)
+{
+ return Vec3(this->x / k, this->y / k, this->z / k);
+}
+
+bool Vec3::operator==(const Vec3& b) const
+{
+ //double epsilon = std::numeric_limits<double>::epsilon();
+ //return (fabs(this->x - b.x) <= epsilon) && (fabs(this->y - b.y) <= epsilon)
+ //&& (fabs(this->z - b.z) <= epsilon);
+ return this->x == b.x && this->y == b.y && this->z == b.z;
+}
+
+bool Vec3::operator!=(const Vec3& b) const
+ {
+ return (!(*this == b));
+ }
+
+bool Vec3::operator<(const Vec3& b) const
+{
+ if(this->x < b.x)
+ return true;
+ else if(this->x > b.x)
+ return false;
+ else
+ {
+ if(this->y < b.y)
+ return true;
+ else if(this->y > b.y)
+ return false;
+ else
+ {
+ if(this->z < b.z)
+ return true;
+ return false;
+ }
+ }
+}
+
+bool Vec3::operator>(const Vec3& b)
+{
+ if(this->x > b.x)
+ return true;
+ else if(this->x < b.x)
+ return false;
+ else
+ {
+ if(this->y > b.y)
+ return true;
+ else if(this->y < b.y)
+ return false;
+ else
+ {
+ if(this->z > b.z)
+ return true;
+ return false;
+ }
+ }
+}
+
+float Vec3::Norm() const
+{
+ return sqrt(this->x * this->x + this->y * this->y + this->z * this->z);
+}
+
+void Vec3::Normalize()
+{
+ float norm = sqrt(this->x * this->x + this->y * this->y + this->z * this->z);
+ this->x = this->x / norm;
+ this->y = this->y / norm;
+ this->z = this->z / norm;
+}
+
+void Vec3::set(const float& x, const float& y, const float& z)
+{
+ this->x = x;
+ this->y = y;
+ this->z = z;
+}
+
+float Vec3::Dot(const Vec3& b) const
+{
+ return (this->x * b.x + this->y * b.y + this->z * b.z);
+}
+
+Vec3 Vec3::Orthogonal(const Vec3& b) const
+{
+ Vec3 u = Vec3(fabs(b.x), fabs(b.y), fabs(b.z));
+ int i = 0;
+ int j = 1;
+ if ((u.x > u.y) && (u.z > u.y))
+ j = 2;
+ else
+ {
+ i = 1;
+ j = 2;
+ if (u.x > u.z)
+ j = 0;
+ }
+ u = Vec3();
+ u[i] = b[j];
+ u[j] = -b[i];
+ return u;
+}
+
+Vec3 Vec3::Cross(const Vec3& b) const
+{
+ return Vec3(this->y * b.z - this->z * b.y, this->z * b.x - this->x * b.z,
+ this->x * b.y - this->y * b.x);
+}
+
+
+float Vec3::EculideanDistance(const Vec3& b) const
+{
+ return sqrt(
+ pow(this->x - b.x, 2.0) + pow(this->y - b.y, 2.0)
+ + pow(this->z - b.z, 2.0));
+}
+
+std::ostream&
+operator<<(std::ostream& os, const Vec3& coord)
+{
+ os << coord[0] << " " << coord[1] << " " << coord[2];
+ return os;
+}
+
--- /dev/null
+/**
+ * Class that models a vector in 3D
+ * Author: Diego Cáceres
+ * Modifications by: Alfredo Morales Pinzón
+ * Date modifications: Octobre 7, 2014
+ */
+#ifndef VEC3_H
+#define VEC3_H
+
+#include <appli/TempAirwaysAppli/MathLib/MathLib_Export.h>
+#include <cmath>
+#include <iostream>
+#include <limits> // std::numeric_limits
+
+//------------------------------------------------------------------------------
+/**
+ * @brief The Vec3 class that represents a vector in 3D
+ */
+class MathLib_EXPORT Vec3
+{
+private:
+ float x, y, z; //!< Coordinates
+public:
+ //----------------------------------------------------------------------------
+ /**
+ * @brief Vec3
+ */
+ Vec3();
+ //----------------------------------------------------------------------------
+ /**
+ * @brief Vec3
+ * @param x
+ * @param y
+ * @param z
+ */
+ Vec3(const float& x, const float& y, const float& z);
+ //----------------------------------------------------------------------------
+ //Getters
+ //----------------------------------------------------------------------------
+ /**
+ * @brief GetVec3
+ * @return
+ */
+ const float* GetVec3() const;
+ //----------------------------------------------------------------------------
+ /**
+ * @brief GetVec3
+ * @return
+ */
+ float* GetVec3();
+ //----------------------------------------------------------------------------
+ /**
+ * @brief Norm
+ * @return
+ */
+ float Norm() const;
+
+ /**
+ * Method that normalizes the vector. v = v/|v|
+ */
+ void Normalize();
+
+ /**
+ * Method that sets all the values
+ */
+ void set(const float& x, const float& y, const float& z);
+
+ /**
+ * @brief Dot
+ * @param b
+ * @return
+ */
+ float Dot(const Vec3& b) const;
+ //----------------------------------------------------------------------------
+ /**
+ * @brief Orthogonal
+ * @param b
+ * @return
+ */
+ Vec3 Orthogonal(const Vec3& b) const;
+ //----------------------------------------------------------------------------
+ /**
+ * @brief Cross
+ * @param b
+ * @return
+ */
+ Vec3 Cross(const Vec3& b) const;
+ //----------------------------------------------------------------------------
+ /**
+ * @brief EculideanDistance between two points
+ * @param b
+ * @return
+ */
+ float EculideanDistance(const Vec3& b) const;
+ //----------------------------------------------------------------------------
+ //Operator Overloads
+ //----------------------------------------------------------------------------
+ /**
+ * @brief The copy operator
+ * @param vec
+ * @return
+ */
+ Vec3& operator=(const Vec3& vec);
+ //----------------------------------------------------------------------------
+ /**
+ * @brief operator []
+ * @param i
+ * @return
+ */
+ float& operator[](unsigned int i);
+ //----------------------------------------------------------------------------
+ /**
+ * @brief operator []
+ * @param i
+ * @return
+ */
+ const float& operator[](unsigned int i) const;
+ //----------------------------------------------------------------------------
+ /**
+ * @brief operator +
+ * @return
+ */
+ Vec3 operator+() const;
+ //----------------------------------------------------------------------------
+ /**
+ * @brief operator -
+ * @return
+ */
+ Vec3 operator-() const;
+ //----------------------------------------------------------------------------
+ /**
+ * @brief operator +=
+ * @param b
+ * @return
+ */
+ Vec3& operator+=(const Vec3& b);
+ //----------------------------------------------------------------------------
+ /**
+ * @brief operator -=
+ * @param b
+ * @return
+ */
+ Vec3& operator-=(const Vec3& b);
+ //----------------------------------------------------------------------------
+ /**
+ * @brief operator *=
+ * @param k
+ * @return
+ */
+ Vec3& operator*=(const float& k);
+ //----------------------------------------------------------------------------
+ /**
+ * @brief operator /=
+ * @param k
+ * @return
+ */
+ Vec3& operator/=(const float& k);
+ //----------------------------------------------------------------------------
+ /**
+ * @brief operator +
+ * @param b
+ * @return
+ */
+ Vec3 operator+(const Vec3& b);
+ //----------------------------------------------------------------------------
+ /**
+ * @brief operator -
+ * @param b
+ * @return
+ */
+ Vec3 operator-(const Vec3& b);
+ //----------------------------------------------------------------------------
+ /**
+ * @brief operator *
+ * @param k
+ * @return
+ */
+ Vec3 operator*(const float& k);
+ //----------------------------------------------------------------------------
+ /**
+ * @brief operator /
+ * @param k
+ * @return
+ */
+ Vec3 operator/(const float& k);
+ //----------------------------------------------------------------------------
+ /**
+ * @brief operator ==
+ * @param b
+ * @return
+ */
+ bool operator==(const Vec3& b) const;
+ //----------------------------------------------------------------------------
+ /**
+ * @brief operator !=
+ * @param b
+ * @return
+ */
+ bool operator!=(const Vec3& b) const;
+ //----------------------------------------------------------------------------
+ /**
+ * @brief operator <
+ * @param b
+ * @return
+ */
+ bool operator<(const Vec3& b) const;
+ //----------------------------------------------------------------------------
+ /**
+ * @brief operator >
+ * @param b
+ * @return
+ */
+ bool operator>(const Vec3& b);
+ //----------------------------------------------------------------------------
+ friend std::ostream& operator<<(std::ostream& os, const Vec3& coord);
+ //----------------------------------------------------------------------------
+
+};
+//------------------------------------------------------------------------------
+
+#endif // VEC3_H
--- /dev/null
+#include <iostream>
+#include <cstdlib>
+#include <sstream>
+
+#include <boost/filesystem.hpp>
+#include <boost/program_options.hpp>
+
+#include "vec3.h"
+#include "airwaysTree.h"
+
+#include <cpPlugins/Interface.h>
+#include <cpPlugins/Workspace.h>
+
+using namespace airways;
+namespace po = boost::program_options;
+
+namespace
+{
+ const size_t ERROR_IN_COMMAND_LINE = 1;
+ const size_t SUCCESS = 0;
+ const size_t ERROR_UNHANDLED_EXCEPTION = 2;
+} // namespace
+
+typedef std::vector< AirwaysTree > AirwaysVector;
+
+// Auxiliar struct to save info for execution.
+typedef void* TImagePointer;
+struct TreeInfo{
+ TImagePointer image;
+ Vec3 seed;
+ std::string folderpath_pigResults;
+ std::string pig_name;
+ std::string image_name;
+
+};
+// Constants
+unsigned char red[3] = { 255, 0, 0 };
+unsigned char green[3] = { 0, 255, 0 };
+unsigned char blue[3] = { 0, 0, 255 };
+unsigned char yellow[3]= { 255, 255, 0 };
+unsigned char purple[3]= { 255, 0, 255 };
+unsigned char cyan[3]= { 0, 255, 255 };
+
+cpPlugins::Interface myPlugins;
+cpPlugins::Workspace myWorkspace;
+
+// -------------------------------------------------------------------------
+void Load_cpPlugins( const std::string& plugins, const std::string& ws );
+void CreateResultDirectory(AirwaysTree tree);
+void DrawVTKLinesFromTree(AirwaysTree& tree, const std::string filename, bool common);
+AirwaysTree& CreateAirwaysTreeFromSegmentation(Vec3 seed, TImagePointer input_image);
+vector<TreeInfo> ReadInputFile(const char* filename);
+void printCommonTreeBetweenTwoTrees(AirwaysTree tree_A, AirwaysTree tree_B, std::vector< std::pair< std::pair<Node*, Node*>, std::pair<Node*, Node*> > > vector_pair_edges_A_to_B, unsigned int Q, unsigned int F);
+void printMatchingResultToFile(AirwaysTree tree_A, AirwaysTree tree_B, std::map< unsigned int, std::vector<Node*> > map_A_to_B, std::map< unsigned int, std::vector<Node*> > map_B_to_A, unsigned int Q, unsigned int F);
+void createLinesAndPointsForVTK(const Node* node, vtkSmartPointer<vtkPoints>& pts,vtkSmartPointer<vtkCellArray>& lines, vtkSmartPointer<vtkUnsignedCharArray>& colors, bool common, bool isRoot);
+
+// -------------------------------------------------------------------------
+int main( int argc, char* argv[] )
+{
+ try
+ {
+ // Define and parse the program options
+ po::options_description desc("Options");
+ desc.add_options()("use_file", po::value<std::string>(), "Adds the filepath containing the file to be used in creaAirwaysTree -- Mandatory")
+ ("build_trees","Creates trees from a given filepath (arg) and stores it in a .vtk file")
+ ("subtree_levels", po::value<int>(),"Get a subtree(s) by levels - result: a vtk and reconstructed img .mhd")
+ ("subtree_length", po::value<float>(),"Get a subtree(s) by length - result: a vtk and reconstructed img .mhd")
+ ("subtree_diameter", po::value<float>(),"Get a subtree(s) by diameter - result: a vtk and reconstructed img .mhd")
+ ("compare_trees", "Compare the trees given in the input file or the subtrees using an option - result: a vtk and reconstructed img .mhd");
+ //TODO: Fix image segmentation. ("segment_images", "Creates a binary image corresponding to the segmentation of of a given original image (arg) and seed (arg), saves it to a .mhd file and uses it for subsequent operations")
+ po::variables_map vm;
+ try
+ {
+ std::string fileName;
+ vector<TreeInfo> treeInfoVector;
+ AirwaysVector aVector;
+ TImagePointer segmentationImage;
+ Vec3 seed;
+ po::store(po::parse_command_line(argc, argv, desc), vm); // can throw
+ if (vm.count("use_file"))
+ {
+ fileName = vm["use_file"].as<std::string>();
+ treeInfoVector = ReadInputFile(fileName.c_str());
+ }
+
+ if(vm.count("segment_images"))
+ {
+ }
+
+ // Build trees option
+ if (vm.count("build_trees"))
+ {
+ for(unsigned int i = 0; i < treeInfoVector.size(); ++ i){
+ TreeInfo info = treeInfoVector[i];
+ AirwaysTree tree = CreateAirwaysTreeFromSegmentation(info.seed, info.image);
+ tree.SetResultPath(info.folderpath_pigResults);
+ tree.SetImageName(info.image_name);
+ aVector.push_back(tree);
+ }
+ for (unsigned int i = 0; i < aVector.size(); ++i)
+ {
+ CreateResultDirectory(aVector[i]);
+ std::string fullPath = aVector[i].GetResultPath() + "/" + aVector[i].GetImageName() + ".vtk";
+ DrawVTKLinesFromTree(aVector[i], fullPath, false);
+
+ } //for
+ } //if
+
+ // Subtree levels option
+ if (vm.count("subtree_levels"))
+ {
+ } //if
+ else if (vm.count("subtree_length"))
+ {
+ } //if
+ else if (vm.count("subtree_diameter"))
+ {
+ } //if
+
+ if (vm.count("compare_trees"))
+ {
+ std::cout << "Option: Compare trees" << std::endl;
+ /**
+ * this piece of code only test in the case of two trees, so it can be replaced
+ * to a code which does a cascade
+ */
+
+ //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ //CODIGO DE COMPARACIÓN DE ARBOLES
+ //std::cout << "Comparing trees ..." << std::endl;
+ //aVector[0].CompareTrees(aVector[1]);
+ //std::cout << "Comparing trees ... OK" << std::endl;
+
+ std::cout << "Tree A ... " << std::endl;
+ //aVector[0].printNodeAndChildrenIds();
+ std::cout << "Tree A ... OK" << std::endl;
+
+ std::cout << "Tree B ... " << std::endl;
+ //aVector[1].printNodeAndChildrenIds();
+ std::cout << "Tree B ... OK" << std::endl;
+
+ std::cout << "Comparing trees Orkisz-Morales..." << std::endl;
+ // Vectors to save the common and uncommon nodes
+ vec_nodes commonA;
+ vec_nodes commonB;
+ vec_nodes nonCommonA;
+ vec_nodes nonCommonB;
+ std::map< unsigned int, std::vector<Node*> > map_A_to_B;
+ std::map< unsigned int, std::vector<Node*> > map_B_to_A;
+ std::vector< std::pair< std::pair<Node*, Node*>, std::pair<Node*, Node*> > > vector_pair_edges_A_to_B;
+
+ // Input parameter for comparison
+ unsigned int Q = 1; // Corresponds to the depth to select "fathers" nodes
+ unsigned int F = 1; // Correspond to the depth to select "family" nodes
+ std::cout << "Q = " << Q << std::endl;
+ std::cout << "F = " << F << std::endl;
+
+ clock_t before_compare = clock();
+ aVector[0].CompareTreesOrkiszMorales(aVector[1], Q, F, commonA, commonB, nonCommonA, nonCommonB, map_A_to_B, map_B_to_A, vector_pair_edges_A_to_B);
+ clock_t after_compare = clock();
+ double compare_time = double(after_compare - before_compare) / CLOCKS_PER_SEC;
+ std::cout << "Matching time: " << compare_time << std::endl;
+
+ // Print the common tree with common edges
+ printCommonTreeBetweenTwoTrees(aVector[0], aVector[1], vector_pair_edges_A_to_B, Q, F);
+ printMatchingResultToFile(aVector[0], aVector[1], map_A_to_B, map_B_to_A, Q, F);
+ //printSeparatedCommonTrees(aVector[0], aVector[1], vector_pair_edges_A_to_B, Q, F);
+ //reconstructAndSaveCommonTrees(aVector[0], aVector[1], vector_pair_edges_A_to_B, Q, F);
+
+ // Print all the maps that have more that two connections
+ std::cout << "Printing multiple relation A to B. Total relations: " << map_A_to_B.size() << std::endl;
+ for(std::map< unsigned int, std::vector<Node*> >::iterator it_map_A_to_B = map_A_to_B.begin(); it_map_A_to_B != map_A_to_B.end(); ++it_map_A_to_B)
+ {
+ if((*it_map_A_to_B).second.size() > 1)
+ {
+ std::cout << "Multiple relation A to B with size:" << (*it_map_A_to_B).second.size() << ", from Id:" << (*it_map_A_to_B).first << std::endl;
+ for(std::vector<Node*>::iterator it_node = (*it_map_A_to_B).second.begin(); it_node != (*it_map_A_to_B).second.end(); ++it_node)
+ {
+ std::cout << "Id: " << (*it_node)->GetId() << std::endl;
+ }
+ }
+ }
+ std::cout << "Printing multiple relation A to B ... OK" << std::endl;
+
+ std::cout << "Printing multiple relation B to A. Total relations: " << map_B_to_A.size() << std::endl;
+ for(std::map< unsigned int, std::vector<Node*> >::iterator it_map_B_to_A = map_B_to_A.begin(); it_map_B_to_A != map_B_to_A.end(); ++it_map_B_to_A)
+ {
+ if((*it_map_B_to_A).second.size() > 1)
+ {
+ std::cout << "Multiple relation B to A with size:" << (*it_map_B_to_A).second.size() << ", from Id:" << (*it_map_B_to_A).first << std::endl;
+ for(std::vector<Node*>::iterator it_node = (*it_map_B_to_A).second.begin(); it_node != (*it_map_B_to_A).second.end(); ++it_node)
+ {
+ std::cout << "Id: " << (*it_node)->GetId() << std::endl;
+ }
+ }
+ }
+ std::cout << "Printing multiple relation B to A ... OK" << std::endl;
+
+ // Mark only the common nodes
+ aVector[0].UnMarkAll();
+ aVector[1].UnMarkAll();
+
+ // --------------------------------------
+ // ------ Get common paths marked -------
+ // --------------------------------------
+
+ /*
+ std::string path_folder = "/run/media/alfredo/Data/Pulmones/results/airwaysGraphs/comparisonResutls/probes/";
+ std::string suffix_vtk = ".vtk";
+ std::string name_tree = "TreeA_dP2P";
+ int iteration = 1;
+ for(vec_nodes::iterator it_A = commonA.begin(); it_A != commonA.end(); ++it_A)
+ {
+ aVector[0].MarkPathFromNodeToNode(aVector[0].GetRoot(), (*it_A));
+
+ std::stringstream filepath_actualIteration;
+ filepath_actualIteration << path_folder << name_tree << "__" << iteration << "_idA_" << (*it_A)->GetId() << suffix_vtk;
+
+ DrawVTKLinesFromTree(aVector[0], filepath_actualIteration.str(), true);
+ ++iteration;
+ }
+
+ iteration = 1;
+ name_tree = "TreeB_dP2P";
+ for(vec_nodes::iterator it_B = commonB.begin(); it_B != commonB.end(); ++it_B)
+ {
+ aVector[1].MarkPathFromNodeToNode(aVector[1].GetRoot(), (*it_B));
+
+ std::stringstream filepath_actualIteration;
+ filepath_actualIteration << path_folder << name_tree << "__" << iteration << "_idB_" << (*it_B)->GetId() << suffix_vtk;
+
+ DrawVTKLinesFromTree(aVector[1], filepath_actualIteration.str(), true);
+ ++iteration;
+ }
+ */
+
+ // --------------------------------------
+ // --------------------------------------
+
+ /*
+ //XXXXXXXXXXXXXXXXXXXXXXXXXXXX
+ //1. PRINT EACH NODE OF EACH TREE AND SAVE IT USING AS NAME ITS ID
+
+ // Tree A
+ for(int i = 2; i < aVector[0].GetWeight(); ++i)
+ {
+ //std::cout << "Writing idA: " << i << std::endl;
+ AirwaysTree* tree_id = aVector[0].GetSingleDetachedTreeNodeById(i);
+ //std::cout << "Got id: " << i << ", numberOfNodes: " << tree_id->GetWeight() << std::endl;
+ if(tree_id)
+ {
+ std::stringstream filepath_actualIteration;
+ filepath_actualIteration << aVector[0].GetResultPath() << "/" << aVector[0].GetImageName() << "_id_" << i << suffix_vtk;
+ DrawVTKLinesFromTree(*tree_id, filepath_actualIteration.str(), false);
+ delete(tree_id);
+ }
+ else
+ std::cout << "Tree NULL" << std::endl;
+ }
+
+ // Tree B
+ for(int i = 2; i < aVector[1].GetWeight(); ++i)
+ {
+ //std::cout << "Writing idA: " << i << std::endl;
+ AirwaysTree* tree_id = aVector[1].GetSingleDetachedTreeNodeById(i);
+ //std::cout << "Got id: " << i << ", numberOfNodes: " << tree_id->GetWeight() << std::endl;
+ if(tree_id)
+ {
+ std::stringstream filepath_actualIteration;
+ filepath_actualIteration << aVector[1].GetResultPath() << "/" << aVector[1].GetImageName() << "_id_" << i << suffix_vtk;
+ DrawVTKLinesFromTree(*tree_id, filepath_actualIteration.str(), false);
+ delete(tree_id);
+ }
+ else
+ std::cout << "Tree NULL" << std::endl;
+ }
+
+ std::cout << "Comparing trees Orkisz-Morales... OK" << std::endl;
+
+ */
+
+ //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ /*// AMP
+ aVector[0].SubIsomorphism(aVector[1]);
+
+ std::cout << "Flag after SubIsomorphism" << std::endl;
+ std::string fullPathA = aVector[0].GetResultPath() + "/" + aVector[0].GetImageName() + "_CMP.vtk";
+ std::string fullPathB = aVector[1].GetResultPath() + "/" + aVector[1].GetImageName() + "_CMP.vtk";
+ DrawVTKLinesFromTree(aVector[0], fullPathA, true);
+ DrawVTKLinesFromTree(aVector[1], fullPathB, true);
+
+ std::cout << "Flag after write vtk" << std::endl;
+ TInputImage::Pointer imgA, imgB;
+ aVector[0].ImageReconstruction(imgA);
+ std::cout << "Flag after Reconstruction A" << std::endl;
+
+ aVector[1].ImageReconstruction(imgB);
+ std::cout << "Flag after Reconstruction B" << std::endl;
+
+ std::string fullPathA_mhd = aVector[0].GetResultPath() + "/" + aVector[0].GetImageName() + "_CMP.vtk";
+ std::string fullPathB_mhd = aVector[1].GetResultPath() + "/" + aVector[1].GetImageName() + "_CMP.vtk";
+ WriteImage(imgA, fullPathA_mhd);
+ WriteImage(imgB, fullPathB_mhd);
+
+ //AMP*/
+ }
+ po::notify(vm); // throws on error, so do after help in case
+ }
+ catch (std::exception& e)
+ {
+ std::cerr << "Unhandled Exception reached the top of main: " << e.what() << ", application will now exit" << std::endl;
+ return ERROR_UNHANDLED_EXCEPTION;
+ } // yrt
+ }
+ catch (std::exception& e)
+ {
+ std::cerr << "Unhandled Exception reached the top of main: " << e.what() << ", application will now exit" << std::endl;
+ return ERROR_UNHANDLED_EXCEPTION;
+ } // yrt
+ return( 0 );
+}
+
+// -------------------------------------------------------------------------
+void Load_cpPlugins( const std::string& plugins, const std::string& ws )
+{
+ myPlugins.LoadConfiguration( plugins );
+ myWorkspace.SetInterface( &myPlugins );
+ std::string err = myWorkspace.LoadWorkspace( ws );
+ if( err != "" )
+ {
+ std::cerr << "Error: " << err << std::endl;
+ std::exit( 1 );
+
+ } // fi
+}
+
+// -------------------------------------------------------------------------
+void CreateResultDirectory(AirwaysTree tree)
+{
+ boost::filesystem::path dir(tree.GetResultPath());
+ if (boost::filesystem::create_directories(dir))
+ {
+ std::cout << "FolderName = " << tree.GetResultPath() << " has been created" << std::endl;
+ } //if
+ else
+ std::cout << "FolderName = " << tree.GetResultPath() << " already exists" << std::endl;
+}
+
+// -------------------------------------------------------------------------
+void DrawVTKLinesFromTree(AirwaysTree& tree, const std::string filepath, bool common)
+{
+ // Create the array of points, lines, and colors
+ vtkSmartPointer<vtkPoints> pts = vtkSmartPointer<vtkPoints>::New();
+ vtkSmartPointer<vtkCellArray> lines = vtkSmartPointer<vtkCellArray>::New();
+ vtkSmartPointer<vtkUnsignedCharArray> colors = vtkSmartPointer<vtkUnsignedCharArray>::New();
+ colors->SetNumberOfComponents(3);
+ colors->SetName("Colors");
+
+ //vtk sphere for sources
+ //vtkSmartPointer<vtkSphereSource> sphereSource = vtkSmartPointer<vtkSphereSource>::New();
+ //Vec3 root = tree.GetRoot()->GetCoords();
+ //sphereSource->SetCenter(root[0], root[1], root[2]);
+ //sphereSource->SetRadius(1);
+ //UndirectedGraph
+
+ srand(time(NULL));
+ unsigned int id = 1;
+
+ // Create and fill the points, lines, and color vectors
+ //CalculateVTKLinesFromEdges(tree.GetRoot(), 0, id, pts, lines, colors, common);
+ pts->SetNumberOfPoints(tree.GetWeight()+1);
+ createLinesAndPointsForVTK(tree.GetRoot(), pts, lines, colors, common, true);
+
+ // Create the polydata
+ vtkSmartPointer<vtkPolyData> linesPolyData = vtkSmartPointer<vtkPolyData>::New();
+
+ // Set the points, lines, and colors
+ linesPolyData->SetPoints(pts);
+ linesPolyData->SetLines(lines);
+ linesPolyData->GetCellData()->SetScalars(colors);
+
+ // Write the file
+ vtkSmartPointer<vtkPolyDataWriter> writer = vtkSmartPointer<vtkPolyDataWriter>::New();
+ writer->SetFileName(filepath.c_str());
+
+#if VTK_MAJOR_VERSION <= 5
+ writer->SetInput(linesPolyData);
+#else
+ writer->SetInputData(linesPolyData);
+#endif
+ writer->Write();
+
+ //The following code is to test the vtk polydata and visualize it
+ /* // Visualize
+ vtkSmartPointer<vtkPolyDataMapper> mapper =
+ vtkSmartPointer<vtkPolyDataMapper>::New();
+ #if VTK_MAJOR_VERSION <= 5
+ mapper->SetInput(linesPolyData);
+ #else
+ mapper->SetInputData(linesPolyData);
+ #endif
+ vtkSmartPointer<vtkActor> actor = vtkSmartPointer<vtkActor>::New();
+ actor->SetMapper(mapper);
+
+ //sphere
+ vtkSmartPointer<vtkPolyDataMapper> mapperSphere = vtkSmartPointer<
+ vtkPolyDataMapper>::New();
+ mapperSphere->SetInputConnection(sphereSource->GetOutputPort());
+
+ vtkSmartPointer<vtkActor> actorSphere = vtkSmartPointer<vtkActor>::New();
+ actorSphere->SetMapper(mapperSphere);
+ //end sphere
+
+ vtkSmartPointer<vtkRenderer> renderer = vtkSmartPointer<vtkRenderer>::New();
+ vtkSmartPointer<vtkRenderWindow> renderWindow = vtkSmartPointer<
+ vtkRenderWindow>::New();
+ renderWindow->AddRenderer(renderer);
+ vtkSmartPointer<vtkRenderWindowInteractor> renderWindowInteractor =
+ vtkSmartPointer<vtkRenderWindowInteractor>::New();
+ renderWindowInteractor->SetRenderWindow(renderWindow);
+ renderer->AddActor(actorSphere);
+ renderer->AddActor(actor);
+
+ renderWindow->Render();
+ renderWindowInteractor->Start();*/
+}
+
+// -------------------------------------------------------------------------
+#include <fpa/Base/ImageSkeleton.h>
+#include <fpa/Image/MinimumSpanningTree.h>
+#include <cpExtensions/DataStructures/ImageIndexesContainer.h>
+#include <itkImage.h>
+
+template< class TImage, class TVertex >
+Node* FAVertexToNode( TVertex vertex, TImage* image )
+{
+ //The FrontAlgorithms Vertex is an TImageType::IndexType
+ typename TImage::PointType point;
+ image->TransformIndexToPhysicalPoint( vertex, point );
+ Vec3 alfPoint(point[0],point[1],point[2]);
+ return new Node(alfPoint);
+}
+
+AirwaysTree& ConvertFilterToAirwaysTree( )
+{
+ /*
+ typedef FilterType TFilter;
+ typedef typename TFilter::TVertex TVertexFA;
+ typedef typename TFilter::TVertices TVerticesFA;
+ typedef typename TFilter::TUniqueVertices TUniqueVerticesFA;
+ typedef typename TFilter::TVertexCompare TVertexCompareFA;
+ typedef typename TFilter::TBranches TBranchesFA;
+ typedef typename TFilter::TMinimumSpanningTree TMst;
+ */
+ typedef
+ fpa::Base::ImageSkeleton< fpa::Image::MinimumSpanningTree< 3 > >
+ TFASkeleton;
+ typedef TFASkeleton::TVertex TVertexFA;
+ typedef TFASkeleton::TVertexCmp TVertexCompareFA;
+ typedef cpExtensions::DataStructures::ImageIndexesContainer< 3 > TVerticesFA;
+
+ typedef Node TVertexAirways;
+ typedef Edge TEdgeAirways;
+ typedef pair_posVox_rad TSkelePoint;
+ typedef vec_pair_posVox_rad TSkelePoints;
+ typedef std::map< TVertexFA, TVertexAirways*, TVertexCompareFA > VertexMap;
+ VertexMap vertexMap;
+ std::time_t start,end;
+ std::time(&start);
+ std::cout << "Starting conversion " << std::endl;
+
+ auto w_filter = myWorkspace.GetFilter( "eb" );
+ auto branches = w_filter->GetOutputData( "Skeleton" )->GetITK< TFASkeleton >( );
+ auto& endpoints = w_filter->GetOutputData( "EndPoints" )->GetITK< TVerticesFA >( )->Get( );
+ auto& bifurcations = w_filter->GetOutputData( "Bifurcations" )->GetITK< TVerticesFA >( )->Get( );
+ auto image = myWorkspace.GetFilter( "reader" )->GetOutputData( "Output" )->GetITK< itk::ImageBase< 3 > >( );
+ auto distance_map = myWorkspace.GetFilter( "dmap" )->GetOutputData( "Output" )->GetITK< itk::Image< float, 3 > >( );
+
+ /*
+ TBranchesFA* branches = filter->GetBranches( );
+ TMst* mst = filter->GetMinimumSpanningTree();
+ */
+ int leoTreeWeigth = endpoints.size()+bifurcations.size()+1;
+ std::cout<< "Creates FA Tree with : "<<leoTreeWeigth << " nodes "<<std::endl;
+
+ auto seed0 = myWorkspace.GetFilter( "seed" )->GetOutputData( "Output" )->GetITK< TVerticesFA >( )->Get( )[ 0 ];
+ //Fill vertex map. Gotta do it with bifurcations, endpoints and the root.
+ //Do the root
+ vertexMap[ seed0 ] = FAVertexToNode( seed0, image );
+ //Do Endpoints
+ auto eIt = endpoints.begin();
+ for( ; eIt != endpoints.end(); ++eIt )
+ {
+ vertexMap[*eIt]=FAVertexToNode( *eIt,image );
+ }
+
+ auto biIt = bifurcations.begin();
+ for(; biIt != bifurcations.end(); ++biIt)
+ {
+ vertexMap[*biIt]=FAVertexToNode(*biIt,image);
+ }
+
+ //Now navigate branches to make the edges of the tree.
+ auto bIt = branches->Get( ).begin();
+ for(; bIt!=branches->Get( ).end(); ++bIt)
+ {
+ auto dVertex = bIt->first;
+ TVertexAirways* destination = vertexMap[dVertex];
+ if( destination == NULL )
+ {
+ vertexMap[dVertex]=FAVertexToNode(dVertex,image);
+ destination = vertexMap[dVertex];
+ }
+ auto brIt = bIt->second.begin( );
+ for( ; brIt != bIt->second.end( ); ++brIt )
+ {
+ auto sVertex = brIt->first;
+ TVertexAirways* source = vertexMap[sVertex];
+ if( source == NULL )
+ {
+ vertexMap[sVertex]=FAVertexToNode(sVertex,image);
+ source = vertexMap[sVertex];
+ }
+ destination->SetFather(source);
+ source->AddChild(destination);
+ TEdgeAirways* edge = new Edge();
+ edge->SetSource(source);
+ edge->SetTarget(destination);
+ //Update path info for the edge. This is why we don't need to call UpdateEdges on the constructor of the tree.
+
+ auto edgePath = brIt->second->GetVertexList( );
+ for( unsigned int pIt = 0; pIt < edgePath->Size( ); ++pIt )
+ {
+ auto cidx = edgePath->GetElement(pIt);
+ itk::ImageBase< 3 >::PointType pnt;
+ image->TransformContinuousIndexToPhysicalPoint(cidx,pnt);
+ TVertexFA idx;
+ image->TransformPhysicalPointToIndex(pnt,idx);
+ Vec3 coords = FAVertexToNode(idx,image)->GetCoords();
+ pair_posVox_rad skPair(coords,distance_map->GetPixel(idx));
+ edge->AddSkeletonPairInfo(skPair);
+ }
+ destination->SetEdge(edge);
+ }
+ }
+ AirwaysTree* tree = new AirwaysTree(dynamic_cast< TInputImage* >( image ), NULL, vertexMap[seed0], false);
+ std::time(&end);
+ std::cout << "Finished conversion. New AlfTree has weight: "<<tree->GetWeight()<<". Takes "<<(end-start)<<" s." << std::endl;
+ return *tree;
+}
+
+// -------------------------------------------------------------------------
+AirwaysTree& CreateAirwaysTreeFromSegmentation(Vec3 seed, TImagePointer input_image)
+{
+ std::string err = myWorkspace.Execute( "eb" );
+ if( err != "" )
+ {
+ std::cerr << "Error: " << err << std::endl;
+ std::exit( 1 );
+
+ } // fi
+ return( ConvertFilterToAirwaysTree( ) );
+}
+
+// -------------------------------------------------------------------------
+vector<TreeInfo> ReadInputFile(const char* filename)
+{
+ // Variables
+ std::ifstream infile(filename);
+ std::string line;
+ std::string folderpath_allResults;
+ vector<TreeInfo> vectorInfo;
+ // Parse the file
+ bool firstLine = false;
+
+ while (std::getline(infile, line))
+ {
+ // First line contains the output folder
+ std::istringstream iss(line);
+ if (!firstLine)
+ {
+ if (!(iss >> folderpath_allResults))
+ {
+ std::cout << "no file" << std::endl;
+ return vectorInfo;
+ } //if
+ firstLine = true;
+ } //if
+ // Other lines, not the first one, contain the information for the airways to be created
+ else
+ {
+ TreeInfo treeInfo;
+ float point_trachea[3];
+ std::string filepath_airwayImage, name_pig, name_image;
+ if (!(iss >> point_trachea[0] >> point_trachea[1] >> point_trachea[2] >> filepath_airwayImage >> name_pig >> name_image))
+ {
+ std::cout << "There is no tree information in the file." << std::endl;
+ break;
+ } //if
+ else
+ {
+ // Print information
+ std::cout << "Point trachea:[" << point_trachea[0] << "," << point_trachea[1] << "," << point_trachea[2] << "]" << std::endl;
+ }
+
+ //Save info.
+
+ //Save seed info.
+ Vec3 seed(point_trachea[0], point_trachea[1], point_trachea[2]); //real coords seed
+ treeInfo.seed = seed;
+ // Create the outputs
+ treeInfo.folderpath_pigResults = folderpath_allResults + "/" + name_pig + "/";
+ treeInfo.pig_name = name_pig;
+ treeInfo.image_name = name_image;
+ // Read the image
+ /*
+ treeInfo.image = ReadImage<TInputImage>(filepath_airwayImage);
+ */
+ Load_cpPlugins( "./plugins.cfg", "./workspace_airwaysappli.wxml" );
+
+ // Execute first pipeline's part
+ std::stringstream seed_stream;
+ seed_stream
+ << point_trachea[0] << " "
+ << point_trachea[1] << " "
+ << point_trachea[2];
+ myWorkspace.SetParameter( "FileNames@reader", filepath_airwayImage );
+ myWorkspace.SetParameter( "Text@seed", seed_stream.str( ) );
+ vectorInfo.push_back(treeInfo);
+
+ } //else
+ } //while
+ return vectorInfo;
+}
+
+// -------------------------------------------------------------------------
+void printCommonTreeBetweenTwoTrees(AirwaysTree tree_A, AirwaysTree tree_B, std::vector< std::pair< std::pair<Node*, Node*>, std::pair<Node*, Node*> > > vector_pair_edges_A_to_B, unsigned int Q, unsigned int F)
+{
+ std::cout << "printCommonTreeBetweenTwoTrees, edges:" << vector_pair_edges_A_to_B.size() << std::endl;
+
+ // Vtk points, cell array, and colors
+ vtkSmartPointer<vtkPoints> points_common = vtkSmartPointer<vtkPoints>::New();
+ vtkSmartPointer<vtkCellArray> lines_common = vtkSmartPointer<vtkCellArray>::New();
+ vtkSmartPointer<vtkUnsignedCharArray> colors_common = vtkSmartPointer<vtkUnsignedCharArray>::New();
+ colors_common->SetNumberOfComponents(3);
+ colors_common->SetName("Colors");
+
+ // Add all the points
+ // 0 - 1000 points for tree A and from 1001 for tree B
+ points_common->SetNumberOfPoints(1000 + tree_B.GetWeight() + 100);
+
+ // Add points for tree A
+ vec_nodes vector_nodesA = tree_A.GetNodes();
+ for(vec_nodes::iterator it_nodesA = vector_nodesA.begin(); it_nodesA != vector_nodesA.end(); ++it_nodesA)
+ {
+ points_common->SetPoint((*it_nodesA)->GetId(),(*it_nodesA)->GetCoords().GetVec3());
+ }
+
+ // Add points for tree B
+ vec_nodes vector_nodesB = tree_B.GetNodes();
+ for(vec_nodes::iterator it_nodesB = vector_nodesB.begin(); it_nodesB != vector_nodesB.end(); ++it_nodesB)
+ {
+ points_common->SetPoint((*it_nodesB)->GetId()+1000,(*it_nodesB)->GetCoords().GetVec3());
+ }
+
+ // Add the edges for both trees
+ std::vector< std::pair< std::pair<Node*, Node*>, std::pair<Node*, Node*> > >::iterator it_edges = vector_pair_edges_A_to_B.begin();
+ int number_pairs = 0;
+ for(; it_edges != vector_pair_edges_A_to_B.end(); ++it_edges)
+ {
+ std::cout << "Pair:" << number_pairs << std::endl;
+ std::pair<Node*, Node*> edge_A = (*it_edges).first;
+ std::pair<Node*, Node*> edge_B = (*it_edges).second;
+
+ vec_nodes path_A;
+ edge_A.first->GetPathToNode(edge_A.second, path_A);
+
+ vec_nodes path_B;
+ edge_B.first->GetPathToNode(edge_B.second, path_B);
+
+ // Set color to be used
+ int numColor = number_pairs % 6;
+ unsigned char color[3];
+ color[0] = green[0];
+ color[1] = green[1];
+ color[2] = green[2];
+ if(numColor == 1)
+ {
+ color[0] = red[0];
+ color[1] = red[1];
+ color[2] = red[2];
+ }
+ else if(numColor == 2)
+ {
+ color[0] = blue[0];
+ color[1] = blue[1];
+ color[2] = blue[2];
+ }
+ else if(numColor == 3)
+ {
+ color[0] = yellow[0];
+ color[1] = yellow[1];
+ color[2] = yellow[2];
+ }
+ else if(numColor == 4)
+ {
+ color[0] = purple[0];
+ color[1] = purple[1];
+ color[2] = purple[2];
+ }
+ else if(numColor == 5)
+ {
+ color[0] = cyan[0];
+ color[1] = cyan[1];
+ color[2] = cyan[2];
+ }
+
+ number_pairs++;
+
+ // Trace line A
+ //if(path_A.size() > 0 && number_pairs < 50)
+ if( path_A.size() )
+ {
+ vtkSmartPointer<vtkLine> line = vtkSmartPointer<vtkLine>::New();
+ lines_common->InsertNextCell(path_A.size());
+ int id_point = 0;
+ for(vec_nodes::iterator it_pathA = path_A.begin(); it_pathA != path_A.end(); ++it_pathA)
+ {
+ lines_common->InsertCellPoint((*it_pathA)->GetId());
+ //line->GetPointIds()->SetId( id_point, (*it_pathA)->GetId() );
+ //id_point++;
+ }
+ //lines_common->InsertNextCell(line);
+ colors_common->InsertNextTupleValue(color);
+ }
+ else
+ {
+ std::cout << "No path A" << std::endl;
+ }
+
+
+ // Trace line B
+ //if(path_B.size() > 0 && number_pairs < 50)
+ if( path_B.size() )
+ {
+ vtkSmartPointer<vtkLine> line = vtkSmartPointer<vtkLine>::New();
+ lines_common->InsertNextCell(path_B.size());
+ int id_point = 0;
+ for(vec_nodes::iterator it_pathB = path_B.begin(); it_pathB != path_B.end(); ++it_pathB)
+ {
+ lines_common->InsertCellPoint((*it_pathB)->GetId()+1000);
+ //line->GetPointIds()->SetId( id_point, (*it_pathB)->GetId()+1000 );
+ id_point++;
+ }
+ //lines_common->InsertNextCell(line);
+ colors_common->InsertNextTupleValue(color);
+ }
+ else
+ {
+ std::cout << "No path B" << std::endl;
+ }
+ }
+
+ // Save the polydata
+ // Create the polydata
+ vtkSmartPointer<vtkPolyData> linesPolyData = vtkSmartPointer<vtkPolyData>::New();
+
+ // Set the points, lines, and colors
+ linesPolyData->SetPoints(points_common);
+ linesPolyData->SetLines(lines_common);
+ linesPolyData->GetCellData()->SetScalars(colors_common);
+
+ // ------------------------------------------------
+ // Write the vtk file
+
+ // Create the pathfile to save
+ std::stringstream filepath_actualIteration;
+ filepath_actualIteration << tree_A.GetResultPath() << "/" << tree_A.GetImageName() << "_" << tree_B.GetImageName() << "_Q_" << Q << "_F_" << F << ".vtk";
+ std::cout << "File to save:" << filepath_actualIteration.str() << std::endl;
+
+ // Create the writer
+ vtkSmartPointer<vtkPolyDataWriter> writer = vtkSmartPointer<vtkPolyDataWriter>::New();
+ writer->SetFileName(filepath_actualIteration.str().c_str());
+
+ // Set input and write
+#if VTK_MAJOR_VERSION <= 5
+ writer->SetInput(linesPolyData);
+#else
+ writer->SetInputData(linesPolyData);
+#endif
+ writer->Write();
+
+
+
+ // ***************************************
+ // Save the links in a file
+ // *******
+
+ // Create the outputfile
+ std::stringstream filepath_evaluation;
+ filepath_evaluation << tree_A.GetResultPath() << "/" << tree_A.GetImageName() << "_" << tree_B.GetImageName() << "_Q_" << Q << "_F_" << F << ".csv";
+
+ std::ofstream file(filepath_evaluation.str().c_str());
+ if(file.is_open())
+ {
+ // Save the header
+ //file << "Node_" << tree_A.GetImageName() << " " << "Node_" << tree_B.GetImageName()<< "\n";
+ file << "TreeName "
+ "idLocalN1 idLocalN2 "
+ "Node1x Node1y Node1z "
+ "Node2x Node2y Node2z "
+ "idGlobalN1 "
+ "idMatch" << std::endl;
+
+
+ std::vector< std::pair< std::pair<Node*, Node*>, std::pair<Node*, Node*> > >::iterator it_edges = vector_pair_edges_A_to_B.begin();
+ for(; it_edges != vector_pair_edges_A_to_B.end(); ++it_edges)
+ {
+
+ std::pair<Node*, Node*> edge_A = (*it_edges).first;
+ std::pair<Node*, Node*> edge_B = (*it_edges).second;
+
+ //file << edge_A.second->GetId() << " " << edge_B.second->GetId()+1000 << "\n";
+ Vec3 coord_EndNodeA = edge_A.second->GetCoords();
+ Vec3 coord_EndNodeB = edge_B.second->GetCoords();
+ file << tree_A.GetImageName() << " " <<
+ edge_A.second->GetId() << " " << edge_B.second->GetId()+1000 << " " <<
+ coord_EndNodeA[0] << " " << coord_EndNodeA[1] << " " << coord_EndNodeA[2] << " " <<
+ coord_EndNodeB[0] << " " << coord_EndNodeB[1] << " " << coord_EndNodeB[2] << " " <<
+ tree_A.GetImageName() << coord_EndNodeA[0] << coord_EndNodeA[1] << coord_EndNodeA[2] << " " <<
+ tree_B.GetImageName() << coord_EndNodeB[0] << coord_EndNodeB[1] << coord_EndNodeB[2] <<"\n";
+
+ file << tree_B.GetImageName() << " " <<
+ edge_B.second->GetId()+1000 << " " << edge_A.second->GetId() << " " <<
+ coord_EndNodeB[0] << " " << coord_EndNodeB[1] << " " << coord_EndNodeB[2] << " " <<
+ coord_EndNodeA[0] << " " << coord_EndNodeA[1] << " " << coord_EndNodeA[2] << " " <<
+ tree_B.GetImageName() << coord_EndNodeB[0] << coord_EndNodeB[1] << coord_EndNodeB[2] << " " <<
+ tree_A.GetImageName() << coord_EndNodeA[0] << coord_EndNodeA[1] << coord_EndNodeA[2] << "\n";
+ }
+ }
+
+ file.close();
+ // *******
+ // ***************************************
+
+
+
+
+ std::cout << "printCommonTreeBetweenTwoTrees ... OK" << std::endl;
+}
+
+// -------------------------------------------------------------------------
+void printMatchingResultToFile(AirwaysTree tree_A, AirwaysTree tree_B, std::map< unsigned int, std::vector<Node*> > map_A_to_B, std::map< unsigned int, std::vector<Node*> > map_B_to_A, unsigned int Q, unsigned int F)
+{
+ std::cout << "Printing matching result ... " << std::endl;
+
+ // Variables and types
+ typedef std::map< unsigned int, vec_nodes > map_id_node;
+
+ // Create the outputfile
+ std::stringstream filepath_evaluation;
+ filepath_evaluation << tree_A.GetResultPath() << "/" << tree_A.GetImageName() << "_" << tree_B.GetImageName() << "_Q_" << Q << "_F_" << F << "_matching.csv";
+
+ std::ofstream file(filepath_evaluation.str().c_str());
+ if(file.is_open())
+ {
+ // Save the header
+ file << "TreeName idLocalN1 idLocalN2 "
+ "Node1x Node1y Node1z "
+ "Node2x Node2y Node2z "
+ "idGlobalN1 "
+ "idMatch "
+ "typeMatch_match_1_nonmatch_0 "
+ "depth "
+ "leaf"<< std::endl;
+
+ // ---------------
+ // From A to B
+
+ // Save the match or non-match for each node from A to B
+ for(int id_a=1; id_a <= tree_A.GetWeight( ); ++id_a)
+ {
+ // Get the actual node in tree A
+ Node* node_a = tree_A.GetNodeById( id_a );
+ Vec3 coords_a = node_a->GetCoords();
+ unsigned int depth_a = tree_A.GetDepthById(id_a);
+ bool leaf_a = node_a->IsLeaf();
+
+ // Check if the node was matched
+ map_id_node::iterator it_a2b = map_A_to_B.find( id_a );
+ if( it_a2b != map_A_to_B.end( ) )
+ {
+ // Get the correspoding matching nodes and print them
+ vec_nodes nodes_B = (*it_a2b).second;
+
+ vec_nodes::iterator it_nodes_b = nodes_B.begin( );
+ for( ; it_nodes_b != nodes_B.end( ); ++it_nodes_b)
+ {
+ Vec3 coords_b = (*it_nodes_b)->GetCoords();
+
+ file << tree_A.GetImageName() << " " << id_a << " " << (*it_nodes_b)->GetId()+1000 << " "
+ << coords_a[0] << " " << coords_a[1] << " " << coords_a[2] << " "
+ << coords_b[0] << " " << coords_b[1] << " " << coords_b[2] << " "
+ << tree_A.GetImageName() << coords_a[0] << coords_a[1] << coords_a[2] << " "
+ << tree_A.GetImageName() << coords_a[0] << coords_a[1] << coords_a[2] << coords_b[0] << coords_b[1] << coords_b[2] << " "
+ << "1" << " "
+ << depth_a << " "
+ << leaf_a << "\n";
+ }
+ }
+ else
+ {
+ file << tree_A.GetImageName() << " " << id_a << " " << "0" << " "
+ << coords_a[0] << " " << coords_a[1] << " " << coords_a[2] << " "
+ << "0" << " " << "0" << " " << "0" << " "
+ << tree_A.GetImageName() << coords_a[0] << coords_a[1] << coords_a[2] << " "
+ << tree_A.GetImageName() << coords_a[0] << coords_a[1] << coords_a[2] << "0" << "0" << "0" << " "
+ << "0" << " "
+ << depth_a << " "
+ << leaf_a << "\n";
+ }
+ }
+
+ // ---------------
+ // From B to A
+
+ // Save the match or non-match for each node from A to B
+ for(int id_b=1; id_b <= tree_B.GetWeight( ); ++id_b)
+ {
+ // Get the actual node in tree B
+ Node* node_b = tree_B.GetNodeById( id_b );
+ Vec3 coords_b = node_b->GetCoords();
+ unsigned int depth_b = tree_B.GetDepthById(id_b);
+ bool leaf_b = node_b->IsLeaf();
+
+ // Check if the node was matched
+ map_id_node::iterator it_b2a = map_B_to_A.find( id_b );
+ if( it_b2a != map_B_to_A.end( ) )
+ {
+ // Get the correspoding matching nodes and print them
+ vec_nodes nodes_A = (*it_b2a).second;
+
+ vec_nodes::iterator it_nodes_a = nodes_A.begin( );
+ for( ; it_nodes_a != nodes_A.end( ); ++it_nodes_a)
+ {
+ Vec3 coords_a = (*it_nodes_a)->GetCoords();
+
+ file << tree_B.GetImageName() << " " << id_b+1000 << " " << (*it_nodes_a)->GetId() << " "
+ << coords_b[0] << " " << coords_b[1] << " " << coords_b[2] << " "
+ << coords_a[0] << " " << coords_a[1] << " " << coords_a[2] << " "
+ << tree_B.GetImageName() << coords_b[0] << coords_b[1] << coords_b[2] << " "
+ << tree_B.GetImageName() << coords_b[0] << coords_b[1] << coords_b[2] << coords_a[0] << coords_a[1] << coords_a[2] << " "
+ << "1" << " "
+ << depth_b << " "
+ << leaf_b << "\n";
+ }
+ }
+ else
+ {
+ file << tree_B.GetImageName() << " " << id_b+1000 << " " << "0" << " "
+ << coords_b[0] << " " << coords_b[1] << " " << coords_b[2] << " "
+ << "0" << " " << "0" << " " << "0" << " "
+ << tree_B.GetImageName() << coords_b[0] << coords_b[1] << coords_b[2] << " "
+ << tree_B.GetImageName() << coords_b[0] << coords_b[1] << coords_b[2] << "0" << "0" << "0" << " "
+ << "0" << " "
+ << depth_b << " "
+ << leaf_b << "\n";
+ }//esle
+ }//rof
+ }//fi
+
+ file.close(); // Close the file
+
+ std::cout << "Printing matching result DONE" << std::endl;
+}
+
+// -------------------------------------------------------------------------
+void createLinesAndPointsForVTK(const Node* node, vtkSmartPointer<vtkPoints>& pts,vtkSmartPointer<vtkCellArray>& lines, vtkSmartPointer<vtkUnsignedCharArray>& colors, bool common, bool isRoot)
+{
+ // Insert the actual point/node
+ //vtkIdType id_father = idNonRoot;
+ //if(isRoot)
+ //id_father = pts->InsertNextPoint(node->GetCoords().GetVec3());
+ vtkIdType id_father = node->GetId();
+ if(isRoot)
+ pts->SetPoint(id_father,node->GetCoords().GetVec3());
+
+ // Iterate over the children
+ const vec_nodes children = node->GetChildren();
+ for (vec_nodes::const_iterator it_child = children.begin(); it_child != children.end(); ++it_child)
+ {
+ if (!(*it_child)->IsMarked() && common)
+ continue;
+
+ //vtkIdType id_son = pts->InsertNextPoint((*it_child)->GetCoords().GetVec3());
+ vtkIdType id_son = (*it_child)->GetId();
+ pts->SetPoint(id_son,(*it_child)->GetCoords().GetVec3());
+
+ // Set color to be used
+ int numColor = (*it_child)->GetLevel() % 4;
+ if(numColor == 0)
+ colors->InsertNextTupleValue(green);
+ else if(numColor == 1)
+ colors->InsertNextTupleValue(red);
+ else if(numColor == 2)
+ colors->InsertNextTupleValue(red);
+ else
+ colors->InsertNextTupleValue(red);
+
+ // Add the line
+ vtkSmartPointer<vtkLine> line = vtkSmartPointer<vtkLine>::New();
+ line->GetPointIds()->SetId(0, id_father);
+ line->GetPointIds()->SetId(1, id_son);
+ lines->InsertNextCell(line);
+
+ createLinesAndPointsForVTK(*it_child, pts, lines, colors, common, false);
+ } //for
+}
+
+// eof - $RCSfile$
git://git.creatis.insa-lyon.fr / FrontAlgorithms.git / commitdiff