Instantly share code, notes, and snippets.

Embed
What would you like to do?
/* *
* Copyright 1993-2012 NVIDIA Corporation. All rights reserved.
*
* Please refer to the NVIDIA end user license agreement (EULA) associated
* with this source code for terms and conditions that govern your use of
* this software. Any use, reproduction, disclosure, or distribution of
* this software and related documentation outside the terms of the EULA
* is strictly prohibited.
*/
#include <stdio.h>
#include <stdlib.h>
#include <vector_types.h>
#include <cuComplex.h>
#include <png.h>
#include <string.h>
/**
* This macro checks return value of the CUDA runtime call and exits
* the application if the call failed.
*/
#define CUDA_CHECK_RETURN(value) { \
cudaError_t _m_cudaStat = value; \
if (_m_cudaStat != cudaSuccess) { \
fprintf(stderr, "Error %s at line %d in file %s\n", \
cudaGetErrorString(_m_cudaStat), __LINE__, __FILE__); \
exit(1); \
} }
static const int DIMX = (4096);
static const int DIMY = (4096);
#define MAXITERATIONS (256)
#define TYPEFLOAT
#ifdef TYPEFLOAT
#define TYPE float
#define cTYPE cuFloatComplex
#define cMakecuComplex(re,i) make_cuFloatComplex(re,i)
#endif
#ifdef TYPEDOUBLE
#define TYPE double
#define cMakecuComplex(re,i) make_cuDoubleComplex(re,i)
#endif
__host__ __device__ static __inline__ cuDoubleComplex cuCexp(cuDoubleComplex x)
{
double factor = exp(x.x);
return make_cuDoubleComplex(factor * cos(x.y), factor * sin(x.y));
}
__host__ __device__ static __inline__ cuFloatComplex cuCexp(cuFloatComplex x)
{
float factor = exp(x.x);
return make_cuFloatComplex(factor * cos(x.y), factor * sin(x.y));
}
cTYPE c0;
__device__ cTYPE juliaFunctor(cTYPE p,cTYPE c){
return cuCaddf(cuCmulf(p,p),c);
//cTYPE a = cMakecuComplex(-1.0,0);
//return cuCaddf(cuCmulf(cuCexp(p),p),c);
}
__device__ int evolveComplexPoint(cTYPE p,cTYPE c){
int it =1;
while(it <= MAXITERATIONS && cuCabsf(p) <=4){
p=juliaFunctor(p,c);
it++;
}
//printf("re=%f,im=%f, it=%i\n",p.y,p.y,it);
return it;
}
//To adjust the vector to be inside of a circle in the complex plane, radius
//of 'scale' and centered on 'center' you need to adjust by the size of the
//picture.
#define moveX (0)
#define moveY (0)
; //you can change these to zoom and change position
__device__ cTYPE convertToComplex(int x, int y,float zoom){
// TYPE factor = sqrt((DIMX/2.0))* sqrt((DIMX/2.0)) + (DIMY/2)*(DIMY/2) ;
TYPE jx = 1.5 * (x - DIMX / 2) / (0.5 * zoom * DIMX) + moveX;
TYPE jy = (y - DIMY / 2) / (0.5 * zoom * DIMY) + moveY;
return cMakecuComplex(jx,jy);
}
__global__ void computeJulia(int* data,cTYPE c,float zoom){
int i = blockIdx.x * blockDim.x + threadIdx.x;
int j = blockIdx.y * blockDim.y + threadIdx.y;
if(i<DIMX && j<DIMY){
cTYPE p = convertToComplex(i,j,zoom);
data[i*DIMY+j] = evolveComplexPoint(p,c);
//printf("i=%i,j=%i: re=%f,im=%f, it=%i\n",i,j,p.y,p.y,data[i*DIMY+j]);
}
}
#define MAX_DWELL (MAXITERATIONS)
///** gets the color, given the dwell (on host) */
//#define CUT_DWELL (MAX_DWELL/4 )
//void dwell_color(int *r, int *g, int *b, int dwell) {
// // black for the Mandelbrot set
//
// *r=dwell%128;
// *g=255;
// *b=(dwell<=MAXITERATIONS)*255;
//
//} // dwell_color
#define CUT_DWELL (MAX_DWELL / 4)
void dwell_color(int *r, int *g, int *b, int dwell) {
// black for the Mandelbrot set
if(dwell >= MAX_DWELL) {
*r = *g = *b = 0;
} else {
// cut at zero
if(dwell < 0)
dwell = 0;
if(dwell <= CUT_DWELL) {
// from black to blue the first half
*r = *g = 0;
*b = 128 + dwell * 127 / (CUT_DWELL);
} else {
// from blue to white for the second half
*b = 255;
*r = *g = (dwell - CUT_DWELL) * 255 / (MAX_DWELL - CUT_DWELL);
}
}
} // dwell_color
/** save the dwell into a PNG file
@remarks: code to save PNG file taken from here
(error handling is removed):
http://www.labbookpages.co.uk/software/imgProc/libPNG.html
*/
void save_image(const char *filename, int *dwells, int w, int h) {
png_bytep row;
FILE *fp = fopen(filename, "wb");
png_structp png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, 0, 0, 0);
png_infop info_ptr = png_create_info_struct(png_ptr);
// exception handling
setjmp(png_jmpbuf(png_ptr));
png_init_io(png_ptr, fp);
// write header (8 bit colour depth)
png_set_IHDR(png_ptr, info_ptr, w, h,
8, PNG_COLOR_TYPE_RGB, PNG_INTERLACE_NONE,
PNG_COMPRESSION_TYPE_BASE, PNG_FILTER_TYPE_BASE);
// set title
png_text title_text;
title_text.compression = PNG_TEXT_COMPRESSION_NONE;
title_text.key = "Title";
title_text.text = "Mandelbrot set, per-pixel";
png_set_text(png_ptr, info_ptr, &title_text, 1);
png_write_info(png_ptr, info_ptr);
// write image data
row = (png_bytep) malloc(3 * w * sizeof(png_byte));
for (int y = 0; y < h; y++) {
for (int x = 0; x < w; x++) {
int r, g, b;
dwell_color(&r, &g, &b, dwells[y * w + x]);
row[3 * x + 0] = (png_byte)r;
row[3 * x + 1] = (png_byte)g;
row[3 * x + 2] = (png_byte)b;
}
png_write_row(png_ptr, row);
}
png_write_end(png_ptr, NULL);
fclose(fp);
png_free_data(png_ptr, info_ptr, PNG_FREE_ALL, -1);
png_destroy_write_struct(&png_ptr, (png_infopp)NULL);
free(row);
} // save_image
/** a useful function to compute the number of threads */
int divup(int x, int y) { return x / y + (x % y ? 1 : 0); }
char* IMAGE_PATH ="./julia";
static const int BLOCKX = 32;
static const int BLOCKY = 32;
/**
* Host function that prepares data array and passes it to the CUDA kernel.
*/
//# c = 1j # dentrite fractal
//# c = -0.123 + 0.745j # douady's rabbit fractal
//# c = -0.750 + 0j # san marco fractal
//# c = -0.391 - 0.587j # siegel disk fractal
//# c = -0.7 - 0.3j # NEAT cauliflower thingy
//# c = -0.75 - 0.2j # galaxies
//# c = -0.75 + 0.15j # groovy
//# c = -0.7 + 0.35j # frost
int main(void) {
char fileName[32];
size_t dataSize= sizeof(int)*DIMX*DIMY;
int* hdata;
CUDA_CHECK_RETURN(cudaMallocHost(&hdata,sizeof(int)*DIMX*DIMY));
int* ddata;
CUDA_CHECK_RETURN(cudaMalloc(&ddata,sizeof(int)*DIMX*DIMY));
dim3 bs(BLOCKX, BLOCKY);
dim3 gs(divup(DIMX, bs.x), divup(DIMY, bs.y));
float incre =0.00000003;
float inci =-0.00009;
float startre=-0.75;
float starti=0.09;
float zoom=1.0;
for(int i=0;i<2500;i++){
//zoom+=0.001;
printf("julia_%f_%f.png\n", startre+i*incre,starti+i*inci);
c0 = cMakecuComplex(startre+i*incre,starti+i*inci);
// sprintf(fileName, "julia_%f_%f.png", startre+i*incre,starti+i*inci);
sprintf(fileName, "julia_%05d.png", i);
computeJulia<<<gs,bs>>>(ddata,c0,zoom);
CUDA_CHECK_RETURN(cudaMemcpy(hdata, ddata, dataSize, cudaMemcpyDeviceToHost));
save_image(fileName,hdata,DIMX,DIMY);
}
// for(int i=0;i<DIMX*DIMY;i++){
// printf("%i, ",hdata[i]);
// }
cudaFreeHost(hdata);
cudaFree(ddata);
return 0;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment