From 6878d573e50cca6ad917b67d596e5bb38c93b42a Mon Sep 17 00:00:00 2001 From: Tom Smeding Date: Tue, 5 May 2020 16:47:42 +0200 Subject: Subsampling for faster image --- Makefile | 2 +- mandel.cu | 102 +++++++++++++++++++++++---- mandel.h | 2 +- sdl.cpp | 238 +++++++++++++++++++++++++++++++++++--------------------------- 4 files changed, 222 insertions(+), 122 deletions(-) diff --git a/Makefile b/Makefile index 2b12c32..f5eb1b4 100644 --- a/Makefile +++ b/Makefile @@ -1,7 +1,7 @@ CXX = g++ NVCC = nvcc CUFLAGS = -O3 -CXXFLAGS = -Wall -Wextra -std=c++11 -O3 +CXXFLAGS = -Wall -Wextra -std=c++17 -O3 LDFLAGS = -L/opt/cuda/lib64 -lcudart CXXFLAGS += $(shell pkg-config --cflags sdl2) diff --git a/mandel.cu b/mandel.cu index 480daed..386f962 100644 --- a/mandel.cu +++ b/mandel.cu @@ -60,34 +60,77 @@ void mandel_free(Mandel *ctx) { delete ctx; } +// These macros may assume that the parameters are simple variable names. +#if 1 +#define MANDEL_FUNCTION_REAL(a, b, a2, b2, x, y) (a2 - b2 + x) +#define MANDEL_FUNCTION_IMAG(a, b, a2, b2, x, y) (2 * a * b + y) +#endif +#if 0 +#define MANDEL_FUNCTION_REAL(a, b, a2, b2, x, y) (a2 - b2 + x + 1) +#define MANDEL_FUNCTION_IMAG(a, b, a2, b2, x, y) (2 * a * b + y) +#endif +#if 0 +#define MANDEL_FUNCTION_REAL(a, b, a2, b2, x, y) (a2 * a - 3 * a * b * b + x) +#define MANDEL_FUNCTION_IMAG(a, b, a2, b2, x, y) (3 * a * a * b - b2 * b + y) +#endif + #define MANDEL_GENERIC(dst, ctx, par, ix, iy, idx) { \ const double x = (par).cx - (par).imgw / 2 + (par).imgw * (ix) / ((ctx).w-1); \ const double y = (par).cy - (ctx).imgh / 2 + (ctx).imgh * (iy) / ((ctx).h-1); \ double a = x, b = y, a2 = a * a, b2 = b * b; \ int16_t iter, maxiter = (par).maxit; \ for (iter = 0; iter < maxiter && a2 + b2 < 4; iter++) { \ - b = 2 * a * b + y; a = a2 - b2 + x; \ + double newa = MANDEL_FUNCTION_REAL(a, b, a2, b2, x, y); \ + double newb = MANDEL_FUNCTION_IMAG(a, b, a2, b2, x, y); \ + a = newa; b = newb; \ a2 = a * a; b2 = b * b; \ } \ (dst)[idx] = iter; \ } -__global__ void mandel_gpu(int16_t *dst, const Mandel *ctx, const Params *par) { - const int idx = blockDim.x * blockIdx.x + threadIdx.x; - const int ix = idx % (int)ctx->w, iy = ctx->h - 1 - idx / (int)ctx->w; +__global__ void mandel_gpu(int16_t *dst, const Mandel *ctx, const Params *par, int subsample) { + const int subidx = blockDim.x * blockIdx.x + threadIdx.x; + const int subw = (ctx->w + subsample - 1) / subsample; + const int subix = subidx % subw, subiy = subidx / subw; + const int ix = subix * subsample, iy = subiy * subsample; + const int idx = ctx->w * iy + ix; + const int riy = ctx->h - 1 - iy; + if (iy >= ctx->h) return; - MANDEL_GENERIC(dst, *ctx, *par, ix, iy, idx); + MANDEL_GENERIC(dst, *ctx, *par, ix, riy, idx); +} + +__global__ void unsubsample_gpu(int16_t *dst, const Mandel *ctx, int subsample) { + const int idx = blockDim.x * blockIdx.x + threadIdx.x; + const int ix = idx % ctx->w, iy = idx / ctx->w; + const int cellix = ix / subsample * subsample, celliy = iy / subsample * subsample; + + dst[ctx->w * iy + ix] = dst[ctx->w * celliy + cellix]; } // Unused in GPU mode __attribute__((unused)) -static inline void mandel_cpu(int16_t *dst, const Mandel *ctx, const Params *par) { - int idx = 0; - for (int iy = 0; iy < ctx->h; iy++) { - for (int ix = 0; ix < ctx->w; ix++) { - MANDEL_GENERIC(dst, *ctx, *par, ix, iy, idx); - idx++; +static inline void mandel_cpu(int16_t *dst, const Mandel *ctx, const Params *par, int subsample) { + for (int iy = 0; iy < ctx->h; iy += subsample) { + for (int ix = 0; ix < ctx->w; ix += subsample) { + int idx = ctx->w * iy + ix; + int riy = ctx->h - 1 - iy; + MANDEL_GENERIC(dst, *ctx, *par, ix, riy, idx); + } + } + + if (subsample == 1) return; + + for (int iy = 0; iy < ctx->h; iy += subsample) { + for (int ix = 0; ix < ctx->w; ix += subsample) { + const int16_t value = dst[ctx->w * iy + ix]; + + for (int subiy = 0; subiy < subsample; subiy++) { + for (int subix = 0; subix < subsample; subix++) { + dst[ctx->w * (iy + subiy) + ix + subix] = value; + } + } } } } @@ -126,20 +169,28 @@ double mandel_imgh(const Mandel *ctx, const Params *par) { return par->imgw * ctx->h / ctx->w; } -void mandel_render(uint8_t *dst, Mandel *ctx, const Params *par) { +void mandel_render(uint8_t *dst, Mandel *ctx, const Params *par, size_t subsample) { ctx->imgh = mandel_imgh(ctx, par); int64_t t1 = gettimestamp(); #ifdef USE_GPU + assert((size_t)(int)subsample == subsample); + + const int subfactor = subsample * subsample; + CUDA_CHECK(cudaMemcpy(ctx->devPar, par, sizeof(Params), cudaMemcpyHostToDevice)); CUDA_CHECK(cudaMemcpy(ctx->devCtx, ctx, sizeof(Mandel), cudaMemcpyHostToDevice)); const int nblocks = (ctx->w * ctx->h + 1023) / 1024; - mandel_gpu<<>>(ctx->devImg, ctx->devCtx, ctx->devPar); + mandel_gpu<<>>(ctx->devImg, ctx->devCtx, ctx->devPar, subsample); + + if (subsample > 1) { + unsubsample_gpu<<>>(ctx->devImg, ctx->devCtx, subsample); + } CUDA_CHECK(cudaMemcpy(ctx->img, ctx->devImg, ctx->w * ctx->h * sizeof(int16_t), cudaMemcpyDeviceToHost)); #else - mandel_cpu(ctx->img, ctx, par); + mandel_cpu(ctx->img, ctx, par, subsample); #endif int64_t t2 = gettimestamp(); @@ -155,20 +206,41 @@ void mandel_render(uint8_t *dst, Mandel *ctx, const Params *par) { int64_t t3 = gettimestamp(); cout << "gpu part: " << (t2 - t1) / 1000000.0 << " sec " - << "cpu part: " << (t3 - t2) / 1000000 << " sec" << endl; + << "cpu part: " << (t3 - t2) / 1000000.0 << " sec" << endl; } #if 0 int main() { Mandel *ctx = mandel_init(1920, 1080); Params par = mandel_default_params(); + +#if 0 par.cx = -0.73844331961137488; par.cy = -0.20921562151741355; par.imgw = 9.6624448855068765e-08; par.maxit = 20000; + uint8_t *img = new uint8_t[3 * 1920 * 1080]; mandel_render(img, ctx, &par); mandel_free(ctx); bmp_rgb_encode_file("out.bmp", img, 1920, 1080); +#else + par.cx = -0.23845305789504126; + par.cy = 0.87200567648858607; + par.imgw = 3.5; // zoom to 7.5815186948616694e-14 + par.maxit = 6144; + + uint8_t *img = new uint8_t[3 * 1920 * 1080]; + + for (int i = 0; i <= 141; i++) { + mandel_render(img, ctx, &par, 1); + char fname[64]; + sprintf(fname, "outdir/out%03d.bmp", i); + bmp_rgb_encode_file(fname, img, 1920, 1080); + par.imgw *= 0.8; + } + + mandel_free(ctx); +#endif } #endif diff --git a/mandel.h b/mandel.h index aed1738..99d4b60 100644 --- a/mandel.h +++ b/mandel.h @@ -15,4 +15,4 @@ void mandel_free(Mandel *ctx); double mandel_imgh(const Mandel *ctx, const Params *par); // dst[3*w*h] -void mandel_render(uint8_t *dst, Mandel *ctx, const Params *par); +void mandel_render(uint8_t *dst, Mandel *ctx, const Params *par, size_t subsample); diff --git a/sdl.cpp b/sdl.cpp index dcc2bb9..ea2f0f9 100644 --- a/sdl.cpp +++ b/sdl.cpp @@ -1,6 +1,7 @@ #include #include #include +#include #include #include "mandel.h" @@ -8,111 +9,119 @@ using namespace std; #if 1 -int main() { - if (SDL_Init(SDL_INIT_VIDEO) < 0) { - cerr << "SDL init error: " << SDL_GetError() << endl; - return 1; - } - - SDL_Window *window = SDL_CreateWindow("Mandel", SDL_WINDOWPOS_UNDEFINED, SDL_WINDOWPOS_UNDEFINED, 960, 640, SDL_WINDOW_SHOWN); - if (!window) { - cerr << "SDL window creation error: " << SDL_GetError() << endl; - return 1; - } +struct State { + SDL_Window *window = nullptr; - int w, h; - SDL_GetWindowSize(window, &w, &h); - Mandel *ctx = mandel_init(w, h); + int w = -1, h = -1; + Mandel *ctx = nullptr; Params par = mandel_default_params(); - uint8_t *imgdata = new uint8_t[3 * w * h]; - mandel_render(imgdata, ctx, &par); + uint8_t *imgdata = nullptr; - SDL_Surface *mandelsurf = SDL_CreateRGBSurfaceFrom(imgdata, w, h, 24, 3 * w, 0x0000000ff, 0x0000ff00, 0x00ff0000, 0); + SDL_Surface *mandelsurf = nullptr; - SDL_Surface *winsurf = SDL_GetWindowSurface(window); + size_t current_resolution = 4; - SDL_BlitSurface(mandelsurf, nullptr, winsurf, nullptr); + State() { + if (SDL_Init(SDL_INIT_VIDEO) < 0) { + cerr << "SDL init error: " << SDL_GetError() << endl; + exit(1); + } - SDL_UpdateWindowSurface(window); + window = SDL_CreateWindow("Mandel", SDL_WINDOWPOS_UNDEFINED, SDL_WINDOWPOS_UNDEFINED, 960, 640, SDL_WINDOW_SHOWN); + if (!window) { + cerr << "SDL window creation error: " << SDL_GetError() << endl; + exit(1); + } - while (true) { - const double min_img_size = min(par.imgw, mandel_imgh(ctx, &par)); - const double scroll_dist = min_img_size * 0.1; + handle_resize(); + mandel_render(imgdata, ctx, &par, current_resolution); - bool need_redraw = false; - - { - stringstream ss; - ss << "Mandel " - << setprecision(17) - << " x=" << par.cx - << " y=" << par.cy - << " w=" << par.imgw - << " i=" << par.maxit; - SDL_SetWindowTitle(window, ss.str().data()); - } + blit(); + SDL_UpdateWindowSurface(window); + } + void blit() { + SDL_Surface *winsurf = SDL_GetWindowSurface(window); + if (SDL_MUSTLOCK(winsurf)) SDL_LockSurface(winsurf); + SDL_BlitSurface(mandelsurf, nullptr, winsurf, nullptr); + if (SDL_MUSTLOCK(winsurf)) SDL_UnlockSurface(winsurf); + } + + void update_title() { + stringstream ss; + ss << "Mandel " + << setprecision(17) + << " x=" << par.cx + << " y=" << par.cy + << " w=" << par.imgw + << " i=" << par.maxit; + SDL_SetWindowTitle(window, ss.str().data()); + } + + SDL_Event wait_event() const { SDL_Event e; if (SDL_WaitEvent(&e) == 0) { cerr << "SDL event wait error: " << SDL_GetError() << endl; - break; + exit(1); } + return e; + } - winsurf = SDL_GetWindowSurface(window); + optional poll_event() const { + SDL_Event e; + if (SDL_PollEvent(&e) == 0) { + return nullopt; + } + return {e}; + } + + void handle_resize() { + if (ctx) mandel_free(ctx); + if (imgdata) delete[] imgdata; + if (mandelsurf) SDL_FreeSurface(mandelsurf); + + SDL_GetWindowSize(window, &w, &h); + cerr << "w=" << w << " h=" << h << endl; + imgdata = new uint8_t[3 * w * h]; + mandelsurf = SDL_CreateRGBSurfaceFrom(imgdata, w, h, 24, 3 * w, 255, 255 << 8, 255 << 16, 0); + ctx = mandel_init(w, h); + } + + enum class action { + quit, + redraw, + idle, + }; + + action handle_event(const SDL_Event &e) { + const double min_img_size = min(par.imgw, mandel_imgh(ctx, &par)); + const double scroll_dist = min_img_size * 0.1; switch (e.type) { case SDL_QUIT: - break; + return action::quit; case SDL_KEYDOWN: switch (e.key.keysym.sym) { - case SDLK_UP: - par.cy += scroll_dist; - need_redraw = true; - break; - - case SDLK_DOWN: - par.cy -= scroll_dist; - need_redraw = true; - break; + case SDLK_UP: par.cy += scroll_dist; return action::redraw; + case SDLK_DOWN: par.cy -= scroll_dist; return action::redraw; + case SDLK_RIGHT: par.cx += scroll_dist; return action::redraw; + case SDLK_LEFT: par.cx -= scroll_dist; return action::redraw; - case SDLK_RIGHT: - par.cx += scroll_dist; - need_redraw = true; - break; + case SDLK_KP_PLUS: case SDLK_PLUS: case SDLK_EQUALS: + par.imgw *= 0.8; return action::redraw; - case SDLK_LEFT: - par.cx -= scroll_dist; - need_redraw = true; - break; - - case SDLK_KP_PLUS: - case SDLK_PLUS: - case SDLK_EQUALS: - par.imgw *= 0.8; - need_redraw = true; - break; - - case SDLK_KP_MINUS: - case SDLK_MINUS: - par.imgw *= 1.25; - need_redraw = true; - break; + case SDLK_KP_MINUS: case SDLK_MINUS: + par.imgw *= 1.25; return action::redraw; case SDLK_PAGEUP: - par.maxit += 128; - need_redraw = true; - break; - + par.maxit += 128; return action::redraw; case SDLK_PAGEDOWN: - if (par.maxit > 256) { - par.maxit -= 128; - need_redraw = true; - } + if (par.maxit > 256) { par.maxit -= 128; return action::redraw; } break; case SDLK_q: - goto cleanup; + return action::quit; } break; @@ -120,43 +129,62 @@ int main() { switch (e.window.event) { case SDL_WINDOWEVENT_FOCUS_GAINED: SDL_UpdateWindowSurface(window); - break; + return action::idle; case SDL_WINDOWEVENT_SIZE_CHANGED: - mandel_free(ctx); - delete[] imgdata; - SDL_FreeSurface(mandelsurf); - - winsurf = SDL_GetWindowSurface(window); - SDL_GetWindowSize(window, &w, &h); - cerr << "w=" << w << " h=" << h << endl; - cerr << "surface: w=" << winsurf->w << " h=" << winsurf->h << endl; - imgdata = new uint8_t[3 * w * h]; - mandelsurf = SDL_CreateRGBSurfaceFrom(imgdata, w, h, 24, 3 * w, 0x0000000ff, 0x0000ff00, 0x00ff0000, 0); - ctx = mandel_init(w, h); - need_redraw = true; - break; + handle_resize(); + return action::redraw; } break; } - if (need_redraw) { - SDL_GetWindowSize(window, &w, &h); - cerr << "w=" << w << " h=" << h << endl; - cerr << "winsurf: w=" << winsurf->w << " h=" << winsurf->h << endl; - cerr << "mandelsurf: w=" << mandelsurf->w << " h=" << mandelsurf->h << endl; - mandel_render(imgdata, ctx, &par); + return action::idle; + } +}; - if (SDL_MUSTLOCK(winsurf)) SDL_LockSurface(winsurf); - SDL_BlitSurface(mandelsurf, nullptr, winsurf, nullptr); - if (SDL_MUSTLOCK(winsurf)) SDL_UnlockSurface(winsurf); +int main() { + State state; - SDL_UpdateWindowSurface(window); + while (true) { + state.update_title(); + + State::action action; + if (state.current_resolution == 1) { + SDL_Event ev = state.wait_event(); + action = state.handle_event(ev); + } else { + if (optional optev = state.poll_event()) { + action = state.handle_event(*optev); + } else { + action = State::action::idle; + } } - } -cleanup: - SDL_DestroyWindow(window); - SDL_Quit(); + switch (action) { + case State::action::quit: + return 0; + + case State::action::redraw: + state.current_resolution = 4; + + SDL_GetWindowSize(state.window, &state.w, &state.h); + cerr << "w=" << state.w << " h=" << state.h << endl; + // cerr << "mandelsurf: w=" << state.mandelsurf->w << " h=" << state.mandelsurf->h << endl; + mandel_render(state.imgdata, state.ctx, &state.par, state.current_resolution); + + state.blit(); + SDL_UpdateWindowSurface(state.window); + break; + + case State::action::idle: + if (state.current_resolution > 1) { + state.current_resolution /= 2; + mandel_render(state.imgdata, state.ctx, &state.par, state.current_resolution); + } + state.blit(); + SDL_UpdateWindowSurface(state.window); + break; + } + } } #endif -- cgit v1.2.3