From 36d560deb506320496c89fdea506cb916ae3dbca Mon Sep 17 00:00:00 2001 From: Leonardo Florez-Valencia Date: Wed, 13 Apr 2016 19:06:52 -0500 Subject: [PATCH] ... --- CMakeLists.txt | 3 +- appli/CMakeLists.txt | 14 +- .../AirwaysLib/CMakeLists.txt | 45 + .../AirwaysLib/airwaysEdge.cxx | 739 ++++++ .../TempAirwaysAppli/AirwaysLib/airwaysEdge.h | 323 +++ .../AirwaysLib/airwaysNode.cxx | 1655 ++++++++++++ .../TempAirwaysAppli/AirwaysLib/airwaysNode.h | 597 +++++ .../AirwaysLib/airwaysTree.cxx | 2251 +++++++++++++++++ .../TempAirwaysAppli/AirwaysLib/airwaysTree.h | 372 +++ .../AirwaysLib/airwaysTreeTypeDefinition.h | 132 + appli/TempAirwaysAppli/CMakeLists.txt | 22 + appli/TempAirwaysAppli/MathLib/CMakeLists.txt | 38 + appli/TempAirwaysAppli/MathLib/Quaternion.cpp | 168 ++ appli/TempAirwaysAppli/MathLib/Quaternion.h | 174 ++ appli/TempAirwaysAppli/MathLib/vec3.cxx | 233 ++ appli/TempAirwaysAppli/MathLib/vec3.h | 221 ++ appli/TempAirwaysAppli/TempAirwaysAppli.cxx | 1025 ++++++++ 17 files changed, 8007 insertions(+), 5 deletions(-) create mode 100644 appli/TempAirwaysAppli/AirwaysLib/CMakeLists.txt create mode 100644 appli/TempAirwaysAppli/AirwaysLib/airwaysEdge.cxx create mode 100644 appli/TempAirwaysAppli/AirwaysLib/airwaysEdge.h create mode 100644 appli/TempAirwaysAppli/AirwaysLib/airwaysNode.cxx create mode 100644 appli/TempAirwaysAppli/AirwaysLib/airwaysNode.h create mode 100644 appli/TempAirwaysAppli/AirwaysLib/airwaysTree.cxx create mode 100644 appli/TempAirwaysAppli/AirwaysLib/airwaysTree.h create mode 100644 appli/TempAirwaysAppli/AirwaysLib/airwaysTreeTypeDefinition.h create mode 100644 appli/TempAirwaysAppli/CMakeLists.txt create mode 100644 appli/TempAirwaysAppli/MathLib/CMakeLists.txt create mode 100644 appli/TempAirwaysAppli/MathLib/Quaternion.cpp create mode 100644 appli/TempAirwaysAppli/MathLib/Quaternion.h create mode 100644 appli/TempAirwaysAppli/MathLib/vec3.cxx create mode 100644 appli/TempAirwaysAppli/MathLib/vec3.h create mode 100644 appli/TempAirwaysAppli/TempAirwaysAppli.cxx diff --git a/CMakeLists.txt b/CMakeLists.txt index c022af8..30f36c0 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -23,6 +23,7 @@ FIND_PACKAGE(cpPlugins REQUIRED) ## ============= OPTION(BUILD_PLUGINS "Build plugins" ON) +OPTION(BUILD_TempAirwaysAppli "Build temporary appli" OFF) ## ========================= ## == Include directories == @@ -47,7 +48,7 @@ SUBDIRS( cmake lib plugins - # appli + appli ) ## eof - $RCSfile$ diff --git a/appli/CMakeLists.txt b/appli/CMakeLists.txt index 6d67748..d9ec026 100644 --- a/appli/CMakeLists.txt +++ b/appli/CMakeLists.txt @@ -1,8 +1,14 @@ -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$ diff --git a/appli/TempAirwaysAppli/AirwaysLib/CMakeLists.txt b/appli/TempAirwaysAppli/AirwaysLib/CMakeLists.txt new file mode 100644 index 0000000..2c05400 --- /dev/null +++ b/appli/TempAirwaysAppli/AirwaysLib/CMakeLists.txt @@ -0,0 +1,45 @@ +## ============================= +## = 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$ diff --git a/appli/TempAirwaysAppli/AirwaysLib/airwaysEdge.cxx b/appli/TempAirwaysAppli/AirwaysLib/airwaysEdge.cxx new file mode 100644 index 0000000..56a9084 --- /dev/null +++ b/appli/TempAirwaysAppli/AirwaysLib/airwaysEdge.cxx @@ -0,0 +1,739 @@ +/* + * airwaysEdge.cxx + * + * Created on: May 12, 2014 + * Author: caceres + */ + +#include "airwaysEdge.h" +#include + +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::max(); + double max = std::numeric_limits::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 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::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 */ diff --git a/appli/TempAirwaysAppli/AirwaysLib/airwaysEdge.h b/appli/TempAirwaysAppli/AirwaysLib/airwaysEdge.h new file mode 100644 index 0000000..fe7adf1 --- /dev/null +++ b/appli/TempAirwaysAppli/AirwaysLib/airwaysEdge.h @@ -0,0 +1,323 @@ +/* + * airwaysEdge.h + * + * Created on: May 12, 2014 + * Author: caceres + */ + +#ifndef AIRWAYSEDGE_H_ +#define AIRWAYSEDGE_H_ + +#include +#include +#include +#include +#include + +#include "airwaysNode.h" +#include "Quaternion.h" +#include + +/** + * @brief Airways project namespace + */ +namespace airways +{ + +/* + * Pair of 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 pair_posVox_rad; + +/* + * A vector stocking all the points between two nodes + */ +typedef std::vector 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 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 + * + */ + 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_ */ diff --git a/appli/TempAirwaysAppli/AirwaysLib/airwaysNode.cxx b/appli/TempAirwaysAppli/AirwaysLib/airwaysNode.cxx new file mode 100644 index 0000000..27f6845 --- /dev/null +++ b/appli/TempAirwaysAppli/AirwaysLib/airwaysNode.cxx @@ -0,0 +1,1655 @@ +/* + * airwaysNode.cxx + * + * Created on: May 12, 2014 + * Author: Diego Cáceres + * + * Modified by: Alfredo Morales Pinzón + */ + +#include + +#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::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 map_dist_pairNodes_A_to_B = this->CalculateDistanceToAllBranches_FatherAndFamily(nodeB, Q, F); + std::multimap 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_A_to_B = this->GetClosestBranch_FatherAndFamily(nodeB); + // 2.1.2 Get best translated branch from B to A + //std::pair pair_B_to_A = nodeB->GetClosestBranch_FatherAndFamily(this); + + // ------------------------ + // ------------------------ + // Making the matching faster + // Get the next best pair + std::pair pair_A_to_B = GetPairWithClosestDistance_FirstNodeNonMarked(map_dist_pairNodes_A_to_B); + std::pair 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 map_dist_pairNodes_A_to_B = this->CalculateDistanceToAllBranches_FatherAndFamily(nodeB); + //std::map 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, 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 best_pair_score(pair_best, distance_bestFinal); + return best_pair_score; +} + +std::multimap< double, std::pair > 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 actual_dist_pair(distance_actual,pair_actual); + map_dist_pairNodeNode.insert(actual_dist_pair); + } + } + + return map_dist_pairNodeNode; +} + + +/*std::pair< std::pair, 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 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::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& 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 <= 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 <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 > Node::GetPairWithClosestDistance_FirstNodeNonMarked(std::multimap map_dist_pairNodes) +{ + // Variables + bool pairFound = false; + std::pair pairNodes_null (NULL, NULL); + std::pair > pair_closest_dist_pairNodes (-1, pairNodes_null); + + std::multimap::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 */ diff --git a/appli/TempAirwaysAppli/AirwaysLib/airwaysNode.h b/appli/TempAirwaysAppli/AirwaysLib/airwaysNode.h new file mode 100644 index 0000000..f002038 --- /dev/null +++ b/appli/TempAirwaysAppli/AirwaysLib/airwaysNode.h @@ -0,0 +1,597 @@ +/* + * airwaysNode.h + * + * Created on: May 12, 2014 + * Author: caceres + * + * Modified by: Alfredo Morales Pinzón + */ + +#ifndef AIRWAYSNODE_H_ +#define AIRWAYSNODE_H_ + +#include +#include + +#include "../MathLib/vec3.h" +#include "airwaysEdge.h" +#include "Quaternion.h" +#include + +// 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 pair_nodes; + + // Vector of nodes + typedef std::vector 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& 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, 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 > 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, 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& 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 + */ + std::pair > GetPairWithClosestDistance_FirstNodeNonMarked(std::multimap 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_ */ diff --git a/appli/TempAirwaysAppli/AirwaysLib/airwaysTree.cxx b/appli/TempAirwaysAppli/AirwaysLib/airwaysTree.cxx new file mode 100644 index 0000000..0f0f821 --- /dev/null +++ b/appli/TempAirwaysAppli/AirwaysLib/airwaysTree.cxx @@ -0,0 +1,2251 @@ +/* + * 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 >& map_A_to_B, std::map< unsigned int, std::vector >& map_B_to_A, std::vector< std::pair< std::pair, std::pair > >& 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 nonCommonNodes_A; + std::map commonNodes_A; + std::map nonCommonNodes_B; + std::map 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::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::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::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::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, 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 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 <m_root->printUpToLevel(level); + cout << "Printing done" <m_root->printNodeAndChildrenIds( ); + cout << "Printing nodes and children DONE" <= 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" <m_root->createQuaternionsUpToLevel(level, level); + } + else + cout << "Quaternions must be created from level 2." <m_root->getStasticsBifurcations(p); + + cout << "Bifurcations:" << endl; + for(int i = 0; i < 6; i++) + { + cout << i+1 << "->" << p[i] << endl; + } + + cout << "Statistics bifurcations done " <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 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 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 colors = vtkSmartPointer::New(); + colors->SetNumberOfComponents(3); + colors->SetName("Colors"); + + //pointer of lines and points + vtkSmartPointer lines = vtkSmartPointer::New(); + vtkSmartPointer pts = vtkSmartPointer::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 linesPolyData = vtkSmartPointer::New(); + linesPolyData->SetPoints(pts); + linesPolyData->SetLines(lines); + + linesPolyData->GetCellData()->SetScalars(colors); + + // Write the file + vtkSmartPointer writer = vtkSmartPointer::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& pts, + vtkSmartPointer& lines, + vtkSmartPointer& 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 line = vtkSmartPointer::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 + diff --git a/appli/TempAirwaysAppli/AirwaysLib/airwaysTree.h b/appli/TempAirwaysAppli/AirwaysLib/airwaysTree.h new file mode 100644 index 0000000..0134708 --- /dev/null +++ b/appli/TempAirwaysAppli/AirwaysLib/airwaysTree.h @@ -0,0 +1,372 @@ +/* + * 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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// ITK +#include "itkImage.h" +#include "itkImageRegionIterator.h" +#include "itkLineIterator.h" +#include +#include + +namespace airways +{ + +// Types + +// Vector of nodes +typedef std::vector vec_nodes; + +// Pair of nodes +typedef std::pair pair_nodes; + +//Pair of a pair of nodes and a score +typedef std::pair Pair_PairNodes_Score; + +/** + * + */ +//typedef std::pair Pair_PairNodes_Score; + +/** + * Vector of pairs of a pair of nodes and a score + */ +typedef std::vector Vector_Pair_PairNodes_Score; + +/*! + * A vector stocking all the points between two nodes + */ +typedef std::vector 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 >& map_A_to_B, std::map< unsigned int, std::vector >& map_B_to_A, std::vector< std::pair< std::pair, std::pair > >& 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& pts, + vtkSmartPointer& lines, + vtkSmartPointer& 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_ */ diff --git a/appli/TempAirwaysAppli/AirwaysLib/airwaysTreeTypeDefinition.h b/appli/TempAirwaysAppli/AirwaysLib/airwaysTreeTypeDefinition.h new file mode 100644 index 0000000..235430a --- /dev/null +++ b/appli/TempAirwaysAppli/AirwaysLib/airwaysTreeTypeDefinition.h @@ -0,0 +1,132 @@ +/* + * airwaysTreeTypeDefinition.h + * + * Created on: May 19, 2014 + * Author: caceres + */ + +#ifndef AIRWAYSTREETYPEDEFINITION_H_ +#define AIRWAYSTREETYPEDEFINITION_H_ + +//std +#include +#include +#include +#include +//boost +#include +//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 pair_posVox_rad; + typedef std::vector vec_pair_posVox_rad; + + /** + * Pair of id nodes + */ + + + /** + * Pair of a pair of nodes and a score + */ + //typedef std::pair Pair_PairNodes_Score; + + //---------------------------------------------------------------------------- + /* + * Typedef declaration + */ + /** + * @typedef Defines a graph using the Boost DiGraph library + */ + typedef boost::adjacency_list DiGraph; + //end typedef Graph + //-------------------------------------------------------------------------- + /** + * @typedef Defines boost::vertex_descriptor of DiGraph + */ + typedef boost::graph_traits::vertex_descriptor DiGraph_VertexDescriptor; + //-------------------------------------------------------------------------- + ///DiGraph edge_descriptor + typedef boost::graph_traits::edge_descriptor DiGraph_EdgeDescriptor; + //-------------------------------------------------------------------------- + ///DiGraph edge pair + typedef std::pair DiGraph_EdgePair; + //-------------------------------------------------------------------------- + ///DiGraph Vertex IndexMap + typedef boost::property_map::type DiGraph_VertexIndexMap; + //-------------------------------------------------------------------------- + ///DiGraph VertexIterator + typedef boost::graph_traits::vertex_iterator DiGraph_VertexIterator; + //-------------------------------------------------------------------------- + ///DiGraph Vertex iterator pair + typedef std::pair DiGraph_VertexIteratorPair; + //-------------------------------------------------------------------------- + ///DiGraph Edge iterator + typedef boost::graph_traits::edge_iterator DiGraph_EdgeIterator; + //-------------------------------------------------------------------------- + ///DiGraph In Edge iterator + typedef boost::graph_traits::in_edge_iterator DiGraph_InEdgeIterator; + //-------------------------------------------------------------------------- + ///DiGraph Out Edge iterator + typedef boost::graph_traits::out_edge_iterator DiGraph_OutEdgeIterator; + //-------------------------------------------------------------------------- + ///DiGraph Adjacency iterator + typedef boost::graph_traits::adjacency_iterator DiGraph_AdjacencyIterator; + //-------------------------------------------------------------------------- + ///DiGraph Vertex descriptor vector + typedef std::vector DiGraph_VertexVector; + //-------------------------------------------------------------------------- + ///DiGraph Vertex descriptor queue + typedef std::queue DiGraph_VertexQueue; + //-------------------------------------------------------------------------- + ///DiGraph Edge descriptor vector + typedef std::vector DiGraph_EdgeVector; + //---------------------------------------------------------------------------- + const unsigned int Dim = 3; + typedef double TPixel; + typedef unsigned char TPixel2; + typedef itk::Image TInputImage; + typedef itk::Image TInputImage2; + typedef itk::Image TInputImage3; + //---------------------------------------------------------------------------- + typedef itk::ImageDuplicator DuplicatorType; + //-------------------------------------------------------------------------- + typedef itk::ConstNeighborhoodIterator ConstNeighborhoodIterator; + typedef itk::ImageRegionIterator RegionIterator; + //---------------------------------------------------------------------------- + typedef itk::SphereSpatialFunction<> SphereType; + //---------------------------------------------------------------------------- + typedef itk::ConnectedThresholdImageFilter ConnectedFilterType; + //---------------------------------------------------------------------------- + typedef itk::SubtractImageFilter SubtractImageFilterType; + //---------------------------------------------------------------------------- + typedef itk::CastImageFilter CastFilterType; + typedef itk::ConnectedComponentImageFilter ConnectedComponentImageFilterType; + //---------------------------------------------------------------------------- + + //---------------------------------------------------------------------------- + typedef std::vector EdgeVector; + //---------------------------------------------------------------------------- + typedef std::queue Vec3Queue; + typedef std::vector Vec3List; +// +/* + * End of typedef declaration + */ +} + +#endif /* AIRWAYSTREETYPEDEFINITION_H_ */ diff --git a/appli/TempAirwaysAppli/CMakeLists.txt b/appli/TempAirwaysAppli/CMakeLists.txt new file mode 100644 index 0000000..4b80e61 --- /dev/null +++ b/appli/TempAirwaysAppli/CMakeLists.txt @@ -0,0 +1,22 @@ + +## 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$ diff --git a/appli/TempAirwaysAppli/MathLib/CMakeLists.txt b/appli/TempAirwaysAppli/MathLib/CMakeLists.txt new file mode 100644 index 0000000..6812792 --- /dev/null +++ b/appli/TempAirwaysAppli/MathLib/CMakeLists.txt @@ -0,0 +1,38 @@ +## ============================= +## = 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$ diff --git a/appli/TempAirwaysAppli/MathLib/Quaternion.cpp b/appli/TempAirwaysAppli/MathLib/Quaternion.cpp new file mode 100644 index 0000000..6beaa48 --- /dev/null +++ b/appli/TempAirwaysAppli/MathLib/Quaternion.cpp @@ -0,0 +1,168 @@ +//---------------------------------------------------------------------------- +// 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 << " " < + +//----------------- +// OWN +//----------------- +#include +#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 diff --git a/appli/TempAirwaysAppli/MathLib/vec3.cxx b/appli/TempAirwaysAppli/MathLib/vec3.cxx new file mode 100644 index 0000000..ca35977 --- /dev/null +++ b/appli/TempAirwaysAppli/MathLib/vec3.cxx @@ -0,0 +1,233 @@ +#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::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; +} + diff --git a/appli/TempAirwaysAppli/MathLib/vec3.h b/appli/TempAirwaysAppli/MathLib/vec3.h new file mode 100644 index 0000000..22f9ee3 --- /dev/null +++ b/appli/TempAirwaysAppli/MathLib/vec3.h @@ -0,0 +1,221 @@ +/** + * 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 +#include +#include +#include // 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 diff --git a/appli/TempAirwaysAppli/TempAirwaysAppli.cxx b/appli/TempAirwaysAppli/TempAirwaysAppli.cxx new file mode 100644 index 0000000..3dfbcc3 --- /dev/null +++ b/appli/TempAirwaysAppli/TempAirwaysAppli.cxx @@ -0,0 +1,1025 @@ +#include +#include +#include + +#include +#include + +#include "vec3.h" +#include "airwaysTree.h" + +#include +#include + +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 ReadInputFile(const char* filename); +void printCommonTreeBetweenTwoTrees(AirwaysTree tree_A, AirwaysTree tree_B, std::vector< std::pair< std::pair, std::pair > > 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 > map_A_to_B, std::map< unsigned int, std::vector > map_B_to_A, unsigned int Q, unsigned int F); +void createLinesAndPointsForVTK(const Node* node, vtkSmartPointer& pts,vtkSmartPointer& lines, vtkSmartPointer& 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(), "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(),"Get a subtree(s) by levels - result: a vtk and reconstructed img .mhd") + ("subtree_length", po::value(),"Get a subtree(s) by length - result: a vtk and reconstructed img .mhd") + ("subtree_diameter", po::value(),"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 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(); + 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 > map_A_to_B; + std::map< unsigned int, std::vector > map_B_to_A; + std::vector< std::pair< std::pair, std::pair > > 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 >::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::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 >::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::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 pts = vtkSmartPointer::New(); + vtkSmartPointer lines = vtkSmartPointer::New(); + vtkSmartPointer colors = vtkSmartPointer::New(); + colors->SetNumberOfComponents(3); + colors->SetName("Colors"); + + //vtk sphere for sources + //vtkSmartPointer sphereSource = vtkSmartPointer::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 linesPolyData = vtkSmartPointer::New(); + + // Set the points, lines, and colors + linesPolyData->SetPoints(pts); + linesPolyData->SetLines(lines); + linesPolyData->GetCellData()->SetScalars(colors); + + // Write the file + vtkSmartPointer writer = vtkSmartPointer::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 mapper = + vtkSmartPointer::New(); + #if VTK_MAJOR_VERSION <= 5 + mapper->SetInput(linesPolyData); + #else + mapper->SetInputData(linesPolyData); + #endif + vtkSmartPointer actor = vtkSmartPointer::New(); + actor->SetMapper(mapper); + + //sphere + vtkSmartPointer mapperSphere = vtkSmartPointer< + vtkPolyDataMapper>::New(); + mapperSphere->SetInputConnection(sphereSource->GetOutputPort()); + + vtkSmartPointer actorSphere = vtkSmartPointer::New(); + actorSphere->SetMapper(mapperSphere); + //end sphere + + vtkSmartPointer renderer = vtkSmartPointer::New(); + vtkSmartPointer renderWindow = vtkSmartPointer< + vtkRenderWindow>::New(); + renderWindow->AddRenderer(renderer); + vtkSmartPointer renderWindowInteractor = + vtkSmartPointer::New(); + renderWindowInteractor->SetRenderWindow(renderWindow); + renderer->AddActor(actorSphere); + renderer->AddActor(actor); + + renderWindow->Render(); + renderWindowInteractor->Start();*/ +} + +// ------------------------------------------------------------------------- +#include +#include +#include +#include + +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 : "<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: "<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 ReadInputFile(const char* filename) +{ + // Variables + std::ifstream infile(filename); + std::string line; + std::string folderpath_allResults; + vector 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(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, std::pair > > 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 points_common = vtkSmartPointer::New(); + vtkSmartPointer lines_common = vtkSmartPointer::New(); + vtkSmartPointer colors_common = vtkSmartPointer::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, std::pair > >::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 edge_A = (*it_edges).first; + std::pair 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 line = vtkSmartPointer::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 line = vtkSmartPointer::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 linesPolyData = vtkSmartPointer::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 writer = vtkSmartPointer::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, std::pair > >::iterator it_edges = vector_pair_edges_A_to_B.begin(); + for(; it_edges != vector_pair_edges_A_to_B.end(); ++it_edges) + { + + std::pair edge_A = (*it_edges).first; + std::pair 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 > map_A_to_B, std::map< unsigned int, std::vector > 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& pts,vtkSmartPointer& lines, vtkSmartPointer& 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 line = vtkSmartPointer::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$ -- 2.47.1