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

Revision 1867, 18.7 KB checked in by bittner, 18 years ago (diff)

merge, global lines, rss sampling updates

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