Commit 7bdfdbb5 authored by Sreekanth Arikatla's avatar Sreekanth Arikatla

Merge branch 'OctreeDataStructure' into 'master'

Octree Data Structure

See merge request iMSTK/iMSTK!342
parents d3f7e50b 5e964bcb
......@@ -142,9 +142,9 @@ main()
}
++mode;
debugPoints->turnDataModifiedFlagON();
debugLines->turnDataModifiedFlagON();
debugTriangles->turnDataModifiedFlagON();
debugPoints->setDataModified(true);
debugLines->setDataModified(true);
debugTriangles->setDataModified(true);
// Must call to update render data
for (auto& delegate : renderer->getDebugRenderDelegates())
......
This diff is collapsed.
###########################################################################
#
# Copyright (c) Kitware, Inc.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0.txt
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#
###########################################################################
project(Example-Octree)
#-----------------------------------------------------------------------------
# Create executable
#-----------------------------------------------------------------------------
add_executable(${PROJECT_NAME} OctreeExample.cpp BunnyMesh.cpp)
#-----------------------------------------------------------------------------
# Add shaders
#-----------------------------------------------------------------------------
include(imstkCopyAndCompileShaders)
CopyAndCompileShaders()
#-----------------------------------------------------------------------------
# Link libraries to executable
#-----------------------------------------------------------------------------
target_link_libraries(${PROJECT_NAME} SimulationManager)
/*=========================================================================
Library: iMSTK
Copyright (c) Kitware, Inc. & Center for Modeling, Simulation,
& Imaging in Medicine, Rensselaer Polytechnic Institute.
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0.txt
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
=========================================================================*/
#include "imstkSimulationManager.h"
#include "imstkSceneObject.h"
#include "imstkDebugGeometry.h"
#include "imstkAPIUtilities.h"
#include "imstkLooseOctree.h"
#include "imstkVTKViewer.h"
#include "imstkTimer.h"
#include "imstkVTKTextStatusManager.h"
#include <thread>
#include <chrono>
#include <sstream>
using namespace imstk;
//#define NUM_MESHES 4u
#define NUM_MESHES 10u
// Load bunny mesh data (vertex positions and triangle faces)
std::pair<StdVectorOfVec3d, std::vector<std::array<size_t, 3>>> getBunny();
static std::pair<StdVectorOfVec3d, std::vector<std::array<size_t, 3>>> g_BunnyData = getBunny();
///
/// \brief Read a mesh, create a visual scene object and add to the scene
///
std::shared_ptr<VisualObject>
createMeshObject(const std::shared_ptr<imstk::Scene>& scene,
const std::string& objectName,
Color color)
{
// Create a surface mesh for the bunny
auto meshObj = std::make_shared<SurfaceMesh>();
meshObj->initialize(g_BunnyData.first, g_BunnyData.second);
// Create a visiual model
auto visualModel = std::make_shared<VisualModel>(meshObj);
auto material = std::make_shared<RenderMaterial>();
material->setDebugColor(color); // Wireframe color
material->setLineWidth(4);
material->setDisplayMode(RenderMaterial::DisplayMode::WIREFRAME_SURFACE);
visualModel->setRenderMaterial(material);
auto visualObject = std::make_shared<VisualObject>(objectName);
visualObject->addVisualModel(visualModel);
scene->addSceneObject(visualObject);
return visualObject;
}
///
/// \brief Generate a random color
///
Color
getRandomColor()
{
Color color(0, 0, 0, 1);
while (true)
{
for (unsigned int i = 0; i < 3; ++i)
{
color.rgba[i] = static_cast<double>(rand()) / static_cast<double>(RAND_MAX);
}
if (color.rgba[0] > 0.95
|| color.rgba[1] > 0.95
|| color.rgba[2] > 0.95)
{
break;
}
}
return color;
}
///
/// \brief This example demonstrates the usage of octree
///
int
main()
{
// SDK and Scene
auto sdk = std::make_shared<SimulationManager>();
auto scene = sdk->createNewScene("Octree Example");
sdk->setActiveScene(scene);
// Get the VTKViewer
auto viewer = std::dynamic_pointer_cast<VTKViewer>(sdk->getViewer());
viewer->getVtkRenderWindow()->SetSize(1920, 1080);
auto statusManager = viewer->getTextStatusManager();
statusManager->setStatusFontSize(VTKTextStatusManager::Custom, 30);
statusManager->setStatusFontColor(VTKTextStatusManager::Custom, Color::Orange);
// Get VTK Renderer
auto renderer = std::dynamic_pointer_cast<VTKRenderer>(viewer->getActiveRenderer());
LOG_IF(FATAL, (!renderer)) << "Invalid renderer: Only VTKRenderer is supported for debug rendering";
// srand(123456); // Deterministic random generation, for random colors
srand(static_cast<unsigned int>(time(nullptr)));
// Create meshes
std::vector<std::shared_ptr<SurfaceMesh>> triMeshes;
for (unsigned int i = 0; i < NUM_MESHES; ++i)
{
const auto sceneObj = createMeshObject(scene, "Mesh-" + std::to_string(triMeshes.size()), getRandomColor());
triMeshes.push_back(std::dynamic_pointer_cast<SurfaceMesh>(sceneObj->getVisualGeometry()));
}
// Compute the scale factor to scale meshes such that meshes with different sizes are still visualized consistently
Vec3d lowerCorner, upperCorner;
const auto pointset = std::dynamic_pointer_cast<PointSet>(triMeshes.front());
ParallelUtils::ParallelReduce::findAABB(pointset->getVertexPositions(), lowerCorner, upperCorner);
const auto scaleFactor = 20.0 / (upperCorner - lowerCorner).norm();
for (const auto& mesh: triMeshes)
{
mesh->scale(scaleFactor, Geometry::TransformType::ApplyToData);
}
StopWatch timer;
timer.start();
// Create octree
LooseOctree octree(Vec3d(0, 0, 0), 100.0, 0.125, 2.0, "TestOctree");
// Add all meshes to the octree
for (const auto& mesh: triMeshes)
{
octree.addTriangleMesh(mesh);
}
// Build octree after adding all geometries
octree.build();
LOG(INFO) << "Build octree time: " << timer.getTimeElapsed() << " ms";
// Always rebuild tree from scratch in each update (default update is incremental update)
// This is significantly slower than incremental update!
// octree.setAlwaysRebuild(true);
// Create debug geometry for the octree (render up to 8 levels, and render all non-emtpy nodes)
const auto debugOctree = octree.getDebugGeometry(8, true);
scene->addDebugGeometry(debugOctree);
// Data for animation
const double translation = 15.0;
StdVectorOfVec3d centers;
StdVectorOfVec3d dirs;
for (unsigned int i = 0; i < NUM_MESHES; ++i)
{
centers.push_back(Vec3d(translation, 0, 0));
dirs.push_back(Vec3d(-1, 0, 0));
}
// Transform the mesh objects
const double angle = 2.0 * PI / NUM_MESHES;
for (unsigned int i = 0; i < NUM_MESHES; ++i)
{
const auto rotation = angle * static_cast<double>(i);
triMeshes[i]->translate(translation, 0, 1, Geometry::TransformType::ApplyToData);
triMeshes[i]->rotate(Vec3d(0, 1, 0), rotation, Geometry::TransformType::ApplyToData);
auto t = centers[i][0];
centers[i][0] = std::cos(rotation) * t;
centers[i][2] = -std::sin(rotation) * t;
t = dirs[i][0];
dirs[i][0] = std::cos(rotation) * t;
dirs[i][2] = -std::sin(rotation) * t;
}
auto updateFunc =
[&](Module*) {
// Move objects
for (size_t i = 0; i < triMeshes.size(); ++i)
{
triMeshes[i]->translate(dirs[i][0], dirs[i][1], dirs[i][2], Geometry::TransformType::ApplyToData);
centers[i] += dirs[i];
}
Vec3d lowerCorners, upperCorner;
ParallelUtils::ParallelReduce::findAABB(centers, lowerCorners, upperCorner);
if ((lowerCorners - upperCorner).norm() > 70.0)
{
for (size_t i = 0; i < dirs.size(); ++i)
{
dirs[i] = -dirs[i]; // Change moving direction to move the objects back if they have moved too far
}
}
StopWatch timer;
timer.start();
octree.update();
const auto updateTime = timer.getTimeElapsed();
const auto numActiveNodes = octree.getNumActiveNodes();
const auto numAllocatedNodes = octree.getNumAllocatedNodes();
const auto maxNumPrimitivesInNodes = octree.getMaxNumPrimitivesInNodes();
std::stringstream ss;
ss << "Octree update time: " << updateTime << " ms\n"
<< "Active nodes: " << numActiveNodes
<< " (" << static_cast<double>(numActiveNodes) / static_cast<double>(numAllocatedNodes) * 100.0
<< " % usage / total allocated nodes: " << numAllocatedNodes << ")\n"
<< "Max number of primitives in tree nodes: " << maxNumPrimitivesInNodes;
statusManager->setCustomStatus(ss.str());
// Update debug rendering data
octree.updateDebugGeometry();
for (auto& delegate : renderer->getDebugRenderDelegates())
{
delegate->updateDataSource();
}
// Pause for a while
std::this_thread::sleep_for(std::chrono::milliseconds(100));
};
sdk->getSceneManager(scene)->setPostUpdateCallback(updateFunc);
// Set Camera configuration
auto cam = scene->getCamera();
cam->setPosition(Vec3d(0, 15, 50));
cam->setFocalPoint(Vec3d(0, 0, 0));
// Light
{
auto light = std::make_shared<DirectionalLight>("Light 1");
light->setFocalPoint(Vec3d(-1, -1, -1));
light->setIntensity(1);
scene->addLight(light);
}
{
auto light = std::make_shared<DirectionalLight>("Light 2");
light->setFocalPoint(Vec3d(1, -1, -1));
light->setIntensity(1);
scene->addLight(light);
}
// Run
sdk->startSimulation(SimulationStatus::PAUSED);
return 0;
}
......@@ -5,6 +5,7 @@ include(imstkAddLibrary)
imstk_add_library( Datastructures
DEPENDS
Core
Geometry
)
#-----------------------------------------------------------------------------
......
/*=========================================================================
Library: iMSTK
Copyright (c) Kitware, Inc. & Center for Modeling, Simulation,
& Imaging in Medicine, Rensselaer Polytechnic Institute.
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0.txt
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
=========================================================================*/
#include "gtest/gtest.h"
#include "gmock/gmock.h"
#include "imstkMath.h"
#include "imstkPointSet.h"
#include "imstkSurfaceMesh.h"
#include "imstkLooseOctree.h"
#include "imstkParallelUtils.h"
#include <queue>
using namespace imstk;
#define BOUND 10.0
#define SPHERE_RADIUS Real(10)
#define SPHERE_CENTER Vec3r(0, 0, 0)
#define PARTICLE_RADIUS Real(2)
#define ITERATIONS 10
///
/// \brief Generate a PointSet
///
std::shared_ptr<PointSet>
generatePointSet()
{
const Vec3r sphereCenter = SPHERE_CENTER;
const auto sphereRadiusSqr = SPHERE_RADIUS * SPHERE_RADIUS;
const auto spacing = Real(2) * PARTICLE_RADIUS;
const int N = int(2 * SPHERE_RADIUS / spacing);
StdVectorOfVec3r particles;
particles.reserve(N * N * N);
const Vec3r corner = sphereCenter - Vec3r(SPHERE_RADIUS, SPHERE_RADIUS, SPHERE_RADIUS);
for (int i = 0; i < N; ++i)
{
for (int j = 0; j < N; ++j)
{
for (int k = 0; k < N; ++k)
{
const Vec3r ppos = corner + Vec3r(spacing * Real(i), spacing * Real(j), spacing * Real(k));
const Vec3r d = ppos - sphereCenter;
if (d.squaredNorm() < sphereRadiusSqr)
{
particles.push_back(ppos);
}
}
}
}
auto pointset = std::make_shared<PointSet>();
pointset->initialize(particles);
return pointset;
}
///
/// \brief Generate a triangle soup (random discrete triangles)
///
std::shared_ptr<SurfaceMesh>
generateMesh()
{
StdVectorOfVec3r vertices;
std::vector<SurfaceMesh::TriangleArray> faces;
auto randD = [] { return (static_cast<double>(rand()) / static_cast<double>(RAND_MAX) * 2.0 - 1.0) * BOUND; };
vertices.reserve(300);
faces.reserve(100);
for (size_t i = 0; i < 100; ++i)
{
faces.push_back({ i*3, i*3 + 1, i*3 + 2 });
vertices.push_back(Vec3d(randD(), randD(), randD()));
vertices.push_back(Vec3d(randD(), randD(), randD()));
vertices.push_back(Vec3d(randD(), randD(), randD()));
}
auto mesh = std::make_shared<SurfaceMesh>();
mesh->initialize(vertices, faces);
return mesh;
}
///
/// \brief Randomize particle positions
///
void
randomizePositions(const std::shared_ptr<PointSet>& pointset)
{
for (size_t i = 0; i < pointset->getNumVertices(); ++i)
{
pointset->setVertexPosition(i, Vec3r(
(static_cast<Real>(rand()) / static_cast<Real>(RAND_MAX) * 2.0 - 1.0) * BOUND,
(static_cast<Real>(rand()) / static_cast<Real>(RAND_MAX) * 2.0 - 1.0) * BOUND,
(static_cast<Real>(rand()) / static_cast<Real>(RAND_MAX) * 2.0 - 1.0) * BOUND
));
}
}
///
/// \brief Randomize triangle positions
///
void
randomizePositions(const std::shared_ptr<SurfaceMesh>& mesh)
{
for (size_t i = 0; i < mesh->getNumTriangles(); ++i)
{
const auto translation = Vec3r(
(static_cast<Real>(rand()) / static_cast<Real>(RAND_MAX) * 2.0 - 1.0) * BOUND,
(static_cast<Real>(rand()) / static_cast<Real>(RAND_MAX) * 2.0 - 1.0) * BOUND,
(static_cast<Real>(rand()) / static_cast<Real>(RAND_MAX) * 2.0 - 1.0) * BOUND
);
const auto face = mesh->getTrianglesVertices()[i];
for (unsigned int j = 0; j < 3; ++j)
{
mesh->setVertexPosition(face[j], mesh->getVertexPosition(face[j]) + translation);
}
}
}
namespace imstk
{
class LooseOctreeTest : public ::testing::Test
{
public:
LooseOctreeTest() {}
void reset()
{
m_Octree.reset(new LooseOctree(Vec3d(0, 0, 0), 100.0, 0.1, 4));
}
void buildExample()
{
reset();
m_PointSet = generatePointSet();
m_Octree->addPointSet(m_PointSet);
EXPECT_EQ(m_Octree->m_vPrimitivePtrs[OctreePrimitiveType::Point].size(), m_PointSet->getNumVertices());
EXPECT_EQ(m_Octree->m_vPrimitivePtrs[OctreePrimitiveType::Triangle].size(), 0);
EXPECT_EQ(m_Octree->m_vPrimitivePtrs[OctreePrimitiveType::AnalyticalGeometry].size(), 0);
m_Mesh = generateMesh();
m_Octree->addTriangleMesh(m_Mesh);
EXPECT_EQ(m_Octree->m_vPrimitivePtrs[OctreePrimitiveType::Point].size(), m_PointSet->getNumVertices());
EXPECT_EQ(m_Octree->m_vPrimitivePtrs[OctreePrimitiveType::Triangle].size(), m_Mesh->getNumTriangles());
EXPECT_EQ(m_Octree->m_vPrimitivePtrs[OctreePrimitiveType::AnalyticalGeometry].size(), 0);
m_Octree->build();
}
void testOctree()
{
m_Octree->update();
// Check node pool
{
unsigned int numInactive = 0;
auto pIter = m_Octree->m_pNodeBlockPoolHead;
while (pIter) {
pIter = pIter->m_NextBlock;
numInactive += 8u;
}
EXPECT_EQ(numInactive, m_Octree->m_NumAvaiableBlocksInPool * 8u);
////////////////////////////////////////////////////////////////////////////////
unsigned int numActives = 0;
std::queue<OctreeNode*> nodeQ;
nodeQ.push(m_Octree->getRootNode());
while (!nodeQ.empty()) {
++numActives;
const auto node = nodeQ.front();
nodeQ.pop();
if (!node->isLeaf())
{
for (unsigned int childIdx = 0; childIdx < 8u; ++childIdx)
{
nodeQ.push(node->getChildNode(childIdx));
}
}
}
EXPECT_EQ(numActives, m_Octree->getNumActiveNodes());
EXPECT_EQ(numActives + numInactive, m_Octree->m_NumAllocatedNodes);
}
// Check primitives in nodes
{
unsigned int primitiveCounts[OctreePrimitiveType::NumPrimitiveTypes];
for (int i = 0; i < OctreePrimitiveType::NumPrimitiveTypes; ++i)
{
primitiveCounts[i] = 0;
}
for (const auto& it: m_Octree->m_sActiveTreeNodeBlocks)
{
for (unsigned int idx = 0; idx < 8u; ++idx)
{
for (int i = 0; i < OctreePrimitiveType::NumPrimitiveTypes; ++i)
{
primitiveCounts[i] += it->m_Nodes[idx].m_PrimitiveCounts[i];
unsigned int numPrimitives = 0;
auto pIter = it->m_Nodes[idx].m_pPrimitiveListHeads[i];
while (pIter) {
pIter = pIter->m_pNext;
++numPrimitives;
}
EXPECT_EQ(numPrimitives, it->m_Nodes[idx].m_PrimitiveCounts[i]);
}
}
}
for (int i = 0; i < OctreePrimitiveType::NumPrimitiveTypes; ++i)
{
EXPECT_EQ(primitiveCounts[i], m_Octree->m_vPrimitivePtrs[i].size());
}
}
// Check validity of primitives
{
for (int i = 0; i < OctreePrimitiveType::NumPrimitiveTypes; ++i)
{
for (const auto& pPrimitive : m_Octree->m_vPrimitivePtrs[i])
{
const auto pNode = pPrimitive->m_pNode;
if (pPrimitive->m_pGeometry->getType() == Geometry::Type::PointSet)
{
EXPECT_EQ(pNode->looselyContains(pPrimitive->m_Position), true);
}
else
{
EXPECT_EQ(pNode->looselyContains(pPrimitive->m_LowerCorner, pPrimitive->m_UpperCorner), true);
}
}
}
}
}
void testDummyPrimitives(bool bRebuild)
{
// Test points
reset();
m_Octree->setAlwaysRebuild(bRebuild);
std::vector<std::shared_ptr<PointSet>> pointsets;
for (int iter = 0; iter < ITERATIONS; ++iter)
{
const auto pointset = std::make_shared<PointSet>();
pointset->initialize({ Vec3d(rand(), rand(), rand() ) * 10.0 / static_cast<double>(RAND_MAX) });
m_Octree->addPointSet(pointset);
pointsets.push_back(std::move(pointset));
m_Octree->build();
EXPECT_EQ(m_Octree->m_MaxDepth, 10);
const auto vPrimitives = m_Octree->m_vPrimitivePtrs[OctreePrimitiveType::Point];
EXPECT_EQ(vPrimitives.size(), iter + 1);
const auto pPrimitive = vPrimitives.back();
const auto pNode = pPrimitive->m_pNode;
EXPECT_EQ(pNode->m_Depth, 10);
EXPECT_EQ(std::abs(pNode->m_HalfWidth * 2.0 - 0.1953125) < 1e-8, true);
EXPECT_EQ(pNode->looselyContains(pPrimitive->m_Position), true);
EXPECT_EQ(pNode->contains(pPrimitive->m_Position), true);
}
auto testTriangle =
[&](bool bHasPoint)
{
std::vector<std::shared_ptr<SurfaceMesh>> meshes;
for (int iter = 0; iter < ITERATIONS; ++iter)
{
StdVectorOfVec3d vertices {
Vec3d(0, 0, 0),
Vec3d(1, 0, 0),
Vec3d(1, 1, 1)
};