Problems in custom ray marcher

[Efoia]

Hello all, I've been struggling with a few problems in my custom ray marcher for quite some time (on and off for months) and was hoping to get some guidance.  The ray marcher works fine for the KIFS fractals.  For other fractals however, I get strange clipping and pinching issues.  I've stripped my code down to the bare minimum needed to reproduce the problem and have included it in this post.  I've tried distance equations from multiple sources and the problems persist, so this leads me to believe that the source of the problems is in the ray marcher itself.

Below are some of the issues I'm experiencing, the image colors correspond the the surface normals.

The julia fractal appears correct at close camera distances, but starts clipping as the camera moves away.  The screenshots below from top to bottom are me moving the camera away from the fractal:









Here is the same problem occurring with the mandelbulb:









The knot fractal appears to work correctly initially:



But at some camera angles I see what looks like pinches in the surface:



I also get these artifacts:



Here is a variant that looks correct from one angle:



But from a side camera angle it appears completely incorrect:



Here is my ray marcher stripped down to only the necessary parts:

Code:

#version 330 core

uniform mat4 M;
uniform vec2 resolution;
uniform float epsilon;

layout(location = 0) out vec3 color;

float DE(vec3 p);
vec3 NDE(vec3 p);

vec3 MultiplyMatrixPoint(vec3 point, mat4 matrix) { return (matrix * vec4(point, 1)).xyz; }
vec3 MultiplyMatrixVector(vec3 vector, mat4 matrix) { return (matrix * vec4(vector, 0)).xyz; }
vec3 MultiplyModelMatrixPoint(vec3 vector) { return (M * vec4(vector, 1)).xyz; }
vec3 MultiplyModelMatrixVector(vec3 vector) { return (M * vec4(vector, 0)).xyz; }

float DE(vec3 pos)
{
    // Implementation
    // ...
}

vec3 NDE(vec3 point)
{
    return normalize(vec3(DE(point + vec3(epsilon, 0, 0)) - DE(point - vec3(epsilon, 0, 0)),
                          DE(point + vec3(0, epsilon, 0)) - DE(point - vec3(0, epsilon, 0)),
                          DE(point + vec3(0, 0, epsilon)) - DE(point - vec3(0, 0, epsilon))));
}

vec3 ray(vec3 ro, vec3 rd)
{
    float t = 0.0;

    for (int i = 0; i < 400; ++i)
    {
        vec3 p = ro + rd * t;

        float d = DE(p);
        t += d;

        if (d < epsilon)
        {
            // Return the normal as the color
            return NDE(ro + rd * t);
        }
    }

    return vec3(0);
}

void main()
{
    float focalLength = .1;
    float zoomFactor = 0.07;

    vec2 pos;

    pos.x = gl_FragCoord.x / resolution.x;
    pos.y = gl_FragCoord.y / resolution.y;

    pos.x = zoomFactor * pos.x - (zoomFactor * .5);
    pos.y = zoomFactor * pos.y - (zoomFactor * .5);

    pos.x *= resolution.x / resolution.y;

    vec3 rd = normalize(vec3(pos.x, pos.y, -focalLength));
    vec3 ro = vec3(pos.x, pos.y, 0.0);

    ro = MultiplyModelMatrixPoint(ro);
    rd = MultiplyModelMatrixVector(rd);

    color = ray(ro, rd);
}

Here is how I'm calculating the model matrix that I use for the camera view:

Code:

glm::mat4 modelMatrix = glm::mat4(1.0f);
    
modelMatrix = glm::translate(modelMatrix, glm::vec3(-pan_x, pan_y, 0.0f));
modelMatrix = glm::rotate(modelMatrix, static_cast<float>(-rotate_x), glm::vec3(0, 1, 0));
modelMatrix = glm::rotate(modelMatrix, static_cast<float>(-rotate_y), glm::vec3(1, 0, 0));
modelMatrix = glm::translate(modelMatrix, glm::vec3(0.0f, 0.0f, -pan_z));

Julia distance equation for reference.  I don't think the distance equations themselves are relevant to the problem but I can post the others if requested:

Code:

float DE(vec3 point)
{
    vec4 c = vec4(0.0, 0.0, 0.0, 0.0);

    vec4 z = vec4(point, 0.0);
    float md2 = 1.0;
    float mz2 = dot(z, z);
    vec4 nz;

    for (int i = 0; i < 11; i++)
    {
        md2 *= 4.0*mz2;
        nz.x = z.x*z.x - dot(z.yzw, z.yzw);
        nz.yzw = 2.0*z.x*z.yzw;
        z = nz + c;

        mz2 = dot(z, z);
        if (mz2 > 10)
        {
            break;
        }
    }

    return 0.25*sqrt(mz2 / md2)*log(mz2);
}

Thanks for your time.

[Buddhi]

It looks like overstepping during ray-marching.
Try to reduce ray-marching step by multiplying estimated distance by some number lower than 1.0, e.g.:

        float d = DE(p);
        t += 0.5 * d;

DE() calculates only estimated distance so it can be much to high in some cases.

[Efoia]

Thank you so much for your reply, overstepping was precisely the source of the problem!  In a single line of code you have ended a months-long headache for me!   

From the complete renderer: