2015-08-15 27 views
1

我刚写完一个呈现Mandelbrot集的图像的cuda程序。我建立的方式是,您将创建图像的函数传递给每个单位的像素的比例以及复平面中图像中心的x和y坐标。我想从许多帧创建一个深度缩放电影,我需要我的程序能够自动确定“有趣”的东西将会发生的中心(而不是放大到只有一种颜色的区域)。我应该如何选择要放大的坐标。如何确定在Mandelbrot集上缩放的好中心

这是我的代码,如果有人感兴趣。

#include <iostream> 
#include <thrust/complex.h> 
#include <cuda.h> 
#include <cassert> 
#include <cstdio> 
#include <algorithm> 

typedef double real; 

inline void cuda_error(cudaError_t code, const char* lbl) 
{ 
    if(code != cudaSuccess) 
    { 
     std::cerr << lbl << " : " << cudaGetErrorString(code) << std::endl; 
     exit(1); 
    } 
} 

__global__ void mandelbrot_kernel(unsigned char* pix, real cx, real cy, real pix_scale, size_t w, size_t h, int iters) 
{ 
    cy = -cy; 
    real sx = cx - (w * pix_scale)/2; 
    real sy = cy - (w * pix_scale)/2; 

    size_t x = (size_t)blockIdx.x * blockDim.x + threadIdx.x; 
    size_t y = (size_t)blockIdx.y * blockDim.y + threadIdx.y; 
    if(x >= w || y >= h) 
     return; 

    thrust::complex<real> c(sx + pix_scale * x, sy + pix_scale * y); 
    thrust::complex<real> z(0, 0); 
    int i = 0; 
    for(; i < iters && thrust::abs(z) < 2; ++i) 
     z = z * z + c; 

    real scale = 255.0/(real)iters; 
    size_t q = 3 * (w * y + x); 
    pix[q] = i * scale; 
    pix[q + 1] = 255 * sinf(z.imag()); 
    pix[q + 2] = 255 * sinf(z.real()); 
} 

void shade_mandelbrot(unsigned char* pix, real* devs, real cx, real cy, real pix_scale, int w, int h, int iters) 
{ 
    dim3 blockDim(16, 16); 
    dim3 gridDim((w + 15)/16, (h + 15)/16); 
    mandelbrot_kernel<<<gridDim, blockDim>>>(pix, cx, cy, pix_scale, w, h, iters); 
} 

void ppm_write(FILE* f, unsigned char* pix, int w, int h) 
{ 
    assert(fprintf(f, "P6 %d %d 255\n", w, h) > 0); 
    size_t sz = 3 * (size_t)w * (size_t)h; 
    assert(fwrite(pix, 1, sz, f) == sz); 
} 

int main() 
{ 
    int dim = 2000; 
    int w = dim; 
    int h = dim; 
    int imgs = 200; 
    int iters = 1024; 

    real cx = -0.7463, cy = 0.1102; 
    cuda_error(cudaSetDevice(0), "Set Device"); 
    unsigned char* pix_buffers[2]; 
    real* dev_buffers[2]; 

    cuda_error(cudaHostAlloc(pix_buffers, 3 * sizeof(unsigned char) * w * h, 0), "Host Alloc 1"); 
    cuda_error(cudaHostAlloc(pix_buffers + 1, 3 * sizeof(unsigned char) * w * h, 0), "Host Alloc 2"); 

    real scale = 8.0/w; 
    shade_mandelbrot(pix_buffers[0], dev_buffers[0], cx, cy, scale, w, h, iters); 
    for(int i = 0; i < imgs; i++) 
    { 
     cuda_error(cudaDeviceSynchronize(), "Sync"); 

     std::cout << scale << std::endl; 
     if(i < (imgs - 1)) 
      shade_mandelbrot(pix_buffers[(i + 1) % 2], dev_buffers[(i + 1) % 2], cx, cy, scale *= 0.97, w, h, 255); 
     char fn[100]; 
     sprintf(fn, "/media/chase/3161D67803D8C5BE/Mandelbroght/image%06d.ppm", i); 
     puts(fn); 
     FILE* f = fopen(fn, "w"); 
     assert(f); 

     ppm_write(f, pix_buffers[i % 2], w, h); 
     fclose(f); 
    } 


    cuda_error(cudaFreeHost(pix_buffers[0]), "Host Free 1"); 
    cuda_error(cudaFreeHost(pix_buffers[1]), "Host Free 2"); 
    return 0; 
} 

回答

3

点具有高迭代计数(但不等于iters)时,因为它们接近设定边界离开内环将呈现有趣的行为。您可以随机选取点,通过算法运行并使用计数最高的点作为中心。如果您随机抽取几个点,以最高的迭代次数计算点,在其周围生成几个点,查看是否获得更高的迭代次数,并重复这些点中的最佳点,则可能会更快地获得结果。

+0

Succinct;写得好。做得好。 – duffymo

0

那么我想出了一个很好的想法。我所做的是计算图像中每个16x16平方的熵。然后,我只是放大该图像的最大熵区域。

__global__ void entropy_kernel(unsigned char* pix, float* entropy, size_t w, size_t h) 
{ 
    __shared__ float probs[256]; 
    __shared__ float e; 

    if(threadIdx.x == 0 && threadIdx.y == 0) 
    { 
     e = 0; 
     for(int i = 0; i < 256; i++) 
      probs[i] = 0; 
    } 

    __syncthreads(); 

    int x = blockIdx.x * ENTROPY_BLOCK_DIM + threadIdx.x; 
    int y = blockIdx.y * ENTROPY_BLOCK_DIM + threadIdx.y; 
    int px = pix[3 * (y * w + x)]; 
    float p = 1.0/(float)(ENTROPY_BLOCK_DIM * ENTROPY_BLOCK_DIM); 
    atomicAdd(probs + px, p); 

    __syncthreads(); 

    p = probs[px]; 
    if(p) atomicAdd(&e, p * log10f(p)); 

    __syncthreads(); 

    if(threadIdx.x == 0 && threadIdx.y == 0) 
    { 
     entropy[blockIdx.y * gridDim.x + blockIdx.x] = -e; 
    } 
} 

结果如何。 https://www.youtube.com/watch?v=mtxbdoJBA0Q