<|start_filename|>mypaint-symmetry.h<|end_filename|>
#ifndef MYPAINTSYMMETRY_H
#define MYPAINTSYMMETRY_H
/* libmypaint - The MyPaint Brush Library
 * Copyright (C) 2017-2019 The MyPaint Team
 *
 * 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 "mypaint-matrix.h"
#include "mypaint-glib-compat.h"

/**
  * MyPaintSymmetryType: Enumeration of different kinds of symmetry
  *
  * Prefix = 'MYPAINT_SYMMETRY_TYPE_'
  * VERTICAL: reflection across the y-axis
  * HORIZONTAL: reflection across the x-axis
  * VERTHORZ: reflection across x-axis and y-axis, special case of SNOWFLAKE
  * ROTATIONAL: rotational symmetry by N symmetry lines around a point
  * SNOWFLAKE: rotational symmetry w. reflection across the N symmetry lines
  */
typedef enum {
    MYPAINT_SYMMETRY_TYPE_VERTICAL,
    MYPAINT_SYMMETRY_TYPE_HORIZONTAL,
    MYPAINT_SYMMETRY_TYPE_VERTHORZ,
    MYPAINT_SYMMETRY_TYPE_ROTATIONAL,
    MYPAINT_SYMMETRY_TYPE_SNOWFLAKE,
    MYPAINT_SYMMETRY_TYPES_COUNT
} MyPaintSymmetryType;


/**
  * MyPaintSymmetryState: Contains the basis for symmetry calculations
  *
  * This is used to calculate the matrices that are
  * used for the actual symmetry calculations, and to
  * determine whether the matrices need to be recalculated.
  */
typedef struct {
  MyPaintSymmetryType type;
  float center_x;
  float center_y;
  float angle;
  float num_lines;
} MyPaintSymmetryState;

/**
  * MyPaintSymmetryData: Contains data used for symmetry calculations
  *
  * Instances contain a current and pending symmetry basis, and the
  * matrices used for the actual symmetry transforms. When the pending
  * state is modified, the "pending_changes" flag should be set.
  * Matrix recalculation should not be performed during draw operations.
  */
typedef struct {
  MyPaintSymmetryState state_current;
  MyPaintSymmetryState state_pending;
  gboolean pending_changes;
  gboolean active;
  int num_symmetry_matrices;
  MyPaintTransform *symmetry_matrices;
} MyPaintSymmetryData;

void mypaint_update_symmetry_state(MyPaintSymmetryData * const symmetry_data);

MyPaintSymmetryData mypaint_default_symmetry_data();

void mypaint_symmetry_data_destroy(MyPaintSymmetryData *);

void mypaint_symmetry_set_pending(
    MyPaintSymmetryData* data, gboolean active, float center_x, float center_y,
    float symmetry_angle, MyPaintSymmetryType symmetry_type, int rot_symmetry_lines);

#endif


<|start_filename|>fastapprox/fastlambertw.h<|end_filename|>
/*=====================================================================*
 *                   Copyright (C) 2011 <NAME>                   *
 * All rights reserved.                                                *
 *                                                                     *
 * Redistribution and use in source and binary forms, with             *
 * or without modification, are permitted provided that the            *
 * following conditions are met:                                       *
 *                                                                     *
 *     * Redistributions of source code must retain the                *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer.                                       *
 *                                                                     *
 *     * Redistributions in binary form must reproduce the             *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer in the documentation and/or            *
 *     other materials provided with the distribution.                 *
 *                                                                     *
 *     * Neither the name of <NAME> nor the names                *
 *     of other contributors may be used to endorse or promote         *
 *     products derived from this software without specific            *
 *     prior written permission.                                       *
 *                                                                     *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND              *
 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,         *
 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES               *
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE             *
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER               *
 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,                 *
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES            *
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE           *
 * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR                *
 * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF          *
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT           *
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY              *
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE             *
 * POSSIBILITY OF SUCH DAMAGE.                                         *
 *                                                                     *
 * Contact: <NAME> <<EMAIL>>                            *
 *=====================================================================*/

#ifndef __FAST_LAMBERT_W_H_
#define __FAST_LAMBERT_W_H_

#include <stdint.h>
#include "fastexp.h"
#include "fastlog.h"
#include "sse.h"

// these functions compute the upper branch aka W_0

static inline float
fastlambertw (float x)
{
  static const float threshold = 2.26445f;

  float c = (x < threshold) ? 1.546865557f : 1.0f;
  float d = (x < threshold) ? 2.250366841f : 0.0f;
  float a = (x < threshold) ? -0.737769969f : 0.0f;

  float logterm = fastlog (c * x + d);
  float loglogterm = fastlog (logterm);

  float minusw = -a - logterm + loglogterm - loglogterm / logterm;
  float expminusw = fastexp (minusw);
  float xexpminusw = x * expminusw;
  float pexpminusw = xexpminusw - minusw;

  return (2.0f * xexpminusw - minusw * (4.0f * xexpminusw - minusw * pexpminusw)) /
         (2.0f + pexpminusw * (2.0f - minusw));
}

static inline float
fasterlambertw (float x)
{
  static const float threshold = 2.26445f;

  float c = (x < threshold) ? 1.546865557f : 1.0f;
  float d = (x < threshold) ? 2.250366841f : 0.0f;
  float a = (x < threshold) ? -0.737769969f : 0.0f;

  float logterm = fasterlog (c * x + d);
  float loglogterm = fasterlog (logterm);

  float w = a + logterm - loglogterm + loglogterm / logterm;
  float expw = fasterexp (-w);

  return (w * w + expw * x) / (1.0f + w);
}

static inline float
fastlambertwexpx (float x)
{
  static const float k = 1.1765631309f;
  static const float a = 0.94537622168f;

  float logarg = fmaxf (x, k);
  float powarg = (x < k) ? a * (x - k) : 0;

  float logterm = fastlog (logarg);
  float powterm = fasterpow2 (powarg);  // don't need accuracy here

  float w = powterm * (logarg - logterm + logterm / logarg);
  float logw = fastlog (w);
  float p = x - logw;

  return w * (2.0f + p + w * (3.0f + 2.0f * p)) /
         (2.0f - p + w * (5.0f + 2.0f * w));
}

static inline float
fasterlambertwexpx (float x)
{
  static const float k = 1.1765631309f;
  static const float a = 0.94537622168f;

  float logarg = fmaxf (x, k);
  float powarg = (x < k) ? a * (x - k) : 0;

  float logterm = fasterlog (logarg);
  float powterm = fasterpow2 (powarg);

  float w = powterm * (logarg - logterm + logterm / logarg);
  float logw = fasterlog (w);

  return w * (1.0f + x - logw) / (1.0f + w);
}

#ifdef __SSE2__

static inline v4sf
vfastlambertw (v4sf x)
{
  const v4sf threshold = v4sfl (2.26445f);

  v4sf under = _mm_cmplt_ps (x, threshold);
  v4sf c = _mm_or_ps (_mm_and_ps (under, v4sfl (1.546865557f)),
                      _mm_andnot_ps (under, v4sfl (1.0f)));
  v4sf d = _mm_and_ps (under, v4sfl (2.250366841f));
  v4sf a = _mm_and_ps (under, v4sfl (-0.737769969f));

  v4sf logterm = vfastlog (c * x + d);
  v4sf loglogterm = vfastlog (logterm);

  v4sf minusw = -a - logterm + loglogterm - loglogterm / logterm;
  v4sf expminusw = vfastexp (minusw);
  v4sf xexpminusw = x * expminusw;
  v4sf pexpminusw = xexpminusw - minusw;

  return (v4sfl (2.0f) * xexpminusw - minusw * (v4sfl (4.0f) * xexpminusw - minusw * pexpminusw)) / 
         (v4sfl (2.0f) + pexpminusw * (v4sfl (2.0f) - minusw));
}

static inline v4sf
vfasterlambertw (v4sf x)
{
  const v4sf threshold = v4sfl (2.26445f);

  v4sf under = _mm_cmplt_ps (x, threshold);
  v4sf c = _mm_or_ps (_mm_and_ps (under, v4sfl (1.546865557f)),
                      _mm_andnot_ps (under, v4sfl (1.0f)));
  v4sf d = _mm_and_ps (under, v4sfl (2.250366841f));
  v4sf a = _mm_and_ps (under, v4sfl (-0.737769969f));

  v4sf logterm = vfasterlog (c * x + d);
  v4sf loglogterm = vfasterlog (logterm);

  v4sf w = a + logterm - loglogterm + loglogterm / logterm;
  v4sf expw = vfasterexp (-w);

  return (w * w + expw * x) / (v4sfl (1.0f) + w);
}

static inline v4sf
vfastlambertwexpx (v4sf x)
{
  const v4sf k = v4sfl (1.1765631309f);
  const v4sf a = v4sfl (0.94537622168f);
  const v4sf two = v4sfl (2.0f);
  const v4sf three = v4sfl (3.0f);
  const v4sf five = v4sfl (5.0f);

  v4sf logarg = _mm_max_ps (x, k);
  v4sf powarg = _mm_and_ps (_mm_cmplt_ps (x, k), a * (x - k));

  v4sf logterm = vfastlog (logarg);
  v4sf powterm = vfasterpow2 (powarg);  // don't need accuracy here

  v4sf w = powterm * (logarg - logterm + logterm / logarg);
  v4sf logw = vfastlog (w);
  v4sf p = x - logw;

  return w * (two + p + w * (three + two * p)) /
         (two - p + w * (five + two * w));
}

static inline v4sf
vfasterlambertwexpx (v4sf x)
{
  const v4sf k = v4sfl (1.1765631309f);
  const v4sf a = v4sfl (0.94537622168f);

  v4sf logarg = _mm_max_ps (x, k);
  v4sf powarg = _mm_and_ps (_mm_cmplt_ps (x, k), a * (x - k));

  v4sf logterm = vfasterlog (logarg);
  v4sf powterm = vfasterpow2 (powarg);

  v4sf w = powterm * (logarg - logterm + logterm / logarg);
  v4sf logw = vfasterlog (w);

  return w * (v4sfl (1.0f) + x - logw) / (v4sfl (1.0f) + w);
}

#endif // __SSE2__

#endif // __FAST_LAMBERT_W_H_


<|start_filename|>fifo.c<|end_filename|>
/*
 *
 * Based on fifo.c from syslog-win32 by <NAME>
 * http://syslog-win32.svn.sourceforge.net/viewvc/syslog-win32/tags/syslog-win32-0-3/daemon/fifo.c
 *
 * THIS SOFTWARE IS NOT COPYRIGHTED
 *
 * This source code is offered for use in the public domain. You may
 * use, modify or distribute it freely.
 *
 * This code is distributed in the hope that it will be useful but
 * WITHOUT ANY WARRANTY. ALL WARRANTIES, EXPRESS OR IMPLIED ARE HEREBY
 * DISCLAIMED. This includes but is not limited to warranties of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 *
 */

#include "config.h"

#include <stdlib.h>
#include "fifo.h"

struct fifo_item
{
    struct fifo_item *next;  /* queue is a single-linked list */
    void *payload;
};

struct fifo
{
    struct fifo_item *first;  /* first pushed item */
    struct fifo_item *last;   /* last pushed item */
    int item_count;
};

/*
 * fifo_new:
 *
 * Allocate and initialize fifo structure. Add an empty item to the fifo.
 */
struct fifo* fifo_new(void)
{
    struct fifo *ret = (struct fifo *) malloc(sizeof(struct fifo));
    ret->first = NULL;
    ret->last = NULL;
    ret->item_count = 0;
    return ret;
}

/*
 * fifo_free:
 *
 * Delete all items and free fifo structure.
 */
void fifo_free(struct fifo* queue, FifoUserFreeFunction user_free)
{
    struct fifo_item *item;

    while ((item = queue->first) != NULL)
    {
        queue->first = item->next;
        user_free(item);
    }
    free(queue);
}

/*
 * fifo_push:
 *
 * Add item to queue.
 */
void fifo_push(struct fifo* queue, void* data)
{
    struct fifo_item *item = (struct fifo_item*) malloc(sizeof(struct fifo_item));
    item->next = NULL;
    item->payload = data;
    if (!queue->last)
        queue->first = item;
    else
        queue->last->next = item;
    queue->last = item;
}

/*
 * fifo_pop:
 *
 * Extract item from queue. Returns NULL on empty queue.
 */
void* fifo_pop(struct fifo* queue)
{
    struct fifo_item *item;
    void *data;

    item = queue->first;
    if (!item)
        return NULL;

    queue->first = item->next;
    if (!queue->first)
        queue->last = NULL;

    data = item->payload;
    free(item);
    queue->item_count--;
    return data;
}

void* fifo_peek_first(struct fifo *queue) {
    return (!queue->first) ? NULL : queue->first->payload;
}

void* fifo_peek_last(struct fifo *queue) {
    return (!queue->last) ? NULL : queue->last->payload;
}


<|start_filename|>examples/minimal.c<|end_filename|>
#include "libmypaint.c"
#include "mypaint-fixed-tiled-surface.h"

void
stroke_to(MyPaintBrush *brush, MyPaintSurface *surf, float x, float y)
{
    float viewzoom = 1.0, viewrotation = 0.0, barrel_rotation = 0.0;
    float pressure = 1.0, ytilt = 0.0, xtilt = 0.0, dtime = 1.0/10;
    gboolean linear = FALSE
    mypaint_brush_stroke_to
      (brush, surf, x, y, pressure, xtilt, ytilt, dtime, viewzoom, viewrotation, barrel_rotation, linear);
}

int
main(int argc, char argv[]) {

    MyPaintBrush *brush = mypaint_brush_new();
    int w = 300;
    int h = 150;
    float wq = (float)w / 5;
    float hq = (float)h / 5;

    /* Create an instance of the simple and naive fixed_tile_surface to draw on. */
    MyPaintFixedTiledSurface *surface = mypaint_fixed_tiled_surface_new(w, h);

    /* Create a brush with default settings for all parameters, then set its color to red. */
    mypaint_brush_from_defaults(brush);
    mypaint_brush_set_base_value(brush, MYPAINT_BRUSH_SETTING_COLOR_H, 0.0);
    mypaint_brush_set_base_value(brush, MYPAINT_BRUSH_SETTING_COLOR_S, 1.0);
    mypaint_brush_set_base_value(brush, MYPAINT_BRUSH_SETTING_COLOR_V, 1.0);

    /* Draw a rectangle on the surface using the brush */
    mypaint_surface_begin_atomic((MyPaintSurface*)surface);
    stroke_to(brush, (MyPaintSurface*)surface, wq, hq);
    stroke_to(brush, (MyPaintSurface*)surface, 4 * wq, hq);
    stroke_to(brush, (MyPaintSurface*)surface, 4 * wq, 4 * hq);
    stroke_to(brush, (MyPaintSurface*)surface, wq, 4 * hq);
    stroke_to(brush, (MyPaintSurface*)surface, wq, hq);

    /* Finalize the surface operation, passing one or more invalidation
       rectangles to get information about which areas were affected by
       the operations between ``surface_begin_atomic`` and ``surface_end_atomic.``
    */
    MyPaintRectangle roi;
    MyPaintRectangles rois;
    rois.num_rectangles = 1;
    rois.rectangles = &roi;
    mypaint_surface_end_atomic((MyPaintSurface *)surface, &rois);

    /* Write the surface pixels to a ppm image file */
    fprintf(stdout, "Writing output\n");
    write_ppm(surface, "output.ppm");

    mypaint_brush_unref(brush);
    mypaint_surface_unref((MyPaintSurface *)surface);
}


<|start_filename|>helpers.c<|end_filename|>
/* libmypaint - The MyPaint Brush Library
 * Copyright (C) 2007-2008 <NAME> <<EMAIL>>
 *
 * 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.
 */

#ifndef HELPERS_C
#define HELPERS_C

#include "config.h"

#include <assert.h>
#include <stdint.h>
#include <math.h>
#include "fastapprox/fastpow.h"

#include "helpers.h"

/*const float T_MATRIX[3][36] = {{0.000578913,0.001952085,0.009886235,0.032720398,0.100474668,0.183366464,0.233267126,0.172815304,0.021160832,-0.170961409,-0.358555623,-0.487793958,-0.674399544,-0.886748322,-0.97045709,-0.872696304,-0.559560624,-0.134497482,0.395369748,0.969077244,1.563646415,1.918490755,2.226446938,2.219830783,1.916051812,1.395620385,0.990444867,0.604042138,0.353697296,0.192706913,0.098266461,0.042521122,0.021860797,0.011569942,0.004800182,0.002704537},*/
/*{-0.000491981,-0.00166858,-0.008527451,-0.028611512,-0.089589024,-0.169698855,-0.232545306,-0.211643919,-0.117700145,0.039996723,0.233957719,0.411776827,0.669587627,1.014305033,1.33449208,1.570104952,1.575060777,1.504833712,1.290156767,1.008658851,0.712494742,0.377174433,0.138783274,-0.025203917,-0.099437546,-0.104503807,-0.088552175,-0.059244144,-0.036402168,-0.020300987,-0.010518378,-0.004600355,-0.002372843,-0.001255839,-0.000521027,-0.000293559},*/
/*{0.00344444,0.011708103,0.059997877,0.202735833,0.644403983,1.282371701,1.953710066,2.206582549,2.085203252,1.555923679,0.970316553,0.491245277,0.242873415,0.070279007,-0.061461896,-0.131571399,-0.164055717,-0.176617196,-0.165922717,-0.144226781,-0.119603243,-0.085357673,-0.061960232,-0.041673571,-0.026318559,-0.01522674,-0.009013144,-0.004849822,-0.002624901,-0.001371407,-0.00067843,-0.000287423,-0.000146799,-7.76941E-05,-3.2234E-05,-1.81614E-05}};*/

/*const float spectral_r[36] = {0.001476566,0.001476571,0.001476697,0.001477428,0.001480081,0.00148896,0.001510625,0.001553174,0.001622652,0.001723776,0.001858242,0.002028291,0.002237804,0.002500589,0.002843823,0.003312684,0.003983853,0.004975295,0.0064965,0.008912057,0.012881757,0.019616497,0.031083884,0.050157206,0.079379607,0.117964522,0.16013148,0.196518985,0.222295104,0.237687584,0.245789155,0.249669549,0.251605188,0.252511963,0.252870747,0.252999473};*/

/*const float spectral_g[36] = {0.004457553,0.004459131,0.004461979,0.004471645,0.004506172,0.004629205,0.004938512,0.005602633,0.006868923,0.009209474,0.013450892,0.021153289,0.034963807,0.059025032,0.094985053,0.129917789,0.136598211,0.112473456,0.080301564,0.055184078,0.038887514,0.028984983,0.022965813,0.019364324,0.01723558,0.016001357,0.015318335,0.014937495,0.014732014,0.014630031,0.014587647,0.014573872,0.014574704,0.014576967,0.014575977,0.014579698};*/

/*const float spectral_b[36] = {0.089623258,0.08963333,0.089677437,0.089887675,0.090595886,0.092619509,0.096261517,0.099409964,0.0953724,0.077775801,0.053270305,0.03267985,0.019295409,0.011424407,0.006963547,0.004447989,0.003004361,0.002142226,0.001608547,0.001266496,0.001041612,0.00089261,0.000792505,0.000726769,0.000685047,0.000659708,0.000644691,0.000636345,0.000631858,0.000629574,0.000628471,0.000627966,0.000627715,0.000627593,0.000627548,0.000627528};*/

static const float T_MATRIX_SMALL[3][10] = {{0.026595621243689,0.049779426257903,0.022449850859496,-0.218453689278271
,-0.256894883201278,0.445881722194840,0.772365886289756,0.194498761382537
,0.014038157587820,0.007687264480513}
,{-0.032601672674412,-0.061021043498478,-0.052490001018404
,0.206659098273522,0.572496335158169,0.317837248815438,-0.021216624031211
,-0.019387668756117,-0.001521339050858,-0.000835181622534}
,{0.339475473216284,0.635401374177222,0.771520797089589,0.113222640692379
,-0.055251113343776,-0.048222578468680,-0.012966666339586
,-0.001523814504223,-0.000094718948810,-0.000051604594741}};

static const float spectral_r_small[10] = {0.009281362787953,0.009732627042016,0.011254252737167,0.015105578649573
,0.024797924177217,0.083622585502406,0.977865045723212,1.000000000000000
,0.999961046144372,0.999999992756822};

static const float spectral_g_small[10] = {0.002854127435775,0.003917589679914,0.012132151699187,0.748259205918013
,1.000000000000000,0.865695937531795,0.037477469241101,0.022816789725717
,0.021747419446456,0.021384940572308};

static const float spectral_b_small[10] = {0.537052150373386,0.546646402401469,0.575501819073983,0.258778829633924
,0.041709923751716,0.012662638828324,0.007485593127390,0.006766900622462
,0.006699764779016,0.006676219883241};


float rand_gauss (RngDouble * rng)
{
  double sum = 0.0;
  sum += rng_double_next(rng);
  sum += rng_double_next(rng);
  sum += rng_double_next(rng);
  sum += rng_double_next(rng);
  return sum * 1.73205080757 - 3.46410161514;
}

// C fmodf function is not "arithmetic modulo"; it doesn't handle negative dividends as you might expect
// if you expect 0 or a positive number when dealing with negatives, use
// this function instead.
float mod_arith(float a, float N)
{
    float ret = a - N * floor (a / N);
    return ret;
}

// Returns the smallest angular difference
float smallest_angular_difference(float angleA, float angleB)
{
    float a;
    a = angleB - angleA;
    a = mod_arith((a + 180), 360) - 180;
    a += (a>180) ? -360 : (a<-180) ? 360 : 0;
    return a;
}

// stolen from GIMP (gimpcolorspace.c)
// (from gimp_rgb_to_hsv)
void
rgb_to_hsv_float (float *r_ /*h*/, float *g_ /*s*/, float *b_ /*v*/)
{
  float max, min, delta;
  float h, s, v;
  float r, g, b;

  h = 0.0; // silence gcc warning

  r = *r_;
  g = *g_;
  b = *b_;

  r = CLAMP(r, 0.0, 1.0);
  g = CLAMP(g, 0.0, 1.0);
  b = CLAMP(b, 0.0, 1.0);

  max = MAX3(r, g, b);
  min = MIN3(r, g, b);

  v = max;
  delta = max - min;

  if (delta > 0.0001)
    {
      s = delta / max;

      if (r == max)
        {
          h = (g - b) / delta;
          if (h < 0.0)
            h += 6.0;
        }
      else if (g == max)
        {
          h = 2.0 + (b - r) / delta;
        }
      else if (b == max)
        {
          h = 4.0 + (r - g) / delta;
        }

      h /= 6.0;
    }
  else
    {
      s = 0.0;
      h = 0.0;
    }

  *r_ = h;
  *g_ = s;
  *b_ = v;
}

// (from gimp_hsv_to_rgb)
void
hsv_to_rgb_float (float *h_, float *s_, float *v_)
{
  int    i;
  double f, w, q, t;
  float h, s, v;
  float r, g, b;
  r = g = b = 0.0; // silence gcc warning

  h = *h_;
  s = *s_;
  v = *v_;

  h = h - floor(h);
  s = CLAMP(s, 0.0, 1.0);
  v = CLAMP(v, 0.0, 1.0);

  double hue;

  if (s == 0.0)
    {
      r = v;
      g = v;
      b = v;
    }
  else
    {
      hue = h;

      if (hue == 1.0)
        hue = 0.0;

      hue *= 6.0;

      i = (int) hue;
      f = hue - i;
      w = v * (1.0 - s);
      q = v * (1.0 - (s * f));
      t = v * (1.0 - (s * (1.0 - f)));

      switch (i)
        {
        case 0:
          r = v;
          g = t;
          b = w;
          break;
        case 1:
          r = q;
          g = v;
          b = w;
          break;
        case 2:
          r = w;
          g = v;
          b = t;
          break;
        case 3:
          r = w;
          g = q;
          b = v;
          break;
        case 4:
          r = t;
          g = w;
          b = v;
          break;
        case 5:
          r = v;
          g = w;
          b = q;
          break;
        }
    }

  *h_ = r;
  *s_ = g;
  *v_ = b;
}

