gistfile1.cu

/* *
 * 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;
}