Ignore:
Timestamp:
01/18/07 10:25:29 (17 years ago)
Author:
bittner
Message:

mutation updates

File:
1 edited

Legend:

Unmodified
Added
Removed

  • GTP/trunk/Lib/Vis/Preprocessing/src/SamplingStrategy.cpp

    r1986 r1989  
    66#include "AxisAlignedBox3.h" 
    77#include "RssTree.h" 
    8 #include "Vector2.h" 
    9 #include "RndGauss.h" 
     8#include "Mutation.h" 
    109 
    1110namespace GtpVisibilityPreprocessor { 
    12  
    13 #define MUTATION_USE_CDF 0 
    14  
    1511 
    1612//HaltonSequence SamplingStrategy::sHalton; 
    … …  
    4844        // tmp changed matt. Q: should one rejected sample  
    4945        // terminate the whole method? 
    50         if (0) 
    51         { 
    52                 for (; i < number; i++) { 
    53                         if (!GenerateSample(ray)) 
    54                                 return i; 
    55                         rays.push_back(ray); 
    56                 } 
    57         } 
    58         else 
    59         { 
    60                 for (; i < number; i++)  
    61                 { 
    62                         int j = 0; 
    63                         bool sampleGenerated = false; 
    64  
    65                         for (j = 0; !sampleGenerated && (j < maxTries); ++ j) 
    66                         { 
    67                                 sampleGenerated = GenerateSample(ray); 
    68  
    69                                 if (sampleGenerated) 
    70                                 {                
    71                                         ++ samples; 
    72                                         rays.push_back(ray); 
    73                                 } 
    74                         } 
    75                 } 
    76         } 
    77          
     46        for (; i < number; i++)  
     47          { 
     48                int j = 0; 
     49                bool sampleGenerated = false; 
     50                 
     51                for (j = 0; !sampleGenerated && (j < maxTries); ++ j) 
     52                  { 
     53                        sampleGenerated = GenerateSample(ray); 
     54                         
     55                        if (sampleGenerated) 
     56                          {              
     57                                ++ samples; 
     58                                rays.push_back(ray); 
     59                          } 
     60                  } 
     61          } 
     62 
     63 
    7864        return samples; 
    7965} 
    … …  
    686672} 
    687673 
    688 void 
    689 MutationBasedDistribution::Update(VssRayContainer &vssRays) 
    690 { 
    691   //  for (int i=0; i < mRays.size(); i++) 
    692   //    cout<<mRays[i].mSamples<<" "; 
    693   //  cout<<endl; 
    694   cerr<<"Muattion update..."<<endl; 
    695   cerr<<"rays = "<<mRays.size()<<endl; 
    696   if (mRays.size()) { 
    697         cerr<<"Oversampling factors = "<< 
    698           GetEntry(0).mSamples<<" "<< 
    699           GetEntry(1).mSamples<<" "<< 
    700           GetEntry(2).mSamples<<" "<< 
    701           GetEntry(3).mSamples<<" "<< 
    702           GetEntry(4).mSamples<<" "<< 
    703           GetEntry(5).mSamples<<" ... "<< 
    704           GetEntry(mRays.size()-6).mSamples<<" "<< 
    705           GetEntry(mRays.size()-5).mSamples<<" "<< 
    706           GetEntry(mRays.size()-4).mSamples<<" "<< 
    707           GetEntry(mRays.size()-3).mSamples<<" "<< 
    708           GetEntry(mRays.size()-2).mSamples<<" "<< 
    709           GetEntry(mRays.size()-1).mSamples<<endl; 
    710   } 
    711   int contributingRays = 0; 
    712   for (int i=0; i < vssRays.size(); i++) { 
    713         if (vssRays[i]->mPvsContribution) { 
    714           contributingRays++; 
    715           if (mRays.size() < mMaxRays) { 
    716                 VssRay *newRay = new VssRay(*vssRays[i]); 
    717                 // add this ray 
    718                 newRay->Ref(); 
    719                 mRays.push_back(RayEntry(newRay)); 
    720           } else { 
    721                 // unref the old ray 
    722                 *mRays[mBufferStart].mRay = *vssRays[i]; 
    723                 mRays[mBufferStart].mSamples = 0; 
    724                 //              mRays[mBufferStart] = RayEntry(newRay); 
    725                 mBufferStart++; 
    726                 if (mBufferStart >= mMaxRays) 
    727                   mBufferStart = 0; 
    728           } 
    729         } 
    730   } 
    731    
    732   float pContributingRays = contributingRays/(float)vssRays.size(); 
    733   float importance = 1.0f; //sqr(1.0f/(pContributingRays + 1e-5)); 
    734   // set this values for last contributingRays 
    735   int index = mBufferStart - 1; 
    736    
    737   for (int i=0; i < contributingRays; i++, index--) { 
    738         if (index < 0) 
    739           index = mRays.size()-1; 
    740         mRays[index].mImportance = importance; 
    741   } 
    742  
    743 #if MUTATION_USE_CDF 
    744   // compute cdf 
    745   mRays[0].mCdf = mRays[0].mImportance/(mRays[0].mSamples+1); 
    746   for (int i=1; i < mRays.size(); i++)  
    747         mRays[i].mCdf = mRays[i-1].mCdf + mRays[i].mImportance/(mRays[i].mSamples+1); 
    748    
    749   float scale = 1.0f/mRays[i-1].mCdf; 
    750   for (i=0; i < mRays.size(); i++) { 
    751         mRays[i].mCdf *= scale; 
    752   } 
    753 #endif 
    754    
    755   cout<<"Importance = "<< 
    756         GetEntry(0).mImportance<<" "<< 
    757         GetEntry(mRays.size()-1).mImportance<<endl; 
    758    
    759   cerr<<"Mutation update done."<<endl; 
    760 } 
    761  
    762  
    763 Vector3 
    764 MutationBasedDistribution::ComputeOriginMutation(const VssRay &ray, 
    765                                                                                                  const Vector3 &U, 
    766                                                                                                  const Vector3 &V, 
    767                                                                                                  const Vector2 vr2, 
    768                                                                                                  const float radius 
    769                                                                                                  ) 
    770 { 
    771 #if 0 
    772   Vector3 v; 
    773   if (d.DrivingAxis() == 0) 
    774         v = Vector3(0, r[0]-0.5f, r[1]-0.5f); 
    775   else 
    776         if (d.DrivingAxis() == 1) 
    777           v = Vector3(r[0]-0.5f, 0, r[1]-0.5f); 
    778         else 
    779           v = Vector3(r[0]-0.5f, r[1]-0.5f, 0); 
    780   return v*(2*radius); 
    781 #endif 
    782 #if 0 
    783   return (U*(r[0] - 0.5f) + V*(r[1] - 0.5f))*(2*radius); 
    784 #endif 
    785  
    786  
    787   // Output random variable 
    788   Vector2 gaussvec2;  
    789    
    790   // Here we apply transform to gaussian, so 2D bivariate 
    791   // normal distribution 
    792   //  float sigma = ComputeSigmaFromRadius(radius); 
    793   float sigma = radius; 
    794   GaussianOn2D(vr2, 
    795                            sigma, // input 
    796                            gaussvec2);  // output 
    797  
    798          
    799     // Here we tranform the point correctly to 3D space using base 
    800     // vectors of the 3D space defined by the direction 
    801     Vector3 shift = gaussvec2.xx * U + gaussvec2.yy * V; 
    802          
    803         //      cout<<shift<<endl; 
    804         return shift; 
    805 } 
    806  
    807 Vector3 
    808 MutationBasedDistribution::ComputeTerminationMutation(const VssRay &ray, 
    809                                                                                                           const Vector3 &U, 
    810                                                                                                           const Vector3 &V, 
    811                                                                                                           const Vector2 vr2, 
    812                                                                                                           const float radius 
    813                                                                                                           ) 
    814 { 
    815 #if 0 
    816   Vector3 v; 
    817   // mutate the termination 
    818   if (d.DrivingAxis() == 0) 
    819         v = Vector3(0, r[2]-0.5f, r[3]-0.5f); 
    820   else 
    821         if (d.DrivingAxis() == 1) 
    822           v = Vector3(r[2]-0.5f, 0, r[3]-0.5f); 
    823         else 
    824           v = Vector3(r[2]-0.5f, r[3]-0.5f, 0); 
    825    
    826   //   Vector3 nv; 
    827    
    828   //   if (Magnitude(v) > Limits::Small) 
    829   //    nv = Normalize(v); 
    830   //   else 
    831   //    nv = v; 
    832    
    833   //  v = nv*size + v*size; 
    834  
    835   return v*(4.0f*radius); 
    836 #endif 
    837 #if 0 
    838   return (U*(vr2.xx - 0.5f) + V*(vr2.yy - 0.5f))*(4.0f*radius); 
    839 #endif 
    840   Vector2 gaussvec2;  
    841 #if 1 
    842   float sigma = radius; 
    843   GaussianOn2D(vr2, 
    844                            sigma, // input 
    845                            gaussvec2);  // output 
    846   Vector3 shift = gaussvec2.xx * U + gaussvec2.yy * V; 
    847   //    cout<<shift<<endl; 
    848   return shift; 
    849 #endif 
    850 #if 0 
    851   // Here we estimate standard deviation (sigma) from radius 
    852   float sigma = 1.1f*ComputeSigmaFromRadius(radius); 
    853   Vector3 vr3(vr2.xx, vr2.yy, RandomValue(0,1)); 
    854   PolarGaussianOnDisk(vr3, 
    855                                           sigma, 
    856                                           radius, // input 
    857                                           gaussvec2); // output 
    858    
    859   // Here we tranform the point correctly to 3D space using base 
    860   // vectors of the 3D space defined by the direction 
    861   Vector3 shift = gaussvec2.xx * U + gaussvec2.yy * V; 
    862    
    863   //    cout<<shift<<endl; 
    864   return shift; 
    865 #endif 
    866 } 
    867  
    868  
    869 bool 
    870 MutationBasedDistribution::GenerateSample(SimpleRay &sray) 
    871 { 
    872   float rr[5]; 
    873  
    874   if (mRays.size() == 0) { 
    875         // use direction based distribution 
    876         Vector3 origin, direction; 
    877         static HaltonSequence halton; 
    878          
    879         halton.GetNext(5, rr); 
    880         mPreprocessor.mViewCellsManager->GetViewPoint(origin, 
    881                                                                                                   Vector3(rr[0], rr[1], rr[2])); 
    882          
    883  
    884         direction = UniformRandomVector(rr[3], rr[4]); 
    885          
    886         const float pdf = 1.0f; 
    887         sray = SimpleRay(origin, direction, MUTATION_BASED_DISTRIBUTION, pdf); 
    888  
    889         return true; 
    890   }  
    891  
    892   int index; 
    893  
    894 #if !MUTATION_USE_CDF 
    895   // get tail of the buffer 
    896   index = (mLastIndex+1)%mRays.size(); 
    897   if (mRays[index].GetSamplingFactor() > 
    898           mRays[mLastIndex].GetSamplingFactor()) { 
    899         // search back for index where this is valid 
    900         index = (mLastIndex - 1 + mRays.size())%mRays.size(); 
    901         for (int i=0; i < mRays.size(); i++) { 
    902            
    903           //      if (mRays[index].mSamples > mRays[mLastIndex].mSamples) 
    904           //            break; 
    905           if (mRays[index].GetSamplingFactor() > 
    906                   mRays[mLastIndex].GetSamplingFactor() ) 
    907                 break; 
    908           index = (index - 1 + mRays.size())%mRays.size(); 
    909         } 
    910         // go one step back 
    911         index = (index+1)%mRays.size(); 
    912   } 
    913 #else 
    914   static HaltonSequence iHalton; 
    915   iHalton.GetNext(1, rr); 
    916   //rr[0] = RandomValue(0,1); 
    917   // use binary search to find index with this cdf 
    918   int l=0, r=mRays.size()-1; 
    919   while(l<r) { 
    920         int i = (l+r)/2; 
    921         if (rr[0] < mRays[i].mCdf ) 
    922           r = i; 
    923         else 
    924           l = i+1; 
    925   } 
    926   index = l; 
    927   //  if (rr[0] >= mRays[r].mCdf) 
    928   //    index = r; 
    929   //  else 
    930   //    index = l; 
    931          
    932          
    933 #endif 
    934   //  cout<<index<<" "<<rr[0]<<" "<<mRays[index].mCdf<<" "<<mRays[(index+1)%mRays.size()].mCdf<<endl; 
    935    
    936   VssRay *ray = mRays[index].mRay; 
    937   mRays[index].mSamples++; 
    938   mLastIndex = index; 
    939  
    940   mRays[index].mHalton.GetNext(4, rr); 
    941  
    942   // mutate the origin 
    943   Vector3 d = ray->GetDir(); 
    944  
    945  
    946   AxisAlignedBox3 box = ray->mTerminationObject->GetBox(); 
    947   float objectRadius = 0.5f*Magnitude(box.Diagonal()); 
    948   //  cout<<objectRadius<<endl; 
    949   if (objectRadius < Limits::Small) 
    950         return false; 
    951    
    952   // Compute right handed coordinate system from direction 
    953   Vector3 U, V; 
    954   Vector3 nd = Normalize(d); 
    955   nd.RightHandedBase(U, V); 
    956  
    957   Vector3 origin = ray->mOrigin; 
    958   Vector3 termination = ray->mTermination; //box.Center(); //ray->mTermination; //box.Center(); 
    959  
    960   float radiusExtension = 1.0f; 
    961   //  + mRays[index].mSamples/50.0f; 
    962    
    963   //  origin += ComputeOriginMutation(*ray, U, V, Vector2(r[0], r[1]), 0.5f*mOriginMutationSize*radiusExtension); 
    964   origin += ComputeOriginMutation(*ray, U, V, 
    965                                                                   Vector2(rr[0], rr[1]), 
    966                                                                   objectRadius*radiusExtension); 
    967   termination += ComputeTerminationMutation(*ray, U, V, 
    968                                                                                         Vector2(rr[2], rr[3]), 
    969                                                                                         objectRadius*radiusExtension); 
    970    
    971   Vector3 direction = termination - origin; 
    972    
    973   if (Magnitude(direction) < Limits::Small) 
    974         return false; 
    975    
    976   // shift the origin a little bit 
    977   origin += direction*0.5f; 
    978    
    979   direction.Normalize(); 
    980    
    981   // $$ jb the pdf is yet not correct for all sampling methods! 
    982   const float pdf = 1.0f; 
    983    
    984   sray = SimpleRay(origin, direction, MUTATION_BASED_DISTRIBUTION, pdf); 
    985   return true; 
    986 } 
    987  
    988 MutationBasedDistribution::MutationBasedDistribution(Preprocessor &preprocessor 
    989                                                                                                          ) : 
    990   SamplingStrategy(preprocessor) 
    991 { 
    992   mType = MUTATION_BASED_DISTRIBUTION; 
    993   mBufferStart = 0; 
    994   mMaxRays = 500000; 
    995   mRays.reserve(mMaxRays); 
    996   mOriginMutationSize = 10.0f; 
    997   mLastIndex = 0; 
    998   //  mOriginMutationSize = Magnitude(preprocessor.mViewCellsManager-> 
    999   //                                                              GetViewSpaceBox().Diagonal())*1e-3; 
    1000    
    1001 } 
    1002674 
    1003675 
Note: See TracChangeset for help on using the changeset viewer.