PMP 라이브러리 Half Edge 테스트

폴리곤 메시 처리 라이브러리중 PMP 라이브러리가 있음..

https://www.pmp-library.org/

OpenMesh 와 비교하여 비교적 심플한 사용법과..상당히 구조적으로 잘 짜여져 있음..
근데 사실.. half edge 의 특성상 half edge 한개가 2개 이상의 face를 공유할수 없는 문제로..
실제 mesh 편집 프로그램에서는 사용하기가 어렵다..

예를 들어 Solid 메시는 한개의 half edge를 두고 2개 이상을 공유해야 하는 문제가 허다하고..
shell 메세중.. 2개이상을 지원해야 하는 경우도 종종있다..

암튼.. 난 이런 사실을 모른체.. 이걸 써야겠다하고.. 맘먹고 테스트해본 코드임..
그냥 이런 라이브러리가 있구나.. 하는 정도만.. 기억할것임.
지금은 직접 만들어 쓰고 있다..

#include <SurfaceMesh.h>

int main(void)
{
	pmp::SurfaceMesh mesh;

	pmp::Vertex v0,v1,v2,v3;

    // add 4 vertices
    v0 = mesh.add_vertex(pmp::Point(0,0,0));
    v1 = mesh.add_vertex(pmp::Point(1,0,0));
    v2 = mesh.add_vertex(pmp::Point(0,1,0));
    v3 = mesh.add_vertex(pmp::Point(0,0,1));

	// temp 노드 셋팅
	pmp::VertexProperty<bool> vtmp;
    vtmp = mesh.add_vertex_property<bool>("v:temp", false);
	vtmp[v0] = true;
	vtmp[v2] = true;
	
	// solver id 셋팅
	pmp::VertexProperty<unsigned int> vsolver;
    vsolver = mesh.add_vertex_property<unsigned int>("v:solver");
	vsolver[v0] = 100;
	vsolver[v1] = 101;
	vsolver[v2] = 102;
	vsolver[v3] = 103;
    //mesh.remove_vertex_property(vsolver);

    auto points = mesh.get_vertex_property<bool>("v:temp");

    for (auto v : mesh.vertices())
    {
        if ( points[v] ) {
			printf("%d : temp node, solver id: %d\n", v.idx(), vsolver[v]);
		} else {
			printf("%d : node, solver id: %d\n", v.idx(), vsolver[v]);
		}
    }

    // add 4 triangular faces
	pmp::Face f0 = mesh.add_triangle(v0,v1,v3);
    mesh.add_triangle(v1,v2,v3);
    mesh.add_triangle(v2,v0,v3);
    mesh.add_triangle(v0,v2,v1);
	
	// face에 대한 solver id 셋팅
	pmp::FaceProperty<unsigned int> fsolver;
    fsolver = mesh.add_face_property<unsigned int>("f:solver");
	fsolver[f0] = 10000;
	
	// face에 대한 part id 셋팅
	pmp::FaceProperty<unsigned int> fpart;
    fpart = mesh.add_face_property<unsigned int>("f:part");
	fpart[f0] = 1001;

	// part 1001 번에 대한 opengl 출력용 tri indices
	std::map<unsigned int, std::vector<unsigned int>> indices;
	std::map<pmp::Face, unsigned int> offset;
	offset[f0].push_back(indices[1001].size());
	indices[1001].push_back(v0.idx());
	indices[1001].push_back(v1.idx());
	indices[1001].push_back(v3.idx());
	offset[f1].push_back(indices[1001].size());
	indices[1001].push_back(v1.idx());
	indices[1001].push_back(v2.idx());
	indices[1001].push_back(v3.idx());

    std::cout << "vertices: " << mesh.n_vertices() << std::endl;
    std::cout << "edges: "    << mesh.n_edges()    << std::endl;
    std::cout << "faces: "    << mesh.n_faces()    << std::endl;
}

robin hood hashing

https://github.com/Tessil/robin-map

std::map과 std::unordered_map 을 속도면에서 발라버리는.. 헤쉬맵..
사용법은 std와 동일하다..

단순히 std::map과 std::unordered_map 을 tsl::robin_map 으로 바꾸어서 컴파일하니..
노드 1488581 개를 가진 메시를 기준으로..
메시 구조를 빌딩하는 타임이 26초 걸리던것이.. 18초가 되었다.. 지쟈스..

clear도 빨라졌네.. ㅎㅎ

여러 해쉬 라이브러리와 비교한 밴치마킹 사이트..
https://martin.ankerl.com/2019/04/01/hashmap-benchmarks-03-01-result-InsertHugeInt/

Line to Quad shader

OpenGL에서 Line을 Quad로 변경하는 OpenGL Geometry 쉐이더..
screen_size 는 viewport width, height 의 값.

아래는 메시의 edge (wireframe) 를 그린것인데..

line_width를 1.0으로 하면.. Line 처럼 보이지만..

10.0으로 변경하면..

요렇게 Quad Primitive 로 라인을 그린걸 볼수 있음.

이렇게 Line을 Quad로 그리는 이유는
1. 속도 때문. 태초에 OpenGL 은 Triangle (Quad는 내부적으로 Triangle 두개로 표현됨) 을 그리는데 빠른 성능을 내기때문임.
2. Line은 offset (glPolygonOffset)을 할수 없어서 Face 위에 정확하게 라인을 그리기가 어려움. 이렇게 Quad로 변경해서 그리면 Offset 이 가능하므로 선을 surface 위에 이쁘게 그릴수가 있음.

#version 400 compatibility
layout(lines) in;
layout(triangle_strip, max_vertices = 4) out;

uniform float line_width;
uniform vec2 screen_size;

void main()
{
  // start/end points of line
  vec4 p0 = gl_in[0].gl_Position;
  vec4 p1 = gl_in[1].gl_Position;
  
  // convert to screen space
  p0.xy = p0.xy / p0.w * screen_size;
  p1.xy = p1.xy / p1.w * screen_size;

  // compute dir and normal
  vec2 lineDir = p1.xy - p0.xy;
  vec2 lineNormal = normalize(vec2(-lineDir.y, lineDir.x));
  
  
  // create screen-aligned quad
  vec2 offset = lineNormal * line_width;


  gl_Position = vec4( (p0.xy + offset) * p0.w / screen_size, p0.z, p0.w);
  EmitVertex();
  
  gl_Position = vec4( (p0.xy - offset) * p0.w / screen_size, p0.z, p0.w);
  EmitVertex();
  
  gl_Position = vec4( (p1.xy + offset) * p1.w / screen_size, p1.z, p1.w);
  EmitVertex();
  
  gl_Position = vec4( (p1.xy - offset) * p1.w / screen_size, p1.z, p1.w);
  EmitVertex();
  
  EndPrimitive();
}

KD-Tree 예제

공간 검색에 탁월하다는 KD-Tree 알고리즘..
http://blog.daum.net/pg365/140 이 싸이트에 간단하게 구현해놓은 코드가 있어 테스트해봄.

kd_tree.h 다운로드

이 코드의 문제점은..
1. insert 속도가 너무 느리다
2. 밸런싱 기능이 없다
3. 삭제 기능이 없다

#include <stdio.h>
#include <assert.h>
#include <time.h>
#include <windows.h>
#include "kd_tree.h"

unsigned int get_msec(void)
{
	return GetTickCount ();
}

class KDTree : public kd_tree<double, 3>
{
public:
	void insert_xyz (double x, double y, double z, void *data)
	{
		double buf[3] = { x, y, z };
		
		insert (buf, data);
	}

	kd_node *nearest_neighbour_search (double x, double y, double z)
	{
		double pos[3] = {x, y, z};
		
		return nn_search (pos);
	}

	list<kd_node *> *bounding_box_search (double lx, double ly, double lz, double hx, double hy, double hz)
	{
		double low[3]  = { lx, ly, lz };
		double high[3] = { hx, hy, hz };

		return range_search (low, high);
	}
};


