RayTracing BVH problém – C / C++ – Fórum – Programujte.com

 Dělám projekt RayTracing a zasekl jsem se na nepříjemném problému ohledně BVH stromu. Jedná se konkrétně o část s procházením, vytvoření se mně zdá být v pohodě. Jedná se mně o to, že pokud použije model zajíce s 947 trojúhelníky (standfort bunny), tak pokud nechám zajíce vykreslit pomocí vytvořeného BVH stromu a jeho procházení, tak se vykreslí shluk nesmyslných obrazců. Pokud ovšem využiji pouze test průniků paprsků s trojúhelníkem (947 volání na každý paprsek), tak se vytvoří obrys zajíce korektně. V této fázi by mně stačil zatím obrys zajíce.

Zdrojový kód vytváření BVH stromu: začíná se na ose Z axis = 2, items_ obsahuje všechny trojúhelníky načtené z modelu zajíce.

Node * BVH::BuildTree(int from, int to, const int depth, char axis)
{
	const int n = to - from + 1;
	Node * node = new Node(from, to);

	for (int i = from; i <= to; ++i)
	{
		node->bounding.Merge(items_[i].Bounds());
	}

	if (n <= max_leaf_items_)
	{
		leafs_++;
		max_depth_ = MAX(depth, max_depth_);
		return node;
	}

	// pokud to neni leaf -> node
	nodes_ += 2;

	const int pivot = n / 2 + from;

	Sort(&items_[from], n, axis);

	// rozdeleni na dva podintervaly
	axis = (axis + 1) % 3;

	// leafy
	node->children[0] = BuildTree(from, pivot - 1, depth + 1, axis);
	node->children[1] = BuildTree(pivot, to, depth + 1, axis);

	return node; // vraci koren
}

Procházení stromu:

void BVH::Traverse(Ray & ray)
{
	Traverse(ray, root_);
}

void BVH::Traverse(Ray & ray, Node * node)
{
	if (RayBoxIntersection(ray, node->bounding))
	{
		if (node->IsLeaf())
		{
			// projdeme trojuhelniky
			for (int i = node->span[0]; i <= node->span[1]; i++) RayTriangleIntersection97(&items_[i], ray);
		}
		else
		{
			Traverse(ray, node->children[0]);
			Traverse(ray, node->children[1]);
		}
	}
}

Procházení má fungovat tak, že se zjistí průnik paprsku s celým AABB boxem kořene BVH stromu, a pokud se nalezne tento průnik, tak se zjišťuje až do jaké pod části stromu BVH průnik spadá a v ní se pak zavolá již průnik jednotlivých trojúhelníků s paprskem. Věděl by někdo co s tím ? Nebo mám přihodit i fotky problému a zbylý zdrojový kód ?

Bez zbytku by ti pomohl asi jen někdo, kdo už to přímo neřešil, takže to sem hoď a uvidíme :)

bool BVH::RayBoxIntersection(Ray & ray, AABB & box)
{
	Vector3 r_dir = Vector3(1 / ray.direction.x, 1 / ray.direction.y, 1 / ray.direction.z);

	float
	l1	= (box.bounds[0].x - ray.origin.x) * r_dir.x,
	l2	= (box.bounds[1].x - ray.origin.x) * r_dir.x,
	lmin	= MIN(l1,l2),
	lmax	= MAX(l1,l2);

	l1	= (box.bounds[0].y - ray.origin.y) * r_dir.y;
	l2	= (box.bounds[1].y - ray.origin.y) * r_dir.y;
	lmin	= MAX(MIN(l1,l2), lmin);
	lmax	= MIN(MAX(l1,l2), lmax);
		
	l1	= (box.bounds[0].z - ray.origin.z) * r_dir.z;
	l2	= (box.bounds[1].z - ray.origin.z) * r_dir.z;
	lmin	= MAX(MIN(l1,l2), lmin);
	lmax	= MIN(MAX(l1,l2), lmax);
		
	//return ((lmax > 0.f) & (lmax >= lmin));
	//return ((lmax > 0.f) & (lmax > lmin));
	return ((lmax >= 0.f) & (lmax >= lmin));
}

void BVH::RayTriangleIntersection97(Triangle * item, Ray & ray)
{
	#define SMALL_NUM  0.00000001 // anything that avoids division overflow

	Vector3 u, v, n; // triangle vectors
	Vector3 dir, w0, w; // ray vectors
	REAL r, a, b; // params to calc ray-plane intersect

	// get triangle edge vectors and plane normal
	u = item->vertices[1] - item->vertices[0];
	v = item->vertices[2] - item->vertices[0];
	n = u.CrossProduct(v); // cross product
	if (n.IsEmpty()) return; // triangle is degenerate, do not deal with this case

	// ray direction vector
	dir = ray.direction - ray.origin;
	w0 = ray.origin - item->vertices[0];
	a = (-1) * n.DotProduct(w0); // -dot(n, w0);
	b = n.DotProduct(dir);

	// ray is parallel to triangle plane
	if (fabs(b) < SMALL_NUM) return;

	// get intersect point of ray with triangle plane
	r = a / b;
	// ray goes away from triangle
	if (r < 0.0) return; // no intersect
	// for a segment, also test if (r > 1.0) => no intersect
	Vector3 output = Vector3(0, 0, 0); // vysledny vektor - I
	output = ray.origin + r * dir;           // intersect point of ray and plane
		
	// is I inside T? je bod output uvnitr trojuhelniku item
	REAL uu, uv, vv, wu, wv, D;
	uu = u.DotProduct(u);
	uv = u.DotProduct(v);
	vv = v.DotProduct(v);
	w = output - item->vertices[0];
	wu = w.DotProduct(u);
	wv = w.DotProduct(v);
	D = uv * uv - uu * vv;

	// get and test parametric coords
	REAL s, t;
	s = (uv * wv - vv * wu) / D;
	// I is outside T
	if (s < 0.0 || s > 1.0) return;
	t = (uv * wu - uu * wv) / D;
	// I is outside T
	if (t < 0.0 || (s + t) > 1.0) return;

	// nastavime nove t
	ray.SetT(r, item);
}

void BVH::Sort(Triangle * items, int n, char axis)
{
	#define MAX_LEVELS  300 

	int beg[MAX_LEVELS];
	int end[MAX_LEVELS];

	int i=0, L, R, swap ;
	Triangle piv;

	beg[0]=0; end[0]=n;
			
	while (i>=0) {
	L=beg[i]; R=end[i]-1;
	if (L<R) {
		piv=items[L];
		while (L<R) {
			while (items[R].Bounds().Centroid().data[axis] >=piv.Bounds().Centroid().data[axis] && L<R) R--; if (L<R) items[L++]=items[R];
			while (items[L].Bounds().Centroid().data[axis] <= piv.Bounds().Centroid().data[axis] && L < R) L++; if (L < R) items[R--] = items[L];
		}
		items[L]=piv; beg[i+1]=L+1; end[i+1]=end[i]; end[i++]=L;
		if (end[i]-beg[i]>end[i-1]-beg[i-1]) {
			swap=beg[i]; beg[i]=beg[i-1]; beg[i-1]=swap;
			swap=end[i]; end[i]=end[i-1]; end[i-1]=swap; }}
	else {
	i--; }}
}

