MIS BPT with Russian Roulette

[shocker_0x15]

Hello.

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.

refPT_1000P.png (241.38 KiB) Viewed 7687 times

BPT_withoutRR_0300P.png (261.71 KiB) Viewed 7687 times

BPT_withRR_0300P.png (319.15 KiB) Viewed 7687 times

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.
        //--------------------------------------------------------------------------------------------------------------
    }
}

[ingenious]

How do you account for RR in the MIS weight? You just need to multiply the pdf of the bounce direction with the survival probability. As I see from the code, you have tested with a constant probability of 0.5.

[shocker_0x15]

Thank you for your prompt reply.
As you say, now I'm testing RR with a constant probability 0.5.
The first few bounces don't consider RR, absolutely continue (but in the above code, RR is considered from the first bounce).

I took RR into account in MIS weight by the below code.
The precomputed PDFs (for surface area) don't contain RR probability, so a path generation probability is multiplied by or divided by RR probability.

Code: Select all

//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;   
}

I assume that it computes correctly, but the result seems incorrect.

Figure_03.png (25.26 KiB) Viewed 7614 times

[spectral]

Hi,

You tell that you don't account for MIS at x(0) (camera vertex), so maybe you have a very small geometric term that does not cancel due to RR !
I say this because fireflies appears in the corner, when the distance between x(i) and y(i) is near 0 !

What you can try is to put check this when you do the connection, if G is too small then use 0.1. Just for test of course ! Something like :

Code: Select all

distance2 = max(0.01, distance2)

Honestly, I'm not sure that it is a good practice to put a RR for small path, maybe use RR only for path length > 3 ! For sure, Ingenious maybe you can confirm ?
3

[shocker_0x15]

Thank you spectral.

As you say, I tested getting rid of very small G factor by your code.
As we expected, fireflies seem to be decreased, but whole of the image of BPT using RR seems to be brighter as before.

Now I'm testing with minimum absolute path length (it is 1) for convenience.
In practice, of course I set it about 3.

decFireflies.png (325.96 KiB) Viewed 7614 times

[apaffy]

It's difficult to find a bug in the code just by looking, but I thought I'd comment on this:

Of course, I know that it is required that the probability for generating a adjacent vertex in MIS weight computations contain RR continuation probability.

Actually, as Dietger noted on another thread, you don't need to take RR into account for MIS weights, any weights are valid as long they sum to 1. You might lose a very small amount of efficiency, but it's completely valid to construct MIS weights ignoring RR, so perhaps this will help to find the bug.

(Edit: in this post viewtopic.php?f=3&t=516&p=1533#p1533)

[shocker_0x15]

My first post contained misleading sentence.

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.

to be exact

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 (strictly, minimum subpath lengths are set 3), the third is the BPT which uses RR.
The first twos look the same, but the third is not, looks slightly brighter.

In "BPT doesn't use RR", minimum subpath lengths are set 3.
So for 2 bounce light transport, RR must not be used.

----------------

Thanks apaffy.

I simply tried to comment out multiplying and dividing by RR in MIS computations. (not in subpath traceing)
But the result yet seems incorrect:?

LuxRender includes BPT implementation and it seems to consider RR in MIS weight computations.

[apaffy]

I simply tried to comment out multiplying and dividing by RR in MIS computations. (not in subpath traceing)
But the result yet seems incorrect:?

LuxRender includes BPT implementation and it seems to consider RR in MIS weight computations.

Then I think you have to conclude that the remaining MIS code is not correct. If you use uniform weights is the image correct? If you accumulate just the weights does it converge to an image of 1s?

[shocker_0x15]

Thank you apaffy.

I tried to check uniform weighting (in this case, there are 3 strategies [s=0, t=4], [s=1, t=3], [s=2, t=2], they are assigned 0.33333... weight).
And I got the result, which shows that BPT and reference will converge to the same.
The minimum subpath lengths are set 1 there.

compUniBPT.png (297.93 KiB) Viewed 7572 times

From this result, I think that my MIS weights computation part is incorrect.

If you accumulate just the weights does it converge to an image of 1s?

MIS weights for certain length sum to 1 only if each strategy exactly traces the same path, don't it ?
So, do you mean that the sum of MIS weights "converges" to an image of 1s if it is averaged by the number of passes ?

I tried to create two weight maps.
The one represents the result that it doesn't use RR.
The other represent the result that it uses RR.
I assigned 100 brightness to 1.0 weights.

weights_dist_withoutRR.png (160.15 KiB) Viewed 7572 times

weights_dist_withRR.png (155.46 KiB) Viewed 7572 times

Each of them seems obviously different from a image of 1s.
Is it caused by subpath pairs that cannot connect each other ?

[apaffy]

I tried to check uniform weighting (in this case, there are 3 strategies [s=0, t=4], [s=1, t=3], [s=2, t=2], they are assigned 0.33333... weight).
And I got the result, which shows that BPT and reference will converge to the same.

From this result, I think that my MIS weights computation part is incorrect.

Yep if the weights produce different results to the uniform case then there must be some error in the weights.

If you accumulate just the weights does it converge to an image of 1s?

MIS weights for certain length sum to 1 only if each strategy exactly traces the same path, don't it ?
So, do you mean that the sum of MIS weights "converges" to an image of 1s if it is averaged by the number of passes ?

Hmm I should have sanity checked this debug suggestion before posting, I remember doing something similar before but plainly this not a useful image as-is. Please ignore this until I regain my sanity, sorry!