Context Navigation

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

Revision 2011, 21.9 KB checked in by bittner, 18 years ago (diff)
mutation strategy updated, importance depending on rsuccess

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	#include "Exporter.h"
12
13	#ifdef GTP_INTERNAL
14	#include "ArchModeler2MLRT.hxx"
15	#endif
16
17	namespace GtpVisibilityPreprocessor {
18
19	#define MUTATION_USE_CDF 0
20	#define USE_SILHOUETTE_MUTATIONS 0
21
22	#define USE_SIL_TERMINATION_MUTATION 1
23	#define MUTATE_ORIGIN 0
24
25	#define EVALUATE_MUTATION_STATS 1
26
27	#define Q_SEARCH_STEPS 3
28
29	#define SORT_RAY_ENTRIES 1
30
31	// use avg ray contribution as importance
32	// if 0 the importance is evaluated from the succ of mutations
33	#define USE_AVG_CONTRIBUTION 1
34
35	MutationBasedDistribution::RayEntry &
36	MutationBasedDistribution::GetEntry(const int index)
37	{
38	#if SORT_RAY_ENTRIES
39	return mRays[index];
40	#else
41	return mRays[(mBufferStart+index)%mRays.size()];
42	#endif
43	}
44
45	void
46	MutationBasedDistribution::Update(VssRayContainer &vssRays)
47	{
48	// for (int i=0; i < mRays.size(); i++)
49	// cout<<mRays[i].mMutations<<" ";
50	// cout<<endl;
51	cerr<<"Muattion update..."<<endl;
52	cerr<<"rays = "<<mRays.size()<<endl;
53	if (mRays.size()) {
54	cerr<<"Oversampling factors = "<<
55	GetEntry(0).mMutations<<" "<<
56	GetEntry(1).mMutations<<" "<<
57	GetEntry(2).mMutations<<" "<<
58	GetEntry(3).mMutations<<" "<<
59	GetEntry(4).mMutations<<" "<<
60	GetEntry(5).mMutations<<" ... "<<
61	GetEntry(mRays.size()-6).mMutations<<" "<<
62	GetEntry(mRays.size()-5).mMutations<<" "<<
63	GetEntry(mRays.size()-4).mMutations<<" "<<
64	GetEntry(mRays.size()-3).mMutations<<" "<<
65	GetEntry(mRays.size()-2).mMutations<<" "<<
66	GetEntry(mRays.size()-1).mMutations<<endl;
67	}
68	int contributingRays = 0;
69
70	int mutationRays = 0;
71	int dummyNcMutations = 0;
72	int dummyCMutations = 0;
73
74	int reverseCandidates = 0;
75
76	#if 0
77	sort(mRays.begin(), mRays.end());
78	// reset the start of the buffer
79	mBufferStart = 0;
80	#endif
81
82	for (int i=0; i < vssRays.size(); i++) {
83	if (vssRays[i]->mPvsContribution) {
84	// reset the counter of unsuccsseful mutation for a generating ray (if it exists)
85	if (vssRays[i]->mDistribution == MUTATION_BASED_DISTRIBUTION &&
86	vssRays[i]->mGeneratorId != -1
87	) {
88	mRays[vssRays[i]->mGeneratorId].mUnsuccessfulMutations = 0;
89	#if EVALUATE_MUTATION_STATS
90	mutationRays++;
91
92	Intersectable *newObject =
93	mPreprocessor.mViewCellsManager->GetIntersectable(
94	*vssRays[i],
95	true);
96
97	Intersectable *oldObject =
98	mPreprocessor.mViewCellsManager->GetIntersectable(
99	*mRays[vssRays[i]->mGeneratorId].mRay,
100	true);
101
102	if (oldObject == newObject)
103	dummyCMutations++;
104	#endif
105	}
106	contributingRays++;
107	if (mRays.size() < mMaxRays) {
108	VssRay newRay = new VssRay(vssRays[i]);
109	// add this ray
110	newRay->Ref();
111	mRays.push_back(RayEntry(newRay));
112	} else {
113	// unref the old ray
114	mRays[mBufferStart].mRay = vssRays[i];
115	mRays[mBufferStart].mMutations = 0;
116	mRays[mBufferStart].mUnsuccessfulMutations = 0;
117	mRays[mBufferStart].ResetReverseMutation();
118	// mRays[mBufferStart] = RayEntry(newRay);
119	mBufferStart++;
120	if (mBufferStart >= mMaxRays)
121	mBufferStart = 0;
122	}
123	} else {
124	if (vssRays[i]->mDistribution == MUTATION_BASED_DISTRIBUTION &&
125	vssRays[i]->mGeneratorId != -1
126	) {
127	// check whether not to store a new backward mutation candidate
128	VssRay *oldRay = mRays[vssRays[i]->mGeneratorId].mRay;
129	VssRay *newRay = vssRays[i];
130
131	#define DIST_THRESHOLD 3.0f
132
133	Intersectable *oldObject =
134	mPreprocessor.mViewCellsManager->GetIntersectable(
135	*oldRay,
136	true);
137
138
139	if (!mRays[newRay->mGeneratorId].HasReverseMutation()) {
140	if (DotProd(oldRay->GetDir(), newRay->GetDir()) > 0.0f) {
141	float oldDist = Magnitude(oldRay->mTermination - newRay->mOrigin);
142	float newDist = Magnitude(newRay->mTermination - newRay->mOrigin);
143
144	if (newDist < oldDist - oldObject->GetBox().Radius()*DIST_THRESHOLD) {
145	Vector3 origin, termination;
146	if (ComputeReverseMutation(oldRay, newRay, origin, termination)) {
147	mRays[newRay->mGeneratorId].SetReverseMutation(origin, termination);
148	}
149
150	reverseCandidates++;
151	//mReverseCandidates
152	}
153	}
154	}
155	#if EVALUATE_MUTATION_STATS
156	mutationRays++;
157
158	Intersectable *newObject =
159	mPreprocessor.mViewCellsManager->GetIntersectable(
160	*vssRays[i],
161	true);
162
163
164	if (oldObject == newObject)
165	dummyNcMutations++;
166	#endif
167	}
168	}
169	}
170
171	if (mutationRays) {
172	cout<<"Mutated rays:"<<mutationRays<<endl;
173	cout<<"Dummy mutations ratio:"<<100.0f*(dummyCMutations + dummyNcMutations)/
174	(float)mutationRays<<"%"<<endl;
175	cout<<"Dummy NC mutations ratio:"<<100.0f*dummyNcMutations/(float)mutationRays<<"%"<<endl;
176	cout<<"Dummy C mutations ratio:"<<100.0f*dummyCMutations/(float)mutationRays<<"%"<<endl;
177	cout<<"Reverse candidates:"<<reverseCandidates<<endl;
178	}
179
180	float pContributingRays = contributingRays/(float)vssRays.size();
181
182	cout<<"Percentage of contributing rays:"<<pContributingRays<<endl;
183
184	#if USE_AVG_CONTRIBUTION
185	float importance = 1.0f/(pContributingRays + 1e-5);
186	// float importance = 1.0f;
187	// set this values for last contributingRays
188	int index = mBufferStart - 1;
189
190	for (int i=0; i < contributingRays; i++, index--) {
191	if (index < 0)
192	index = mRays.size()-1;
193	mRays[index].mImportance = importance;
194	}
195	#else
196	// use unsucc mutation samples as feedback on importance
197	for (int i=0; i < mRays.size(); i++) {
198	const float minImportance = 0.1f;
199	const int minImportanceSamples = 20;
200	mRays[i].mImportance = minImportance +
201	(1-minImportance)exp(-3.0fmRays[i].mUnsuccessfulMutations/minImportanceSamples);
202
203	//mRays[i].mImportance = 1.0f/(mRays[i].mUnsuccessfulMutations+3);
204	// mRays[i].mImportance = 1.0f;
205	}
206	#endif
207
208	#if SORT_RAY_ENTRIES
209	long t1 = GetTime();
210	sort(mRays.begin(), mRays.end());
211	// reset the start of the buffer
212	mBufferStart = 0;
213	mLastIndex = mRays.size();
214	cout<<"Mutation candidates sorted in "<<TimeDiff(t1, GetTime())<<" ms."<<endl;
215	#endif
216
217	#if MUTATION_USE_CDF
218	// compute cdf
219	mRays[0].mCdf = mRays[0].mImportance/(mRays[0].mMutations+1);
220	for (int i=1; i < mRays.size(); i++)
221	mRays[i].mCdf = mRays[i-1].mCdf + mRays[i].mImportance/(mRays[i].mMutations+1);
222
223	float scale = 1.0f/mRays[i-1].mCdf;
224	for (i=0; i < mRays.size(); i++) {
225	mRays[i].mCdf *= scale;
226	}
227	#endif
228
229	cout<<"Importance = "<<
230	GetEntry(0).mImportance<<" "<<
231	GetEntry(mRays.size()-1).mImportance<<endl;
232
233	cout<<"Sampling factor = "<<
234	GetEntry(0).GetSamplingFactor()<<" "<<
235	GetEntry(mRays.size()-1).GetSamplingFactor()<<endl;
236
237	cerr<<"Mutation update done."<<endl;
238	}
239
240
241	Vector3
242	MutationBasedDistribution::ComputeOriginMutation(const VssRay &ray,
243	const Vector3 &U,
244	const Vector3 &V,
245	const Vector2 vr2,
246	const float radius
247	)
248	{
249	#if 0
250	Vector3 v;
251	if (d.DrivingAxis() == 0)
252	v = Vector3(0, r[0]-0.5f, r[1]-0.5f);
253	else
254	if (d.DrivingAxis() == 1)
255	v = Vector3(r[0]-0.5f, 0, r[1]-0.5f);
256	else
257	v = Vector3(r[0]-0.5f, r[1]-0.5f, 0);
258	return v(2radius);
259	#endif
260	#if 0
261	return (U(r[0] - 0.5f) + V(r[1] - 0.5f))(2radius);
262	#endif
263
264
265	// Output random variable
266	Vector2 gaussvec2;
267
268	// Here we apply transform to gaussian, so 2D bivariate
269	// normal distribution
270	// float sigma = ComputeSigmaFromRadius(radius);
271	float sigma = radius;
272	GaussianOn2D(vr2,
273	sigma, // input
274	gaussvec2); // output
275
276
277	// Here we tranform the point correctly to 3D space using base
278	// vectors of the 3D space defined by the direction
279	Vector3 shift = gaussvec2.xx * U + gaussvec2.yy * V;
280
281	// cout<<shift<<endl;
282	return shift;
283	}
284
285	Vector3
286	MutationBasedDistribution::ComputeTerminationMutation(const VssRay &ray,
287	const Vector3 &U,
288	const Vector3 &V,
289	const Vector2 vr2,
290	const float radius
291	)
292	{
293	#if 0
294	Vector3 v;
295	// mutate the termination
296	if (d.DrivingAxis() == 0)
297	v = Vector3(0, r[2]-0.5f, r[3]-0.5f);
298	else
299	if (d.DrivingAxis() == 1)
300	v = Vector3(r[2]-0.5f, 0, r[3]-0.5f);
301	else
302	v = Vector3(r[2]-0.5f, r[3]-0.5f, 0);
303
304	// Vector3 nv;
305
306	// if (Magnitude(v) > Limits::Small)
307	// nv = Normalize(v);
308	// else
309	// nv = v;
310
311	// v = nvsize + vsize;
312
313	return v(4.0fradius);
314	#endif
315	#if 0
316	return (U(vr2.xx - 0.5f) + V(vr2.yy - 0.5f))(4.0fradius);
317	#endif
318	Vector2 gaussvec2;
319	#if 1
320	float sigma = radius;
321	GaussianOn2D(vr2,
322	sigma, // input
323	gaussvec2); // output
324	Vector3 shift = gaussvec2.xx * U + gaussvec2.yy * V;
325	// cout<<shift<<endl;
326	return shift;
327	#endif
328	#if 0
329	// Here we estimate standard deviation (sigma) from radius
330	float sigma = 1.1f*ComputeSigmaFromRadius(radius);
331	Vector3 vr3(vr2.xx, vr2.yy, RandomValue(0,1));
332	PolarGaussianOnDisk(vr3,
333	sigma,
334	radius, // input
335	gaussvec2); // output
336
337	// Here we tranform the point correctly to 3D space using base
338	// vectors of the 3D space defined by the direction
339	Vector3 shift = gaussvec2.xx * U + gaussvec2.yy * V;
340
341	// cout<<shift<<endl;
342	return shift;
343	#endif
344	}
345
346	bool
347	MutationBasedDistribution::ComputeReverseMutation(
348	const VssRay &oldRay,
349	const VssRay &newRay,
350	Vector3 &origin,
351	Vector3 &termination
352	)
353	{
354	// first reconstruct the termination point
355	Vector3 oldDir = Normalize(oldRay.GetDir());
356	Plane3 oldPlane(oldDir, oldRay.mTermination);
357
358	termination = oldPlane.FindIntersection(newRay.mOrigin,
359	newRay.mTermination);
360
361	// now find the new origin of the ray by casting ray backward from the termination and termining
362	// silhouette point with respect to the occluding object (object containing the newRay termination)
363
364	Plane3 newPlane(oldDir, newRay.mTermination);
365
366	Vector3 oldPivot = newPlane.FindIntersection(oldRay.mOrigin,
367	oldRay.mTermination);
368
369	Vector3 newPivot = newRay.mTermination;
370	Vector3 line = 2.0f*(oldPivot - newPivot);
371
372	Intersectable *occluder = mPreprocessor.mViewCellsManager->GetIntersectable(
373	newRay,
374	true);
375
376	AxisAlignedBox3 box = occluder->GetBox();
377	box.Scale(2.0f);
378
379	const int packetSize = 4;
380	static int hit_triangles[packetSize];
381	static float dist[packetSize];
382	static Vector3 dirs[packetSize];
383	static Vector3 shifts[packetSize];
384	// now find the silhouette along the line
385	int i;
386	float left = 0.0f;
387	float right = 1.0f;
388	// cast rays to find silhouette ray
389	for (int j=0; j < Q_SEARCH_STEPS; j++) {
390	for (i=0; i < packetSize; i++) {
391	float r = left + (i+1)*(right-left)/(packetSize+1);
392	shifts[i] = r*line;
393	dirs[i] = Normalize(newPivot + shifts[i] - termination );
394	mlrtaStoreRayASEye4(&termination.x,
395	&dirs[i].x,
396	i);
397	}
398
399	mlrtaTraverseGroupASEye4(&box.Min().x,
400	&box.Max().x,
401	hit_triangles,
402	dist);
403
404	for (i=0; i < packetSize; i++) {
405	if (hit_triangles[i] == -1) {
406	// break on first passing ray
407	break;
408	}
409	}
410	float rr = left + (i+1)*(right-left)/(packetSize+1);
411	float rl = left + i*(right-left)/(packetSize+1);
412	left = rl;
413	right = rr;
414	}
415
416	float t = right;
417	if (right==1.0f)
418	return false;
419
420	if (i == packetSize)
421	origin = newPivot + right*line;
422	else
423	origin = newPivot + shifts[i];
424
425	if (0) {
426
427	static VssRayContainer rRays;
428	static int counter = 0;
429	char filename[256];
430
431	if (counter < 50) {
432	sprintf(filename, "reverse_rays_%03d.x3d", counter++);
433
434	VssRay tRay(origin, termination, NULL, NULL);
435	rRays.push_back((VssRay *)&oldRay);
436	rRays.push_back((VssRay *)&newRay);
437	rRays.push_back(&tRay);
438
439	Exporter *exporter = NULL;
440	exporter = Exporter::GetExporter(filename);
441
442	exporter->SetFilled();
443
444	Intersectable *occludee = mPreprocessor.mViewCellsManager->GetIntersectable(
445	oldRay,
446	true);
447
448	exporter->SetForcedMaterial(RgbColor(0,0,1));
449	exporter->ExportIntersectable(occluder);
450	exporter->SetForcedMaterial(RgbColor(0,1,0));
451	exporter->ExportIntersectable(occludee);
452	exporter->ResetForcedMaterial();
453
454	exporter->SetWireframe();
455
456
457	exporter->ExportRays(rRays, RgbColor(1, 0, 0));
458	delete exporter;
459	rRays.clear();
460	}
461	}
462
463
464
465	return true;
466
467	// now the origin and termination is swapped compred to the generator ray
468	// swap(origin, termination);???
469	// -> perhaps not neccessary as the reverse mutation wil only be used once!
470	}
471
472	Vector3
473	MutationBasedDistribution::ComputeSilhouetteTerminationMutation(const VssRay &ray,
474	const Vector3 &origin,
475	const AxisAlignedBox3 &box,
476	const Vector3 &U,
477	const Vector3 &V,
478	const Vector2 vr2,
479	const float radius
480	)
481	{
482	const int packetSize = 4;
483	static int hit_triangles[packetSize];
484	static float dist[packetSize];
485	static Vector3 dirs[packetSize];
486	static Vector3 shifts[packetSize];
487	// mutate the
488	float alpha = RandomValue(0.0f, 2.0f*M_PI);
489	//float alpha = vr2.x2.0fM_PI;
490
491	// direction along which we will mutate the ray
492	Vector3 line = sin(alpha)U + cos(alpha)V;
493
494	// cout<<line<<endl;
495	// create 16 rays along the selected dir
496	int i;
497	float left = 0.0f;
498	float right = radius;
499	// cast rays to find silhouette ray
500	for (int j=0; j < Q_SEARCH_STEPS; j++) {
501	for (i=0; i < packetSize; i++) {
502	float r = left + (i+1)*(right-left)/(packetSize+1);
503	shifts[i] = r*line;
504	dirs[i] = Normalize(ray.mTermination + shifts[i] - origin );
505	mlrtaStoreRayASEye4(&origin.x,
506	&dirs[i].x,
507	i);
508	}
509
510	mlrtaTraverseGroupASEye4(&box.Min().x,
511	&box.Max().x,
512	hit_triangles,
513	dist);
514
515	for (i=0; i < packetSize; i++) {
516	if (hit_triangles[i] == -1) {
517	// if (hit_triangles[i] == -1 \|\| !box.IsInside(origin + dist[i]*dirs[i])) {
518	// break on first passing ray
519	break;
520	}
521	}
522	float rr = left + (i+1)*(right-left)/(packetSize+1);
523	float rl = left + i*(right-left)/(packetSize+1);
524	left = rl;
525	right = rr;
526	}
527
528	if (i == packetSize) {
529	// cerr<<"Warning: hit the same box here should never happen!"<<endl;
530	// shift the ray even a bit more
531	//cout<<"W"<<i<<endl;
532	// return (RandomValue(1.0f, 1.5f)radius)line;
533	return right*line;
534	}
535
536	// cout<<i<<endl;
537	return shifts[i];
538	}
539
540
541	bool
542	MutationBasedDistribution::GenerateSample(SimpleRay &sray)
543	{
544
545	if (mRays.size() == 0) {
546	float rr[5];
547	// use direction based distribution
548	Vector3 origin, direction;
549	static HaltonSequence halton;
550
551	halton.GetNext(5, rr);
552	mPreprocessor.mViewCellsManager->GetViewPoint(origin,
553	Vector3(rr[0], rr[1], rr[2]));
554
555
556	direction = UniformRandomVector(rr[3], rr[4]);
557
558	const float pdf = 1.0f;
559	sray = SimpleRay(origin, direction, MUTATION_BASED_DISTRIBUTION, pdf);
560	sray.mGeneratorId = -1;
561
562	return true;
563	}
564
565	int index;
566
567	#if !MUTATION_USE_CDF
568	#if SORT_RAY_ENTRIES
569	index = mLastIndex - 1;
570	if (index < 0 \|\| index >= mRays.size()-1) {
571	index = mRays.size() - 1;
572	} else
573	if (
574	mRays[index].GetSamplingFactor() >= mRays[mLastIndex].GetSamplingFactor()) {
575	// make another round
576
577	// cout<<"R2"<<endl;
578	// cout<<mLastIndex<<endl;
579	// cout<<index<<endl;
580	index = mRays.size() - 1;
581	}
582	#else
583	// get tail of the buffer
584	index = (mLastIndex+1)%mRays.size();
585	if (mRays[index].GetSamplingFactor() >
586	mRays[mLastIndex].GetSamplingFactor()) {
587	// search back for index where this is valid
588	index = (mLastIndex - 1 + mRays.size())%mRays.size();
589	for (int i=0; i < mRays.size(); i++) {
590
591	// if (mRays[index].mMutations > mRays[mLastIndex].mMutations)
592	// break;
593	if (mRays[index].GetSamplingFactor() >
594	mRays[mLastIndex].GetSamplingFactor() )
595	break;
596	index = (index - 1 + mRays.size())%mRays.size();
597	}
598	// go one step back
599	index = (index+1)%mRays.size();
600	}
601	#endif
602	#else
603	static HaltonSequence iHalton;
604	iHalton.GetNext(1, rr);
605	//rr[0] = RandomValue(0,1);
606	// use binary search to find index with this cdf
607	int l=0, r=mRays.size()-1;
608	while(l<r) {
609	int i = (l+r)/2;
610	if (rr[0] < mRays[i].mCdf )
611	r = i;
612	else
613	l = i+1;
614	}
615	index = l;
616	// if (rr[0] >= mRays[r].mCdf)
617	// index = r;
618	// else
619	// index = l;
620
621
622	#endif
623	// cout<<index<<" "<<rr[0]<<" "<<mRays[index].mCdf<<" "<<mRays[(index+1)%mRays.size()].mCdf<<endl;
624
625	mLastIndex = index;
626	// Debug<<index<<" "<<mRays[index].GetSamplingFactor()<<endl;
627
628	if (mRays[index].HasReverseMutation()) {
629	//cout<<"R "<<mRays[index].mutatedOrigin<<" "<<mRays[index].mutatedTermination<<endl;
630	sray = SimpleRay(mRays[index].mutatedOrigin,
631	Normalize(mRays[index].mutatedTermination - mRays[index].mutatedOrigin),
632	MUTATION_BASED_DISTRIBUTION, 1.0f);
633	sray.mGeneratorId = index;
634	mRays[index].ResetReverseMutation();
635	mRays[index].mMutations++;
636	mRays[index].mUnsuccessfulMutations++;
637
638	return true;
639	}
640
641	#if USE_SILHOUETTE_MUTATIONS
642	return GenerateSilhouetteMutation(index, sray);
643	#else
644	return GenerateMutation(index, sray);
645	#endif
646	}
647
648
649
650
651
652	bool
653	MutationBasedDistribution::GenerateMutationCandidate(const int index,
654	SimpleRay &sray,
655	Intersectable *object,
656	const AxisAlignedBox3 &box
657	)
658	{
659	float rr[4];
660
661	VssRay *ray = mRays[index].mRay;
662
663	mRays[index].mHalton.GetNext(4, rr);
664	// rr[0] = RandomValue(0.0f,0.99999f);
665	// rr[1] = RandomValue(0.0f,0.99999f);
666	// rr[2] = RandomValue(0.0f,0.99999f);
667	// rr[3] = RandomValue(0.0f,0.99999f);
668
669	// mutate the origin
670	Vector3 d = ray->GetDir();
671
672	float objectRadius = 0.5f*Magnitude(box.Diagonal());
673	// cout<<objectRadius<<endl;
674	if (objectRadius < Limits::Small)
675	return false;
676
677	// Compute right handed coordinate system from direction
678	Vector3 U, V;
679	Vector3 nd = Normalize(d);
680	nd.RightHandedBase(U, V);
681
682	Vector3 origin = ray->mOrigin;
683	Vector3 termination = ray->mTermination; //box.Center(); //ray->mTermination; //box.Center();
684
685	// optimal for Pompeii 0.1f
686	// optimal for Vienna 0.5f
687
688	float radiusExtension = 0.1f;
689	// + mRays[index].mMutations/50.0f;
690
691	float mutationRadius = objectRadius*radiusExtension;
692
693	// tmp for pompeii
694	// mutationRadius = 0.22f;
695
696	#if MUTATE_ORIGIN
697	origin += ComputeOriginMutation(*ray, U, V,
698	Vector2(rr[0], rr[1]),
699	mutationRadius);
700	#endif
701
702	#if USE_SIL_TERMINATION_MUTATION
703	termination += ComputeSilhouetteTerminationMutation(*ray,
704	origin,
705	box,
706	U, V,
707	Vector2(rr[2], rr[3]),
708	2.0f*objectRadius);
709	#else
710	termination += ComputeTerminationMutation(*ray, U, V,
711	Vector2(rr[2], rr[3]),
712	mutationRadius);
713	#endif
714	Vector3 direction = termination - origin;
715
716	if (Magnitude(direction) < Limits::Small)
717	return false;
718
719	// shift the origin a little bit
720	origin += direction*0.5f;
721
722	direction.Normalize();
723
724	// $$ jb the pdf is yet not correct for all sampling methods!
725	const float pdf = 1.0f;
726
727	sray = SimpleRay(origin, direction, MUTATION_BASED_DISTRIBUTION, pdf);
728	sray.mGeneratorId = index;
729	}
730
731	bool
732	MutationBasedDistribution::GenerateMutation(const int index, SimpleRay &sray)
733	{
734	VssRay *ray = mRays[index].mRay;
735
736	Intersectable *object = mPreprocessor.mViewCellsManager->GetIntersectable(
737	*ray,
738	true);
739
740	AxisAlignedBox3 box = object->GetBox();
741
742	if (GenerateMutationCandidate(index, sray, object, box)) {
743	mRays[index].mMutations++;
744	mRays[index].mUnsuccessfulMutations++;
745
746	return true;
747	}
748	return false;
749	}
750
751	bool
752	MutationBasedDistribution::GenerateSilhouetteMutation(const int index, SimpleRay &sray)
753	{
754	#ifndef GTP_INTERNAL
755	return GenerateMutation(index, sray);
756	#else
757	const int packetSize = 4;
758	const int maxTries = 8;
759
760	static int hit_triangles[16];
761	static float dist[16];
762
763	SimpleRay mutationCandidates[packetSize];
764	int candidates = 0;
765
766	VssRay *ray = mRays[index].mRay;
767
768	Intersectable *object = mPreprocessor.mViewCellsManager->GetIntersectable(
769	*ray,
770	true);
771
772	AxisAlignedBox3 box = object->GetBox();
773
774	int id = 0;
775	int silhouetteRays = 0;
776	int tries = 0;
777	while (silhouetteRays == 0 && tries < maxTries) {
778	for (candidates = 0; candidates < packetSize && tries < maxTries; tries++)
779	if (GenerateMutationCandidate(index, mutationCandidates[candidates], object, box))
780	candidates++;
781
782	if (candidates < packetSize)
783	break;
784
785	// cout<<candidates<<endl;
786	// cast rays to find silhouette edge
787	for (int i=0; i < packetSize; i++)
788	mlrtaStoreRayAS4(&mutationCandidates[i].mOrigin.x,
789	&mutationCandidates[i].mDirection.x,
790	i);
791
792	mlrtaTraverseGroupAS4(&box.Min().x,
793	&box.Max().x,
794	hit_triangles,
795	dist);
796
797	for (int i=0; i < packetSize; i++)
798	if (hit_triangles[i] == -1) {
799	silhouetteRays++;
800	id = i;
801	break;
802	}
803	}
804
805	if (candidates == 0)
806	return false;
807
808	// cout<<id<<endl;
809	// cout<<tries<<endl;
810	sray = mutationCandidates[id];
811	mRays[index].mMutations++;
812	mRays[index].mUnsuccessfulMutations++;
813
814	return true;
815	#endif
816	}
817
818
819
820
821	MutationBasedDistribution::MutationBasedDistribution(Preprocessor &preprocessor
822	) :
823	SamplingStrategy(preprocessor)
824	{
825	mType = MUTATION_BASED_DISTRIBUTION;
826	mBufferStart = 0;
827	mMaxRays = 500000;
828	mRays.reserve(mMaxRays);
829	mOriginMutationSize = 10.0f;
830	mLastIndex = 0;
831	// mOriginMutationSize = Magnitude(preprocessor.mViewCellsManager->
832	// GetViewSpaceBox().Diagonal())*1e-3;
833
834	}
835
836
837
838	}

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

Download in other formats: