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
00032 using namespace lux;
00033
00034 #define SAMPLE_FLOATS 7
00035
00036
00037 static float mutate(const float x)
00038 {
00039 static const float s1 = 1.f / 1024.f, s2 = 1.f / 16.f;
00040 float randomValue = lux::random::floatValue();
00041 float dx = s2 * powf(s1 / s2, fabsf(2.f * randomValue - 1.f));
00042 if (randomValue < 0.5f) {
00043 float x1 = x + dx;
00044 return (x1 > 1.f) ? x1 - 1.f : x1;
00045 } else {
00046 float x1 = x - dx;
00047 return (x1 < 0.f) ? x1 + 1.f : x1;
00048 }
00049 }
00050
00051
00052 static float mutateScaled(const float x, const float mini, const float maxi, const float range)
00053 {
00054 float randomValue = lux::random::floatValue();
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 xStart, int xEnd, int yStart, int yEnd,
00067 int totMutations, float rng, int sw) :
00068 Sampler(xStart, xEnd, yStart, yEnd, 1), LY(0.), gain(0.f),
00069 totalSamples(0), totalTimes(0), totalMutations(totMutations), chain(0),
00070 numChains(0), mutation(0), consecRejects(0), stamp(0),
00071 range(rng), weight(0.), alpha(0.), baseImage(NULL), sampleImage(NULL),
00072 timeImage(NULL), strataWidth(sw)
00073 {
00074
00075 strataSamples = (float *)AllocAligned(2 * sw * sw * sizeof(float));
00076 strataSqr = sw*sw;
00077 currentStrata = strataSqr;
00078 }
00079
00080 ERPTSampler::~ERPTSampler() {
00081 FreeAligned(sampleImage);
00082 FreeAligned(baseImage);
00083 FreeAligned(timeImage);
00084 }
00085
00086
00087 ERPTSampler* ERPTSampler::clone() const
00088 {
00089 ERPTSampler *newSampler = new ERPTSampler(*this);
00090 newSampler->totalSamples = 0;
00091 newSampler->sampleImage = NULL;
00092 return newSampler;
00093 }
00094
00095 static void initERPT(ERPTSampler *sampler, const Sample *sample)
00096 {
00097 u_int i;
00098 sampler->normalSamples = SAMPLE_FLOATS;
00099 for (i = 0; i < sample->n1D.size(); ++i)
00100 sampler->normalSamples += sample->n1D[i];
00101 for (i = 0; i < sample->n2D.size(); ++i)
00102 sampler->normalSamples += 2 * sample->n2D[i];
00103 sampler->totalSamples = sampler->normalSamples;
00104 sampler->offset = new int[sample->nxD.size()];
00105 sampler->totalTimes = 0;
00106 for (i = 0; i < sample->nxD.size(); ++i) {
00107 sampler->offset[i] = sampler->totalSamples;
00108 sampler->totalTimes += sample->nxD[i];
00109 sampler->totalSamples += sample->dxD[i] * sample->nxD[i];
00110 }
00111 sampler->sampleImage = (float *)AllocAligned(sampler->totalSamples * sizeof(float));
00112 sampler->baseImage = (float *)AllocAligned(sampler->totalSamples * sizeof(float));
00113 sampler->timeImage = (int *)AllocAligned(sampler->totalTimes * sizeof(int));
00114 }
00115
00116
00117 bool ERPTSampler::GetNextSample(Sample *sample, u_int *use_pos)
00118 {
00119 sample->sampler = this;
00120 if (sampleImage == NULL) {
00121 initERPT(this, sample);
00122 }
00123
00124 if ((chain == 0 && mutation == 0) || initCount < initSamples) {
00125
00126
00127 if ((film->haltSamplePerPixel > 0) && film->enoughSamplePerPixel)
00128 return false;
00129
00130 if(currentStrata == strataSqr) {
00131
00132 StratifiedSample2D(strataSamples, strataWidth, strataWidth, true);
00133 Shuffle(strataSamples, strataSqr, 2);
00134 currentStrata = 0;
00135 }
00136
00137 sample->imageX = strataSamples[currentStrata*2] * (xPixelEnd - xPixelStart) + xPixelStart;
00138 sample->imageY = strataSamples[(currentStrata*2)+1] * (yPixelEnd - yPixelStart) + yPixelStart;
00139 currentStrata++;
00140
00141 sample->lensU = lux::random::floatValue();
00142 sample->lensV = lux::random::floatValue();
00143 sample->time = lux::random::floatValue();
00144 sample->wavelengths = lux::random::floatValue();
00145 sample->singleWavelength = lux::random::floatValue();
00146 for (int i = SAMPLE_FLOATS; i < normalSamples; ++i)
00147 sample->oneD[0][i - SAMPLE_FLOATS] = lux::random::floatValue();
00148 for (int i = 0; i < totalTimes; ++i)
00149 sample->timexD[0][i] = -1;
00150 sample->stamp = 0;
00151 } else {
00152 if (mutation == 0) {
00153
00154
00155 sample->imageX = baseImage[0];
00156 sample->imageY = baseImage[1];
00157 sample->lensU = baseImage[2];
00158 sample->lensV = baseImage[3];
00159 sample->time = baseImage[4];
00160 sample->wavelengths = baseImage[5];
00161 sample->singleWavelength = lux::random::floatValue();
00162 for (int i = SAMPLE_FLOATS; i < totalSamples; ++i)
00163 sample->oneD[0][i - SAMPLE_FLOATS] = baseImage[i];
00164 for (int i = 0; i < totalTimes; ++i) {
00165 if (sample->timexD[0][i] != -1)
00166 sample->timexD[0][i] = 0;
00167 }
00168 sample->stamp = 0;
00169 }
00170
00171
00172 sample->imageX = mutateScaled(sampleImage[0], xPixelStart, xPixelEnd, range);
00173 sample->imageY = mutateScaled(sampleImage[1], yPixelStart, yPixelEnd, range);
00174 sample->lensU = mutate(sampleImage[2]);
00175 sample->lensV = mutate(sampleImage[3]);
00176 sample->time = mutate(sampleImage[4]);
00177 sample->wavelengths = mutate(sampleImage[5]);
00178 sample->singleWavelength = lux::random::floatValue();
00179 for (int i = SAMPLE_FLOATS; i < normalSamples; ++i)
00180 sample->oneD[0][i - SAMPLE_FLOATS] = mutate(sampleImage[i]);
00181 ++(sample->stamp);
00182 }
00183
00184 return true;
00185 }
00186
00187 float *ERPTSampler::GetLazyValues(Sample *sample, u_int num, u_int pos)
00188 {
00189 float *data = sample->xD[num] + pos * sample->dxD[num];
00190 if (sample->timexD[num][pos] != sample->stamp) {
00191 if (sample->timexD[num][pos] == -1) {
00192 float *image = baseImage + offset[num] + pos * sample->dxD[num];
00193 for (u_int i = 0; i < sample->dxD[num]; ++i) {
00194 data[i] = lux::random::floatValue();
00195 image[i] = data[i];
00196 }
00197 sample->timexD[num][pos] = 0;
00198 } else {
00199 float *image = sampleImage + offset[num] + pos * sample->dxD[num];
00200 for (u_int i = 0; i < sample->dxD[num]; ++i){
00201 data[i] = image[i];
00202 }
00203 }
00204 for (int &time = sample->timexD[num][pos]; time < sample->stamp; ++time) {
00205 for (u_int i = 0; i < sample->dxD[num]; ++i)
00206 data[i] = mutate(data[i]);
00207 }
00208 }
00209 return data;
00210 }
00211
00212
00213 void ERPTSampler::AddSample(float imageX, float imageY, const Sample &sample, const Ray &ray, const XYZColor &newL, float newAlpha, int id)
00214 {
00215 if (!isSampleEnd) {
00216 sample.AddContribution(imageX, imageY, newL, newAlpha, id);
00217 } else {
00218 sample.contributions.clear();
00219 sample.AddContribution(imageX, imageY, newL, newAlpha, id);
00220 }
00221 }
00222
00223 void ERPTSampler::AddSample(const Sample &sample)
00224 {
00225 vector<Sample::Contribution> &newContributions(sample.contributions);
00226 float newLY = 0.0f;
00227 for(u_int i = 0; i < newContributions.size(); ++i)
00228 newLY += newContributions[i].color.y();
00229
00230 if (initCount < initSamples) {
00231 meanIntensity += newLY;
00232 ++(initCount);
00233 if (initCount < initSamples)
00234 return;
00235 if (meanIntensity == 0.) meanIntensity = 1.;
00236 meanIntensity /= initSamples * totalMutations;
00237 }
00238
00239 if (chain == 0 && mutation == 0) {
00240 gain = newLY / (meanIntensity * totalSamples);
00241 if (gain < 1.f)
00242 numChains = Floor2Int(lux::random::floatValue() + gain);
00243 else
00244 numChains = Floor2Int(gain);
00245 if (numChains == 0)
00246 return;
00247 gain /= numChains;
00248 film->AddSampleCount(1.f);
00249 }
00250
00251 float accProb = min(1.0f, newLY / LY);
00252 if (mutation == 0)
00253 accProb = 1.f;
00254 float newWeight = accProb;
00255 weight += 1.f - accProb;
00256
00257
00258 if (accProb == 1.f || lux::random::floatValue() < accProb) {
00259
00260 weight *= gain * meanIntensity / LY;
00261 if (!isinf(weight) && LY > 0.f) {
00262 for(u_int i = 0; i < oldContributions.size(); ++i) {
00263 XYZColor color = oldContributions[i].color;
00264 color *= weight;
00265 film->AddSample(oldContributions[i].imageX, oldContributions[i].imageY,
00266 color, oldContributions[i].alpha,
00267 oldContributions[i].buffer, oldContributions[i].bufferGroup);
00268 }
00269 }
00270
00271 weight = newWeight;
00272 LY = newLY;
00273 oldContributions = newContributions;
00274 sampleImage[0] = sample.imageX;
00275 sampleImage[1] = sample.imageY;
00276 sampleImage[2] = sample.lensU;
00277 sampleImage[3] = sample.lensV;
00278 sampleImage[4] = sample.time;
00279 sampleImage[5] = sample.wavelengths;
00280 sampleImage[6] = sample.singleWavelength;
00281 for (int i = SAMPLE_FLOATS; i < totalSamples; ++i)
00282 sampleImage[i] = sample.oneD[0][i - SAMPLE_FLOATS];
00283 for (int i = 0 ; i < totalTimes; ++i)
00284 timeImage[i] = sample.timexD[0][i];
00285 stamp = sample.stamp;
00286 if (chain == 0 && mutation == 0) {
00287 for (int i = 0; i < totalSamples; ++i)
00288 baseImage[i] = sampleImage[i];
00289 }
00290 consecRejects = 0;
00291 } else {
00292
00293 newWeight *= gain * meanIntensity / newLY;
00294 if (!isinf(newWeight) && newLY > 0.f) {
00295 for(u_int i = 0; i < newContributions.size(); ++i) {
00296 XYZColor color = newContributions[i].color;
00297 color *= newWeight;
00298 film->AddSample(newContributions[i].imageX, newContributions[i].imageY,
00299 color, newContributions[i].alpha,
00300 newContributions[i].buffer, newContributions[i].bufferGroup);
00301 }
00302 }
00303
00304 for (int i = 0; i < totalTimes; ++i)
00305 sample.timexD[0][i] = timeImage[i];
00306 sample.stamp = stamp;
00307 ++consecRejects;
00308 }
00309 if (++mutation >= totalMutations) {
00310 mutation = 0;
00311 if (++chain >= numChains)
00312 chain = 0;
00313 }
00314 }
00315
00316 Sampler* ERPTSampler::CreateSampler(const ParamSet ¶ms, const Film *film)
00317 {
00318 int xStart, xEnd, yStart, yEnd;
00319 film->GetSampleExtent(&xStart, &xEnd, &yStart, &yEnd);
00320 initSamples = params.FindOneInt("initsamples", 100000);
00321 initCount = 0;
00322 meanIntensity = 0.;
00323 int totMutations = params.FindOneInt("chainlength", 2000);
00324 float range = params.FindOneFloat("mutationrange", (xEnd - xStart + yEnd - yStart) / 50.);
00325 int stratawidth = params.FindOneInt("stratawidth", 256);
00326
00327 return new ERPTSampler(xStart, xEnd, yStart, yEnd, totMutations, range, stratawidth);
00328 }
00329
00330 int ERPTSampler::initCount, ERPTSampler::initSamples;
00331 float ERPTSampler::meanIntensity = 0.f;
00332