Context Navigation

source: GTP/trunk/Lib/Vis/Preprocessing/src/Mutation.cpp @ 1997

Revision 1997, 14.5 KB checked in by bittner, 17 years ago (diff)
sil mutation

Line
1	#include "SamplingStrategy.h"
2	#include "Ray.h"
3	#include "Intersectable.h"
4	#include "Preprocessor.h"
5	#include "ViewCellsManager.h"
6	#include "AxisAlignedBox3.h"
7	#include "RssTree.h"
8	#include "Vector2.h"
9	#include "RndGauss.h"
10	#include "Mutation.h"
11
12	#ifdef GTP_INTERNAL
13	#include "ArchModeler2MLRT.hxx"
14	#endif
15
16	namespace GtpVisibilityPreprocessor {
17
18	#define MUTATION_USE_CDF 0
19	#define USE_SILHOUETTE_MUTATIONS 0
20
21	#define USE_SIL_TERMINATION_MUTATION 1
22
23	#define EVALUATE_MUTATION_STATS 1
24
25	void
26	MutationBasedDistribution::Update(VssRayContainer &vssRays)
27	{
28	// for (int i=0; i < mRays.size(); i++)
29	// cout<<mRays[i].mMutations<<" ";
30	// cout<<endl;
31	cerr<<"Muattion update..."<<endl;
32	cerr<<"rays = "<<mRays.size()<<endl;
33	if (mRays.size()) {
34	cerr<<"Oversampling factors = "<<
35	GetEntry(0).mMutations<<" "<<
36	GetEntry(1).mMutations<<" "<<
37	GetEntry(2).mMutations<<" "<<
38	GetEntry(3).mMutations<<" "<<
39	GetEntry(4).mMutations<<" "<<
40	GetEntry(5).mMutations<<" ... "<<
41	GetEntry(mRays.size()-6).mMutations<<" "<<
42	GetEntry(mRays.size()-5).mMutations<<" "<<
43	GetEntry(mRays.size()-4).mMutations<<" "<<
44	GetEntry(mRays.size()-3).mMutations<<" "<<
45	GetEntry(mRays.size()-2).mMutations<<" "<<
46	GetEntry(mRays.size()-1).mMutations<<endl;
47	}
48	int contributingRays = 0;
49
50	int mutationRays = 0;
51	int dummyNcMutations = 0;
52	int dummyCMutations = 0;
53
54	for (int i=0; i < vssRays.size(); i++) {
55	if (vssRays[i]->mPvsContribution) {
56	// reset the counter of unsuccsseful mutation for a generating ray (if it exists)
57	if (vssRays[i]->mDistribution == MUTATION_BASED_DISTRIBUTION &&
58	vssRays[i]->mGeneratorId != -1
59	) {
60	mRays[vssRays[i]->mGeneratorId].mUnsuccessfulMutations = 0;
61	#if EVALUATE_MUTATION_STATS
62	mutationRays++;
63
64	Intersectable *newObject =
65	mPreprocessor.mViewCellsManager->GetIntersectable(
66	*vssRays[i],
67	true);
68
69	Intersectable *oldObject =
70	mPreprocessor.mViewCellsManager->GetIntersectable(
71	*mRays[vssRays[i]->mGeneratorId].mRay,
72	true);
73
74	if (oldObject == newObject)
75	dummyCMutations++;
76	#endif
77	}
78	contributingRays++;
79	if (mRays.size() < mMaxRays) {
80	VssRay newRay = new VssRay(vssRays[i]);
81	// add this ray
82	newRay->Ref();
83	mRays.push_back(RayEntry(newRay));
84	} else {
85	// unref the old ray
86	mRays[mBufferStart].mRay = vssRays[i];
87	mRays[mBufferStart].mMutations = 0;
88	// mRays[mBufferStart] = RayEntry(newRay);
89	mBufferStart++;
90	if (mBufferStart >= mMaxRays)
91	mBufferStart = 0;
92	}
93	} else {
94	#if EVALUATE_MUTATION_STATS
95	if (vssRays[i]->mDistribution == MUTATION_BASED_DISTRIBUTION &&
96	vssRays[i]->mGeneratorId != -1
97	) {
98	mutationRays++;
99
100	Intersectable *newObject =
101	mPreprocessor.mViewCellsManager->GetIntersectable(
102	*vssRays[i],
103	true);
104
105	Intersectable *oldObject =
106	mPreprocessor.mViewCellsManager->GetIntersectable(
107	*mRays[vssRays[i]->mGeneratorId].mRay,
108	true);
109
110	if (oldObject == newObject)
111	dummyNcMutations++;
112	}
113	#endif
114	}
115	}
116
117	if (mutationRays) {
118	cout<<"Mutated rays:"<<mutationRays<<endl;
119	cout<<"Dummy mutations ratio:"<<100.0f*(dummyCMutations + dummyNcMutations)/
120	(float)mutationRays<<"%"<<endl;
121	cout<<"Dummy NC mutations ratio:"<<100.0f*dummyNcMutations/(float)mutationRays<<"%"<<endl;
122	cout<<"Dummy C mutations ratio:"<<100.0f*dummyCMutations/(float)mutationRays<<"%"<<endl;
123	}
124
125	float pContributingRays = contributingRays/(float)vssRays.size();
126	float importance = 1.0f/(pContributingRays + 1e-5);
127	// set this values for last contributingRays
128	int index = mBufferStart - 1;
129
130	for (int i=0; i < contributingRays; i++, index--) {
131	if (index < 0)
132	index = mRays.size()-1;
133	mRays[index].mImportance = importance;
134	}
135
136	#if MUTATION_USE_CDF
137	// compute cdf
138	mRays[0].mCdf = mRays[0].mImportance/(mRays[0].mMutations+1);
139	for (int i=1; i < mRays.size(); i++)
140	mRays[i].mCdf = mRays[i-1].mCdf + mRays[i].mImportance/(mRays[i].mMutations+1);
141
142	float scale = 1.0f/mRays[i-1].mCdf;
143	for (i=0; i < mRays.size(); i++) {
144	mRays[i].mCdf *= scale;
145	}
146	#endif
147
148	cout<<"Importance = "<<
149	GetEntry(0).mImportance<<" "<<
150	GetEntry(mRays.size()-1).mImportance<<endl;
151
152	cerr<<"Mutation update done."<<endl;
153	}
154
155
156	Vector3
157	MutationBasedDistribution::ComputeOriginMutation(const VssRay &ray,
158	const Vector3 &U,
159	const Vector3 &V,
160	const Vector2 vr2,
161	const float radius
162	)
163	{
164	#if 0
165	Vector3 v;
166	if (d.DrivingAxis() == 0)
167	v = Vector3(0, r[0]-0.5f, r[1]-0.5f);
168	else
169	if (d.DrivingAxis() == 1)
170	v = Vector3(r[0]-0.5f, 0, r[1]-0.5f);
171	else
172	v = Vector3(r[0]-0.5f, r[1]-0.5f, 0);
173	return v(2radius);
174	#endif
175	#if 0
176	return (U(r[0] - 0.5f) + V(r[1] - 0.5f))(2radius);
177	#endif
178
179
180	// Output random variable
181	Vector2 gaussvec2;
182
183	// Here we apply transform to gaussian, so 2D bivariate
184	// normal distribution
185	// float sigma = ComputeSigmaFromRadius(radius);
186	float sigma = radius;
187	GaussianOn2D(vr2,
188	sigma, // input
189	gaussvec2); // output
190
191
192	// Here we tranform the point correctly to 3D space using base
193	// vectors of the 3D space defined by the direction
194	Vector3 shift = gaussvec2.xx * U + gaussvec2.yy * V;
195
196	// cout<<shift<<endl;
197	return shift;
198	}
199
200	Vector3
201	MutationBasedDistribution::ComputeTerminationMutation(const VssRay &ray,
202	const Vector3 &U,
203	const Vector3 &V,
204	const Vector2 vr2,
205	const float radius
206	)
207	{
208	#if 0
209	Vector3 v;
210	// mutate the termination
211	if (d.DrivingAxis() == 0)
212	v = Vector3(0, r[2]-0.5f, r[3]-0.5f);
213	else
214	if (d.DrivingAxis() == 1)
215	v = Vector3(r[2]-0.5f, 0, r[3]-0.5f);
216	else
217	v = Vector3(r[2]-0.5f, r[3]-0.5f, 0);
218
219	// Vector3 nv;
220
221	// if (Magnitude(v) > Limits::Small)
222	// nv = Normalize(v);
223	// else
224	// nv = v;
225
226	// v = nvsize + vsize;
227
228	return v(4.0fradius);
229	#endif
230	#if 0
231	return (U(vr2.xx - 0.5f) + V(vr2.yy - 0.5f))(4.0fradius);
232	#endif
233	Vector2 gaussvec2;
234	#if 1
235	float sigma = radius;
236	GaussianOn2D(vr2,
237	sigma, // input
238	gaussvec2); // output
239	Vector3 shift = gaussvec2.xx * U + gaussvec2.yy * V;
240	// cout<<shift<<endl;
241	return shift;
242	#endif
243	#if 0
244	// Here we estimate standard deviation (sigma) from radius
245	float sigma = 1.1f*ComputeSigmaFromRadius(radius);
246	Vector3 vr3(vr2.xx, vr2.yy, RandomValue(0,1));
247	PolarGaussianOnDisk(vr3,
248	sigma,
249	radius, // input
250	gaussvec2); // output
251
252	// Here we tranform the point correctly to 3D space using base
253	// vectors of the 3D space defined by the direction
254	Vector3 shift = gaussvec2.xx * U + gaussvec2.yy * V;
255
256	// cout<<shift<<endl;
257	return shift;
258	#endif
259	}
260
261
262	Vector3
263	MutationBasedDistribution::ComputeSilhouetteTerminationMutation(const VssRay &ray,
264	const Vector3 &origin,
265	const AxisAlignedBox3 &box,
266	const Vector3 &U,
267	const Vector3 &V,
268	const Vector2 vr2,
269	const float radius
270	)
271	{
272	const int packetSize = 16;
273
274	static int hit_triangles[packetSize];
275	static float dist[packetSize];
276	static Vector3 dirs[packetSize];
277	// mutate the
278	float alpha = RandomValue(0.0f, 2.0f*M_PI);
279
280	// direction along which we will mutate the ray
281	Vector3 line = sin(alpha)U + cos(alpha)V;
282
283	// cout<<line<<endl;
284	// create 16 rays along the selected dir
285
286	// cast rays to find silhouette ray
287	for (int i=0; i < packetSize; i++) {
288	dirs[i] = Normalize(ray.mTermination + ((radius/(packetSize - i))*line) - origin );
289	mlrtaStoreRayAS16(&origin.x,
290	&dirs[i].x,
291	i);
292	}
293
294	mlrtaTraverseGroupAS16(&box.Min().x,
295	&box.Max().x,
296	hit_triangles,
297	dist);
298
299	for (int i=0; i < packetSize; i++) {
300	if (hit_triangles[i] == -1 \|\| !box.IsInside(origin + dist[i]*dirs[i])) {
301	// break on first passing ray
302	break;
303	}
304	}
305
306	if (i == packetSize) {
307	// cerr<<"Warning: hit the same box here should never happen!"<<endl;
308	// shift the ray even a bit more
309	// cout<<"W"<<i<<endl;
310	return ray.mTermination + (RandomValue(1.0f, 2.0f)radius)line;
311	}
312
313	// cout<<i<<endl;
314	return dirs[i];
315	}
316
317
318	bool
319	MutationBasedDistribution::GenerateSample(SimpleRay &sray)
320	{
321
322	if (mRays.size() == 0) {
323	float rr[5];
324	// use direction based distribution
325	Vector3 origin, direction;
326	static HaltonSequence halton;
327
328	halton.GetNext(5, rr);
329	mPreprocessor.mViewCellsManager->GetViewPoint(origin,
330	Vector3(rr[0], rr[1], rr[2]));
331
332
333	direction = UniformRandomVector(rr[3], rr[4]);
334
335	const float pdf = 1.0f;
336	sray = SimpleRay(origin, direction, MUTATION_BASED_DISTRIBUTION, pdf);
337	sray.mGeneratorId = -1;
338
339	return true;
340	}
341
342	int index;
343
344	#if !MUTATION_USE_CDF
345	// get tail of the buffer
346	index = (mLastIndex+1)%mRays.size();
347	if (mRays[index].GetSamplingFactor() >
348	mRays[mLastIndex].GetSamplingFactor()) {
349	// search back for index where this is valid
350	index = (mLastIndex - 1 + mRays.size())%mRays.size();
351	for (int i=0; i < mRays.size(); i++) {
352
353	// if (mRays[index].mMutations > mRays[mLastIndex].mMutations)
354	// break;
355	if (mRays[index].GetSamplingFactor() >
356	mRays[mLastIndex].GetSamplingFactor() )
357	break;
358	index = (index - 1 + mRays.size())%mRays.size();
359	}
360	// go one step back
361	index = (index+1)%mRays.size();
362	}
363	#else
364	static HaltonSequence iHalton;
365	iHalton.GetNext(1, rr);
366	//rr[0] = RandomValue(0,1);
367	// use binary search to find index with this cdf
368	int l=0, r=mRays.size()-1;
369	while(l<r) {
370	int i = (l+r)/2;
371	if (rr[0] < mRays[i].mCdf )
372	r = i;
373	else
374	l = i+1;
375	}
376	index = l;
377	// if (rr[0] >= mRays[r].mCdf)
378	// index = r;
379	// else
380	// index = l;
381
382
383	#endif
384	// cout<<index<<" "<<rr[0]<<" "<<mRays[index].mCdf<<" "<<mRays[(index+1)%mRays.size()].mCdf<<endl;
385
386	mLastIndex = index;
387
388	#if USE_SILHOUETTE_MUTATIONS
389	return GenerateSilhouetteMutation(index, sray);
390	#else
391	return GenerateMutation(index, sray);
392	#endif
393
394	}
395
396
397
398
399
400	bool
401	MutationBasedDistribution::GenerateMutationCandidate(const int index,
402	SimpleRay &sray,
403	Intersectable *object,
404	const AxisAlignedBox3 &box
405	)
406	{
407	float rr[4];
408
409	VssRay *ray = mRays[index].mRay;
410
411	mRays[index].mHalton.GetNext(4, rr);
412
413	// mutate the origin
414	Vector3 d = ray->GetDir();
415
416	float objectRadius = 0.5f*Magnitude(box.Diagonal());
417	// cout<<objectRadius<<endl;
418	if (objectRadius < Limits::Small)
419	return false;
420
421	// Compute right handed coordinate system from direction
422	Vector3 U, V;
423	Vector3 nd = Normalize(d);
424	nd.RightHandedBase(U, V);
425
426	Vector3 origin = ray->mOrigin;
427	Vector3 termination = ray->mTermination; //box.Center(); //ray->mTermination; //box.Center();
428
429	// optimal for Pompeii 0.1f
430	// optimal for Vienna 0.5f
431
432	float radiusExtension = 0.05f;
433	// + mRays[index].mMutations/50.0f;
434
435	float mutationRadius = objectRadius*radiusExtension;
436
437	// tmp for pompeii
438	mutationRadius = 0.22f;
439
440	origin += ComputeOriginMutation(*ray, U, V,
441	Vector2(rr[0], rr[1]),
442	mutationRadius);
443
444	#if USE_SIL_TERMINATION_MUTATION
445	termination += ComputeSilhouetteTerminationMutation(*ray,
446	origin,
447	box,
448	U, V,
449	Vector2(rr[2], rr[3]),
450	3.0f*objectRadius);
451	#else
452	termination += ComputeTerminationMutation(*ray, U, V,
453	Vector2(rr[2], rr[3]),
454	mutationRadius);
455	#endif
456	Vector3 direction = termination - origin;
457
458	if (Magnitude(direction) < Limits::Small)
459	return false;
460
461	// shift the origin a little bit
462	origin += direction*0.5f;
463
464	direction.Normalize();
465
466	// $$ jb the pdf is yet not correct for all sampling methods!
467	const float pdf = 1.0f;
468
469	sray = SimpleRay(origin, direction, MUTATION_BASED_DISTRIBUTION, pdf);
470	sray.mGeneratorId = index;
471	}
472
473	bool
474	MutationBasedDistribution::GenerateMutation(const int index, SimpleRay &sray)
475	{
476	VssRay *ray = mRays[index].mRay;
477
478	Intersectable *object = mPreprocessor.mViewCellsManager->GetIntersectable(
479	*ray,
480	true);
481
482	AxisAlignedBox3 box = object->GetBox();
483
484	if (GenerateMutationCandidate(index, sray, object, box)) {
485	mRays[index].mMutations++;
486	return true;
487	}
488	return false;
489	}
490
491	bool
492	MutationBasedDistribution::GenerateSilhouetteMutation(const int index, SimpleRay &sray)
493	{
494	#ifndef GTP_INTERNAL
495	return GenerateMutation(index, sray);
496	#else
497	const int packetSize = 4;
498	const int maxTries = 8;
499
500	static int hit_triangles[16];
501	static float dist[16];
502
503	SimpleRay mutationCandidates[packetSize];
504	int candidates = 0;
505
506	VssRay *ray = mRays[index].mRay;
507
508	Intersectable *object = mPreprocessor.mViewCellsManager->GetIntersectable(
509	*ray,
510	true);
511
512	AxisAlignedBox3 box = object->GetBox();
513
514	int id = 0;
515	int silhouetteRays = 0;
516	int tries = 0;
517	while (silhouetteRays == 0 && tries < maxTries) {
518	for (candidates = 0; candidates < packetSize && tries < maxTries; tries++)
519	if (GenerateMutationCandidate(index, mutationCandidates[candidates], object, box))
520	candidates++;
521
522	if (candidates < packetSize)
523	break;
524
525	// cout<<candidates<<endl;
526	// cast rays to find silhouette edge
527	for (int i=0; i < packetSize; i++)
528	mlrtaStoreRayAS4(&mutationCandidates[i].mOrigin.x,
529	&mutationCandidates[i].mDirection.x,
530	i);
531
532	mlrtaTraverseGroupAS4(&box.Min().x,
533	&box.Max().x,
534	hit_triangles,
535	dist);
536
537	for (int i=0; i < packetSize; i++)
538	if (hit_triangles[i] == -1) {
539	silhouetteRays++;
540	id = i;
541	break;
542	}
543	}
544
545	if (candidates == 0)
546	return false;
547
548	// cout<<id<<endl;
549	// cout<<tries<<endl;
550	sray = mutationCandidates[id];
551	mRays[index].mMutations++;
552
553	return true;
554	#endif
555	}
556
557
558
559
560	MutationBasedDistribution::MutationBasedDistribution(Preprocessor &preprocessor
561	) :
562	SamplingStrategy(preprocessor)
563	{
564	mType = MUTATION_BASED_DISTRIBUTION;
565	mBufferStart = 0;
566	mMaxRays = 500000;
567	mRays.reserve(mMaxRays);
568	mOriginMutationSize = 10.0f;
569	mLastIndex = 0;
570	// mOriginMutationSize = Magnitude(preprocessor.mViewCellsManager->
571	// GetViewSpaceBox().Diagonal())*1e-3;
572
573	}
574
575
576	}

Note: See TracBrowser for help on using the repository browser.

Download in other formats: