/* * Josh Holtrop * 2007-11-28 * Newton Fractal Renderer using MPI and SDL */ #include #include #include #include "complex.h" #define LEEWAY 0.001 #define MAX_ITERATIONS 32 #define DEFAULT_WIN_WIDTH 800 #define DEFAULT_WIN_HEIGHT 600 #define PROGNAME "Josh's MPI/pthread newton fractal renderer" /* Prototypes */ SDL_Surface * sdl_init(unsigned int width, unsigned int height); void calculateRow(int winWidth, int winHeight, int row, double viewWidth, double viewHeight, double x_center, double y_center, Uint32 * rowVals); Uint32 computePoint(double x, double y); void draw(SDL_Surface * screen, int winWidth, int winHeight, double x_center, double y_center, double viewWidth, double viewHeight); void main_loop(SDL_Surface * screen, int winWidth, int winHeight, double x_center, double y_center, double viewWidth, double viewHeight); int main(int argc, char * argv[]) { SDL_Surface * screen; int winWidth = DEFAULT_WIN_WIDTH; int winHeight = DEFAULT_WIN_HEIGHT; double viewWidth = 2.0, x_center = 0.0, y_center = 0.0; double viewHeight = ((double)winHeight/(double)winWidth) * viewWidth; if ( !(screen = sdl_init(winWidth, winHeight)) ) { fprintf(stderr, "SDL initialization error!\n"); return -1; } main_loop(screen, winWidth, winHeight, x_center, y_center, viewWidth, viewHeight); return 0; } void main_loop(SDL_Surface * screen, int winWidth, int winHeight, double x_center, double y_center, double viewWidth, double viewHeight) { SDL_Event event; for (;;) { draw(screen, winWidth, winHeight, x_center, y_center, viewWidth, viewHeight); while(SDL_WaitEvent(&event)) if (event.type == SDL_QUIT || event.type == SDL_MOUSEBUTTONDOWN) break; switch (event.type) { case SDL_QUIT: return; case SDL_MOUSEBUTTONDOWN: { int button = event.button.button; int x = event.motion.x; int y = event.motion.y; switch (button) { case 1: /* re-center on click point, zoom in */ x_center += (viewWidth / winWidth) * (x - (winWidth >> 1)); y_center += (viewHeight / winHeight) * ((winHeight >> 1) - y); viewWidth /= 2.0; viewHeight /= 2.0; break; case 2: /* rectangular selection to zoom into */ { int ox = x, oy = y; while (SDL_WaitEvent(&event)) if (event.type == SDL_MOUSEBUTTONUP) break; x = event.motion.x; y = event.motion.y; if (x < ox) { int t = ox; ox = x; x = ox; } if (y < oy) { int t = oy; oy = y; y = oy; } x_center += (viewWidth / winWidth) * (ox + ((x-ox) >> 1) - (winWidth >> 1)); y_center += (viewHeight / winHeight) * ((winHeight >> 1) - (oy + ((y-oy) >> 1))); viewWidth = (viewWidth / winWidth) * (x-ox+1); viewHeight = (viewHeight / winHeight) * (y-oy+1); } break; case 3: /* reset view */ viewWidth = 2.0; x_center = 0.0; y_center = 0.0; viewHeight = ((double)winHeight/(double)winWidth) * viewWidth; break; case 4: /* zoom in */ viewWidth /= 2.0; viewHeight /= 2.0; break; case 5: /* zoom out */ viewWidth *= 2.0; viewHeight *= 2.0; break; } } } } } void draw(SDL_Surface * screen, int winWidth, int winHeight, double x_center, double y_center, double viewWidth, double viewHeight) { Uint32 * pixels = (Uint32 *) screen->pixels; int ix, iy; for (iy = 0; iy < winHeight; iy++) { calculateRow(winWidth, winHeight, iy, viewWidth, viewHeight, x_center, y_center, pixels); pixels += winWidth; } SDL_Flip(screen); } SDL_Surface * sdl_init(unsigned int width, unsigned int height) { SDL_Surface * screen; if (SDL_Init(SDL_INIT_VIDEO)) { fprintf(stderr, "Failed to initialize SDL!\n"); return NULL; } atexit(SDL_Quit); if (!(screen = SDL_SetVideoMode(width, height, 32, SDL_DOUBLEBUF | SDL_HWSURFACE))) { fprintf(stderr, "Failed to set video mode!\n"); return NULL; } SDL_WM_SetCaption(PROGNAME, PROGNAME); return screen; } void calculateRow(int winWidth, int winHeight, int row, double viewWidth, double viewHeight, double x_center, double y_center, Uint32 * rowVals) { int i; double xspacing = viewWidth / winWidth; double yspacing = viewHeight / winHeight; double y = ((winHeight >> 1) - row) * yspacing + y_center; double x = x_center - (viewWidth / 2); for (i = 0; i < winWidth; i++, x += xspacing) *rowVals++ = computePoint(x, y); } Uint32 computePoint(double x, double y) { int n; complex_t rootGuess = {y, x}; complex_t t1, t2; for (n = 0; n < MAX_ITERATIONS; n++) { /* percentage of color to illuminate based on current iteration */ double pIlum = 1.0 - (n * 0.03); /* * These if statements check to see if the complex number (the root * guess) is within LEEWAY distance of a real root. If so, a unique * color is returned reflecting which root it is close to and how many * iterations it took for the root guess to get to that root */ if ((fabs(rootGuess.a - 1) < LEEWAY) && (fabs(rootGuess.b) < LEEWAY)) return ((Uint32)(0xFF * pIlum)) << 16; if ((fabs(rootGuess.a + 1) < LEEWAY) && (fabs(rootGuess.b) < LEEWAY)) return ((Uint32)(0xFF * pIlum) << 16) + ((Uint32)(0xFF * pIlum) << 8); if ((fabs(rootGuess.a - .5) < LEEWAY) && (fabs(rootGuess.b - sqrt(3)/2) < LEEWAY)) return (Uint32)(0xFF * pIlum) << 8; if ((fabs(rootGuess.a + .5) < LEEWAY) && (fabs(rootGuess.b - sqrt(3)/2) < LEEWAY)) return (Uint32)(0x88 * pIlum); if ((fabs(rootGuess.a - .5) < LEEWAY) && (fabs(rootGuess.b + sqrt(3)/2) < LEEWAY)) return (Uint32)(0xFF * pIlum); if ((fabs(rootGuess.a + .5) < LEEWAY) && (fabs(rootGuess.b + sqrt(3)/2) < LEEWAY)) return ((Uint32)(0xFF * pIlum) << 16) + (Uint32)(0xFF * pIlum); /* This expression evaluates the next complex number to be tested * for being close to a root. It uses Newton's method for finding * roots of equations according to the following recursive equation: * x_n+1 = x_n - y_n / dy_n * --> x_n+1 = x_n - (x_n^6 - 1) / 6x_n^5 * More information can be found at: * http://www.willamette.edu/~sekino/fractal/fractal.htm * http://www.math.iastate.edu/danwell/Fexplain/newt1.html * and Thomas' CALCULUS 10th edition by Finney, Weir, & Giordano pg. 302 */ complex_pow(&rootGuess, 6, &t1); complex_subs(&t1, 1, &t1); complex_pow(&rootGuess, 5, &t2); complex_muls(&t2, 6, &t2); complex_div(&t1, &t2, &t1); complex_sub(&rootGuess, &t1, &rootGuess); } return 0; }