BVH::BVH(Geometry * geometry, int max_leaf_items_p_)
{
	max_leaf_items_ = max_leaf_items_p_;
	leafs_ = nodes_ = max_depth_ = 0;
	int items_count_ = geometry->number_of_faces();
	items_ = new Triangle[items_count_];
	for (int i = 0; i < items_count_; i++) items_[i] = geometry->GetTriangle(i);
	Triangle * tri = &items_[0];
	printf("Build tree\n");
	printf("Start...\n");
	root_ = BuildTree(0, items_count_ - 1, 0, 2);
	nodes_++; // jeden si pripocist za root samotny
	// vypis
	printf("Leafs : %d\n", leafs_);
	printf("Nodes : %d\n", nodes_);
	printf("Items : %d\n", items_count_);
	printf("Max d.: %d\n", max_depth_);
	printf("Leaf items: %d\n", max_leaf_items_);
	// bounds=(-0.095, -0.062, 0.033) x (0.061, 0.059, 0.187)
	printf("Bounds=(%.3f, %.3f, %.3f) x (%.3f, %.3f, %.3f)\n",
		 root_->bounding.bounds[0].x, root_->bounding.bounds[0].y, root_->bounding.bounds[0].z,
		 root_->bounding.bounds[1].x, root_->bounding.bounds[1].y, root_->bounding.bounds[1].z);
}

struct Node
{
	AABB bounding; // omezujici kvadr vsech itemu v uzlu
	int span[2]; // uzavreny interval vsech indexu
	Node * children[2]; // dva potomci max. leva a prava cast

	bool IsLeaf()
	{
		return ((children[0] == NULL) && (children[1] == NULL));
	}

	Node(const int from, const int to)
	{
		span[0] = from; 
		span[1] = to;
		children[0] = children[1] = NULL;
	}

	// destruktor
	~Node()
	{
		// smazani obou potomku, pokud je to mozne
		if (children[0] != NULL)
		{
			delete children[0];
			children[0] = NULL;
		}
		if (children[1] != NULL)
		{
			delete children[1];
			children[1] = NULL;
		}
	}
};

struct AABB
{
	Vector3 bounds[2];

	AABB()
	{
		bounds[0] = Vector3(MAX_REAL, MAX_REAL, MAX_REAL); // minimum nastaveno na maximum
		bounds[1] = Vector3(-MAX_REAL, -MAX_REAL, -MAX_REAL); // maximum nastaveno na minimum
	}

	AABB(const Vector3 & p0, const Vector3 & p1)
	{
		bounds[0] = p0; // min
		bounds[1] = p1; // max
	}

	void Merge(AABB & aaBB)
	{
		// minima z nimin
		bounds[0].x = MIN(bounds[0].x, aaBB.bounds[0].x);
		bounds[0].y = MIN(bounds[0].y, aaBB.bounds[0].y);
		bounds[0].z = MIN(bounds[0].z, aaBB.bounds[0].z);

		// maxima z maxim
		bounds[1].x = MAX(bounds[1].x, aaBB.bounds[1].x);
		bounds[1].y = MAX(bounds[1].y, aaBB.bounds[1].y);
		bounds[1].z = MAX(bounds[1].z, aaBB.bounds[1].z);
	}

	Vector3 Centroid()
	{
		// secteme a podelime 2, vratime
		return Vector3((bounds[0].x + bounds[1].x) / 2, (bounds[0].y + bounds[1].y) / 2, (bounds[0].z + bounds[1].z) / 2);
	}
};

struct Triangle
{
	Vector3 vertices[3];

	Triangle()
	{
		vertices[0] = Vector3(0, 0, 0);
		vertices[1] = Vector3(0, 0, 0);
		vertices[2] = Vector3(0, 0, 0);
	}
	// x, y, z -> u modelu prohodit
	Triangle( Vector3 & p0, Vector3 & p1, Vector3 & p2 )
	{
		vertices[0] = p0;
		vertices[1] = p2;
		vertices[2] = p1;
	}

	AABB Bounds()
	{
		AABB output;

		// min vektor
		output.bounds[0].x = MIN(MIN(vertices[0].x, vertices[1].x), vertices[2].x);
		output.bounds[0].y = MIN(MIN(vertices[0].y, vertices[1].y), vertices[2].y);
		output.bounds[0].z = MIN(MIN(vertices[0].z, vertices[1].z), vertices[2].z);

		// max vektor
		output.bounds[1].x = MAX(MAX(vertices[0].x, vertices[1].x), vertices[2].x);
		output.bounds[1].y = MAX(MAX(vertices[0].y, vertices[1].y), vertices[2].y);
		output.bounds[1].z = MAX(MAX(vertices[0].z, vertices[1].z), vertices[2].z);

		return output;
	}