int main(int argc, char **argv)
{
	KDTree kd;

	unsigned int start = get_msec();

	int point_id = 0;

	printf("-- insert points ... \n"); fflush(stdout);

	for(int i=0; i<1000; i++) {

		for (int x=0; x < 10; x++) {
			for (int y=0; y < 10; y++) {
				for (int z=0; z < 10; z++) {
					kd.insert_xyz (x, y, z, (void *)point_id);
					point_id++;
				}
			}
		}

		for (int x=20; x < 30; x++) {
			for (int y=20; y < 30; y++) {
				for (int z=0; z < 10; z++) {
					kd.insert_xyz (x, y, z, (void *)point_id);
					point_id++;
				}
			}
		}

		for (int x=0; x < 10; x++) {
			for (int y=20; y < 30; y++) {
				for (int z=0; z < 10; z++) {
					kd.insert_xyz (x, y, z, (void *)point_id);
					point_id++;
				}
			}
		}

		for (int x=20; x < 30; x++) {
			for (int y=0; y < 10; y++) {
				for (int z=0; z < 10; z++) {
					kd.insert_xyz (x, y, z, (void *)point_id);
					point_id++;
				}
			}
		}
	}

	unsigned int msec = get_msec() - start;
	printf("inserting %d points... %.3f sec\n", point_id, msec/1000.0);

	start = get_msec();

	list<kd_tree<double,3>::kd_node *> *ret = kd.bounding_box_search(0, 0, 0, 10, 10, 10);
#if 1
	for (list<kd_tree<double,3>::kd_node *>::iterator it = ret->begin (); it != ret->end(); it++) {
		kd_tree<double,3>::kd_node *n = *it;
		printf ("(%g, %g, %g) \n", n->_pos[0], n->_pos[1], n->_pos[2]);
	}
#endif

	msec = get_msec() - start;
	printf("range query returned %d items in %.5f sec\n", ret->size(), msec/1000.0);
	return 0;
}

OpenGL 로 만드는 해석 모델 편집기..

요즘 짜투리 시간을 이용해서..
OpenGL로 해석 모델을 편집하는 프로그램을 만들고 있다..
이래저래 문제도 많았는데..
하나둘씩 해결이 되고 있다..

투명문제도.. 텍스트 출력 속도도.. 해결이 되었음..

파트가 너무 많으면.. 인덱스 트랜스퍼가 너무 많다보니.. 느림..
이건 꼭 해결하고 넘어가야할 산…

아래는 인터넷상에 굴러다니는 모델을 불러온 화면..

가장 근접 거리의 포인트 검색

가장 근접거리의 포인트를 단시간내에 찾는 알고리즘중에..
octree, kd-tree 가 제일 유명한듯하다..
은근히 이 알고리즘이 구현된 쓸만한 라이브러리가 없다..
찾아보다 결국 PCL (Point cloud library) 이라는 라이브러리가 제일 유명하더라는..
모든 포인트를 클라우드라는 영역에 집어넣고 특정 포인트에서 쿼리를 때리면..
가장 근접한 몇개의 포인트를 찾아준다..

보통 octree나 kd-tree는 방향성을 주어 검색을 할 수 없다..
PCL은 다행히 octree 에 한정하여 원점과 방향을 지정하여 근접 포인트 쿼리가 가능하다.

아래는 테스트한 코드..

pcl::PointCloud<pcl::PointXYZL>::Ptr m_point_cloud (new pcl::PointCloud<pcl::PointXYZL> ());
pcl::octree::OctreePointCloudSearch<pcl::PointXYZL> m_octree_search (0.2f);

// 포인트와 해당 ID를 집어넣는다.
for(int i=0; i<10000; i++) {
	...
	m_point_cloud->push_back (pcl::PointXYZL ((float)px, (float)py, (float)pz, id));
}

// 원점 지정
Eigen::Vector3f _origin (origin[0],origin[1],origin[2]);
// 방향 지정
Eigen::Vector3f _ray (ray[0],ray[1],ray[2]);

// resolution 지정 (반지름 30을 갖는 sphre 정도라고 생각하면 될려나?)
float resolution = 30.f;

printf("octree ray intersection with resolution %f\n", resolution); 
fflush(stdout);

// 쿼리
std::vector<int> indicesInRay;
pcl::octree::OctreePointCloudSearch<pcl::PointXYZL> octree(resolution);
octree.setInputCloud (m_point_cloud);
octree.addPointsFromInputCloud ();
octree.getIntersectedVoxelIndices (_origin, _ray, indicesInRay);

// 탐색된 포인트 확인
for(int i=0; i<indicesInRay.size(); i++) {
	int idx = indicesInRay[i];
	pcl::PointXYZL pt = m_point_cloud->points[idx];
	//printf("[%d]: %d: %f, %f, %f\n", idx, pt.label, pt.x, pt.y, pt.z);
}

