We are hiring ! See our job offers.
Revision 7d1d9c1c9d76fee0d09beff2396006f8754e8613 authored by Gert Meijer on 06 May 2021, 03:23:18 UTC, committed by Gert Meijer on 06 May 2021, 03:23:18 UTC
1 parent 328a518
prickliness.cpp
``````#include "precomp.h"
#include "prickliness.h"

#include "DrawPrickliness.h"

namespace Prickliness
{
//Get a vector to a point inside a convex DCEL plane
Vector_3 PricklinessArrangement::GetVectorInPlane(const Arrangement_2::Face_const_handle face) const
{
std::vector<Point_2> edgePoints;

if (face->has_outer_ccb())
{
Arrangement_2::Ccb_halfedge_const_circulator currentHalfEdge = face->outer_ccb();

//Select a point on 2 of the (non-fictitious) edges to calculate a point inside the face
do
{
//Fictitious: https://doc.cgal.org/latest/Arrangement_on_surface_2/index.html#title32
if (!currentHalfEdge->is_fictitious())
{
if (currentHalfEdge->curve().is_ray())
{
//Edge is an infinite ray originating from a vertex on an unbounded face, push an arbitrary point on the ray
auto rayStart = currentHalfEdge->curve().ray().start();
auto rayVector = currentHalfEdge->curve().ray().direction().to_vector();

edgePoints.push_back(rayStart + (rayVector * Kernel::FT(0.5)));
}
else if (currentHalfEdge->curve().is_segment())
{
//Edge is a line segment bounded by two vertices, push the point in the middle of the line segment
auto heSegment = currentHalfEdge->curve().segment();
auto segmentVector = heSegment.end() - heSegment.start();

edgePoints.push_back(heSegment.start() + (segmentVector * Kernel::FT(0.5)));
}
}

} while (++currentHalfEdge != face->outer_ccb() && edgePoints.size() < 2);
}

//Get the point in between the two edge points and convert it to 3D
//the arrangement plane is at z = 10
Point_2 startFacePoint = edgePoints.at(0) + ((edgePoints.at(1) - edgePoints.at(0)) * Kernel::FT(0.5));
Vector_3 faceVector(startFacePoint.x(), startFacePoint.y(), 10.0);

return faceVector;
}

bool PricklinessArrangement::IsBoundaryVert(const Delaunay& DT, Delaunay::Face_circulator adjacentFaces) const
{

do
{
{
return true;
}

return false;
}

PricklinessArrangement::PricklinessArrangement(const Delaunay& DT)
{
Plane_3 arrangementPlane(Point_3(0.0, 0.0, 10.0), Vector_3(0.0, 0.0, 1.0));

std::vector<Arr_X_monotone_curve_2> arrangementLines;
std::vector<std::vector<Plane_3>> perpPlanes;
std::vector<int> vertexDegrees;
int currentVertIndex = 0;

std::cout << "\t Constructing perpendicular planes" << std::endl;

size_t DTVertCount = DT.number_of_vertices();
//size_t percent = DTVertCount / 100;

//Loop through all the finite vertices of the DT, the underlying structures adds a point at infinity to connect all the boundary verts together

for (Del_finite_vertices_iterator vertIterator = DT.finite_vertices_begin(); vertIterator != DT.finite_vertices_end(); vertIterator++)
{
if (IsBoundaryVert(DT, vertIterator->incident_faces()))
{
continue;
}

//Add a perp plane vector for this vertex
perpPlanes.push_back(std::vector<Plane_3>());

//Store current vertex index for the result plane traversal
PerpPlaneIndices vnpTriple;
vnpTriple.vertexIndex = currentVertIndex;

//Get next vertex and its neighbours iterator
Point_3 currentVertex = vertIterator->point();
Del_vertex_circulator neighbourCirculator = vertIterator->incident_vertices(), done(neighbourCirculator);

int neighbourIndex = 0;
if (neighbourCirculator != 0)
{
do
{
//Skip the vertex at infinity that the triangulation DS uses to connect the boundary verts
if (!DT.is_infinite(neighbourCirculator))
{
//Get neighbour vertex and subtract current current vertex to move it towards the origin
Point_3 neighbour = neighbourCirculator->point();
Vector_3 neighbourToOrigin = neighbour - currentVertex;

//Construct a line going through the neighbour and current vector at the origin
//Flip the direction so it goes towards the current vector
//(it's important to make sure it's oriented this way because we need to check if we are on the local maxima side of the plane with a dot on the normal)
Line_3 lineThroughOrigin = Line_3(CGAL::ORIGIN, neighbourToOrigin).opposite();

//Store the index of the neighbour so we can use it to get it's perp plane
vnpTriple.planeIndex = perpPlanes.at(currentVertIndex).size();

//Construct a plane perpendiculat to the neighbour line at the origin
Plane_3 perpPlane = lineThroughOrigin.perpendicular_plane(CGAL::ORIGIN);
perpPlanes.at(currentVertIndex).push_back(perpPlane);

//Intersect the perp plane with the arrangement plane, store the resulting 2D line so we can mass add them in the arrangement
auto intersection = CGAL::intersection(arrangementPlane, perpPlane);

if (intersection)
{
if (Line_3* arrangementLine = boost::get<Line_3>(&*intersection))
{
Point_3 p3 = arrangementLine->point();
Vector_3 v3 = arrangementLine->to_vector();

Arr_data_line_2 arrangementLine2(Point_2(p3.x(), p3.y()), Vector_2(v3.x(), v3.y()));

arrangementLines.push_back(Arr_X_monotone_curve_2(arrangementLine2, vnpTriple));

//CGAL::insert(resultArrangement, Arr_X_monotone_curve_2(arrangementLine2, vnpTriple));
}
else
{
//This will never happen because the perp planes go through the origin and the arrangement plane lies above it.
std::cout << "Result of intersection between perpendicular plane and arrangement plane is a plane!" << std::endl;
}
}
else
{
//Parallel planes degenerate case: Ignore? Should never happen unless two vertices are right on top of each other
std::cout << "Empty intersection between perpendicular plane and arrangement plane!" << std::endl;
}

neighbourIndex++;
}
} while (++neighbourCirculator != done);

vertexDegrees.push_back(neighbourIndex); //Save degree this way because vert->degree() counts fictitious neighbours..
}

currentVertIndex++;

//if (currentVertIndex % percent == 0)
//{
//    std::cout << currentVertIndex / percent << std::endl;
//}

}

std::cout << "\t Constructing 2d arrangement (DCEL)" << std::endl;
//Construct an arrangement (DCEL) that keeps track of the original lines the half-edges are created from
//Arrangement_2 resultArrangement;
CGAL::insert(resultArrangement, arrangementLines.begin(), arrangementLines.end());
arrangementLines.clear();
//FIFO queues for breadth first traversal of the arrangement
std::queue<Arrangement_2::Face_handle> faceQueue;
std::queue<std::vector<int>> abovePlaneCountsQueue;

std::cout << "\t Determining local maxima for DCEL faces" << std::endl;
if (resultArrangement.number_of_faces() > 1)
{
//Calculate the local maxima for the (random) starting face
auto startingFace = resultArrangement.faces_begin();

//Get a vector to a point in the face so we can check if we are above the perp planes
Vector_3 startFaceVector = GetVectorInPlane(startingFace);

//Construct the data structure keeping track of planes and the local maxima
std::vector<int> abovePlanesCount;
abovePlanesCount.reserve(perpPlanes.size());

for (std::vector<Plane_3>& planeVector : perpPlanes)
{
int aboveCount = 0;
for (Plane_3& plane : planeVector)
{
//Count as above if perpendicular plane normal dot face > 0
if (CGAL::scalar_product(plane.orthogonal_vector(), startFaceVector) > Kernel::FT(0))
{
aboveCount++;
}
}
abovePlanesCount.push_back(aboveCount);
}

//Count for how many terrain vertices we are in the local maxima cone, this is the # local maxima for this arrangement face
//We are in the cone if we are above all of its planes
int localMaximaCount = 0;
for (size_t i = 0; i < abovePlanesCount.size(); i++)
{
if (abovePlanesCount.at(i) == vertexDegrees.at(i))
{
localMaximaCount++;
}
}

startingFace->data().localMaximaCount = localMaximaCount;

//Put the starting face and the counters for the planes on the queues for the breadth first iteration
faceQueue.push(startingFace);
abovePlaneCountsQueue.push(abovePlanesCount);

//Breadth first search over the terrain:
//Everytime we cross an edge check its corresponding perp plane(s) with a vector to the point in the next face:
//if we are above a plane increment its corresponding above counter, decrement otherwise.
//if the above counter is equal to the degree of the vertex after an increment we are within its local maxima cone and increment # local maxima
//if the above counter is equal to the degree of the vertex before an decrement we have left the local maxima cone and decrement # local maxima
while (!faceQueue.empty())
{
auto currentFace = faceQueue.front();
faceQueue.pop();

std::vector<int> currentAbovePlaneCounts = abovePlaneCountsQueue.front();
abovePlaneCountsQueue.pop();

//Iterate through the incident edges
Arrangement_2::Ccb_halfedge_circulator currentHalfEdge = currentFace->outer_ccb();
do
{
//Fictitious: https://doc.cgal.org/latest/Arrangement_on_surface_2/index.html#title32
if (!currentHalfEdge->is_fictitious())
{
//Skip neighbour face if already visited
if (currentHalfEdge->twin()->face()->data().localMaximaCount != -1)
{
continue;
}

std::vector<int> neighbourAbovePlanesCounts = currentAbovePlaneCounts;
int localMaximaCount = currentFace->data().localMaximaCount;

Vector_3 faceVector = GetVectorInPlane(currentHalfEdge->twin()->face());

//This is a for loop because the DCEL stores the data in a list in case of the degenerate case where multiple lines overlap
for (PerpPlaneIndices ppIndices : currentHalfEdge->curve().data())
{
Plane_3 perpPlane = perpPlanes.at(ppIndices.vertexIndex).at(ppIndices.planeIndex);

//Which side of the perp plane did we go to?
if (CGAL::scalar_product(perpPlane.orthogonal_vector(), faceVector) > Kernel::FT(0))
{
//Above, increment above counter and check if inside local maxima cone
neighbourAbovePlanesCounts.at(ppIndices.vertexIndex)++;

if (neighbourAbovePlanesCounts.at(ppIndices.vertexIndex) == vertexDegrees.at(ppIndices.vertexIndex))
{
localMaximaCount++;
}
}
else
{
//Below, check if we were inside the local maxima cone and decrement above counter
if (neighbourAbovePlanesCounts.at(ppIndices.vertexIndex) == vertexDegrees.at(ppIndices.vertexIndex))
{
localMaximaCount--;
}
neighbourAbovePlanesCounts.at(ppIndices.vertexIndex)--;
}
}

currentHalfEdge->twin()->face()->data().localMaximaCount = localMaximaCount;

//Local maxima count for this face current most?
if (prickliness < localMaximaCount)
{
prickliness = localMaximaCount;
pricklyDirection = faceVector;
}

faceQueue.push(currentHalfEdge->twin()->face());
abovePlaneCountsQueue.push(neighbourAbovePlanesCounts);
}
} while (++currentHalfEdge != currentFace->outer_ccb());
}
}

//CGAL::draw(resultArrangement);
}

void PricklinessArrangement::Draw()
{
Arrangement_2::Vertex_const_iterator vit = resultArrangement.vertices_begin();

Kernel::FT minX = vit->point().x();
Kernel::FT minY = vit->point().y();
Kernel::FT maxX = vit->point().x();
Kernel::FT maxY = vit->point().y();

for (; vit != resultArrangement.vertices_end(); ++vit) {

if (vit->point().x() < minX) { minX = vit->point().x(); }
if (vit->point().y() < minY) { minY = vit->point().y(); }
if (vit->point().x() > maxX) { maxX = vit->point().x(); }
if (vit->point().y() > maxY) { maxY = vit->point().y(); }
}

Point_2 minPoint(minX, minY);
Point_2 maxPoint(maxX, maxY);

minPoint -= Vector_2(10, 10);
maxPoint += Vector_2(10, 10);
std::cout << "Min:" << minPoint << " Max:" << maxPoint << std::endl;

Point_2 topLeft(minPoint.x(), maxPoint.y());
Point_2 bottomRight(maxPoint.x(), minPoint.y());

std::vector<Arr_traits_2::Curve_2> boundarybox;
boundarybox.push_back(Arr_traits_2::Curve_2(topLeft, minPoint));
boundarybox.push_back(Arr_traits_2::Curve_2(minPoint, bottomRight));
boundarybox.push_back(Arr_traits_2::Curve_2(bottomRight, maxPoint));
boundarybox.push_back(Arr_traits_2::Curve_2(maxPoint, topLeft));

CGAL::insert(resultArrangement, boundarybox.begin(), boundarybox.end());

std::cout << "Drawing Prickliness (" << prickliness << "): " << resultArrangement.number_of_faces() << " Faces, " << resultArrangement.number_of_edges() << " Edges, " << resultArrangement.number_of_vertices() << " Vertices." << std::endl;

DrawPrickliness(resultArrangement, prickliness, pricklyDirection);
}

void PricklinessArrangement::WriteToIpe(const std::filesystem::path& filePath, bool centerOnOrigin)
{
//Merge faces
Arrangement_2::Edge_iterator edgeIt;
std::vector<Arrangement_2::Halfedge_handle> redundantEdges;
std::vector<int> pricklyValues;
for (edgeIt = resultArrangement.edges_begin(); edgeIt != resultArrangement.edges_end(); ++edgeIt)
{
Arrangement_2::Halfedge_handle hEdge = edgeIt;

if (hEdge->face()->data().localMaximaCount == hEdge->twin()->face()->data().localMaximaCount)
{
int prickly = hEdge->face()->data().localMaximaCount;
redundantEdges.push_back(hEdge);
pricklyValues.push_back(hEdge->face()->data().localMaximaCount);
}
}

int currIndex = 0;
for (Arrangement_2::Halfedge_handle hEdge : redundantEdges)
{
Arrangement_2::Face_handle fh = resultArrangement.remove_edge(hEdge, true, true);
fh->data().localMaximaCount = pricklyValues.at(currIndex);
++currIndex;
}

//Remove redundant vertices
std::vector<Arrangement_2::Vertex_handle> vertexHandles;
for (Arrangement_2::Vertex_handle vertexHandle : resultArrangement.vertex_handles())
{
vertexHandles.push_back(vertexHandle);
}
for (Arrangement_2::Vertex_handle vertexHandle : vertexHandles)
{
CGAL::remove_vertex(resultArrangement, vertexHandle);
}
vertexHandles.clear();

#if 0  //Boundary box variable size
Arrangement_2::Vertex_const_iterator vit = resultArrangement.vertices_begin();

Kernel::FT minX = vit->point().x();
Kernel::FT minY = vit->point().y();
Kernel::FT maxX = vit->point().x();
Kernel::FT maxY = vit->point().y();

for (; vit != resultArrangement.vertices_end(); ++vit) {

if (vit->point().x() < minX) { minX = vit->point().x(); }
if (vit->point().y() < minY) { minY = vit->point().y(); }
if (vit->point().x() > maxX) { maxX = vit->point().x(); }
if (vit->point().y() > maxY) { maxY = vit->point().y(); }
}

Point_2 minPoint(minX, minY);
Point_2 maxPoint(maxX, maxY);

minPoint -= Vector_2(1, 1);
maxPoint += Vector_2(1, 1);
#else //Static size

Point_2 minPoint;
Point_2 maxPoint;

if (centerOnOrigin)
{
//square centered on prickly direction
minPoint = CGAL::ORIGIN;
maxPoint = CGAL::ORIGIN;
}
else
{
//square centered on prickly direction
minPoint = Point_2(pricklyDirection.x(), pricklyDirection.y());
maxPoint = Point_2(pricklyDirection.x(), pricklyDirection.y());
}

minPoint -= Vector_2(4, 4);
maxPoint += Vector_2(4, 4);

#endif // boundary box

std::cout << "Min:" << minPoint << " Max:" << maxPoint << std::endl;

Point_2 topLeft(minPoint.x(), maxPoint.y());
Point_2 bottomRight(maxPoint.x(), minPoint.y());

std::vector<Arr_traits_2::Curve_2> boundarybox;
boundarybox.push_back(Arr_traits_2::Curve_2(topLeft, minPoint));
boundarybox.push_back(Arr_traits_2::Curve_2(minPoint, bottomRight));
boundarybox.push_back(Arr_traits_2::Curve_2(bottomRight, maxPoint));
boundarybox.push_back(Arr_traits_2::Curve_2(maxPoint, topLeft));

CGAL::insert(resultArrangement, boundarybox.begin(), boundarybox.end());

const std::string ipeFileStart = "<?xml version=\"1.0\"?>\n"
"<!DOCTYPE ipe SYSTEM \"ipe.dtd\">\n"
"<ipe version = \"70206\" creator=\"Ipe 7.2.7\">\n"
"<page>\n"
"<layer name = \"alpha\"/>\n"
"<view layers = \"alpha\" active=\"alpha\"/>\n";

const std::string ipeFileEnd = "</page>\n</ipe>";

Arrangement_2::Face_const_iterator fit;

//Find lowest and highest values in current view (This double loop stuff is kind of ugly but it's only for some basic i/o so whatever)
int currentViewMaximaMax = 0;
int currentViewMaximaMin = prickliness;

for (fit = resultArrangement.faces_begin(); fit != resultArrangement.faces_end(); ++fit)
{
//Only export finite faces
if (!fit->is_unbounded() && !fit->is_fictitious() && fit->has_outer_ccb())
{
if (fit->data().localMaximaCount != -1)
{
auto edgeCirculator = fit->outer_ccb();
auto startEdgeCirc = fit->outer_ccb();

bool inside = true;

do
{
if (edgeCirculator->source()->point().x() < minPoint.x() ||
edgeCirculator->source()->point().x() > maxPoint.x() ||
edgeCirculator->source()->point().y() < minPoint.y() ||
edgeCirculator->source()->point().y() > maxPoint.y())
{
inside = false;
break;
}

} while (++edgeCirculator != startEdgeCirc);

//Skip faces outside of bounding box
if (!inside)
{
continue;
}

if (currentViewMaximaMax < fit->data().localMaximaCount)
{
currentViewMaximaMax = fit->data().localMaximaCount;
}

if (currentViewMaximaMin > fit->data().localMaximaCount)
{
currentViewMaximaMin = fit->data().localMaximaCount;
}
}
}
}

std::cout << "Min prickly in view: " << currentViewMaximaMin << std::endl;
std::cout << "Max prickly in view: " << currentViewMaximaMax << std::endl;

std::ostringstream fileName;
fileName << "prickliness_";

if (centerOnOrigin)
{
fileName << "O";
}
else
{
fileName << "P";
}

fileName << "_" << currentViewMaximaMin << "_to_" << currentViewMaximaMax << ".ipe";

std::ofstream ipeFileStream(filePath / fileName.str());

if (!ipeFileStream.is_open())
{
std::cout << "\t Ipe file path invalid or not openable.." << std::endl;
return;
}

ipeFileStream << ipeFileStart;

//Export the arrangement faces.
for (fit = resultArrangement.faces_begin(); fit != resultArrangement.faces_end(); ++fit)
{
//Only export finite faces
if (!fit->is_unbounded() && !fit->is_fictitious() && fit->has_outer_ccb())
{

auto edgeCirculator = fit->outer_ccb();
auto startEdgeCirc = fit->outer_ccb();

bool inside = true;

do
{
if (edgeCirculator->source()->point().x() < minPoint.x() ||
edgeCirculator->source()->point().x() > maxPoint.x() ||
edgeCirculator->source()->point().y() < minPoint.y() ||
edgeCirculator->source()->point().y() > maxPoint.y())
{
inside = false;
break;
}

} while (++edgeCirculator != startEdgeCirc);

//Skip faces outside of bounding box
if (!inside)
{
continue;
}

//Determine face color
if (fit->data().localMaximaCount != -1)
{
//double shade = (((double)fit->data().localMaximaCount / (double)prickliness));
//double shade = (((double)(fit->data().localMaximaCount - currentViewMaximaMin) / (double)(currentViewMaximaMax - currentViewMaximaMin)));

////Green -> Red (https://stackoverflow.com/a/7947812)
//double r = std::fmin(2.0 * shade, 1.0);
//double g = std::fmin(2.0 * (1.0 - shade), 1.0);
//double b = 0.0;

////Yellow -> Red
//double r = 1.0;
//double g = 1.0 - shade;
//double b = 0;

//PricklyColor color = GetColorForPercentage((double)(fit->data().localMaximaCount - currentViewMaximaMin) / (double)(currentViewMaximaMax - currentViewMaximaMin));
PricklyColor color = GetColorForPercentage((double)(fit->data().localMaximaCount) / (double)(currentViewMaximaMax));

//Polygon opening tag, set both stroke and fill to pickly shade so we don't get the wierd lines in between the faces
ipeFileStream << "<path stroke = \"" << color.r << " " << color.g << " " << color.b << "\" fill=\"" << color.r << " " << color.g << " " << color.b << "\">\n";
}
else
{
//White face for prickliless (I just made that word)faces
ipeFileStream << "<path stroke = \"1 1 1\" fill=\"1 1 1\">\n";
}

edgeCirculator = fit->outer_ccb();
startEdgeCirc = fit->outer_ccb();

//Circle along the edges and add the vertices to the ipe polygon
bool first = true;
do
{
ipeFileStream << (edgeCirculator->source()->point().x() * Kernel::FT(10)) << " " << (edgeCirculator->source()->point().y() * Kernel::FT(10));

if (first)
{
ipeFileStream << " m\n";
first = false;
}
else
{
ipeFileStream << " l\n";
}
} while (++edgeCirculator != startEdgeCirc);

ipeFileStream << "h\n";

//Polygon closing tag
ipeFileStream << "</path>\n";

}
}

ipeFileStream << ipeFileEnd;

}

//Thank you https://stackoverflow.com/a/7128796
PricklinessArrangement::PricklyColor PricklinessArrangement::GetColorForPercentage(double percent) const
{
static std::vector<std::pair<double, PricklyColor>> percentColors
{
{ 0.0,   { 1.0,	        1.0,	     1.0 }},
{ 0.125, { 0.929411765,	0.97254902,	 0.694117647 }},
{ 0.25,	 { 0.780392157,	0.91372549,	 0.705882353}},
{ 0.375, { 0.498039216,	0.803921569, 0.733333333}},
{ 0.5,	 { 0.254901961,	0.71372549,	 0.768627451}},
{ 0.625, { 0.11372549,	0.568627451, 0.752941176}},
{ 0.75,	 { 0.133333333,	0.368627451, 0.658823529}},
{ 0.875, { 0.145098039,	0.203921569, 0.580392157}},
{ 1.0,	 { 0.031372549,	0.11372549,	 0.345098039}}
};

int i = 1;
for (; i < percentColors.size() - 1; i++)
{
if (percent < percentColors.at(i).first)
{
break;
}
}

std::pair<double, PricklyColor> lower = percentColors.at(i - 1);
std::pair<double, PricklyColor> upper = percentColors.at(i);
double range = upper.first - lower.first;
double rangePercent = (percent - lower.first) / range;
double percentLower = 1 - rangePercent;
double percentUpper = rangePercent;

PricklyColor color;
color.r = lower.second.r * percentLower + upper.second.r * percentUpper;
color.g = lower.second.g * percentLower + upper.second.g * percentUpper;
color.b = lower.second.b * percentLower + upper.second.b * percentUpper;

return color;
}
}``````

Computing file changes ...