Welcome to Fractal Forums

using System;
using System.Collections.Generic;
using System.Text;
using System.Drawing;
using System.Drawing.Imaging;
using System.Runtime.InteropServices;

using OpenTK;
using OpenTK.Graphics;
using OpenTK.Graphics.OpenGL;
using OpenTK.Input;



using Cloo;


namespace TKBuddhabrot
{
    class TKBuddhabrot : GameWindow
    {

private static string kernelSource = @"
__kernel void mandelbrot(
  const float deltaReal,
  const float deltaImaginary,
  const float realMin,
  const float imaginaryMin,
  const unsigned int maxIter,
  const unsigned int escapeOrbit,
  const unsigned int hRes,
  __global int* outputi
) {

  int xId = get_global_id(0);
  int yId = get_global_id(1);

  float realPos = realMin + (xId * deltaReal);
  float imaginaryPos = imaginaryMin + (yId * deltaImaginary);
  float real = realPos;
  float imaginary = imaginaryPos;
  float realSquared = real * real;
  float imaginarySquared = imaginary * imaginary;

  int iter = 0;
  while ( (iter < maxIter) && ((realSquared + imaginarySquared) < escapeOrbit) )
  {
    imaginary = (2 * (real * imaginary)) + imaginaryPos;
    real = realSquared - imaginarySquared + realPos;
    realSquared = real * real;
    imaginarySquared = imaginary * imaginary;
    iter++;
  }
  if(iter >= maxIter){
        iter = 0;
  }
  outputi[(yId * hRes) + xId] = iter;
}



";

        ComputePlatform platform;
        ComputeContextPropertyList properties; 
        ComputeContext context;

        Bitmap bmp;

        float realMin, realMax, imaginaryMin, imaginaryMax, deltaReal, deltaImaginary;
        int maxiter, screenSizeInPixel, escapeOrbit;

        static int initialScreenWidth = 800;
        static int initialScreenHeight = 800;



        /// <summary>Creates a window with the specified title.</summary>
        public TKBuddhabrot() : base(initialScreenWidth, initialScreenHeight, GraphicsMode.Default, "TKBuddhabrot")
        {
            VSync = VSyncMode.On;
        }

        /// <summary>Load resources here.</summary>
        /// <param name="e">Not used.</param>
        protected override void OnLoad(EventArgs e)
        {
            base.OnLoad(e);

            //Create Bitmap
            bmp = new Bitmap(ClientRectangle.Width, ClientRectangle.Height);

            //OpenGL Stuff
            GL.ClearColor(0.1f, 0.2f, 0.5f, 0.0f);

            //OpenCL initialisation
            platform = ComputePlatform.Platforms[0];
            Console.WriteLine("Compute platform : " + platform.ToString());

            properties = new ComputeContextPropertyList(platform);            
            context = new ComputeContext(platform.Devices, properties, null, IntPtr.Zero);
            Console.WriteLine("Compute context : " + context.ToString());

            //Mandelbrot Specific
            realMin = -2.25f;
            realMax = 0.75f;
            imaginaryMin = -1.5f;
            imaginaryMax = 1.5f;
            maxiter = 64;
            escapeOrbit = 4;

            deltaReal = (realMax - realMin) / (ClientRectangle.Width - 1);
            deltaImaginary = (imaginaryMax - imaginaryMin) / (ClientRectangle.Height - 1);
            screenSizeInPixel = ClientRectangle.Width * ClientRectangle.Height;

            //OpenCL Buffer
            ComputeBuffer<float> kernelOutput = new ComputeBuffer<float>(context, ComputeMemoryFlags.WriteOnly, screenSizeInPixel);

            //Build OpenCL kernel
            ComputeProgram program = new ComputeProgram(context, new string[] { kernelSource });
            program.Build(null, null, null, IntPtr.Zero);
            ComputeKernel kernel = program.CreateKernel("mandelbrot");

            //OpenCL args
            //  const float deltaReal,
            //  const float deltaImaginary,
            //  const float realMin,
            //  const float imaginaryMin,
            //  const unsigned int maxIter,
            //  const unsigned int escapeOrbit,
            //  const unsigned int hRes,
            //  __global int* outputi

            kernel.SetValueArgument<float>(0, deltaReal);
            kernel.SetValueArgument<float>(1, deltaImaginary);
            kernel.SetValueArgument<float>(2, realMin);
            kernel.SetValueArgument<float>(3, imaginaryMin);
            kernel.SetValueArgument<int>(4, maxiter);
            kernel.SetValueArgument<int>(5, escapeOrbit);
            kernel.SetValueArgument<int>(6, ClientRectangle.Width);
            kernel.SetMemoryArgument(7, kernelOutput);

            //Execute
            ComputeCommandQueue commands = new ComputeCommandQueue(context, context.Devices[0], ComputeCommandQueueFlags.None);
            ComputeEventList events = new ComputeEventList();
            commands.Execute(kernel, null, new long[] { ClientRectangle.Width, ClientRectangle.Height }, null, events);

            //Get result
            int[] kernelResult = new int[screenSizeInPixel];
            GCHandle kernelResultHandle = GCHandle.Alloc(kernelResult, GCHandleType.Pinned);

            commands.Read(kernelOutput, false, 0, screenSizeInPixel, kernelResultHandle.AddrOfPinnedObject(), events);
            commands.Finish();

            //Finish openCL stuff
            kernelResultHandle.Free();
            
            int maxfound = 0;
            foreach (int iter in kernelResult)
            {
                if (iter > maxfound) maxfound = iter;
            }

            //Use the result
            int x, y;
            for (x = 0; x < bmp.Width; x++)
            {
                for (y = 0; y < bmp.Height; y++)
                {
                    Color c = Color.FromArgb((int)((float)kernelResult[x + y * bmp.Width] / maxfound * 255.0),
                        (int)((float)kernelResult[x + y * bmp.Width] / maxfound * 255.0),
                        (int)((float)kernelResult[x + y * bmp.Width] / maxfound * 255.0)
                        );

                    bmp.SetPixel(x, y, c);
                }
            }

            BitmapData bmp_data = bmp.LockBits(new Rectangle(0, 0, bmp.Width, bmp.Height), ImageLockMode.ReadOnly, System.Drawing.Imaging.PixelFormat.Format32bppArgb);

            GL.TexImage2D(TextureTarget.Texture2D, 0, PixelInternalFormat.Rgba, bmp_data.Width, bmp_data.Height, 0,
                OpenTK.Graphics.OpenGL.PixelFormat.Bgra, PixelType.UnsignedByte, bmp_data.Scan0);

            bmp.UnlockBits(bmp_data);
            GL.TexParameter(TextureTarget.Texture2D, TextureParameterName.TextureMinFilter, (int)TextureMinFilter.Linear);
            GL.TexParameter(TextureTarget.Texture2D, TextureParameterName.TextureMagFilter, (int)TextureMagFilter.Linear);

            Console.WriteLine("done");

        }

        /// <summary>
        /// Called when your window is resized. Set your viewport here. It is also
        /// a good place to set up your projection matrix (which probably changes
        /// along when the aspect ratio of your window).
        /// </summary>
        /// <param name="e">Not used.</param>
        protected override void OnResize(EventArgs e)
        {
            base.OnResize(e);

            GL.Viewport(ClientRectangle.X, ClientRectangle.Y, ClientRectangle.Width, ClientRectangle.Height);

            Matrix4 projection = Matrix4.CreatePerspectiveFieldOfView((float)Math.PI / 4, Width / (float)Height, 1.0f, 64.0f);
            GL.MatrixMode(MatrixMode.Projection);
            GL.LoadMatrix(ref projection);
        }

        /// <summary>
        /// Called when it is time to setup the next frame. Add you game logic here.
        /// </summary>
        /// <param name="e">Contains timing information for framerate independent logic.</param>
        protected override void OnUpdateFrame(FrameEventArgs e)
        {
            base.OnUpdateFrame(e);

            if (Keyboard[Key.Escape])
                Exit();
        }

        /// <summary>
        /// Called when it is time to render the next frame. Add your rendering code here.
        /// </summary>
        /// <param name="e">Contains timing information.</param>
        protected override void OnRenderFrame(FrameEventArgs e)
        {
            base.OnRenderFrame(e);

            GL.Clear(ClearBufferMask.ColorBufferBit | ClearBufferMask.DepthBufferBit);

            GL.Disable(EnableCap.DepthTest);
            GL.Enable(EnableCap.Texture2D);
            GL.Enable(EnableCap.Blend);
            GL.BlendFunc(BlendingFactorSrc.SrcAlpha, BlendingFactorDest.OneMinusSrcAlpha);


            GL.MatrixMode(MatrixMode.Projection);

            GL.LoadIdentity();
            GL.Ortho(0, ClientRectangle.Width, ClientRectangle.Height, 0, 0, 1);
          
  
            GL.Begin(BeginMode.Polygon);
                GL.TexCoord2(0.0, 1.0);
                GL.Vertex2(0, 0);

                GL.TexCoord2(1.0, 1.0);
                GL.Vertex2(ClientRectangle.Width, 0);

                GL.TexCoord2(1.0, 0.0);
                GL.Vertex2(ClientRectangle.Width, ClientRectangle.Height);

                GL.TexCoord2(0.0, 0.0);
                GL.Vertex2(0,ClientRectangle.Height);
            GL.End();
            
            SwapBuffers();
        }

        /// <summary>
        /// The main entry point for the application.
        /// </summary>
        [STAThread]
        static void Main()
        {
            // The 'using' idiom guarantees proper resource cleanup.
            // We request 30 UpdateFrame events per second, and unlimited
            // RenderFrame events (as fast as the computer can handle).
            using (TKBuddhabrot game = new TKBuddhabrot())
            {
                Console.WriteLine("Display device List (may be usefull for debug)");
                foreach (DisplayDevice device in DisplayDevice.AvailableDisplays)
                {

                    Console.WriteLine("-------------");
                    Console.WriteLine("is primary : " + device.IsPrimary);
                    Console.WriteLine("bound : " + device.Bounds);
                    Console.WriteLine("Refresh rate : " + device.RefreshRate);
                    Console.WriteLine("bpp : " + device.BitsPerPixel);
                    //foreach (DisplayResolution res in device.AvailableResolutions) { Console.WriteLine(res); }

                }
                Console.WriteLine("-------------");
                game.Run(30.0);
            }
        }


    }
}

bool isInMSet(
  float cr, 
  float ci, 
  const unsigned int maxIter,
  const float escapeOrbit
)
{
    int iter = 0;
    float zr = 0.0;
    float zi = 0.0;
    float zr2 = zr * zr;
    float zi2 = zi * zi;
    float temp = 0.0;

    //Check if c is in the 2nd order period bulb.
    if( sqrt( ((cr+1.0) * (cr+1.0)) + (ci * ci) ) < 0.25 )
    {
        return true;
    }

    //Check if c is in the main cardiod
    //IF ((ABS( 1.0 - SQRT(1-(4*c)) ))  < 1.0 ) THEN RETURN TRUE (main cardioid)

     while ( (iter < maxIter) && ((zr2 + zi2) < escapeOrbit) )
    {
        temp = zr * zi;
        zr2 = zr * zr;
        zi2 = zi * zi;
        zr = zr2 - zi2 + cr;
        zi = temp + temp + ci;
        iter++;
    }

    if ( iter < maxIter )
    {
        return false;
    } else {
        return true;
    }

}    

__kernel void mandelbrot(
  const float realMax,
  const float imaginaryMax,
  const float realMin,
  const float imaginaryMin,
  const unsigned int maxIter,
  const unsigned int escapeOrbit,
  const unsigned int hRes,
  __global int* outputi
) {

  const int xId = get_global_id(0);
  const int yId = get_global_id(1);
  const int maxX = get_global_size(0);
  const int maxY = get_global_size(1);

  float deltaReal = (realMax - realMin) / (maxX - 1);
  float deltaImaginary = (imaginaryMax - imaginaryMin) / (maxY - 1);

  float realPos = realMin + (xId * deltaReal);
  float imaginaryPos = imaginaryMin + (yId * deltaImaginary);

  float real = realPos;
  float imaginary = imaginaryPos;

  if(isInMSet(real, imaginary, maxIter, escapeOrbit) == false)
  {

    int iter = 0;
    float zr = 0.0;
    float zi = 0.0;
    float zr2 = zr * zr;
    float zi2 = zi * zi;
    float cr = real;
    float ci = imaginary;
    float temp = 0.0;
    
    while ((iter < maxIter))
    {
        temp = zr * zi;
        zr2 = zr * zr;
        zi2 = zi * zi;
        zr = zr2 - zi2 + cr;
        zi = temp + temp + ci;
        int x = (int)(maxX * (zr - realMin) / (realMax - realMin));
        int y = (int)(maxY * (zi - imaginaryMin) / (imaginaryMax - imaginaryMin));
        
        if( (x > 0) && (y > 0) && (x < maxX) && (y < maxY) && (iter > 2))
        {
            outputi[(y * hRes) + x] += 1;
        }
        iter++;
    }

  
  }

}

//A function to check if the choosen point is in the mandelbrot set
bool isInMSet(
  float cr, 
  float ci, 
  const unsigned int maxIter,
  const float escapeOrbit
)
{
    int iter = 0;
    float zr = 0.0;
    float zi = 0.0;
    float zr2 = zr * zr;
    float zi2 = zi * zi;
    float temp = 0.0;

    //Quick rejection check if c is in the 2nd order period bulb.
    if( sqrt( ((cr+1.0) * (cr+1.0)) + (ci * ci) ) < 0.25 )
    {
        return true;
    }

    //Quick rejection check if c is in the main cardiod
    //IF ((ABS( 1.0 - SQRT(1-(4*c)) ))  < 1.0 ) THEN RETURN TRUE (main cardioid)
    float tempi = ci*(-4.0);
    float tempr = 1.0 - cr*4.0;
    float theta = atan2(tempi,tempr)/2.0;
    float r = pow((tempr*tempr + tempi*tempi),0.25);
    tempr = 1.0 - r * cos(theta);
    tempi = -r * sin(theta);
    if( (tempr * tempr + tempi * tempi) < 1.0)
    {
        return true;
    }


    //Bruteforce check if c is escaping escapeOrbit (with a good old iteration up to maxIter)
    while( (iter < maxIter) && ((zr2 + zi2) < escapeOrbit) )
    {
        temp = zr * zi;
        zr2 = zr * zr;
        zi2 = zi * zi;
        zr = zr2 - zi2 + cr;
        zi = temp + temp + ci;
        iter++;
    }

    if ( iter < maxIter )
    {
        return false;
    } else {
        return true;
    }

}    

//Main kernel

__kernel void buddhabrot(
  const float realMax,
  const float imaginaryMax,
  const float realMin,
  const float imaginaryMin,
  const unsigned int maxIter,
  const unsigned int escapeOrbit,
  const unsigned int hRes,
  const float offset,
  __global int* outputi
) {

  const int xId = get_global_id(0);
  const int yId = get_global_id(1);
  const int offsetStep = get_global_id(2);

  const int maxX = get_global_size(0);
  const int maxY = get_global_size(1);

  const float deltaReal = (realMax - realMin) / (maxX - 1);
  const float deltaImaginary = (imaginaryMax - imaginaryMin) / (maxY - 1);

  float cr = realMin + (xId * deltaReal) + (offsetStep * offset);
  float ci = imaginaryMin + (yId * deltaImaginary);

  int iter = 0;
  float zr = 0.0;
  float zi = 0.0;
  float zr2 = zr * zr;
  float zi2 = zi * zi;
  float temp = 0.0;
  int x, y;
    
  if(isInMSet(cr, ci, maxIter, escapeOrbit) == false)
  {
    iter = 0;
    zr = 0.0;
    zi = 0.0;
    zr2 = zr * zr;
    zi2 = zi * zi;
    temp = 0.0;
    while ((iter < maxIter) && ((zr2 + zi2) < escapeOrbit) )
    {
      temp = zr * zi;
      zr2 = zr * zr;
      zi2 = zi * zi;
      zr = zr2 - zi2 + cr;
      zi = temp + temp + ci;
      x = (maxX * (zr - realMin) / (realMax - realMin));
      y = (maxY * (zi - imaginaryMin) / (imaginaryMax - imaginaryMin));
        
      if( (x > 0) && (y > 0) && (x < maxX) && (y < maxY) && (iter > 2))
      {
        outputi[(y * hRes) + x] += 1;
      }
      iter++;
    } //EndWhile
  } //EndIf
} //EndKernel

__kernel void xorshift(
    uint s1,
    uint s2,
    uint s3,
    uint s4,
    const int bufferSize,
    __global uint* outputUInt,
    __global float* outputFloat
)
{
    uint st;

    for(int i=0; i < bufferSize; i++)
    {
        st = s1 ^ (s1 << 11);
        s1 = s2;
        s2 = s3;
        s3 = s4;
        s4 = s4 ^ (s4 >> 19) ^ ( st ^ (st >> 18));
        outputUInt[i] = s4;
        outputFloat[i] = s4 / 4294967295.0; 
    }
}

float zx = -0.0423594f, zy = -0.985749f, zoom = 120.f;

// the formula for computing the real and imaginary bounds is like this, ih/iw is the height to width ratio.
float zrmin = zx - 1.0f / zoom, zrmax = zx + 1.0f / zoom;
float zimin = zy - 1.0f / zoom * ih/iw, zimax = zy + 1.0f / zoom * ih/iw;

The following error has occurred during XML parsing:

