apps/plugins/fractals/mandelbrot_set.c

   1 /***************************************************************************
   2  *             __________               __   ___.
   3  *   Open      \______   \ ____   ____ |  | _\_ |__   _______  ___
   4  *   Source     |       _//  _ \_/ ___\|  |/ /| __ \ /  _ \  \/  /
   5  *   Jukebox    |    |   (  <_> )  \___|    < | \_\ (  <_> > <  <
   6  *   Firmware   |____|_  /\____/ \___  >__|_ \|___  /\____/__/\_ \
   7  *                     \/            \/     \/    \/            \/
   8  * $Id$
   9  *
  10  * Copyright (C) 2004 Matthias Wientapper
  11  * Heavily extended 2005 Jens Arnold
  12  *
  13  *
  14  * This program is free software; you can redistribute it and/or
  15  * modify it under the terms of the GNU General Public License
  16  * as published by the Free Software Foundation; either version 2
  17  * of the License, or (at your option) any later version.
  18  *
  19  * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY OF ANY
  20  * KIND, either express or implied.
  21  *
  22  ****************************************************************************/
  23 #include "fractal_sets.h"
  24 #include "mandelbrot_set.h"
  25
  26 #define BUTTON_YIELD_TIMEOUT (HZ / 4)
  27
  28 #ifdef USEGSLIB
  29 static unsigned char imgbuffer[LCD_HEIGHT];
  30 #else
  31 static fb_data imgbuffer[LCD_HEIGHT];
  32 #endif
  33
  34 #ifdef USEGSLIB
  35 #define LCOLOR(iter) ((iter ^ 7) << 5)
  36 #else
  37 /*
  38  * Spread iter's colors over color range.
  39  * 345 (=15*26-45) is max_iter maximal value
  40  * This implementation ignores pixel format, thus it is not uniformly spread
  41  */
  42 #define LCOLOR(iter) ((iter << LCD_DEPTH) / 345)
  43 #endif
  44
  45 #ifdef HAVE_LCD_COLOR
  46 #define COLOR(iter) (fb_data)LCOLOR(iter)
  47 #define CONVERGENCE_COLOR LCD_RGBPACK(0, 0, 0)
  48 #else /* greyscale */
  49 #define COLOR(iter) (unsigned char)LCOLOR(iter)
  50 #define CONVERGENCE_COLOR 0
  51 #endif
  52
  53 #if CONFIG_LCD == LCD_SSD1815
  54 /* Recorder, Ondio: pixel_height == 1.25 * pixel_width */
  55 #define MB_HEIGHT (LCD_HEIGHT*5/4)
  56 #else
  57 /* square pixels */
  58 #define MB_HEIGHT LCD_HEIGHT
  59 #endif
  60
  61 #define MB_XOFS (-0x03000000L)             /* -0.75 (s5.26) */
  62 #if (3000 * MB_HEIGHT / LCD_WIDTH) >= 2400 /* width is limiting factor */
  63 #define MB_XFAC (0x06000000LL)             /* 1.5 (s5.26) */
  64 #define MB_YFAC (MB_XFAC*MB_HEIGHT / LCD_WIDTH)
  65 #else                                      /* height is limiting factor */
  66 #define MB_YFAC (0x04cccccdLL)             /* 1.2 (s5.26) */
  67 #define MB_XFAC (MB_YFAC*LCD_WIDTH / MB_HEIGHT)
  68 #endif
  69
  70 #ifndef USEGSLIB
  71 #define UPDATE_FREQ (HZ/50)
  72 #endif
  73
  74 /* fixed point format s5.26: sign, 5 bits integer part, 26 bits fractional part */
  75 struct fractal_ops *ops;
  76 long x_min;
  77 long x_max;
  78 long x_step;
  79 long x_delta;
  80 long y_min;
  81 long y_max;
  82 long y_step;
  83 long y_delta;
  84 int step_log2;
  85 unsigned max_iter;
  86
  87 static void mandelbrot_init(void);
  88
  89 static int mandelbrot_calc_low_prec(struct fractal_rect *rect,
  90         int (*button_yield_cb)(void *), void *button_yield_ctx);
  91
  92 static int mandelbrot_calc_high_prec(struct fractal_rect *rect,
  93         int (*button_yield_cb)(void *), void *button_yield_ctx);
  94
  95 static void mandelbrot_move(int dx, int dy);
  96
  97 static int mandelbrot_zoom(int factor);
  98
  99 static int mandelbrot_precision(int d);
 100
 101 struct fractal_ops mandelbrot_ops =
 102 {
 103     .init = mandelbrot_init,
 104     .calc = NULL,
 105     .move = mandelbrot_move,
 106     .zoom = mandelbrot_zoom,
 107     .precision = mandelbrot_precision,
 108 };
 109
 110 #define LOG2_OUT_OF_BOUNDS -32767
 111
 112 static int ilog2_fp(long value) /* calculate integer log2(value_fp_6.26) */
 113 {
 114     int i = 0;
 115
 116     if (value <= 0)
 117     {
 118         return LOG2_OUT_OF_BOUNDS;
 119     }
 120     else if (value > (1L << 26))
 121     {
 122         while (value >= (2L << 26))
 123         {
 124             value >>= 1;
 125             i++;
 126         }
 127     }
 128     else
 129     {
 130         while (value < (1L << 26))
 131         {
 132             value <<= 1;
 133             i--;
 134         }
 135     }
 136     return i;
 137 }
 138
 139 static int recalc_parameters(void)
 140 {
 141     x_step = (x_max - x_min) / LCD_WIDTH;
 142     y_step = (y_max - y_min) / LCD_HEIGHT;
 143     step_log2 = ilog2_fp(MIN(x_step, y_step));
 144
 145     if (step_log2 == LOG2_OUT_OF_BOUNDS)
 146         return 1; /* out of bounds */
 147
 148     x_delta = X_DELTA(x_step);
 149     y_delta = Y_DELTA(y_step);
 150     max_iter = MAX(15, -15 * step_log2 - 45);
 151
 152     ops->calc = (step_log2 <= -10) ?
 153         mandelbrot_calc_high_prec : mandelbrot_calc_low_prec;
 154
 155     return 0;
 156 }
 157
 158 static void mandelbrot_init(void)
 159 {
 160     ops = &mandelbrot_ops;
 161
 162     x_min = MB_XOFS - MB_XFAC;
 163     x_max = MB_XOFS + MB_XFAC;
 164     y_min = -MB_YFAC;
 165     y_max = MB_YFAC;
 166
 167     recalc_parameters();
 168 }
 169
 170 static int mandelbrot_calc_low_prec(struct fractal_rect *rect,
 171         int (*button_yield_cb)(void *), void *button_yield_ctx)
 172 {
 173 #ifndef USEGSLIB
 174     long next_update = *rb->current_tick;
 175     int last_px = rect->px_min;
 176 #endif
 177     unsigned n_iter;
 178     long a32, b32;
 179     short x, x2, y, y2, a, b;
 180     int p_x, p_y;
 181     unsigned long last_yield = *rb->current_tick;
 182     unsigned long last_button_yield = *rb->current_tick;
 183
 184     a32 = x_min + x_step * rect->px_min;
 185
 186     for (p_x = rect->px_min; p_x < rect->px_max; p_x++)
 187     {
 188         a = a32 >> 16;
 189
 190         b32 = y_min + y_step * (LCD_HEIGHT - rect->py_max);
 191
 192         for (p_y = rect->py_max - 1; p_y >= rect->py_min; p_y--)
 193         {
 194             b = b32 >> 16;
 195             x = a;
 196             y = b;
 197             n_iter = 0;
 198
 199             while (++n_iter <= max_iter)
 200             {
 201                 x2 = MULS16_ASR10(x, x);
 202                 y2 = MULS16_ASR10(y, y);
 203
 204                 if (x2 + y2 > (4<<10)) break;
 205
 206                 y = 2 * MULS16_ASR10(x, y) + b;
 207                 x = x2 - y2 + a;
 208             }
 209
 210             if  (n_iter > max_iter)
 211                 imgbuffer[p_y] = CONVERGENCE_COLOR;
 212             else
 213                 imgbuffer[p_y] = COLOR(n_iter);
 214
 215             /* be nice to other threads:
 216              * if at least one tick has passed, yield */
 217             if  (TIME_AFTER(*rb->current_tick, last_yield))
 218             {
 219                 rb->yield();
 220                 last_yield = *rb->current_tick;
 221             }
 222
 223             if (TIME_AFTER(*rb->current_tick, last_button_yield))
 224             {
 225                 if (button_yield_cb(button_yield_ctx))
 226                 {
 227 #ifndef USEGSLIB
 228                     /* update screen part that was changed since last yield */
 229                     rb->lcd_update_rect(last_px, rect->py_min,
 230                             p_x - last_px + 1, rect->py_max - rect->py_min);
 231 #endif
 232                     rect->px_min = (p_x == 0) ? 0 : p_x - 1;
 233
 234                     return 1;
 235                 }
 236
 237                 last_button_yield = *rb->current_tick + BUTTON_YIELD_TIMEOUT;
 238             }
 239
 240             b32 += y_step;
 241         }
 242 #ifdef USEGSLIB
 243         grey_ub_gray_bitmap_part(imgbuffer, 0, rect->py_min, 1,
 244                 p_x, rect->py_min, 1, rect->py_max - rect->py_min);
 245 #else
 246         rb->lcd_bitmap_part(imgbuffer, 0, rect->py_min, 1,
 247                 p_x, rect->py_min, 1, rect->py_max - rect->py_min);
 248
 249         if ((p_x == rect->px_max - 1) ||
 250                 TIME_AFTER(*rb->current_tick, next_update))
 251         {
 252             next_update = *rb->current_tick + UPDATE_FREQ;
 253
 254             /* update screen part that was changed since last yield */
 255             rb->lcd_update_rect(last_px, rect->py_min,
 256                     p_x - last_px + 1, rect->py_max - rect->py_min);
 257             last_px = p_x;
 258         }
 259 #endif
 260
 261         a32 += x_step;
 262     }
 263
 264     rect->valid = 0;
 265
 266     return 0;
 267 }
 268
 269 static int mandelbrot_calc_high_prec(struct fractal_rect *rect,
 270         int (*button_yield_cb)(void *), void *button_yield_ctx)
 271 {
 272 #ifndef USEGSLIB
 273     long next_update = *rb->current_tick;
 274     int last_px = rect->px_min;
 275 #endif
 276     unsigned n_iter;
 277     long x, x2, y, y2, a, b;
 278     int p_x, p_y;
 279     unsigned long last_yield = *rb->current_tick;
 280     unsigned long last_button_yield = *rb->current_tick;
 281
 282     MULS32_INIT();
 283
 284     a = x_min + x_step * rect->px_min;
 285
 286     for (p_x = rect->px_min; p_x < rect->px_max; p_x++)
 287     {
 288         b = y_min + y_step * (LCD_HEIGHT - rect->py_max);
 289
 290         for (p_y = rect->py_max - 1; p_y >= rect->py_min; p_y--)
 291         {
 292             x = a;
 293             y = b;
 294             n_iter = 0;
 295
 296             while (++n_iter <= max_iter)
 297             {
 298                 x2 = MULS32_ASR26(x, x);
 299                 y2 = MULS32_ASR26(y, y);
 300
 301                 if (x2 + y2 > (4L<<26)) break;
 302
 303                 y = 2 * MULS32_ASR26(x, y) + b;
 304                 x = x2 - y2 + a;
 305             }
 306
 307             if  (n_iter > max_iter)
 308                 imgbuffer[p_y] = CONVERGENCE_COLOR;
 309             else
 310                 imgbuffer[p_y] = COLOR(n_iter);
 311
 312             /* be nice to other threads:
 313              * if at least one tick has passed, yield */
 314             if  (TIME_AFTER(*rb->current_tick, last_yield))
 315             {
 316                 rb->yield();
 317                 last_yield = *rb->current_tick;
 318             }
 319
 320             if (TIME_AFTER(*rb->current_tick, last_button_yield))
 321             {
 322                 if (button_yield_cb(button_yield_ctx))
 323                 {
 324 #ifndef USEGSLIB
 325                     /* update screen part that was changed since last yield */
 326                     rb->lcd_update_rect(last_px, rect->py_min,
 327                             p_x - last_px + 1, rect->py_max - rect->py_min);
 328 #endif
 329                     rect->px_min = (p_x == 0) ? 0 : p_x - 1;
 330
 331                     return 1;
 332                 }
 333
 334                 last_button_yield = *rb->current_tick + BUTTON_YIELD_TIMEOUT;
 335             }
 336
 337             b += y_step;
 338         }
 339 #ifdef USEGSLIB
 340         grey_ub_gray_bitmap_part(imgbuffer, 0, rect->py_min, 1,
 341                 p_x, rect->py_min, 1, rect->py_max - rect->py_min);
 342 #else
 343         rb->lcd_bitmap_part(imgbuffer, 0, rect->py_min, 1,
 344                 p_x, rect->py_min, 1, rect->py_max - rect->py_min);
 345
 346         if ((p_x == rect->px_max - 1) ||
 347                 TIME_AFTER(*rb->current_tick, next_update))
 348         {
 349             next_update = *rb->current_tick + UPDATE_FREQ;
 350
 351             /* update screen part that was changed since last yield */
 352             rb->lcd_update_rect(last_px, rect->py_min,
 353                     p_x - last_px + 1, rect->py_max - rect->py_min);
 354             last_px = p_x;
 355         }
 356 #endif
 357         a += x_step;
 358     }
 359
 360     rect->valid = 0;
 361
 362     return 0;
 363 }
 364
 365 static void mandelbrot_move(int x_factor, int y_factor)
 366 {
 367     long dx = (long)x_factor * x_delta;
 368     long dy = (long)y_factor * y_delta;
 369
 370     x_min += dx;
 371     x_max += dx;
 372     y_min += dy;
 373     y_max += dy;
 374 }
 375
 376 static int mandelbrot_zoom(int factor)
 377 {
 378     int res;
 379     long factor_x = (long)factor * x_delta;
 380     long factor_y = (long)factor * y_delta;
 381
 382     x_min += factor_x;
 383     x_max -= factor_x;
 384     y_min += factor_y;
 385     y_max -= factor_y;
 386
 387     res = recalc_parameters();
 388     if (res) /* zoom not possible, revert */
 389     {
 390         mandelbrot_zoom(-factor);
 391     }
 392
 393     return res;
 394 }
 395
 396 static int mandelbrot_precision(int d)
 397 {
 398     int changed = 0;
 399
 400     /* Increase precision */
 401     for (; d > 0; d--)
 402     {
 403         max_iter += max_iter / 2;
 404         changed = 1;
 405     }
 406
 407     /* Decrease precision */
 408     for (; d < 0 && max_iter >= 15; d++)
 409     {
 410         max_iter -= max_iter / 3;
 411         changed = 1;
 412     }
 413
 414     return changed;
 415 }
 416