GPU实现julia集曲线

         julia集是满足某个复数计算函数的所有点构成的边界。对于函数参数的所有取值，生成的边界将形成一种不规则的碎片形状，这是数学中最有趣和最漂亮的形状之一。

         生成julia的算法非常简单，Julia集的基本算法是,通过一个简单的迭代等式对复平面中的点求值，如果在计算某个点时，迭代等式的计算结果是发散的，那么这个点就不属于Julia集合，相反，如果在迭代等式中计算得到的一系列值都位于某个边界范围之内，那么这个点就属于Julia集合。迭代等式为：

       Z(n+1)=Z(n)*Z(n)+C

       基于GPU的Julia集的算法如下：

           

#include<stdlib.h>
#include<stdio.h>
#include<string.h>
#include<math.h>
#include<cutil.h>
#include<cpu_bitmap.h>// 该头文件打不开
#define DIM 128
//基于GPU的Julia集
/*__device__修饰符表示代码将在GPU上而不是主机上运行。只能从其他__device__
函数或者从__gloabl__函数中调用它们。
*/
struct cuComplex
{
      float r;
      float i;
       cuComplex(float a,float b):r(a),i(b){}

      //成员函数，模的平方
      __device__ float magnitude2(void)

     {        return r*r+i*i;     }

         //成员函数，复数乘法的重载

    __device__ cuComplex operator*(const cuComplex&a)
      {          return cuComplex(r*a.r-i*a.i,i*a.i+r*a.i);     }

   //成员函数，复数加法的重载
     __device__ cuComplex operator+(const cuComplex& a)
    {         return cuComplex(r=a.r,i+a.i);     }
};

//函数功能，将像素转换为得数

__device__ int julia(int x,int y)
{
    const float scale=1.5;//放大系数

    float jx=scale*(float)(DIM/2-x)/(DIM/2);
    float jy=scale*(float)(DIM/2-y)/(DIM/2);

    cuComplex c(-0.8,0.156);//迭代函数中常量的定义

    cuComplex a(jx,jy);//迭代函数中Z变量的初值

    int i=0;
    for(i=0;i<200;i++)
    {
      a=a*a+c;
       if(a.magnitude2()>1000)//模的平方大于1000的话，则说明发散，不属于julia集

        return 0;
    }
    return 1;//Z属于julia集，则返加1

}
/*第一步：将kernel声明为一个__global__ 类型的函数。线程格每一维的
大小与图像每一维大小是相等的，因些在（0，0）
和（DIM-1，DIM-1）之间的每个像素点都能获得一个线程块
第二步：得到输出 缓冲区ptr中的线性偏移，这个偏移是通中另一个
内置变量girdDim来计算的，对所有的线程块来说，gidDim是一个常数，
用来保存线程格第一维的大小.在示例中，gridDim的值是(DIM,DIM), 因此，
将行索引乘以线程格的宽，再加上列索引，就得到了ptr中的唯一索引，
其取值范围为（DIM*DIM-1）。
最后分析判断某个点是否属于Julia集的代码。
*/
__global__ void kernel(unsigned char *ptr)
{
   //将threadIdx/BlockIdx映身到像素位置
    int x=blockIdx.x;
    int y=blockIdx.y;
    int offset=x+y*gridDim.x;

    //计算这个位置上的值
    int juliaValue=julia(x,y);
 
     ptr[offset*4+0]=255*juliaValue;
     ptr[offset*4+1]=0;
     ptr[offset*4+2]=0;
     ptr[offset*4+3]=255;
}
int main(void)
{
       CPUBitmap bitmap(DIM,DIM);//定义位图对象
       unsigned char *dev_bitmap;//用来保存设备上数据的副本。

       //申请内存空间
       HANDLE_ERROR(cudaMalloc((void**)&dev_bitmap,bitmap.image_size()));

       dim3 grid(DIM,DIM);//声明二维的线程格，每个线程都要执行kernel函数的一个副本。
       kernel<<<grid,1>>>(dev_bitmap);//将dim3的变量grid传递给CUDA运行时
 
      //将处理结果复制回主机，复制方向中指定为从设备到主机
      HANDLE_ERROR(cudaMemcpy(bitmap.get_ptr(),dev_bitmap,bitmap.image_size(),cudaMemcpyDeviceToHost));
 
      bitmap.display_and_exit();
      cudaFree(dev_bitmap);
}

        

 

         程序写完了，但在编译时，由于不认cpu_bitmap头文件，所以本程序没能运行。