Welcome to Fractal Forums

    static float eye = 0.0f; eye = glm::clamp(eye+0.005f,0.0f,1.0f); // animate
    const glm::mat4 projection = glm::perspective(35.0f, (float)width/height, 0.1f, 10.0f),
            modelview = glm::lookAt(glm::vec3(cos(eye),sin(eye),-1),glm::vec3(0,0,0),glm::vec3(0,0,1));     
    const float epsilon = sqrt(1.0f/std::max(width,height))/2.0f;
    for(int y=0; y<height; y++) {
            for(int x=0; x<width; x++) {
                    glm::vec3 p = glm::unProject(glm::vec3(x,y,0),modelview,projection,glm::vec4(0,0,width,height)),
                            dir = glm::unProject(glm::vec3(x,y,1),modelview,projection,glm::vec4(0,0,width,height))-p;
                    //std::cout << x << "," << y << " " << p.x << "," << p.y << "," << p.z << " " << dir.x << "," << dir.y << "," << dir.z << std::endl;
                    float D = 0;
                    for(int i=0; i<MAX_ITER; i++) {
                            const float d = DE(dir*D + p);
                            D += d;
                            if(d<epsilon) {
                                    depth_bmp[y*width+x] = 255.0f/(i+1);
                                    break;
                            }
                    }
            }
    }

    float DE(glm::vec3 p) {
    	const float Scale = -1.77f, fixedRadius2 = 1.0f, minRadius2 = (0.5f*0.5f);
    	const glm::vec3 p0 = p;
    	float dr = 1.0f;
    	for(int n = 0; n < 13; n++) {
    		// Reflect
    		p = (glm::clamp(p,-1.0f,1.0f) * 2.0f) - p;
    		// Sphere Inversion
    		const float r2 = glm::dot(p,p);
    		if(r2<minRadius2) {
    			const float t = (fixedRadius2/minRadius2);
    			p *= t;
    			dr *= t;
    		} else if(r2<fixedRadius2) {
    			const float t = (fixedRadius2/r2);
    			p *= t;
    			dr *= t;
    		}
    		// Scale & Translate
                    p = p * Scale + p0;
                    dr = dr * abs(Scale) + 1.0f;
    	}
    	return glm::length(p)/abs(dr);
    }

enum { DEBUG = 0 };
const float rotation = (double)now_nanosecs()/1000000000, distance = 8.0f;
const glm::vec3 eye = glm::vec3(cos(rotation),0.5f,sin(rotation))*distance;
const glm::mat4
	projection = glm::perspective(35.0f,(float)depth_w/depth_h,0.1f,100.0f),
	modelview = glm::lookAt(eye,glm::vec3(0,0,0),glm::vec3(0.0f,1.0f,0.0f));	
for(int y=0; y<depth_h; y++) {
	for(int x=0; x<depth_w; x++) {
		glm::vec3
			p = glm::unProject(glm::vec3(x,y,0),modelview,projection,glm::vec4(0,0,depth_w,depth_h)),
			dir = glm::unProject(glm::vec3(x,y,1),modelview,projection,glm::vec4(0,0,depth_w,depth_h))-p;
		const float len = glm::length(dir),
			epsilon = 0.001; // how best to compute scaled epsilon too?
		dir = glm::normalize(dir);
		if(DEBUG) std::cout << x << "," << y << " " << p.x << "," << p.y << "," << p.z << " " << dir.x << "," << dir.y << "," << dir.z << std::endl;
		float D = 0.01f;
		for(int i=0; i<MAX_ITER; i++) {
			const float d = DE(p + dir*D);
			if(DEBUG) std::cout << " " << x << "," << y << " " << i << " " << D << "," << d << "," << (D+d) << "," << (D+d>len) << std::endl;
			D += d;
			if(D>len)
				break;
			if(d<epsilon) {
				depth_bmp[y*depth_w+x] = 255.0f/(i+1); // colour by ray march steps
				break;
			}
		}
	}
}

	const float Scale = -1.77f, fixedRadius2 = 1.0f, minRadius2 = (0.5f*0.5f);
	const glm::vec3 p0 = p;
	float dr = Scale;
	for(int n = 0; n < 13; n++) {
		// Reflect
		p = (glm::clamp(p,-1.0f,1.0f) * 2.0f) - p;
		// Sphere Inversion
		const float r2 = glm::dot(p,p);
		if(r2<minRadius2) {
			const float t = (fixedRadius2/minRadius2);
			p *= t;
			dr *= t;
		} else if(r2<fixedRadius2) {
			const float t = (fixedRadius2/r2);
			p *= t;
			dr *= t;
		}
		// Scale & Translate
                p = p * Scale + p0;
                dr *= Scale;
	}
	return glm::length(p)/abs(dr);

float AbsScalem1 = abs(Scale - 1.0);
float AbsScaleRaisedTo1mIters = pow(abs(Scale), float(1-Iterations));

float DE(vec3 pos) {
	vec4 p = vec4(pos,1), p0 = p;  // p.w is the distance estimate
	for (int i=0; i<Iterations; i++) {
		p.xyz = clamp(p.xyz, -1.0, 1.0) * 2.0 - p.xyz;  
		float r2 = dot(p.xyz, p.xyz);
		p *= clamp(max(MinRad2/r2, MinRad2), 0.0, 1.0);  
		p = p*scale + p0;
	}
	return ((length(p.xyz) - AbsScalem1) / p.w - AbsScaleRaisedTo1mIters);
}

	const int Iterations = 13;
	const float Scale = 2.0f,
		MinRad2 = (0.5f*0.5f),
		AbsScalem1 = abs(Scale - 1.0f),
		AbsScaleRaisedTo1mIters = pow(abs(Scale),float(1-Iterations));
	const vec3 p0 = p;
	float p_w = 1.0f;
	for(int i=0; i<Iterations; i++) {
		p = (clamp(p,-1.0f,1.0f) * 2.0f) - p;
		const float r2 = dot(p,p),
			adj = clamp(max(MinRad2/r2,MinRad2),0.0f,1.0f)*Scale;
		p = p*adj + p0;
		p_w = p_w*adj + 1.0f;
	}
	return ((length(p)-AbsScalem1)/p_w) - AbsScaleRaisedTo1mIters;

	const Vec3 pos = p;
	const int Iterations = 20, Power = 8;
	const Float DEPTH_OF_FIELD = 1000.0;
	Float dr = 1.0f, r = 0.0f;
	for(int i=0; i<Iterations; i++) {
		r = length(p);
		if(r>DEPTH_OF_FIELD) break;
		const Float theta = acos(p.z/r)*Power,
        		phi = atan2(p.y,p.x)*Power;
		dr =  pow(r,Power-1)*Power*dr + 1.0f;
		p.x = sin(theta)*cos(phi);
		p.y = sin(phi)*sin(theta);
		p.z = cos(theta);
		p *= pow(r,Power);
		p += pos;
	}
	return 0.5f*log(r)*r/dr;