00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #include "metrosampler.h"
00027 #include "memory.h"
00028 #include "scene.h"
00029 #include "dynload.h"
00030
00031 using namespace lux;
00032
00033 #define SAMPLE_FLOATS 7
00034
00035
00036 static float mutate(const float x, const float randomValue)
00037 {
00038 static const float s1 = 1.f / 512.f, s2 = 1.f / 16.f;
00039 const float dx = s1 / (s1 / s2 + fabsf(2.f * randomValue - 1.f)) -
00040 s1 / (s1 / s2 + 1.f);
00041 if (randomValue < 0.5f) {
00042 float x1 = x + dx;
00043 return (x1 < 1.f) ? x1 : x1 - 1.f;
00044 } else {
00045 float x1 = x - dx;
00046 return (x1 < 0.f) ? x1 + 1.f : x1;
00047 }
00048 }
00049
00050
00051 static float mutateScaled(const float x, const float randomValue, const float mini, const float maxi, const float range)
00052 {
00053 static const float s1 = 32.f;
00054 const float dx = range / (1.f + s1 * fabsf(2.f * randomValue - 1.f)) -
00055 range / (1.f + s1);
00056 if (randomValue < 0.5f) {
00057 float x1 = x + dx;
00058 return (x1 < maxi) ? x1 : x1 - maxi + mini;
00059 } else {
00060 float x1 = x - dx;
00061 return (x1 < mini) ? x1 - mini + maxi : x1;
00062 }
00063 }
00064
00065 static const u_int rngN = 8191;
00066 static const u_int rngA = 884;
00067 static float rngDummy;
00068 #define rngGet(__pos) (modff(rngSamples[(rngBase + (__pos)) % rngN] + rngRotation[(__pos)], &rngDummy))
00069 #define rngGet2(__pos,__off) (modff(rngSamples[(rngBase + (__pos) + (__off)) % rngN] + rngRotation[(__pos)], &rngDummy))
00070
00071 MetropolisSampler::MetropolisSampler(int xStart, int xEnd, int yStart, int yEnd,
00072 int maxRej, float largeProb, float microProb, float rng,
00073 bool useV) :
00074 Sampler(xStart, xEnd, yStart, yEnd, 1), normalSamples(0), totalSamples(0),
00075 totalTimes(0), maxRejects(maxRej), consecRejects(0), pLarge(largeProb),
00076 pMicro(microProb), range(rng), useVariance(useV), sampleImage(NULL),
00077 timeImage(NULL), offset(NULL), rngRotation(NULL), rngBase(0),
00078 rngOffset(0), large(true), stamp(0), numMicro(-1), posMicro(-1), weight(0.f),
00079 LY(0.f), alpha(0.f), totalLY(0.f), sampleCount(0.f)
00080 {
00081
00082 rngSamples = AllocAligned<float>(rngN);
00083 rngSamples[0] = 0.f;
00084 rngSamples[1] = 1.f / rngN;
00085 u_int rngI = 1;
00086 for (u_int i = 2; i < rngN; ++i) {
00087 rngI = (rngI * rngA) % rngN;
00088 rngSamples[i] = static_cast<float>(rngI) / rngN;
00089 }
00090 }
00091
00092 MetropolisSampler::~MetropolisSampler() {
00093 FreeAligned(sampleImage);
00094 FreeAligned(timeImage);
00095 FreeAligned(rngSamples);
00096 FreeAligned(rngRotation);
00097 delete[] offset;
00098 }
00099
00100
00101 MetropolisSampler* MetropolisSampler::clone() const
00102 {
00103 MetropolisSampler *newSampler = new MetropolisSampler(*this);
00104 newSampler->totalSamples = 0;
00105 newSampler->sampleImage = NULL;
00106
00107 return newSampler;
00108 }
00109
00110 static void initMetropolis(MetropolisSampler *sampler, const Sample *sample)
00111 {
00112 u_int i;
00113
00114 sampler->normalSamples = SAMPLE_FLOATS;
00115 for (i = 0; i < sample->n1D.size(); ++i)
00116 sampler->normalSamples += sample->n1D[i];
00117 for (i = 0; i < sample->n2D.size(); ++i)
00118 sampler->normalSamples += 2 * sample->n2D[i];
00119
00120
00121 sampler->totalSamples = sampler->normalSamples;
00122 sampler->offset = new int[sample->nxD.size()];
00123 sampler->totalTimes = 0;
00124 for (i = 0; i < sample->nxD.size(); ++i) {
00125 sampler->offset[i] = sampler->totalSamples;
00126 sampler->totalTimes += sample->nxD[i];
00127 sampler->totalSamples += sample->dxD[i] * sample->nxD[i];
00128 }
00129
00130
00131 sampler->sampleImage = AllocAligned<float>(sampler->totalSamples);
00132 sampler->timeImage = AllocAligned<int>(sampler->totalTimes);
00133
00134
00135 sampler->contribBuffer = sampler->film->scene->contribPool->Next(NULL);
00136
00137
00138
00139
00140 sampler->rngOffset = sampler->totalSamples;
00141 if (sampler->rngOffset >= rngN)
00142 sampler->rngOffset = sampler->rngOffset % (rngN - 1) + 1;
00143
00144
00145 sampler->rngBase = rngN - sampler->rngOffset;
00146
00147 sampler->rngRotation = AllocAligned<float>(sampler->totalSamples);
00148 }
00149
00150
00151 bool MetropolisSampler::GetNextSample(Sample *sample, u_int *use_pos)
00152 {
00153 sample->sampler = this;
00154 const float mutationSelector = tspack->rng->floatValue();
00155 large = (mutationSelector < pLarge);
00156 if (sampleImage == NULL) {
00157 initMetropolis(this, sample);
00158 large = true;
00159 }
00160
00161
00162
00163 rngBase += rngOffset;
00164 if (rngBase > rngN - 1)
00165 rngBase -= rngN;
00166
00167
00168
00169 if (rngBase == 0) {
00170
00171
00172 if (film->enoughSamplePerPixel)
00173 return false;
00174 for (int i = 0; i < totalSamples; ++i)
00175 rngRotation[i] = tspack->rng->floatValue();
00176 }
00177 if (large) {
00178
00179
00180 sample->imageX = rngGet(0) * (xPixelEnd - xPixelStart) + xPixelStart;
00181 sample->imageY = rngGet(1) * (yPixelEnd - yPixelStart) + yPixelStart;
00182 sample->lensU = rngGet(2);
00183 sample->lensV = rngGet(3);
00184 sample->time = rngGet(4);
00185 sample->wavelengths = rngGet(5);
00186 sample->singleWavelength = rngGet(6);
00187 for (int i = SAMPLE_FLOATS; i < normalSamples; ++i)
00188 sample->oneD[0][i - SAMPLE_FLOATS] = rngGet(i);
00189
00190 for (int i = 0; i < totalTimes; ++i)
00191 sample->timexD[0][i] = -1;
00192
00193 sample->stamp = 0;
00194 numMicro = -1;
00195 } else {
00196
00197
00198 if (1.f - mutationSelector < pMicro)
00199 numMicro = min<int>(sample->nxD.size(), Float2Int((1.f - mutationSelector) / pMicro * (sample->nxD.size() + 1)));
00200 else
00201 numMicro = -1;
00202 if (numMicro > 0) {
00203
00204
00205 u_int maxPos = 0;
00206 for (; maxPos < sample->nxD[numMicro - 1]; ++maxPos)
00207 if (sample->timexD[numMicro - 1][maxPos] < sample->stamp)
00208 break;
00209 posMicro = min<int>(sample->nxD[numMicro - 1] - 1, Float2Int(tspack->rng->floatValue() * maxPos));
00210 } else {
00211
00212
00213 posMicro = -1;
00214 sample->imageX = mutateScaled(sampleImage[0],
00215 rngGet(0), xPixelStart, xPixelEnd, range);
00216 sample->imageY = mutateScaled(sampleImage[1],
00217 rngGet(1), yPixelStart, yPixelEnd, range);
00218 sample->lensU = mutate(sampleImage[2], rngGet(2));
00219 sample->lensV = mutate(sampleImage[3], rngGet(3));
00220 sample->time = mutate(sampleImage[4], rngGet(4));
00221 sample->wavelengths = mutate(sampleImage[5], rngGet(5));
00222 sample->singleWavelength = mutateScaled(sampleImage[6],
00223 rngGet(6), 0.f, 1.f, 1.f);
00224 for (int i = SAMPLE_FLOATS; i < normalSamples; ++i)
00225 sample->oneD[0][i - SAMPLE_FLOATS] =
00226 mutate(sampleImage[i], rngGet(i));
00227 }
00228
00229 ++(sample->stamp);
00230 }
00231 return true;
00232 }
00233
00234 float *MetropolisSampler::GetLazyValues(Sample *sample, u_int num, u_int pos)
00235 {
00236
00237 const u_int size = sample->dxD[num];
00238 float *data = sample->xD[num] + pos * size;
00239
00240 int stampLimit = sample->stamp;
00241
00242 if (numMicro >= 0 && num != u_int(numMicro - 1) && pos != u_int(posMicro))
00243 --stampLimit;
00244
00245 if (sample->timexD[num][pos] != stampLimit) {
00246
00247
00248 if (sample->timexD[num][pos] == -1) {
00249 for (u_int i = 0; i < size; ++i)
00250 data[i] = rngGet(normalSamples + pos * size + i);
00251 sample->timexD[num][pos] = 0;
00252 } else {
00253 float *image = sampleImage + offset[num] + pos * size;
00254 for (u_int i = 0; i < size; ++i)
00255 data[i] = image[i];
00256 }
00257
00258 for (; sample->timexD[num][pos] < stampLimit; ++(sample->timexD[num][pos])) {
00259 for (u_int i = 0; i < size; ++i)
00260 data[i] = mutate(data[i], rngGet2(normalSamples + pos * size + i, rngOffset * (stampLimit - sample->timexD[num][pos] + 1)));
00261 }
00262 }
00263 return data;
00264 }
00265
00266 void MetropolisSampler::AddSample(const Sample &sample)
00267 {
00268 vector<Contribution> &newContributions(sample.contributions);
00269 float newLY = 0.f;
00270 for(u_int i = 0; i < newContributions.size(); ++i) {
00271 const float ly = newContributions[i].color.Y();
00272 if (ly > 0.f) {
00273 if (isinf(ly))
00274 newContributions[i].color = XYZColor(0.f);
00275 else {
00276 if (useVariance && newContributions[i].variance > 0.f)
00277 newLY += ly * newContributions[i].variance;
00278 else
00279 newLY += ly;
00280 }
00281 }
00282 }
00283
00284 if (large) {
00285 totalLY += newLY;
00286 ++sampleCount;
00287 }
00288 const float meanIntensity = totalLY > 0.f ? totalLY / sampleCount : 1.f;
00289
00290 contribBuffer->AddSampleCount(1.f);
00291
00292
00293 float accProb;
00294 if (LY > 0.f && consecRejects < maxRejects)
00295 accProb = min(1.f, newLY / LY);
00296 else
00297 accProb = 1.f;
00298 const float newWeight = accProb + (large ? 1.f : 0.f);
00299 weight += 1.f - accProb;
00300
00301 if (accProb == 1.f || tspack->rng->floatValue() < accProb) {
00302
00303 const float norm = weight / (LY / meanIntensity + pLarge);
00304 if (norm > 0.f) {
00305 for(u_int i = 0; i < oldContributions.size(); ++i) {
00306
00307 if(&oldContributions && !contribBuffer->Add(&oldContributions[i], norm)) {
00308 contribBuffer = film->scene->contribPool->Next(contribBuffer);
00309 contribBuffer->Add(&oldContributions[i], norm);
00310 }
00311 }
00312 }
00313
00314 weight = newWeight;
00315 LY = newLY;
00316 oldContributions = newContributions;
00317 sampleImage[0] = sample.imageX;
00318 sampleImage[1] = sample.imageY;
00319 sampleImage[2] = sample.lensU;
00320 sampleImage[3] = sample.lensV;
00321 sampleImage[4] = sample.time;
00322 sampleImage[5] = sample.wavelengths;
00323 sampleImage[6] = sample.singleWavelength;
00324 for (int i = SAMPLE_FLOATS; i < totalSamples; ++i)
00325 sampleImage[i] = sample.oneD[0][i - SAMPLE_FLOATS];
00326 for (int i = 0 ; i < totalTimes; ++i)
00327 timeImage[i] = sample.timexD[0][i];
00328 stamp = sample.stamp;
00329
00330 consecRejects = 0;
00331 } else {
00332
00333 const float norm = newWeight / (newLY / meanIntensity + pLarge);
00334 if (norm > 0.f) {
00335 for(u_int i = 0; i < newContributions.size(); ++i) {
00336
00337 if(!contribBuffer->Add(&newContributions[i], norm)) {
00338 contribBuffer = film->scene->contribPool->Next(contribBuffer);
00339 contribBuffer->Add(&newContributions[i], norm);
00340 }
00341 }
00342 }
00343
00344 for (int i = 0; i < totalTimes; ++i)
00345 sample.timexD[0][i] = timeImage[i];
00346 sample.stamp = stamp;
00347
00348 consecRejects++;
00349 }
00350 newContributions.clear();
00351 }
00352
00353 Sampler* MetropolisSampler::CreateSampler(const ParamSet ¶ms, const Film *film)
00354 {
00355 int xStart, xEnd, yStart, yEnd;
00356 film->GetSampleExtent(&xStart, &xEnd, &yStart, &yEnd);
00357 int maxConsecRejects = params.FindOneInt("maxconsecrejects", 512);
00358 float largeMutationProb = params.FindOneFloat("largemutationprob", 0.4f);
00359 float microMutationProb = params.FindOneFloat("micromutationprob", 0.f);
00360 float range = params.FindOneFloat("mutationrange", (xEnd - xStart + yEnd - yStart) / 32.);
00361 bool useVariance = params.FindOneBool("usevariance", false);
00362
00363 return new MetropolisSampler(xStart, xEnd, yStart, yEnd, maxConsecRejects,
00364 largeMutationProb, microMutationProb, range, useVariance);
00365 }
00366
00367 static DynamicLoader::RegisterSampler<MetropolisSampler> r("metropolis");