// (from gimp_rgb_to_hsl)
void
rgb_to_hsl_float (float *r_, float *g_, float *b_)
{
  double max, min, delta;

  float h, s, l;
  float r, g, b;

  // silence gcc warnings
  h=0;

  r = *r_;
  g = *g_;
  b = *b_;

  r = CLAMP(r, 0.0, 1.0);
  g = CLAMP(g, 0.0, 1.0);
  b = CLAMP(b, 0.0, 1.0);

  max = MAX3(r, g, b);
  min = MIN3(r, g, b);

  l = (max + min) / 2.0;

  if (max == min)
    {
      s = 0.0;
      h = 0.0; //GIMP_HSL_UNDEFINED;
    }
  else
    {
      if (l <= 0.5)
        s = (max - min) / (max + min);
      else
        s = (max - min) / (2.0 - max - min);

      delta = max - min;

      if (delta == 0.0)
        delta = 1.0;

      if (r == max)
        {
          h = (g - b) / delta;
        }
      else if (g == max)
        {
          h = 2.0 + (b - r) / delta;
        }
      else if (b == max)
        {
          h = 4.0 + (r - g) / delta;
        }

      h /= 6.0;

      if (h < 0.0)
        h += 1.0;
    }

  *r_ = h;
  *g_ = s;
  *b_ = l;
}

static double
hsl_value (double n1,
           double n2,
           double hue)
{
  double val;

  if (hue > 6.0)
    hue -= 6.0;
  else if (hue < 0.0)
    hue += 6.0;

  if (hue < 1.0)
    val = n1 + (n2 - n1) * hue;
  else if (hue < 3.0)
    val = n2;
  else if (hue < 4.0)
    val = n1 + (n2 - n1) * (4.0 - hue);
  else
    val = n1;

  return val;
}


/**
 * gimp_hsl_to_rgb:
 * @hsl: A color value in the HSL colorspace
 * @rgb: The value converted to a value in the RGB colorspace
 *
 * Convert a HSL color value to an RGB color value.
 **/
void
hsl_to_rgb_float (float *h_, float *s_, float *l_)
{
  float h, s, l;
  float r, g, b;

  h = *h_;
  s = *s_;
  l = *l_;

  h = h - floor(h);
  s = CLAMP(s, 0.0, 1.0);
  l = CLAMP(l, 0.0, 1.0);

  if (s == 0)
    {
      /*  achromatic case  */
      r = l;
      g = l;
      b = l;
    }
  else
    {
      double m1, m2;

      if (l <= 0.5)
        m2 = l * (1.0 + s);
      else
        m2 = l + s - l * s;

      m1 = 2.0 * l - m2;

      r = hsl_value (m1, m2, h * 6.0 + 2.0);
      g = hsl_value (m1, m2, h * 6.0);
      b = hsl_value (m1, m2, h * 6.0 - 2.0);
    }

  *h_ = r;
  *s_ = g;
  *l_ = b;
}

void
rgb_to_hcy_float (float *r_, float *g_, float *b_) {
	
	float _HCY_RED_LUMA = 0.2162;
	float _HCY_GREEN_LUMA = 0.7152;
	float _HCY_BLUE_LUMA = 0.0722;
	float h, c, y;
	float r, g, b;
	float p, n, d;

	r = *r_;
	g = *g_;
	b = *b_;

	// Luma is just a weighted sum of the three components.
	y = _HCY_RED_LUMA*r + _HCY_GREEN_LUMA*g + _HCY_BLUE_LUMA*b;

	// Hue. First pick a sector based on the greatest RGB component, then add
	// the scaled difference of the other two RGB components.
	p = MAX3(r, g, b);
	n = MIN3(r, g, b);
	d = p - n; // An absolute measure of chroma: only used for scaling

	if (n == p){
		h = 0.0;
	} else if (p == r){
		h = (g - b)/d;
		if (h < 0){
			h += 6.0;
		}
	} else if (p == g){
		h = ((b - r)/d) + 2.0;
	} else {  // p==b
		h = ((r - g)/d) + 4.0;
	}
	h /= 6.0;
	h = fmod(h,1.0);

	// Chroma, relative to the RGB gamut envelope.
	if ((r == g) && (g == b)){
		// Avoid a division by zero for the achromatic case.
		c = 0.0;
	} else {
		// For the derivation, see the GLHS paper.
		c = MAX((y-n)/y, (p-y)/(1-y));
	}

	*r_ = h;
	*g_ = c;
	*b_ = y;
}

void
hcy_to_rgb_float (float *h_, float *c_, float *y_) {
	
	float _HCY_RED_LUMA = 0.2162;
	float _HCY_GREEN_LUMA = 0.7152;
	float _HCY_BLUE_LUMA = 0.0722;
	float h, c, y;
	float r, g, b;
	float th, tm;

	h = *h_;
	c = *c_;
	y = *y_;

	h = h - floor(h);
	c = CLAMP(c, 0.0f, 1.0f);
	y = CLAMP(y, 0.0f, 1.0f);

	if (c == 0)	{
	  /*  achromatic case  */
	  r = y;
	  g = y;
	  b = y;
	}

	h = fmod(h, 1.0);
	h *= 6.0;

	if (h < 1){
		// implies (p==r and h==(g-b)/d and g>=b)
		th = h;
		tm = _HCY_RED_LUMA + _HCY_GREEN_LUMA * th;
	} else if (h < 2) {
		// implies (p==g and h==((b-r)/d)+2.0 and b<r)
		th = 2.0 - h;
		tm = _HCY_GREEN_LUMA + _HCY_RED_LUMA * th;
	} else if (h < 3){
		// implies (p==g and h==((b-r)/d)+2.0 and b>=g)
		th = h - 2.0;
		tm = _HCY_GREEN_LUMA + _HCY_BLUE_LUMA * th;
	} else if (h < 4) {
		// implies (p==b and h==((r-g)/d)+4.0 and r<g)
		th = 4.0 - h;
		tm = _HCY_BLUE_LUMA + _HCY_GREEN_LUMA * th;
	} else if (h < 5){
		// implies (p==b and h==((r-g)/d)+4.0 and r>=g)
		th = h - 4.0;
		tm = _HCY_BLUE_LUMA + _HCY_RED_LUMA * th;
	} else {
		// implies (p==r and h==(g-b)/d and g<b)
		th = 6.0 - h;
		tm = _HCY_RED_LUMA + _HCY_BLUE_LUMA * th;
	}

	float n,p,o;
	// Calculate the RGB components in sorted order
	if (tm >= y){
		p = y + y*c*(1-tm)/tm;
		o = y + y*c*(th-tm)/tm;
		n = y - (y*c);
	}else{
		p = y + (1-y)*c;
		o = y + (1-y)*c*(th-tm)/(1-tm);
		n = y - (1-y)*c*tm/(1-tm);
	}

	// Back to RGB order
	if (h < 1){
		r = p;
		g = o;
		b = n;
	} else if (h < 2){
		r = o;
		g = p;
		b = n;
	} else if (h < 3){
		r = n;
		g = p;
		b = o;
	} else if (h < 4){
		r = n;
		g = o;
		b = p;
	} else if (h < 5){
		r = o;
		g = n;
		b = p;
	}else{ 
		r = p;
		g = n;
		b = o;
	}

	*h_ = r;
	*c_ = g;
	*y_ = b;
}


void
rgb_to_spectral (float r, float g, float b, float *spectral_) {
  float offset = 1.0 - WGM_EPSILON;
  r = r * offset + WGM_EPSILON;
  g = g * offset + WGM_EPSILON;
  b = b * offset + WGM_EPSILON;
  //upsample rgb to spectral primaries
  float spec_r[10] = {0};
  for (int i=0; i < 10; i++) {
    spec_r[i] = spectral_r_small[i] * r;
  }
  float spec_g[10] = {0};
  for (int i=0; i < 10; i++) {
    spec_g[i] = spectral_g_small[i] * g;
  }
  float spec_b[10] = {0};
  for (int i=0; i < 10; i++) {
    spec_b[i] = spectral_b_small[i] * b;
  }
  //collapse into one spd
  for (int i=0; i<10; i++) {
    spectral_[i] += spec_r[i] + spec_g[i] + spec_b[i];
  }

}

void
spectral_to_rgb (float *spectral, float *rgb_) {
  float offset = 1.0 - WGM_EPSILON;
  // We need this tmp. array to allow auto vectorization.
  float tmp[3] = {0};
  for (int i=0; i<10; i++) {
    tmp[0] += T_MATRIX_SMALL[0][i] * spectral[i];
    tmp[1] += T_MATRIX_SMALL[1][i] * spectral[i];
    tmp[2] += T_MATRIX_SMALL[2][i] * spectral[i];
  }
  for (int i=0; i<3; i++) {
    rgb_[i] = CLAMP((tmp[i] - WGM_EPSILON) / offset, 0.0f, 1.0f);
  }
}


//function to make it easy to blend two spectral colors via weighted geometric mean
//a is the current smudge state, b is the get_color or brush color
float * mix_colors(float *a, float *b, float fac, float paint_mode)
{
  static float result[4] = {0};  
  
  float opa_a = fac;
  float opa_b = 1.0-opa_a;
  result[3] = CLAMP(opa_a * a[3] + opa_b * b[3], 0.0f, 1.0f);
  // Guard against NaN from division by zero
  float sfac_a = a[3] == 0 ? 0.0 : opa_a * a[3] / (a[3] + b[3] * opa_b);
  float sfac_b = 1 - sfac_a;

  if (paint_mode > 0.0) { 
    float spec_a[10] = {0};
    float spec_b[10] = {0};

    rgb_to_spectral(a[0], a[1], a[2], spec_a);
    rgb_to_spectral(b[0], b[1], b[2], spec_b);
    
    //blend spectral primaries subtractive WGM
    float spectralmix[10] = {0};
    for (int i=0; i < 10; i++) {
      // 'mix_colors' is called infrequently enough that we
      // can afford to not use the faster approximations here.
      spectralmix[i] = powf(spec_a[i], sfac_a) * powf(spec_b[i], sfac_b);
    }
    
    //convert to RGB
    float rgb_result[3] = {0};
    spectral_to_rgb(spectralmix, rgb_result);

    for (int i=0; i < 3; i++) {
      result[i] = rgb_result[i];
    }
  }
  
  if (paint_mode < 1.0) {
    for (int i=0; i<3; i++) {
      result[i] = result[i] * paint_mode + (1-paint_mode) * (a[i] * opa_a + b[i] * opa_b);
    }
  }

  return result;
}

#endif //HELPERS_C


