I agree, it is difficult. And I had read the two papers during literature survey but there always seems to be a gap in explaining how they get from the equation as stated using forward and backward differences and this code implementation. However, we know what we really need to compute. Given the distance/time values of 2 adjacent neighbours, how do we compute the distance/time it takes to reach the pixel in question ? (Distance and time are interchangeable since the speed is taken to be 1 in a direction normal to the curve at all points.)

Let's break it down a bit:

In `_fast_marching_method` I always pass the adjacent neighbours, in the sense:

    0 A 0
    B ? D
    0 C 0

I'll pass A-B or A-D or B-C or C-D not A-C or B-D. So that's where the `2` comes into the picture in `r` computation. Essentially, since the distance between them is `sqrt(2)`, that is the max difference between there distance/time values. 

So, in what case will it be lesser ? I can imagine a curve inclined at -45 deg hitting B and C. In this case, intuitively, the values should be the same. Assume the values for both is `1`. So my `r` would simply be `sqrt(2)`. Now, since my `s` value according to L210 would be less than `1`. I go to L214. And `s = (u1 + u2 + r) / 2`. This would be essentially be `1 / sqrt(2)` times greater than `u1` or `u2` i.e. `1`.

Now if we imagine a isosceles right angle triangle as:

    B -- ?
     \   |
       \ |
         C
And the curve parallel to side BC, then yes the distance value given as `1 + 1/sqrt(2)` for `?` given `1` for `B` and `C`, does make sense.

Does this train of thought make sense ?