Context Navigation

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

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

Rev	Line
[1991]	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
[1997]	20
	21	#define USE_SIL_TERMINATION_MUTATION 1
	22
[1991]	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
[1997]	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;
[1991]	273
[1997]	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	}
[1991]	312
[1997]	313	// cout<<i<<endl;
	314	return dirs[i];
	315	}
	316
	317
[1991]	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();
[1995]	428
	429	// optimal for Pompeii 0.1f
	430	// optimal for Vienna 0.5f
[1991]	431
	432	float radiusExtension = 0.05f;
	433	// + mRays[index].mMutations/50.0f;
[1995]	434
	435	float mutationRadius = objectRadius*radiusExtension;
	436
	437	// tmp for pompeii
	438	mutationRadius = 0.22f;
[1991]	439
	440	origin += ComputeOriginMutation(*ray, U, V,
	441	Vector2(rr[0], rr[1]),
[1995]	442	mutationRadius);
[1991]	443
[1997]	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
[1991]	452	termination += ComputeTerminationMutation(*ray, U, V,
	453	Vector2(rr[2], rr[3]),
[1995]	454	mutationRadius);
[1997]	455	#endif
[1991]	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: