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 "metropolis.h"
00030
00031 using namespace lux;
00032
00033
00034 float MetroSample::mutate (const float x) const {
00035 static const float s1 = 1/1024., s2 = 1/16.;
00036 float dx = s2 * exp(-log(s2/s1) * lux::random::floatValue());
00037 if (lux::random::floatValue() < 0.5) {
00038 float x1 = x + dx;
00039 return (x1 > 1) ? x1 - 1 : x1;
00040 } else {
00041 float x1 = x - dx;
00042 return (x1 < 0) ? x1 + 1 : x1;
00043 }
00044 }
00045
00046
00047
00048 float MetroSample::value (const int i, float defval) const {
00049 if (i >= (int) values.size()) {
00050 int oldsize = values.size();
00051 int newsize = i + 1;
00052 values.resize(newsize);
00053 modify.resize(newsize);
00054 for (int j = oldsize; j < newsize; j++) {
00055 if(defval == -1.)
00056 values[j] = lux::random::floatValue();
00057 else
00058 values[j] = defval;
00059 modify[j] = time;
00060 }
00061 }
00062 for (; modify[i] < time; modify[i]++)
00063 values[i] = mutate(values[i]);
00064 return values[i];
00065 }
00066
00067
00068 MetroSample MetroSample::next () const {
00069 MetroSample newsamp(*this);
00070 newsamp.time = time + 1;
00071 return newsamp;
00072 }
00073
00074
00075 void Metropolis::GetNext(float& bs1, float& bs2, float& bcs, int pathLength)
00076 {
00077
00078 bs1 = newsamp.value(5+(3*pathLength), bs1);
00079 bs2 = newsamp.value(5+(3*pathLength)+1, bs2);
00080 bcs = newsamp.value(5+(3*pathLength)+2, bcs);
00081 }
00082
00083
00084 bool Metropolis::GetNextSample(Sampler *sampler, Sample *sample, u_int *use_pos)
00085 {
00086 large = (lux::random::floatValue() < pLarge);
00087 if(large) {
00088
00089 newsamp = MetroSample();
00090
00091
00092 if(!sampler->GetNextSample(sample, use_pos))
00093 return false;
00094
00095
00096 newsamp.value(0, (sample->imageX - xStart) / (xEnd - xStart));
00097 newsamp.value(1, (sample->imageY - yStart) / (yEnd - yStart));
00098 newsamp.value(2, sample->lensU);
00099 newsamp.value(3, sample->lensV);
00100 newsamp.value(4, sample->time);
00101 } else {
00102
00103 newsamp = msamp.next();
00104
00105
00106 sample->imageX = newsamp.value(0, -1.) * (xEnd - xStart) + xStart;
00107 sample->imageY = newsamp.value(1, -1.) * (yEnd - yStart) + yStart;
00108 sample->lensU = newsamp.value(2, -1.);
00109 sample->lensV = newsamp.value(3, -1.);
00110 sample->time = newsamp.value(4, -1.);
00111 }
00112
00113 return true;
00114 }
00115
00116
00117 void Metropolis::AddSample(const Sample &sample, const Ray &ray,
00118 const Spectrum &newL, float alpha, Film *film)
00119 {
00120
00121 float accprob = min(1.0f, newL.y()/L.y());
00122
00123
00124 if (L.y() > 0.f && L != Spectrum(0.f))
00125 film->AddSample(msamp.value(0, -1.) * (xEnd - xStart) + xStart,
00126 msamp.value(1, -1.) * (yEnd - yStart) + yStart,
00127 L * (1. / L.y()) * (1. - accprob), alpha);
00128
00129
00130 if (newL.y() > 0.f && newL != Spectrum(0.f))
00131 film->AddSample(newsamp.value(0, -1.) * (xEnd - xStart) + xStart,
00132 newsamp.value(1, -1.) * (yEnd - yStart) + yStart,
00133 newL * (1. / newL.y()) * accprob, alpha);
00134
00135
00136 if (lux::random::floatValue() < accprob || consec_rejects > maxReject) {
00137 msamp = newsamp;
00138 L = newL;
00139 consec_rejects = 0;
00140 } else {
00141 consec_rejects++;
00142 }
00143 }
00144