source: GTP/trunk/App/Demos/Vis/FriendlyCulling/src/SampleGenerator.cpp @ 3355

#include "SampleGenerator.h"
#include "common.h"


using namespace std;
using namespace CHCDemoEngine;


SampleGenerator::SampleGenerator(int numSamples, float radius):
mNumSamples(numSamples), mRadius(radius)
{
        mHalton = new HaltonSequence(2);
}


SampleGenerator::~SampleGenerator()
{
        DEL_PTR(mHalton);
}


PoissonDiscSampleGenerator2D::PoissonDiscSampleGenerator2D(int numSamples, float radius):
SampleGenerator(numSamples, radius)
{}


void PoissonDiscSampleGenerator2D::Generate(float *samples) const
{
        // Poisson disc sampling generator using relaxation 
        // dart-throwing as proposed by McCool et al.
        // the min distance requirement is relaxed if we are not 
        // able to place any dart for a number of tries
        // the solution is a possion sampling with respect 
        // to the adjusted min distance
        // this sampling scheme has the benefit that it is hierarchical 
        int maxTries = 1000;
        const float f_reduction = 0.95f;
        
        float r[2];

        // the maximal possible radius: our radius is a fraction of this radius
        // this is used as a measure of the quality of distribution of the point samples
        //const float rmax = 2.0f * mRadius * sqrt(1.0f / (2.0f * sqrt(3.0f) * mNumSamples));
        float rmax = 0.5f; //mRadius * sqrt(1.0f / (2.0f * sqrt(3.0f) * mNumSamples));

        // generates poisson distribution on disc
        // start with some thresholds: all samples lie on the circumference of circle
         
        //float minDist = 2.0f * rmax;
        float minDist = rmax;
        float sqrMinDist = minDist * minDist;

        int tries = 0;

        //cout << "minDist before= " << minDist / rmax << endl;
        cout << "minDist before= " << rmax << endl;

        Sample2 *s = (Sample2 *)samples;

        // check if on disc
        for (int i = 0; i < mNumSamples; ++ i)
        {
                // repeat until valid sample was found
                while (1)
                {
                        // q: should we use halton or does it conflict with the poisson disc properties?
                        //r[0] = RandomValue(0, 1); r[1] = RandomValue(0, 1);
                        mHalton->GetNext(r);

#if 0
                        // scale to -1 .. 1
                        const float rx = r[0] * 2.0f - 1.0f;
                        const float ry = r[1] * 2.0f - 1.0f;

                        // check if in disk, else exit early
                        const float distanceSquared = rx * rx + ry * ry;
                        if (rx * rx + ry * ry > mRadius * mRadius) 
                                // also avoid case that sample exactly in center
                                //|| (distanceSquared <= 1e-3f)
                        {
                                continue;
                        }
#else
                        const float rx = r[0];
                        const float ry = r[1];
#endif

                        bool sampleValid = true;

                        // check poisson property
                        for (int j = 0; ((j < i) && sampleValid); ++ j)
                        {
                                const float dist = 
                                        (s[j].x - rx) * (s[j].x - rx) +
                                        (s[j].y - ry) * (s[j].y - ry);
                        
                                if (dist < sqrMinDist) sampleValid = false;
                        }

                        if (sampleValid)
                        {
                                s[i].x = rx;
                                s[i].y = ry;
                                break;
                        }

                        ++ tries;

                        if (tries > maxTries)
                        {
                                minDist *= f_reduction;
                                sqrMinDist = minDist * minDist;

                                maxTries += 1000;
                        }
                }
        }

        for (int i = 0; i < mNumSamples; ++ i)
        {
                const float a = 2.0f * M_PI * s[i].x;
                const float r = sqrt(s[i].y);
                
                const float rad = mRadius * r;

                s[i].x = rad * cos(a);
                s[i].y = rad * sin(a);
        }

        rmax = mRadius * sqrt(1.0f / (2.0f * sqrt(3.0f) * mNumSamples));
        cout << "minDist after= " << (float)minDist * mNumSamples << " #tries: " << tries << " samples: " << mNumSamples << endl;
        cout << "minDist after= " << (float)minDist / rmax << " #tries: " << tries << endl;
}


RandomSampleGenerator2D::RandomSampleGenerator2D(int numSamples, float radius):
SampleGenerator(numSamples, radius)
{}


#if 0
void RandomSampleGenerator2D::Generate(float *samples) const
{
        Sample2 *s = (Sample2 *)samples;

        int numSamples = 0;

        float r[2];

        while (numSamples < mNumSamples)
        {
                mHalton->GetNext(r);
                
                const float rx = r[0] * 2.0f - 1.0f;
                const float ry = r[1] * 2.0f - 1.0f;

                // check if in disk, else exit early
                if (rx * rx + ry * ry > mRadius * mRadius)
                        continue;

                s[numSamples].x = rx;
                s[numSamples].y = ry;

                ++ numSamples;
        }
}

#else

void RandomSampleGenerator2D::Generate(float *samples) const
{
        Sample2 *s = (Sample2 *)samples;

        int numSamples = 0;
        float x[2];

        static float total1 = 0;
        static float total2 = 0;
        static int totalSamples = 0;

        for (int i = 0; i < mNumSamples; ++ i)
        {
                //x[0] = RandomValue(0, 1); x[1] = RandomValue(0, 1);
                mHalton->GetNext(x);
                
                const float a = 2.0f * M_PI * x[0];
                const float r = mRadius * sqrt(x[1]);

                s[i].x = r * cos(a);
                s[i].y = r * sin(a);

                /*total1 += x[0];
                total2 += x[1];
                totalSamples ++;

                if (totalSamples % 1000 == 1)
                {
                        float n1 = (float)total1 / totalSamples;
                        float n2 = (float)total2 / totalSamples;

                        cout << "here3 " << n1 << " " << n2 << endl;
                }*/
        }
}

#endif

SphericalSampleGenerator3D::SphericalSampleGenerator3D(int numSamples, float radius):
SampleGenerator(numSamples, radius)
{}


void SphericalSampleGenerator3D::Generate(float *samples) const
{
        float r[2];
        Sample3 *s = (Sample3 *)samples;

        for (int i = 0; i < mNumSamples; ++ i)
        {
                r[0] = RandomValue(0, 1);
                r[1] = RandomValue(0, 1);

                const float theta = 2.0f * acos(sqrt(1.0f - r[0]));
                const float phi = 2.0f * M_PI * r[1];
 
                s[i].x = mRadius * sin(theta) * cos(phi);
                s[i].y = mRadius * sin(theta) * sin(phi);
                s[i].z = mRadius * cos(theta);
        }
}


QuadraticDiscSampleGenerator2D::QuadraticDiscSampleGenerator2D(int numSamples, 
                                                                                                                           float radius):
//PoissonDiscSampleGenerator2D(numSamples, radius)
RandomSampleGenerator2D(numSamples, radius)
{
}


void QuadraticDiscSampleGenerator2D::Generate(float *samples) const
{
        Sample2 *s = (Sample2 *)samples;

        int numSamples = 0;
        float x[2];

        for (int i = 0; i < mNumSamples; ++ i)
        {
                x[0] = RandomValue(0, 1); x[1] = RandomValue(0, 1);
                //mHalton->GetNext(x);
                
                const float a = 2.0f * M_PI * x[0];
                const float r = sqrt(x[1]);
                //const float rad = mRadius * r * r * r * r;
                const float rad = mRadius * r;// * r;

                s[i].x = rad * cos(a);
                s[i].y = rad * sin(a);
        }
}