Julia集是一个著名的分形集,它是复数经过迭代得到的,是满足某个复数计算函数的所有点构成的边界。
算法思想:
通过一个简单的迭代等式对复平面中的点求值,如果在计算某个点时迭代的结果是发散的,那么这个点就不属于Julia集。相反,如果迭代计算得到的一系列值都位于某个边界范围之内,那么这个点就属于Julia集。迭代计算的公式为:
Z(n+1) = Z(n)^2 + C (C为常数)
基于GPU的Julia集通过DIM*DIM个thread计算平面中各点的值,最终用OpenGL显示计算的结果如下:
主要程序代码:
#include <stdio.h>
#include <stdlib.h>
#include <cuda_runtime.h>
#include <cutil_inline.h>
#include "CPUBitmap.h"
//#include <device_launch_parameters.h>
#define DIM 512
struct myComplex{
float r;
float i;
__device__ myComplex(float a, float b): r(a), i(b){}
__device__ float magnitude2(){
return r*r+i*i;
}
__device__ myComplex operator*(const myComplex& a){
return myComplex(r*a.r - i*a.i, i*a.r + r*a.i);
}
__device__ myComplex operator+(const myComplex& a){
return myComplex(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);
myComplex c(-0.8, 0.156);
myComplex a(jx, jy);
int i = 0;
for (i=0; i<200; i++)
{
a = a * a + c;
if (a.magnitude2() > 1000)
{
return 0;
}
}
return 1;
}
__global__ void kernel(unsigned char* ptr)
{
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] = 255 *juliaValue;
ptr[offset*4 + 2] = 255 *juliaValue;
ptr[offset*4 + 3] = 255;
}
int main(int argc, char* argv[])
{
if(!InitCUDA()) {
return 0;
}
CPUBitmap bitmap(DIM, DIM);
unsigned char *dev_bitmap;
cutilSafeCall(cudaMalloc((void**)&dev_bitmap, bitmap.image_size()));
dim3 grid(DIM, DIM);
kernel<<<grid, 1>>>(dev_bitmap);
cutilSafeCall(cudaMemcpy(bitmap.get_ptr(), dev_bitmap, bitmap.image_size(), cudaMemcpyDeviceToHost));
bitmap.display_and_exit();
cutilSafeCall(cudaFree(dev_bitmap));
return 0;
}