Welcome to Fractal Forums

float mandelBoxDistance(float4 v,float scale,int maxIter)
{
	int iter = 0;
	v = fabs(v);
	float4 c = v;
	float distance = scale;
	float fixedRadius = 1.0f;

	float fR2 = fixedRadius * fixedRadius;
	float minRadius = 0.5f;
	float mR2 = minRadius * minRadius;

	while (iter < maxIter)
	{
		if (v.x > 1) v.x = 2-v.x;
		else if (v.x < -1) -2-v.x;
		if (v.y > 1) v.y = 2-v.y;
		else if (v.y < -1) -2-v.y;
		if (v.z > 1) v.z = 2-v.z;
		else if (v.z < -1) -2-v.z;		
		
		float r2 = v.x*v.x+v.y*v.y+v.z*v.z;
		if (r2 > 150)
			break;
		if (r2 < mR2)
		{	
			v *= fR2 / mR2;
			distance *= fR2 / mR2;
		}
		else if (r2 < fR2) 
		{
			v *= fR2 / r2;
			distance *= fR2 / r2;
		}
		distance *= scale;
		v = v * scale + c;
		iter++;
	}
	return length(v)/fabs(distance);
}
kernel void MandelBox(const float4 camera,
					  const float4 middle,
					  const float4 u,
					  const float4 v,
					  const int maxIter,
					  const int2 size,
					  const float fieldOfView,
					  const float farPlane,
					  const float scale,
					  __global write_only int* output,
					  __global write_only int* zBuffer)
{
	int x = get_global_id(0);
	int y = get_global_id(1);

	int index = 3*(x + y * size.x);

	float relx = (-0.5f+(float)x/(float)size.x);
	float rely = (-0.5f+(float)y/(float)size.y);

	float4 g = middle + relx * u * fieldOfView + rely * v * fieldOfView;

	float4 direction = normalize(g - camera);

	float t = 0.0f;
	int steps = 0;
	while (t < farPlane)
	{
		steps++;
		float4 c = camera + t * direction;
		float distance = 0.5f*mandelBoxDistance(c,scale,maxIter);
		if (distance < 0.001f)
			break;
		t += distance;
	}
	float4 pos = camera + t * direction;
	float deltaD = 0.001f;
	float4 deltaX = (float4)(deltaD,0,0,0);
	float4 deltaY = (float4)(0,deltaD,0,0);
	float4 deltaZ = (float4)(0,0,deltaD,0);
	float4 normal;
	float4 lightSource = camera - (float4)(0,0,0,0);
	normal.x = mandelBoxDistance(pos+deltaX,scale,maxIter) - mandelBoxDistance(pos-deltaX,scale,maxIter);
	normal.y = mandelBoxDistance(pos+deltaY,scale,maxIter) - mandelBoxDistance(pos-deltaY,scale,maxIter);
	normal.z = mandelBoxDistance(pos+deltaZ,scale,maxIter) - mandelBoxDistance(pos-deltaZ,scale,maxIter);
	normal.w = 0;
	normal = normalize(normal);

	float4 lightVector = normalize(lightSource - pos);
	float4 halfVector = normalize(direction + lightVector);
	int specularPower = 5;
	float color = fmax(0.0f,dot(normal,lightVector));
	// pown(fmax(0.0f,dot(normal,halfVector)),specularPower)
        // Blinn Phong shading looks rather ugly atm
	output[index] = (uchar)(255 * color);
	output[index+1] = (uchar)(255 * color);
	output[index+2] = (uchar)(255 * color);
	if (t >= farPlane)
	{
		output[index] = 0;
		output[index+1] = 0;
		output[index+2] = 0;
	}
    zBuffer[index/3] = (int)((t / farPlane) * 256 * 256 * 256 * 128);
}

__kernel void computeMandelbulb(
    __global float *DE,
    __global float *cx,
    __global float *cy,
    __global float *cz,
    int maxStep,
    int maxIter,
    __global int *cstep,
    const float p
    )
{

    const int gid = get_global_id(0);

    float x = cx[gid];
    float y = cy[gid];
    float z = -4.0;

    float nx = x;
    float ny = y;
    float nz = z;

    int step = 0;
    int iter = 0;
    float r = 0.0;
    float dr = 1.0;
    r=sqrt(x*x+y*y+z*z);
    float th = atan(y/x)*p;
    float ph = asin(z/r)*p;
    float r2p = pow(r,p);

    for(step = 0; step < maxStep && z < 4.0; step++) {
        //r2p = pow(r,p);
        iter = 0;
        r = 0.0;
        dr = 1.0;
        nx = x;
        ny = y;
        nz = z;
        r=sqrt(x*x+y*y+z*z);
        //th = atan(y/x)*p;
        //ph = asin(z/r)*p;
        while(iter < maxIter && r<2.0  ) {
            r2p = pow(r,p);
            th=atan(ny/nx)*p;
            ph=asin(nz/r)*p;
            nx=r2p*cos(ph)*cos(th)+x;
            ny=r2p*cos(ph)*sin(th)+y;
            nz=r2p*sin(ph)+z;
            dr=dr*p*pow(r,p-1.0)+1.0;
            r=sqrt(nx*nx+ny*ny+nz*nz);
            iter++; 
        }

        if(0.5*log(r)*r/dr < 0.00001) {
            break;
        } else {
            z+=0.5*log(r)*r/dr;
        }
    }

    if(0.5*log(r)*r/dr < 0.00001 && 0.5*log(r)*r/dr > -0.000001) {
//      DE[gid] = r;//0.5*log(r)*r/dr;
//      cx[gid] = x;
//      cy[gid] = y;
        cz[gid] = z;
        cstep[gid] = step;
    } else {
        cstep[gid] = 0;

    }

}