helpers.c

/* libmypaint - The MyPaint Brush Library
 * Copyright (C) 2007-2008 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.
 */

#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