Context Navigation

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

Revision 1991, 12.5 KB checked in by bittner, 17 years ago (diff)

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

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

Download in other formats: