brushlib/brush.hpp

   1 /* brushlib - The MyPaint Brush Library
   2  * Copyright (C) 2007-2011 Martin Renold <martinxyz@gmx.ch>
   3  *
   4  * Permission to use, copy, modify, and/or distribute this software for any
   5  * purpose with or without fee is hereby granted, provided that the above
   6  * copyright notice and this permission notice appear in all copies.
   7  *
   8  * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
   9  * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
  10  * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
  11  * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
  12  * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
  13  * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
  14  * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
  15  */
  16
  17 #include <stdio.h>
  18 #include <string.h>
  19 #include <glib.h>
  20 #include <math.h>
  21 #include "Python.h"
  22
  23 #include "brushsettings.hpp"
  24 #include "mapping.hpp"
  25
  26 #define ACTUAL_RADIUS_MIN 0.2
  27 #define ACTUAL_RADIUS_MAX 800 // safety guard against radius like 1e20 and against rendering overload with unexpected brush dynamics
  28
  29 /* The Brush class stores two things:
  30    b) settings: constant during a stroke (eg. size, spacing, dynamics, color selected by the user)
  31    a) states: modified during a stroke (eg. speed, smudge colors, time/distance to next dab, position filter states)
  32
  33    FIXME: Actually those are two orthogonal things. Should separate them:
  34           a) brush settings class that is saved/loaded/selected  (without states)
  35           b) brush core class to draw the dabs (using an instance of the above)
  36
  37    In python, there are two kinds of instances from this: a "global
  38    brush" which does the cursor tracking, and the "brushlist" where
  39    the states are ignored. When a brush is selected, its settings are
  40    copied into the global one, leaving the state intact.
  41  */
  42
  43
  44 class Brush {
  45 public:
  46   bool print_inputs; // debug menu
  47   // for stroke splitting (undo/redo)
  48   double stroke_total_painting_time;
  49   double stroke_current_idling_time;
  50
  51 private:
  52   // see also brushsettings.py
  53
  54   // the states (get_state, set_state, reset) that change during a stroke
  55   float states[STATE_COUNT];
  56   GRand * rng;
  57
  58   // Those mappings describe how to calculate the current value for each setting.
  59   // Most of settings will be constant (eg. only their base_value is used).
  60   Mapping * settings[BRUSH_SETTINGS_COUNT];
  61
  62   // the current value of all settings (calculated using the current state)
  63   float settings_value[BRUSH_SETTINGS_COUNT];
  64
  65   // cached calculation results
  66   float speed_mapping_gamma[2], speed_mapping_m[2], speed_mapping_q[2];
  67
  68   bool reset_requested;
  69
  70 public:
  71   Brush() {
  72     for (int i=0; i<BRUSH_SETTINGS_COUNT; i++) {
  73       settings[i] = new Mapping(INPUT_COUNT);
  74     }
  75     rng = g_rand_new();
  76     print_inputs = false;
  77
  78     for (int i=0; i<STATE_COUNT; i++) {
  79       states[i] = 0;
  80     }
  81     new_stroke();
  82
  83     settings_base_values_have_changed();
  84
  85     reset_requested = true;
  86   }
  87
  88   ~Brush() {
  89     for (int i=0; i<BRUSH_SETTINGS_COUNT; i++) {
  90       delete settings[i];
  91     }
  92     g_rand_free (rng); rng = NULL;
  93   }
  94
  95   void reset()
  96   {
  97     reset_requested = true;
  98   }
  99
 100   void new_stroke()
 101   {
 102     stroke_current_idling_time = 0;
 103     stroke_total_painting_time = 0;
 104   }
 105
 106   void set_base_value (int id, float value) {
 107     assert (id >= 0 && id < BRUSH_SETTINGS_COUNT);
 108     settings[id]->base_value = value;
 109
 110     settings_base_values_have_changed ();
 111   }
 112
 113   void set_mapping_n (int id, int input, int n) {
 114     assert (id >= 0 && id < BRUSH_SETTINGS_COUNT);
 115     settings[id]->set_n (input, n);
 116   }
 117
 118   void set_mapping_point (int id, int input, int index, float x, float y) {
 119     assert (id >= 0 && id < BRUSH_SETTINGS_COUNT);
 120     settings[id]->set_point (input, index, x, y);
 121   }
 122
 123   float get_state (int i)
 124   {
 125     assert (i >= 0 && i < STATE_COUNT);
 126     return states[i];
 127   }
 128
 129   void set_state (int i, float value)
 130   {
 131     assert (i >= 0 && i < STATE_COUNT);
 132     states[i] = value;
 133   }
 134
 135 private:
 136   // returns the fraction still left after t seconds
 137   float exp_decay (float T_const, float t)
 138   {
 139     // the argument might not make mathematical sense (whatever.)
 140     if (T_const <= 0.001) {
 141       return 0.0;
 142     } else {
 143       return exp(- t / T_const);
 144     }
 145   }
 146
 147
 148   void settings_base_values_have_changed ()
 149   {
 150     // precalculate stuff that does not change dynamically
 151
 152     // Precalculate how the physical speed will be mapped to the speed input value.
 153     // The forumla for this mapping is:
 154     //
 155     // y = log(gamma+x)*m + q;
 156     //
 157     // x: the physical speed (pixels per basic dab radius)
 158     // y: the speed input that will be reported
 159     // gamma: parameter set by ths user (small means a logarithmic mapping, big linear)
 160     // m, q: parameters to scale and translate the curve
 161     //
 162     // The code below calculates m and q given gamma and two hardcoded constraints.
 163     //
 164     for (int i=0; i<2; i++) {
 165       float gamma;
 166       gamma = settings[(i==0)?BRUSH_SPEED1_GAMMA:BRUSH_SPEED2_GAMMA]->base_value;
 167       gamma = exp(gamma);
 168
 169       float fix1_x, fix1_y, fix2_x, fix2_dy;
 170       fix1_x = 45.0;
 171       fix1_y = 0.5;
 172       fix2_x = 45.0;
 173       fix2_dy = 0.015;
 174
 175       float m, q;
 176       float c1;
 177       c1 = log(fix1_x+gamma);
 178       m = fix2_dy * (fix2_x + gamma);
 179       q = fix1_y - m*c1;
 180
 181       speed_mapping_gamma[i] = gamma;
 182       speed_mapping_m[i] = m;
 183       speed_mapping_q[i] = q;
 184     }
 185   }
 186
 187   // This function runs a brush "simulation" step. Usually it is
 188   // called once or twice per dab. In theory the precision of the
 189   // "simulation" gets better when it is called more often. In
 190   // practice this only matters if there are some highly nonlinear
 191   // mappings in critical places or extremely few events per second.
 192   //
 193   // note: parameters are is dx/ddab, ..., dtime/ddab (dab is the number, 5.0 = 5th dab)
 194   void update_states_and_setting_values (float step_dx, float step_dy, float step_dpressure, float step_declination, float step_ascension, float step_dtime)
 195   {
 196     float pressure;
 197     float inputs[INPUT_COUNT];
 198
 199     if (step_dtime < 0.0) {
 200       printf("Time is running backwards!\n");
 201       step_dtime = 0.001;
 202     } else if (step_dtime == 0.0) {
 203       // FIXME: happens about every 10th start, workaround (against division by zero)
 204       step_dtime = 0.001;
 205     }
 206
 207     states[STATE_X]        += step_dx;
 208     states[STATE_Y]        += step_dy;
 209     states[STATE_PRESSURE] += step_dpressure;
 210
 211     states[STATE_DECLINATION] += step_declination;
 212     states[STATE_ASCENSION] += step_ascension;
 213
 214     float base_radius = expf(settings[BRUSH_RADIUS_LOGARITHMIC]->base_value);
 215
 216     // FIXME: does happen (interpolation problem?)
 217     states[STATE_PRESSURE] = CLAMP(states[STATE_PRESSURE], 0.0, 1.0);
 218     pressure = states[STATE_PRESSURE];
 219
 220     { // start / end stroke (for "stroke" input only)
 221       if (!states[STATE_STROKE_STARTED]) {
 222         if (pressure > settings[BRUSH_STROKE_THRESHOLD]->base_value + 0.0001) {
 223           // start new stroke
 224           //printf("stroke start %f\n", pressure);
 225           states[STATE_STROKE_STARTED] = 1;
 226           states[STATE_STROKE] = 0.0;
 227         }
 228       } else {
 229         if (pressure <= settings[BRUSH_STROKE_THRESHOLD]->base_value * 0.9 + 0.0001) {
 230           // end stroke
 231           //printf("stroke end\n");
 232           states[STATE_STROKE_STARTED] = 0;
 233         }
 234       }
 235     }
 236
 237     // now follows input handling
 238
 239     float norm_dx, norm_dy, norm_dist, norm_speed;
 240     norm_dx = step_dx / step_dtime / base_radius;
 241     norm_dy = step_dy / step_dtime / base_radius;
 242     norm_speed = sqrt(SQR(norm_dx) + SQR(norm_dy));
 243     norm_dist = norm_speed * step_dtime;
 244
 245     inputs[INPUT_PRESSURE] = pressure;
 246     inputs[INPUT_SPEED1] = log(speed_mapping_gamma[0] + states[STATE_NORM_SPEED1_SLOW])*speed_mapping_m[0] + speed_mapping_q[0];
 247     inputs[INPUT_SPEED2] = log(speed_mapping_gamma[1] + states[STATE_NORM_SPEED2_SLOW])*speed_mapping_m[1] + speed_mapping_q[1];
 248     inputs[INPUT_RANDOM] = g_rand_double (rng);
 249     inputs[INPUT_STROKE] = MIN(states[STATE_STROKE], 1.0);
 250     inputs[INPUT_DIRECTION] = fmodf (atan2f (states[STATE_DIRECTION_DY], states[STATE_DIRECTION_DX])/(2*M_PI)*360 + 180.0, 180.0);
 251     inputs[INPUT_TILT_DECLINATION] = states[STATE_DECLINATION];
 252     inputs[INPUT_TILT_ASCENSION] = states[STATE_ASCENSION];
 253     inputs[INPUT_CUSTOM] = states[STATE_CUSTOM_INPUT];
 254     if (print_inputs) {
 255       g_print("press=% 4.3f, speed1=% 4.4f\tspeed2=% 4.4f\tstroke=% 4.3f\tcustom=% 4.3f\n", (double)inputs[INPUT_PRESSURE], (double)inputs[INPUT_SPEED1], (double)inputs[INPUT_SPEED2], (double)inputs[INPUT_STROKE], (double)inputs[INPUT_CUSTOM]);
 256     }
 257     // FIXME: this one fails!!!
 258     //assert(inputs[INPUT_SPEED1] >= 0.0 && inputs[INPUT_SPEED1] < 1e8); // checking for inf
 259
 260     for (int i=0; i<BRUSH_SETTINGS_COUNT; i++) {
 261       settings_value[i] = settings[i]->calculate (inputs);
 262     }
 263
 264     {
 265       float fac = 1.0 - exp_decay (settings_value[BRUSH_SLOW_TRACKING_PER_DAB], 1.0);
 266       states[STATE_ACTUAL_X] += (states[STATE_X] - states[STATE_ACTUAL_X]) * fac; // FIXME: should this depend on base radius?
 267       states[STATE_ACTUAL_Y] += (states[STATE_Y] - states[STATE_ACTUAL_Y]) * fac;
 268     }
 269
 270     { // slow speed
 271       float fac;
 272       fac = 1.0 - exp_decay (settings_value[BRUSH_SPEED1_SLOWNESS], step_dtime);
 273       states[STATE_NORM_SPEED1_SLOW] += (norm_speed - states[STATE_NORM_SPEED1_SLOW]) * fac;
 274       fac = 1.0 - exp_decay (settings_value[BRUSH_SPEED2_SLOWNESS], step_dtime);
 275       states[STATE_NORM_SPEED2_SLOW] += (norm_speed - states[STATE_NORM_SPEED2_SLOW]) * fac;
 276     }
 277
 278     { // slow speed, but as vector this time
 279
 280       // FIXME: offset_by_speed should be removed.
 281       //   Is it broken, non-smooth, system-dependent math?!
 282       //   A replacement could be a directed random offset.
 283
 284       float time_constant = exp(settings_value[BRUSH_OFFSET_BY_SPEED_SLOWNESS]*0.01)-1.0;
 285       // Workaround for a bug that happens mainly on Windows, causing
 286       // individual dabs to be placed far far away. Using the speed
 287       // with zero filtering is just asking for trouble anyway.
 288       if (time_constant < 0.002) time_constant = 0.002;
 289       float fac = 1.0 - exp_decay (time_constant, step_dtime);
 290       states[STATE_NORM_DX_SLOW] += (norm_dx - states[STATE_NORM_DX_SLOW]) * fac;
 291       states[STATE_NORM_DY_SLOW] += (norm_dy - states[STATE_NORM_DY_SLOW]) * fac;
 292     }
 293
 294     { // orientation (similar lowpass filter as above, but use dabtime instead of wallclock time)
 295       float dx = step_dx / base_radius;
 296       float dy = step_dy / base_radius;
 297       float step_in_dabtime = hypotf(dx, dy); // FIXME: are we recalculating something here that we already have?
 298       float fac = 1.0 - exp_decay (exp(settings_value[BRUSH_DIRECTION_FILTER]*0.5)-1.0, step_in_dabtime);
 299
 300       float dx_old = states[STATE_DIRECTION_DX];
 301       float dy_old = states[STATE_DIRECTION_DY];
 302       // use the opposite speed vector if it is closer (we don't care about 180 degree turns)
 303       if (SQR(dx_old-dx) + SQR(dy_old-dy) > SQR(dx_old-(-dx)) + SQR(dy_old-(-dy))) {
 304         dx = -dx;
 305         dy = -dy;
 306       }
 307       states[STATE_DIRECTION_DX] += (dx - states[STATE_DIRECTION_DX]) * fac;
 308       states[STATE_DIRECTION_DY] += (dy - states[STATE_DIRECTION_DY]) * fac;
 309     }
 310
 311     { // custom input
 312       float fac;
 313       fac = 1.0 - exp_decay (settings_value[BRUSH_CUSTOM_INPUT_SLOWNESS], 0.1);
 314       states[STATE_CUSTOM_INPUT] += (settings_value[BRUSH_CUSTOM_INPUT] - states[STATE_CUSTOM_INPUT]) * fac;
 315     }
 316
 317     { // stroke length
 318       float frequency;
 319       float wrap;
 320       frequency = expf(-settings_value[BRUSH_STROKE_DURATION_LOGARITHMIC]);
 321       states[STATE_STROKE] += norm_dist * frequency;
 322       // can happen, probably caused by rounding
 323       if (states[STATE_STROKE] < 0) states[STATE_STROKE] = 0;
 324       wrap = 1.0 + settings_value[BRUSH_STROKE_HOLDTIME];
 325       if (states[STATE_STROKE] > wrap) {
 326         if (wrap > 9.9 + 1.0) {
 327           // "inifinity", just hold stroke somewhere >= 1.0
 328           states[STATE_STROKE] = 1.0;
 329         } else {
 330           states[STATE_STROKE] = fmodf(states[STATE_STROKE], wrap);
 331           // just in case
 332           if (states[STATE_STROKE] < 0) states[STATE_STROKE] = 0;
 333         }
 334       }
 335     }
 336
 337     // calculate final radius
 338     float radius_log;
 339     radius_log = settings_value[BRUSH_RADIUS_LOGARITHMIC];
 340     states[STATE_ACTUAL_RADIUS] = expf(radius_log);
 341     if (states[STATE_ACTUAL_RADIUS] < ACTUAL_RADIUS_MIN) states[STATE_ACTUAL_RADIUS] = ACTUAL_RADIUS_MIN;
 342     if (states[STATE_ACTUAL_RADIUS] > ACTUAL_RADIUS_MAX) states[STATE_ACTUAL_RADIUS] = ACTUAL_RADIUS_MAX;
 343
 344     // aspect ratio (needs to be caluclated here because it can affect the dab spacing)
 345     states[STATE_ACTUAL_ELLIPTICAL_DAB_RATIO] = settings_value[BRUSH_ELLIPTICAL_DAB_RATIO];
 346     states[STATE_ACTUAL_ELLIPTICAL_DAB_ANGLE] = settings_value[BRUSH_ELLIPTICAL_DAB_ANGLE];
 347   }
 348
 349   // Called only from stroke_to(). Calculate everything needed to
 350   // draw the dab, then let the surface do the actual drawing.
 351   //
 352   // This is only gets called right after update_states_and_setting_values().
 353   // Returns true if the surface was modified.
 354   bool prepare_and_draw_dab (Surface * surface)
 355   {
 356     float x, y, opaque;
 357     float radius;
 358
 359     // ensure we don't get a positive result with two negative opaque values
 360     if (settings_value[BRUSH_OPAQUE] < 0) settings_value[BRUSH_OPAQUE] = 0;
 361     opaque = settings_value[BRUSH_OPAQUE] * settings_value[BRUSH_OPAQUE_MULTIPLY];
 362     opaque = CLAMP(opaque, 0.0, 1.0);
 363     //if (opaque == 0.0) return false; <-- cannot do that, since we need to update smudge state.
 364     if (settings_value[BRUSH_OPAQUE_LINEARIZE]) {
 365       // OPTIMIZE: no need to recalculate this for each dab
 366       float alpha, beta, alpha_dab, beta_dab;
 367       float dabs_per_pixel;
 368       // dabs_per_pixel is just estimated roughly, I didn't think hard
 369       // about the case when the radius changes during the stroke
 370       dabs_per_pixel = (
 371                         settings[BRUSH_DABS_PER_ACTUAL_RADIUS]->base_value +
 372                         settings[BRUSH_DABS_PER_BASIC_RADIUS]->base_value
 373                         ) * 2.0;
 374
 375       // the correction is probably not wanted if the dabs don't overlap
 376       if (dabs_per_pixel < 1.0) dabs_per_pixel = 1.0;
 377
 378       // interpret the user-setting smoothly
 379       dabs_per_pixel = 1.0 + settings[BRUSH_OPAQUE_LINEARIZE]->base_value*(dabs_per_pixel-1.0);
 380
 381       // see doc/brushdab_saturation.png
 382       //      beta = beta_dab^dabs_per_pixel
 383       // <==> beta_dab = beta^(1/dabs_per_pixel)
 384       alpha = opaque;
 385       beta = 1.0-alpha;
 386       beta_dab = powf(beta, 1.0/dabs_per_pixel);
 387       alpha_dab = 1.0-beta_dab;
 388       opaque = alpha_dab;
 389     }
 390
 391     x = states[STATE_ACTUAL_X];
 392     y = states[STATE_ACTUAL_Y];
 393
 394     float base_radius = expf(settings[BRUSH_RADIUS_LOGARITHMIC]->base_value);
 395
 396     if (settings_value[BRUSH_OFFSET_BY_SPEED]) {
 397       x += states[STATE_NORM_DX_SLOW] * settings_value[BRUSH_OFFSET_BY_SPEED] * 0.1 * base_radius;
 398       y += states[STATE_NORM_DY_SLOW] * settings_value[BRUSH_OFFSET_BY_SPEED] * 0.1 * base_radius;
 399     }
 400
 401     if (settings_value[BRUSH_OFFSET_BY_RANDOM]) {
 402       float amp = settings_value[BRUSH_OFFSET_BY_RANDOM];
 403       if (amp < 0.0) amp = 0.0;
 404       x += rand_gauss (rng) * amp * base_radius;
 405       y += rand_gauss (rng) * amp * base_radius;
 406     }
 407
 408
 409     radius = states[STATE_ACTUAL_RADIUS];
 410     if (settings_value[BRUSH_RADIUS_BY_RANDOM]) {
 411       float radius_log, alpha_correction;
 412       // go back to logarithmic radius to add the noise
 413       radius_log  = settings_value[BRUSH_RADIUS_LOGARITHMIC];
 414       radius_log += rand_gauss (rng) * settings_value[BRUSH_RADIUS_BY_RANDOM];
 415       radius = expf(radius_log);
 416       radius = CLAMP(radius, ACTUAL_RADIUS_MIN, ACTUAL_RADIUS_MAX);
 417       alpha_correction = states[STATE_ACTUAL_RADIUS] / radius;
 418       alpha_correction = SQR(alpha_correction);
 419       if (alpha_correction <= 1.0) {
 420         opaque *= alpha_correction;
 421       }
 422     }
 423
 424     // color part
 425
 426     float color_h = settings[BRUSH_COLOR_H]->base_value;
 427     float color_s = settings[BRUSH_COLOR_S]->base_value;
 428     float color_v = settings[BRUSH_COLOR_V]->base_value;
 429     float eraser_target_alpha = 1.0;
 430     if (settings_value[BRUSH_SMUDGE] > 0.0) {
 431       // mix (in RGB) the smudge color with the brush color
 432       hsv_to_rgb_float (&color_h, &color_s, &color_v);
 433       float fac = settings_value[BRUSH_SMUDGE];
 434       if (fac > 1.0) fac = 1.0;
 435       // If the smudge color somewhat transparent, then the resulting
 436       // dab will do erasing towards that transparency level.
 437       // see also ../doc/smudge_math.png
 438       eraser_target_alpha = (1-fac)*1.0 + fac*states[STATE_SMUDGE_A];
 439       // fix rounding errors (they really seem to happen in the previous line)
 440       eraser_target_alpha = CLAMP(eraser_target_alpha, 0.0, 1.0);
 441       if (eraser_target_alpha > 0) {
 442         color_h = (fac*states[STATE_SMUDGE_RA] + (1-fac)*color_h) / eraser_target_alpha;
 443         color_s = (fac*states[STATE_SMUDGE_GA] + (1-fac)*color_s) / eraser_target_alpha;
 444         color_v = (fac*states[STATE_SMUDGE_BA] + (1-fac)*color_v) / eraser_target_alpha;
 445       } else {
 446         // we are only erasing; the color does not matter
 447         color_h = 1.0;
 448         color_s = 0.0;
 449         color_v = 0.0;
 450       }
 451       rgb_to_hsv_float (&color_h, &color_s, &color_v);
 452     }
 453
 454     if (settings_value[BRUSH_SMUDGE_LENGTH] < 1.0 and
 455         // optimization, since normal brushes have smudge_length == 0.5 without actually smudging
 456         (settings_value[BRUSH_SMUDGE] != 0.0 or not settings[BRUSH_SMUDGE]->is_constant())) {
 457
 458       float fac = settings_value[BRUSH_SMUDGE_LENGTH];
 459       if (fac < 0.0) fac = 0;
 460       int px, py;
 461       px = ROUND(x);
 462       py = ROUND(y);
 463
 464       // Calling get_color() is almost as expensive as rendering a
 465       // dab. Because of this we use the previous value if it is not
 466       // expected to hurt quality too much.
 467       float r, g, b, a;
 468       states[STATE_LAST_GETCOLOR_RECENTNESS] *= fac;
 469       if (states[STATE_LAST_GETCOLOR_RECENTNESS] < 0.5*fac) {
 470         states[STATE_LAST_GETCOLOR_RECENTNESS] = 1.0;
 471
 472         float smudge_radius = radius * expf(settings_value[BRUSH_SMUDGE_RADIUS_LOG]);
 473         smudge_radius = CLAMP(smudge_radius, ACTUAL_RADIUS_MIN, ACTUAL_RADIUS_MAX);
 474         surface->get_color (px, py, smudge_radius, &r, &g, &b, &a);
 475
 476         states[STATE_LAST_GETCOLOR_R] = r;
 477         states[STATE_LAST_GETCOLOR_G] = g;
 478         states[STATE_LAST_GETCOLOR_B] = b;
 479         states[STATE_LAST_GETCOLOR_A] = a;
 480       } else {
 481         r = states[STATE_LAST_GETCOLOR_R];
 482         g = states[STATE_LAST_GETCOLOR_G];
 483         b = states[STATE_LAST_GETCOLOR_B];
 484         a = states[STATE_LAST_GETCOLOR_A];
 485       }
 486
 487       // updated the smudge color (stored with premultiplied alpha)
 488       states[STATE_SMUDGE_A ] = fac*states[STATE_SMUDGE_A ] + (1-fac)*a;
 489       // fix rounding errors
 490       states[STATE_SMUDGE_A ] = CLAMP(states[STATE_SMUDGE_A], 0.0, 1.0);
 491
 492       states[STATE_SMUDGE_RA] = fac*states[STATE_SMUDGE_RA] + (1-fac)*r*a;
 493       states[STATE_SMUDGE_GA] = fac*states[STATE_SMUDGE_GA] + (1-fac)*g*a;
 494       states[STATE_SMUDGE_BA] = fac*states[STATE_SMUDGE_BA] + (1-fac)*b*a;
 495     }
 496
 497     // eraser
 498     if (settings_value[BRUSH_ERASER]) {
 499       eraser_target_alpha *= (1.0-settings_value[BRUSH_ERASER]);
 500     }
 501
 502     // HSV color change
 503     color_h += settings_value[BRUSH_CHANGE_COLOR_H];
 504     color_s += settings_value[BRUSH_CHANGE_COLOR_HSV_S];
 505     color_v += settings_value[BRUSH_CHANGE_COLOR_V];
 506
 507     // HSL color change
 508     if (settings_value[BRUSH_CHANGE_COLOR_L] || settings_value[BRUSH_CHANGE_COLOR_HSL_S]) {
 509       // (calculating way too much here, can be optimized if neccessary)
 510       // this function will CLAMP the inputs
 511       hsv_to_rgb_float (&color_h, &color_s, &color_v);
 512       rgb_to_hsl_float (&color_h, &color_s, &color_v);
 513       color_v += settings_value[BRUSH_CHANGE_COLOR_L];
 514       color_s += settings_value[BRUSH_CHANGE_COLOR_HSL_S];
 515       hsl_to_rgb_float (&color_h, &color_s, &color_v);
 516       rgb_to_hsv_float (&color_h, &color_s, &color_v);
 517     }
 518
 519     float hardness = CLAMP(settings_value[BRUSH_HARDNESS], 0.0, 1.0);
 520
 521     // anti-aliasing attempt (works surprisingly well for ink brushes)
 522     float current_fadeout_in_pixels = radius * (1.0 - hardness);
 523     float min_fadeout_in_pixels = settings_value[BRUSH_ANTI_ALIASING];
 524     if (current_fadeout_in_pixels < min_fadeout_in_pixels) {
 525       // need to soften the brush (decrease hardness), but keep optical radius
 526       // so we tune both radius and hardness, to get the desired fadeout_in_pixels
 527       float current_optical_radius = radius - (1.0-hardness)*radius/2.0;
 528
 529       // Equation 1: (new fadeout must be equal to min_fadeout)
 530       //   min_fadeout_in_pixels = radius_new*(1.0 - hardness_new)
 531       // Equation 2: (optical radius must remain unchanged)
 532       //   current_optical_radius = radius_new - (1.0-hardness_new)*radius_new/2.0
 533       //
 534       // Solved Equation 1 for hardness_new, using Equation 2: (thanks to mathomatic)
 535       float hardness_new = ((current_optical_radius - (min_fadeout_in_pixels/2.0))/(current_optical_radius + (min_fadeout_in_pixels/2.0)));
 536       // Using Equation 1:
 537       float radius_new = (min_fadeout_in_pixels/(1.0 - hardness_new));
 538
 539       hardness = hardness_new;
 540       radius = radius_new;
 541     }
 542
 543     // the functions below will CLAMP most inputs
 544     hsv_to_rgb_float (&color_h, &color_s, &color_v);
 545     return surface->draw_dab (x, y, radius, color_h, color_s, color_v, opaque, hardness, eraser_target_alpha,
 546                               states[STATE_ACTUAL_ELLIPTICAL_DAB_RATIO], states[STATE_ACTUAL_ELLIPTICAL_DAB_ANGLE],
 547                               settings_value[BRUSH_LOCK_ALPHA]);
 548   }
 549
 550   // How many dabs will be drawn between the current and the next (x, y, pressure, +dt) position?
 551   float count_dabs_to (float x, float y, float pressure, float dt)
 552   {
 553     float xx, yy;
 554     float res1, res2, res3;
 555     float dist;
 556
 557     if (states[STATE_ACTUAL_RADIUS] == 0.0) states[STATE_ACTUAL_RADIUS] = expf(settings[BRUSH_RADIUS_LOGARITHMIC]->base_value);
 558     if (states[STATE_ACTUAL_RADIUS] < ACTUAL_RADIUS_MIN) states[STATE_ACTUAL_RADIUS] = ACTUAL_RADIUS_MIN;
 559     if (states[STATE_ACTUAL_RADIUS] > ACTUAL_RADIUS_MAX) states[STATE_ACTUAL_RADIUS] = ACTUAL_RADIUS_MAX;
 560
 561
 562     // OPTIMIZE: expf() called too often
 563     float base_radius = expf(settings[BRUSH_RADIUS_LOGARITHMIC]->base_value);
 564     if (base_radius < ACTUAL_RADIUS_MIN) base_radius = ACTUAL_RADIUS_MIN;
 565     if (base_radius > ACTUAL_RADIUS_MAX) base_radius = ACTUAL_RADIUS_MAX;
 566     //if (base_radius < 0.5) base_radius = 0.5;
 567     //if (base_radius > 500.0) base_radius = 500.0;
 568
 569     xx = x - states[STATE_X];
 570     yy = y - states[STATE_Y];
 571     //dp = pressure - pressure; // Not useful?
 572     // TODO: control rate with pressure (dabs per pressure) (dpressure is useless)
 573
 574     if (states[STATE_ACTUAL_ELLIPTICAL_DAB_RATIO] > 1.0) {
 575       // code duplication, see tiledsurface::draw_dab()
 576       float angle_rad=states[STATE_ACTUAL_ELLIPTICAL_DAB_ANGLE]/360*2*M_PI;
 577       float cs=cos(angle_rad);
 578       float sn=sin(angle_rad);
 579       float yyr=(yy*cs-xx*sn)*states[STATE_ACTUAL_ELLIPTICAL_DAB_RATIO];
 580       float xxr=yy*sn+xx*cs;
 581       dist = sqrt(yyr*yyr + xxr*xxr);
 582     } else {
 583       dist = hypotf(xx, yy);
 584     }
 585
 586     // FIXME: no need for base_value or for the range checks above IF always the interpolation
 587     //        function will be called before this one
 588     res1 = dist / states[STATE_ACTUAL_RADIUS] * settings[BRUSH_DABS_PER_ACTUAL_RADIUS]->base_value;
 589     res2 = dist / base_radius   * settings[BRUSH_DABS_PER_BASIC_RADIUS]->base_value;
 590     res3 = dt * settings[BRUSH_DABS_PER_SECOND]->base_value;
 591     return res1 + res2 + res3;
 592   }
 593
 594 public:
 595   // This function:
 596   // - is called once for each motion event
 597   // - does motion event interpolation
 598   // - paints zero, one or several dabs
 599   // - decides whether the stroke is finished (for undo/redo)
 600   // returns true if the stroke is finished or empty
 601   bool stroke_to (Surface * surface, float x, float y, float pressure, float xtilt, float ytilt, double dtime)
 602   {
 603     //printf("%f %f %f %f\n", (double)dtime, (double)x, (double)y, (double)pressure);
 604
 605     float tilt_ascension = 0.0;
 606     float tilt_declination = 90.0;
 607     if (xtilt != 0 || ytilt != 0) {
 608       // shield us from insane tilt input
 609       xtilt = CLAMP(xtilt, -1.0, 1.0);
 610       ytilt = CLAMP(ytilt, -1.0, 1.0);
 611       assert(isfinite(xtilt) && isfinite(ytilt));
 612
 613       tilt_ascension = 180.0*atan2(-xtilt, ytilt)/M_PI;
 614       float e;
 615       if (abs(xtilt) > abs(ytilt)) {
 616         e = sqrt(1+ytilt*ytilt);
 617       } else {
 618         e = sqrt(1+xtilt*xtilt);
 619       }
 620       float rad = hypot(xtilt, ytilt);
 621       float cos_alpha = rad/e;
 622       if (cos_alpha >= 1.0) cos_alpha = 1.0; // fix numerical inaccuracy
 623       tilt_declination = 180.0*acos(cos_alpha)/M_PI;
 624
 625       assert(isfinite(tilt_ascension));
 626       assert(isfinite(tilt_declination));
 627     }
 628
 629     // printf("xtilt %f, ytilt %f\n", (double)xtilt, (double)ytilt);
 630     // printf("ascension %f, declination %f\n", (double)tilt_ascension, (double)tilt_declination);
 631
 632     pressure = CLAMP(pressure, 0.0, 1.0);
 633     if (!isfinite(x) || !isfinite(y) ||
 634         (x > 1e10 || y > 1e10 || x < -1e10 || y < -1e10)) {
 635       // workaround attempt for https://gna.org/bugs/?14372
 636       g_print("Warning: ignoring brush::stroke_to with insane inputs (x = %f, y = %f)\n", (double)x, (double)y);
 637       x = 0.0;
 638       y = 0.0;
 639       pressure = 0.0;
 640     }
 641     // the assertion below is better than out-of-memory later at save time
 642     assert(x < 1e8 && y < 1e8 && x > -1e8 && y > -1e8);
 643
 644     if (dtime < 0) g_print("Time jumped backwards by dtime=%f seconds!\n", dtime);
 645     if (dtime <= 0) dtime = 0.0001; // protect against possible division by zero bugs
 646
 647     if (dtime > 0.100 && pressure && states[STATE_PRESSURE] == 0) {
 648       // Workaround for tablets that don't report motion events without pressure.
 649       // This is to avoid linear interpolation of the pressure between two events.
 650       stroke_to (surface, x, y, 0.0, 90.0, 0.0, dtime-0.0001);
 651       dtime = 0.0001;
 652     }
 653
 654     g_rand_set_seed (rng, states[STATE_RNG_SEED]);
 655
 656     { // calculate the actual "virtual" cursor position
 657
 658       // noise first
 659       if (settings[BRUSH_TRACKING_NOISE]->base_value) {
 660         // OPTIMIZE: expf() called too often
 661         float base_radius = expf(settings[BRUSH_RADIUS_LOGARITHMIC]->base_value);
 662
 663         x += rand_gauss (rng) * settings[BRUSH_TRACKING_NOISE]->base_value * base_radius;
 664         y += rand_gauss (rng) * settings[BRUSH_TRACKING_NOISE]->base_value * base_radius;
 665       }
 666
 667       float fac = 1.0 - exp_decay (settings[BRUSH_SLOW_TRACKING]->base_value, 100.0*dtime);
 668       x = states[STATE_X] + (x - states[STATE_X]) * fac;
 669       y = states[STATE_Y] + (y - states[STATE_Y]) * fac;
 670     }
 671
 672     // draw many (or zero) dabs to the next position
 673
 674     // see doc/stroke2dabs.png
 675     float dist_moved = states[STATE_DIST];
 676     float dist_todo = count_dabs_to (x, y, pressure, dtime);
 677
 678     //if (dtime > 5 || dist_todo > 300) {
 679     if (dtime > 5 || reset_requested) {
 680       reset_requested = false;
 681
 682       /*
 683         TODO:
 684         if (dist_todo > 300) {
 685         // this happens quite often, eg when moving the cursor back into the window
 686         // FIXME: bad to hardcode a distance treshold here - might look at zoomed image
 687         //        better detect leaving/entering the window and reset then.
 688         g_print ("Warning: NOT drawing %f dabs.\n", dist_todo);
 689         g_print ("dtime=%f, dx=%f\n", dtime, x-states[STATE_X]);
 690         //must_reset = 1;
 691         }
 692       */
 693
 694       //printf("Brush reset.\n");
 695       for (int i=0; i<STATE_COUNT; i++) {
 696         states[i] = 0;
 697       }
 698
 699       states[STATE_X] = x;
 700       states[STATE_Y] = y;
 701       states[STATE_PRESSURE] = pressure;
 702
 703       // not resetting, because they will get overwritten below:
 704       //dx, dy, dpress, dtime
 705
 706       states[STATE_ACTUAL_X] = states[STATE_X];
 707       states[STATE_ACTUAL_Y] = states[STATE_Y];
 708       states[STATE_STROKE] = 1.0; // start in a state as if the stroke was long finished
 709
 710       return true;
 711     }
 712
 713     //g_print("dist = %f\n", states[STATE_DIST]);
 714     enum { UNKNOWN, YES, NO } painted = UNKNOWN;
 715     double dtime_left = dtime;
 716
 717     float step_dx, step_dy, step_dpressure, step_dtime;
 718     float step_declination, step_ascension;
 719     while (dist_moved + dist_todo >= 1.0) { // there are dabs pending
 720       { // linear interpolation (nonlinear variant was too slow, see SVN log)
 721         float frac; // fraction of the remaining distance to move
 722         if (dist_moved > 0) {
 723           // "move" the brush exactly to the first dab (moving less than one dab)
 724           frac = (1.0 - dist_moved) / dist_todo;
 725           dist_moved = 0;
 726         } else {
 727           // "move" the brush from one dab to the next
 728           frac = 1.0 / dist_todo;
 729         }
 730         step_dx        = frac * (x - states[STATE_X]);
 731         step_dy        = frac * (y - states[STATE_Y]);
 732         step_dpressure = frac * (pressure - states[STATE_PRESSURE]);
 733         step_dtime     = frac * (dtime_left - 0.0);
 734         step_declination = frac * (tilt_declination - states[STATE_DECLINATION]);
 735         step_ascension   = frac * (tilt_ascension - states[STATE_ASCENSION]);
 736         // Though it looks different, time is interpolated exactly like x/y/pressure.
 737       }
 738
 739       update_states_and_setting_values (step_dx, step_dy, step_dpressure, step_declination, step_ascension, step_dtime);
 740       bool painted_now = prepare_and_draw_dab (surface);
 741       if (painted_now) {
 742         painted = YES;
 743       } else if (painted == UNKNOWN) {
 744         painted = NO;
 745       }
 746
 747       dtime_left   -= step_dtime;
 748       dist_todo  = count_dabs_to (x, y, pressure, dtime_left);
 749     }
 750
 751     {
 752       // "move" the brush to the current time (no more dab will happen)
 753       // Important to do this at least once every event, because
 754       // brush_count_dabs_to depends on the radius and the radius can
 755       // depend on something that changes much faster than only every
 756       // dab (eg speed).
 757
 758       step_dx        = x - states[STATE_X];
 759       step_dy        = y - states[STATE_Y];
 760       step_dpressure = pressure - states[STATE_PRESSURE];
 761       step_declination = tilt_declination - states[STATE_DECLINATION];
 762       step_ascension = tilt_ascension - states[STATE_ASCENSION];
 763       step_dtime     = dtime_left;
 764
 765       //dtime_left = 0; but that value is not used any more
 766
 767       update_states_and_setting_values (step_dx, step_dy, step_dpressure, step_declination, step_ascension, step_dtime);
 768     }
 769
 770     // save the fraction of a dab that is already done now
 771     states[STATE_DIST] = dist_moved + dist_todo;
 772     //g_print("dist_final = %f\n", states[STATE_DIST]);
 773
 774     // next seed for the RNG (GRand has no get_state() and states[] must always contain our full state)
 775     states[STATE_RNG_SEED] = g_rand_int(rng);
 776
 777     // stroke separation logic (for undo/redo)
 778
 779     if (painted == UNKNOWN) {
 780       if (stroke_current_idling_time > 0 || stroke_total_painting_time == 0) {
 781         // still idling
 782         painted = NO;
 783       } else {
 784         // probably still painting (we get more events than brushdabs)
 785         painted = YES;
 786         //if (pressure == 0) g_print ("info: assuming 'still painting' while there is no pressure\n");
 787       }
 788     }
 789     if (painted == YES) {
 790       //if (stroke_current_idling_time > 0) g_print ("idling ==> painting\n");
 791       stroke_total_painting_time += dtime;
 792       stroke_current_idling_time = 0;
 793       // force a stroke split after some time
 794       if (stroke_total_painting_time > 4 + 3*pressure) {
 795         // but only if pressure is not being released
 796         // FIXME: use some smoothed state for dpressure, not the output of the interpolation code
 797         //        (which might easily wrongly give dpressure == 0)
 798         if (step_dpressure >= 0) {
 799           return true;
 800         }
 801       }
 802     } else if (painted == NO) {
 803       //if (stroke_current_idling_time == 0) g_print ("painting ==> idling\n");
 804       stroke_current_idling_time += dtime;
 805       if (stroke_total_painting_time == 0) {
 806         // not yet painted, start a new stroke if we have accumulated a lot of irrelevant motion events
 807         if (stroke_current_idling_time > 1.0) {
 808           return true;
 809         }
 810       } else {
 811         // Usually we have pressure==0 here. But some brushes can paint
 812         // nothing at full pressure (eg gappy lines, or a stroke that
 813         // fades out). In either case this is the prefered moment to split.
 814 //        if (stroke_total_painting_time+stroke_current_idling_time > 0.9 + 5*pressure) {
 815           return true;
 816 //        }
 817       }
 818     }
 819     return false;
 820   }
 821
 822 };