File: D:\projets\Buddhabrot-cuda\Buddhabrot\Buddhabrot_vc90.vcproj
Line: 22
Column: 4
Error Message:
Custom build rules file 'd:\projets\common\Cuda.rules' was not found or failed to load.
The file 'D:\projets\Buddhabrot-cuda\Buddhabrot\Buddhabrot_vc90.vcproj' has failed to load.

        while( (iter < maxIter) && ((zr*zr+zi*zi) < escapeOrbit) )
        {
            temp = zr * zi;
            zr = zr*zr - zi*zi + cr;
            zi = temp + temp + ci;
            //etc ....
        }

        while( (iter < maxIter) && ((zr2+zi2) < escapeOrbit) )
        {
            temp = zr * zi;
            zr2 = zr * zr;
            zi2 = zi * zi;
            zr = zr2 - zi2 + cr;
            zi = temp + temp + ci;
            //etc ....
        }

//Check if choosen point is in MSet
bool isInMSet(
    const float2 c,
    const uint minIter,
    const uint maxIter,
    const float escapeOrbit)
{
    int iter = 0;
    float2 z = 0.0;

    if( !(((c.x-0.25)*(c.x-0.25) + (c.y * c.y))*(((c.x-0.25)*(c.x-0.25) + (c.y * c.y))+(c.x-0.25)) < 0.25* c.y * c.y))  //main cardioid
    {
        if( !((c.x+1.0) * (c.x+1.0) + (c.y * c.y) < 0.0625))            //2nd order period bulb
        {
            if (!(( ((c.x+1.309)*(c.x+1.309)) + c.y*c.y) < 0.00345))    //smaller bulb left of the period-2 bulb
            {
                if (!((((c.x+0.125)*(c.x+0.125)) + (c.y-0.744)*(c.y-0.744)) < 0.0088))      // smaller bulb bottom of the main cardioid
                {
                    if (!((((c.x+0.125)*(c.x+0.125)) + (c.y+0.744)*(c.y+0.744)) < 0.0088))  //smaller bulb top of the main cardioid
                    {
                        while( (iter < maxIter) && (z.x*z.x + z.y*z.y < escapeOrbit) )      //Bruteforce check  
                        {
                            z = (float2)(z.x * z.x - z.y * z.y, (z.x * z.y * 2.0)) + c;
                            iter++;
                        }
                        if( (iter > minIter) && (iter < maxIter))
                        {
                            return false;
                        }
                    }
                }
            }
        }
    }
    return true;
}

//Main kernel
__kernel void buddhabrot(
    const float realMin,
    const float realMax,
    const float imaginaryMin,
    const float imaginaryMax,
    const uint  minIter,
    const uint  maxIter,
    const uint  width,
    const uint  height,
    const float escapeOrbit,
    const uint4 minColor,
    const uint4 maxColor,
    __global float2* randomXYBuffer,
    __global uint4*  outputBuffer)
{

    float2 rand = randomXYBuffer[get_global_id(0)];    

    const float deltaReal = (realMax - realMin);
    const float deltaImaginary = (imaginaryMax - imaginaryMin);

    //mix(a,b,c) = a + (b-a)*c //(c must be in the range 0.0 ... 1.0
    float2 c = (float2)(mix(realMin, realMax, rand.x) , mix(imaginaryMin, imaginaryMax, rand.y));

    if( isInMSet(c, minIter, maxIter, escapeOrbit) == false)
    {
        int x, y;
        int iter = 0;
        float2 z = 0.0;
    
        while( (iter < maxIter) && ((z.x*z.x+z.y*z.y) < escapeOrbit) )
        {
            z = (float2)(z.x * z.x - z.y * z.y, (z.x * z.y * 2.0)) + c;
            x = (width * (z.x - realMin) / deltaReal);
            y = (height * (z.y - imaginaryMin) / deltaImaginary);

            if( (iter > minIter) && (x>0) && (y>0) && (x<width) && (y<height) )
            {
                if( (iter > minColor.x) && (iter < maxColor.x) ) { outputBuffer[x + (y * width)].x++; }
                if( (iter > minColor.y) && (iter < maxColor.y) ) { outputBuffer[x + (y * width)].y++; }
                if( (iter > minColor.z) && (iter < maxColor.z) ) { outputBuffer[x + (y * width)].z++; }
            }
            iter++;
        }
    }
}

__kernel void xorshift(
    uint s1,
    uint s2,
    const int bufferSize,
    __global float2* randomXYBuffer
)
{
    uint st;

    for(int i=0; i < bufferSize; i++)
    {
        st = s1 ^ (s1 << 11);
        s1 = s2;
        s2 = s2 ^ (s2 >> 19) ^ ( st ^ (st >> 18));
        randomXYBuffer[i] = (float2)((float)st / UINT_MAX,(float)s1 / UINT_MAX);
    }
}