GPU编程:学习Julia曲线绘制cpu和GPU实现方式

Julia 集介绍

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

Julia 集也是一个典型的分形，只是在表达上相当复杂，难以用古典的数学方法描述。文章源自菜鸟学院-https://www.cainiaoxueyuan.com/suanfa/6939.html

Julia 集集合由一个复变函数f(z) = z*z + c  生成，其中c为常数。文章源自菜鸟学院-https://www.cainiaoxueyuan.com/suanfa/6939.html

她是满足某个复数计算函数的所有点构成的边界。对于函数参数的所有取值，生成的边界将形成一种不规则的碎片形状，这是数学中最有趣和最漂亮的形状之一。文章源自菜鸟学院-https://www.cainiaoxueyuan.com/suanfa/6939.html

Julia 集的基本算法原理

通过一个简单的迭代等式对复平面中的点求值，如果在计算某个点时，迭代等式的计算结果是发散的，那么这个点就不属于Julia集合，文章源自菜鸟学院-https://www.cainiaoxueyuan.com/suanfa/6939.html

相反，如果在迭代等式中计算得到的一系列值都位于某个边界范围之内，那么这个点就属于Julia集合。文章源自菜鸟学院-https://www.cainiaoxueyuan.com/suanfa/6939.html

由于c可以是任意值，所以当c取不同的值时，制出的图形也不相同。文章源自菜鸟学院-https://www.cainiaoxueyuan.com/suanfa/6939.html

初学GPU编程，学习Julia曲线绘制的cpu和GPU实现方式，并将实现代码贴下文章源自菜鸟学院-https://www.cainiaoxueyuan.com/suanfa/6939.html

具体代码段已添加相应的备注，可对比学习文章源自菜鸟学院-https://www.cainiaoxueyuan.com/suanfa/6939.html

cpu版本文章源自菜鸟学院-https://www.cainiaoxueyuan.com/suanfa/6939.html

#include <stdio.h>
#include "cuda_runtime.h"
#include <opencv2/core/cuda_devptrs.hpp>
#include "cpu_bitmap.h"
using namespace cv;
using namespace cv::gpu;文章源自菜鸟学院-https://www.cainiaoxueyuan.com/suanfa/6939.html

//懒得用book.h了……
#define DIM 1000文章源自菜鸟学院-https://www.cainiaoxueyuan.com/suanfa/6939.html

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);
}
};文章源自菜鸟学院-https://www.cainiaoxueyuan.com/suanfa/6939.html

int julia(int x, int y){文章源自菜鸟学院-https://www.cainiaoxueyuan.com/suanfa/6939.html

//将点坐标换成复数平面上的坐标，复数平面的远点在图的中心（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);文章源自菜鸟学院-https://www.cainiaoxueyuan.com/suanfa/6939.html

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;
}文章源自菜鸟学院-https://www.cainiaoxueyuan.com/suanfa/6939.html

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;
}
}
}文章源自菜鸟学院-https://www.cainiaoxueyuan.com/suanfa/6939.html

int main()
{
CPUBitmap bitmap(DIM, DIM); //create the appropriate size bitmap image
unsigned char *ptr = bitmap.get_ptr(); //
kernel(ptr); //处理每个点
bitmap.display_and_exit(); //显示处理
}文章源自菜鸟学院-https://www.cainiaoxueyuan.com/suanfa/6939.html

效果图为：文章源自菜鸟学院-https://www.cainiaoxueyuan.com/suanfa/6939.html

GPU版本文章源自菜鸟学院-https://www.cainiaoxueyuan.com/suanfa/6939.html

#include <stdio.h>
#include "cuda_runtime.h"
#include <opencv2/core/cuda_devptrs.hpp>
#include "cpu_bitmap.h"
using namespace cv;
using namespace cv::gpu;文章源自菜鸟学院-https://www.cainiaoxueyuan.com/suanfa/6939.html

//懒得用book.h了……
#define DIM 200文章源自菜鸟学院-https://www.cainiaoxueyuan.com/suanfa/6939.html

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);
}
};文章源自菜鸟学院-https://www.cainiaoxueyuan.com/suanfa/6939.html

//这个函数是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);文章源自菜鸟学院-https://www.cainiaoxueyuan.com/suanfa/6939.html

cuComplex c(-0.8, 0.156); //这里的常数C在书上写着：这个值刚好能生成一张有趣的图片
cuComplex a(jx, jy);文章源自菜鸟学院-https://www.cainiaoxueyuan.com/suanfa/6939.html

int i = 0;
for (i=0; i<200; i++) { //执行迭代运算，判断该像素点是不是属于Julia集，是的返回值是1,
a = a * a + c;
if (a.magnitude2() > 1000)
return 0;
}文章源自菜鸟学院-https://www.cainiaoxueyuan.com/suanfa/6939.html

return 1;
}文章源自菜鸟学院-https://www.cainiaoxueyuan.com/suanfa/6939.html

__global__ void kernel( unsigned char *ptr ){
int x = blockIdx.x;
int y = blockIdx.y;
int offset = x + y * DIM;文章源自菜鸟学院-https://www.cainiaoxueyuan.com/suanfa/6939.html

//循环不见啦！因为都并行了嘛。
//这里是没个像素有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; //代表黑色
}文章源自菜鸟学院-https://www.cainiaoxueyuan.com/suanfa/6939.html

struct DataBlock {
unsigned char *dev_bitmap;
};文章源自菜鸟学院-https://www.cainiaoxueyuan.com/suanfa/6939.html

int main( void ) {文章源自菜鸟学院-https://www.cainiaoxueyuan.com/suanfa/6939.html

DataBlock data;
CPUBitmap bitmap( DIM, DIM, &data );
unsigned char *dev_bitmap;文章源自菜鸟学院-https://www.cainiaoxueyuan.com/suanfa/6939.html

//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);文章源自菜鸟学院-https://www.cainiaoxueyuan.com/suanfa/6939.html

/*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);
}文章源自菜鸟学院-https://www.cainiaoxueyuan.com/suanfa/6939.html

这里DIM设为1000的话，可能会因为内存不够，而不能显示图片，所以我将DIM设为200文章源自菜鸟学院-https://www.cainiaoxueyuan.com/suanfa/6939.html

效果为：文章源自菜鸟学院-https://www.cainiaoxueyuan.com/suanfa/6939.html

具体demo，需要的可 下载文章源自菜鸟学院-https://www.cainiaoxueyuan.com/suanfa/6939.html

---------------------
作者：Nani_xiao
来源：CSDN文章源自菜鸟学院-https://www.cainiaoxueyuan.com/suanfa/6939.html