void ShowNPoints_Tools::SeparateSplines()
{
WidgetShowNPoints* wsp = bbGetInputWidgetShowNPoints();
- if((wsp->GetLstModelShowNPointsSize() > 0) && (bbGetInputParams().size() == 3))
+ if((wsp->GetLstModelShowNPointsSize() > 0) && (bbGetInputParams().size() == 6))
{
- std::vector<double> refPoint = bbGetInputParams();
- std::vector<double> pointsX, pointsY, pointsZ;
+ std::vector<double> refPoint = {bbGetInputParams()[0], bbGetInputParams()[1], bbGetInputParams()[2]};
+ std::vector<double> pointsX, pointsY, pointsZ, refPlaneNorm, planeNormal;
std::vector<std::vector<std::vector<double>>> newGroups;
std::vector<std::vector<double>> planesNorm;
- double rotAxis[3] = {0,0,1};
+ std::vector<double> rotAxis = {bbGetInputParams()[3], bbGetInputParams()[4], bbGetInputParams()[5]}; ;
- double vector1[3], vector2[3], planeNormal[3], rotatedNormal[3], axisPoint[3], refPlaneNorm[3];
- int idPoint2, idPoint3;
+ std::vector<double> intersectionTop, intersectionBottom;
int numGroups = wsp->GetLstModelShowNPointsSize();
for(int i = 0; i < numGroups; i++){
pointsY = wsp->GetModelShowNPoints(i)->GetLstPointsY();
pointsZ = wsp->GetModelShowNPoints(i)->GetLstPointsZ();
- idPoint2 = ((int)pointsX.size()/3)-1;
- idPoint3 = idPoint2*2;
-
- vector1[0] = pointsX[idPoint2]-refPoint[0];
- vector1[1] = pointsY[idPoint2]-refPoint[1];
- vector1[2] = pointsZ[idPoint2]-refPoint[2];
-
- vector2[0] = pointsX[idPoint3]-refPoint[0];
- vector2[1] = pointsY[idPoint3]-refPoint[1];
- vector2[2] = pointsZ[idPoint3]-refPoint[2];
-
- vtkMath::Cross(vector1, vector2, planeNormal);
- vtkMath::Normalize(planeNormal);
- vtkMath::Cross(planeNormal, rotAxis, rotatedNormal);
- vtkMath::Normalize(rotatedNormal);
-
+ planeNormal = GetPlaneNormalFromPointsRefPoint(pointsX, pointsY, pointsZ, refPoint);
+
if(i == 0){
- std::copy(std::begin(planeNormal), std::end(planeNormal), std::begin(refPlaneNorm));
+ refPlaneNorm = planeNormal;
}
- std::vector<std::vector<double>> pointsOver, pointsUnder, pointsOn;
- std::vector<double> point(3);
- for(int j = 0; j < pointsX.size(); j++){
- point[0] = pointsX[j];
- point[1] = pointsY[j];
- point[2] = pointsZ[j];
- vtkMath::Subtract(point.data(), refPoint.data(), vector1);
- if(vtkMath::Dot(vector1, rotatedNormal) > 0){
- pointsOver.push_back(point);
+ //Get points over, under and on plane, all groups ordered in the direction of the axis.
+ std::vector<std::vector<std::vector<double>>> pointsAroundPlane;
+ pointsAroundPlane = GetOrderedPointsAroundPlane(refPoint, planeNormal, rotAxis, pointsX, pointsY, pointsZ);
+ std::vector<std::vector<double>> pointsOver = pointsAroundPlane[0];
+ std::vector<std::vector<double>> pointsUnder = pointsAroundPlane[1];
+ std::vector<std::vector<double>> pointsOn = pointsAroundPlane[2];
+
+ //Calculate points where splines should intersect
+ double middlePoint[3];
+ if(intersectionTop.empty() && intersectionBottom.empty()){
+ if(!pointsOver.empty() && !pointsUnder.empty()){
+ vtkMath::Add(pointsOver.back().data(), pointsUnder.back().data(), middlePoint);
+ vtkMath::MultiplyScalar(middlePoint, 0.5);
+ intersectionTop = GetProjectionPointOnAxis(middlePoint, refPoint.data(), rotAxis.data());
+
+ vtkMath::Add(pointsOver.front().data(), pointsUnder.front().data(), middlePoint);
+ vtkMath::MultiplyScalar(middlePoint, 0.5);
+ intersectionBottom = GetProjectionPointOnAxis(middlePoint, refPoint.data(), rotAxis.data());
+ }
+ else if(!pointsOn.empty()){
+ intersectionTop = pointsOn.back();
+ intersectionBottom = pointsOn.front();
}
- else if (vtkMath::Dot(vector1, rotatedNormal) < 0){
- pointsUnder.push_back(point);
+ else{
+ intersectionTop = pointsOver.empty()?pointsUnder.back():pointsOver.back();
+ intersectionBottom = pointsOver.empty()?pointsUnder.front():pointsOver.front();
}
- //else{
- // pointsOn.push_back(j);
- //}
}
-
- std::vector<double> normalToSave(std::begin(planeNormal), std::end(planeNormal));
- if(!pointsOver.empty()){
+
+ //Save groups and their normals
+ if(!pointsOver.empty()){
+ pointsOver.push_back(intersectionTop);
+ pointsOver.insert(pointsOver.begin(), intersectionBottom);
newGroups.push_back(pointsOver);
- planesNorm.push_back(normalToSave);
+ planesNorm.push_back(planeNormal);
}
- if(!pointsUnder.empty()){
- std::reverse(pointsUnder.begin(), pointsUnder.end());
+ if(!pointsUnder.empty()){
+ pointsUnder.push_back(intersectionTop);
+ pointsUnder.insert(pointsUnder.begin(), intersectionBottom);
newGroups.push_back(pointsUnder);
- vtkMath::MultiplyScalar(planeNormal, -1);
- normalToSave[0] = planeNormal[0]; normalToSave[1] = planeNormal[1]; normalToSave[2] = planeNormal[2];
- planesNorm.push_back(normalToSave);
+ double tempPlaneNorm[3] = {planeNormal[0], planeNormal[1], planeNormal[2]};
+ vtkMath::MultiplyScalar(tempPlaneNorm, -1);
+ planeNormal[0] = tempPlaneNorm[0]; planeNormal[1] = tempPlaneNorm[1]; planeNormal[2] = tempPlaneNorm[2];
+ planesNorm.push_back(planeNormal);
}
-
-
}
+ //Sort Groups by their plane angles, from the first plane around the rotation axis
auto comp = [&](int id1,int id2)-> bool {
- //double dot = vtkMath::Dot(planesNorm[id1].data(), refPlaneNorm);
- //double det = vtkMath::Determinant3x3(planesNorm[id1].data(), refPlaneNorm, rotAxis);
- //double angle1 = atan2(det, dot);
- double angle1 = vtkMath::SignedAngleBetweenVectors(planesNorm[id1].data(), refPlaneNorm, rotAxis);
+ double angle1 = vtkMath::SignedAngleBetweenVectors(planesNorm[id1].data(), refPlaneNorm.data(), rotAxis.data());
if(angle1 < 0) angle1 += 2*vtkMath::Pi();
- double angle2 = vtkMath::SignedAngleBetweenVectors(planesNorm[id2].data(), refPlaneNorm, rotAxis);
+ double angle2 = vtkMath::SignedAngleBetweenVectors(planesNorm[id2].data(), refPlaneNorm.data(), rotAxis.data());
if(angle2 < 0) angle2 += 2*vtkMath::Pi();
return angle1 < angle2;
};
std::iota(std::begin(indexVect), std::end(indexVect), 0);
std::sort(indexVect.begin(), indexVect.end(), comp);
+ //Build Groups
wsp->ResetCollections_();
for(int i = 0; i < newGroups.size(); i++){
std::vector<std::vector<double>> group = newGroups[indexVect[i]];
wsp->InsertCollectionAfter_();
}
}
+ wsp->InvertLstPoints_();
wsp->SetOutputBox();
wsp->UndoRedo_SaveCollection();
}
}
+std::vector<std::vector<std::vector<double>>> ShowNPoints_Tools::GetOrderedPointsAroundPlane(std::vector<double> planeOrigin, std::vector<double> planeNormal, std::vector<double> rotAxis, std::vector<double> pX, std::vector<double> pY, std::vector<double> pZ){
+
+ std::vector<std::vector<double>> pointsOver, pointsUnder, pointsOn;
+ std::vector<double> point(3);
+ double rotatedNormal[3], vector1[3];
+ vtkMath::Cross(planeNormal.data(), rotAxis.data(), rotatedNormal);
+ vtkMath::Normalize(rotatedNormal);
+
+ bool changedSide = false;
+ int startingSide = -1; //-1 not started, 0 over, 1 under
+ int placedElements = 0;
+ for(int j = 0; j < pX.size(); j++){
+ point[0] = pX[j];
+ point[1] = pY[j];
+ point[2] = pZ[j];
+ vtkMath::Subtract(point.data(), planeOrigin.data(), vector1);
+ if(vtkMath::Dot(vector1, rotatedNormal) > 0){
+ if(!pointsUnder.empty()){
+ changedSide = true;
+ }
+ else {
+ startingSide = 0;
+ }
+
+ if(changedSide && startingSide == 0){
+ pointsOver.insert(pointsOver.begin()+placedElements, point);
+ placedElements++;
+ }
+ else{
+ pointsOver.push_back(point);
+ }
+ }
+ else if (vtkMath::Dot(vector1, rotatedNormal) < 0){
+ if(!pointsOver.empty())
+ {
+ changedSide = true;
+ }
+ else startingSide = 1;
+
+ if(changedSide && startingSide == 1){
+ pointsUnder.insert(pointsUnder.begin()+placedElements, point);
+ placedElements++;
+ }
+ else{
+ pointsUnder.push_back(point);
+ }
+ }
+ else{
+ pointsOn.push_back(point);
+ }
+ }
+ ////////
+ /// Check direction of points and reverse if needed.
+ /// Always follows axis vector direction
+ ////////
+ if(pointsUnder.size() > 2){
+ vtkMath::Subtract(pointsUnder.front(), pointsUnder.back(), vector1);
+ double dot = vtkMath::Dot(vector1, rotAxis.data());
+ if(dot < 0){
+ std::reverse(pointsUnder.begin(), pointsUnder.end());
+ if(startingSide == 1 && !pointsOn.empty()){
+ std::reverse(pointsOn.begin(), pointsOn.end());
+ }
+ }
+ }
+ if(pointsOver.size() > 2){
+ vtkMath::Subtract(pointsOver.front(), pointsOver.back(), vector1);
+ double dot = vtkMath::Dot(vector1, rotAxis.data());
+ if(dot < 0){
+ std::reverse(pointsOver.begin(), pointsOver.end());
+ if(startingSide == 0 && !pointsOn.empty()){
+ std::reverse(pointsOn.begin(), pointsOn.end());
+ }
+ }
+ }
+
+ std::vector<std::vector<std::vector<double>>> results = {pointsOver, pointsUnder, pointsOn};
+ return results;
+}
+
+std::vector<double> ShowNPoints_Tools::GetPlaneNormalFromPointsRefPoint(std::vector<double> pointsX, std::vector<double> pointsY, std::vector<double> pointsZ, std::vector<double> refPoint)
+{
+ double vector1[3], vector2[3], planeNormal[3];
+ int idPoint2, idPoint3;
+ idPoint2 = ((int)pointsX.size()/3)-1;
+ idPoint3 = ((int)pointsX.size()/2)+1;
+
+ vector1[0] = pointsX[idPoint2]-refPoint[0];
+ vector1[1] = pointsY[idPoint2]-refPoint[1];
+ vector1[2] = pointsZ[idPoint2]-refPoint[2];
+
+ vector2[0] = pointsX[idPoint3]-refPoint[0];
+ vector2[1] = pointsY[idPoint3]-refPoint[1];
+ vector2[2] = pointsZ[idPoint3]-refPoint[2];
+
+ vtkMath::Cross(vector1, vector2, planeNormal);
+ vtkMath::Normalize(planeNormal);
+
+ std::vector<double> result(std::begin(planeNormal), std::end(planeNormal));
+ return result;
+}
+
+std::vector<double> ShowNPoints_Tools::GetProjectionPointOnAxis(double pointToProject[3], double originAxis[3], double axisDir[3])
+{
+ double projectionVect[3], vectorToPoint[3], pointProjected[3];
+ vtkMath::Subtract(pointToProject, originAxis, vectorToPoint);
+
+ vtkMath::ProjectVector(vectorToPoint, axisDir, projectionVect);
+ vtkMath::Add(originAxis, projectionVect, pointProjected);
+
+ std::vector<double> result(std::begin(pointProjected), std::end(pointProjected));
+ return result;
+}
+
+
//=====
// Before editing this file, make sure it's a file of your own (i.e.: it wasn't generated from xml description; if so : your modifications will be lost)
//=====
git://git.creatis.insa-lyon.fr / creaMaracasVisu.git / blobdiff