octree.deleteTree();

Outline 그리기

opengl 에서 outline을 그리는 몇가지 방법중의 하나..
바로 stencil buffer 를 이용함.

glPushAttrib(GL_ALL_ATTRIB_BITS);

glClearStencil(0);
glClear(GL_STENCIL_BUFFER_BIT);
glEnable(GL_STENCIL_TEST);
// Set the stencil buffer to write a 1 in every time
// a pixel is written to the screen
glStencilFunc(GL_ALWAYS, 1, 0xFFFF);
glStencilOp(GL_KEEP, GL_KEEP, GL_REPLACE);
glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
glColor3f(1.0f, 0.0f, 0.0f);

...draw primitive...

glStencilFunc(GL_NOTEQUAL, 1, 0xFFFF);
glStencilOp(GL_KEEP, GL_KEEP, GL_REPLACE);
// Draw the object with thick lines
glLineWidth(4.0f);
glDisable (GL_LINE_SMOOTH);
glDisable(GL_LIGHTING);

glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
glColor3f(1.0f, 1.0f, 1.0f);

...draw primitive...

glPopAttrib();

스트링 offset 파싱 속도 비교

다음과 같은.. 스트링이 있을때..

char *s = "10000001   1201712009430120094781202507512025075";

8,8,8,8,8,8,8 자리로 끊어서 읽고 싶은데..
sscanf는 공백을 무시고 읽어.. 원하는데로 동작되지 않음..
하여 찾다보니.. boost에 offset_separator 라는것이 있음..
테스트 하는 김에 속도 비교를 해봄..

sscanf의 경우.. (대략 88.637946초)

std::chrono::system_clock::time_point start = std::chrono::system_clock::now();
for(int i=0; i<100000000; i++) {
	char *s = "10000001   1201712009430120094781202507512025075";
	unsigned int eid, pid, n1, n2, n3, n4;
	int cnt = sscanf(s, "%8u%8u%8u%8u%8u%8u", &eid, &pid, &n1, &n2, &n3, &n4);
}
std::chrono::system_clock::time_point end = std::chrono::system_clock::now();
std::chrono::duration<double> sec = end - start;
printf("%f(sec)}", sec.count());

boost의 offset_separator (대략 72.378287초)

std::chrono::system_clock::time_point start = std::chrono::system_clock::now();
for(int i=0; i<100000000; i++) {
	std::string s = "10000001   1201712009430120094781202507512025075";

   int offsets[] = {8,8,8,8,8,8};
   boost::offset_separator f(offsets, offsets+6);
   boost::tokenizer<boost::offset_separator> tok(s,f);
   //for(boost::tokenizer<boost::offset_separator>::iterator beg=tok.begin(); beg!=tok.end();++beg){
	 //  std::cout << *beg << "\n";
  // }
   const size_t count = std::distance(tok.begin(), tok.end());
   //std::cout << count << "\n";
   boost::tokenizer<boost::offset_separator>::iterator it1 = tok.begin();
   //std::cout << *it1 << "\n";
   ++it1;
   //std::cout << *it1 << "\n";
   ++it1;
   //std::cout << *it1 << "\n";
   ++it1;
   //std::cout << *it1 << "\n";
   ++it1;
   //std::cout << *it1 << "\n";
   ++it1;
   //std::cout << *it1 << "\n";
}
std::chrono::system_clock::time_point end = std::chrono::system_clock::now();
std::chrono::duration<double> sec = end - start;
printf("%f(sec)}", sec.count());

역시.. 구글링을 해보면 sscanf의 속도가 느리다고 하는데.. 맞는가봄..
일단 boost의 경우 자리수를 지정해서 그 만큼 짤라내는게 아주 잘 되고..
sscanf보다 빠르니.. 요걸 쓰기로 결정..

gmsh를 이용한 Meshing..

자.. 우선 gmsh 의 입력 geometry를 작성..

이건 cube 형태의 geometry..

//test.geo
Point(1) = {1, 0, 0, 1.0};
Point(2) = {1, 1, 0, 1.0};
Point(3) = {0, 1, 0, 1.0};
Point(4) = {0, 0, 1, 1.0};
Point(5) = {1, 0, 1, 1.0};
Point(6) = {1, 1, 1, 1.0};
Point(7) = {0, 1, 1, 1.0};
Point(8) = {0, 0, 0, 1.0};
Line(1) = {7, 6};
Line(2) = {6, 5};
Line(3) = {5, 1};
Line(4) = {1, 8};
Line(5) = {8, 3};
Line(6) = {3, 7};
Line(7) = {7, 4};
Line(8) = {4, 8};
Line(9) = {4, 5};
Line(10) = {2, 1};
Line(11) = {2, 6};
Line(12) = {2, 3};
Line Loop(1) = {6, 1, -11, 12};
Plane Surface(1) = {1};
Line Loop(2) = {11, 2, 3, -10};
Plane Surface(2) = {2};
Line Loop(3) = {2, -9, -7, 1};
Plane Surface(3) = {-3};
Line Loop(4) = {6, 7, 8, 5};
Plane Surface(4) = {-4};
Line Loop(5) = {8, -4, -3, -9};
Plane Surface(5) = {5};
Line Loop(6) = {10, 4, 5, -12};
Plane Surface(6) = {6};
Physical Surface(1) = {4, 3, 2, 6};
Physical Surface(2) = {1};
Physical Surface(3) = {5};
Surface Loop(1) = {6, 2, 1, 4, 3, 5};
Volume(1) = {1};

이제 gmsh 옵션으로.. 2d (-2), 알고리즘은 meshadapt를.. 아웃풋 포맷은 LsDyna의 key 포맷으로..
요소 최소/최대 사이즈를 지정하면..

gmsh.exe test.geo -2 -algo meshadapt -o test.key -format key -clmin 0.01 -clmax 0.05
...
Info    : Meshing 2D...
Info    : Meshing surface 1 (Plane, Frontal Quad)
Info    : Meshing surface 2 (Plane, Frontal Quad)
Info    : Meshing surface 3 (Plane, Frontal Quad)
Info    : Meshing surface 4 (Plane, Frontal Quad)
Info    : Meshing surface 5 (Plane, Frontal Quad)
Info    : Meshing surface 6 (Plane, Frontal Quad)
Info    : Done meshing 2D (0.109375 s)
Info    : 2402 nodes 5048 elements
Info    : Writing 'test.key'...
Info    : Done writing 'test.key'

짜잔~

요렇게 생성이 됨.

자.. 이제 quad 메싱은 맨 위의 입력 geom 파일 끝 줄에 아래를 추가..

Recombine Surface "*";

다음.. 아래와 같이..하고 algorithm은 delquad를 사용하여 매싱을 하면..

gmsh.exe test.geo -2 -algo delquad -o test.key -format key -clmin 0.01 -clmax 0.05
...
Info    : Meshing surface 6 (Plane, Frontal Quad)
Info    : Blossom: 1160 internal 76 closed
Info    : Blossom recombination completed (0 s): 400 quads, 0 triangles, 0 invalid quads, 0 quads with Q < 0.1, avg Q = 1, min Q = 1
Info    : Done meshing 2D (0.09375 s)
Info    : 2402 nodes 2648 elements
Info    : Writing 'test.key'...
Info    : Done writing 'test.key'

오… 꽤 쓸만함..

근데.. 문제는 gmsh의 경우 tria + quad 의 mixed 는 안된다고함..
단 한 surface는 tria를.. 다른 surface는 quad 는 가능함..

무슨 말이냐하면.. 아래와 같이 특정 1번 surface만 recombine을 적용하면..

//Recombine Surface "*";
Recombine Surface{1};

한면만 quad로 meshing이된다..

다음 tetra 는 직접 step 캐드 파일로부터 해봄..
다음과 같은 test.step 파일이 있을때..

아래와 같이 geom 파일을 작성하고..

// stl2msh.geo
Merge "test.step";
Surface Loop(1) = 1;
Volume(1) = 1;
Physical Volume("obj") = {1};

//Mesh.Algorithm3D = 4; //(1=tetgen, 4=netgen, 7=MMG3D, 9=R-tree)
Mesh.Recombine3DAll=1;
Mesh.Smoothing=0;
Mesh.Optimize=1;
Mesh.OptimizeNetgen=1;

gmsh를 아래와 같이 돌리면..

gmsh.exe stl2msh.geo -3 -o test.key -format key -clmax 1.0 -clmin 0.1

