Resulting Probabilty Density in Path Tracer for paths using Next Event Estimation I am trying to implement my own Gradient Domain Path Tracer by following the code of this guy who already implemented it: 

https://gist.github.com/BachiLi/4f5c6e5a4fef5773dab1

I already managed to go through various steps but I wanted to do something more. I expanded the code of in the reference by implementing next event estimation and here are some results.

Normal Path Tracer image:



Gradient Domain resulting image:



Results are already good.. but as said before, I wanted something more. So I implemented Next Event Estimation and here is the result of the basic Path Tracer:



Here is my code:

private Vector3 SampleWithNEE( Ray ray )
{
  // prepare
  Vector3 T = (1,1,1), E = (0,0,0), NL = (0,-1,0);
  int depth = 0;
  // random walk
  while (depth++ &lt; MAXDEPTH)
  {
    // find nearest ray/scene intersection
    Scene.Intersect( ray );
    if (ray.objIdx == -1) break; //if there is no intersection
    Vector3 I = ray.O + ray.t * ray.D; //go to the Hit Point on the scene
    Material material = scene.GetMaterial( ray.objIdx, I );
    if (material.emissive) //case of a light
    {
        E += material.diffuse;
        break;
    }
    // next event estimation
    Vector3 BRDF = material.diffuse * 1 / PI;
    float f = RTTools.RandomFloat();
    Vector3 L = Scene.RandomPointOnLight() - I;
    float dist = L.Length();
    L = Vector3.Normalize( L );
    float NLdotL = Math.Abs( Vector3.Dot( NL, -L ) );
    float NdotL = Vector3.Dot( ray.N, L );
    if (NdotL &gt; 0)
    {
        Ray r = new Ray( I + L * EPSILON, L, dist - 2 * EPSILON ); //make it a tiny bit shorter otherwise I risk to hit my starting and destination point
        Scene.Intersect( r );
        if (r.objIdx == -1) //no occlusion towards the light
        {
            float solidAngle= (nldotl * light.getArea()) / (dist * dist);
            E += T * (NdotL) * solidAngle * BRDF * light.emission;
        }
    }
    // sample random direction on hemisphere
    Vector3 R = DiffuseReflectionCosWeighted( ray.N );
    float hemi_PDF = Vector3.Dot( R, ray.N ) / PI;
    T *= (Vector3.Dot( R, ray.N ) / hemiPDF) * BRDF;
    ray = new Ray( I + R * EPSILON, R, 1e34f );
  }
  return E;
}


Things work really fine and the results are shown with the picture above. One more thing: I have only diffuse surfaces in my scene. 

Now, the problem is that in this method I use 2 kind of PDF's:


one is given by randomly sampling the light in the Direct Lighting of the Next Event Estimation, indeed that SolidAngle is our PDF, or better 1/PDF. 
while the second PDF is the one led by the use of DiffuseReflectionCosWeighted which brings a PDF equals to CosTheta/PI. 


Everything fine so far and for any implementation detail you can just look at my code but problems come with my Gradient Domain Path Tracer. Indeed, there, as also in the reference link above implemented by Tzu-Mao Li, I need the final probability density of a whole path to compute the final gradient image. How did I calculate it in case without Next Event Estimation (NEE) ? In that case (since I have only Diffuse surfaces) this probability is the product of CosTheta / PI at each bounce in the scene. Everything is fine and the resulting gradient image is shown above.

Instead, in case I use NEE things don't work anymore because the Probability Density of my whole path changes and I don't manage to understand how it is. The resulting Gradient Domain Image with Next Event Estimation is:



I need to understand how to calculate the final density probability of a path. Can you help me doing it? Thanks in advance!
I don't have any experience with Gradient Domain Path Tracing, but here are my thoughts:

There seems to be a different problem

If you look carefully at the little spikes of distortion in the final image, you will see that they are all lit from the same direction - on their top left side at a consistent 45 degrees. The sphere also appears to be lit from this angle, rather than from above by the light source.

This is unlikely to be explained by an incorrect probability estimation for a path. I would expect there to be a different problem with the code, that these distortions are hinting at.

I will therefore address these two separate points:


You want to know how to calculate the Probability Density of a path when using Next Event Estimation.
There is evidence of some problem unrelated to this.


I'll also review the code for non-essential points - but I'll leave this until after the essentials.

Probability Density of a path when using Next Event Estimation

