903 lines
28 KiB
C++
903 lines
28 KiB
C++
////////////////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Paella
|
|
// Copyright (C) 2015 - Thomas FORGIONE, Emilie JALRAS, Marion LENFANT, Thierry MALON, Amandine PAILLOUX
|
|
// Authors :
|
|
// Thomas FORGIONE
|
|
// Emilie JALRAS
|
|
// Marion LENFANT
|
|
// Thierry MALON
|
|
// Amandine PAILLOUX
|
|
//
|
|
// This file is part of the project Paella
|
|
// This software is provided 'as-is', without any express or implied warranty.
|
|
// In no event will the authors be held liable for any damages arising from the use of this software.
|
|
//
|
|
// Permission is granted to anyone to use this software for any purpose,
|
|
// including commercial applications, and to alter it and redistribute it freely,
|
|
// subject to the following restrictions:
|
|
//
|
|
// 1. The origin of this software must not be misrepresented;
|
|
// you must not claim that you wrote the original software.
|
|
// If you use this software in a product, an acknowledgment
|
|
// in the product documentation would be appreciated but is not required.
|
|
//
|
|
// 2. Altered source versions must be plainly marked as such,
|
|
// and must not be misrepresented as being the original software.
|
|
//
|
|
// 3. This notice may not be removed or altered from any source distribution.
|
|
////////////////////////////////////////////////////////////////////////////////
|
|
#include <iostream>
|
|
#include <fstream>
|
|
#include <sstream>
|
|
#include <cmath>
|
|
#include <list>
|
|
#include <algorithm>
|
|
|
|
#include "Meshing/Skeleton3D.hpp"
|
|
#include "Utility/TupleFunctions.hpp"
|
|
#include "Geometry/MathFunctions.hpp"
|
|
#include "Geometry/Point.hpp"
|
|
|
|
namespace pae
|
|
{
|
|
|
|
namespace detail
|
|
{
|
|
|
|
// template<typename T>
|
|
// std::ostream& operator<<(std::ostream& out, std::pair<std::vector<std::list<geo::Point<T>>>, std::vector<std::list<geo::Point<T>>>> const& pair)
|
|
// {
|
|
// unsigned int i = 0;
|
|
// out << "up_points_vector_list" << std::endl;
|
|
// for (auto const& l : pair.first)
|
|
// {
|
|
// out << "circle_" << i << std::endl;
|
|
// for (auto const& p : l)
|
|
// {
|
|
// out << p.first << " ";
|
|
// }
|
|
// out << std::endl;
|
|
// i++;
|
|
// }
|
|
//
|
|
// out << std::endl << "Second list" << std::endl;
|
|
// i = 0;
|
|
// for (auto const& l : pair.second)
|
|
// {
|
|
//
|
|
// out << "circle_" << i << std::endl;
|
|
// for (auto const& p : l)
|
|
// {
|
|
// out << p.first << " ";
|
|
// }
|
|
// out << std::endl;
|
|
// i++;
|
|
// }
|
|
//
|
|
//
|
|
// return out;
|
|
// }
|
|
|
|
} // namespace detail
|
|
|
|
Skeleton3D::Skeleton3D() : splines{}
|
|
{
|
|
|
|
}
|
|
|
|
bool Skeleton3D::loadFromFile(std::string const& path)
|
|
{
|
|
std::fstream file{path};
|
|
|
|
if (!file)
|
|
{
|
|
std::cerr << "Warning : couldn't load Skeleton from file " << path << std::endl;
|
|
return false;
|
|
}
|
|
|
|
std::stringstream stream;
|
|
std::string line;
|
|
bool alreadyFoundL = false;
|
|
geo::Spline<geo::Vector3f, float> currentSpline;
|
|
|
|
while (std::getline(file, line))
|
|
{
|
|
if (line.size() == 0)
|
|
continue;
|
|
|
|
char first;
|
|
stream.clear();
|
|
stream.str(line);
|
|
|
|
stream >> first;
|
|
|
|
switch (first)
|
|
{
|
|
case 'd':
|
|
{
|
|
// Add the current spline
|
|
splines.push_back(currentSpline);
|
|
|
|
// Prepare the new one
|
|
currentSpline = geo::Spline<geo::Vector3f, float>{};
|
|
|
|
// Set the degree of the spline
|
|
int i;
|
|
stream >> i;
|
|
currentSpline.degree = i;
|
|
|
|
break;
|
|
}
|
|
|
|
case 'n':
|
|
{
|
|
float x;
|
|
while (stream >> x)
|
|
{
|
|
currentSpline.nodes.push_back(x);
|
|
}
|
|
break;
|
|
}
|
|
|
|
case 'p':
|
|
{
|
|
// Add control point to the spline
|
|
currentSpline.addControlPoint(line.substr(1));
|
|
break;
|
|
}
|
|
|
|
case 'l':
|
|
{
|
|
if (!alreadyFoundL)
|
|
{
|
|
splines.push_back(currentSpline);
|
|
alreadyFoundL = true;
|
|
}
|
|
|
|
float size_junction;
|
|
stream >> size_junction;
|
|
|
|
unsigned int x;
|
|
|
|
junctions.push_back(Junction{});
|
|
Junction& j = junctions[junctions.size()-1];
|
|
|
|
for (unsigned int i = 0; i<size_junction; i++)
|
|
{
|
|
stream >> x;
|
|
j.push_back(std::make_pair(x-1,0.0f));
|
|
}
|
|
|
|
float y;
|
|
for (unsigned int i = 0; i<size_junction; i++)
|
|
{
|
|
stream >> y;
|
|
j[i].second = y;
|
|
}
|
|
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
|
|
if (!alreadyFoundL)
|
|
{
|
|
splines.push_back(currentSpline);
|
|
}
|
|
|
|
// Remove the first spline that is empty
|
|
splines.erase(splines.begin());
|
|
|
|
return true;
|
|
|
|
}
|
|
|
|
void Skeleton3D::meshJunctions(geo::Mesh& mesh, std::vector<std::vector<geo::Circle<float>>> const& splines_circles) const
|
|
{
|
|
for(unsigned int i = 0; i < junctions.size(); i++)
|
|
{
|
|
geo::Vector3f junction_point; // sphere's center of the junction point
|
|
float junction_radius; // sphere's radius
|
|
|
|
// Get the junction elements
|
|
auto spline = splines[junctions[i][0].first];
|
|
auto tuple = spline(junctions[i][0].second);
|
|
|
|
junction_point = std::get<0>(tuple);
|
|
junction_radius = std::get<1>(tuple);
|
|
|
|
|
|
std::vector<geo::Circle<float>> junction_circles; // characteristic circles from the splines of the junction
|
|
for(unsigned int j = 0; j<junctions[i].size(); j++)
|
|
{
|
|
unsigned int indSpline = junctions[i][j].first;
|
|
float t = junctions[i][j].second;
|
|
auto spline_circles = splines_circles[indSpline];
|
|
|
|
if (t==0.0f)
|
|
{
|
|
junction_circles.push_back(spline_circles[0]);
|
|
}
|
|
else
|
|
{
|
|
junction_circles.push_back(spline_circles[spline_circles.size()-1]);
|
|
}
|
|
}
|
|
meshJunction(mesh,junction_point,junction_radius, junction_circles);
|
|
}
|
|
|
|
}
|
|
|
|
geo::Vector4f meshJunction(geo::Mesh& mesh, geo::Vector3f const& junction_point, float junction_radius, std::vector<geo::Circle<float>> const& junction_circles)
|
|
{
|
|
std::vector<geo::Vector3f> centerPoints;
|
|
|
|
for (auto const& circle : junction_circles)
|
|
{
|
|
centerPoints.push_back(circle.center);
|
|
}
|
|
// Find the medium plane that cut the circles in two
|
|
geo::Vector4f plane = geo::closestPlane<float>(centerPoints);
|
|
|
|
// Compute list of points upside the plane and downside
|
|
auto sorted_points = detail::sortedPoints(plane, junction_circles);
|
|
|
|
geo::Vector3f direction{plane.x(),plane.y(),plane.z()};
|
|
direction /= direction.norm();
|
|
geo::Vector3f up_point_proj = junction_point - direction*junction_radius;
|
|
geo::Vector3f down_point_proj = junction_point + direction*junction_radius;
|
|
|
|
mesh.vertices.push_back(up_point_proj);
|
|
mesh.vertices.push_back(down_point_proj);
|
|
|
|
// choose the indice of the two new points
|
|
unsigned int ind_up = mesh.vertices.size()-2;
|
|
unsigned int ind_down = mesh.vertices.size()-1;
|
|
|
|
// mesh the up part of the sphere
|
|
for (auto const& up_points : sorted_points.first)
|
|
{
|
|
auto last = std::end(up_points);
|
|
last--;
|
|
|
|
for(auto it = std::begin(up_points), end = std::end(up_points); it != end; ++it)
|
|
{
|
|
if (it != last)
|
|
{
|
|
auto cp = it;
|
|
cp ++;
|
|
detail::pushFaceOriented(mesh,
|
|
geo::Vector3u{it->first,cp->first,ind_up},
|
|
junction_point);
|
|
}
|
|
}
|
|
}
|
|
|
|
// mesh the down part of the sphere
|
|
for (auto const& down_points : sorted_points.second)
|
|
{
|
|
auto last = std::end(down_points);
|
|
last--;
|
|
|
|
for(auto it = std::begin(down_points), end = std::end(down_points); it != end; ++it)
|
|
{
|
|
if (it != last)
|
|
{
|
|
auto cp = it;
|
|
cp ++;
|
|
detail::pushFaceOriented(mesh,
|
|
geo::Vector3u{it->first,cp->first,ind_down},
|
|
junction_point);
|
|
}
|
|
}
|
|
}
|
|
|
|
// connect the four points around the plane
|
|
|
|
auto points_c1_up = sorted_points.first.at(0);
|
|
auto points_c2_up = sorted_points.first.at(1);
|
|
|
|
auto points_c1_down = sorted_points.second.at(0);
|
|
auto points_c2_down = sorted_points.second.at(1);
|
|
|
|
auto p_c1_up = points_c1_up.front();
|
|
auto p_c1_up1 = points_c1_up.back();
|
|
|
|
auto p_c2_up = points_c2_up.front();
|
|
auto p_c2_up1 = points_c2_up.back();
|
|
|
|
// if (points_c1_up.size() == 0 || points_c2_up.size() == 0 || points_c1_down.size() == 0 || points_c2_down.size() == 0)
|
|
// {
|
|
// throw std::runtime_error("Splitting the circles didn't work : a circle had no points");
|
|
// }
|
|
|
|
std::vector<std::pair<geo::Point<float>, geo::Point<float>>> pairs{
|
|
{p_c1_up,p_c2_up},
|
|
{p_c1_up,p_c2_up1},
|
|
{p_c1_up1,p_c2_up},
|
|
{p_c1_up1,p_c2_up1}
|
|
};
|
|
|
|
auto min_pair = std::min_element(std::begin(pairs), std::end(pairs),
|
|
[] (std::pair<geo::Point<float>, geo::Point<float>> const& points,
|
|
std::pair<geo::Point<float>, geo::Point<float>> const& points2) {
|
|
|
|
return (points.first.second-points.second.second).norm2() <
|
|
(points2.first.second-points2.second.second).norm2();
|
|
|
|
});
|
|
|
|
unsigned int ind_i = std::distance(std::begin(pairs), min_pair);
|
|
|
|
std::pair<geo::Point<float>, geo::Point<float>> next_points;
|
|
|
|
unsigned int ind_p_c1_up = 0, ind_p_c2_up = 0, ind_p_c2_up1 = 0, ind_p_c1_up1 = 0, ind_p_c1_down = 0, ind_p_c1_down1 = 0, ind_p_c2_down = 0, ind_p_c2_down1 = 0;
|
|
|
|
std::cout << " ind_i = " << ind_i <<std::endl;
|
|
switch(ind_i)
|
|
{
|
|
case 0 :
|
|
{
|
|
ind_p_c1_up = p_c1_up.first;
|
|
ind_p_c2_up = p_c2_up.first;
|
|
ind_p_c1_down = points_c1_down.back().first;
|
|
ind_p_c2_down = points_c2_down.back().first;
|
|
|
|
|
|
detail::pushFaceOriented(mesh,
|
|
geo::Vector3u{ind_p_c1_up, ind_p_c2_up, ind_p_c1_down},
|
|
junction_point);
|
|
detail::pushFaceOriented(mesh,
|
|
geo::Vector3u{ind_p_c1_down, ind_p_c2_up, ind_p_c2_down},
|
|
junction_point);
|
|
detail::pushFaceOriented(mesh,
|
|
geo::Vector3u{ind_p_c1_up, ind_up, ind_p_c2_up},
|
|
junction_point);
|
|
detail::pushFaceOriented(mesh,
|
|
geo::Vector3u{ind_p_c1_down, ind_p_c2_down, ind_down},
|
|
junction_point);
|
|
|
|
geo::Point<float> p_c2_down = points_c2_down.front();
|
|
next_points = std::make_pair(p_c2_up1, p_c2_down);
|
|
|
|
break;
|
|
}
|
|
case 1 :
|
|
{
|
|
ind_p_c1_up = p_c1_up.first;
|
|
ind_p_c1_down1 = points_c1_down.back().first;
|
|
ind_p_c2_up1 = p_c2_up1.first;
|
|
ind_p_c2_down = points_c2_down.front().first;
|
|
|
|
detail::pushFaceOriented(mesh,
|
|
geo::Vector3u{ind_p_c1_up, ind_p_c2_up1, ind_p_c1_down1},
|
|
junction_point);
|
|
detail::pushFaceOriented(mesh,
|
|
geo::Vector3u{ind_p_c1_down1, ind_p_c2_up1, ind_p_c2_down},
|
|
junction_point);
|
|
detail::pushFaceOriented(mesh,
|
|
geo::Vector3u{ind_p_c1_up, ind_up, ind_p_c2_up1},
|
|
junction_point);
|
|
detail::pushFaceOriented(mesh,
|
|
geo::Vector3u{ind_p_c2_down, ind_p_c1_down1, ind_down},
|
|
junction_point);
|
|
|
|
geo::Point<float> p_c2_down1 = points_c2_down.back();
|
|
next_points = std::make_pair(p_c2_up,p_c2_down1);
|
|
|
|
break;
|
|
}
|
|
case 2 :
|
|
{
|
|
|
|
ind_p_c1_up1 = p_c1_up1.first;
|
|
ind_p_c1_down = points_c1_down.front().first;
|
|
ind_p_c2_up = p_c2_up.first;
|
|
ind_p_c2_down1 = points_c2_down.back().first;
|
|
|
|
detail::pushFaceOriented(mesh,
|
|
geo::Vector3u{ind_p_c1_up1, ind_p_c2_up, ind_p_c1_down},
|
|
junction_point);
|
|
detail::pushFaceOriented(mesh,
|
|
geo::Vector3u{ind_p_c1_up1, ind_p_c1_down, ind_p_c2_down1},
|
|
junction_point);
|
|
detail::pushFaceOriented(mesh,
|
|
geo::Vector3u{ind_p_c1_up1, ind_up, ind_p_c2_up},
|
|
junction_point);
|
|
detail::pushFaceOriented(mesh,
|
|
geo::Vector3u{ind_p_c1_down, ind_p_c2_down1, ind_down},
|
|
junction_point);
|
|
|
|
geo::Point<float> p_c2_down = points_c2_down.front();
|
|
next_points = std::make_pair(p_c2_up1, p_c2_down);
|
|
|
|
break;
|
|
}
|
|
case 3 :
|
|
{
|
|
|
|
|
|
ind_p_c2_up = p_c2_up1.first;
|
|
ind_p_c1_up = p_c1_up1.first;
|
|
ind_p_c1_down = points_c1_down.front().first;
|
|
ind_p_c2_down = points_c2_down.front().first;
|
|
|
|
detail::pushFaceOriented(mesh,
|
|
geo::Vector3u{ind_p_c1_up, ind_p_c2_up, ind_p_c1_down},
|
|
junction_point);
|
|
detail::pushFaceOriented(mesh,
|
|
geo::Vector3u{ind_p_c1_down, ind_p_c2_up, ind_p_c2_down},
|
|
junction_point);
|
|
detail::pushFaceOriented(mesh,
|
|
geo::Vector3u{ind_p_c1_up, ind_up, ind_p_c2_up},
|
|
junction_point);
|
|
detail::pushFaceOriented(mesh,
|
|
geo::Vector3u{ind_p_c1_down, ind_p_c2_down, ind_down},
|
|
junction_point);
|
|
|
|
auto p_c2_down1 = points_c2_down.back();
|
|
next_points = std::make_pair(p_c2_up,p_c2_down1);
|
|
|
|
|
|
break;
|
|
}
|
|
|
|
}
|
|
|
|
for (unsigned int i = 1; i < sorted_points.first.size(); i++)
|
|
{
|
|
auto points_c3_up = sorted_points.first[(i+1)%sorted_points.first.size()];
|
|
auto points_c3_down = sorted_points.second[(i+1)%sorted_points.first.size()];
|
|
|
|
auto p_c3_up = points_c3_up.front();
|
|
auto p_c3_up1 = points_c3_up.back();
|
|
|
|
auto p_c3_down = points_c3_down.front();
|
|
auto p_c3_down1 = points_c3_down.back(); ind_p_c1_up1 = p_c1_up1.first;
|
|
ind_p_c1_down = points_c1_down.back().first;
|
|
|
|
|
|
auto p_c2_up = next_points.first;
|
|
auto p_c2_down = next_points.second;
|
|
|
|
float dist_1 = (p_c2_up.second - p_c3_up.second).norm2();
|
|
float dist_2 = (p_c2_up.second - p_c3_up1.second).norm2();
|
|
|
|
if (dist_1<dist_2)
|
|
{
|
|
unsigned int ind_p_c2_up = p_c2_up.first;
|
|
unsigned int ind_p_c2_down = p_c2_down.first;
|
|
|
|
unsigned int ind_p_c3_up = p_c3_up.first;
|
|
unsigned int ind_p_c3_down = p_c3_down1.first;
|
|
|
|
detail::pushFaceOriented(mesh,
|
|
geo::Vector3u{ind_p_c2_up, ind_p_c3_up, ind_p_c2_down},
|
|
junction_point);
|
|
detail::pushFaceOriented(mesh,
|
|
geo::Vector3u{ind_p_c2_down, ind_p_c3_up, ind_p_c3_down},
|
|
junction_point);
|
|
detail::pushFaceOriented(mesh,
|
|
geo::Vector3u{ind_p_c2_up, ind_up, ind_p_c3_up},
|
|
junction_point);
|
|
detail::pushFaceOriented(mesh,
|
|
geo::Vector3u{ind_p_c2_down, ind_p_c3_down, ind_down},
|
|
junction_point);
|
|
|
|
next_points.first = p_c3_up1;
|
|
next_points.second = p_c3_down;
|
|
|
|
}
|
|
|
|
else
|
|
{
|
|
unsigned int ind_p_c2_up = p_c2_up.first;
|
|
unsigned int ind_p_c2_down = p_c2_down.first;
|
|
|
|
unsigned int ind_p_c3_up = p_c3_up1.first;
|
|
unsigned int ind_p_c3_down = p_c3_down.first;
|
|
|
|
detail::pushFaceOriented(mesh,
|
|
geo::Vector3u{ind_p_c2_up, ind_p_c3_up, ind_p_c2_down},
|
|
junction_point);
|
|
detail::pushFaceOriented(mesh,
|
|
geo::Vector3u{ind_p_c2_down, ind_p_c3_up, ind_p_c3_down},
|
|
junction_point);
|
|
detail::pushFaceOriented(mesh,
|
|
geo::Vector3u{ind_p_c2_up, ind_up, ind_p_c3_up},
|
|
junction_point);
|
|
detail::pushFaceOriented(mesh,
|
|
geo::Vector3u{ind_p_c2_down, ind_p_c3_down, ind_down},
|
|
junction_point);
|
|
|
|
next_points.first = p_c3_up;
|
|
next_points.second = p_c3_down1;
|
|
|
|
}
|
|
|
|
|
|
|
|
}
|
|
|
|
return plane;
|
|
|
|
}
|
|
|
|
std::vector<std::pair<unsigned int, float>> Skeleton3D::findExtremities() const
|
|
{
|
|
std::vector<std::pair<unsigned int, float>> extremities;
|
|
std::vector<unsigned int> indexJunctions;
|
|
|
|
//construction of a vector containing indexes of a spline *2 + t
|
|
if (junctions.empty())
|
|
{
|
|
for (unsigned int i=0;i<splines.size();i++)
|
|
{
|
|
extremities.push_back(std::make_pair(i,0));
|
|
extremities.push_back(std::make_pair(i,1));
|
|
}
|
|
return extremities;
|
|
}
|
|
|
|
for (unsigned int i=0;i<junctions.size();i++)
|
|
{
|
|
for (unsigned int j=0;j<junctions[i].size();j++)
|
|
{
|
|
std::pair<unsigned int, float> extremSpline = junctions[i][j];
|
|
indexJunctions.push_back(extremSpline.first*2 + extremSpline.second);
|
|
}
|
|
}
|
|
|
|
std::sort(indexJunctions.begin(),indexJunctions.end());
|
|
|
|
unsigned int index=0;
|
|
for (unsigned int i=0;i<splines.size();i++)
|
|
{
|
|
if (i*2==indexJunctions[index])
|
|
{
|
|
index++;
|
|
}
|
|
else
|
|
{
|
|
extremities.push_back(std::make_pair(i,0));
|
|
}
|
|
if (((i*2)+1)==indexJunctions[index])
|
|
{
|
|
index++;
|
|
}
|
|
else
|
|
{
|
|
extremities.push_back(std::make_pair(i,1));
|
|
}
|
|
}
|
|
return extremities;
|
|
}
|
|
|
|
void Skeleton3D::meshExtremities(geo::Mesh& mesh, std::vector<std::vector<geo::Circle<float>>> const& splines_circles, unsigned int profondeur) const
|
|
{
|
|
std::vector<std::pair<unsigned int, float>> extremities = findExtremities();
|
|
for (unsigned int i=0;i<extremities.size();i++)
|
|
{
|
|
meshExtremity(mesh, splines_circles, profondeur, extremities[i]);
|
|
}
|
|
}
|
|
|
|
void Skeleton3D::meshExtremity(geo::Mesh& mesh, std::vector<std::vector<geo::Circle<float>>> const& splines_circles, unsigned int profondeur, std::pair<unsigned int,float> extremity) const
|
|
{
|
|
// Get the junction elements
|
|
auto tuple = splines[extremity.first](extremity.second);
|
|
auto C_t = std::get<0>(tuple);
|
|
auto r_t = std::get<1>(tuple);
|
|
|
|
auto S_prime_t = splines[extremity.first].prime(extremity.second);
|
|
auto C_prime_t = std::get<0>(S_prime_t);
|
|
|
|
C_prime_t /= C_prime_t.norm();
|
|
|
|
auto intersection = C_t;
|
|
|
|
if (extremity.second == 0)
|
|
{
|
|
intersection -= r_t * C_prime_t;
|
|
}
|
|
else
|
|
{
|
|
intersection += r_t * C_prime_t;
|
|
}
|
|
|
|
geo::Circle<float> extremCircle;
|
|
if (extremity.second == 0)
|
|
{
|
|
extremCircle = splines_circles[extremity.first][0];
|
|
}
|
|
else
|
|
{
|
|
extremCircle = splines_circles[extremity.first][splines_circles[extremity.first].size()-1];
|
|
}
|
|
|
|
|
|
unsigned int intersectionIndex = mesh.vertices.size();
|
|
|
|
if (profondeur==0)
|
|
{
|
|
mesh.vertices.push_back(intersection);
|
|
for (unsigned int i=0;i<extremCircle.points.size();i++)
|
|
{
|
|
auto ip1 = extremCircle.points[i].first;
|
|
if (i==extremCircle.points.size()-1)
|
|
{
|
|
auto ip2 = extremCircle.points[0].first;
|
|
mesh.faces.push_back(sortedIndexes(ip1,intersectionIndex,ip2,extremity.second));
|
|
}
|
|
else
|
|
{
|
|
auto ip2 = extremCircle.points[i+1].first;
|
|
mesh.faces.push_back(sortedIndexes(ip1,intersectionIndex,ip2,extremity.second));
|
|
}
|
|
}
|
|
}
|
|
else
|
|
{
|
|
subdivision(mesh, profondeur, extremCircle, intersectionIndex, intersection, C_t, r_t, extremity.second);
|
|
}
|
|
}
|
|
|
|
void Skeleton3D::subdivision(geo::Mesh& mesh, unsigned int profondeur, geo::Circle<float> extremCircle, unsigned int intersectionIndex, geo::Vector3f intersection, geo::Vector3f center, float radius, float t) const
|
|
{
|
|
std::vector<std::vector<geo::Vector3f>> arcs;
|
|
for (unsigned int i=0;i<extremCircle.points.size();i++)
|
|
{
|
|
// compute the points build by subdivision on one arc (between a point on the extrem
|
|
// circle and the extremity)
|
|
std::vector<geo::Vector3f> arc;
|
|
arc.push_back(intersection);
|
|
arcSubdivise(intersection, extremCircle.points[i].second, profondeur, arc, center, radius);
|
|
arc.push_back(extremCircle.points[i].second);
|
|
arcs.push_back(arc);
|
|
}
|
|
|
|
// Compute the mesh
|
|
|
|
//push the new vertices
|
|
unsigned int indexVertices = mesh.vertices.size();
|
|
unsigned int indexVerticesInit = indexVertices;
|
|
unsigned int nbPoints = arcs[0].size()-2;
|
|
|
|
//push the intersection
|
|
mesh.vertices.push_back(arcs[0][0]);
|
|
|
|
// push all the new vertices
|
|
for (unsigned int i=0;i<arcs.size();i++)
|
|
{
|
|
for (unsigned int j=1;j<arcs[i].size()-1;j++)
|
|
{
|
|
mesh.vertices.push_back(arcs[i][j]);
|
|
}
|
|
}
|
|
|
|
//push the faces
|
|
for (unsigned int i=0;i<arcs.size()-1;i++)
|
|
{
|
|
auto courant = arcs[i];
|
|
|
|
for (unsigned int j=0;j<courant.size()-2;j++)
|
|
{
|
|
if (j==0)
|
|
{
|
|
mesh.faces.push_back(sortedIndexes(intersectionIndex, indexVertices + 1 + nbPoints, indexVertices + 1,t));
|
|
}
|
|
else
|
|
{
|
|
mesh.faces.push_back(sortedIndexes(indexVertices + j, indexVertices + j + 1 + nbPoints, indexVertices + j + 1, t));
|
|
mesh.faces.push_back(sortedIndexes(indexVertices + j, indexVertices + j + nbPoints, indexVertices + j + 1 + nbPoints, t));
|
|
}
|
|
}
|
|
|
|
// before last point
|
|
unsigned int i1 = extremCircle.points[i].first;
|
|
unsigned int i2 = extremCircle.points[i+1].first;
|
|
|
|
mesh.faces.push_back(sortedIndexes(indexVertices + nbPoints, i2, i1, t));
|
|
mesh.faces.push_back(sortedIndexes(indexVertices + nbPoints , indexVertices + 2*nbPoints , i2, t));
|
|
|
|
// increment the indexVertices
|
|
indexVertices += nbPoints;
|
|
|
|
}
|
|
|
|
// the last mesh (last arc)
|
|
auto courant = arcs[arcs.size()-1];
|
|
|
|
for (unsigned int j=0;j<courant.size()-2;j++)
|
|
{
|
|
if (j==0)
|
|
{
|
|
mesh.faces.push_back(sortedIndexes(intersectionIndex, indexVerticesInit + 1, indexVertices + 1, t));
|
|
}
|
|
else
|
|
{
|
|
mesh.faces.push_back(sortedIndexes(indexVertices + j, indexVerticesInit + j + 1, indexVertices + j + 1, t));
|
|
mesh.faces.push_back(sortedIndexes(indexVertices + j, indexVerticesInit + j, indexVerticesInit + j + 1 , t));
|
|
}
|
|
}
|
|
|
|
// before last point
|
|
unsigned int i1 = extremCircle.points[arcs.size()-1].first;
|
|
unsigned int i2 = extremCircle.points[0].first;
|
|
|
|
mesh.faces.push_back(sortedIndexes(indexVertices + nbPoints, i2, i1, t));
|
|
mesh.faces.push_back(sortedIndexes(indexVertices + nbPoints, indexVerticesInit + nbPoints, i2, t));
|
|
}
|
|
|
|
geo::Vector3u Skeleton3D::sortedIndexes(unsigned int a, unsigned int b, unsigned int c, float t) const
|
|
{
|
|
if(t==0)
|
|
{
|
|
return geo::Vector3u {a,b,c};
|
|
}
|
|
else
|
|
{
|
|
return geo::Vector3u {a,c,b};
|
|
}
|
|
}
|
|
|
|
//pt1 the extremity
|
|
void Skeleton3D::arcSubdivise(geo::Vector3f pt1, geo::Vector3f pt2, unsigned int profondeur, std::vector<geo::Vector3f>& arc, geo::Vector3f center, float radius) const
|
|
{
|
|
if (profondeur!=0)
|
|
{
|
|
// compute middle of arc
|
|
geo::Vector3f middleSegment;
|
|
middleSegment = (pt1+pt2)/2;
|
|
// projection on the sphere
|
|
geo::Vector3f direction = (middleSegment - center);
|
|
direction /= direction.norm();
|
|
direction *= radius;
|
|
geo::Vector3f middle = center + direction;
|
|
|
|
arcSubdivise(pt1, middle, profondeur-1, arc, center, radius);
|
|
arc.push_back(middle);
|
|
arcSubdivise(middle, pt2, profondeur-1, arc, center, radius);
|
|
}
|
|
}
|
|
|
|
geo::Mesh Skeleton3D::computeMesh(float pointsProportion, unsigned int nbPointsPerCircle, unsigned int depth,bool withExtremities, bool withJunctions) const
|
|
{
|
|
geo::Mesh mesh;
|
|
std::vector<std::vector<geo::Circle<float>>> splines_circles;
|
|
unsigned int globalOffset = 0;
|
|
|
|
for (auto const& spline : splines)
|
|
{
|
|
float l = geo::splineLenght(spline);
|
|
unsigned int nbCirclesOnTheSpline = std::max(static_cast<unsigned int>(std::round(pointsProportion*l)),2u);
|
|
auto spl = spline.computeCircles(nbCirclesOnTheSpline,nbPointsPerCircle,globalOffset);
|
|
splines_circles.push_back(spl);
|
|
mesh += geo::splineMeshing(spl);
|
|
globalOffset+=nbPointsPerCircle*spl.size();
|
|
}
|
|
|
|
if (withExtremities)
|
|
meshExtremities(mesh, splines_circles, depth);
|
|
|
|
if (withJunctions)
|
|
meshJunctions(mesh, splines_circles);
|
|
|
|
return mesh;
|
|
}
|
|
|
|
std::ostream& operator<<(std::ostream& out, Skeleton3D const& s)
|
|
{
|
|
for (auto const& spline:s.splines)
|
|
{
|
|
out << spline << std::endl;
|
|
}
|
|
|
|
for (auto const& junction:s.junctions)
|
|
{
|
|
out << "l " << junction.size() << " ";
|
|
for (auto const& pair:junction)
|
|
{
|
|
out << pair.first << " ";
|
|
}
|
|
for (auto const& pair:junction)
|
|
{
|
|
out << pair.second << " ";
|
|
}
|
|
out << std::endl;
|
|
}
|
|
return out;
|
|
}
|
|
|
|
namespace detail
|
|
{
|
|
|
|
std::pair<std::vector<std::list<geo::Point<float>>>, std::vector<std::list<geo::Point<float>>>> sortedPoints(geo::Vector4f const& plane, std::vector<geo::Circle<float>> const& junction_circles)
|
|
{
|
|
std::vector<std::list<geo::Point<float>>> up_points;
|
|
std::vector<std::list<geo::Point<float>>> down_points;
|
|
|
|
for (unsigned int i = 0; i < junction_circles.size(); i++)
|
|
{
|
|
std::list<geo::Point<float>> first_points; // points located on one side of the circle, the first side reached.The two sides are separated by the median plane computed earlier.
|
|
std::list<geo::Point<float>> second_points; // points located on the other side of the circle
|
|
const bool up_side = geo::whichSide(plane, junction_circles[i].points[0].second); // constant boolean initialized at the beginning to know with which side we begin the way on the circle.
|
|
bool side_current_points = up_side; // boolean to know on which side we must had the points.
|
|
|
|
for (unsigned int j = 0; j < junction_circles[i].points.size(); j++)
|
|
{
|
|
if (up_side == side_current_points)
|
|
{
|
|
if (geo::whichSide(plane, junction_circles[i].points[j].second) == up_side)
|
|
{
|
|
first_points.push_back(junction_circles[i].points[j]);
|
|
}
|
|
else
|
|
{
|
|
side_current_points = geo::whichSide(plane, junction_circles[i].points[j].second);
|
|
second_points.push_back(junction_circles[i].points[j]);
|
|
}
|
|
}
|
|
else
|
|
{
|
|
if (geo::whichSide(plane, junction_circles[i].points[j].second) == side_current_points)
|
|
{
|
|
second_points.push_back(junction_circles[i].points[j]);
|
|
}
|
|
else
|
|
{
|
|
break;
|
|
}
|
|
}
|
|
|
|
}
|
|
|
|
for (int j = static_cast<int>(junction_circles[i].points.size()-1); j >= 0; j--)
|
|
{
|
|
if (geo::whichSide(plane, junction_circles[i].points[j].second) == up_side)
|
|
{
|
|
first_points.push_front(junction_circles[i].points[j]);
|
|
}
|
|
else
|
|
{
|
|
break;
|
|
}
|
|
}
|
|
|
|
if (up_side)
|
|
{
|
|
up_points.push_back(first_points);
|
|
down_points.push_back(second_points);
|
|
}
|
|
else
|
|
{
|
|
up_points.push_back(second_points);
|
|
down_points.push_back(first_points);
|
|
|
|
}
|
|
}
|
|
|
|
return std::make_pair(up_points, down_points);
|
|
}
|
|
|
|
void pushFaceOriented(geo::Mesh& mesh, geo::Vector3u const& face, geo::Vector3f const& sphere_center)
|
|
{
|
|
geo::Vector3f v1 = mesh.vertices[face.x()] - sphere_center;
|
|
geo::Vector3f v2 = geo::crossProduct(mesh.vertices[face.z()] - mesh.vertices[face.y()],
|
|
mesh.vertices[face.z()] - mesh.vertices[face.x()]);
|
|
|
|
if (v1*v2 < 0)
|
|
{
|
|
mesh.faces.push_back(face);
|
|
}
|
|
else
|
|
{
|
|
mesh.faces.push_back(geo::Vector3u{face.x(),face.z(),face.y()});
|
|
}
|
|
}
|
|
|
|
} // namespace detail
|
|
|
|
} // namespace pae
|