GPU编程：学习Julia曲线绘制CPU和GPU实现方法

Julia集介绍

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

Julia集也是一种典型的分形，但其表达相当复杂，很难用经典数学方法来描述。

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

她是所有满足某种复杂算术函数的点所形成的边界。对于函数参数的所有值，得到的边界将形成不规则的剪切形状，这是数学中最有趣、最美丽的形状之一。

基本Julia集算法原理

通过简单的迭代方程评估复杂平面中的点。如果迭代方程的计算结果在计算某个点时发散，则说明该点不属于Julia集，

反之，如果迭代方程计算的一系列值在一定限度内范围，则该点属于Julia集。

由于c可以是任意值，所以当c取不同的值时，得到的图形也会不同。

GPU编程初学者，学习Julia曲线绘制的CPU和GPU实现，下面贴出实现代码

具体代码段添加了相应注释，可以对比学习

cpu版本

#include
#include "cuda_runtime.h"
#include
#include "cpu_bitmap。
使用命名空间 cv::gpu ;

//懒得用book.h...
#define DIM 1000

struct cuComplex {
float r;
float i;
cuComplex( float a, float b ) : r(a), i(b) {}
r 浮点数 2 r + i * i ; }
cuComplex 运算符*(const cuComplex& a) {
return cuComplex(r*a.r - i*a.i, i*a.r + r*a.i);
}♸+ cu♝ ​​{
return cuComplex(r+a.r, i+a.i);
}
};

int julia(int x, int y){

//将用复数上的坐标替换点 d飞机。复平面最外点位于图像中心(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); //这是任意的，但它恰好产生一个有趣的
cuComplex a(jx, jy);
//计算两百次。每次计算结束后，判断是否超过阈值（这里是1000）。如果超过，就会发散，返回0
int i = 0； for(i = 0; i a = a * a + c;
if(a.magnitude2() > 1000)
返回 0;
} 返回 1;
}

void kernel(unsigned char * ptr){
for(int y = 0; y for(int x = 0; 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 位图(DIM, DIM); //创建适当大小的位图图像
unsigned char *ptr = bitmap .get_ptr(); //
核心(ptr); //处理每个点
bitmap.display_and_exit(); //屏幕处理
}

渲染为：

GPU版本

#include
#include "cuda_runtime.h"include
#include "cpu_bitmap.h"
using namespace cv ;
using namespace cv::gpu;

//懒得用book..定义DIM 200

struct cuComplex {
float r;
float i;
__device__ cuComplex( float a, float b ) : r(a), i(b) {} //这里记得添加 __device！ ！书上没有添加！！
__device__ floatitude2( 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 运算符+(const cuComplex& a) {
return cuComplex(r+a.r, i+a.i);
}
/ }
/ };在设备上执行，因此声明为 __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 =比例 * (浮点)(DIM/2 - y)/(DIM/2);

cuComplex c(-0.8, 0.156); //这里的常量C书上是这样写的：这个值刚好足以生成一个有趣的图像
cuComplex a(jx, jy);

int i = 0;
for (i=0; i 1000)
return 0;
}❝ }

__global__ void core( unsigned char *ptr ){
int x = blockIdx.x;
int y = blockIdx.y;
int offset = x + y/; 循环消失了！因为它们都是并行运行的。
//这里每个像素有4个字节，所以索引后需要乘以4才能找到对应像素的位置，
//那么这四个字节代表RGBA模型，第一个它代表红色
int juliaValue = julia( x, y );
ptr[偏移*4 + 0] = 255 * juliaValue;
ptr[偏移*4 + 1] = 0; *4 + 2] = 0 ;
ptr[偏移*4 + 3] = 255; //代表黑色
}

struct DataBlock {
unsigned char *dev_bitmap;
} ;

int main( void ) {it CPU;位图( DIM, DIM, &data ) .image_size());
dim3 网格(DIM, DIM); //dim3是cuda的内置类型，包括三个int类型，x,y,z
core(dev_bitmap);

/*HANDLE_ERROR(cudaMemcpy(bitmap.get_ptr(),
dev_bitmap,
bitmap.image_size(),
cudaMemcpyDeviceToHost)❀/get.) ,
dev_bitmap ,
位图.image_size( ) ),
cudaMemcpyDeviceToHost);
bitmap.display_and_exit();
cudaFree(dev_bitmap);如果是1000，可能是内存不够显示图像，所以我把DIM设置为200

效果是：

具体demo，需要的可以下载

---------------- -----
作者：Nani_xiao
来源：CSDN