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

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


using namespace std;

HaltonSequence SphericalSampleGenerator::sHalton;
HaltonSequence PoissonDiscSampleGenerator::sHalton;
HaltonSequence GaussianSampleGenerator::sHalton;


SampleGenerator::SampleGenerator(int numSamples, float radius):
mNumSamples(numSamples), mRadius(radius)
{}


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


void PoissonDiscSampleGenerator::Generate(float *samples) const
{
        // this is a hacky poisson sampling generator which does random dart-throwing
        // until it is not able to place any dart for a number of tries
        // in this case, the required min distance is reduced
        // the solution is a possion sampling with respect to the adjusted min distance
        // better solutions have been proposed, i.e., using hierarchical sampling

        const float maxTries = 1000;
        const float f_reduction = 0.9f;

        static HaltonSequence halton;
        float r[2];

        // generates poisson distribution on disc
        float minDist = 2.0f / sqrt((float)mNumSamples);

        //cout << "minDist before= " << minDist << endl;
        Sample2 *s = (Sample2 *)samples;

        for (int i = 0; i < mNumSamples; ++ i)
        {
                int tries = 0, totalTries = 0;

                // repeat until valid sample was found
                while (1)
                {
                        ++ tries;
                        ++ totalTries;

                        halton.GetNext(2, 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 > 1)
                                continue;

                        bool sampleValid = true;

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

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

                        if (tries > maxTries)
                        {
                                minDist *= f_reduction;
                                tries = 0;
                        }
                }
                
                //cout << "sample: " << samples[i].x << " " << i << " " << samples[i].y << " r: " << sqrt(samples[i].x * samples[i].x  + samples[i].y * samples[i].y) << " tries: " << totalTries << endl;
        }

        //cout << "minDist after= " << minDist << endl;
}



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


void GaussianSampleGenerator::Generate(float *samples) const
{
        float r[2];

        const float sigma = mRadius;

        Sample2 *s = (Sample2 *)samples;

        const float p0 = 1.0f / (sigma * sqrt(2.0f * M_PI));
        const float p1 = 1.0f / (sigma * sqrt(2.0f * M_PI));

        for (int i = 0; i < mNumSamples; ++ i)
        {
                sHalton.GetNext(2, r);

                float exp_x = -(r[0] * r[0]) / (2.0f * sigma * sigma);
                float exp_y = -(r[1] * r[1]) / (2.0f * sigma * sigma);

                s[i].x = p0 * pow(M_E, exp_x);
                s[i].y = p1 * pow(M_E, exp_y);
        }

        //cout << "minDist after= " << minDist << endl;
}


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


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

        // idea: use poisson distribution to generate samples
        PoissonDiscSampleGenerator poisson(mNumSamples, 1.0f);
        Sample2 *pSamples = new Sample2[mNumSamples];

        poisson.Generate((float *)pSamples);

        for (int i = 0; i < mNumSamples; ++ i)
        {
                //sHalton.GetNext(2, r);
                r[0] = pSamples[i].x; r[1] = pSamples[i].y;

                // create stratified samples over sphere
                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);
        }

        delete [] pSamples;
}