Looking at the paper on which the code you are following was based it seems that the novel shift mapping described in section 5.2 is defined in terms of the reflective properties of the surfaces found at the vertices of the path. I must emphasise that I don't have a full understanding of this, but it suggests that Next Event Estimation may not require a change in this approach, as the surfaces encountered will be the same. Hopefully once the other problems are cleared up it will be easier to judge whether the image looks correct.

Note that section 5.2 of the paper already mentions (just below Figure 10) that they take into account sampling the emitter "either using BSDF or area sampling".

The difference with Next Event Estimation is that the area sampling happens at every vertex of the path, but it isn't obvious to me that this should cause a problem.

The fact that your scene uses only diffuse surfaces means that the offset path should in most cases rejoin the base path at the second vertex, so you would only need to recalculate the area sampling for the first vertex of the offset path.

The cause of the incorrect lighting direction

In reading through the code to familiarise myself with how it works, I noticed that NLdotL is calculated, but then not used. A text search revealed that the only other occurrence has a different case: nldotl. Here are the two variables in context (first and ninth lines of this excerpt):

float NLdotL = Math.Abs( Vector3.Dot( NL, -L ) );
float NdotL = Vector3.Dot( ray.N, L );
if (NdotL &gt; 0)
{
    Ray r = new Ray( I + L * EPSILON, L, dist - 2 * EPSILON ); //make it a tiny bit shorter otherwise I risk to hit my starting and destination point
    Scene.Intersect( r );
    if (r.objIdx == -1) //no occlusion towards the light
    {
        float solidAngle= (nldotl * light.getArea()) / (dist * dist);
        E += T * (NdotL) * solidAngle * BRDF * light.emission;
    }
}


Since nldotl is not defined, the result of the code is undefined behaviour. In practice, the program is likely to act as if nldotl is zero, or for some compilers perhaps a constant arbitrary value, or even a different arbitrary value at each iteration. For your particular compiler it appears to be a constant value, and I strongly suspect this is the cause of the distinctive alignment of lighting angle on all the speckles and the sphere. If there is also another contributing problem it will be easier to analyse it in a separate question once this initial problem has been fixed.

It may be worth considering using a compiler and/or flag setting that gives an error or at the very least a warning for undefined variables, as this type of mistake is both very easy to make, and very easy to overlook later.

Additional contribution of light source

There appears to be another problem that will cause results to be incorrect in a more subtle way, with no obvious distortion. Due to Next Event Estimation, the light source is contributing to each step along the path. This means it should not make a contribution if the path itself hits the light source directly. Otherwise, for that path the light source will be contributing twice. You can correct this by changing:

if (material.emissive) //case of a light
{
    E += material.diffuse;
    break;
}


to:

if (material.emissive) //case of a light
{
    break;
}


This will give a zero contribution from the intersection with the light.

Note that since I can only see this one function, I can't guess whether this will also make the light appear black in the image. You may or may not need to adjust for rays that start at the camera and hit the light directly.



Code review

Doubly finite rays

I'm used to defining a ray as a half infinite line segment - having a starting point but no ending point. I notice that this code gives a ray both a starting point and a length. The only place I can see a reason for this is when testing a shadow ray against the light source: the code checks that there are no intersections on the way to the light, so intersections behind the light (or on the light itself) must be excluded. In all other places, the ray is defined with a pseudo-infinite length (1e34f).

The following suggestion won't affect the correctness of your code, but it may be more readable and avoid you having to work around the need for infinity and having to account for epsilon twice.

If the ray is simply a starting point and a direction, then shadow rays can simply check that the first intersection is the light, rather than checking that there is no intersection. For example, by replacing:

if (r.objIdx == -1) //no occlusion towards the light


with:

if (r.objIdx == LIGHT) // There is no object intersected before the light


Here I use LIGHT as a placeholder for the id of the light source, since that part of the code is not included in the question.


This will always give false if the light is occluded by a nearer object.
This will always give true if the ray hits the light before any more distant objects.


This is therefore equivalent to the current code, but does not require the ray to store a length.

Lights that don't reflect

The code currently models lights and surfaces separately. This means that if an object is a light then it is only emissive, and reflects no light from other objects.

This causes a negligible difference in the example scene in the question, which has a single bright light. However, if used with a number of dimmer lights, it would be more noticeable that they do not light each other. In many cases the difference will not be noticeable, so this is not a problem if it works for the scenes you wish to render.