1 /*
2
3 LaplaceRelaxation.c : Laplacian relaxation demo
4
5 Demo Program that
6 - loads an image
7 - converts it to a brightness map
8 - binarize brightness map
9 - performs a Laplace relaxation on it
10 - displays the solution using a sinusoidally modulated lookup table
11
12 Very slow! :-) TODO: use imagefilter MMX routines
13
14 (C) A. Schiffler, Sep 2010, zlib License
15
16 */
17
18 #include <stdlib.h>
19 #include <stdio.h>
20 #include <string.h>
21
22 #include "SDL.h"
23
24 #ifdef WIN32
25 #include <windows.h>
26 #include "SDL_gfxPrimitives.h"
27 #include "SDL_rotozoom.h"
28 #include "SDL_imageFilter.h"
29 #ifndef bcmp
30 #define bcmp(s1, s2, n) memcmp ((s1), (s2), (n))
31 #endif
32 #else
33 #include "SDL/SDL_gfxPrimitives.h"
34 #include "SDL/SDL_rotozoom.h"
35 #include "SDL/SDL_imageFilter.h"
36 #endif
37
38 int iterations = 400;
39 int contours = 20;
40 int threshold = 4;
41 int w = 640, h = 480;
42
HandleEvent()43 void HandleEvent()
44 {
45 SDL_Event event;
46
47 /* Check for events */
48 SDL_PollEvent(&event);
49 switch (event.type) {
50 case SDL_KEYDOWN:
51 case SDL_QUIT:
52 exit(0);
53 break;
54 }
55 }
56
ClearScreen(SDL_Surface * screen)57 void ClearScreen(SDL_Surface *screen)
58 {
59 int i;
60 /* Set the screen to black */
61 if ( SDL_LockSurface(screen) == 0 ) {
62 Uint8 *pixels;
63 pixels = (Uint8 *)screen->pixels;
64 for ( i=0; i<screen->h; ++i ) {
65 memset(pixels, 0,
66 screen->w*screen->format->BytesPerPixel);
67 pixels += screen->pitch;
68 }
69 SDL_UnlockSurface(screen);
70 }
71 }
72
73 /*
74 Returns a potential value from either the static or the dynamic charge map based
75 on the binarized map at coordinates x,y for a grid of size w,h clamping to 0 at the edges.
76 */
GetPotential(Uint8 * map,Uint8 * greyscale,double * relaxed,int x,int y,int w,int h)77 double GetPotential(Uint8 *map, Uint8 *greyscale, double *relaxed, int x, int y, int w, int h)
78 {
79 int offset;
80 double v;
81
82 if ((x<0) || (y<0) || (x>(w-1)) || (y>(h-1))) {
83 v = 0.0;
84 } else {
85 offset = x+w*y;
86 if ((Uint8)map[offset]!=0) {
87 v = (double)greyscale[offset];
88 } else {
89 v = relaxed[offset];
90 }
91 }
92
93 return v;
94 }
95
96 /*
97 * Return the pixel value at (x, y)
98 */
getPixel(SDL_Surface * surface,int x,int y)99 Uint32 getPixel(SDL_Surface *surface, int x, int y)
100 {
101 int bpp = surface->format->BytesPerPixel;
102
103 /* Here p is the address to the pixel we want to retrieve */
104 Uint8 *p = (Uint8 *)surface->pixels + y * surface->pitch + x * bpp;
105
106 switch (bpp) {
107 case 1:
108 return *p;
109
110 case 2:
111 return *(Uint16 *)p;
112
113 case 3:
114 if (SDL_BYTEORDER == SDL_BIG_ENDIAN)
115 return p[0] << 16 | p[1] << 8 | p[2];
116 else
117 return p[0] | p[1] << 8 | p[2] << 16;
118
119 case 4:
120 return *(Uint32 *)p;
121
122 default:
123 return 0; /* shouldn't happen, but avoids warnings */
124 } // switch
125 }
126
Draw(SDL_Surface * screen,char * bmpfile)127 void Draw (SDL_Surface *screen, char *bmpfile)
128 {
129 char messageText[128];
130 Uint32 rmask, gmask, bmask, amask;
131 SDL_Surface *picture;
132 SDL_Surface *mapped_picture;
133 SDL_Surface *rotozoom_picture;
134 SDL_PixelFormat *pixelFormat;
135 Uint8 *grayscale, *map, *curmap;
136 double *unrelaxed, *relaxed, *currelax;
137 int mapsize, relaxsize;
138 int rowskip;
139 char *pixel;
140 Uint32 p;
141 int x, y;
142 Uint8 r, g, b, a;
143 double dy;
144 double t, s;
145 int i;
146 double c1, c2, c3, c4, ca;
147 Uint8 lookupTable[256];
148 double zw, zh, zf;
149
150 /* SDL interprets each pixel as a 32-bit number, so our masks must depend
151 on the endianness (byte order) of the machine */
152 #if SDL_BYTEORDER == SDL_BIG_ENDIAN
153 rmask = 0xff000000;
154 gmask = 0x00ff0000;
155 bmask = 0x0000ff00;
156 amask = 0x000000ff;
157 #else
158 rmask = 0x000000ff;
159 gmask = 0x0000ff00;
160 bmask = 0x00ff0000;
161 amask = 0xff000000;
162 #endif
163
164 /* Load the image into a surface */
165 printf("Loading picture: %s\n", bmpfile);
166 picture = SDL_LoadBMP(bmpfile);
167 if ( picture == NULL ) {
168 fprintf(stderr, "Couldn't load %s: %s\n", bmpfile, SDL_GetError());
169 return;
170 }
171
172 /* Convert image to a brightness map */
173 printf("Allocating workbuffers\n");
174 mapsize = picture->w * picture->h;
175 relaxsize = mapsize * sizeof(double);
176 grayscale = (Uint8 *)malloc(mapsize);
177 map = (Uint8 *)malloc(mapsize);
178 unrelaxed = (double *)malloc(relaxsize);
179 relaxed = (double *)malloc(relaxsize);
180 if ((grayscale == NULL) || (map == NULL) || (unrelaxed == NULL) || (relaxed == NULL))
181 {
182 fprintf(stderr, "Out of memory\n");
183 return;
184 }
185 memset(grayscale, 0, mapsize);
186 memset(map, 0, mapsize);
187 memset(unrelaxed, 0, relaxsize);
188 memset(relaxed, 0, relaxsize);
189
190 printf("Converting image to brightness map\n");
191 pixel = picture->pixels;
192 pixelFormat = picture->format;
193 rowskip = (picture->pitch - picture->w * picture->format->BytesPerPixel);
194 curmap = grayscale;
195 for (y=0; y < picture->h; y++) {
196 for (x=0; x < picture->w; x++) {
197 // Get RGB
198 p = getPixel(picture, x, y);
199 SDL_GetRGBA(p, pixelFormat, &r, &g, &b, &a);
200
201 // Calculate luminance (Y = 0.3R + 0.59G + 0.11B;)
202 dy = (0.30 * r) + (0.59 * g) + (0.11 * b);
203 if (dy<0.0) {
204 dy=0.0;
205 } else if (dy>255.0) {
206 dy=255.0;
207 }
208 *curmap = (Uint8)dy;
209
210 // Next pixel
211 pixel += picture->format->BytesPerPixel;
212 curmap++;
213 }
214
215 // Next row
216 pixel += rowskip;
217 }
218
219 /* --- Prepare relaxation loop --- */
220
221 /* Binarize luminance map */
222 SDL_imageFilterBinarizeUsingThreshold(grayscale, map, mapsize, threshold);
223
224 /* Create cos^5 based lookup table */
225 t = 0.0;
226 for (i = 0; i < 256; i++)
227 {
228 s = 1.0 - 0.5 * (1.0 + cos(t));
229 s = 255.0 * s * s * s * s * s;
230 lookupTable[i] = (Uint8)s;
231 t += ((double)contours*2.0*3.141592653589/128.0);
232 }
233
234 /* Create new surface for relaxed image */
235 mapped_picture = SDL_CreateRGBSurface(SDL_SWSURFACE, picture->w, picture->h, 32,
236 rmask, gmask, bmask, amask);
237 if (mapped_picture == NULL) {
238 fprintf(stderr, "CreateRGBSurface failed: %s\n", SDL_GetError());
239 return;
240 }
241
242 /* Apply relaxation algorithm */
243 printf("Applying Laplacian relaxation: %i iterations\n", iterations);
244 // Iterate to relax
245 for (i = 0; i <= iterations; i++)
246 {
247 // Backup original
248 memcpy(unrelaxed, relaxed, relaxsize);
249 // Average into relaxed
250 for (x=0; x < picture->w; x++) {
251 for (y=0; y < picture->h; y++) {
252 // up
253 c1 = GetPotential(map, grayscale, unrelaxed, x, y-1, picture->w, picture->h);
254 // down
255 c2 = GetPotential(map, grayscale, unrelaxed, x, y+1, picture->w, picture->h);
256 // left
257 c3 = GetPotential(map, grayscale, unrelaxed, x-1, y, picture->w, picture->h);
258 // right
259 c4 = GetPotential(map, grayscale, unrelaxed, x+1, y, picture->w, picture->h);
260 // average and store
261 ca = ((c1 + c2 + c3 + c4)/4.0);
262 relaxed[x+y*picture->w] = ca;
263 }
264 }
265
266 // Draw only sometimes
267 if (((i % 10)==0) || (i==iterations)) {
268
269 /* --- Create image with contour map --- */
270
271 /* Fill output surface via colormap */
272 currelax = relaxed;
273 for (y=0; y<mapped_picture->h; y++) {
274 for (x=0; x<mapped_picture->w; x++) {
275 if (map[x+y*picture->w]!=0) {
276 r = g = b = grayscale[x+y*picture->w];
277 } else {
278 r = g = b = lookupTable[(Uint8)*currelax];
279 }
280 pixelRGBA(mapped_picture, x, y, r, g, b, 255);
281 currelax++;
282 }
283 }
284
285 /* --- Scale and draw to screen --- */
286
287 /* Scale to screen size */
288 zw = (double)screen->w/(double)mapped_picture->w;
289 zh = (double)screen->h/(double)mapped_picture->h;
290 zf = (zw < zh) ? zw : zh;
291 if ((rotozoom_picture=zoomSurface(mapped_picture, zf, zf, 1))==NULL) {
292 fprintf(stderr, "Rotozoom failed: %s\n", SDL_GetError());
293 return;
294 }
295
296 /* Draw surface to screen */
297 if ( SDL_BlitSurface(rotozoom_picture, NULL, screen, NULL) < 0 ) {
298 fprintf(stderr, "Blit failed: %s\n", SDL_GetError());
299 return;
300 }
301 SDL_FreeSurface(rotozoom_picture);
302
303 /* Info */
304 if (i != iterations) {
305 sprintf(messageText,"%i", i);
306 stringRGBA(screen, 8, 8, messageText, 255, 255, 255, 255);
307 }
308
309 /* Display by flipping screens */
310 SDL_Flip(screen);
311 }
312
313 /* Maybe quit */
314 HandleEvent();
315 }
316
317 /* Save final picture */
318 if (SDL_SaveBMP(mapped_picture, "result.bmp") <0) {
319 fprintf(stderr, "Save BMP failed: %s\n", SDL_GetError());
320 return;
321 }
322 free(map);
323 free(grayscale);
324 free(unrelaxed);
325 free(relaxed);
326 SDL_FreeSurface(picture);
327 SDL_FreeSurface(mapped_picture);
328
329 return;
330 }
331
main(int argc,char * argv[])332 int main ( int argc, char *argv[] )
333 {
334 SDL_Surface *screen;
335 int desired_bpp;
336 Uint32 video_flags;
337 char *bmpfile = NULL;
338
339 /* Title */
340 fprintf(stderr,"Laplace relaxation demo\n");
341
342 /* Set default options and check command-line */
343 desired_bpp = 32;
344 video_flags = SDL_SWSURFACE | SDL_SRCALPHA | SDL_RESIZABLE | SDL_DOUBLEBUF;
345 while ( argc > 1 ) {
346 if ( strcmp(argv[1], "-iterations") == 0 ) {
347 if ( argv[2] && ((iterations = atoi(argv[2])) > 0) ) {
348 argv += 2;
349 argc -= 2;
350 } else {
351 fprintf(stderr,
352 "The -iterations option requires an argument\n");
353 exit(1);
354 }
355 } else if ( strcmp(argv[1], "-contours") == 0 ) {
356 if ( argv[2] && ((contours = atoi(argv[2])) > 0) ) {
357 argv += 2;
358 argc -= 2;
359 } else {
360 fprintf(stderr,
361 "The -contours option requires an argument\n");
362 exit(1);
363 }
364 } else if ( strcmp(argv[1], "-threshold") == 0 ) {
365 if ( argv[2] && ((threshold = atoi(argv[2])) > 0) ) {
366 argv += 2;
367 argc -= 2;
368 } else {
369 fprintf(stderr,
370 "The -threshold option requires an argument\n");
371 exit(1);
372 }
373 } else if (( strcmp(argv[1], "-help") == 0 ) || (strcmp(argv[1], "--help") == 0)) {
374 printf ("Usage:\n%s [options] filename\n");
375 printf ("Options:\n");
376 printf (" -iterations # Number of relaxation iterations to perform, default: 400\n");
377 printf (" -contours # Number of contours over range, default: 8\n");
378 printf (" -threshold # Binarization threshold for Y, default: 20\n");
379 exit(0);
380 } else {
381 bmpfile = argv[1];
382 break;
383 }
384 }
385
386 /* Initialize SDL */
387 if ( SDL_Init(SDL_INIT_VIDEO) < 0 ) {
388 fprintf(stderr,
389 "Couldn't initialize SDL: %s\n", SDL_GetError());
390 exit(1);
391 }
392 atexit(SDL_Quit); /* Clean up on exit */
393
394 /* Initialize the display */
395 screen = SDL_SetVideoMode(w, h, desired_bpp, video_flags);
396 if ( screen == NULL ) {
397 fprintf(stderr, "Couldn't set %dx%dx%d video mode: %s\n",
398 w, h, desired_bpp, SDL_GetError());
399 exit(1);
400 }
401
402 /* Show some info */
403 printf("Set %dx%dx%d mode\n",
404 screen->w, screen->h, screen->format->BitsPerPixel);
405 printf("Video surface located in %s memory.\n",
406 (screen->flags & SDL_HWSURFACE) ? "video" : "system");
407
408 /* Check for double buffering */
409 if ( screen->flags & SDL_DOUBLEBUF ) {
410 printf("Double-buffering enabled - good!\n");
411 }
412
413 /* Set the window manager title bar */
414 SDL_WM_SetCaption("Laplace relaxation demo", "laplacerelaxation");
415
416 /* Do all the drawing work */
417 Draw(screen, bmpfile);
418
419 /* Wait for keypress */
420 while(1) {
421 HandleEvent();
422 SDL_Delay(100);
423 }
424
425 return(0);
426 }
427