text: /* libmypaint - The MyPaint Brush Library
   * Copyright (C) 2007-2011 Martin Renold <martinxyz@gmx.ch>
   *
   * Permission to use, copy, modify, and/or distribute this software for any
   * purpose with or without fee is hereby granted, provided that the above
   * copyright notice and this permission notice appear in all copies.
   *
   * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
   * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
   * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
   * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
   * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
   * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
   * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
   */
  
  #include "config.h"
  
  #include <stdlib.h>
  #include <stdio.h>
  #include <string.h>
  #include <math.h>
  #include <assert.h>
  
  #if MYPAINT_CONFIG_USE_GLIB
  #include <glib.h>
  #include "glib/mypaint-brush.h"
  #endif
  
  #include "mypaint-brush.h"
  
  #include "mypaint-brush-settings.h"
  #include "mypaint-mapping.h"
  #include "helpers.h"
  #include "rng-double.h"
  
  #include <json.h>
  
  #ifdef _MSC_VER
  #if _MSC_VER < 1700     // Visual Studio 2012 and later has isfinite and roundf
    #include <float.h>
    static inline int    isfinite(double x) { return _finite(x); }
    static inline float  roundf  (float  x) { return x >= 0.0f ? floorf(x + 0.5f) : ceilf(x - 0.5f); }
  #endif
  #endif
  
  // Conversion from degrees to radians
  #define RADIANS(x) ((x) * M_PI / 180.0)
  // Conversion from radians to degrees
  #define DEGREES(x) (((x) / (2 * M_PI)) * 360.0)
  
  #define ACTUAL_RADIUS_MIN 0.2
  #define ACTUAL_RADIUS_MAX 1000 // safety guard against radius like 1e20 and against rendering overload with unexpected brush dynamics
  
  #define GRID_SIZE 256.0
  
  /* Named indices for the smudge bucket arrays */
  enum {
    SMUDGE_R, SMUDGE_G, SMUDGE_B, SMUDGE_A,
    PREV_COL_R, PREV_COL_G, PREV_COL_B, PREV_COL_A,
    PREV_COL_RECENTNESS,
    SMUDGE_BUCKET_SIZE
  };
  
  /* The Brush class stores two things:
     b) settings: constant during a stroke (eg. size, spacing, dynamics, color selected by the user)
     a) states: modified during a stroke (eg. speed, smudge colors, time/distance to next dab, position filter states)
  
     FIXME: Actually those are two orthogonal things. Should separate them:
            a) brush settings class that is saved/loaded/selected  (without states)
            b) brush core class to draw the dabs (using an instance of the above)
  
     In python, there are two kinds of instances from this: a "global
     brush" which does the cursor tracking, and the "brushlist" where
     the states are ignored. When a brush is selected, its settings are
     copied into the global one, leaving the state intact.
   */
  
  
  /**
    * MyPaintBrush:
    *
    * The MyPaint brush engine class.
    */
  struct MyPaintBrush {
  
      gboolean print_inputs; // debug menu
      // for stroke splitting (undo/redo)
      double stroke_total_painting_time;
      double stroke_current_idling_time;
  
      // the states (get_state, set_state, reset) that change during a stroke
      float states[MYPAINT_BRUSH_STATES_COUNT];
      // smudge bucket array: part of the state, but stored separately.
      // Usually used for brushes with multiple offset dabs, where each
      // dab is assigned its own bucket containing a smudge state.
      float *smudge_buckets;
      int num_buckets;
      int min_bucket_used;
      int max_bucket_used;
  
      double random_input;
      float skip;
      float skip_last_x;
      float skip_last_y;
      float skipped_dtime;
      RngDouble * rng;
  
      // Those mappings describe how to calculate the current value for each setting.
      // Most of settings will be constant (eg. only their base_value is used).
      MyPaintMapping * settings[MYPAINT_BRUSH_SETTINGS_COUNT];
  
      // the current value of all settings (calculated using the current state)
      float settings_value[MYPAINT_BRUSH_SETTINGS_COUNT];
  
      // see also brushsettings.py
  
      // cached calculation results
      float speed_mapping_gamma[2];
      float speed_mapping_m[2];
      float speed_mapping_q[2];
  
      gboolean reset_requested;
      json_object *brush_json;
      int refcount;
  };
  
  /* Macros for accessing states and setting values
     Although macros are never nice, simple file-local macros are warranted
     here since it massively improves code readability.
     The INPUT macro relies on 'inputs' being a float array of size MYPAINT_BRUSH_INPUTS_COUNT,
     accessible in the scope where the macro is used.
  */
  #define STATE(self, state_name) ((self)->states[MYPAINT_BRUSH_STATE_##state_name])
  #define SETTING(self, setting_name) ((self)->settings_value[MYPAINT_BRUSH_SETTING_##setting_name])
  #define BASEVAL(self, setting_name) (mypaint_mapping_get_base_value((self)->settings[MYPAINT_BRUSH_SETTING_##setting_name]))
  #define INPUT(input_name) (inputs[MYPAINT_BRUSH_INPUT_##input_name])
  
  void settings_base_values_have_changed (MyPaintBrush *self);
  
  
  #include "glib/mypaint-brush.c"
  
  void
  brush_reset(MyPaintBrush *self)
  {
      self->skip = 0;
      self->skip_last_x = 0;
      self->skip_last_y = 0;
      self->skipped_dtime = 0;
      // Clear states
      memset(self->states, 0, sizeof(self->states));
      // Set the flip state such that it will be at "1" for the first
      // dab, and then switch between -1 and 1 for the subsequent dabs.
      STATE(self, FLIP) = -1;
      // Clear smudge buckets
      if (self->smudge_buckets) {
        int min_index = self->min_bucket_used;
        if (min_index != -1) {
          int max_index = self->max_bucket_used;
          size_t num_bytes = (max_index - min_index) * sizeof(self->smudge_buckets[0]) * SMUDGE_BUCKET_SIZE;
          memset(self->smudge_buckets + min_index, 0, num_bytes);
          self->min_bucket_used = -1;
          self->max_bucket_used = -1;
        }
      }
  }
  
  /**
    * mypaint_brush_new:
    *
    * Create a new MyPaint brush engine instance.
    * Initial reference count is 1. Release references using mypaint_brush_unref()
    */
  MyPaintBrush *
  mypaint_brush_new(void)
  {
    return mypaint_brush_new_with_buckets(0);
  }
  
  /**
    * mypaint_brush_new_with_buckets:
    *
    * Create a new MyPaint brush engine instance with smudge bucket support.
    * MyPaint 2.0 supports up to 255 buckets per brush.
    * Initial reference count is 1. Release references using mypaint_brush_unref()
    */
  MyPaintBrush *
  mypaint_brush_new_with_buckets(int num_smudge_buckets)
  {
      MyPaintBrush *self = (MyPaintBrush *)malloc(sizeof(MyPaintBrush));
  
      if (!self) {
        return NULL;
      }
  
      if (num_smudge_buckets > 0) {
        float *bucket_array = malloc(num_smudge_buckets * SMUDGE_BUCKET_SIZE * sizeof(float));
        if (!bucket_array) {
          free(self);
          return NULL;
        }
        self->smudge_buckets = bucket_array;
        self->num_buckets = num_smudge_buckets;
        // Set up min/max to initialize (clear) the array in the call to brush_reset.
        self->min_bucket_used = 0;
        self->max_bucket_used = self->num_buckets - 1;
      } else {
        self->smudge_buckets = NULL;
        self->num_buckets = 0;
      }
  
      self->refcount = 1;
      for (int i = 0; i < MYPAINT_BRUSH_SETTINGS_COUNT; i++) {
        self->settings[i] = mypaint_mapping_new(MYPAINT_BRUSH_INPUTS_COUNT);
      }
      self->rng = rng_double_new(1000);
      self->random_input = 0;
      self->print_inputs = FALSE;
  
      brush_reset(self);
  
      mypaint_brush_new_stroke(self);
      settings_base_values_have_changed(self);
  
      self->reset_requested = TRUE;
  
      self->brush_json = json_object_new_object();
  
      return self;
  }
  
  void
  brush_free(MyPaintBrush *self)
  {
      for (int i = 0; i < MYPAINT_BRUSH_SETTINGS_COUNT; i++) {
          mypaint_mapping_free(self->settings[i]);
      }
      rng_double_free (self->rng);
      self->rng = NULL;
  
      if (self->brush_json) {
          json_object_put(self->brush_json);
      }
  
      free(self->smudge_buckets);
      free(self);
  }
  
  /**
    * mypaint_brush_unref: (skip)
    *
    * Decrease the reference count. Will be freed when it hits 0.
    */
  void
  mypaint_brush_unref(MyPaintBrush *self)
  {
      self->refcount--;
      if (self->refcount == 0) {
          brush_free(self);
      }
  }
  /**
    * mypaint_brush_ref: (skip)
    *
    * Increase the reference count.
    */
  void
  mypaint_brush_ref(MyPaintBrush *self)
  {
      self->refcount++;
  }
  
  /**
    * mypaint_brush_get_total_stroke_painting_time:
    *
    * Return the total amount of painting time for the current stroke.
    */
  double
  mypaint_brush_get_total_stroke_painting_time(MyPaintBrush *self)
  {
      return self->stroke_total_painting_time;
  }
  
  /**
    * mypaint_brush_set_print_inputs:
    *
    * Enable/Disable printing of brush engine inputs on stderr. Intended for debugging only.
    */
  void
  mypaint_brush_set_print_inputs(MyPaintBrush *self, gboolean enabled)
  {
      self->print_inputs = enabled;
  }
  
  /**
    * mypaint_brush_reset:
    *
    * Reset the current brush engine state.
    * Used when the next mypaint_brush_stroke_to() call is not related to the current state.
    * Note that the reset request is queued and changes in state will only happen on next stroke_to()
    */
  void
  mypaint_brush_reset(MyPaintBrush *self)
  {
      self->reset_requested = TRUE;
  }
  
  /**
    * mypaint_brush_new_stroke:
    *
    * Start a new stroke.
    */
  void
  mypaint_brush_new_stroke(MyPaintBrush *self)
  {
      self->stroke_current_idling_time = 0;
      self->stroke_total_painting_time = 0;
  }
  
  /**
    * mypaint_brush_set_base_value:
    *
    * Set the base value of a brush setting.
    */
  void
  mypaint_brush_set_base_value(MyPaintBrush *self, MyPaintBrushSetting id, float value)
  {
      assert (id < MYPAINT_BRUSH_SETTINGS_COUNT);
      mypaint_mapping_set_base_value(self->settings[id], value);
  
      settings_base_values_have_changed (self);
  }
  
  /**
    * mypaint_brush_get_base_value:
    *
    * Get the base value of a brush setting.
    */
  float
  mypaint_brush_get_base_value(MyPaintBrush *self, MyPaintBrushSetting id)
  {
      assert (id < MYPAINT_BRUSH_SETTINGS_COUNT);
      return mypaint_mapping_get_base_value(self->settings[id]);
  }
  
  /**
    * mypaint_brush_set_mapping_n:
    *
    * Set the number of points used for the dynamics mapping between a #MyPaintBrushInput and #MyPaintBrushSetting.
    */
  void
  mypaint_brush_set_mapping_n(MyPaintBrush *self, MyPaintBrushSetting id, MyPaintBrushInput input, int n)
  {
      assert (id < MYPAINT_BRUSH_SETTINGS_COUNT);
      mypaint_mapping_set_n(self->settings[id], input, n);
  }
  
  /**
    * mypaint_brush_get_mapping_n:
    *
    * Get the number of points used for the dynamics mapping between a #MyPaintBrushInput and #MyPaintBrushSetting.
    */
  int
  mypaint_brush_get_mapping_n(MyPaintBrush *self, MyPaintBrushSetting id, MyPaintBrushInput input)
  {
      return mypaint_mapping_get_n(self->settings[id], input);
  }
  
  /**
    * mypaint_brush_is_constant:
    *
    * Returns TRUE if the brush has no dynamics for the given #MyPaintBrushSetting
    */
  gboolean
  mypaint_brush_is_constant(MyPaintBrush *self, MyPaintBrushSetting id)
  {
      assert (id < MYPAINT_BRUSH_SETTINGS_COUNT);
      return mypaint_mapping_is_constant(self->settings[id]);
  }
  
  /**
    * mypaint_brush_get_inputs_used_n:
    *
    * Returns how many inputs are used for the dynamics of a #MyPaintBrushSetting
    */
  int
  mypaint_brush_get_inputs_used_n(MyPaintBrush *self, MyPaintBrushSetting id)
  {
      assert (id < MYPAINT_BRUSH_SETTINGS_COUNT);
      return mypaint_mapping_get_inputs_used_n(self->settings[id]);
  }
  
  /**
    * mypaint_brush_set_mapping_point:
    *
    * Set a X,Y point of a dynamics mapping.
    * The index must be within the number of points set using mypaint_brush_set_mapping_n()
    */
  void
  mypaint_brush_set_mapping_point(MyPaintBrush *self, MyPaintBrushSetting id, MyPaintBrushInput input, int index, float x, float y)
  {
      assert (id < MYPAINT_BRUSH_SETTINGS_COUNT);
      mypaint_mapping_set_point(self->settings[id], input, index, x, y);
  }
  
  /**
   * mypaint_brush_get_mapping_point:
   * @x: (out): Location to return the X value
   * @y: (out): Location to return the Y value
   *
   * Get a X,Y point of a dynamics mapping.
   **/
  void
  mypaint_brush_get_mapping_point(MyPaintBrush *self, MyPaintBrushSetting id, MyPaintBrushInput input, int index, float *x, float *y)
  {
      assert (id < MYPAINT_BRUSH_SETTINGS_COUNT);
      mypaint_mapping_get_point(self->settings[id], input, index, x, y);
  }
  
  /**
   * mypaint_brush_get_state:
   *
   * Get an internal brush engine state.
   * Normally used for debugging, but can be used to implement record & replay functionality.
   **/
  float
  mypaint_brush_get_state(MyPaintBrush *self, MyPaintBrushState i)
  {
      assert (i < MYPAINT_BRUSH_STATES_COUNT);
      return self->states[i];
  }
  
  /**
   * mypaint_brush_set_state:
   *
   * Set an internal brush engine state.
   * Normally used for debugging, but can be used to implement record & replay functionality.
   **/
  void
  mypaint_brush_set_state(MyPaintBrush *self, MyPaintBrushState i, float value)
  {
      assert (i < MYPAINT_BRUSH_STATES_COUNT);
      self->states[i] = value;
  }
  /**
   * mypaint_brush_set_smudge_bucket_state: (skip)
   *
   * Set an internal smudge bucket state for the bucket of the given index.
   * Returns true if the state could be updated, false otherwise.
   * Setting the state can only fail if the bucket index lies outside the
   * capacity of the brush.
   *
   **/
  gboolean
  mypaint_brush_set_smudge_bucket_state(
      MyPaintBrush* self, int bucket_index,
      float r, float g, float b, float a,
      float prev_r, float prev_g, float prev_b, float prev_a,
      float prev_color_recentness)
  {
    if (self->num_buckets > bucket_index) {
      float* bucket = &self->smudge_buckets[bucket_index];
      bucket[SMUDGE_R] = r;
      bucket[SMUDGE_G] = g;
      bucket[SMUDGE_B] = b;
      bucket[SMUDGE_A] = a;
      bucket[PREV_COL_R] = prev_r;
      bucket[PREV_COL_G] = prev_g;
      bucket[PREV_COL_B] = prev_b;
      bucket[PREV_COL_A] = prev_a;
      bucket[PREV_COL_RECENTNESS] = prev_color_recentness;
      return TRUE;
    }
    return FALSE;
  }
  
  /**
   * mypaint_brush_get_smudge_bucket_state: (skip)
   *
   * Retrieve the internal state for the smudge bucket of the given index.
   * If the state could be retrieved, populating the out variables, true is returned.
   * Getting the state can only fail if the bucket index lies outside the
   * capacity of the brush.
   *
   **/
  gboolean
  mypaint_brush_get_smudge_bucket_state(
      const MyPaintBrush* self, int bucket_index,
      float* r, float* g, float* b, float* a,
      float* prev_r, float* prev_g, float* prev_b, float* prev_a,
      float* prev_color_recentness)
  {
    if (self->num_buckets > bucket_index) {
      float* bucket = &self->smudge_buckets[bucket_index];
      *r = bucket[SMUDGE_R];
      *g = bucket[SMUDGE_G];
      *b = bucket[SMUDGE_B];
      *a = bucket[SMUDGE_A];
      *prev_r = bucket[PREV_COL_R];
      *prev_g = bucket[PREV_COL_G];
      *prev_b = bucket[PREV_COL_B];
      *prev_a = bucket[PREV_COL_A];
      *prev_color_recentness = bucket[PREV_COL_RECENTNESS];
      return TRUE;
    }
    return FALSE;
  }
  
  /**
   * mypaint_brush_get_min_smudge_bucket_used: (skip)
   *
   * Get the lowest index of the indices of buckets that have been written to.
   * If no smudge bucket has been written to, -1 is returned.
   *
   */
  int mypaint_brush_get_min_smudge_bucket_used(const MyPaintBrush* self)
  {
    return self->min_bucket_used;
  }
  
  /**
   * mypaint_brush_get_max_smudge_bucket_used: (skip)
   *
   * Get the largest index of the indices of buckets that have been written to.
   * If no smudge bucket has been written to, -1 is returned.
   *
   */
  int mypaint_brush_get_max_smudge_bucket_used(const MyPaintBrush* self)
  {
    return self->max_bucket_used;
  }
  
    // returns the fraction still left after t seconds
    float exp_decay (float T_const, float t)
    {
      // the argument might not make mathematical sense (whatever.)
      if (T_const <= 0.001) {
        return 0.0;
      }
  
      const float arg = -t / T_const;
      return expf(arg);
    }
  
  
    void settings_base_values_have_changed (MyPaintBrush *self)
    {
      // precalculate stuff that does not change dynamically
  
      // Precalculate how the physical speed will be mapped to the speed input value.
      // The formula for this mapping is:
      //
      // y = log(gamma+x)*m + q;
      //
      // x: the physical speed (pixels per basic dab radius)
      // y: the speed input that will be reported
      // gamma: parameter set by the user (small means a logarithmic mapping, big linear)
      // m, q: parameters to scale and translate the curve
      //
      // The code below calculates m and q given gamma and two hardcoded constraints.
      //
      for (int i = 0; i < 2; i++) {
        const float gamma = expf(i == 0 ? BASEVAL(self, SPEED1_GAMMA) : BASEVAL(self, SPEED2_GAMMA));
  
        const float fix1_x = 45.0;
        const float fix1_y = 0.5;
        const float fix2_x = 45.0;
        const float fix2_dy = 0.015;
  
        const float c1 = log(fix1_x + gamma);
        const float m = fix2_dy * (fix2_x + gamma);
        const float q = fix1_y - m * c1;
  
        self->speed_mapping_gamma[i] = gamma;
        self->speed_mapping_m[i] = m;
        self->speed_mapping_q[i] = q;
      }
    }
  
  typedef struct {
    float x;
    float y;
  } Offsets;
  
  Offsets
  directional_offsets(const MyPaintBrush* const self, const float base_radius, const int brush_flip)
  {
      const float offset_mult = expf(SETTING(self, OFFSET_MULTIPLIER));
      // Sanity check - it is easy to reach infinite multipliers w. logarithmic parameters
      if (!isfinite(offset_mult)) {
          Offsets offs = {0.0f, 0.0f};
          return offs;
      }
  
      float dx = SETTING(self, OFFSET_X);
      float dy = SETTING(self, OFFSET_Y);
  
      //Anti_Art offsets tweaked by BrienD.  Adjusted with ANGLE_ADJ and OFFSET_MULTIPLIER
      const float offset_angle_adj = SETTING(self, OFFSET_ANGLE_ADJ);
      const float dir_angle_dy = STATE(self, DIRECTION_ANGLE_DY);
      const float dir_angle_dx = STATE(self, DIRECTION_ANGLE_DX);
      const float angle_deg = fmodf(DEGREES(atan2f(dir_angle_dy, dir_angle_dx)) - 90, 360);
  
      //offset to one side of direction
      const float offset_angle = SETTING(self, OFFSET_ANGLE);
      if (offset_angle) {
          const float dir_angle = RADIANS(angle_deg + offset_angle_adj);
          dx += cos(dir_angle) * offset_angle;
          dy += sin(dir_angle) * offset_angle;
      }
  
      //offset to one side of ascension angle
      const float view_rotation = STATE(self, VIEWROTATION);
      const float offset_angle_asc = SETTING(self, OFFSET_ANGLE_ASC);
      if (offset_angle_asc) {
          const float ascension = STATE(self, ASCENSION);
          const float asc_angle = RADIANS(ascension - view_rotation + offset_angle_adj);
          dx += cos(asc_angle) * offset_angle_asc;
          dy += sin(asc_angle) * offset_angle_asc;
        }
  
      //offset to one side of view orientation
      const float view_offset = SETTING(self, OFFSET_ANGLE_VIEW);
      if (view_offset) {
          const float view_angle = RADIANS(view_rotation + offset_angle_adj);
          dx += cos(-view_angle) * view_offset;
          dy += sin(-view_angle) * view_offset;
      }
  
      //offset mirrored to sides of direction
      const float offset_dir_mirror = MAX(0.0, SETTING(self, OFFSET_ANGLE_2));
      if (offset_dir_mirror) {
          const float dir_mirror_angle = RADIANS(angle_deg + offset_angle_adj * brush_flip);
          const float offset_factor = offset_dir_mirror * brush_flip;
          dx += cos(dir_mirror_angle) * offset_factor;
          dy += sin(dir_mirror_angle) * offset_factor;
      }
  
      //offset mirrored to sides of ascension angle
      const float offset_asc_mirror = MAX(0.0, SETTING(self, OFFSET_ANGLE_2_ASC));
      if (offset_asc_mirror) {
          const float ascension = STATE(self, ASCENSION);
          const float asc_angle = RADIANS(ascension - view_rotation + offset_angle_adj * brush_flip);
          const float offset_factor = brush_flip * offset_asc_mirror;
          dx += cos(asc_angle) * offset_factor;
          dy += sin(asc_angle) * offset_factor;
      }
  
      //offset mirrored to sides of view orientation
      const float offset_view_mirror = MAX(0.0, SETTING(self, OFFSET_ANGLE_2_VIEW));
      if (offset_view_mirror) {
          const float offset_factor = brush_flip * offset_view_mirror;
          const float offset_angle_rad = RADIANS(view_rotation + offset_angle_adj);
          dx += cos(-offset_angle_rad) * offset_factor;
          dy += sin(-offset_angle_rad) * offset_factor;
      }
      // Clamp the final offsets to avoid potential memory issues (extreme memory use from redraws)
      // Allow offsets up to the 1080 * 3 pixels. Unlikely to hamper anyone artistically.
      const float lim = 3240;
      const float base_mul = base_radius * offset_mult;
      Offsets offs = {CLAMP(dx * base_mul, -lim, lim),  CLAMP(dy * base_mul, -lim, lim)};
      return offs;
  }
  
  // Debugging: print brush inputs/states (not all of them)
  void print_inputs(MyPaintBrush *self, float* inputs)
  {
      printf(
          "press=% 4.3f, speed1=% 4.4f\tspeed2=% 4.4f",
          INPUT(PRESSURE),
          INPUT(SPEED1),
          INPUT(SPEED2)
          );
      printf(
          "\tstroke=% 4.3f\tcustom=% 4.3f",
          INPUT(STROKE),
          INPUT(CUSTOM)
          );
      printf(
          "\tviewzoom=% 4.3f\tviewrotation=% 4.3f",
          INPUT(VIEWZOOM),
          STATE(self, VIEWROTATION)
          );
      printf(
          "\tasc=% 4.3f\tdir=% 4.3f\tdec=% 4.3f\tdabang=% 4.3f",
          INPUT(TILT_ASCENSION),
          INPUT(DIRECTION),
          INPUT(TILT_DECLINATION),
          STATE(self, ACTUAL_ELLIPTICAL_DAB_ANGLE)
          );
      printf(
          "\txtilt=% 4.3f\tytilt=% 4.3fattack=% 4.3f",
          INPUT(TILT_DECLINATIONX),
          INPUT(TILT_DECLINATIONY),
          INPUT(ATTACK_ANGLE)
          );
      printf("\n");
  }
  
    // This function runs a brush "simulation" step. Usually it is
    // called once or twice per dab. In theory the precision of the
    // "simulation" gets better when it is called more often. In
    // practice this only matters if there are some highly nonlinear
    // mappings in critical places or extremely few events per second.
    //
    // note: parameters are is dx/ddab, ..., dtime/ddab (dab is the number, 5.0 = 5th dab)
    void update_states_and_setting_values (MyPaintBrush *self, float step_ddab, float step_dx, float step_dy, float step_dpressure, float step_declination, float step_ascension, float step_dtime, float step_viewzoom, float step_viewrotation, float step_declinationx, float step_declinationy, float step_barrel_rotation)
    {
      if (step_dtime < 0.0) {
        printf("Time is running backwards!\n");
        step_dtime = 0.001;
      } else if (step_dtime == 0.0) {
        // FIXME: happens about every 10th start, workaround (against division by zero)
        step_dtime = 0.001;
      }
  
      STATE(self, X) += step_dx;
      STATE(self, Y) += step_dy;
      STATE(self, PRESSURE) += step_dpressure;
  
      STATE(self, DECLINATION) += step_declination;
      STATE(self, ASCENSION) += step_ascension;
      STATE(self, DECLINATIONX) += step_declinationx;
      STATE(self, DECLINATIONY) += step_declinationy;
  
      STATE(self, VIEWZOOM) = step_viewzoom;
      const float viewrotation = mod_arith(DEGREES(step_viewrotation) + 180.0, 360.0) - 180.0;
      STATE(self, VIEWROTATION) = viewrotation;
  
      { // Gridmap state update
          const float x = STATE(self, ACTUAL_X);
          const float y = STATE(self, ACTUAL_Y);
          const float scale = expf(SETTING(self, GRIDMAP_SCALE));
          const float scale_x = SETTING(self, GRIDMAP_SCALE_X);
          const float scale_y = SETTING(self, GRIDMAP_SCALE_Y);
          const float scaled_size = scale * GRID_SIZE;
          STATE(self, GRIDMAP_X) = mod_arith(fabsf(x * scale_x), scaled_size) / scaled_size * GRID_SIZE;
          STATE(self, GRIDMAP_Y) = mod_arith(fabsf(y * scale_y), scaled_size) / scaled_size * GRID_SIZE;
          if (x < 0.0) {
              STATE(self, GRIDMAP_X) = GRID_SIZE - STATE(self, GRIDMAP_X);
          }
          if (y < 0.0) {
              STATE(self, GRIDMAP_Y) = GRID_SIZE - STATE(self, GRIDMAP_Y);
          }
      }
  
      float base_radius = expf(BASEVAL(self, RADIUS_LOGARITHMIC));
      STATE(self, BARREL_ROTATION) += step_barrel_rotation;
  
      // FIXME: does happen (interpolation problem?)
      if (STATE(self, PRESSURE) <= 0.0) STATE(self, PRESSURE) = 0.0;
      const float pressure = STATE(self, PRESSURE);
  
      { // start / end stroke (for "stroke" input only)
        const float lim = 0.0001;
        const float threshold = BASEVAL(self, STROKE_THRESHOLD);
        const float started = STATE(self, STROKE_STARTED);
        if (!started && pressure > threshold + lim) {
          // start new stroke
          STATE(self, STROKE_STARTED) = 1;
          STATE(self, STROKE) = 0.0;
        } else if (started && pressure <= threshold * 0.9 + lim) {
          // end stroke
          STATE(self, STROKE_STARTED) = 0;
        }
      }
  
      // now follows input handling
  
      //adjust speed with viewzoom
      const float norm_dx = step_dx / step_dtime * STATE(self, VIEWZOOM);
      const float norm_dy = step_dy / step_dtime * STATE(self, VIEWZOOM);
  
      const float norm_speed = hypotf(norm_dx, norm_dy);
      //norm_dist should relate to brush size, whereas norm_speed should not
      const float norm_dist = hypotf(step_dx / step_dtime / base_radius, step_dy / step_dtime / base_radius) * step_dtime;
  
      float inputs[MYPAINT_BRUSH_INPUTS_COUNT];
  
      INPUT(PRESSURE) = pressure * expf(BASEVAL(self, PRESSURE_GAIN_LOG));
  
      const float m0 = self->speed_mapping_m[0];
      const float q0 = self->speed_mapping_q[0];
      const float m1 = self->speed_mapping_m[1];
      const float q1 = self->speed_mapping_q[1];
      INPUT(SPEED1) = log(self->speed_mapping_gamma[0] + STATE(self, NORM_SPEED1_SLOW)) * m0 + q0;
      INPUT(SPEED2) = log(self->speed_mapping_gamma[1] + STATE(self, NORM_SPEED2_SLOW)) * m1 + q1;
  
      INPUT(RANDOM) = self->random_input;
      INPUT(STROKE) = MIN(STATE(self, STROKE), 1.0);
  
      //correct direction for varying view rotation
      const float dir_angle = atan2f(STATE(self, DIRECTION_DY), STATE(self, DIRECTION_DX));
      INPUT(DIRECTION) = mod_arith(DEGREES(dir_angle) + viewrotation + 180.0, 180.0);
      const float dir_angle_360 = atan2f(STATE(self, DIRECTION_ANGLE_DY), STATE(self, DIRECTION_ANGLE_DX));
      INPUT(DIRECTION_ANGLE) = fmodf(DEGREES(dir_angle_360) + viewrotation + 360.0, 360.0) ;
      INPUT(TILT_DECLINATION) = STATE(self, DECLINATION);
      //correct ascension for varying view rotation, use custom mod
      INPUT(TILT_ASCENSION) = mod_arith(STATE(self, ASCENSION) + viewrotation + 180.0, 360.0) - 180.0;
      INPUT(VIEWZOOM) = BASEVAL(self, RADIUS_LOGARITHMIC) - logf(base_radius / STATE(self, VIEWZOOM));
      INPUT(ATTACK_ANGLE) = smallest_angular_difference(STATE(self, ASCENSION), mod_arith(DEGREES(dir_angle_360) + 90, 360));
      INPUT(BRUSH_RADIUS) = BASEVAL(self, RADIUS_LOGARITHMIC);
  
      INPUT(GRIDMAP_X) = CLAMP(STATE(self, GRIDMAP_X), 0.0, GRID_SIZE);
      INPUT(GRIDMAP_Y) = CLAMP(STATE(self, GRIDMAP_Y), 0.0, GRID_SIZE);
  
      INPUT(TILT_DECLINATIONX) = STATE(self, DECLINATIONX);
      INPUT(TILT_DECLINATIONY) = STATE(self, DECLINATIONY);
  
      INPUT(CUSTOM) = STATE(self, CUSTOM_INPUT);
      INPUT(BARREL_ROTATION) = mod_arith(STATE(self, BARREL_ROTATION), 360);
  
      if (self->print_inputs) {
          print_inputs(self, inputs);
      }
  
      for (int i = 0; i < MYPAINT_BRUSH_SETTINGS_COUNT; i++) {
        self->settings_value[i] = mypaint_mapping_calculate(self->settings[i], (inputs));
      }
  
      STATE(self, DABS_PER_BASIC_RADIUS) = SETTING(self, DABS_PER_BASIC_RADIUS);
      STATE(self, DABS_PER_ACTUAL_RADIUS) = SETTING(self, DABS_PER_ACTUAL_RADIUS);
      STATE(self, DABS_PER_SECOND) = SETTING(self, DABS_PER_SECOND);
  
      {
        const float fac = 1.0 - exp_decay(SETTING(self, SLOW_TRACKING_PER_DAB), step_ddab);
        STATE(self, ACTUAL_X) += (STATE(self, X) - STATE(self, ACTUAL_X)) * fac;
        STATE(self, ACTUAL_Y) += (STATE(self, Y) - STATE(self, ACTUAL_Y)) * fac;
      }
  
      { // slow speed
        const float fac1 = 1.0 - exp_decay(SETTING(self, SPEED1_SLOWNESS), step_dtime);
        STATE(self, NORM_SPEED1_SLOW) += (norm_speed - STATE(self, NORM_SPEED1_SLOW)) * fac1;
        const float fac2 = 1.0 - exp_decay (SETTING(self, SPEED2_SLOWNESS), step_dtime);
        STATE(self, NORM_SPEED2_SLOW) += (norm_speed - STATE(self, NORM_SPEED2_SLOW)) * fac2;
      }
  
      { // slow speed, but as vector this time
        float time_constant = expf(SETTING(self, OFFSET_BY_SPEED_SLOWNESS)*0.01)-1.0;
        // Workaround for a bug that happens mainly on Windows, causing
        // individual dabs to be placed far far away. Using the speed
        // with zero filtering is just asking for trouble anyway.
        if (time_constant < 0.002) time_constant = 0.002;
        const float fac = 1.0 - exp_decay (time_constant, step_dtime);
        STATE(self, NORM_DX_SLOW) += (norm_dx - STATE(self, NORM_DX_SLOW)) * fac;
        STATE(self, NORM_DY_SLOW) += (norm_dy - STATE(self, NORM_DY_SLOW)) * fac;
      }
  
      { // orientation (similar lowpass filter as above, but use dabtime instead of wallclock time)
        // adjust speed with viewzoom
        float dx = step_dx * STATE(self, VIEWZOOM);
        float dy = step_dy * STATE(self, VIEWZOOM);
  
        const float step_in_dabtime = hypotf(dx, dy);
        const float fac = 1.0 - exp_decay(expf(SETTING(self, DIRECTION_FILTER) * 0.5) - 1.0, step_in_dabtime);
  
        const float dx_old = STATE(self, DIRECTION_DX);
        const float dy_old = STATE(self, DIRECTION_DY);
  
        // 360 Direction
        STATE(self, DIRECTION_ANGLE_DX) += (dx - STATE(self, DIRECTION_ANGLE_DX)) * fac;
        STATE(self, DIRECTION_ANGLE_DY) += (dy - STATE(self, DIRECTION_ANGLE_DY)) * fac;
  
        // use the opposite speed vector if it is closer (we don't care about 180 degree turns)
        if (SQR(dx_old-dx) + SQR(dy_old-dy) > SQR(dx_old-(-dx)) + SQR(dy_old-(-dy))) {
          dx = -dx;
          dy = -dy;
        }
        STATE(self, DIRECTION_DX) += (dx - STATE(self, DIRECTION_DX)) * fac;
        STATE(self, DIRECTION_DY) += (dy - STATE(self, DIRECTION_DY)) * fac;
      }
  
      { // custom input
        const float fac = 1.0 - exp_decay (SETTING(self, CUSTOM_INPUT_SLOWNESS), 0.1);
        STATE(self, CUSTOM_INPUT) += (SETTING(self, CUSTOM_INPUT) - STATE(self, CUSTOM_INPUT)) * fac;
      }
  
      { // stroke length
        const float frequency = expf(-SETTING(self, STROKE_DURATION_LOGARITHMIC));
        const float stroke = MAX(0, STATE(self, STROKE) + norm_dist * frequency);
        const float wrap = 1.0 + MAX(0, SETTING(self, STROKE_HOLDTIME));
        // If the hold time is above 9.9, it is considered infinite, and if the stroke value has reached
        // that threshold it is no longer updated (until the stroke is reset, or the hold time changes).
        if (stroke >= wrap && wrap > 9.9 + 1.0) {
          STATE(self, STROKE) = 1.0;
        } else if (stroke >= wrap) {
          STATE(self, STROKE) = fmodf(stroke, wrap);
        } else {
          STATE(self, STROKE) = stroke;
        }
      }
  
      // calculate final radius
      const float radius_log = SETTING(self, RADIUS_LOGARITHMIC);
      STATE(self, ACTUAL_RADIUS) = expf(radius_log);
      if (STATE(self, ACTUAL_RADIUS) < ACTUAL_RADIUS_MIN) STATE(self, ACTUAL_RADIUS) = ACTUAL_RADIUS_MIN;
      if (STATE(self, ACTUAL_RADIUS) > ACTUAL_RADIUS_MAX) STATE(self, ACTUAL_RADIUS) = ACTUAL_RADIUS_MAX;
  
      // aspect ratio (needs to be calculated here because it can affect the dab spacing)
      STATE(self, ACTUAL_ELLIPTICAL_DAB_RATIO) = SETTING(self, ELLIPTICAL_DAB_RATIO);
      // correct dab angle for view rotation
      STATE(self, ACTUAL_ELLIPTICAL_DAB_ANGLE) = mod_arith(SETTING(self, ELLIPTICAL_DAB_ANGLE) - viewrotation + 180.0, 180.0) - 180.0;
    }
  
    float *fetch_smudge_bucket(MyPaintBrush *self) {
      if (!self->smudge_buckets || !self->num_buckets) {
        return &STATE(self, SMUDGE_RA);
      }
      const int bucket_index = CLAMP(roundf(SETTING(self, SMUDGE_BUCKET)), 0, self->num_buckets - 1);
      if (self->min_bucket_used == -1 || self->min_bucket_used > bucket_index) {
        self->min_bucket_used = bucket_index;
      }
      if (self->max_bucket_used < bucket_index) {
        self->max_bucket_used = bucket_index;
      }
      return &self->smudge_buckets[bucket_index * SMUDGE_BUCKET_SIZE];
    }
  
    gboolean
    update_smudge_color(
        const MyPaintBrush* self, MyPaintSurface* surface, float* const smudge_bucket, const float smudge_length, int px,
        int py, const float radius, const float legacy_smudge, const float paint_factor)
    {
  
        // Value between 0.01 and 1.0 that determines how often the canvas should be resampled
        float update_factor = MAX(0.01, smudge_length);
  
  
        // Calling get_color() is almost as expensive as rendering a
        // dab. Because of this we use the previous value if it is not
        // expected to hurt quality too much. We call it at most every
        // second dab.
        float r, g, b, a;
        const float smudge_length_log = SETTING(self, SMUDGE_LENGTH_LOG);
  
        const float recentness = smudge_bucket[PREV_COL_RECENTNESS] * update_factor;
        smudge_bucket[PREV_COL_RECENTNESS] = recentness;
  
        const float margin = 0.0000000000000001;
        if (recentness < MIN(1.0, powf(0.5 * update_factor, smudge_length_log) + margin)) {
            if (recentness == 0.0) {
                // first initialization of smudge color (initiate with color sampled from canvas)
                update_factor = 0.0;
            }
            smudge_bucket[PREV_COL_RECENTNESS] = 1.0;
  
            const float radius_log = SETTING(self, SMUDGE_RADIUS_LOG);
            const float smudge_radius = CLAMP(radius * expf(radius_log), ACTUAL_RADIUS_MIN, ACTUAL_RADIUS_MAX);
  
            // Sample colors on the canvas, using a negative value for the paint factor
            // means that the old sampling method is used, instead of weighted spectral.
            mypaint_surface_get_color(
                surface, px, py, smudge_radius, &r, &g, &b, &a, legacy_smudge ? -1.0 : paint_factor);
  
            // don't draw unless the picked-up alpha is above a certain level
            // this is sort of like lock_alpha but for smudge
            // negative values reverse this idea
            const float smudge_op_lim = SETTING(self, SMUDGE_TRANSPARENCY);
            if ((smudge_op_lim > 0.0 && a < smudge_op_lim) || (smudge_op_lim < 0.0 && a > -smudge_op_lim)) {
                return TRUE; // signals the caller to return early
            }
            smudge_bucket[PREV_COL_R] = r;
            smudge_bucket[PREV_COL_G] = g;
            smudge_bucket[PREV_COL_B] = b;
            smudge_bucket[PREV_COL_A] = a;
        } else {
            r = smudge_bucket[PREV_COL_R];
            g = smudge_bucket[PREV_COL_G];
            b = smudge_bucket[PREV_COL_B];
            a = smudge_bucket[PREV_COL_A];
        }
  
        if (legacy_smudge) {
            const float fac_old = update_factor;
            const float fac_new = (1.0 - update_factor) * a;
            smudge_bucket[SMUDGE_R] = fac_old * smudge_bucket[SMUDGE_R] + fac_new * r;
            smudge_bucket[SMUDGE_G] = fac_old * smudge_bucket[SMUDGE_G] + fac_new * g;
            smudge_bucket[SMUDGE_B] = fac_old * smudge_bucket[SMUDGE_B] + fac_new * b;
            smudge_bucket[SMUDGE_A] = CLAMP((fac_old * smudge_bucket[SMUDGE_A] + fac_new), 0.0, 1.0);
        } else if (a > WGM_EPSILON * 10) {
            float prev_smudge_color[4] = {smudge_bucket[SMUDGE_R], smudge_bucket[SMUDGE_G], smudge_bucket[SMUDGE_B],
                                          smudge_bucket[SMUDGE_A]};
            float sampled_color[4] = {r, g, b, a};
  
            float* smudge_new = mix_colors(prev_smudge_color, sampled_color, update_factor, paint_factor);
            smudge_bucket[SMUDGE_R] = smudge_new[SMUDGE_R];
            smudge_bucket[SMUDGE_G] = smudge_new[SMUDGE_G];
            smudge_bucket[SMUDGE_B] = smudge_new[SMUDGE_B];
            smudge_bucket[SMUDGE_A] = smudge_new[SMUDGE_A];
        } else {
            // To avoid color noise from spectral mixing with a low alpha,
            // we'll just decrease the alpha of the existing smudge color.
            smudge_bucket[SMUDGE_A] = (smudge_bucket[SMUDGE_A] + a) / 2;
        }
        return FALSE; // signals the caller to not return early (the default)
    }
  
    float
    apply_smudge(
        const float* const smudge_bucket, const float smudge_value, const gboolean legacy_smudge,
        const float paint_factor, float* color_r, float* color_g, float* color_b)
    {
        float smudge_factor = MIN(1.0, smudge_value);
  
        // If the smudge color somewhat transparent, then the resulting
        // dab will do erasing towards that transparency level.
        // see also ../doc/smudge_math.png
        const float eraser_target_alpha =
            CLAMP((1.0 - smudge_factor) + smudge_factor * smudge_bucket[SMUDGE_A], 0.0, 1.0);
  
        if (eraser_target_alpha > 0) {
            if (legacy_smudge) {
                const float col_factor = 1.0 - smudge_factor;
                *color_r = (smudge_factor * smudge_bucket[SMUDGE_R] + col_factor * *color_r) / eraser_target_alpha;
                *color_g = (smudge_factor * smudge_bucket[SMUDGE_G] + col_factor * *color_g) / eraser_target_alpha;
                *color_b = (smudge_factor * smudge_bucket[SMUDGE_B] + col_factor * *color_b) / eraser_target_alpha;
            } else {
                float smudge_color[4] = {smudge_bucket[SMUDGE_R], smudge_bucket[SMUDGE_G], smudge_bucket[SMUDGE_B],
                                         smudge_bucket[SMUDGE_A]};
                float brush_color[4] = {*color_r, *color_g, *color_b, 1.0};
                float* color_new = mix_colors(smudge_color, brush_color, smudge_factor, paint_factor);
                *color_r = color_new[SMUDGE_R];
                *color_g = color_new[SMUDGE_G];
                *color_b = color_new[SMUDGE_B];
            }
        } else {
            // we are only erasing; the color does (should) not matter
            // we set the color to a clear red to see bugs easier.
            *color_r = 1.0;
            *color_g = 0.0;
            *color_b = 0.0;
        }
        return eraser_target_alpha;
    }
  
    // Called only from stroke_to(). Calculate everything needed to
    // draw the dab, then let the surface do the actual drawing.
    //
    // This is only gets called right after update_states_and_setting_values().
    // Returns TRUE if the surface was modified.
  gboolean prepare_and_draw_dab (MyPaintBrush *self, MyPaintSurface * surface, gboolean linear)
    {
      const float opaque_fac = SETTING(self, OPAQUE_MULTIPLY);
      // ensure we don't get a positive result with two negative opaque values
      float opaque = MAX(0.0, SETTING(self, OPAQUE));
      opaque = CLAMP(opaque * opaque_fac, 0.0, 1.0);
  
      const float opaque_linearize = BASEVAL(self, OPAQUE_LINEARIZE);
  
      if (opaque_linearize) {
        // OPTIMIZE: no need to recalculate this for each dab
        float alpha, beta, alpha_dab, beta_dab;
        float dabs_per_pixel;
        // dabs_per_pixel is just estimated roughly, I didn't think hard
        // about the case when the radius changes during the stroke
        dabs_per_pixel = (STATE(self, DABS_PER_ACTUAL_RADIUS) + STATE(self, DABS_PER_BASIC_RADIUS)) * 2.0;
  
        // the correction is probably not wanted if the dabs don't overlap
        if (dabs_per_pixel < 1.0) dabs_per_pixel = 1.0;
  
        // interpret the user-setting smoothly
        dabs_per_pixel = 1.0 + opaque_linearize * (dabs_per_pixel-1.0);
  
        // see doc/brushdab_saturation.png
        //      beta = beta_dab^dabs_per_pixel
        // <==> beta_dab = beta^(1/dabs_per_pixel)
        alpha = opaque;
        beta = 1.0-alpha;
        beta_dab = powf(beta, 1.0/dabs_per_pixel);
        alpha_dab = 1.0-beta_dab;
        opaque = alpha_dab;
      }
  
      float x = STATE(self, ACTUAL_X);
      float y = STATE(self, ACTUAL_Y);
  
      float base_radius = expf(BASEVAL(self, RADIUS_LOGARITHMIC));
  
      // Directional offsets
      Offsets offs = directional_offsets(self, base_radius, (int)STATE(self, FLIP));
      x += offs.x;
      y += offs.y;
  
      const float view_zoom = STATE(self, VIEWZOOM);
      const float offset_by_speed = SETTING(self, OFFSET_BY_SPEED);
      if (offset_by_speed) {
          x += STATE(self, NORM_DX_SLOW) * offset_by_speed * 0.1 / view_zoom;
          y += STATE(self, NORM_DY_SLOW) * offset_by_speed * 0.1 / view_zoom;
      }
  
      const float offset_by_random = SETTING(self, OFFSET_BY_RANDOM);
      if (offset_by_random) {
          float amp = MAX(0.0, offset_by_random);
          x += rand_gauss (self->rng) * amp * base_radius;
          y += rand_gauss (self->rng) * amp * base_radius;
      }
  
      float radius = STATE(self, ACTUAL_RADIUS);
      const float radius_by_random = SETTING(self, RADIUS_BY_RANDOM);
      if (radius_by_random) {
          const float noise = rand_gauss(self->rng) * radius_by_random;
          float radius_log = SETTING(self, RADIUS_LOGARITHMIC) + noise;
          radius = CLAMP(expf(radius_log), ACTUAL_RADIUS_MIN, ACTUAL_RADIUS_MAX);
          float alpha_correction = SQR(STATE(self, ACTUAL_RADIUS) / radius);
          if (alpha_correction <= 1.0) {
              opaque *= alpha_correction;
          }
      }
  
      const float paint_factor = SETTING(self, PAINT_MODE);
      const gboolean paint_setting_constant = mypaint_mapping_is_constant(self->settings[MYPAINT_BRUSH_SETTING_PAINT_MODE]);
      const gboolean legacy_smudge = paint_factor <= 0.0 && paint_setting_constant;
  
      //convert to RGB here instead of later
      // color part
      float color_h = BASEVAL(self, COLOR_H);
      float color_s = BASEVAL(self, COLOR_S);
      float color_v = BASEVAL(self, COLOR_V);
      hsv_to_rgb_float (&color_h, &color_s, &color_v);
  
      // update smudge color
      const float smudge_length = SETTING(self, SMUDGE_LENGTH);
      if (smudge_length < 1.0 && // default smudge length is 0.5, so the smudge factor is checked as well
          (SETTING(self, SMUDGE) != 0.0 || !mypaint_mapping_is_constant(self->settings[MYPAINT_BRUSH_SETTING_SMUDGE]))) {
          float* const bucket = fetch_smudge_bucket(self);
          gboolean return_early = update_smudge_color(
              self, surface, bucket, smudge_length, ROUND(x), ROUND(y), radius, legacy_smudge, paint_factor);
          if (return_early) {
            return FALSE;
          }
      }
  
      float eraser_target_alpha = 1.0;
      const float smudge_value = SETTING(self, SMUDGE);
  
      if (smudge_value > 0.0) {
        float* const bucket = fetch_smudge_bucket(self);
        eraser_target_alpha =
            apply_smudge(bucket, smudge_value, legacy_smudge, paint_factor, &color_h, &color_s, &color_v);
      }
  
      // eraser
      if (SETTING(self, ERASER)) {
        eraser_target_alpha *= (1.0-SETTING(self, ERASER));
      }
  
      /*
        If the colors are stored as linear sRGB, they need to be transformed to
        the form compatible with the hsv/hsl conversion functions, and then
        transformed back after the adjustments.
      */
      gboolean using_hsv_dynamics =
          SETTING(self, CHANGE_COLOR_H) || SETTING(self, CHANGE_COLOR_HSV_S) || SETTING(self, CHANGE_COLOR_V);
      gboolean using_hsl_dynamics = SETTING(self, CHANGE_COLOR_L) || SETTING(self, CHANGE_COLOR_HSL_S);
      gboolean using_color_dynamics = using_hsv_dynamics || using_hsl_dynamics;
  
      // delinearize
      if (linear && using_color_dynamics) {
        color_h = powf(color_h, 1 / 2.2);
        color_s = powf(color_s, 1 / 2.2);
        color_v = powf(color_v, 1 / 2.2);
      }
  
      // HSV color change
      if (using_hsv_dynamics) {
          rgb_to_hsv_float(&color_h, &color_s, &color_v);
          color_h += SETTING(self, CHANGE_COLOR_H);
          color_s += color_s * color_v * SETTING(self, CHANGE_COLOR_HSV_S);
          color_v += SETTING(self, CHANGE_COLOR_V);
          hsv_to_rgb_float(&color_h, &color_s, &color_v);
      }
  
      // HSL color change
      if (using_hsl_dynamics) {
        // (calculating way too much here, can be optimized if necessary)
        // this function will CLAMP the inputs
        rgb_to_hsl_float (&color_h, &color_s, &color_v);
        color_v += SETTING(self, CHANGE_COLOR_L);
        color_s += color_s * MIN(fabsf(1.0f - color_v), fabsf(color_v)) * 2.0f
          * SETTING(self, CHANGE_COLOR_HSL_S);
        hsl_to_rgb_float (&color_h, &color_s, &color_v);
      }
  
      // linearize
      if (linear && using_color_dynamics) {
        color_h = powf(color_h, 2.2);
        color_s = powf(color_s, 2.2);
        color_v = powf(color_v, 2.2);
      }
  
      float hardness = CLAMP(SETTING(self, HARDNESS), 0.0f, 1.0f);
      float softness = CLAMP(SETTING(self, SOFTNESS), 0.0f, 1.0f);
  
      // anti-aliasing attempt (works surprisingly well for ink brushes)
      float current_fadeout_in_pixels = radius * (1.0 - hardness);
      float min_fadeout_in_pixels = SETTING(self, ANTI_ALIASING);
      if (current_fadeout_in_pixels < min_fadeout_in_pixels) {
        // need to soften the brush (decrease hardness), but keep optical radius
        // so we tune both radius and hardness, to get the desired fadeout_in_pixels
        float current_optical_radius = radius - (1.0 - hardness)*radius/2.0;
  
        // Equation 1: (new fadeout must be equal to min_fadeout)
        //   min_fadeout_in_pixels = radius_new*(1.0 - hardness_new)
        // Equation 2: (optical radius must remain unchanged)
        //   current_optical_radius = radius_new - (1.0-hardness_new)*radius_new/2.0
        //
        // Solved Equation 1 for hardness_new, using Equation 2: (thanks to mathomatic)
        float hardness_new = ((current_optical_radius - (min_fadeout_in_pixels/2.0))/(current_optical_radius + (min_fadeout_in_pixels/2.0)));
        // Using Equation 1:
        float radius_new = (min_fadeout_in_pixels/(1.0 - hardness_new));
  
        hardness = hardness_new;
        radius = radius_new;
      }
  
      // snap to pixel
      float snapToPixel = SETTING(self, SNAP_TO_PIXEL);
      if (snapToPixel > 0.0)
      {
        // linear interpolation between non-snapped and snapped
        const float snapped_x = floor(x) + 0.5;
        const float snapped_y = floor(y) + 0.5;
        x = x + (snapped_x - x) * snapToPixel;
        y = y + (snapped_y - y) * snapToPixel;
  
        float snapped_radius = roundf(radius * 2.0) / 2.0;
        if (snapped_radius < 0.5)
          snapped_radius = 0.5;
  
        if (snapToPixel > 0.9999 )
        {
          snapped_radius -= 0.0001; // this fixes precision issues where
                                    // neighbor pixels could be wrongly painted
        }
  
        radius = radius + (snapped_radius - radius) * snapToPixel;
      }
  
      const float dab_ratio = STATE(self, ACTUAL_ELLIPTICAL_DAB_RATIO);
      const float dab_angle = STATE(self, ACTUAL_ELLIPTICAL_DAB_ANGLE);
      const float lock_alpha = SETTING(self, LOCK_ALPHA);
      const float colorize = SETTING(self, COLORIZE);
      const float posterize = SETTING(self, POSTERIZE);
      const float posterize_num = SETTING(self, POSTERIZE_NUM);
  
      return mypaint_surface_draw_dab (
          surface, x, y, radius, color_h, color_s, color_v, opaque, hardness, softness, eraser_target_alpha,
          dab_ratio, dab_angle, lock_alpha, colorize, posterize, posterize_num, paint_factor);
    }
  
    // How many dabs will be drawn between the current and the next (x, y, +dt) position?
    float count_dabs_to (MyPaintBrush *self, float x, float y, float dt)
    {
      const float base_radius_log = BASEVAL(self, RADIUS_LOGARITHMIC);
      const float base_radius = CLAMP(expf(base_radius_log), ACTUAL_RADIUS_MIN, ACTUAL_RADIUS_MAX);
  
      if (STATE(self, ACTUAL_RADIUS) == 0.0) {
        STATE(self, ACTUAL_RADIUS) = base_radius;
      }
  
      const float dx = x - STATE(self, X);
      const float dy = y - STATE(self, Y);
  
      float dist;
  
      if (STATE(self, ACTUAL_ELLIPTICAL_DAB_RATIO) > 1.0) {
        // code duplication, see calculate_rr in mypaint-tiled-surface.c
        float angle_rad = RADIANS(STATE(self, ACTUAL_ELLIPTICAL_DAB_ANGLE));
        float cs = cos(angle_rad);
        float sn = sin(angle_rad);
        float yyr = (dy * cs - dx * sn) * STATE(self, ACTUAL_ELLIPTICAL_DAB_RATIO);
        float xxr = dy * sn + dx * cs;
        dist = sqrt(yyr * yyr + xxr * xxr);
      } else {
        dist = hypotf(dx, dy);
      }
  
      const float res1 = dist / STATE(self, ACTUAL_RADIUS) * STATE(self, DABS_PER_ACTUAL_RADIUS);
      const float res2 = dist / base_radius * STATE(self, DABS_PER_BASIC_RADIUS);
      const float res3 = dt * STATE(self, DABS_PER_SECOND);
      //on first load if isnan the engine messes up and won't paint
      //until you switch modes
      float res4 = res1 + res2 + res3;
      if (isnan(res4) || res4 < 0.0) { res4 = 0.0; }
      return res4;
    }
  
    /**
     * mypaint_brush_stroke_to:
     * @dtime: Time since last motion event, in seconds.
     * @viewzoom: Canvas zoom; 1.0 = 100% zoom. Zoom value v *must* be in range:
     * 0.0 < v < FLOAT_MAX (reasonable max is probably always below 100).
     *
     * Should be called once for each motion event.
     *
     *
     * Returns: non-0 if the stroke is finished or empty, else 0.
     */
    int mypaint_brush_stroke_to (MyPaintBrush *self, MyPaintSurface *surface,
                                  float x, float y, float pressure,
                                 float xtilt, float ytilt, double dtime, float viewzoom, float viewrotation, float barrel_rotation, gboolean linear)
    {
      const float max_dtime = 5;
  
      float tilt_ascension = 0.0;
      float tilt_declination = 90.0;
      float tilt_declinationx = 90.0;
      float tilt_declinationy = 90.0;
      if (xtilt != 0 || ytilt != 0) {
        // shield us from insane tilt input
        xtilt = CLAMP(xtilt, -1.0, 1.0);
        ytilt = CLAMP(ytilt, -1.0, 1.0);
        assert(isfinite(xtilt) && isfinite(ytilt));
  
        tilt_ascension = DEGREES(atan2(-xtilt, ytilt));
        const float rad = hypot(xtilt, ytilt);
        tilt_declination = 90-(rad*60);
        tilt_declinationx = (xtilt * 60);
        tilt_declinationy = (ytilt * 60);
  
        assert(isfinite(tilt_ascension));
        assert(isfinite(tilt_declination));
        assert(isfinite(tilt_declinationx));
        assert(isfinite(tilt_declinationy));
      }
  
      if (pressure <= 0.0) pressure = 0.0;
      if (!isfinite(x) || !isfinite(y) ||
          (x > 1e10 || y > 1e10 || x < -1e10 || y < -1e10)) {
        // workaround attempt for https://gna.org/bugs/?14372
        printf("Warning: ignoring brush::stroke_to with insane inputs (x = %f, y = %f)\n", (double)x, (double)y);
        x = 0.0;
        y = 0.0;
        pressure = 0.0;
        viewzoom = 0.0;
        viewrotation = 0.0;
        barrel_rotation = 0.0;
      }
      // the assertion below is better than out-of-memory later at save time
      assert(x < 1e8 && y < 1e8 && x > -1e8 && y > -1e8);
  
      if (dtime < 0) printf("Time jumped backwards by dtime=%f seconds!\n", dtime);
      if (dtime <= 0) dtime = 0.0001; // protect against possible division by zero bugs
  
      if (dtime > 0.100 && pressure && STATE(self, PRESSURE) == 0) {
        // Workaround for tablets that don't report motion events without pressure.
        // This is to avoid linear interpolation of the pressure between two events.
        mypaint_brush_stroke_to (self, surface, x, y, 0.0, 90.0, 0.0, dtime-0.0001, viewzoom, viewrotation, 0.0, linear);
        dtime = 0.0001;
      }
  
      // skip some length of input if requested (for stable tracking noise)
      if (self->skip > 0.001) {
        float dist = hypotf(self->skip_last_x-x, self->skip_last_y-y);
        self->skip_last_x = x;
        self->skip_last_y = y;
        self->skipped_dtime += dtime;
        self->skip -= dist;
        dtime = self->skipped_dtime;
  
        if (self->skip > 0.001 && !(dtime > max_dtime || self->reset_requested))
          return FALSE;
  
        // skipped
        self->skip = 0;
        self->skip_last_x = 0;
        self->skip_last_y = 0;
        self->skipped_dtime = 0;
      }
  
      { // calculate the actual "virtual" cursor position
  
        // noise first
        if (BASEVAL(self, TRACKING_NOISE)) {
          // OPTIMIZE: expf() called too often
          const float base_radius = expf(BASEVAL(self, RADIUS_LOGARITHMIC));
          const float noise = base_radius * BASEVAL(self, TRACKING_NOISE);
  
          if (noise > 0.001) {
            // we need to skip some length of input to make
            // tracking noise independent from input frequency
            self->skip = 0.5*noise;
            self->skip_last_x = x;
            self->skip_last_y = y;
  
            // add noise
            x += noise * rand_gauss(self->rng);
            y += noise * rand_gauss(self->rng);
          }
        }
  
        const float fac = 1.0 - exp_decay(BASEVAL(self, SLOW_TRACKING), 100.0 * dtime);
        x = STATE(self, X) + (x - STATE(self, X)) * fac;
        y = STATE(self, Y) + (y - STATE(self, Y)) * fac;
      }
  
      if (dtime > max_dtime || self->reset_requested) {
        self->reset_requested = FALSE;
  
        brush_reset(self);
  
        // reset value of random input
        self->random_input = rng_double_next(self->rng);
  
        STATE(self, X) = x;
        STATE(self, Y) = y;
        STATE(self, PRESSURE) = pressure;
  
        // not resetting, because they will get overwritten below:
        //dx, dy, dpress, dtime
  
        STATE(self, ACTUAL_X) = STATE(self, X);
        STATE(self, ACTUAL_Y) = STATE(self, Y);
        STATE(self, STROKE) = 1.0; // start in a state as if the stroke was long finished
  
        return TRUE;
      }
  
      enum { UNKNOWN, YES, NO } painted = UNKNOWN;
      double dtime_left = dtime;
  
      float step_ddab, step_dx, step_dy, step_dpressure, step_dtime;
      float step_declination, step_ascension, step_declinationx, step_declinationy, step_barrel_rotation;
  
      // draw many (or zero) dabs to the next position
      // see doc/images/stroke2dabs.png
      float dabs_moved = STATE(self, PARTIAL_DABS);
      float dabs_todo = count_dabs_to (self, x, y, dtime);
      while (dabs_moved + dabs_todo >= 1.0) { // there are dabs pending
        { // linear interpolation (nonlinear variant was too slow, see SVN log)
          float frac; // fraction of the remaining distance to move
          if (dabs_moved > 0) {
            // "move" the brush exactly to the first dab
            step_ddab = 1.0 - dabs_moved; // the step "moves" the brush by a fraction of one dab
            dabs_moved = 0;
          } else {
            step_ddab = 1.0; // the step "moves" the brush by exactly one dab
          }
          frac = step_ddab / dabs_todo;
          step_dx        = frac * (x - STATE(self, X));
          step_dy        = frac * (y - STATE(self, Y));
          step_dpressure = frac * (pressure - STATE(self, PRESSURE));
          step_dtime     = frac * (dtime_left - 0.0);
          // Though it looks different, time is interpolated exactly like x/y/pressure.
          step_declination = frac * (tilt_declination - STATE(self, DECLINATION));
          step_declinationx = frac * (tilt_declinationx - STATE(self, DECLINATIONX));
          step_declinationy = frac * (tilt_declinationy - STATE(self, DECLINATIONY));
          step_ascension   = frac * smallest_angular_difference(STATE(self, ASCENSION), tilt_ascension);
          //converts barrel_ration to degrees,
          step_barrel_rotation = frac * smallest_angular_difference(STATE(self, BARREL_ROTATION),barrel_rotation * 360);
  
          update_states_and_setting_values (self, step_ddab, step_dx, step_dy,
                                            step_dpressure, step_declination,
                                            step_ascension, step_dtime, viewzoom,
                                            viewrotation, step_declinationx,
                                            step_declinationy, step_barrel_rotation);
        }
  
        // Flips between 1 and -1, used for "mirrored" offsets.
        STATE(self, FLIP) *= -1;
        gboolean painted_now = prepare_and_draw_dab (self, surface, linear);
        if (painted_now) {
          painted = YES;
        } else if (painted == UNKNOWN) {
          painted = NO;
        }
  
        // update value of random input only when drawing the dab
        self->random_input = rng_double_next(self->rng);
  
        dtime_left -= step_dtime;
        dabs_todo = count_dabs_to(self, x, y, dtime_left);
      }
  
      {
        // "move" the brush to the current time (no more dab will happen)
        // Important to do this at least once every event, because
        // brush_count_dabs_to depends on the radius and the radius can
        // depend on something that changes much faster than just every
        // dab.
        step_ddab = dabs_todo; // the step "moves" the brush by a fraction of one dab
        step_dx        = x - STATE(self, X);
        step_dy        = y - STATE(self, Y);
        step_dpressure = pressure - STATE(self, PRESSURE);
        step_declination = tilt_declination - STATE(self, DECLINATION);
        step_declinationx = tilt_declinationx - STATE(self, DECLINATIONX);
        step_declinationy = tilt_declinationy - STATE(self, DECLINATIONY);
        step_ascension = smallest_angular_difference(STATE(self, ASCENSION), tilt_ascension);
        step_dtime     = dtime_left;
        step_barrel_rotation = smallest_angular_difference(STATE(self, BARREL_ROTATION), barrel_rotation * 360);
  
        update_states_and_setting_values (self, step_ddab, step_dx, step_dy, step_dpressure, step_declination, step_ascension, step_dtime, viewzoom, viewrotation, step_declinationx, step_declinationy, step_barrel_rotation);
      }
  
      // save the fraction of a dab that is already done now
      STATE(self, PARTIAL_DABS) = dabs_moved + dabs_todo;
  
      /* not working any more with the new rng...
      // next seed for the RNG (GRand has no get_state() and states[] must always contain our full state)
      STATE(self, RNG_SEED) = rng_double_next(self->rng);
      */
  
      // stroke separation logic (for undo/redo)
  
      if (painted == UNKNOWN) {
        if (self->stroke_current_idling_time > 0 || self->stroke_total_painting_time == 0) {
          // still idling
          painted = NO;
        } else {
          // probably still painting (we get more events than brushdabs)
          painted = YES;
          //if (pressure == 0) g_print ("info: assuming 'still painting' while there is no pressure\n");
        }
      }
      if (painted == YES) {
        //if (stroke_current_idling_time > 0) g_print ("idling ==> painting\n");
        self->stroke_total_painting_time += dtime;
        self->stroke_current_idling_time = 0;
        // force a stroke split after some time
        if (self->stroke_total_painting_time > 4 + 3*pressure) {
          // but only if pressure is not being released
          // FIXME: use some smoothed state for dpressure, not the output of the interpolation code
          //        (which might easily wrongly give dpressure == 0)
          if (step_dpressure >= 0) {
            return TRUE;
          }
        }
      } else if (painted == NO) {
        //if (stroke_current_idling_time == 0) g_print ("painting ==> idling\n");
        self->stroke_current_idling_time += dtime;
        if (self->stroke_total_painting_time == 0) {
          // not yet painted, start a new stroke if we have accumulated a lot of irrelevant motion events
          if (self->stroke_current_idling_time > 1.0) {
            return TRUE;
          }
        } else {
          // Usually we have pressure==0 here. But some brushes can paint
          // nothing at full pressure (eg gappy lines, or a stroke that
          // fades out). In either case this is the preferred moment to split.
          if (self->stroke_total_painting_time+self->stroke_current_idling_time > 0.9 + 5*pressure) {
            return TRUE;
          }
        }
      }
      return FALSE;
    }
  
  // Compat wrapper, for supporting libjson
  static gboolean
  obj_get(json_object *self, const gchar *key, json_object **obj_out) {
  #if JSON_C_MINOR_VERSION >= 10
      return json_object_object_get_ex(self, key, obj_out)
          && (obj_out ? (*obj_out != NULL) : TRUE);
  #else
      json_object *o = json_object_object_get(self, key);
      if (obj_out) {
          *obj_out = o;
      }
      return (o != NULL);
  #endif
  }
  
  static gboolean
  update_brush_setting_from_json_object(MyPaintBrush *self,
                                        char *setting_name,
                                        json_object *setting_obj)
  {
      MyPaintBrushSetting setting_id = mypaint_brush_setting_from_cname(setting_name);
  
      if (setting_id >= MYPAINT_BRUSH_SETTINGS_COUNT) {
          fprintf(stderr, "Warning: Unknown setting_id: %d for setting: %s\n",
                  setting_id, setting_name);
          return FALSE;
      }
  
      if (!json_object_is_type(setting_obj, json_type_object)) {
          fprintf(stderr, "Warning: Wrong type for setting: %s\n", setting_name);
          return FALSE;
      }
  
      // Base value
      json_object *base_value_obj = NULL;
      if (! obj_get(setting_obj, "base_value", &base_value_obj)) {
          fprintf(stderr, "Warning: No 'base_value' field for setting: %s\n", setting_name);
          return FALSE;
      }
      const double base_value = json_object_get_double(base_value_obj);
      mypaint_brush_set_base_value(self, setting_id, base_value);
  
      // Inputs
      json_object *inputs = NULL;
      if (! obj_get(setting_obj, "inputs", &inputs)) {
          fprintf(stderr, "Warning: No 'inputs' field for setting: %s\n", setting_name);
          return FALSE;
      }
      json_object_object_foreach(inputs, input_name, input_obj) {
          MyPaintBrushInput input_id = mypaint_brush_input_from_cname(input_name);
  
          if (!json_object_is_type(input_obj, json_type_array)) {
              fprintf(stderr, "Warning: Wrong inputs type for setting: %s\n", setting_name);
              return FALSE;
          }
  
          if (input_id >= MYPAINT_BRUSH_INPUTS_COUNT) {
              fprintf(stderr, "Warning: Unknown input_id: %d for input: %s\n",
                      input_id, input_name);
              return FALSE;
          }
  
          const int number_of_mapping_points = json_object_array_length(input_obj);
  
          mypaint_brush_set_mapping_n(self, setting_id, input_id, number_of_mapping_points);
  
          for (int i=0; i<number_of_mapping_points; i++) {
              json_object *mapping_point = json_object_array_get_idx(input_obj, i);
  
              json_object *x_obj = json_object_array_get_idx(mapping_point, 0);
              const float x = json_object_get_double(x_obj);
              json_object *y_obj = json_object_array_get_idx(mapping_point, 1);
              const float y = json_object_get_double(y_obj);
  
              mypaint_brush_set_mapping_point(self, setting_id, input_id, i, x, y);
          }
      }
  
      return TRUE;
  }
  
  static gboolean
  update_brush_from_json_object(MyPaintBrush *self)
  {
      // Check version
      json_object *version_object = NULL;
      if (! obj_get(self->brush_json, "version", &version_object)) {
          fprintf(stderr, "Error: No 'version' field for brush\n");
          return FALSE;
      }
      const int version = json_object_get_int(version_object);
      if (version != 3) {
          fprintf(stderr, "Error: Unsupported brush setting version: %d\n", version);
          return FALSE;
      }
  
      // Set settings
      json_object *settings = NULL;
      if (! obj_get(self->brush_json, "settings", &settings)) {
          fprintf(stderr, "Error: No 'settings' field for brush\n");
          return FALSE;
      }
  
      gboolean updated_any = FALSE;
      json_object_object_foreach(settings, setting_name, setting_obj) {
        updated_any |= update_brush_setting_from_json_object(self, setting_name, setting_obj);
      }
      return updated_any;
  }
  
  gboolean
  mypaint_brush_from_string(MyPaintBrush *self, const char *string)
  {
      json_object *brush_json = NULL;
  
      if (self->brush_json) {
          // Free
          json_object_put(self->brush_json);
          self->brush_json = NULL;
      }
      if (string) {
          brush_json = json_tokener_parse(string);
      }
  
      if (brush_json) {
          self->brush_json = brush_json;
          return update_brush_from_json_object(self);
      }
      else {
          self->brush_json = json_object_new_object();
          return FALSE;
      }
  }
  
  
  void
  mypaint_brush_from_defaults(MyPaintBrush *self) {
      for (int s = 0; s < MYPAINT_BRUSH_SETTINGS_COUNT; s++) {
          for (int i = 0; i < MYPAINT_BRUSH_INPUTS_COUNT; i++) {
              mypaint_brush_set_mapping_n(self, s, i, 0);
          }
  
          const float def = mypaint_brush_setting_info(s)->def;
          mypaint_brush_set_base_value(self, s, def);
      }
  
      mypaint_brush_set_mapping_n(self, MYPAINT_BRUSH_SETTING_OPAQUE_MULTIPLY, MYPAINT_BRUSH_INPUT_PRESSURE, 2);
      mypaint_brush_set_mapping_point(self, MYPAINT_BRUSH_SETTING_OPAQUE_MULTIPLY, MYPAINT_BRUSH_INPUT_PRESSURE, 0, 0.0, 0.0);
      mypaint_brush_set_mapping_point(self, MYPAINT_BRUSH_SETTING_OPAQUE_MULTIPLY, MYPAINT_BRUSH_INPUT_PRESSURE, 1, 1.0, 1.0);
  }
  