4 * Created on: May 12, 2014
8 #include "airwaysEdge.h"
14 Edge::Edge() : m_id(0), m_source(NULL), m_target(NULL), m_mark(false), m_angle(0), m_length(
15 0), m_eDistance(0), m_aRadius(0), m_minRadius(0), m_maxRadius(0)
20 Edge::Edge(Edge* edge) : m_id(0), m_source(NULL), m_target(NULL), m_mark(false), m_angle(0), m_length(
21 0), m_eDistance(0), m_aRadius(0), m_minRadius(0), m_maxRadius(0)
23 this->m_id = edge->m_id;
24 this->m_source = edge->m_source;
25 this->m_target = edge->m_target;
26 this->m_angle = edge->m_angle;
27 this->m_length = edge->m_length;
28 this->m_eDistance = edge->m_eDistance;
29 this->m_vec_pair_posVox_rad = edge->m_vec_pair_posVox_rad;
30 this->m_aRadius = edge->m_aRadius;
31 this->m_minRadius = edge->m_minRadius;
32 this->m_maxRadius = edge->m_maxRadius;
39 int Edge::GetAngle() const
44 double Edge::GetARadius() const
46 return this->m_aRadius;
49 double Edge::GetEDistance() const
51 return this->m_eDistance;
54 double Edge::GetLength() const
56 return this->m_length;
59 double Edge::GetMaxRadius() const
61 return this->m_maxRadius;
64 double Edge::GetMinRadius() const
66 return this->m_minRadius;
69 Node* Edge::GetSource() const
71 return this->m_source;
74 Node* Edge::GetTarget() const
76 return this->m_target;
79 const vec_pair_posVox_rad& Edge::GetEdgeInfo() const
81 return m_vec_pair_posVox_rad;
84 bool Edge::IsMarked() const
89 bool Edge::IsPointInfluencedByEdge(float point_x, float point_y, float point_z) const
92 float minDistance = -1.0;
94 float distanceTemp = 0.0;
95 bool influenced = false;
96 Vec3 vector(point_x, point_y, point_z);
98 // Iterate over all voxel in the edge
99 vec_pair_posVox_rad::const_iterator it = this->m_vec_pair_posVox_rad.begin();
100 for (; it != this->m_vec_pair_posVox_rad.end() && !influenced; ++it)
102 // Get the vector from the compared source to each voxel
103 Vec3 voxel_actual = ((*it).first);
104 //Vec3 source_actual = this->m_source->GetCoords( );
105 //Vec3 vector_actual(voxel_actual - source_actual);
107 // Get the difference vector
108 Vec3 vector_diff = vector - voxel_actual;
111 distanceTemp = vector_diff.Norm();
116 minDistance = distanceTemp;
118 if(minDistance < ((*it).second*1.8))
121 else if(distanceTemp < minDistance)
123 minDistance = distanceTemp;
124 if(minDistance < ((*it).second)*1.8)
132 void Edge::SetAngle(const double& angle)
134 this->m_angle = angle;
137 void Edge::SetARadius(const double& aRadius)
139 this->m_aRadius = aRadius;
142 void Edge::SetEDistance(const double& eDistance)
144 this->m_eDistance = eDistance;
147 void Edge::SetLength(const double& length)
149 this->m_length = length;
152 void Edge::SetMaxRadius(const unsigned int& maxRadius)
154 this->m_maxRadius = maxRadius;
157 void Edge::SetMinRadius(const unsigned int& minRadius)
159 this->m_minRadius = minRadius;
162 void Edge::SetSource(Node* source)
164 this->m_source = source;
167 void Edge::SetTarget(Node* target)
169 this->m_target = target;
172 void Edge::AddSkeletonPairInfo(const pair_posVox_rad& skPairInfo)
174 this->m_vec_pair_posVox_rad.push_back(skPairInfo);
177 void Edge::SetSkeletonPairVector(const vec_pair_posVox_rad& skPInfoVector)
179 this->m_vec_pair_posVox_rad = skPInfoVector;
183 void Edge::UpdateEdgeInfo()
185 if (this->m_vec_pair_posVox_rad.empty())
187 std::cout << " ------------------------ This edge has no information" << std::endl;
188 std::cout << "Source:" << this->m_source->GetCoords() << std::endl;
189 std::cout << "Target:" << this->m_target->GetCoords() << std::endl;
192 double min = std::numeric_limits<unsigned int>::max();
193 double max = std::numeric_limits<unsigned int>::min();
195 vec_pair_posVox_rad::iterator it = this->m_vec_pair_posVox_rad.begin();
196 for (; it != this->m_vec_pair_posVox_rad.end(); ++it)
198 pair_posVox_rad info = *it;
199 average += info.second;
200 if (info.second < min)
202 else if (info.second > max)
205 this->m_aRadius = average / this->m_vec_pair_posVox_rad.size();
206 this->m_minRadius = min;
207 this->m_maxRadius = max;
210 double Edge::CompareWith(Edge* edge)
212 //std::cout << "Comparing ... " << std::endl;
213 // To compared the edges, the distance in each component
214 // is calculated and a distance function between edges is calculated.
215 // The components are:
216 // - Each component of the quaternion (son-parent quaterion). (4 components)
218 // - Edge average radius
220 // Difference in each quaternion component
230 /*const Edge* edgeActualParent = this->m_source->GetEdge( );
231 Vec3 sourceActualParent = edgeActualParent->m_source->GetCoords( );
232 Vec3 targetActualParent = edgeActualParent->m_target->GetCoords( );
233 Vec3 sourceActual = this->m_source->GetCoords( );
234 Vec3 targetActual = this->m_target->GetCoords( );
236 Vec3 sp_actual( sourceActualParent - targetActualParent );
237 Vec3 sg_actual( targetActual - sourceActual );
239 const Edge* edgeComparedParent = edge.m_source->GetEdge( );
240 Vec3 sourceComparedParent = edgeComparedParent->m_source->GetCoords( );
241 Vec3 targetComparedParent = edgeComparedParent->m_target->GetCoords( );
242 Vec3 sourceCompared = edge.m_source->GetCoords( );
243 Vec3 targetCompared = edge.m_target->GetCoords( );
245 Vec3 sp_compared( sourceComparedParent - targetComparedParent );
246 Vec3 sg_compared( targetCompared - sourceCompared );
248 Quaternion q_actual(sp_actual, sg_actual);
249 Quaternion q_compared(sp_compared, sg_compared);
251 float dif_r = q_actual.getR() - q_compared.getR();
252 float dif_i = q_actual.getI() - q_compared.getI();
253 float dif_j = q_actual.getJ() - q_compared.getJ();
254 float dif_k = q_actual.getK() - q_compared.getK();
255 float dif_distance = this->m_length - edge.m_length;
256 float dif_radius = this->m_aRadius - edge.m_aRadius;
258 double distanceEdge = GetDistanceToEdge(edge);
260 //std::cout << "Comparing ... OK" << std::endl;
265 float Edge::GetDistanceToEdge(Edge* edge)
267 // a. The parent edges must be aligned before calculate the distance.
268 // b. Parent end-points are aligned also.
269 // Steps a and b represent a rotation and a translation respectively.
270 // c. For each point in the actual Edge (this) the closest point in the
271 // compared edge is found. This distance is calculated and then averaged
272 // for all the points. Must be done in both ways?
273 // We can take the maximum of both distances max(D(a,b),D(b,a))
276 float averageDistance = 0.0;
279 // Get the parent Edges
280 /*const Edge* edgeActualParent = this->m_source->GetEdge( );
281 Vec3 sourceActualParent = edgeActualParent->m_source->GetCoords( );
282 Vec3 targetActualParent = edgeActualParent->m_target->GetCoords( );
283 Vec3 sourceActual = this->m_source->GetCoords( );
284 Vec3 targetActual = this->m_target->GetCoords( );
286 Vec3 sp_actual( sourceActualParent - targetActualParent );
287 Vec3 sg_actual( targetActual - sourceActual );
289 const Edge* edgeComparedParent = edge.m_source->GetEdge( );
290 Vec3 sourceComparedParent = edgeComparedParent->m_source->GetCoords( );
291 Vec3 targetComparedParent = edgeComparedParent->m_target->GetCoords( );
292 Vec3 sourceCompared = edge.m_source->GetCoords( );
293 Vec3 targetCompared = edge.m_target->GetCoords( );
295 Vec3 sp_compared( sourceComparedParent - targetComparedParent );
296 Vec3 sg_compared( targetCompared - sourceCompared );
298 //Get the quaternion between parents
299 Quaternion q_parent(sp_actual, sp_compared);
301 // Get the translation vector to translate the compared parent to the actual parent
302 const Vec3 transParents( targetActualParent - targetComparedParent);
305 // Iterate over the voxels in the compared edge
306 for(vec_pair_posVox_rad::iterator it = edge->m_vec_pair_posVox_rad.begin(); it != edge->m_vec_pair_posVox_rad.end(); ++it)
308 // Get the vector from the compared source to each voxel in the compared edge
309 Vec3 voxel_actual = (*it).first;
310 /*Vec3 vector_compared(voxel_actual - sourceCompared);
312 // Rotate the vector using the parental quaternion
313 Vec3 vector_compared_rotated = q_parent.rotateVector(vector_compared);
316 // Find the distance to the closest point in the actual edge
317 //averageDistance += getDistanceToClosestPoint(sourceActual+vector_compared_rotated);
318 averageDistance += this->GetSmallestDistanceToPoint(voxel_actual);
323 // Get the average value
325 averageDistance = averageDistance / numPoints;
327 return averageDistance;
330 float Edge::GetDistanceToTranslatedEdge(Edge* edge, Vec3 vector_translation)
332 // a. The parent edges must be aligned before calculate the distance.
333 // b. Parent end-points are aligned also.
334 // Steps a and b represent a rotation and a translation respectively.
335 // c. For each point in the actual Edge (this) the closest point in the
336 // compared edge is found. This distance is calculated and then averaged
337 // for all the points. Must be done in both ways?
338 // We can take the maximum of both distances max(D(a,b),D(b,a))
341 float averageDistance = 0.0;
344 // Iterate over the voxels in the compared edge
345 for(vec_pair_posVox_rad::iterator it = edge->m_vec_pair_posVox_rad.begin(); it != edge->m_vec_pair_posVox_rad.end(); ++it)
347 // Get the vector from the compared source to each voxel in the compared edge
348 Vec3 voxel_actual = (*it).first;
350 // Translate the voxel
351 Vec3 voxel_translated = voxel_actual + vector_translation;
353 // Find the distance to the closest point in the actual edge
354 averageDistance += this->GetSmallestDistanceToPoint(voxel_translated);
359 // Get the average value
361 averageDistance = averageDistance / numPoints;
363 return averageDistance;
366 float Edge::GetDistanceWeigthedToTranslatedEdge(Edge* edge, Vec3 vector_translation)
368 // a. The parent edges must be aligned before calculate the distance.
369 // b. Parent end-points are aligned also.
370 // Steps a and b represent a rotation and a translation respectively.
371 // c. For each point in the actual Edge (this) the closest point in the
372 // compared edge is found. This distance is calculated and then averaged
373 // for all the points. Must be done in both ways?
374 // We can take the maximum of both distances max(D(a,b),D(b,a))
377 double alfa = 4; // Weight for duplicated closest points
378 float distance_average = 0.0;
379 float distance_actual = 0.0;
381 std::map<Vec3, vec_pair_posVox_rad > map_vector_pair_vector_distance; // Map to save the correspondences between the given edge voxels and this voxels, and the distance between them.
383 // Iterate over the voxels in the compared edge
384 for(vec_pair_posVox_rad::iterator it = edge->m_vec_pair_posVox_rad.begin(); it != edge->m_vec_pair_posVox_rad.end(); ++it)
386 // Get the vector from the compared source to each voxel in the compared edge
387 Vec3 voxel_actual = (*it).first;
389 // Translate the voxel
390 Vec3 voxel_translated = voxel_actual + vector_translation;
392 // Find the distance to the closest point in the actual edge
393 Vec3 point_closest = this->GetClosestPoint(voxel_translated, distance_actual);
394 distance_average += distance_actual;
397 pair_posVox_rad pair_vector_distance(voxel_actual, distance_actual);
398 map_vector_pair_vector_distance[point_closest].push_back(pair_vector_distance);
403 // Add the weight to shared closest points
404 float distance_weighted = 0.0;
405 for(std::map<Vec3, vec_pair_posVox_rad>::iterator it = map_vector_pair_vector_distance.begin(); it != map_vector_pair_vector_distance.end(); ++it)
407 vec_pair_posVox_rad vector_pair_vector_distance = (*it).second;
408 if(vector_pair_vector_distance.size() > 1)
410 for(vec_pair_posVox_rad::iterator it_vector = vector_pair_vector_distance.begin(); it_vector != vector_pair_vector_distance.end(); ++it_vector)
411 distance_weighted += (*it_vector).second * alfa;
415 for(vec_pair_posVox_rad::iterator it_vector = vector_pair_vector_distance.begin(); it_vector != vector_pair_vector_distance.end(); ++it_vector)
416 distance_weighted += (*it_vector).second;
420 // Get the average value
423 distance_average = distance_average / numPoints;
424 distance_weighted = distance_weighted / numPoints;
427 //std::cout << "[DistanceW;DistAvg][" << distance_weighted << ";" << distance_average << "];;";
428 //return distance_average;
429 return distance_weighted;
432 float Edge::GetDistanceToPoint(Vec3 point)
436 float distance = 0.0;
438 // Iterate over all voxel in the edge
439 vec_pair_posVox_rad::iterator it = this->m_vec_pair_posVox_rad.begin();
440 for (; it != this->m_vec_pair_posVox_rad.end(); ++it)
442 // Get the vector from the compared source to each voxel
443 Vec3 voxel_actual = ((*it).first);
445 // Get the difference vector
446 Vec3 vector_diff = point - voxel_actual;
449 distance += vector_diff.Norm();
456 float Edge::GetSmallestDistanceToPoint(Vec3 vector)
459 float minDistance = -1.0;
461 float distanceTemp = 0.0;
463 // Iterate over all voxel in the edge
464 vec_pair_posVox_rad::iterator it = this->m_vec_pair_posVox_rad.begin();
465 for (; it != this->m_vec_pair_posVox_rad.end(); ++it)
467 // Get the vector from the compared source to each voxel
468 Vec3 voxel_actual = ((*it).first);
469 //Vec3 source_actual = this->m_source->GetCoords( );
470 //Vec3 vector_actual(voxel_actual - source_actual);
472 // Get the difference vector
473 Vec3 vector_diff = vector - voxel_actual;
476 distanceTemp = vector_diff.Norm();
481 minDistance = distanceTemp;
484 else if(distanceTemp < minDistance)
485 minDistance = distanceTemp;
491 Vec3 Edge::GetClosestPoint(Vec3 vector, float& distance)
496 float distanceTemp = 0.0;
498 // Iterate over all voxel in the edge
499 vec_pair_posVox_rad::iterator it = this->m_vec_pair_posVox_rad.begin();
500 for (; it != this->m_vec_pair_posVox_rad.end(); ++it)
502 // Get the vector from the compared source to each voxel
503 Vec3 voxel_actual = ((*it).first);
504 //Vec3 source_actual = this->m_source->GetCoords( );
505 //Vec3 vector_actual(voxel_actual - source_actual);
507 // Get the difference vector
508 Vec3 vector_diff = vector - voxel_actual;
511 distanceTemp = vector_diff.Norm();
516 point_closest = voxel_actual;
517 distance = distanceTemp;
520 else if(distanceTemp < distance)
522 point_closest = voxel_actual;
523 distance = distanceTemp;
527 return point_closest;
530 Edge* Edge::ConcatenateToSuperior(Edge* superior)
533 return new Edge(this);
535 // Check concatenation condition in the nodes
536 if(superior->GetTarget()->GetCoords() != this->GetSource()->GetCoords())
538 std::cout << "Edge::ConcatenateToSuperior - superior.target != actualEdge.source. Actual idNode:" << this->m_id << std::endl;
539 std::cout << "[superior.target;actualEdge.source]: [" << superior->GetTarget()->GetCoords() << ";" << this->GetSource()->GetCoords() << "]" << std::endl;
543 // Check concatenation condition in the edge information
544 vec_pair_posVox_rad vectorInfo_thisEdge = this->GetEdgeInfo(); // Actual edge information
545 vec_pair_posVox_rad vectorInfo_superiorEdge = superior->GetEdgeInfo(); // Superior edge information
547 // Check that last superior edge position is equal to first initial actual edge position
548 if(vectorInfo_superiorEdge[vectorInfo_superiorEdge.size()-1].first != vectorInfo_thisEdge[0].first)
550 std::cout << "Edge::ConcatenateToSuperior - superior.info[end] != actualEdge.info[begin]. Actual idNode:" << this->m_id << std::endl;
554 // Change the superior edge target, the superior length, and the euclidian distance
555 superior->SetTarget(this->GetTarget());
556 superior->SetLength(superior->GetLength()+this->GetLength()-1);
557 superior->SetEDistance(superior->GetEDistance()+this->GetEDistance());
559 // Add the position information
560 vec_pair_posVox_rad::iterator it_actualInfo = vectorInfo_thisEdge.begin();
561 ++it_actualInfo; // Skip the first position because is it shared between both edges
562 for(; it_actualInfo != vectorInfo_thisEdge.end(); ++it_actualInfo)
563 vectorInfo_superiorEdge.push_back((*it_actualInfo));
565 // Set the new information
566 superior->SetSkeletonPairVector(vectorInfo_superiorEdge);
572 Edge& Edge::operator=(Edge& edge)
574 this->m_mark = edge.m_mark;
575 this->m_angle = edge.m_angle;
576 this->m_length = edge.m_length;
577 this->m_eDistance = edge.m_eDistance;
578 this->m_aRadius = edge.m_aRadius;
579 this->m_minRadius = edge.m_minRadius;
580 this->m_maxRadius = edge.m_maxRadius;
581 this->m_source = edge.m_source;
582 this->m_target = edge.m_target;
583 this->m_vec_pair_posVox_rad = edge.m_vec_pair_posVox_rad;
587 const Edge& Edge::operator=(const Edge& edge)
589 this->m_mark = edge.m_mark;
590 this->m_angle = edge.m_angle;
591 this->m_length = edge.m_length;
592 this->m_eDistance = edge.m_eDistance;
593 this->m_aRadius = edge.m_aRadius;
594 this->m_minRadius = edge.m_minRadius;
595 this->m_maxRadius = edge.m_maxRadius;
596 this->m_source = edge.m_source;
597 this->m_target = edge.m_target;
598 this->m_vec_pair_posVox_rad = edge.m_vec_pair_posVox_rad;
602 bool Edge::operator==(const Edge& edge)
610 /*if (this->m_angle != edge.m_angle)
613 this->m_angle >= edge.m_angle ? this->m_angle : edge.m_angle;
614 double tol = 0.3 * max;
615 double diff = this->m_angle - edge.m_angle;
616 if (!((-tol <= diff) && (diff <= tol)))
619 if (this->m_length != edge.m_length)
622 this->m_length >= edge.m_length ? this->m_length : edge.m_length;
623 double tol = 0.3 * max;
624 double diff = abs(this->m_length - edge.m_length);
628 if (this->m_eDistance != edge.m_eDistance)
631 this->m_eDistance >= edge.m_eDistance ?
632 this->m_eDistance : edge.m_eDistance;
633 double tol = 0.3 * max;
634 double diff = abs(this->m_eDistance - edge.m_eDistance);
638 /*if (this->m_aRadius != edge.m_aRadius)
641 this->m_aRadius >= edge.m_aRadius ?
642 this->m_aRadius : edge.m_aRadius;
643 double tol = 0.3 * max;
644 double diff = abs(this->m_aRadius - edge.m_aRadius);
648 if (this->m_minRadius != edge.m_minRadius)
651 this->m_minRadius >= edge.m_minRadius ?
652 this->m_minRadius : edge.m_minRadius;
653 double tol = 0.3 * max;
654 double diff = abs(this->m_minRadius - edge.m_minRadius);
658 if (this->m_maxRadius != edge.m_maxRadius)
661 this->m_maxRadius >= edge.m_maxRadius ?
662 this->m_maxRadius : edge.m_maxRadius;
663 double tol = 0.3 * max;
664 double diff = abs(this->m_maxRadius - edge.m_maxRadius);
668 return (cmp1 && cmp2 && cmp3 && cmp4 && cmp5 && cmp6);
671 bool Edge::operator!=(const Edge& edge)
673 return !(*this == edge);
676 bool Edge::operator<(const Edge& edge)
680 if (this->m_length != edge.m_length)
683 this->m_length >= edge.m_length ? this->m_length : edge.m_length;
684 double tol = 0.3 * max;
685 double diff = abs(this->m_length - edge.m_length);
689 if (this->m_eDistance != edge.m_eDistance)
692 this->m_eDistance >= edge.m_eDistance ?
693 this->m_eDistance : edge.m_eDistance;
694 double tol = 0.3 * max;
695 double diff = abs(this->m_eDistance - edge.m_eDistance);
701 if ((this->m_length < edge.m_length)
702 || (this->m_eDistance < edge.m_eDistance))
708 bool Edge::operator>(const Edge& edge)
712 if (this->m_length != edge.m_length)
715 this->m_length >= edge.m_length ? this->m_length : edge.m_length;
716 double tol = 0.3 * max;
717 double diff = abs(this->m_length - edge.m_length);
721 if (this->m_eDistance != edge.m_eDistance)
724 this->m_eDistance >= edge.m_eDistance ?
725 this->m_eDistance : edge.m_eDistance;
726 double tol = 0.3 * max;
727 double diff = abs(this->m_eDistance - edge.m_eDistance);
733 if ((this->m_length > edge.m_length)
734 || (this->m_eDistance > edge.m_eDistance))
739 } /* namespace airways */