218 lines
6.3 KiB
C
218 lines
6.3 KiB
C
/*
|
|
* Josh Holtrop
|
|
* 2007-11-28
|
|
* Newton Fractal Renderer using MPI and SDL
|
|
*/
|
|
|
|
#include <stdio.h>
|
|
#include <math.h>
|
|
#include <SDL/SDL.h>
|
|
#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;
|
|
}
|