<|start_filename|>mypaint-matrix.c<|end_filename|>
/* libmypaint - The MyPaint Brush Library
 * Copyright (C) 2019 The MyPaint Team
 *
 * 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 "mypaint-matrix.h"
#include <math.h>

MyPaintTransform
mypaint_matrix_multiply(const MyPaintTransform m1, const MyPaintTransform m2)
{
    MyPaintTransform result;
    for (int row = 0; row < 3; ++row) {
      for (int col = 0; col < 3; ++col) {
	  result.rows[row][col] =
	    m1.rows[0][col] * m2.rows[row][0] +
	    m1.rows[1][col] * m2.rows[row][1] +
	    m1.rows[2][col] * m2.rows[row][2];
        }
    }
    return result;
}

MyPaintTransform
mypaint_transform_unit()
{
    MyPaintTransform m = {{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}};
    return m;
}

MyPaintTransform
mypaint_transform_rotate_cw(const MyPaintTransform transform, const float angle_radians)
{
    const float a = angle_radians;
    MyPaintTransform factor = {{{cos(a), sin(a), 0}, {-sin(a), cos(a), 0}, {0, 0, 1}}};
    return mypaint_matrix_multiply(transform, factor);
}

MyPaintTransform
mypaint_transform_rotate_ccw(const MyPaintTransform transform, const float angle_radians)
{
    const float a = angle_radians;
    MyPaintTransform factor = {{
        {cos(a), -sin(a), 0},
        {sin(a), cos(a), 0},
        {0, 0, 1},
    }};
    return mypaint_matrix_multiply(transform, factor);
}

MyPaintTransform
mypaint_transform_reflect(const MyPaintTransform transform, const float angle_radians)
{
    float x = cos(angle_radians);
    float y = sin(angle_radians);
    MyPaintTransform factor = {{
        {x * x - y * y, 2.0 * x * y, 0},
        {2.0 * x * y, y * y - x * x, 0},
        {0, 0, 1},
    }};
    return mypaint_matrix_multiply(transform, factor);
}


MyPaintTransform
mypaint_transform_translate(const MyPaintTransform transform, const float x, const float y)
{
    MyPaintTransform factor = {{
	{1, 0, x},
	{0, 1, y},
	{0, 0, 1},
      }};
    return mypaint_matrix_multiply(transform, factor);
}

void
mypaint_transform_point(const MyPaintTransform* const t, float x, float y, float* xout, float* yout)
{
    *xout = t->rows[0][0] * x + t->rows[0][1] * y + t->rows[0][2];
    *yout = t->rows[1][0] * x + t->rows[1][1] * y + t->rows[1][2];
}


<|start_filename|>fastapprox/fasterf.h<|end_filename|>
/*=====================================================================*
 *                   Copyright (C) 2011 <NAME>                   *
 * All rights reserved.                                                *
 *                                                                     *
 * Redistribution and use in source and binary forms, with             *
 * or without modification, are permitted provided that the            *
 * following conditions are met:                                       *
 *                                                                     *
 *     * Redistributions of source code must retain the                *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer.                                       *
 *                                                                     *
 *     * Redistributions in binary form must reproduce the             *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer in the documentation and/or            *
 *     other materials provided with the distribution.                 *
 *                                                                     *
 *     * Neither the name of <NAME> nor the names                *
 *     of other contributors may be used to endorse or promote         *
 *     products derived from this software without specific            *
 *     prior written permission.                                       *
 *                                                                     *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND              *
 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,         *
 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES               *
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE             *
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER               *
 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,                 *
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES            *
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE           *
 * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR                *
 * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF          *
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT           *
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY              *
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE             *
 * POSSIBILITY OF SUCH DAMAGE.                                         *
 *                                                                     *
 * Contact: <NAME> <<EMAIL>>                            *
 *=====================================================================*/

#ifndef __FAST_ERF_H_
#define __FAST_ERF_H_

#include <math.h>
#include <stdint.h>
#include "sse.h"
#include "fastexp.h"
#include "fastlog.h"

// fasterfc: not actually faster than erfcf(3) on newer machines!
// ... although vectorized version is interesting
//     and fastererfc is very fast

static inline float
fasterfc (float x)
{
  static const float k = 3.3509633149424609f;
  static const float a = 0.07219054755431126f;
  static const float b = 15.418191568719577f;
  static const float c = 5.609846028328545f;

  union { float f; uint32_t i; } vc = { c * x };
  float xsq = x * x;
  float xquad = xsq * xsq;

  vc.i |= 0x80000000;

  return 2.0f / (1.0f + fastpow2 (k * x)) - a * x * (b * xquad - 1.0f) * fasterpow2 (vc.f);
}

static inline float
fastererfc (float x)
{
  static const float k = 3.3509633149424609f;

  return 2.0f / (1.0f + fasterpow2 (k * x));
}

// fasterf: not actually faster than erff(3) on newer machines! 
// ... although vectorized version is interesting
//     and fastererf is very fast

static inline float
fasterf (float x)
{
  return 1.0f - fasterfc (x);
}

static inline float
fastererf (float x)
{
  return 1.0f - fastererfc (x);
}

static inline float
fastinverseerf (float x)
{
  static const float invk = 0.30004578719350504f;
  static const float a = 0.020287853348211326f;
  static const float b = 0.07236892874789555f;
  static const float c = 0.9913030456864257f;
  static const float d = 0.8059775923760193f;

  float xsq = x * x;

  return invk * fastlog2 ((1.0f + x) / (1.0f - x)) 
       + x * (a - b * xsq) / (c - d * xsq);
}

static inline float
fasterinverseerf (float x)
{
  static const float invk = 0.30004578719350504f;

  return invk * fasterlog2 ((1.0f + x) / (1.0f - x));
}

#ifdef __SSE2__

static inline v4sf
vfasterfc (v4sf x)
{
  const v4sf k = v4sfl (3.3509633149424609f);
  const v4sf a = v4sfl (0.07219054755431126f);
  const v4sf b = v4sfl (15.418191568719577f);
  const v4sf c = v4sfl (5.609846028328545f);

  union { v4sf f; v4si i; } vc; vc.f = c * x;
  vc.i |= v4sil (0x80000000);

  v4sf xsq = x * x;
  v4sf xquad = xsq * xsq;

  return v4sfl (2.0f) / (v4sfl (1.0f) + vfastpow2 (k * x)) - a * x * (b * xquad - v4sfl (1.0f)) * vfasterpow2 (vc.f);
}

static inline v4sf
vfastererfc (const v4sf x)
{
  const v4sf k = v4sfl (3.3509633149424609f);

  return v4sfl (2.0f) / (v4sfl (1.0f) + vfasterpow2 (k * x));
}

static inline v4sf
vfasterf (v4sf x)
{
  return v4sfl (1.0f) - vfasterfc (x);
}

static inline v4sf
vfastererf (const v4sf x)
{
  return v4sfl (1.0f) - vfastererfc (x);
}

static inline v4sf
vfastinverseerf (v4sf x)
{
  const v4sf invk = v4sfl (0.30004578719350504f);
  const v4sf a = v4sfl (0.020287853348211326f);
  const v4sf b = v4sfl (0.07236892874789555f);
  const v4sf c = v4sfl (0.9913030456864257f);
  const v4sf d = v4sfl (0.8059775923760193f);

  v4sf xsq = x * x;

  return invk * vfastlog2 ((v4sfl (1.0f) + x) / (v4sfl (1.0f) - x)) 
       + x * (a - b * xsq) / (c - d * xsq);
}

static inline v4sf
vfasterinverseerf (v4sf x)
{
  const v4sf invk = v4sfl (0.30004578719350504f);

  return invk * vfasterlog2 ((v4sfl (1.0f) + x) / (v4sfl (1.0f) - x));
}

#endif //__SSE2__

#endif // __FAST_ERF_H_


<|start_filename|>glib/mypaint-brush.h<|end_filename|>
#ifndef MYPAINTBRUSHGLIB_H
#define MYPAINTBRUSHGLIB_H

#include <glib-object.h>

#define MYPAINT_TYPE_BRUSH               (mypaint_brush_get_type ())
#define MYPAINT_VALUE_HOLDS_BRUSH(value) (G_TYPE_CHECK_VALUE_TYPE ((value), MYPAINT_TYPE_BRUSH))
GType mypaint_brush_get_type(void);

#define MYPAINT_TYPE_SURFACE               (mypaint_surface_get_type ())
#define MYPAINT_VALUE_HOLDS_SURFACE(value) (G_TYPE_CHECK_VALUE_TYPE ((value), MYPAINT_TYPE_SURFACE))
GType mypaint_surface_get_type(void);

#define MYPAINT_TYPE_FIXED_TILED_SURFACE               (mypaint_fixed_tiled_surface_get_type ())
#define MYPAINT_VALUE_HOLDS_FIXED_TILED_SURFACE(value) (G_TYPE_CHECK_VALUE_TYPE ((value), MYPAINT_TYPE_FIXED_TILED_SURFACE))
GType mypaint_fixed_tiled_surface_get_type(void);

#define MYPAINT_TYPE_RECTANGLE               (mypaint_rectangle_get_type ())
#define MYPAINT_VALUE_HOLDS_RECTANGLE(value) (G_TYPE_CHECK_VALUE_TYPE ((value), MYPAINT_TYPE_RECTANGLE))
GType mypaint_rectangle_get_type(void);

#endif // MYPAINTBRUSHGLIB_H


<|start_filename|>mypaint-symmetry.c<|end_filename|>
#include "mypaint-symmetry.h"
#include "helpers.h"

#include <stdio.h>
#include <stdlib.h>

#define DEFAULT_NUM_MATRICES 16

void
allocation_failure_warning(int num)
{
    fprintf(stderr, "Critical: failed to allocate memory for %d transformation matrices!\n", num);
}

gboolean
allocate_symmetry_matrices(MyPaintSymmetryData* data, int num_matrices)
{
    int bytes = num_matrices * sizeof(MyPaintTransform);
    void* allocated = realloc(data->symmetry_matrices, bytes);
    if (!allocated) {
        allocation_failure_warning(num_matrices);
        data->num_symmetry_matrices = 0;
        return FALSE;
    } else {
        data->symmetry_matrices = allocated;
        data->num_symmetry_matrices = num_matrices;
        return TRUE;
    }
}

gboolean
symmetry_states_equal(const MyPaintSymmetryState* const s1, const MyPaintSymmetryState* const s2)
{
    return s1->type == s2->type && s1->center_x == s2->center_x && s1->center_y == s2->center_y &&
           s1->angle == s2->angle && s1->num_lines == s2->num_lines;
}

int
num_matrices_required(const MyPaintSymmetryState* state)
{
    switch (state->type) {
    case MYPAINT_SYMMETRY_TYPE_VERTICAL:
    case MYPAINT_SYMMETRY_TYPE_HORIZONTAL:
        return 1;
    case MYPAINT_SYMMETRY_TYPE_VERTHORZ:
        return 3;
    case MYPAINT_SYMMETRY_TYPE_ROTATIONAL:
        return state->num_lines - 1;
    case MYPAINT_SYMMETRY_TYPE_SNOWFLAKE:
        return 2 * state->num_lines - 1;
    default:
        return 0;
    }
}

/* Public functions */

MyPaintSymmetryState
symmetry_state_default()
{
    MyPaintSymmetryState base = {
        .type = MYPAINT_SYMMETRY_TYPE_VERTICAL, .center_x = 0.0, .center_y = 0.0, .angle = 0.0, .num_lines = 2};
    return base;
}

/* If the symmetry state has changed since last, recalculate matrices */
void
mypaint_update_symmetry_state(MyPaintSymmetryData* const self)
{
    if (!self->pending_changes || symmetry_states_equal(&self->state_current, &self->state_pending)) return;
    // Need to recalculate matrices
    // Check if we need to allocate more space
    const int required = num_matrices_required(&self->state_pending);
    if (self->num_symmetry_matrices < required) {
        // Try to allocate space for matrices, skip recalculations if it fails
        if (!allocate_symmetry_matrices(self, required)) return;
    }
    const MyPaintSymmetryState symm = self->state_pending;
    self->state_current = symm;
    float cx = symm.center_x;
    float cy = symm.center_y;
    // Convert angle to radians
    float angle = symm.angle * (M_PI / 180.0);
    float rot_angle = (2.0 * M_PI) / symm.num_lines;
    MyPaintTransform* matrices = self->symmetry_matrices;
    const MyPaintTransform m = mypaint_transform_translate(mypaint_transform_unit(), -cx, -cy);
    switch (symm.type) {
    case MYPAINT_SYMMETRY_TYPE_HORIZONTAL:
    case MYPAINT_SYMMETRY_TYPE_VERTICAL: {
        if (symm.type == MYPAINT_SYMMETRY_TYPE_VERTICAL) {
            angle += M_PI / 2.0;
        }
        matrices[0] = mypaint_transform_reflect(m, -angle);
    } break;
    case MYPAINT_SYMMETRY_TYPE_VERTHORZ: {
        float v_angle = angle + M_PI / 2.0;
        matrices[0] = mypaint_transform_reflect(m, -angle);
        matrices[1] = mypaint_transform_reflect(matrices[0], -v_angle);
        matrices[2] = mypaint_transform_reflect(matrices[1], -angle);
    } break;
    case MYPAINT_SYMMETRY_TYPE_SNOWFLAKE: {
        int base_idx = symm.num_lines - 1;
        for (int i = 0; i < symm.num_lines; ++i) {
            matrices[base_idx + i] =
                mypaint_transform_reflect(mypaint_transform_rotate_cw(m, rot_angle * i), -i * rot_angle - angle);
        }
    }
    case MYPAINT_SYMMETRY_TYPE_ROTATIONAL: {
        for (int i = 1; i < symm.num_lines; ++i) {
            matrices[i - 1] = mypaint_transform_rotate_cw(m, rot_angle * i);
        }
    } break;
    default:
        fprintf(stderr, "Warning: Unhandled symmetry type: %d\n", symm.type);
        return;
    }
    for (int i = 0; i < required; ++i) {
        matrices[i] = mypaint_transform_translate(matrices[i], cx, cy);
    }
    self->pending_changes = FALSE;
}

MyPaintSymmetryData
mypaint_default_symmetry_data()
{
    MyPaintSymmetryData symm_data = {
        .state_current = {.type = -1},
        .state_pending = symmetry_state_default(),
        .pending_changes = TRUE,
        .active = FALSE,
        .num_symmetry_matrices = DEFAULT_NUM_MATRICES,
        .symmetry_matrices = NULL,
    };
    if (allocate_symmetry_matrices(&symm_data, DEFAULT_NUM_MATRICES)) {
        mypaint_update_symmetry_state(&symm_data);
    }
    return symm_data;
}

void
mypaint_symmetry_data_destroy(MyPaintSymmetryData* data)
{
    if (data->symmetry_matrices != NULL) {
        free(data->symmetry_matrices);
    }
}

void
mypaint_symmetry_set_pending(
    MyPaintSymmetryData* data, gboolean active, float center_x, float center_y, float symmetry_angle,
    MyPaintSymmetryType symmetry_type, int rot_symmetry_lines)
{
    data->active = active;
    data->state_pending.center_x = center_x;
    data->state_pending.center_y = center_y;
    data->state_pending.type = symmetry_type;
    data->state_pending.num_lines = MAX(2, rot_symmetry_lines);
    data->state_pending.angle = symmetry_angle;

    data->pending_changes = TRUE;
}


<|start_filename|>tests/test-fixed-tiled-surface.c<|end_filename|>

#include <stddef.h>

#include "mypaint-fixed-tiled-surface.h"
#include "mypaint-test-surface.h"

MyPaintSurface *
fixed_surface_factory(gpointer user_data)
{
    MyPaintFixedTiledSurface * surface = mypaint_fixed_tiled_surface_new(1000, 1000);
    return (MyPaintSurface *)surface;
}

int
main(int argc, char **argv)
{
    return mypaint_test_surface_run(argc, argv, fixed_surface_factory, "MyPaintFixedSurface", NULL);
}



<|start_filename|>tiled-surface-private.h<|end_filename|>


void render_dab_mask (uint16_t * mask,
                        float x, float y,
                        float radius,
                        float hardness,
                        float softness,
                        float aspect_ratio, float angle
                        );


<|start_filename|>mypaint-tiled-surface.c<|end_filename|>
/* libmypaint - The MyPaint Brush Library
 * Copyright (C) 2007-2014 <NAME> <<EMAIL>> et. al.
 *
 * 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 <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>

#ifdef _OPENMP
#include <omp.h>
#endif

#include "mypaint-config.h"
#include "mypaint-tiled-surface.h"
#include "tiled-surface-private.h"
#include "helpers.h"
#include "brushmodes.h"
#include "operationqueue.h"

void process_tile(MyPaintTiledSurface *self, int tx, int ty);

static void
begin_atomic_default(MyPaintSurface *surface)
{
    mypaint_tiled_surface_begin_atomic((MyPaintTiledSurface *)surface);
}

static void
end_atomic_default(MyPaintSurface *surface, MyPaintRectangles *roi)
{
    mypaint_tiled_surface_end_atomic((MyPaintTiledSurface *)surface, roi);
}

void
prepare_bounding_boxes(MyPaintTiledSurface *self) {
    MyPaintSymmetryState symm_state = self->symmetry_data.state_current;
    const gboolean snowflake = symm_state.type == MYPAINT_SYMMETRY_TYPE_SNOWFLAKE;
    const int num_bboxes_desired = symm_state.num_lines * (snowflake ? 2 : 1);
    // If the bounding box array cannot fit one rectangle per symmetry dab,
    // try to allocate enough space for that to be possible.
    // Failure is ok, as the bounding box assignments will be functional anyway.
    if (num_bboxes_desired > self->num_bboxes) {
        const int margin = 10; // Add margin to avoid unnecessary reallocations.
        const int num_to_allocate = num_bboxes_desired + margin;
        int bytes_to_allocate = num_to_allocate * sizeof(MyPaintRectangle);
        MyPaintRectangle* new_bboxes = malloc(bytes_to_allocate);
        if (new_bboxes) {
            if (self->num_bboxes > NUM_BBOXES_DEFAULT) {
                // Free previous allocation
                free(self->bboxes);
            }
            // Initialize memory
            memset(new_bboxes, 0, bytes_to_allocate);
            self->bboxes = new_bboxes;
            self->num_bboxes = num_to_allocate;
            // No need to clear anything after the memset, so reset counter
            self->num_bboxes_dirtied = 0;
        }
    }
    // Clean up any previously populated bounding boxes and reset the counter
    for (int i = 0; i < MIN(self->num_bboxes, self->num_bboxes_dirtied); ++i) {
        self->bboxes[i].height = 0;
        self->bboxes[i].width = 0;
        self->bboxes[i].x = 0;
        self->bboxes[i].y = 0;
    }
    self->num_bboxes_dirtied = 0;
}

/**
 * mypaint_tiled_surface_begin_atomic: (skip)
 *
 * Implementation of #MyPaintSurface::being_atomic vfunc
 * Note: Only intended to be used from #MyPaintTiledSurface subclasses, which should chain up to this
 * if implementing their own #MyPaintSurface::begin_atomic vfunc.
 * Application code should only use mypaint_surface_being_atomic()
 */
void
mypaint_tiled_surface_begin_atomic(MyPaintTiledSurface *self)
{
    mypaint_update_symmetry_state(&self->symmetry_data);
    prepare_bounding_boxes(self);
}

/**
 * mypaint_tiled_surface_end_atomic: (skip)
 *
 * Implementation of #MyPaintSurface::end_atomic vfunc
 * Note: Only intended to be used from #MyPaintTiledSurface subclasses, which should chain up to this
 * if implementing their own #MyPaintSurface::end_atomic vfunc.
 * Application code should only use mypaint_surface_end_atomic().
 */
void
mypaint_tiled_surface_end_atomic(MyPaintTiledSurface *self, MyPaintRectangles *roi)
{
    // Process tiles
    TileIndex *tiles;
    int tiles_n = operation_queue_get_dirty_tiles(self->operation_queue, &tiles);

    #pragma omp parallel for schedule(static) if(self->threadsafe_tile_requests && tiles_n > 3)
    for (int i = 0; i < tiles_n; i++) {
        process_tile(self, tiles[i].x, tiles[i].y);
    }

    operation_queue_clear_dirty_tiles(self->operation_queue);

    if (roi) {
        const int roi_rects = roi->num_rectangles;
        const int num_dirty = self->num_bboxes_dirtied;
        // Clear out the input rectangles that will be overwritten
        for (int i = 0; i < MIN(roi_rects, num_dirty); ++i) {
            roi->rectangles[i].x = 0;
            roi->rectangles[i].y = 0;
            roi->rectangles[i].width = 0;
            roi->rectangles[i].height = 0;
        }
        // Write bounding box rectangles to the output array
        const float bboxes_per_output = MAX(1, (float)num_dirty / roi_rects);
        for (int i = 0; i < num_dirty; ++i) {
            int out_index;
            // If there is not enough space for all rectangles in the output,
            // merge some of the rectangles with their list-adjacent neighbours.
            if (num_dirty > roi_rects) {
                out_index = (int)MIN(roi_rects - 1, roundf((float)i / bboxes_per_output));
            } else {
                out_index = i;
            }
            mypaint_rectangle_expand_to_include_rect(&(roi->rectangles[out_index]), &(self->bboxes[i]));
        }
        // Set the number of rectangles written to, so the caller knows which ones to act on.
        roi->num_rectangles = MIN(roi_rects, num_dirty);
    }
}

/**
 * mypaint_tiled_surface_tile_request_start:
 *
 * Fetch a tile out from the underlying tile store.
 * When successful, request->data will be set to point to the fetched tile.
 * Consumers must *always* call mypaint_tiled_surface_tile_request_end() with the same
 * request to complete the transaction.
 */
void mypaint_tiled_surface_tile_request_start(MyPaintTiledSurface *self, MyPaintTileRequest *request)
{
    assert(self->tile_request_start);
    self->tile_request_start(self, request);
}

/**
 * mypaint_tiled_surface_tile_request_end:
 *
 * Put a (potentially modified) tile back into the underlying tile store.
 *
 * Consumers must *always* call mypaint_tiled_surface_tile_request_start() with the same
 * request to start the transaction before calling this function.
 */
void mypaint_tiled_surface_tile_request_end(MyPaintTiledSurface *self, MyPaintTileRequest *request)
{
    assert(self->tile_request_end);
    self->tile_request_end(self, request);
}

/* FIXME: either expose this through MyPaintSurface, or move it into the brush engine */
/**
 * mypaint_tiled_surface_set_symmetry_state:
 * @active: TRUE to enable, FALSE to disable.
 * @center_x: X axis to mirror events across.
 * @center_y: Y axis to mirror events across.
 * @symmetry_angle: Angle to rotate the symmetry lines
 * @symmetry_type: Symmetry type to activate.
 * @rot_symmetry_lines: Number of rotational symmetry lines.
 *
 * Enable/Disable symmetric brush painting across an X axis.
 *
 */
void
mypaint_tiled_surface_set_symmetry_state(MyPaintTiledSurface *self, gboolean active,
                                         float center_x, float center_y,
                                         float symmetry_angle,
                                         MyPaintSymmetryType symmetry_type,
                                         int rot_symmetry_lines)
{
    mypaint_symmetry_set_pending( // Only write to the pending new state, nothing gets recalculated here
        &self->symmetry_data, active, center_x, center_y, symmetry_angle, symmetry_type, rot_symmetry_lines);
}

/**
 * mypaint_tile_request_init:
 *
 * Initialize a request for use with mypaint_tiled_surface_tile_request_start()
 * and mypaint_tiled_surface_tile_request_end()
 */
void
mypaint_tile_request_init(MyPaintTileRequest *data, int level,
                          int tx, int ty, gboolean readonly)
{
    data->tx = tx;
    data->ty = ty;
    data->readonly = readonly;
    data->buffer = NULL;
    data->context = NULL;
#ifdef _OPENMP
    data->thread_id = omp_get_thread_num();
#else
    data->thread_id = -1;
#endif
    data->mipmap_level = level;
}

