Julia曲线绘制-——CPU版本与GPU版本对比

Julia 集

Julia 集是由法国数学家 Gaston Julia 和 Pierre Faton 在发展了复变函数迭代的基础理论后获得的。

Julia 集也是一个典型的分形，只是在表达上相当复杂，难以用古典的数学方法描述。

Julia 集集合由一个复变函数f(z) = z*z + c  生成，其中c为常数。

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

 

Julia 集的基本算法原理

通过一个简单的迭代等式对复平面中的点求值，如果在计算某个点时，迭代等式的计算结果是发散的，那么这个点就不属于Julia集合，

相反，如果在迭代等式中计算得到的一系列值都位于某个边界范围之内，那么这个点就属于Julia集合。

由于c可以是任意值，所以当c取不同的值时，制出的图形也不相同。

初学GPU编程，学习Julia曲线绘制的cpu和GPU实现方式，并将实现代码贴下

具体代码段已添加相应的备注，可对比学习

cpu版本

#include <stdio.h>
#include "cuda_runtime.h"
#include <opencv2/core/cuda_devptrs.hpp>
#include "cpu_bitmap.h"
using namespace cv;
using namespace cv::gpu;
 
//懒得用book.h了……
#define DIM 1000

struct cuComplex {
	float   r;
	float   i;
	cuComplex( float a, float b ) : r(a), i(b)  {}
	float magnitude2( void ) { return r * r + i * i; }
	cuComplex operator*(const cuComplex& a) {
		return cuComplex(r*a.r - i*a.i, i*a.r + r*a.i);
	}
	cuComplex operator+(const cuComplex& a) {
		return cuComplex(r+a.r, i+a.i);
	}
};

int julia(int x, int y){

	//将点坐标换成复数平面上的坐标，复数平面的远点在图的中心（DIM/2, DIM/2）
	//并且使其范围在[-1.0, 1.0]
	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);                       //this is arbitrary,but it happens to yield an interesting
	cuComplex a(jx, jy);                           
	//计算两百次，每次计算后判断有没有超过阈值（这里是1000），超过就发散，返回0
	int i = 0; for(i = 0; i < 200; i++){
		a = a * a + c;
		if(a.magnitude2() > 1000)
			return 0;
	}
	return 1;
}

void kernel(unsigned char * ptr){
	for(int y = 0; y < DIM; y++){
		for(int x = 0; x < DIM; x++){
			int offset = x + y * DIM;
			int juliaValue = julia(x,y);          //julia()返回bool，是集合中的则返回1，否则返回0
			ptr[offset*4 + 0] = 255 * juliaValue; //设置颜色
			ptr[offset*4 + 1] = 0;
			ptr[offset*4 + 2] = 0;
			ptr[offset*4 + 3] = 255;
		}
	}
}



int main()
{
	CPUBitmap bitmap(DIM, DIM);            //create the appropriate size bitmap image
	unsigned char *ptr = bitmap.get_ptr(); //
	kernel(ptr);                           //处理每个点
	bitmap.display_and_exit();             //显示处理
}

效果图为：

GPU版本

#include <stdio.h>
#include "cuda_runtime.h"
#include <opencv2/core/cuda_devptrs.hpp>
#include "cpu_bitmap.h"
using namespace cv;
using namespace cv::gpu;
 
//懒得用book.h了……
#define DIM 200

struct cuComplex {
	float   r;
	float   i;
	__device__ cuComplex( float a, float b ) : r(a), i(b)  {}   //这里一定要记得加__device__!!书里没有加!!
	__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.r + r*a.i);
	}
	__device__ cuComplex operator+(const cuComplex& a) {
		return cuComplex(r+a.r, i+a.i);
	}
};

//这个函数是kernel函数调用的，在device上执行，所以用__device__声明
__device__ int julia( int x, int y ) { 
	const float scale = 1.5; //为了让图像更大，便于查看，选择了放大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); //这里的常数C在书上写着：这个值刚好能生成一张有趣的图片 
	cuComplex a(jx, jy);

	int i = 0;
	for (i=0; i<200; i++) { //执行迭代运算，判断该像素点是不是属于Julia集，是的返回值是1,  
		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 * DIM;

	//循环不见啦！因为都并行了嘛。
	//这里是没个像素有4个字节，所以索引之后要乘以4才能找到对应像素的位置，
	//然后这四个字节表示的是RGBA模型，第一个是代表红色  
	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; //代表黑色
}

struct DataBlock {  
    unsigned char   *dev_bitmap;  
};  


int main( void ) {

	DataBlock data;
	CPUBitmap bitmap( DIM, DIM, &data );
	unsigned char *dev_bitmap;

	//HANDLE_ERROR(cudaMalloc((void**)&dev_bitmap, bitmap.image_size()));
	cudaMalloc((void**)&dev_bitmap, bitmap.image_size());
	dim3 grid(DIM, DIM);                //dim3是cuda内置类型，包含三个int型，x,y,z
	kernel<<<grid, 1>>>(dev_bitmap);

	/*HANDLE_ERROR(cudaMemcpy(bitmap.get_ptr(),
		dev_bitmap,
		bitmap.image_size(),
		cudaMemcpyDeviceToHost));*/
	cudaMemcpy(bitmap.get_ptr(),
		dev_bitmap,
		bitmap.image_size(),
		cudaMemcpyDeviceToHost);
	bitmap.display_and_exit();
	cudaFree(dev_bitmap);
}

这里DIM设为1000的话，可能会因为内存不够，而不能显示图片，所以我将DIM设为200

效果为：

具体demo，需要的可 下载