double mandelbulb

#pragma kernel CSMain
struct Ray
{
    float3 origin;
    float3 direction;
};

RWTexture2D<float4> Result;
Texture2D<float4> Source;
float maxDst;
float minStepSize;
StructuredBuffer<double> power;
StructuredBuffer<double> spread;
StructuredBuffer<double> scale;
int maxIterations;
float3 lightDir;

float4x4 _CameraToWorld;
float4x4 _CameraInverseProjection;

double length(double3 f)
{
    return sqrt(f.x*f.x + f.y*f.y + f.z*f.z);
}
double dot(double3 f1, double3 f2)
{
    return f1.x * f2.x + f1.y * f2.y + f1.z * f2.z;
}

// Mandelbulb distance estimation:
// http://blog.hvidtfeldts.net/index.php/2011/09/distance-estimated-3d-fractals-v-the-mandelbulb-different-de-approximations/
double2 DstAndIterations(double3 position) {
    position *= scale[0];
    double3 z = position * spread[0];
	double dr = 1.0;
	double r = 0.0;
    int iterations = 0;
    double powe = power[0];

	for (int i = 0; i < maxIterations ; i++) {
        iterations = i;
		r = length(z);

		if (r>2) {
            break;
        }

		// convert to polar coordinates
		double theta = acos(z.z/r);
		double phi = atan2(z.y,z.x);
		dr =  pow( r, powe-1.0)*powe*dr + 1.0;

		// scale and rotate the point
		double zr = pow( r,powe) ;
		theta = theta*powe;
		phi = phi*powe;

		// convert back to cartesian coordinates
		z = zr*float3(sin(theta)*cos(phi), sin(phi)*sin(theta), cos(theta));
		z+=position;
	}
    double dst = 0.5*log(r)*r/dr;
	return float2(iterations,dst*1);
}

Ray CreateCameraRay(float2 uv) {
    float3 origin = mul(_CameraToWorld, float4(0,0,0,1)).xyz;
    float3 direction = mul(_CameraInverseProjection, float4(uv,0,1)).xyz;
    direction = mul(_CameraToWorld, float4(direction,0)).xyz;
    direction = normalize(direction);

    Ray ray;
    ray.origin = origin;
    ray.direction = direction;
    return ray;
}

float3 CalculateNormal(float3 pos)
{
    const float eps = 0.001;

    float3 f1 = float3(
        DstAndIterations(pos + float3(eps, 0.0, 0.0)).y,
        DstAndIterations(pos + float3(0.0, eps, 0.0)).y,
        DstAndIterations(pos + float3(0.0, 0.0, eps)).y
    );
    float3 f2 = float3(
        DstAndIterations(pos - float3(eps, 0.0, 0.0)).y,
        DstAndIterations(pos - float3(0.0, eps, 0.0)).y,
        DstAndIterations(pos - float3(0.0, 0.0, eps)).y
    );

    return normalize(f1 - f2);
}

[numthreads(8,8,1)]
void CSMain (uint3 id : SV_DispatchThreadID)
{
    uint width,height;
    Result.GetDimensions(width, height);

    Result[id.xy] = Source[id.xy];

    float2 uv = id.xy / float2(width,height) * 2 - 1;
    Ray ray = CreateCameraRay(uv);

    float rayDst = 0;
    while(rayDst < maxDst)
    {
        float2 data = DstAndIterations(ray.origin);

        ray.origin += ray.direction * data.y;
        rayDst += data.y;

        if(data.y < minStepSize)
        {
            float3 normal = CalculateNormal(ray.origin);
            float brightness = dot(normal, lightDir);
            Result[id.xy] = float4(brightness, brightness, brightness, 1);
            break;
        }
    }
}