	// nebo 3 normaly v kazdem rohu a aritmeticky prumer deleny 3
	Vector3 Normal(Vector3 & p)
	{
		Vector3 v0  = vertices[2] - vertices[0];
		Vector3 v1 = vertices[1] - vertices[0];
		Vector3 v2 = p - vertices[0];

		REAL dot00 = v0.DotProduct(v0);
		REAL dot01 = v0.DotProduct(v1);
		REAL dot02 = v0.DotProduct(v2);
		REAL dot11 = v1.DotProduct(v1);
		REAL dot12 = v1.DotProduct(v2);

		REAL inv_denom = 1 / (dot00 * dot11 - dot01 * dot01);
		REAL u = (dot11 * dot02 - dot01 * dot12) * inv_denom;
		REAL v = (dot00 * dot12 - dot01 * dot02) * inv_denom;

		v1 *= u;
		v2 *= v;

		// vracime
		return v1.Normal(v2);
	}
};

struct Ray
{
	Vector3 origin;
	Vector3 direction;
	Vector3 backDirection; // vektor pozadi
	REAL t;
	bool changed;
	Triangle * triangle; // zasazeny trojuhelnik

	Ray( Vector3 & origin, Vector3 & direction )
	{
		this->origin = origin;
		this->direction = direction;
		t = MAX_REAL;
		changed = false;
		triangle = NULL;
	}

	// nastavi t, ale nove t musi byt mensi jak posledni jinak se nic nestane
	void SetT(float t)
	{
		if (t < this->t)
		{
			changed = true; // jednou se jiz zmenilo, zasah
			this->t = t;
		}
	}

	void SetT(float t, Triangle * tri)
	{
		if (t < this->t)
		{
			changed = true; // jednou se jiz zmenilo, zasah
			this->t = t;
			triangle = tri;
		}
	}

	Vector3 Target( const REAL t )
	{
		return origin + t * direction;
	}	
};

Mně to funguje špatně i když testuju všechny trojúhelníky a fakt se mi nechce hledat, jestli tam máš chybu ty nebo já :)

#4 zlz
průsečík trojúhelníků vím že funguje správně, to mám odzkoušené. V kombinaci s Traverse a AABB už se ale vykresluje špatně.

Rozhodl jsem se zaplatit převodem na účet 500 Kč, tomu kdo mně tam najde tu chybu. Poslal bych Vám celé solution ve Visual Studiu 2010, jazyk C++ s knihovnou OpenCl. Vážní zájemci piště na email Roman.Makudera@seznam.cz

Jde o to, že to je spousta kódu a další hromadu jsem musel dopsat, abych to zprovoznil. Hoď to sem komplet (jako odkaz), platit nic nemusíš.

Tady je celé solution : http://frasier.ic.cz/…GI_Final.zip

Sakra, taková kravina, ale taky jsem to nemohl najít.

bool BVH::RayBoxIntersection(Ray & _ray, AABB & box)
{
    Ray ray = _ray;
    ray.direction -= ray.origin;

    Vector3 r_dir = Vector3(1 / ray.direction.x, 1 / ray.direction.y, 1 / ray.direction.z);
    ...
}

Díky moc zlz, to byla ta chyba. Nejhorší je, že já k tomu nemám žádné materiály a od těch verzí tohoto algoritmu co jsem si našel na internetu, nikdo neříkal že se tyto dvě části musí odečíst. Díky moc, už to funguje.

Nemusí se odečíst, pokud 'direction' je směrový vektor. Ale když máš ray definovaný jinak, tak to holt musíš přizpůsobit.