// Must be threadsafe
static inline float
calculate_r_sample(float x, float y, float aspect_ratio,
                      float sn, float cs)
{
    const float yyr=(y*cs-x*sn)*aspect_ratio;
    const float xxr=y*sn+x*cs;
    const float r = (yyr*yyr + xxr*xxr);
    return r;
}

static inline float
calculate_rr(int xp, int yp, float x, float y, float aspect_ratio,
                      float sn, float cs, float one_over_radius2)
{
    // code duplication, see brush::count_dabs_to()
    const float yy = (yp + 0.5f - y);
    const float xx = (xp + 0.5f - x);
    const float yyr=(yy*cs-xx*sn)*aspect_ratio;
    const float xxr=yy*sn+xx*cs;
    const float rr = (yyr*yyr + xxr*xxr) * one_over_radius2;
    // rr is in range 0.0..1.0*sqrt(2)
    return rr;
}

static inline float
sign_point_in_line( float px, float py, float vx, float vy )
{
    return (px - vx) * (-vy) - (vx) * (py - vy);
}

static inline void
closest_point_to_line( float lx, float ly, float px, float py, float *ox, float *oy )
{
    const float l2 = lx*lx + ly*ly;
    const float ltp_dot = px*lx + py*ly;
    const float t = ltp_dot / l2;
    *ox = lx * t;
    *oy = ly * t;
}

// Must be threadsafe
//
// This works by taking the visibility at the nearest point
// and dividing by 1.0 + delta.
//
// - nearest point: point where the dab has more influence
// - farthest point: point at a fixed distance away from
//                   the nearest point
// - delta: how much occluded is the farthest point relative
//          to the nearest point
static inline float
calculate_rr_antialiased(int xp, int yp, float x, float y, float aspect_ratio,
                      float sn, float cs, float one_over_radius2,
                      float r_aa_start)
{
    // calculate pixel position and borders in a way
    // that the dab's center is always at zero
    float pixel_right = x - (float)xp;
    float pixel_bottom = y - (float)yp;
    float pixel_center_x = pixel_right - 0.5f;
    float pixel_center_y = pixel_bottom - 0.5f;
    float pixel_left = pixel_right - 1.0f;
    float pixel_top = pixel_bottom - 1.0f;

    float nearest_x, nearest_y; // nearest to origin, but still inside pixel
    float farthest_x, farthest_y; // farthest from origin, but still inside pixel
    float r_near, r_far, rr_near, rr_far;
    // Dab's center is inside pixel?
    if( pixel_left<0 && pixel_right>0 &&
        pixel_top<0 && pixel_bottom>0 )
    {
        nearest_x = 0;
        nearest_y = 0;
        r_near = rr_near = 0;
    }
    else
    {
        closest_point_to_line( cs, sn, pixel_center_x, pixel_center_y, &nearest_x, &nearest_y );
        nearest_x = CLAMP( nearest_x, pixel_left, pixel_right );
        nearest_y = CLAMP( nearest_y, pixel_top, pixel_bottom );
        // XXX: precision of "nearest" values could be improved
        // by intersecting the line that goes from nearest_x/Y to 0
        // with the pixel's borders here, however the improvements
        // would probably not justify the perdormance cost.
        r_near = calculate_r_sample( nearest_x, nearest_y, aspect_ratio, sn, cs );
        rr_near = r_near * one_over_radius2;
    }

    // out of dab's reach?
    if( rr_near > 1.0f )
        return rr_near;

    // check on which side of the dab's line is the pixel center
    float center_sign = sign_point_in_line( pixel_center_x, pixel_center_y, cs, -sn );

    // radius of a circle with area=1
    //   A = pi * r * r
    //   r = sqrt(1/pi)
    const float rad_area_1 = sqrtf( 1.0f / M_PI );

    // center is below dab
    if( center_sign < 0 )
    {
        farthest_x = nearest_x - sn*rad_area_1;
        farthest_y = nearest_y + cs*rad_area_1;
    }
    // above dab
    else
    {
        farthest_x = nearest_x + sn*rad_area_1;
        farthest_y = nearest_y - cs*rad_area_1;
    }

    r_far = calculate_r_sample( farthest_x, farthest_y, aspect_ratio, sn, cs );
    rr_far = r_far * one_over_radius2;

    // check if we can skip heavier AA
    if( r_far < r_aa_start )
        return (rr_far+rr_near) * 0.5f;

    // calculate AA approximate
    float visibilityNear = 1.0f - rr_near;
    float delta = rr_far - rr_near;
    float delta2 = 1.0f + delta;
    visibilityNear /= delta2;

    return 1.0f - visibilityNear;
}

static inline float
calculate_opa(float rr, float hardness,
              float segment1_offset, float segment1_slope,
              float segment2_offset, float segment2_slope) {

    const float fac = rr <= hardness ? segment1_slope : segment2_slope;
    float opa = rr <= hardness ? segment1_offset : segment2_offset;
    opa += rr*fac;

    if (rr > 1.0f) {
        opa = 0.0f;
    }
    #ifdef HEAVY_DEBUG
    assert(isfinite(opa));
    assert(opa >= 0.0f && opa <= 1.0f);
    #endif
    return opa;
}

// Must be threadsafe
void render_dab_mask (uint16_t * mask,
                        float x, float y,
                        float radius,
                        float hardness,
                        float softness,
                        float aspect_ratio, float angle
                        )
{

    hardness = CLAMP(hardness, 0.0, 1.0);
    if (aspect_ratio<1.0) aspect_ratio=1.0;
    assert(hardness != 0.0); // assured by caller

    // For a graphical explanation, see:
    // http://wiki.mypaint.info/Development/Documentation/Brushlib
    //
    // The hardness calculation is explained below:
    //
    // Dab opacity gradually fades out from the center (rr=0) to
    // fringe (rr=1) of the dab. How exactly depends on the hardness.
    // We use two linear segments, for which we pre-calculate slope
    // and offset here.
    //
    // opa
    // ^
    // *   .
    // |        *
    // |          .
    // +-----------*> rr = (distance_from_center/radius)^2
    // 0           1
    //

    float segment1_offset = (1.f)*(1.f-softness);
    float segment1_slope  = -(1.0f/hardness - 1.0f)*(1.f-softness);
    float segment2_offset = hardness/(1.0f-hardness)*(1.f-softness);
    float segment2_slope  = -hardness/(1.0f-hardness)*(1.f-softness);
    // for hardness == 1.0, segment2 will never be used

    float angle_rad=angle/360*2*M_PI;
    float cs=cos(angle_rad);
    float sn=sin(angle_rad);

    const float r_fringe = radius + 1.0f; // +1.0 should not be required, only to be sure
    int x0 = floor (x - r_fringe);
    int y0 = floor (y - r_fringe);
    int x1 = floor (x + r_fringe);
    int y1 = floor (y + r_fringe);
    if (x0 < 0) x0 = 0;
    if (y0 < 0) y0 = 0;
    if (x1 > MYPAINT_TILE_SIZE-1) x1 = MYPAINT_TILE_SIZE-1;
    if (y1 > MYPAINT_TILE_SIZE-1) y1 = MYPAINT_TILE_SIZE-1;
    const float one_over_radius2 = 1.0f/(radius*radius);

    // Pre-calculate rr and put it in the mask.
    // This an optimization that makes use of auto-vectorization
    // OPTIMIZE: if using floats for the brush engine, store these directly in the mask
    float rr_mask[MYPAINT_TILE_SIZE*MYPAINT_TILE_SIZE+2*MYPAINT_TILE_SIZE];

    if (radius < 3.0f)
    {
      const float aa_border = 1.0f;
      float r_aa_start = ((radius>aa_border) ? (radius-aa_border) : 0);
      r_aa_start *= r_aa_start / aspect_ratio;

      for (int yp = y0; yp <= y1; yp++) {
        for (int xp = x0; xp <= x1; xp++) {
          const float rr = calculate_rr_antialiased(xp, yp,
                                  x, y, aspect_ratio,
                                  sn, cs, one_over_radius2,
                                  r_aa_start);
          rr_mask[(yp*MYPAINT_TILE_SIZE)+xp] = rr;
        }
      }
    }
    else
    {
      for (int yp = y0; yp <= y1; yp++) {
        for (int xp = x0; xp <= x1; xp++) {
          const float rr = calculate_rr(xp, yp,
                                  x, y, aspect_ratio,
                                  sn, cs, one_over_radius2);
          rr_mask[(yp*MYPAINT_TILE_SIZE)+xp] = rr;
        }
      }
    }

    // we do run length encoding: if opacity is zero, the next
    // value in the mask is the number of pixels that can be skipped.
    uint16_t * mask_p = mask;
    int skip=0;

    skip += y0*MYPAINT_TILE_SIZE;
    for (int yp = y0; yp <= y1; yp++) {
      skip += x0;

      int xp;
      for (xp = x0; xp <= x1; xp++) {
        const float rr = rr_mask[(yp*MYPAINT_TILE_SIZE)+xp];
        const float opa = calculate_opa(rr, hardness,
                                  segment1_offset, segment1_slope,
                                  segment2_offset, segment2_slope);
        const uint16_t opa_ = opa * (1<<15);
        if (!opa_) {
          skip++;
        } else {
          if (skip) {
            *mask_p++ = 0;
            *mask_p++ = skip*4;
            skip = 0;
          }
          *mask_p++ = opa_;
        }
      }
      skip += MYPAINT_TILE_SIZE-xp;
    }
    *mask_p++ = 0;
    *mask_p++ = 0;
  }

// Must be threadsafe
void
process_op(uint16_t *rgba_p, uint16_t *mask,
           int tx, int ty, OperationDataDrawDab *op)
{

    // first, we calculate the mask (opacity for each pixel)
    render_dab_mask(mask,
                    op->x - tx*MYPAINT_TILE_SIZE,
                    op->y - ty*MYPAINT_TILE_SIZE,
                    op->radius,
                    op->hardness,
                    op->softness,
                    op->aspect_ratio, op->angle
                    );

    // second, we use the mask to stamp a dab for each activated blend mode
    if (op->paint < 1.0) {
      if (op->normal) {
        if (op->color_a == 1.0) {
          draw_dab_pixels_BlendMode_Normal(mask, rgba_p,
                                           op->color_r, op->color_g, op->color_b, op->normal*op->opaque*(1 - op->paint)*(1<<15));
        } else {
          // normal case for brushes that use smudging (eg. watercolor)
          draw_dab_pixels_BlendMode_Normal_and_Eraser(mask, rgba_p,
                                                      op->color_r, op->color_g, op->color_b, op->color_a*(1<<15),
                                                      op->normal*op->opaque*(1 - op->paint)*(1<<15));
        }
      }

      if (op->lock_alpha && op->color_a != 0) {
        draw_dab_pixels_BlendMode_LockAlpha(mask, rgba_p,
                                            op->color_r, op->color_g, op->color_b,
                                            op->lock_alpha*op->opaque*(1 - op->colorize)*(1 - op->posterize)*(1 - op->paint)*(1<<15));
      }
    }
    
    if (op->paint > 0.0) {
      if (op->normal) {
        if (op->color_a == 1.0) {
          draw_dab_pixels_BlendMode_Normal_Paint(mask, rgba_p,
                                           op->color_r, op->color_g, op->color_b, op->normal*op->opaque*op->paint*(1<<15));
        } else {
          // normal case for brushes that use smudging (eg. watercolor)
          draw_dab_pixels_BlendMode_Normal_and_Eraser_Paint(mask, rgba_p,
                                                      op->color_r, op->color_g, op->color_b, op->color_a*(1<<15),
                                                      op->normal*op->opaque*op->paint*(1<<15));
        }
      }

      if (op->lock_alpha && op->color_a != 0) {
        draw_dab_pixels_BlendMode_LockAlpha_Paint(mask, rgba_p,
                                            op->color_r, op->color_g, op->color_b,
                                            op->lock_alpha*op->opaque*(1 - op->colorize)*(1 - op->posterize)*op->paint*(1<<15));
      }
    }
    
    if (op->colorize) {
      draw_dab_pixels_BlendMode_Color(mask, rgba_p,
                                      op->color_r, op->color_g, op->color_b,
                                      op->colorize*op->opaque*(1<<15));
    }
    if (op->posterize) {
      draw_dab_pixels_BlendMode_Posterize(mask, rgba_p,
                                      op->posterize*op->opaque*(1<<15),
                                      op->posterize_num);
    }
}

