Now, I'm implementing Veach-style BPT.
And I could get a result which converge to the reference image generated by using path tracing.
But this is only in case that I don't use Russian roulette in BPT.
For debugging BPT, of course I checked the equality between PT and BPT's unweighted contributions with respect to each path length.
For example, the result of each
- 1 bounce light transport by [reference PT] and [s=0, t=3], [s=1, t=2] are the same, (Now I consider only the pinhole camera model, t<=1 paths are not considered)
- 2 bounce light transport by [reference PT] and [s=0, t=4], [s=1, t=3], [s=2, t=2] are the same,
- 3 bounce ...
...
But I say again, these results are gotten only if I don't use Russian roulette.
Of course, I know that it is required that the probability for generating a adjacent vertex in MIS weight computations contain RR continuation probability.
Would anyone give me some suggestions ?
I attached three images that only represent 2 bounce transport, the first is the reference by PT, the second is the BPT which doesn't use RR, the third is the BPT which uses RR (strictly, minimum subpath lengths are set 3).
The first twos look the same, but the third is not, looks slightly brighter.
I also attached the current code fragments just to make sure.
I'm sorry my poor English.
Thank you.
NOTICE:
consider only the pinhole camera model.
consider only diffuse surfaces.
Now I don't distinguish the shading normal from the geometry normal, but in the test scene, they are the same vectors.
Code: Select all
struct ScatteringEvent {
Object* obj;//----reference of an object on which a subpath vertex is
objID oid;//------collision object, triangle identifier
BSDF* bsdf;//-----reference of a BSDF
Vector3f pos;//---position
Vector3f gn;//----geometry normal
Vector3f sn;//----shading normal
Color weight;//---subpath weight
float dirIPDF;//--directional PDF for generating adjacent vertex toward reverse direction of light flow (solid angle)
float dirOPDF;//--directional PDF for generating adjacent vertex toward direction of light flow (solid angle)
float dAIPDF;//---PDF for generating adjacent vertex toward reverse direction of light flow (surface area)
float dAOPDF;//---PDF for generating adjacent vertex toward direction of light flow (surface area)
float rrIPDF;//---path continuation probability for generating adjacent vertex toward reverse direction of light flow
float rrOPDF;//---path continuation probability for generating adjacent vertex toward direction of light flow
Color bsdfValue;//BSDF value
ScatteringEvent() { };
ScatteringEvent(Object* obj_, BSDF* bsdf_, const objID &oid_,
const Vector3f &pos_, const Vector3f &gn_, const Vector3f &sn_,
const Color &weight_,
const float dirIPDF_, const float dirOPDF_) :
obj(obj_), bsdf(bsdf_), oid(oid_),
pos(pos_), gn(gn_), sn(sn_),
weight(weight_),
dirIPDF(dirIPDF_), dirOPDF(dirOPDF_)
{
dAIPDF = dAOPDF = rrIPDF = rrOPDF = -1;
};
};
Code: Select all
void BidirectionalPathTracingWorker::job() {
BidirectionalPathTracingTaskData& pttask = (BidirectionalPathTracingTaskData&)*_task;
BidirectionalPathTracing* op = (BidirectionalPathTracing*)pttask.p;
Scene* sceneRef = op->_TargetScene;
ICamera* camera = sceneRef->Camera;
for (int curImgPX = 0; curImgPX < op->_width; ++curImgPX) {
//------------------------------------------------------------------------
//calculate the index of the current pixel, ray direction.
int cIdx = pttask.iRow * op->_width + curImgPX;
float curImgX = (curImgPX + op->_rand->GetRandom32bitReal0cTo1c()) / op->_width;
float curImgY = (pttask.iRow + op->_rand->GetRandom32bitReal0cTo1c()) / op->_height;
float pixelPDF = 1.0f / (1.0f * 1.0f);
Vector3f rayDir = camera->GetRayDirection(op->_cliprect.orgX + curImgX * op->_cliprect.sizeW,
op->_cliprect.orgY + curImgY * op->_cliprect.sizeH);
Vector3f camDir = camera->GetDirection();
//------------------------------------------------------------------------
static const int LimitPathLength = 50;
ScatteringEvent lightPath[LimitPathLength + 2];
ScatteringEvent eyePath[LimitPathLength + 2];
//set the minimum path lengths. if it is traced more than this values, a subpath is terminated rondomly by russian roulette.
static const int numAbsLengthLight = 1;
static const int numAbsLengthEye = 1;
int nL = 0;
int nE = 0;
float dirPDF;
float costh;
Color weight;
Intersection isect;
objID collisionID;
Ray ray;
//--------------------------------------------------------------------------------------------------------------
//generate light subpath.
//PA(yi) : probability density for generating i-th vertex with respect to surface area
//Ps(y' -> yi) : probability density for generating i-th vertex from y', measured with respect to solid angle
//Psp(y' -> yi) : probability density for generating i-th vertex from y', measured with respect to projected solid angle
//choose a light and generate a ray.
float sumLightPower = 0;
int whichLight = 0;
for (int i = 0; i < sceneRef->NumOfLights(); ++i) {
sumLightPower += sceneRef->GetLight(i)->emissionPower();
}
float rLight = sumLightPower * op->_rand->GetRandom32bitReal0cTo1o();
float acmLightPower = 0;
for (int i = 0; i < sceneRef->NumOfLights(); ++i) {
acmLightPower += sceneRef->GetLight(i)->emissionPower();
if (acmLightPower >= rLight) {
whichLight = i;
break;
}
}
Object* light = sceneRef->GetLight(whichLight);
float lightPDF = light->emissionPower() / sumLightPower;
float areaPDF;
Vector3f lightGN;
Vector3f lightSN;
Color radiosity;
Color Le = light->sampleRadiance(sampledWLs, &ray, &lightGN, &lightSN, &radiosity, &areaPDF, &dirPDF);
//light - 0th vertex
weight = Color(1.0f, 1.0f, 1.0f);// a0 = 1
lightPath[nL++] = ScatteringEvent(NULL, NULL, objID(),
vz, vz, vz,
weight,
-1, -1);
//light - 1st vertex
weight = radiosity / (areaPDF * lightPDF);// a1 = Le^(0)(y0) / PA(y0)
lightPath[nL++] = ScatteringEvent(light, light->materials[0]->bsdf, objID(sceneRef->LightID(whichLight), 0),
ray.org, lightGN, lightSN,
weight,
-1, dirPDF);
lightPath[1].bsdfValue = Le / radiosity;// Le^(1)(y0 -> y1)
//trace light subpath.
weight = Le / (dirPDF * areaPDF * lightPDF) * fabsf(Vector3f::Dot(ray.dir, lightSN));// a2 = Le^(1)(y0 -> y1) / Psp(y0 -> y1) * a1
collisionID = lightPath[1].oid;
while (FindNearestIntersectionWithRay(&isect, ray, collisionID, sceneRef)) {
Vector3f hitPos = ray.org + ray.dir * isect.t;
Vector3f vin = -ray.dir;
Material* mat = isect.obj->materials[0];
collisionID = isect.objID;
lightPath[nL] = ScatteringEvent(isect.obj, mat->bsdf, isect.objID,
hitPos, isect.gn, isect.sn,
weight,
-1, -1);
//------------------------------------------------------------------------
//reach the depth limit.
if (nL > LimitPathLength) {
break;
}
//------------------------------------------------------------------------
bool pathContinue = false;
//------------------------------------------------------------------------
//sample transport direction.
//terminate path by russian roulette, except for the first few vertices.
Vector3f vout;
BSDFType sampledType;
BSDFSample bsdfSample = BSDFSample(op->_rand->GetRandom32bitReal0cTo1o(),
op->_rand->GetRandom32bitReal0cTo1o());
Color r = mat->bsdf->sample_f(vin, isect.gn, isect.sn, bsdfSample, &vout, &dirPDF, BSDF_ALL, &sampledType);
costh = fabsf(Vector3f::Dot(vout, isect.sn));
Color fraction = r * costh / dirPDF;
float p = 0.5f;//fminf(1.0f, fraction.maxComp());
if (costh != 0) {
if (ray.depth < numAbsLengthLight) {
weight *= fraction;// a_{i} = fs(y_{i-1} -> y_{i} -> y_{i+1}) / Psp(y_{i} -> y_{i+1}) * a_{i-1}
if (weight.hasNonZero()) {
ray = Ray(hitPos, vout, ray.depth + 1);
pathContinue = true;
}
}
else {
if (op->_rand->GetRandom32bitReal0cTo1o() < p) {
weight *= fraction / p;// a_{i} = fs(y_{i-1} -> y_{i} -> y_{i+1}) / (q_{i} * Psp(y_{i} -> y_{i+1})) * a_{i-1}
if (weight.hasNonZero()) {
ray = Ray(hitPos, vout, ray.depth + 1);
pathContinue = true;
}
}
}
}
//------------------------------------------------------------------------
//------------------------------------------------------------------------
//if not terminated, calculate a more per-vertex information.
if (pathContinue) {
lightPath[nL].bsdfValue = r;
lightPath[nL].dirOPDF = dirPDF;
lightPath[nL].dirIPDF = mat->bsdf->pdf(vout, isect.gn, isect.sn, vin, sampledWLs);
lightPath[nL].rrOPDF = p;
lightPath[nL].rrIPDF = 0.5f;//fminf(1.0f, r.maxComp() * fabsf(Vector3f::Dot(vin, isect.sn)) / lightPath[nL].dirIPDF);
++nL;
}
else {
break;
}
//------------------------------------------------------------------------
}
//------------------------------------------------------------------------
//calculate the probability density with respect to surface area.
lightPath[0].dAOPDF = lightPDF * areaPDF;
lightPath[0].rrOPDF = 1.0f;
Vector3f adN = lightPath[2].pos - lightPath[1].pos;
float distN2 = adN.SqLength();
adN = adN.normalize();
lightPath[1].dAOPDF = lightPath[1].dirOPDF * fabsf(Vector3f::Dot(adN, lightPath[2].sn)) / distN2;
lightPath[1].rrOPDF = 1.0f;
for (int i = 2; i < nL; ++i) {
Vector3f adP = lightPath[i - 1].pos - lightPath[i].pos;
float distP2 = adP.SqLength();
adP = adP.normalize();
adN = lightPath[i + 1].pos - lightPath[i].pos;
distN2 = adN.SqLength();
adN = adN.normalize();
lightPath[i].dAIPDF = lightPath[i].dirIPDF * fabsf(Vector3f::Dot(adP, lightPath[i - 1].sn)) / distP2;
lightPath[i].dAOPDF = lightPath[i].dirOPDF * fabsf(Vector3f::Dot(adN, lightPath[i + 1].sn)) / distN2;
}
//------------------------------------------------------------------------
//END : generate light subpath.
//--------------------------------------------------------------------------------------------------------------
//--------------------------------------------------------------------------------------------------------------
//generate eye subpath.
isect = Intersection();
ray = Ray(camera->position, rayDir, 1);
weight = Color(1.0f, 1.0f, 1.0f);
dirPDF = 1.0f;
//eye - 0th vertex
weight = Color(1.0f, 1.0f, 1.0f);
eyePath[nE++] = ScatteringEvent(NULL, NULL, objID(),
vz, vz, vz,
weight,
-1, -1);
//eye - 1st vertex
weight = Color(1.0f, 1.0f, 1.0f);
eyePath[nE++] = ScatteringEvent(NULL, NULL, objID(),
ray.org, camDir, camDir,
weight,
dirPDF, -1);
eyePath[1].bsdfValue = Color(1.0f, 1.0f, 1.0f);
//trace eye subpath.
weight *= 1.0f;
collisionID = objID();
while (FindNearestIntersectionWithRay(&isect, ray, collisionID, sceneRef)) {
Vector3f hitPos = ray.org + ray.dir * isect.t;
Vector3f vout = -ray.dir;
Material* mat = isect.obj->materials[0];
collisionID = isect.objID;
eyePath[nE] = ScatteringEvent(isect.obj, mat->bsdf, collisionID,
hitPos, isect.gn, isect.sn,
weight,
-1, -1);
//------------------------------------------------------------------------
//reach the depth limit.
if (nE > LimitPathLength - 2) {
break;
}
//------------------------------------------------------------------------
bool pathContinue = false;
//------------------------------------------------------------------------
//sample transport direction.
//terminate path by russian roulette, except for the first few vertices.
Vector3f vin;
BSDFType sampledType;
BSDFSample bsdfSample = BSDFSample(op->_rand->GetRandom32bitReal0cTo1o(),
op->_rand->GetRandom32bitReal0cTo1o());
Color r = mat->bsdf->sample_f(vout, isect.gn, isect.sn, bsdfSample, &vin, &dirPDF, BSDF_ALL, &sampledType);
costh = fabsf(Vector3f::Dot(vin, isect.sn));
Color fraction = r * costh / dirPDF;
float p = 0.5f;//fminf(1.0f, fraction.maxComp());
if (costh != 0) {
if (ray.depth < numAbsLengthEye) {
weight *= fraction;
if (weight.hasNonZero()) {
ray = Ray(hitPos, vin, ray.depth + 1);
pathContinue = true;
}
}
else {
if (op->_rand->GetRandom32bitReal0cTo1o() < p) {
weight *= fraction / p;
if (weight.hasNonZero()) {
ray = Ray(hitPos, vin, ray.depth + 1);
pathContinue = true;
}
}
}
}
//------------------------------------------------------------------------
//------------------------------------------------------------------------
//if not terminated, calculate a more per-vertex information.
if (pathContinue) {
eyePath[nE].bsdfValue = r;
eyePath[nE].dirIPDF = dirPDF;
eyePath[nE].dirOPDF = mat->bsdf->pdf(vin, isect.gn, isect.sn, vout, sampledWLs);
eyePath[nE].rrIPDF = p;
eyePath[nE].rrOPDF = 0.5f;//fminf(1.0f, r.maxComp() * fabsf(Vector3f::Dot(vout, isect.sn)) / eyePath[nE].dirOPDF);
nE++;
}
else {
break;
}
//------------------------------------------------------------------------
}
//------------------------------------------------------------------------
//calculate the probability density with respect to surface area.
eyePath[0].dAIPDF = 1.0f;
eyePath[0].rrIPDF = 1.0f;
eyePath[1].dAIPDF = 1.0f;
eyePath[1].rrIPDF = 1.0f;
for (int i = 2; i < nE; ++i) {
Vector3f adP = eyePath[i - 1].pos - eyePath[i].pos;
float distP2 = adP.SqLength();
adP = adP.normalize();
Vector3f adN = eyePath[i + 1].pos - eyePath[i].pos;
float distN2 = adN.SqLength();
adN = adN.normalize();
eyePath[i].dAOPDF = eyePath[i].dirOPDF * fabsf(Vector3f::Dot(adP, eyePath[i - 1].sn)) / distP2;
eyePath[i].dAIPDF = eyePath[i].dirIPDF * fabsf(Vector3f::Dot(adN, eyePath[i + 1].sn)) / distN2;
}
//------------------------------------------------------------------------
//END : generate eye subpath.
//--------------------------------------------------------------------------------------------------------------
//--------------------------------------------------------------------------------------------------------------
//connect subpaths.
Vector3f connect;
Ray visRay;
float distance2;
double pPath;
double invWeight;
double pdfRatio;
Color uwCst;
//s = 0, eye subpath directly hit a light source.
for (int t = 2; t <= nE; ++t) {
if (eyePath[t].obj->materials[0]->edf != NULL) {
//for debugging, consider only 2 bounce light transport.
if (t != 4) {
continue;
}
//---------------------------------------------------------------------------------------------
//calculate MIS weight.
invWeight = 1.0f;
pPath = 1.0f;
float thisAreaPDF;
float thisDirPDF;
Ray dhRay = Ray(eyePath[t - 1].pos, (eyePath[t].pos - eyePath[t - 1].pos).normalize());//a ray which directly hits a light source.
eyePath[t].obj->pdf_radiance(dhRay, 1.0f, &thisAreaPDF, &thisDirPDF);
//t (- 2) -->
//not consider paths that a light subpath directly hits a lens surface and connect a vertex on a lens surface and a light subpath.
int ad = t + 1;
int rm;
for (int i = 0; i < t - 2; ++i) {
if (i >= nL) {
break;
}
rm = ad - 2;
if (i == 0) {
pdfRatio = (eyePath[t].obj->emissionPower() / sumLightPower * thisAreaPDF) / eyePath[rm].dAIPDF;
//consider russian roulette
if (rm > numAbsLengthEye) {
pdfRatio /= eyePath[rm].rrIPDF;
}
}
else if (i == 1) {
Vector3f temp = (eyePath[t].pos - eyePath[t - 1].pos);
pdfRatio = (thisDirPDF * fabsf(Vector3f::Dot(temp.normalize(), eyePath[t - 1].sn)) / temp.SqLength()) / eyePath[rm].dAIPDF;
//consider russian roulette
if (rm > numAbsLengthEye) {
pdfRatio /= eyePath[rm].rrIPDF;
}
}
else {
pdfRatio = eyePath[ad].dAOPDF / eyePath[rm].dAIPDF;
//consider russian roulette
if (t - ad >= numAbsLengthLight) {
pdfRatio *= eyePath[ad].rrOPDF;
}
if (rm > numAbsLengthEye) {
pdfRatio /= eyePath[rm].rrIPDF;
}
}
pPath *= pdfRatio;
invWeight += pPath * pPath;
--ad;
}
float Wst = 1.0f / invWeight;
//---------------------------------------------------------------------------------------------
//calculate unweighted contribution.
uwCst = eyePath[t].obj->radiance(dhRay, 1.0f, sampledWLs) * eyePath[t].weight;
op->_pixelValues[cIdx] += Wst * uwCst;
}
}
//s >= 1, t >= 2
int vIdx;
ScatteringEvent* uniDirPath[(LimitPathLength + 2) * 2];
for (int s = 1; s <= nL; ++s) {
for (int t = 2; t <= nE; ++t) {
//for debugging, consider only 2 bounce light transport.
if (s + t != 4) {
continue;
}
isect = Intersection();
connect = eyePath[t].pos - lightPath[s].pos;
distance2 = connect.SqLength();
connect = connect.normalize();
visRay = Ray(lightPath[s].pos, connect);
//if there are no objects that lie between the two vertices.
if (FindNearestIntersectionWithRay(&isect, visRay, lightPath[s].oid, sceneRef)) {
if (isect.objID.ID != eyePath[t].oid.ID || isect.objID.primitiveID != eyePath[t].oid.primitiveID) {
continue;
}
//---------------------------------------------------------------------------------------------
//calculate MIS weight. for convenience, two subpaths are incorporated into unidirectional path.
Vector3f lastLightSegment = lightPath[s - 1].pos - lightPath[s].pos;
float dist2LLS = lastLightSegment.SqLength();
lastLightSegment = lastLightSegment.normalize();
Vector3f lastEyeSegment = eyePath[t - 1].pos - eyePath[t].pos;
float dist2LES = lastEyeSegment.SqLength();
lastEyeSegment = lastEyeSegment.normalize();
vIdx = 0;
for (int i = 0; i < s; ++i) {
uniDirPath[vIdx++] = &(lightPath[i]);
}
ScatteringEvent lEnd = ScatteringEvent(lightPath[s].obj, lightPath[s].bsdf, lightPath[s].oid,
lightPath[s].pos, lightPath[s].gn, lightPath[s].sn,
lightPath[s].weight,
-1, -1);
ScatteringEvent eEnd = ScatteringEvent(eyePath[t].obj, eyePath[t].bsdf, eyePath[t].oid,
eyePath[t].pos, eyePath[t].gn, eyePath[t].sn,
eyePath[t].weight,
-1, -1);
uniDirPath[vIdx++] = &lEnd;
uniDirPath[vIdx++] = &eEnd;
float cosL = fabsf(Vector3f::Dot(connect, lEnd.sn));
float cosE = fabsf(Vector3f::Dot(connect, eEnd.sn));
Color lBSDF_O;
if (s > 1) {
lBSDF_O = lEnd.bsdf->f(lastLightSegment, lEnd.gn, lEnd.sn, connect, sampledWLs);
lEnd.dirIPDF = lEnd.bsdf->pdf(connect, lEnd.gn, lEnd.sn, lastLightSegment, sampledWLs);
lEnd.dirOPDF = lEnd.bsdf->pdf(lastLightSegment, lEnd.gn, lEnd.sn, connect, sampledWLs);
}
else {
lBSDF_O = lEnd.obj->radiance(Ray(lEnd.pos + connect, -connect), 1.0f, sampledWLs) / radiosity;
lEnd.dirIPDF = -1;
lEnd.dirOPDF = lEnd.obj->materials[0]->edf->pdf_radiance(lEnd.sn, connect);
}
if (lBSDF_O.hasNonZero() == false) {
continue;
}
lEnd.dAIPDF = lEnd.dirIPDF * fabsf(Vector3f::Dot(lastLightSegment, lightPath[s - 1].sn)) / dist2LLS;
lEnd.dAOPDF = lEnd.dirOPDF * cosE / distance2;
lEnd.rrIPDF = 0.5f;//lBSDF_O.maxComp() * fabsf(Vector3f::Dot(lastLightSegment, lEnd.sn)) / lEnd.dirIPDF;
lEnd.rrOPDF = 0.5f;//lBSDF_O.maxComp() * cosL / lEnd.dirOPDF;
Color eBSDF_O = eEnd.bsdf->f(-connect, eEnd.gn, eEnd.sn, lastEyeSegment, sampledWLs);
if (eBSDF_O.hasNonZero() == false) {
continue;
}
eEnd.dirIPDF = eEnd.bsdf->pdf(lastEyeSegment, eEnd.gn, eEnd.sn, -connect, sampledWLs);
eEnd.dirOPDF = eEnd.bsdf->pdf(-connect, eEnd.gn, eEnd.sn, lastEyeSegment, sampledWLs);
eEnd.dAIPDF = eEnd.dirIPDF * cosL / distance2;
eEnd.dAOPDF = eEnd.dirOPDF * fabsf(Vector3f::Dot(lastEyeSegment, eyePath[t - 1].sn)) / dist2LES;
eEnd.rrIPDF = 0.5f;//eBSDF_O.maxComp() * cosE / eEnd.dirIPDF;
eEnd.rrOPDF = 0.5f;//eBSDF_O.maxComp() * fabsf(Vector3f::Dot(lastEyeSegment, eEnd.sn)) / eEnd.dirOPDF;
for (int i = t - 1; i > -1; --i) {
uniDirPath[vIdx++] = &(eyePath[i]);
}
invWeight = 1.0f;
pPath = 1.0f;
int ad;
int rm;
//calculate p_{s+1} ~ p_{s+t-2}.
for (int i = s; i < s + t - 2; ++i) {
if (i >= nL) {
break;
}
ad = i;
rm = ad + 2;
pdfRatio = uniDirPath[ad]->dAOPDF / uniDirPath[rm]->dAIPDF;
//consider russian roulette
if (ad > numAbsLengthLight) {
pdfRatio *= uniDirPath[ad]->rrOPDF;
}
if (s + t - rm >= numAbsLengthEye) {
pdfRatio /= uniDirPath[rm]->rrIPDF;
}
pPath *= pdfRatio;
invWeight += pPath * pPath;
}
//calculate p_{s-1} ~ p{0}.
pPath = 1.0f;
for (int i = t; i < s + t; ++i) {
if (i >= nE) {
break;
}
ad = s + t - i + 1;
rm = ad - 2;
pdfRatio = uniDirPath[ad]->dAIPDF / uniDirPath[rm]->dAOPDF;
//consider russian roulette
if (rm > numAbsLengthLight) {
pdfRatio /= uniDirPath[rm]->rrOPDF;
}
if (s + t - ad >= numAbsLengthEye) {
pdfRatio *= uniDirPath[ad]->rrIPDF;
}
pPath *= pdfRatio;
invWeight += pPath * pPath;
}
float Wst = 1.0f / invWeight;
//END : calculate MIS weight.
//---------------------------------------------------------------------------------------------
//calculate unweighted contribution.
uwCst = lightPath[s].weight * lBSDF_O * cosL * cosE / distance2 * eBSDF_O * eyePath[t].weight;
op->_pixelValues[cIdx] += Wst * uwCst;
}
}
}
//END : connect subpaths.
//--------------------------------------------------------------------------------------------------------------
}
}