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
00027
00028
00029 #include "erpt.h"
00030 #include "dynload.h"
00031 #include "scene.h"
00032 #include "error.h"
00033
00034 using namespace lux;
00035
00036 #define SAMPLE_FLOATS 7
00037
00038
00039 static float mutate(const float x, const float randomValue)
00040 {
00041 static const float s1 = 1.f / 1024.f, s2 = 1.f / 16.f;
00042 float dx = s2 * powf(s1 / s2, fabsf(2.f * randomValue - 1.f));
00043 if (randomValue < 0.5f) {
00044 float x1 = x + dx;
00045 return (x1 > 1.f) ? x1 - 1.f : x1;
00046 } else {
00047 float x1 = x - dx;
00048 return (x1 < 0.f) ? x1 + 1.f : x1;
00049 }
00050 }
00051
00052
00053 static float mutateScaled(const float x, const float randomValue, const float mini, const float maxi, const float range)
00054 {
00055 float dx = range * exp(-log(2.f * range) * fabsf(2.f * randomValue - 1.f));
00056 if (randomValue < 0.5f) {
00057 float x1 = x + dx;
00058 return (x1 > maxi) ? x1 - maxi + mini : x1;
00059 } else {
00060 float x1 = x - dx;
00061 return (x1 < mini) ? x1 - mini + maxi : x1;
00062 }
00063 }
00064
00065
00066 ERPTSampler::ERPTSampler(int totMutations, float microProb, float rng,
00067 Sampler *sampler) :
00068 Sampler(sampler->xPixelStart, sampler->xPixelEnd,
00069 sampler->yPixelStart, sampler->yPixelEnd, sampler->samplesPerPixel),
00070 normalSamples(0), totalSamples(0), totalTimes(0), totalMutations(totMutations),
00071 pMicro(microProb), range(rng), baseSampler(sampler),
00072 baseImage(NULL), sampleImage(NULL), currentImage(NULL),
00073 baseTimeImage(NULL), timeImage(NULL), currentTimeImage(NULL), offset(NULL),
00074 numChains(0), chain(0), mutation(-1), stamp(0), numMicro(-1), posMicro(-1),
00075 baseLY(0.f), quantum(0.f), weight(0.f), LY(0.f), alpha(0.f),
00076 totalLY(0.), sampleCount(0.)
00077 {
00078 }
00079
00080 ERPTSampler::~ERPTSampler() {
00081 FreeAligned(sampleImage);
00082 FreeAligned(baseImage);
00083 FreeAligned(timeImage);
00084 FreeAligned(baseTimeImage);
00085 delete baseSampler;
00086 }
00087
00088
00089 ERPTSampler* ERPTSampler::clone() const
00090 {
00091 ERPTSampler *newSampler = new ERPTSampler(*this);
00092 newSampler->baseSampler = baseSampler->clone();
00093 newSampler->totalSamples = 0;
00094 newSampler->sampleImage = NULL;
00095 return newSampler;
00096 }
00097
00098 static void initERPT(ERPTSampler *sampler, const Sample *sample)
00099 {
00100 u_int i;
00101 sampler->normalSamples = SAMPLE_FLOATS;
00102 for (i = 0; i < sample->n1D.size(); ++i)
00103 sampler->normalSamples += sample->n1D[i];
00104 for (i = 0; i < sample->n2D.size(); ++i)
00105 sampler->normalSamples += 2 * sample->n2D[i];
00106 sampler->totalSamples = sampler->normalSamples;
00107 sampler->offset = new int[sample->nxD.size()];
00108 sampler->totalTimes = 0;
00109 for (i = 0; i < sample->nxD.size(); ++i) {
00110 sampler->offset[i] = sampler->totalSamples;
00111 sampler->totalTimes += sample->nxD[i];
00112 sampler->totalSamples += sample->dxD[i] * sample->nxD[i];
00113 }
00114 sampler->sampleImage = AllocAligned<float>(sampler->totalSamples);
00115 sampler->baseImage = AllocAligned<float>(sampler->totalSamples);
00116 sampler->timeImage = AllocAligned<int>(sampler->totalTimes);
00117 sampler->baseTimeImage = AllocAligned<int>(sampler->totalTimes);
00118 sampler->baseSampler->SetTsPack(sampler->tspack);
00119 sampler->baseSampler->SetFilm(sampler->film);
00120 sampler->mutation = -1;
00121
00122
00123 sampler->contribBuffer = sampler->film->scene->contribPool->Next(NULL);
00124 }
00125
00126
00127 bool ERPTSampler::GetNextSample(Sample *sample, u_int *use_pos)
00128 {
00129 sample->sampler = this;
00130 if (sampleImage == NULL) {
00131 initERPT(this, sample);
00132 }
00133
00134 if (mutation == -1) {
00135
00136
00137 if (film->enoughSamplePerPixel)
00138 return false;
00139
00140 const bool ret = baseSampler->GetNextSample(sample, use_pos);
00141 sample->sampler = this;
00142 for (int i = 0; i < totalTimes; ++i)
00143 sample->timexD[0][i] = -1;
00144 sample->stamp = 0;
00145 currentImage = baseImage;
00146 currentTimeImage = baseTimeImage;
00147 stamp = 0;
00148 numMicro = -1;
00149 return ret;
00150 } else {
00151 if (mutation == 0) {
00152
00153 for (int i = 0; i < totalTimes; ++i)
00154 sample->timexD[0][i] = baseTimeImage[i];
00155 sample->stamp = 0;
00156 currentImage = baseImage;
00157 currentTimeImage = baseTimeImage;
00158 stamp = 0;
00159 }
00160
00161
00162 float mutationSelector = tspack->rng->floatValue();
00163 if (mutationSelector < pMicro)
00164 numMicro = min<int>(sample->nxD.size(), Float2Int(mutationSelector / pMicro * (sample->nxD.size() + 1)));
00165 else
00166 numMicro = -1;
00167 if (numMicro > 0) {
00168 u_int maxPos = 0;
00169 for(; maxPos < sample->nxD[numMicro - 1]; ++maxPos)
00170 if (sample->timexD[numMicro - 1][maxPos] < sample->stamp)
00171 break;
00172 posMicro = min<int>(sample->nxD[numMicro - 1] - 1, Float2Int(tspack->rng->floatValue() * maxPos));
00173 } else {
00174 posMicro = -1;
00175 sample->imageX = mutateScaled(currentImage[0], tspack->rng->floatValue(), xPixelStart, xPixelEnd, range);
00176 sample->imageY = mutateScaled(currentImage[1], tspack->rng->floatValue(), yPixelStart, yPixelEnd, range);
00177 sample->lensU = mutate(currentImage[2], tspack->rng->floatValue());
00178 sample->lensV = mutate(currentImage[3], tspack->rng->floatValue());
00179 sample->time = mutate(currentImage[4], tspack->rng->floatValue());
00180 sample->wavelengths = mutate(currentImage[5], tspack->rng->floatValue());
00181 sample->singleWavelength = mutateScaled(currentImage[6], tspack->rng->floatValue(), 0.f, 1.f, 1.f);
00182 for (int i = SAMPLE_FLOATS; i < normalSamples; ++i)
00183 sample->oneD[0][i - SAMPLE_FLOATS] = mutate(currentImage[i], tspack->rng->floatValue());
00184 }
00185 ++(sample->stamp);
00186 }
00187
00188 return true;
00189 }
00190
00191 float *ERPTSampler::GetLazyValues(Sample *sample, u_int num, u_int pos)
00192 {
00193 const u_int size = sample->dxD[num];
00194 float *data = sample->xD[num] + pos * size;
00195 int stampLimit = sample->stamp;
00196 if (numMicro >= 0 && num != u_int(numMicro - 1) && pos != u_int(posMicro))
00197 --stampLimit;
00198 if (sample->timexD[num][pos] != stampLimit) {
00199 if (sample->timexD[num][pos] == -1) {
00200 baseSampler->GetLazyValues(sample, num, pos);
00201 sample->timexD[num][pos] = 0;
00202 } else {
00203 int start = offset[num] + pos * size;
00204 float *image = currentImage + start;
00205 for (u_int i = 0; i < size; ++i)
00206 data[i] = image[i];
00207 }
00208 for (int &time = sample->timexD[num][pos]; time < stampLimit; ++time) {
00209 for (u_int i = 0; i < size; ++i)
00210 data[i] = mutate(data[i], tspack->rng->floatValue());
00211 }
00212 }
00213 return data;
00214 }
00215
00216
00217 void ERPTSampler::AddSample(const Sample &sample)
00218 {
00219 vector<Contribution> &newContributions(sample.contributions);
00220 float newLY = 0.0f;
00221 for(u_int i = 0; i < newContributions.size(); ++i)
00222 newLY += newContributions[i].color.Y();
00223 if (mutation <= 0) {
00224 if (weight > 0.f) {
00225
00226 weight *= quantum / LY;
00227 if (!isinf(weight) && LY > 0.f) {
00228 for(u_int i = 0; i < oldContributions.size(); ++i) {
00229
00230 if(&oldContributions && !contribBuffer->Add(&oldContributions[i], weight)) {
00231 contribBuffer = film->scene->contribPool->Next(contribBuffer);
00232 contribBuffer->Add(&oldContributions[i], weight);
00233 }
00234 }
00235 }
00236 weight = 0.f;
00237 }
00238 if (mutation == -1) {
00239 contribBuffer->AddSampleCount(1.f);
00240 ++sampleCount;
00241 if (!(newLY > 0.f)) {
00242 newContributions.clear();
00243 return;
00244 }
00245 totalLY += newLY;
00246 const float meanIntensity = totalLY > 0.f ? totalLY / sampleCount : 1.f;
00247
00248 quantum = newLY / meanIntensity;
00249 numChains = max(1, Floor2Int(quantum + .5f));
00250 if (numChains > 100) {
00251 printf("%d chains -> %d\n", numChains, totalMutations);
00252 numChains = totalMutations;
00253 }
00254 quantum /= (numChains * totalSamples);
00255 baseLY = newLY;
00256 baseContributions = newContributions;
00257 baseImage[0] = sample.imageX;
00258 baseImage[1] = sample.imageY;
00259 baseImage[2] = sample.lensU;
00260 baseImage[3] = sample.lensV;
00261 baseImage[4] = sample.time;
00262 baseImage[5] = sample.wavelengths;
00263 baseImage[6] = sample.singleWavelength;
00264 for (int i = SAMPLE_FLOATS; i < totalSamples; ++i)
00265 baseImage[i] = sample.oneD[0][i - SAMPLE_FLOATS];
00266 for (int i = 0 ; i < totalTimes; ++i)
00267 baseTimeImage[i] = sample.timexD[0][i];
00268 mutation = 0;
00269 newContributions.clear();
00270 return;
00271 }
00272 LY = baseLY;
00273 oldContributions = baseContributions;
00274 }
00275
00276 float accProb;
00277 if (LY > 0.f)
00278 accProb = min(1.f, newLY / LY);
00279 else
00280 accProb = 1.f;
00281 float newWeight = accProb;
00282 weight += 1.f - accProb;
00283
00284
00285 if (accProb == 1.f || tspack->rng->floatValue() < accProb) {
00286
00287 weight *= quantum / LY;
00288 if (!isinf(weight) && LY > 0.f) {
00289 for(u_int i = 0; i < oldContributions.size(); ++i) {
00290
00291 if(&oldContributions && !contribBuffer->Add(&oldContributions[i], weight)) {
00292 contribBuffer = film->scene->contribPool->Next(contribBuffer);
00293 contribBuffer->Add(&oldContributions[i], weight);
00294 }
00295 }
00296 }
00297 weight = newWeight;
00298 LY = newLY;
00299 oldContributions = newContributions;
00300
00301
00302 sampleImage[0] = sample.imageX;
00303 sampleImage[1] = sample.imageY;
00304 sampleImage[2] = sample.lensU;
00305 sampleImage[3] = sample.lensV;
00306 sampleImage[4] = sample.time;
00307 sampleImage[5] = sample.wavelengths;
00308 sampleImage[6] = sample.singleWavelength;
00309 for (int i = SAMPLE_FLOATS; i < totalSamples; ++i)
00310 sampleImage[i] = sample.oneD[0][i - SAMPLE_FLOATS];
00311 for (int i = 0 ; i < totalTimes; ++i)
00312 timeImage[i] = sample.timexD[0][i];
00313 stamp = sample.stamp;
00314 currentImage = sampleImage;
00315 currentTimeImage = timeImage;
00316 } else {
00317
00318 newWeight *= quantum / newLY;
00319 if (!isinf(newWeight) && newLY > 0.f) {
00320 for(u_int i = 0; i < newContributions.size(); ++i) {
00321
00322 if(!contribBuffer->Add(&newContributions[i], newWeight)) {
00323 contribBuffer = film->scene->contribPool->Next(contribBuffer);
00324 contribBuffer->Add(&newContributions[i], newWeight);
00325 }
00326 }
00327 }
00328
00329
00330 for (int i = 0; i < totalTimes; ++i)
00331 sample.timexD[0][i] = currentTimeImage[i];
00332 sample.stamp = stamp;
00333 }
00334 if (++mutation >= totalMutations) {
00335 mutation = 0;
00336 if (++chain >= numChains) {
00337 chain = 0;
00338 mutation = -1;
00339 }
00340 }
00341 newContributions.clear();
00342 }
00343
00344 Sampler* ERPTSampler::CreateSampler(const ParamSet ¶ms, const Film *film)
00345 {
00346 int xStart, xEnd, yStart, yEnd;
00347 film->GetSampleExtent(&xStart, &xEnd, &yStart, &yEnd);
00348 int totMutations = params.FindOneInt("chainlength", 100);
00349 float microMutationProb = params.FindOneFloat("micromutationprob", .5f);
00350 float range = params.FindOneFloat("mutationrange", (xEnd - xStart + yEnd - yStart) / 50.);
00351 string base = params.FindOneString("basesampler", "random");
00352 Sampler *sampler = MakeSampler(base, params, film);
00353 if (sampler == NULL) {
00354 luxError(LUX_SYSTEM, LUX_SEVERE, "ERPTSampler: Could not obtain a valid sampler");
00355 return NULL;
00356 }
00357 return new ERPTSampler(totMutations, microMutationProb, range, sampler);
00358 }
00359
00360 static DynamicLoader::RegisterSampler<ERPTSampler> r("erpt");