짜잔.. 잘 된것 같은데.. 문제는 형상이 몇개 빠진듯.. -_-;;
암튼.. 메시 사이즈를 0.4 정도로 주고 다시 해보면..

잘되는듯하다..

gmsh의 경우 hexa가 안되는게 제일 문제..

종합해보면..

1. quad, tria 지원 (mixed 미지원)
2. tetra 지원
3. hexa 미지원

아래는 각 메시 라이브러리별 지원 항목..

ReactPhysics3D 의 광선 교차 테스트

ReactPhysics3D 의 광선 교차 테스트..

#include "reactphysics3d.h"

#include <iostream>
#include <stdexcept>
#include <vector>
#include <map>
#include <sstream>
#include <iomanip>
#include <memory>
#include <algorithm>
#include <set>
#include <string>
#include <fstream>

#define GLM_ENABLE_EXPERIMENTAL
#include <glm/gtx/string_cast.hpp>
#include <glm/glm.hpp>
#include <glm/gtc/matrix_transform.hpp>
#include <glm/gtc/type_ptr.hpp>
#include <glm/gtx/transform.hpp> 
#include <glm/gtx/quaternion.hpp>
#include <glm/gtx/normal.hpp>


// https://github.com/DanielChappuis/reactphysics3d/blob/master/test/tests/collision/TestAABB.h
using namespace reactphysics3d;


class DynamicTreeRaycastCallback : public DynamicAABBTreeRaycastCallback {

    public:

        std::vector<int> mHitNodes;

        // Called when the AABB of a leaf node is hit by a ray
        virtual decimal raycastBroadPhaseShape(int32 nodeId, const Ray& ray) override {
			printf("--hit-- %d\n", nodeId);
            mHitNodes.push_back(nodeId);
            return 1.0;
        }

        void reset() {
            mHitNodes.clear();
        }

        bool isHit(int nodeId) const {
            return std::find(mHitNodes.begin(), mHitNodes.end(), nodeId) != mHitNodes.end();
        }
};


int main(void)
{
	DynamicTreeRaycastCallback mRaycastCallback;

	DynamicAABBTree tree(MemoryManager::getBaseAllocator());

	AABB aabb;

	std::vector<glm::vec3> vertices;

	int solverId = 1;
	vertices.push_back(glm::vec3(0, -2.5, 0));
	vertices.push_back(glm::vec3(-2.5, 2.5, 0));
	vertices.push_back(glm::vec3(2.5, 2.5, 0));
	aabb = AABB::createAABBForTriangle((Vector3*)&vertices[0]);
	int objectId = tree.addObject(aabb, &solverId);
	printf("%d\n", objectId);

	solverId = 2;
	glm::vec3 points2[] = {
		glm::vec3(0, -2.5, 2), glm::vec3(-2.5, 2.5, 2), glm::vec3(2.5, 2.5, 2)
	};
	aabb = AABB::createAABBForTriangle((Vector3*)points2);
	objectId = tree.addObject(aabb, &solverId);
	printf("%d\n", objectId);
	
	solverId = 3;
	vertices.clear();
	vertices.push_back(glm::vec3(0, -2.5, 4));
	vertices.push_back(glm::vec3(-2.5, 2.5, 4));
	vertices.push_back(glm::vec3(2.5, 2.5, 4));
	aabb = AABB::createAABBForTriangle((Vector3*)&vertices[0]);
	objectId = tree.addObject(aabb, &solverId);
	printf("%d\n", objectId);


	printf("Press any key to ray hit test..\n"); fflush(stdout);
	getchar();

	Ray ray(Vector3(0, 0, 10), Vector3(0, 0, -10));

	// hit test 
	printf("Hit test..\n");
	mRaycastCallback.reset();
	tree.raycast(ray, mRaycastCallback);

	// hit test after remove
	printf("Hit test after remove..\n");
	tree.removeObject(1);
	mRaycastCallback.reset();
	tree.raycast(ray, mRaycastCallback);

#if 0
	printf("%f\n", tree.getMin().x);
	printf("%f\n", tree.getMin().y);
	printf("%f\n", tree.getMin().z);
	printf("%f\n", tree.getMax().x);
	printf("%f\n", tree.getMax().y);
	printf("%f\n", tree.getMax().z);
#endif

#if 0
	printf("%d\n", aabb.testRayIntersect(ray));
#endif

	return 0;
}