// Must be threadsafe
void
process_tile(MyPaintTiledSurface *self, int tx, int ty)
{
    TileIndex tile_index = {tx, ty};
    OperationDataDrawDab *op = operation_queue_pop(self->operation_queue, tile_index);
    if (!op) {
        return;
    }

    MyPaintTileRequest request_data;
    const int mipmap_level = 0;
    mypaint_tile_request_init(&request_data, mipmap_level, tx, ty, FALSE);

    mypaint_tiled_surface_tile_request_start(self, &request_data);
    uint16_t * rgba_p = request_data.buffer;
    if (!rgba_p) {
        printf("Warning: Unable to get tile!\n");
        return;
    }

    uint16_t mask[MYPAINT_TILE_SIZE*MYPAINT_TILE_SIZE+2*MYPAINT_TILE_SIZE];

    while (op) {
        process_op(rgba_p, mask, tile_index.x, tile_index.y, op);
        free(op);
        op = operation_queue_pop(self->operation_queue, tile_index);
    }

    mypaint_tiled_surface_tile_request_end(self, &request_data);
}

void
update_dirty_bbox(MyPaintRectangle *bbox, OperationDataDrawDab *op)
{
    int bb_x, bb_y, bb_w, bb_h;
    float r_fringe = op->radius + 1.0f; // +1.0 should not be required, only to be sure
    bb_x = floor (op->x - r_fringe);
    bb_y = floor (op->y - r_fringe);
    bb_w = floor (op->x + r_fringe) - bb_x + 1;
    bb_h = floor (op->y + r_fringe) - bb_y + 1;

    mypaint_rectangle_expand_to_include_point(bbox, bb_x, bb_y);
    mypaint_rectangle_expand_to_include_point(bbox, bb_x+bb_w-1, bb_y+bb_h-1);
}

// returns TRUE if the surface was modified
gboolean draw_dab_internal (MyPaintTiledSurface *self, float x, float y,
               float radius,
               float color_r, float color_g, float color_b,
               float opaque, float hardness, float softness,
               float color_a,
               float aspect_ratio, float angle,
               float lock_alpha,
               float colorize,
               float posterize,
               float posterize_num,
               float paint,
               int bbox_index
               )

{
    OperationDataDrawDab op_struct;
    OperationDataDrawDab *op = &op_struct;

    op->x = x;
    op->y = y;
    op->radius = radius;
    op->aspect_ratio = aspect_ratio;
    op->angle = angle;
    op->opaque = CLAMP(opaque, 0.0f, 1.0f);
    op->hardness = CLAMP(hardness, 0.0f, 1.0f);
    op->softness = CLAMP(softness, 0.0f, 1.0f);
    op->lock_alpha = CLAMP(lock_alpha, 0.0f, 1.0f);
    op->colorize = CLAMP(colorize, 0.0f, 1.0f);
    op->posterize = CLAMP(posterize, 0.0f, 1.0f);
    op->posterize_num= CLAMP(ROUND(posterize_num * 100.0), 1, 128);
    op->paint = CLAMP(paint, 0.0f, 1.0f);
    if (op->radius < 0.1f) return FALSE; // don't bother with dabs smaller than 0.1 pixel
    if (op->hardness == 0.0f) return FALSE; // infintly small center point, fully transparent outside
    if (op->softness == 1.0f) return FALSE;
    if (op->opaque == 0.0f) return FALSE;

    color_r = CLAMP(color_r, 0.0f, 1.0f);
    color_g = CLAMP(color_g, 0.0f, 1.0f);
    color_b = CLAMP(color_b, 0.0f, 1.0f);
    color_a = CLAMP(color_a, 0.0f, 1.0f);

    op->color_r = color_r * (1<<15);
    op->color_g = color_g * (1<<15);
    op->color_b = color_b * (1<<15);
    op->color_a = color_a;

    // blending mode preparation
    op->normal = 1.0f;

    op->normal *= 1.0f-op->lock_alpha;
    op->normal *= 1.0f-op->colorize;
    op->normal *= 1.0f-op->posterize;

    if (op->aspect_ratio<1.0f) op->aspect_ratio=1.0f;

    // Determine the tiles influenced by operation, and queue it for processing for each tile
    float r_fringe = radius + 1.0f; // +1.0 should not be required, only to be sure
      
    int tx1 = floor(floor(x - r_fringe) / MYPAINT_TILE_SIZE);
    int tx2 = floor(floor(x + r_fringe) / MYPAINT_TILE_SIZE);
    int ty1 = floor(floor(y - r_fringe) / MYPAINT_TILE_SIZE);
    int ty2 = floor(floor(y + r_fringe) / MYPAINT_TILE_SIZE);

    for (int ty = ty1; ty <= ty2; ty++) {
        for (int tx = tx1; tx <= tx2; tx++) {
            const TileIndex tile_index = {tx, ty};
            OperationDataDrawDab *op_copy = (OperationDataDrawDab *)malloc(sizeof(OperationDataDrawDab));
            *op_copy = *op;
            operation_queue_add(self->operation_queue, tile_index, op_copy);
        }
    }

    update_dirty_bbox(&self->bboxes[bbox_index], op);

    return TRUE;
}


// returns TRUE if the surface was modified
int draw_dab (MyPaintSurface *surface, float x, float y,
               float radius,
               float color_r, float color_g, float color_b,
               float opaque, float hardness, float softness,
               float color_a,
               float aspect_ratio, float angle,
               float lock_alpha,
               float colorize,
               float posterize,
               float posterize_num,
               float paint)
{
    MyPaintTiledSurface* self = (MyPaintTiledSurface*)surface;

    // These calls are repeated enough to warrant a local macro, for both readability and correctness.
#define DDI(x, y, angle, bb_idx) (draw_dab_internal(\
        self, (x), (y), radius, color_r, color_g, color_b, opaque, \
        hardness, softness, color_a, aspect_ratio, (angle), \
        lock_alpha, colorize, posterize, posterize_num, paint, (bb_idx)))

    // Normal pass
    gboolean surface_modified = DDI(x, y, angle, 0);

    int num_bboxes_used = surface_modified ? 1 : 0;

    // Symmetry pass

    // OPTIMIZATION: skip the symmetry pass if surface was not modified by the initial dab;
    // at current if the initial dab does not modify the surface, none of the symmetry dabs
    // will either. If/when selection masks are added, this optimization _must_ be removed,
    // and `surface_modified` must be or'ed with the result of each call to draw_dab_internal.
    MyPaintSymmetryData *symm_data = &self->symmetry_data;
    if (surface_modified && symm_data->active && symm_data->num_symmetry_matrices) {
        const MyPaintSymmetryState symm = symm_data->state_current;
        const int num_bboxes = self->num_bboxes;
        const float rot_angle = 360.0 / symm.num_lines;
        const MyPaintTransform* const matrices = symm_data->symmetry_matrices;
        float x_out, y_out;

        switch (symm.type) {
        case MYPAINT_SYMMETRY_TYPE_VERTICAL: {
            mypaint_transform_point(&matrices[0], x, y, &x_out, &y_out);
            DDI(x_out, y_out, -2.0 * (90 + symm.angle) - angle, 1);
            num_bboxes_used = 2;
            break;
        }
        case MYPAINT_SYMMETRY_TYPE_HORIZONTAL: {
            mypaint_transform_point(&matrices[0], x, y, &x_out, &y_out);
            DDI(x_out, y_out, -2.0 * symm.angle - angle, 1);
            num_bboxes_used = 2;
            break;
        }
        case MYPAINT_SYMMETRY_TYPE_VERTHORZ: {
            // Reflect across horizontal line
            mypaint_transform_point(&matrices[0], x, y, &x_out, &y_out);
            DDI(x_out, y_out, -2.0 * symm.angle - angle, 1);
            // Then across the vertical line (diagonal)
            mypaint_transform_point(&matrices[1], x, y, &x_out, &y_out);
            DDI(x_out, y_out, angle, 2);
            // Then back across the horizontal line
            mypaint_transform_point(&matrices[2], x, y, &x_out, &y_out);
            DDI(x_out, y_out, -2.0 * symm.angle - angle, 3);
            num_bboxes_used = 4;
            break;
        }
        case MYPAINT_SYMMETRY_TYPE_SNOWFLAKE: {

            // These dabs will occupy the bboxes after the last bbox used by the rotational dabs.
            const int offset = MIN(num_bboxes / 2, symm.num_lines);
            const float dabs_per_bbox = MAX(1, (float)symm.num_lines * 2.0 / num_bboxes);
            const int base_idx = symm.num_lines - 1;
            const float base_angle = -2 * symm.angle - angle;
            // draw snowflake dabs for _all_ symmetry lines as we need to reflect the initial dab.
            for (int dab_count = 0; dab_count < symm.num_lines; dab_count++) {
                // If the number of bboxes cannot fit all snowflake dabs, use half for the rotational dabs
                // and the other half for the reflected dabs. This is not always optimal, but seldom bad.
                const int bbox_idx = offset + MIN(roundf(dab_count / dabs_per_bbox), num_bboxes - 1);
                mypaint_transform_point(&matrices[base_idx + dab_count], x, y, &x_out, &y_out);
                DDI(x_out, y_out, base_angle - dab_count * rot_angle, bbox_idx);
            }
            num_bboxes_used = MIN(self->num_bboxes, symm.num_lines * 2);
            // fall through to rotational to finish the process
        }
        case MYPAINT_SYMMETRY_TYPE_ROTATIONAL: {

            // Set the dab bbox distribution factor based on whether the pass is only
            // rotational, or following a snowflake pass. For the latter, we compress
            // the available range (unimportant if there are enough bboxes to go around).
            const gboolean snowflake = symm.type == MYPAINT_SYMMETRY_TYPE_SNOWFLAKE;
            float dabs_per_bbox = MAX(1, (float)(symm.num_lines * (snowflake ? 2 : 1)) / num_bboxes);

            // draw self->rot_symmetry_lines - 1 rotational dabs since initial pass handles the first dab
            for (int dab_count = 1; dab_count < symm.num_lines; dab_count++) {
                const int bbox_index = MIN(roundf(dab_count / dabs_per_bbox), num_bboxes - 1);
                mypaint_transform_point(&matrices[dab_count - 1], x, y, &x_out, &y_out);
                DDI(x_out, y_out, angle - dab_count * rot_angle, bbox_index);
            }

            // Use existing (larger) number of bboxes if it was set (in a snowflake pass)
            num_bboxes_used = MIN(self->num_bboxes, MAX(symm.num_lines, num_bboxes_used));
            break;
        }
        default:
            fprintf(stderr, "Warning: Unhandled symmetry type: %d\n", symm.type);
            break;
        }
    }
    self->num_bboxes_dirtied = MIN(self->num_bboxes, num_bboxes_used);
    return surface_modified;
#undef DDI
}


