2015-08-15 28 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