00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023 #include "memory.h"
00024
00025
00026
00027
00028
00029
00030
00031
00032 #include "metrosampler.h"
00033 #include "dynload.h"
00034
00035 using namespace lux;
00036
00037 #define SAMPLE_FLOATS 7
00038
00039
00040 static float mutate(const float x)
00041 {
00042 static const float s1 = 1/1024., s2 = 1/16.;
00043 float randomValue = lux::random::floatValue();
00044 float dx = s2 * exp(-log(s2/s1) * fabsf(2.f * randomValue - 1.f));
00045 if (randomValue < 0.5) {
00046 float x1 = x + dx;
00047 return (x1 > 1) ? x1 - 1 : x1;
00048 } else {
00049 float x1 = x - dx;
00050 return (x1 < 0) ? x1 + 1 : x1;
00051 }
00052 }
00053
00054
00055 static float mutateScaled(const float x, const float mini, const float maxi, const float range)
00056 {
00057 static const float s1 = 64.;
00058 float randomValue = lux::random::floatValue();
00059 float dx = range * exp(-log(s1) * fabsf(2.f * randomValue - 1.f));
00060 if (randomValue < 0.5) {
00061 float x1 = x + dx;
00062 return (x1 > maxi) ? x1 - maxi + mini : x1;
00063 } else {
00064 float x1 = x - dx;
00065 return (x1 < mini) ? x1 - mini + maxi : x1;
00066 }
00067 }
00068
00069
00070 MetropolisSampler::MetropolisSampler(int xStart, int xEnd, int yStart, int yEnd,
00071 int maxRej, float largeProb, float rng, int sw, bool useV) :
00072 Sampler(xStart, xEnd, yStart, yEnd, 1), large(true), LY(0.), V(0.),
00073 totalSamples(0), totalTimes(0), maxRejects(maxRej), consecRejects(0), stamp(0),
00074 pLarge(largeProb), range(rng), weight(0.), alpha(0.), sampleImage(NULL),
00075 timeImage(NULL), strataWidth(sw), useVariance(useV)
00076 {
00077
00078 strataSamples = (float *)AllocAligned(2 * sw * sw * sizeof(float));
00079 strataSqr = sw*sw;
00080 currentStrata = strataSqr;
00081 }
00082
00083 MetropolisSampler::~MetropolisSampler() {
00084 FreeAligned(sampleImage);
00085 FreeAligned(strataSamples);
00086 }
00087
00088
00089 MetropolisSampler* MetropolisSampler::clone() const
00090 {
00091 MetropolisSampler *newSampler = new MetropolisSampler(*this);
00092 newSampler->totalSamples = 0;
00093 newSampler->sampleImage = NULL;
00094
00095 return newSampler;
00096 }
00097
00098 static void initMetropolis(MetropolisSampler *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 = (float *)AllocAligned(sampler->totalSamples * sizeof(float));
00115 sampler->timeImage = (int *)AllocAligned(sampler->totalTimes * sizeof(int));
00116 }
00117
00118
00119 bool MetropolisSampler::GetNextSample(Sample *sample, u_int *use_pos)
00120 {
00121 sample->sampler = this;
00122 large = (lux::random::floatValue() < pLarge) || initCount < initSamples;
00123 if (sampleImage == NULL) {
00124 initMetropolis(this, sample);
00125 large = true;
00126 }
00127
00128
00129
00130 if ((film->haltSamplePerPixel > 0) && film->enoughSamplePerPixel)
00131 return false;
00132
00133 if (large) {
00134 if(currentStrata == strataSqr) {
00135
00136 StratifiedSample2D(strataSamples, strataWidth, strataWidth, true);
00137 Shuffle(strataSamples, strataSqr, 2);
00138 currentStrata = 0;
00139 }
00140
00141
00142 sample->imageX = strataSamples[currentStrata*2] * (xPixelEnd - xPixelStart) + xPixelStart;
00143 sample->imageY = strataSamples[(currentStrata*2)+1] * (yPixelEnd - yPixelStart) + yPixelStart;
00144 currentStrata++;
00145
00146 sample->lensU = lux::random::floatValue();
00147 sample->lensV = lux::random::floatValue();
00148 sample->time = lux::random::floatValue();
00149 sample->wavelengths = lux::random::floatValue();
00150 sample->singleWavelength = lux::random::floatValue();
00151 for (int i = SAMPLE_FLOATS; i < normalSamples; ++i)
00152 sample->oneD[0][i - SAMPLE_FLOATS] = lux::random::floatValue();
00153 for (int i = 0; i < totalTimes; ++i)
00154 sample->timexD[0][i] = -1;
00155 sample->stamp = 0;
00156 } else {
00157
00158
00159 sample->imageX = mutateScaled(sampleImage[0], xPixelStart, xPixelEnd, range);
00160 sample->imageY = mutateScaled(sampleImage[1], yPixelStart, yPixelEnd, range);
00161 sample->lensU = mutate(sampleImage[2]);
00162 sample->lensV = mutate(sampleImage[3]);
00163 sample->time = mutate(sampleImage[4]);
00164 sample->wavelengths = mutate(sampleImage[5]);
00165 sample->singleWavelength = lux::random::floatValue();
00166 for (int i = SAMPLE_FLOATS; i < normalSamples; ++i)
00167 sample->oneD[0][i - SAMPLE_FLOATS] = mutate(sampleImage[i]);
00168 ++(sample->stamp);
00169 }
00170
00171 return true;
00172 }
00173
00174 float *MetropolisSampler::GetLazyValues(Sample *sample, u_int num, u_int pos)
00175 {
00176 float *data = sample->xD[num] + pos * sample->dxD[num];
00177 if (sample->timexD[num][pos] != sample->stamp) {
00178 if (sample->timexD[num][pos] == -1) {
00179 for (u_int i = 0; i < sample->dxD[num]; ++i)
00180 data[i] = lux::random::floatValue();
00181 sample->timexD[num][pos] = 0;
00182 } else {
00183 for (u_int i = 0; i < sample->dxD[num]; ++i){
00184 data[i] = sampleImage[offset[num] +
00185 pos * sample->dxD[num] + i];
00186 }
00187 }
00188 for (; sample->timexD[num][pos] < sample->stamp; ++(sample->timexD[num][pos])) {
00189 for (u_int i = 0; i < sample->dxD[num]; ++i)
00190 data[i] = mutate(data[i]);
00191 }
00192 }
00193 return data;
00194 }
00195
00196
00197 void MetropolisSampler::AddSample(float imageX, float imageY, const Sample &sample, const Ray &ray, const XYZColor &newL, float newAlpha, int id)
00198 {
00199 if (!isSampleEnd) {
00200 sample.AddContribution(imageX, imageY, newL, newAlpha, id);
00201 } else {
00202 sample.contributions.clear();
00203 sample.AddContribution(imageX, imageY, newL, newAlpha, id);
00204 }
00205 }
00206
00207 void MetropolisSampler::AddSample(const Sample &sample)
00208 {
00209 vector<Sample::Contribution> &newContributions(sample.contributions);
00210 float newLY = 0.f, newV = 0.f;
00211 for(u_int i = 0; i < newContributions.size(); ++i) {
00212 newLY += newContributions[i].color.y();
00213 if (newContributions[i].color.y() > 0.f)
00214 newV += newContributions[i].color.y() * newContributions[i].variance;
00215 }
00216
00217 if (initCount < initSamples) {
00218 meanIntensity += newLY / initSamples;
00219 ++(initCount);
00220 if (initCount < initSamples)
00221 return;
00222 if (meanIntensity <= 0.f) meanIntensity = 1.f;
00223 }
00224
00225 film->AddSampleCount(1.f);
00226
00227 float accProb, accProb2, factor;
00228 if (LY > 0.f) {
00229 accProb = min(1.f, newLY / LY);
00230 if (useVariance && V > 0.f && newV > 0.f && newLY > 0.f)
00231 factor = newV / (V * accProb);
00232 else
00233 factor = 1.f;
00234 } else {
00235 accProb = 1.f;
00236 factor = 1.f;
00237 }
00238 accProb2 = accProb * factor;
00239 float newWeight = (accProb + (large ? 1.f : 0.f)) / (factor * newLY / meanIntensity + pLarge);
00240 weight += (1.f - accProb) / (LY / (factor * meanIntensity) + pLarge);
00241
00242 if (accProb2 == 1.f || consecRejects >= maxRejects || lux::random::floatValue() < accProb2) {
00243
00244 for(u_int i = 0; i < oldContributions.size(); ++i) {
00245 XYZColor color = oldContributions[i].color;
00246 color *= weight;
00247 film->AddSample(oldContributions[i].imageX, oldContributions[i].imageY,
00248 color, oldContributions[i].alpha,
00249 oldContributions[i].buffer, oldContributions[i].bufferGroup);
00250 }
00251
00252 weight = newWeight;
00253 LY = newLY;
00254 V = newV;
00255 oldContributions = newContributions;
00256 sampleImage[0] = sample.imageX;
00257 sampleImage[1] = sample.imageY;
00258 sampleImage[2] = sample.lensU;
00259 sampleImage[3] = sample.lensV;
00260 sampleImage[4] = sample.time;
00261 sampleImage[5] = sample.wavelengths;
00262 sampleImage[6] = sample.singleWavelength;
00263 for (int i = SAMPLE_FLOATS; i < totalSamples; ++i)
00264 sampleImage[i] = sample.oneD[0][i - SAMPLE_FLOATS];
00265 for (int i = 0 ; i < totalTimes; ++i)
00266 timeImage[i] = sample.timexD[0][i];
00267 stamp = sample.stamp;
00268
00269 consecRejects = 0;
00270 } else {
00271
00272 for(u_int i = 0; i < newContributions.size(); ++i) {
00273 XYZColor color = newContributions[i].color;
00274 color *= newWeight;
00275 film->AddSample(newContributions[i].imageX, newContributions[i].imageY,
00276 color, newContributions[i].alpha,
00277 newContributions[i].buffer, newContributions[i].bufferGroup);
00278 }
00279
00280 for (int i = 0; i < totalTimes; ++i)
00281 sample.timexD[0][i] = timeImage[i];
00282 sample.stamp = stamp;
00283
00284 consecRejects++;
00285 }
00286 }
00287
00288 Sampler* MetropolisSampler::CreateSampler(const ParamSet ¶ms, const Film *film)
00289 {
00290 int xStart, xEnd, yStart, yEnd;
00291 film->GetSampleExtent(&xStart, &xEnd, &yStart, &yEnd);
00292 initSamples = params.FindOneInt("initsamples", 100000);
00293 initCount = 0;
00294 meanIntensity = 0.;
00295 int maxConsecRejects = params.FindOneInt("maxconsecrejects", 512);
00296 float largeMutationProb = params.FindOneFloat("largemutationprob", 0.4f);
00297 float range = params.FindOneFloat("mutationrange", (xEnd - xStart + yEnd - yStart) / 32.);
00298 int stratawidth = params.FindOneInt("stratawidth", 256);
00299 bool useVariance = params.FindOneBool("usevariance", false);
00300
00301 return new MetropolisSampler(xStart, xEnd, yStart, yEnd, maxConsecRejects,
00302 largeMutationProb, range, stratawidth, useVariance);
00303 }
00304
00305 int MetropolisSampler::initCount, MetropolisSampler::initSamples;
00306 float MetropolisSampler::meanIntensity;