void get_color (MyPaintSurface *surface, float x, float y,
                  float radius,
                  float * color_r, float * color_g, float * color_b, float * color_a,
                  float paint
                  )
{
    MyPaintTiledSurface *self = (MyPaintTiledSurface *)surface;

    if (radius < 1.0f) radius = 1.0f;
    const float hardness = 0.5f;
    const float softness = 0.5f;
    const float aspect_ratio = 1.0f;
    const float angle = 0.0f;

    float sum_weight, sum_r, sum_g, sum_b, sum_a;
    sum_weight = sum_r = sum_g = sum_b = sum_a = 0.0f;

    // in case we return with an error
    *color_r = 0.0f;
    *color_g = 1.0f;
    *color_b = 0.0f;

    // WARNING: some code duplication with draw_dab

    float r_fringe = radius + 1.0f; // +1 should not be required, only to be sure

    int tx1 = floor(floor(x - r_fringe) / MYPAINT_TILE_SIZE);
    int tx2 = floor(floor(x + r_fringe) / MYPAINT_TILE_SIZE);
    int ty1 = floor(floor(y - r_fringe) / MYPAINT_TILE_SIZE);
    int ty2 = floor(floor(y + r_fringe) / MYPAINT_TILE_SIZE);
    #ifdef _OPENMP
    int tiles_n = (tx2 - tx1) * (ty2 - ty1);
    #endif

    // Calculate the `guaranteed sample` interval and
    // the percentage of pixels to sample for the dab.
    // The basic idea is to have larger intervals and
    // lower percentages for really large dabs, to
    // avoid accumulated rounding errors and heavier
    // calculations.
    //
    // The values are set so that the number of pixels
    // sampled is _bounded_ linearly by the radius.
    //
    // The constant factor 7 is chosen through manual
    // evaluation of results and gives us a total sample
    // rate bounded by '1/(r * 3.5)'
    // Other models may have better properties, some
    // more thinking needed here.
    //
    // For really small radii we'll sample every pixel
    // in the dab to avoid biasing.
    const int sample_interval = radius <= 2.0f ? 1 : (int)(radius * 7);
    const float random_sample_rate = 1.0f / (7 * radius);

    #pragma omp parallel for schedule(static) if(self->threadsafe_tile_requests && tiles_n > 3)
    for (int ty = ty1; ty <= ty2; ty++) {
      for (int tx = tx1; tx <= tx2; tx++) {

        // Flush queued draw_dab operations
        process_tile(self, tx, ty);

        MyPaintTileRequest request_data;
        const int mipmap_level = 0;
        mypaint_tile_request_init(&request_data, mipmap_level, tx, ty, TRUE);

        mypaint_tiled_surface_tile_request_start(self, &request_data);
        uint16_t * rgba_p = request_data.buffer;
        if (!rgba_p) {
          printf("Warning: Unable to get tile!\n");
          break;
        }

        // first, we calculate the mask (opacity for each pixel)
        uint16_t mask[MYPAINT_TILE_SIZE*MYPAINT_TILE_SIZE+2*MYPAINT_TILE_SIZE];

        render_dab_mask(mask,
                        x - tx*MYPAINT_TILE_SIZE,
                        y - ty*MYPAINT_TILE_SIZE,
                        radius,
                        hardness,
                        softness,
                        aspect_ratio, angle
                        );

        // TODO: try atomic operations instead
        #pragma omp critical
        {
        get_color_pixels_accumulate (
          mask, rgba_p, &sum_weight, &sum_r, &sum_g, &sum_b, &sum_a, paint,
          sample_interval, random_sample_rate);
        }

        mypaint_tiled_surface_tile_request_end(self, &request_data);
      }
    }

    assert(sum_weight > 0.0f);
    sum_a /= sum_weight;

    // For legacy sampling, we need to divide
    // by the total after the accumulation.
    if (paint < 0.0) {
        sum_r /= sum_weight;
        sum_g /= sum_weight;
        sum_b /= sum_weight;
    }

    *color_a = CLAMP(sum_a, 0.0f, 1.0f);
    if (sum_a > 0.0f) {
      // Straighten the color channels if using legacy sampling.
      // Clamp to guard against rounding errors.
      const float demul = paint < 0.0 ? sum_a : 1.0;
      *color_r = CLAMP(sum_r / demul, 0.0f, 1.0f);
      *color_g = CLAMP(sum_g / demul, 0.0f, 1.0f);
      *color_b = CLAMP(sum_b / demul, 0.0f, 1.0f);
    } else {
      // it is all transparent, so don't care about the colors
      // (let's make them ugly so bugs will be visible)
      *color_r = 0.0f;
      *color_g = 1.0f;
      *color_b = 0.0f;
    }
}

/**
 * mypaint_tiled_surface_init: (skip)
 *
 * Initialize the surface, passing in implementations of the tile backend.
 * Note: Only intended to be called from subclasses of #MyPaintTiledSurface
 **/
void
mypaint_tiled_surface_init(MyPaintTiledSurface *self,
                           MyPaintTileRequestStartFunction tile_request_start,
                           MyPaintTileRequestEndFunction tile_request_end)
{
    mypaint_surface_init(&self->parent);
    self->parent.draw_dab = draw_dab;
    self->parent.get_color = get_color;
    self->parent.begin_atomic = begin_atomic_default;
    self->parent.end_atomic = end_atomic_default;

    self->tile_request_end = tile_request_end;
    self->tile_request_start = tile_request_start;

    self->tile_size = MYPAINT_TILE_SIZE;
    self->threadsafe_tile_requests = FALSE;

    self->num_bboxes = NUM_BBOXES_DEFAULT;
    self->bboxes = self->default_bboxes;
    memset(self->bboxes, 0, sizeof(MyPaintRectangle) * NUM_BBOXES_DEFAULT);

    self->symmetry_data = mypaint_default_symmetry_data();
    self->operation_queue = operation_queue_new();
}

/**
 * mypaint_tiled_surface_destroy: (skip)
 *
 * Deallocate resources set up by mypaint_tiled_surface_init()
 * Does not free the #MyPaintTiledSurface itself.
 * Note: Only intended to be called from subclasses of #MyPaintTiledSurface
 */
void
mypaint_tiled_surface_destroy(MyPaintTiledSurface *self)
{
    operation_queue_free(self->operation_queue);
    if (self->bboxes != self->default_bboxes) {
      free(self->bboxes);
    }
    mypaint_symmetry_data_destroy(&self->symmetry_data);
}


<|start_filename|>gegl/mypaint-gegl-surface.h<|end_filename|>
#ifndef MYPAINTGEGLSURFACE_H
#define MYPAINTGEGLSURFACE_H

/* libmypaint - The MyPaint Brush Library
 * Copyright (C) 2012 <NAME> <<EMAIL>>
 *
 * 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 <gegl.h>

G_BEGIN_DECLS

#include "mypaint-config.h"
#include "glib/mypaint-gegl-glib.h"
#include "mypaint-tiled-surface.h"

typedef struct MyPaintGeglTiledSurface MyPaintGeglTiledSurface;

GeglBuffer *
mypaint_gegl_tiled_surface_get_buffer(MyPaintGeglTiledSurface *self);

void
mypaint_gegl_tiled_surface_set_buffer(MyPaintGeglTiledSurface *self, GeglBuffer *buffer);

MyPaintSurface *
mypaint_gegl_tiled_surface_interface(MyPaintGeglTiledSurface *self);

MyPaintGeglTiledSurface *
mypaint_gegl_tiled_surface_new(void);

G_END_DECLS

#endif // MYPAINTGEGLSURFACE_H


<|start_filename|>tests/mypaint-utils-stroke-player.c<|end_filename|>
/* libmypaint - The MyPaint Brush Library
 * Copyright (C) 2012 <NAME> <<EMAIL>>
 *
 * 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 <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>

#include "mypaint-utils-stroke-player.h"

int lines_in_string(const char *str) {
    int lines = 0;
    while (*str) {
        if (*str == '\n') {
            lines++;
        }
        str++;
    }
    return lines;
}

typedef struct {
    gboolean valid;
    float time;
    float x;
    float y;
    float pressure;
    float xtilt;
    float ytilt;
    float viewzoom;
    float viewrotation;
    float barrel_rotation;
} MotionEvent;

struct MyPaintUtilsStrokePlayer {
    MyPaintSurface *surface;
    MyPaintBrush *brush;
    MotionEvent *events;
    int current_event_index;
    int number_of_events;
    gboolean transaction_on_stroke; /* If MyPaintBrush::stroke_to should be done between MyPaintSurface::begin_atomic() end_atomic() calls.*/
    float scale;
};

MyPaintUtilsStrokePlayer *
mypaint_utils_stroke_player_new(void)
{
    MyPaintUtilsStrokePlayer *self = (MyPaintUtilsStrokePlayer *)malloc(sizeof(MyPaintUtilsStrokePlayer));

    self->surface = NULL;
    self->brush = NULL;
    self->events = NULL;
    self->number_of_events = 0;
    self->current_event_index = 0;
    self->transaction_on_stroke = TRUE;
    self->scale = 1.0;

    return self;
}

void
mypaint_utils_stroke_player_free(MyPaintUtilsStrokePlayer *self)
{
    if (self->events) {
        free(self->events);
    }
    free(self);
}

void
mypaint_utils_stroke_player_set_brush(MyPaintUtilsStrokePlayer *self, MyPaintBrush *brush)
{
    self->brush = brush;
}

void mypaint_utils_stroke_player_set_surface(MyPaintUtilsStrokePlayer *self, MyPaintSurface *surface)
{
    self->surface = surface;
}

void
mypaint_utils_stroke_player_set_source_data(MyPaintUtilsStrokePlayer *self, const char *data)
{
    self->number_of_events = lines_in_string(data);
    self->events = (MotionEvent *)malloc(sizeof(MotionEvent) * self->number_of_events);

    char * data_copy = strdup(data);

    assert(data_copy);

    char * line = strtok(data_copy, "\n");;
    for (int i=0; i<self->number_of_events; i++) {
        MotionEvent *event = &self->events[i];

        int matches = sscanf(line, "%f %f %f %f %f",
                             &event->time, &event->x, &event->y, &event->pressure, &event->barrel_rotation);
        if (matches != 4) {
            event->valid = FALSE;
            fprintf(stderr, "Error: Unable to parse line '%s'\n", line);
        } else {
            event->valid = TRUE;
        }
        event->xtilt = 0.0;
        event->ytilt = 0.0;
        event->viewzoom = 1.0;
        event->viewrotation = 0.0;
        event->barrel_rotation = 0.0;

        line = strtok(NULL, "\n");
    }

    free(data_copy);

    mypaint_utils_stroke_player_reset(self);
}

gboolean
mypaint_utils_stroke_player_iterate(MyPaintUtilsStrokePlayer *self)
{
    const MotionEvent *event = &self->events[self->current_event_index];
    const int last_event_index = self->current_event_index-1;
    float last_event_time = 0.0;
    if (last_event_index >= 0) {
        const MotionEvent *last_event = &self->events[last_event_index];
        last_event_time = last_event->time;
    }
    const float dtime = event->time - last_event_time;
    const float linear = FALSE;
    if (event->valid) {
        if (self->transaction_on_stroke) {
            mypaint_surface_begin_atomic(self->surface);
        }

        mypaint_brush_stroke_to(self->brush, self->surface,
                                event->x*self->scale, event->y*self->scale,
                                event->pressure,
                                event->xtilt, event->ytilt, dtime,
                                event->viewzoom, event->viewrotation,
                                event->barrel_rotation, linear);

        if (self->transaction_on_stroke) {
            mypaint_surface_end_atomic(self->surface, NULL);
        }
    }
    self->current_event_index++;

    if (self->current_event_index < self->number_of_events) {
        return TRUE;
    } else {
        mypaint_utils_stroke_player_reset(self);
        return FALSE;
    }
}

void
mypaint_utils_stroke_player_reset(MyPaintUtilsStrokePlayer *self)
{
    self->current_event_index = 0;
}

void
mypaint_utils_stroke_player_run_sync(MyPaintUtilsStrokePlayer *self)
{
    while(mypaint_utils_stroke_player_iterate(self)) {
        ;
    }
}

void
mypaint_utils_stroke_player_set_transactions_on_stroke_to(MyPaintUtilsStrokePlayer *self, gboolean value)
{
    self->transaction_on_stroke = value;
}

void
mypaint_utils_stroke_player_set_scale(MyPaintUtilsStrokePlayer *self, float scale)
{
    self->scale = scale;
}
