source: GTP/trunk/Lib/Vis/Preprocessing/src/MutualVisibility.cpp @ 2601

Revision 2601, 18.6 KB checked in by mattausch, 17 years ago (diff)
Line 
1#include <assert.h>
2#include <stack>
3#include "KdTree.h"
4#include "AxisAlignedBox3.h"
5#include "Ray.h"
6#include "MutualVisibility.h"
7#include "Exporter.h"
8#include "Mesh.h"
9#include "Triangle3.h"
10#include "SceneGraph.h"
11
12using namespace std;
13
14namespace GtpVisibilityPreprocessor {
15
16void
17RayShaft::Init(
18               const Rectangle3 &source,
19               const Rectangle3 &target)
20{
21  mDepth = 0;
22  mSource = source;
23  mTarget = target;
24}
25
26
27Vector3
28RayShaft::GetIntersectionPoint(const int rayIndex,
29                               const int depth) const
30{
31  Vector3 origin, direction;
32  Ray ray;
33  GetRaySegment(rayIndex, ray);
34  if (depth >= mSamples[rayIndex].mIntersections.size()) {
35    cerr<<"depth of sample out of limits"<<endl;
36    exit(1);
37  }
38  return ray.Extrap(mSamples[rayIndex].mIntersections[depth].mT);
39}
40
41void
42RayShaft::GetRaySegment(const int i, Ray &ray) const
43{
44  Vector3 origin, direction;
45  GetRay(i, origin, direction);
46  ray.Init(origin, direction, Ray::LINE_SEGMENT);
47  if ( mSamples[i].IsValid() ) {
48    origin = ray.Extrap(mSamples[i].mMinT);
49    direction = ray.Extrap(mSamples[i].mMaxT) - origin;
50    ray.Init(origin, direction, Ray::LINE_SEGMENT);
51  }
52}
53
54void
55RayShaft::GetRay(const int rayIndex,
56                 Vector3 &origin,
57                 Vector3 &direction) const
58{
59 
60  assert(rayIndex < 4);
61 
62  origin = mSource.mVertices[rayIndex];
63  direction = mTarget.mVertices[rayIndex] - origin;
64}
65
66void
67MutualVisibilitySampler::PerformSplit(
68                                      const RayShaft &sample,
69                                      const bool splitSource,
70                                      const int axis,
71                                      RayShaft &sample1,
72                                      RayShaft &sample2
73                                      )
74{
75 
76
77  // split the triangles
78
79  if (splitSource) {
80    sample1.mTarget = sample.mTarget;
81    sample2.mTarget = sample.mTarget;
82    sample.mSource.Split(
83                         axis,
84                         sample1.mSource,
85                         sample2.mSource);
86  } else {
87   
88    sample1.mSource = sample.mSource;
89    sample2.mSource = sample.mSource;
90   
91    sample.mTarget.Split(
92                         axis,
93                         sample1.mTarget,
94                         sample2.mTarget);
95  }
96 
97  // split the intersections
98  switch (axis) {
99  case 0:
100    sample1.mSamples[0].mIntersections = sample.mSamples[0].mIntersections;
101    sample1.mSamples[3].mIntersections = sample.mSamples[3].mIntersections;
102   
103    sample2.mSamples[1].mIntersections = sample.mSamples[1].mIntersections;
104    sample2.mSamples[2].mIntersections = sample.mSamples[2].mIntersections;
105    break;
106
107  case 1:
108    sample1.mSamples[0].mIntersections = sample.mSamples[0].mIntersections;
109    sample1.mSamples[1].mIntersections = sample.mSamples[1].mIntersections;
110
111    sample2.mSamples[2].mIntersections = sample.mSamples[2].mIntersections;
112    sample2.mSamples[3].mIntersections = sample.mSamples[3].mIntersections;
113    break;
114  }
115
116  // the intersections for the new shaft rays will be established
117  // later
118  sample1.mDepth = sample2.mDepth = sample.mDepth+1;
119}
120
121float
122MutualVisibilitySampler::GetSpatialAngle(const RayShaft &shaft,
123                                         const Vector3 &point
124                                         )
125{
126  const int sampleIndices[]={
127    0,1,2,
128    0,2,3
129  };
130 
131  float sum = 0.0f;
132  int i, j=0;
133
134  if (!shaft.IsValid())
135    return 2.0f*mSolidAngleThreshold;
136
137  for (i=0; i < 2; i++) {
138    Triangle3 triangle(shaft.GetIntersectionPoint(sampleIndices[j++], 0),
139                       shaft.GetIntersectionPoint(sampleIndices[j++], 0),
140                       shaft.GetIntersectionPoint(sampleIndices[j++], 0));
141    sum += triangle.GetSpatialAngle(point);
142  }
143
144  return sum;
145}
146
147
148void
149MutualVisibilitySampler::ComputeError(RayShaft &shaft)
150{
151  // evaluate minimal error which can be achieved by more precise evaluation
152  // if this is above the threshold do not proceed further
153
154  // check whether the samples pierce the same mesh
155  Intersectable::NewMail();
156  int i, j;
157  for (i=0; i < 4; i++) {
158    RaySample *sample = &shaft.mSamples[i];
159    for (j=0; j < sample->mIntersections.size(); j++)
160      sample->mIntersections[j].mObject->Mail();
161  }
162
163  //  int id = Mailable::GetThreadId();
164  for (i=0; i < 4; i++) {
165    RaySample *sample = &shaft.mSamples[i];
166    for (j=0; j < sample->mIntersections.size(); j++) {
167      if (sample->mIntersections[j].mObject->IncMail() == 4) {
168                cerr<<"T";
169                shaft.mError = 0.0f;
170                return;
171      }
172    }
173  }
174   
175 
176  float maxAngle = GetSpatialAngle(shaft,
177                                   shaft.mTarget.GetCenter());
178 
179  for (i=0; i < 3; i++) {
180    float angle = GetSpatialAngle(shaft,
181                                  shaft.mTarget.mVertices[i]);
182    if (angle > maxAngle)
183      maxAngle = angle;
184  }
185
186  maxAngle = MAX_FLOAT;
187  shaft.mError = maxAngle;
188}
189
190
191void
192MutualVisibilitySampler::ConstructInitialSamples2(
193                                                  const AxisAlignedBox3 &source,
194                                                  const AxisAlignedBox3 &target,
195                                                  vector<RayShaft *> &samples
196                                                  )
197{
198  // get all rectangles potentially visible from the source box
199  int i;
200  int sourceMask = 0;
201  for (i=0; i < 8; i++)
202    sourceMask |= source.GetFaceVisibilityMask(target.GetVertex(i));
203
204  // now for each visble source face find all visible target faces
205  for (i=0; i < 6; i++) {
206    Rectangle3 sourceFace = source.GetFace(i);
207    if ( sourceMask &(1<<i) ) {
208      int mask = target.GetFaceVisibilityMask(sourceFace);
209      // construct triangle samples for all visible rectangles
210      for (int j=0; j < 6; j++)
211        if (mask & (1<<j)) {
212          AddInitialSamples2(sourceFace, target.GetFace(j), samples);
213        }
214    }
215  }
216}
217
218void
219MutualVisibilitySampler::ConstructInitialSamples3(
220                                                  const AxisAlignedBox3 &source,
221                                                  const AxisAlignedBox3 &target,
222                                                  vector<RayShaft *> &samples
223                                                  )
224{
225  // get all rectangles potentially visible from the source box
226  int i;
227  int sourceMask = 0;
228  for (i=0; i < 8; i++)
229    sourceMask |= source.GetFaceVisibilityMask(target.GetVertex(i));
230
231  // now for each visble source face find all visible target faces
232  int face;
233  for (face=0; face < 6; face++) {
234    if ( sourceMask & (1<<face) ) {
235      Rectangle3 sourceRect = source.GetFace(face);
236
237      Vector3 targetCenter = target.Center();
238     
239      Vector3 normal = sourceRect.GetNormal();
240     
241      Plane3 sourcePlane(normal, sourceRect.GetVertex(0));
242     
243      Plane3 targetPlane(normal, target.GetVertex(0));
244
245      int i;
246      for (i=1; i < 8; i++) {
247        Vector3 v = target.GetVertex(i);
248        if (targetPlane.Distance(v) < 0)
249          targetPlane = Plane3(normal, v);
250      }
251 
252      Vector3 xBasis = Normalize(sourceRect.GetVertex(1) - sourceRect.GetVertex(0));
253      Vector3 yBasis = Normalize(sourceRect.GetVertex(3) - sourceRect.GetVertex(0));
254
255      // cast rays between the centers of the boxes
256      Vector3 targetRCenter = targetPlane.FindIntersection(sourceRect.GetCenter(),
257                                                           targetCenter);
258     
259      Rectangle3 targetRect;
260     
261      float targetDist[4];
262      targetDist[0] = targetDist[1] = targetDist[2] = targetDist[3] = 0.0f;
263     
264      // cast rays between corresponding vertices of the boxes
265      int j;
266      int intersections=0;
267      for (i=0; i < 4; i++)
268        for (j=0; j < 8; j++) {
269          Vector3 v;
270          Vector3 diff;
271          bool coplanar;
272          float dist;
273         
274          v = targetPlane.FindIntersection(sourceRect.GetVertex(i),
275                                           target.GetVertex(j),
276                                           NULL,
277                                           &coplanar);
278          if (!coplanar) {
279            // evaluate target and
280            diff = targetRCenter - v;
281            dist = DotProd(diff, xBasis);
282            if (dist < targetDist[0])
283              targetDist[0] = dist;
284            if (dist > targetDist[1])
285              targetDist[1] = dist;
286
287            dist = DotProd(diff, yBasis);
288
289            if (dist < targetDist[2])
290              targetDist[2] = dist;
291            if (dist > targetDist[3])
292              targetDist[3] = dist;
293            intersections++;
294          }
295        }
296     
297      if (intersections>=4) {
298        targetRect.mVertices[0] = targetRCenter + targetDist[0]*xBasis + targetDist[2]*yBasis;
299        targetRect.mVertices[1] = targetRCenter + targetDist[1]*xBasis + targetDist[2]*yBasis;
300        targetRect.mVertices[2] = targetRCenter + targetDist[1]*xBasis + targetDist[3]*yBasis;
301        targetRect.mVertices[3] = targetRCenter + targetDist[0]*xBasis + targetDist[3]*yBasis;
302       
303        //      cout<<sourceRect<<targetRect<<endl;
304        AddInitialSamples(sourceRect, targetRect, samples);
305      }
306    }
307  }
308}
309
310void
311MutualVisibilitySampler::ConstructInitialSamples(
312                                                 const AxisAlignedBox3 &source,
313                                                 const AxisAlignedBox3 &target,
314                                                 vector<RayShaft *> &samples
315                                                 )
316{
317  Vector3 sourceCenter = source.Center();
318  Vector3 targetCenter = target.Center();
319  Vector3 normal = Normalize(sourceCenter - targetCenter );
320 
321  Plane3 sourcePlane(normal, source.GetVertex(0));
322  int i;
323  for (i=1; i < 8; i++) {
324    Vector3 v = source.GetVertex(i);
325    if (sourcePlane.Distance(v) > 0)
326      sourcePlane = Plane3(normal, v);
327  }
328 
329  Plane3 targetPlane(normal, target.GetVertex(0));
330  for (i=1; i < 8; i++) {
331    Vector3 v = target.GetVertex(i);
332    if (targetPlane.Distance(v) < 0)
333      targetPlane = Plane3(normal, v);
334  }
335
336 
337  Vector3 xBasis = CrossProd(Vector3(0,1,0), normal);
338
339  if (Magnitude(xBasis) > 1e-6)
340    xBasis.Normalize();
341  else {
342    xBasis = Normalize(CrossProd(Vector3(0,0,1), normal));
343  }
344
345  Vector3 yBasis = Normalize( CrossProd(normal, xBasis) );
346  // cast rays between the centers of the boxes
347  Vector3 targetRCenter = targetPlane.FindIntersection(sourceCenter,
348                                                       targetCenter);
349 
350  Vector3 sourceRCenter = sourcePlane.FindIntersection(sourceCenter,
351                                                       targetCenter);
352 
353
354  Rectangle3 sourceRect;
355  Rectangle3 targetRect;
356
357
358  if (0) {
359  float scale = Magnitude(source.Size())*0.7f;
360  sourceRect.mVertices[0] = sourceRCenter - scale*xBasis - scale*yBasis;
361  sourceRect.mVertices[1] = sourceRCenter + scale*xBasis - scale*yBasis;
362  sourceRect.mVertices[2] = sourceRCenter + scale*xBasis + scale*yBasis;
363  sourceRect.mVertices[3] = sourceRCenter - scale*xBasis + scale*yBasis;
364
365  scale = Magnitude(target.Size())*0.7f;
366  targetRect.mVertices[0] = targetRCenter - scale*xBasis - scale*yBasis;
367  targetRect.mVertices[1] = targetRCenter + scale*xBasis - scale*yBasis;
368  targetRect.mVertices[2] = targetRCenter + scale*xBasis + scale*yBasis;
369  targetRect.mVertices[3] = targetRCenter - scale*xBasis + scale*yBasis;
370  }
371
372
373  float sourceDist[4];
374  float targetDist[4];
375  sourceDist[0] = sourceDist[1] = sourceDist[2] = sourceDist[3] = 0.0f;
376  targetDist[0] = targetDist[1] = targetDist[2] = targetDist[3] = 0.0f;
377
378  // cast rays between corresponding vertices of the boxes
379  int j;
380  for (i=0; i < 8; i++)
381    for (j=0; j < 8; j++) {
382      Vector3 v;
383      Vector3 diff;
384      bool coplanar;
385      float dist;
386      v = sourcePlane.FindIntersection(source.GetVertex(i),
387                                       target.GetVertex(j),
388                                       NULL,
389                                       &coplanar);
390      if (!coplanar) {
391        // evaluate source and
392        diff = sourceRCenter - v;
393        dist = DotProd(diff, xBasis);
394        if (dist < sourceDist[0])
395          sourceDist[0] = dist;
396        if (dist > sourceDist[1])
397          sourceDist[1] = dist;
398
399        dist = DotProd(diff, yBasis);
400        if (dist < sourceDist[2])
401          sourceDist[2] = dist;
402        if (dist > sourceDist[3])
403          sourceDist[3] = dist;
404      }
405     
406      v = targetPlane.FindIntersection(source.GetVertex(i),
407                                       target.GetVertex(j),
408                                       NULL,
409                                       &coplanar);
410      if (!coplanar) {
411        // evaluate target and
412        diff = targetRCenter - v;
413        dist = DotProd(diff, xBasis);
414        if (dist < targetDist[0])
415          targetDist[0] = dist;
416        if (dist > targetDist[1])
417          targetDist[1] = dist;
418        dist = DotProd(diff, yBasis);
419        if (dist < targetDist[2])
420          targetDist[2] = dist;
421        if (dist > targetDist[3])
422          targetDist[3] = dist;
423      }
424    }
425
426  sourceRect.mVertices[0] = sourceRCenter + sourceDist[0]*xBasis + sourceDist[2]*yBasis;
427  sourceRect.mVertices[1] = sourceRCenter + sourceDist[1]*xBasis + sourceDist[2]*yBasis;
428  sourceRect.mVertices[2] = sourceRCenter + sourceDist[1]*xBasis + sourceDist[3]*yBasis;
429  sourceRect.mVertices[3] = sourceRCenter + sourceDist[0]*xBasis + sourceDist[3]*yBasis;
430
431  targetRect.mVertices[0] = targetRCenter + targetDist[0]*xBasis + targetDist[2]*yBasis;
432  targetRect.mVertices[1] = targetRCenter + targetDist[1]*xBasis + targetDist[2]*yBasis;
433  targetRect.mVertices[2] = targetRCenter + targetDist[1]*xBasis + targetDist[3]*yBasis;
434  targetRect.mVertices[3] = targetRCenter + targetDist[0]*xBasis + targetDist[3]*yBasis;
435
436
437  cout<<sourceRect<<targetRect<<endl;
438  AddInitialSamples(sourceRect, targetRect, samples);
439}
440
441void
442MutualVisibilitySampler::AddInitialSamples(
443                                           const Rectangle3 &sourceRect,
444                                           const Rectangle3 &targetRect,
445                                           vector<RayShaft *> &samples
446                                           )
447{
448
449  RayShaft *sample = new RayShaft(sourceRect,
450                                  targetRect);
451  samples.push_back(sample);
452}
453
454void
455MutualVisibilitySampler::AddInitialSamples2(
456                                            const Rectangle3 &sourceRect,
457                                            const Rectangle3 &targetRect,
458                                            vector<RayShaft *> &samples
459                                            )
460{
461  // align the rectangles properly
462  Rectangle3 sRect, tRect;
463  if (DotProd(sourceRect.GetNormal(), targetRect.GetNormal()) < -0.99f) {
464    int i;
465    for (i=0; i < 4; i++) {
466      tRect.mVertices[i] = targetRect.mVertices[( 7 - i )%4];
467    }
468   
469    RayShaft *sample = new RayShaft(sourceRect,
470                                    tRect);
471    samples.push_back(sample);
472  }
473 
474  return;
475 
476  float minDist = MAX_FLOAT;
477  int startI;
478  int startJ;
479  int i, j;
480
481 
482  for (i=0; i < 4; i++) {
483    for (j=0; j < 4; j++) {
484      float dist = Distance(sourceRect.mVertices[i], targetRect.mVertices[j]);
485      if (dist < minDist) {
486        minDist = dist;
487        startI = i;
488        startJ = j;
489      }
490    }
491  }
492  for (i=0; i < 4; i++) {
493    sRect.mVertices[i] = sourceRect.mVertices[(startI + i )%4];
494    tRect.mVertices[i] = targetRect.mVertices[(4 + startJ - i )%4];
495  }
496 
497  RayShaft *sample = new RayShaft(sRect,
498                                  tRect);
499  samples.push_back(sample);
500}
501
502
503void
504MutualVisibilitySampler::ExportShafts(vector<RayShaft *> &shafts,
505                                      const bool singleFile)
506{
507  static int id = 0;
508  char filename[64];
509  if (id > 20)
510    return;
511 
512  Exporter *exporter = NULL;
513  for (int i=0; i < shafts.size(); i++) {
514    if (!exporter) {
515      if (singleFile)
516        sprintf(filename, "shafts-single-%04d.x3d", id++);
517      else
518        sprintf(filename, "shafts%04d-%02d.x3d", id++, i);
519      exporter = Exporter::GetExporter(filename);
520    }
521   
522    exporter->SetWireframe();
523    // exporter->ExportScene(mSceneGraph->GetRoot());
524
525    exporter->ExportBox(mSource);
526    exporter->ExportBox(mTarget);
527
528    exporter->SetFilled();
529
530
531    RayShaft *shaft = shafts[i];
532    Mesh *mesh = new Mesh;
533    mesh->AddRectangle(shaft->mSource);
534    mesh->AddRectangle(shaft->mTarget);
535    vector<Ray *> rays;
536    for (int j=0; j < 4; j++) {
537      Ray *ray = new Ray;
538      shaft->GetRaySegment(j, *ray);
539      rays.push_back(ray);
540    }
541   
542    Material m = RandomMaterial();
543    exporter->SetForcedMaterial(m);
544    MeshInstance mi(mesh);
545    exporter->ExportIntersectable(&mi);
546    exporter->ExportRays(rays, -1.0f, m.mDiffuseColor);
547    if (!singleFile) {
548      delete exporter;
549      exporter = NULL;
550    }
551  }
552  if (exporter)
553    delete exporter;
554 
555}
556
557int
558MutualVisibilitySampler::CastRays(RayShaft &shaft)
559{
560  Ray ray;
561  int i;
562
563  for (i=0; i < 4; i++)
564    if (!shaft.mSamples[i].IsProcessed()) {
565      Vector3 origin, direction;
566      shaft.GetRay(i, origin, direction);
567      // determine intersections with the boxes
568      ray.Init(origin, direction, Ray::LINE_SEGMENT);
569      float stmin, stmax = 0.0f, ttmin=1.0f, ttmax;
570      bool valid = true;
571     
572      if (mUseBoxes) {
573        if (mSource.GetMinMaxT(ray, &stmin, &stmax) &&
574            mTarget.GetMinMaxT(ray, &ttmin, &ttmax)) {
575          shaft.mSamples[i].mMinT = stmax;
576          shaft.mSamples[i].mMaxT = ttmin;
577          origin = ray.Extrap(stmax);
578          direction = ray.Extrap(ttmin) - origin;
579          // reinit the ray
580          ray.Init(origin, direction, Ray::LINE_SEGMENT);
581        } else
582          valid = false;
583      } else {
584        shaft.mSamples[i].mMinT = 0.0f;
585        shaft.mSamples[i].mMaxT = 1.0f;
586      }
587      if (valid) {
588        if (!mKdTree->CastRay(ray)) {
589          cerr<<"V"<<endl;
590          return VISIBLE;
591        }
592        shaft.mSamples[i].mIntersections = ray.intersections;
593        cerr<<"I";
594      } else {
595        cerr<<"X";
596        shaft.mSamples[i].SetInvalid();
597      }
598    }
599  return INVISIBLE;
600}
601
602int
603MutualVisibilitySampler::ComputeVisibility()
604{
605  int result = INVISIBLE;
606
607  vector<RayShaft *> shafts;
608  ConstructInitialSamples3(mSource, mTarget, shafts);
609       
610  if (1)
611    ExportShafts(shafts, false);
612
613  stack<RayShaft *> shaftStack;
614  int i;
615  for (i=0; i < shafts.size(); i++)
616    shaftStack.push(shafts[i]);
617
618  shafts.clear();
619  Ray ray;
620 
621// now process the shafts as long as we have something to do
622  while (!shaftStack.empty()) {
623    RayShaft *shaft = shaftStack.top();
624    shaftStack.pop();
625   
626//      // cast a new ray
627//      int triangleSplitEdge = SetupExtremalRay(sample, source, ray);
628   
629    if (CastRays(*shaft) == VISIBLE) {
630      result = VISIBLE;
631      break;
632    }
633               
634    // compute error ....
635    ComputeError(*shaft);
636    cout<<shaft->mDepth<<"|";
637    if (shaft->IsValid())
638      shafts.push_back(shaft);
639
640    if (shaft->mDepth < 10 &&
641                                shaft->mError > mSolidAngleThreshold) {
642     
643      // generate 2 new samples
644      RayShaft *newSamples[2];
645      newSamples[0] = new RayShaft;
646      newSamples[1] = new RayShaft;
647     
648      // chose what to split
649      bool splitSource = shaft->mSource.GetArea() > shaft->mTarget.GetArea();
650      int axis;
651      if (splitSource) {
652                                axis = shaft->mSource.DominantAxis();
653      } else {
654                                axis = shaft->mTarget.DominantAxis();
655      }
656     
657      PerformSplit(*shaft, splitSource, axis, *newSamples[0], *newSamples[1]);
658      //      delete shaft;
659      shaftStack.push(newSamples[0]);
660      shaftStack.push(newSamples[1]);
661    } else {
662      // store a terminal shaft
663    }
664
665  }
666 
667  while (!shaftStack.empty()) {
668    RayShaft *shaft = shaftStack.top();
669    shaftStack.pop();
670    delete shaft;
671  }
672
673  if (0)
674    ExportShafts(shafts, true);
675
676  for (i=0; i < shafts.size(); i++) {
677    delete shafts[i];
678  }
679
680  return result;
681}
682
683MutualVisibilitySampler::MutualVisibilitySampler(SceneGraph *sceneGraph,
684                                                                                                 KdTree *kdTree,
685                                                                                                 const AxisAlignedBox3 &source,
686                                                                                                 const AxisAlignedBox3 &target,
687                                                                                                 const float solidAngleThreshold)
688{
689        mSceneGraph = sceneGraph;
690        mKdTree = kdTree;
691        mSource = source;
692        mTarget = target;
693        mUseBoxes = true;
694        mSolidAngleThreshold = solidAngleThreshold;
695}
696
697
698int
699ComputeBoxVisibility(SceneGraph *sceneGraph,
700                                         KdTree *kdTree,
701                                         const AxisAlignedBox3 &source,
702                                         const AxisAlignedBox3 &target,
703                                         const float solidAngleThreshold)
704{
705        MutualVisibilitySampler
706                sampler(sceneGraph, kdTree, source, target, solidAngleThreshold);
707
708        int visibility = sampler.ComputeVisibility();
709
710        return visibility;
711
712}
713
714}
Note: See TracBrowser for help on using the repository browser.