Based on Parametric fuzzy sets for automatic color naming

automatic_fuzzy_color_naming.rs

/* DOES NOT RETURN CORRECT RESULTS,
 * current implementation violates the unity-sum constraint
 * (sum of all memberships sometimes exceeds 1) and skips
 * achromatic category memberships. */

extern crate nalgebra;

use nalgebra::{U1, U3, Matrix, MatrixArray};
use std::ops::{Mul, Div};

type Matrix3x3 = Matrix<f32, U3, U3, MatrixArray<f32, U3, U3>>;
type RowMatrix3 = Matrix<f32, U1, U3, MatrixArray<f32, U1, U3>>;
type ColumnMatrix3 = Matrix<f32, U3, U1, MatrixArray<f32, U3, U1>>;

type Point3d = ColumnMatrix3; // (a, b, 1)^T
type Point2d = (f32, f32); // (a, b)
type Translation2d = (f32, f32); // (t_a, t_b)
type Rotation2d = (Rotation, Rotation); // (alpha_a, alpha_b)
type Slope2d = (Slope, Slope); // (beta_a, beta_b)

type Semiaxes = (f32, f32);

type Rotation = f32; // alpha
type Slope = f32; // beta


macro_rules! scalar {
    ($one_by_one_matrix: expr) => { $one_by_one_matrix.as_slice()[0] }
}

fn sigmoid_2d(p: Point3d, t: Translation2d, a: Rotation, b: Slope, axis_vector: RowMatrix3) -> f32 {
    let transform = axis_vector
        .map(|e| e * -b)
        .mul(rotation_matrix(a))
        .mul(translation_matrix(t));

    (1.0 + scalar!(transform.mul(p)).exp()).recip()
}

fn double_sigmoid((p_a, p_b): Point2d, t: Translation2d, (a_a, a_b): Rotation2d, (b_a, b_b): Slope2d) -> f32 {
    let point_3d = ColumnMatrix3::new(p_a, p_b, 1.0);
    let axis_a = RowMatrix3::new(1.0, 0.0, 0.0);
    let axis_b = RowMatrix3::new(0.0, 1.0, 0.0);

    sigmoid_2d(point_3d, t, a_b, b_b, axis_a) * sigmoid_2d(point_3d, t, a_a, b_a, axis_b)
}

fn elliptic_sigmoid((p_a, p_b): Point2d, t: Translation2d, (e_a, e_b): Semiaxes, phi: Rotation, b_e: Slope) -> f32 {
    let point_3d = ColumnMatrix3::new(p_a, p_b, 1.0);
    let axis_a = RowMatrix3::new(1.0, 0.0, 0.0);
    let axis_b = RowMatrix3::new(0.0, 1.0, 0.0);
    
    let curve = {
        let transform_matrix = rotation_matrix(phi).mul(translation_matrix(t));
        (scalar!(axis_a.mul(transform_matrix).mul(point_3d)).div(e_a).powi(2) +
         scalar!(axis_b.mul(transform_matrix).mul(point_3d)).div(e_b).powi(2) - 1.0)
    };

    (1.0 + (-b_e * curve).exp()).recip()
}

fn triple_sigmoid_with_elliptical_center(p: Point2d, t: Translation2d, /* DS */ a2d: Rotation2d, b2d: Slope2d, /* ES */ e: Semiaxes, phi: Rotation, b_e: Slope) -> f32 {
    double_sigmoid(p, t, a2d, b2d) * elliptic_sigmoid(p, t, e, phi, b_e)
}

type ColorSample = (f32, f32, f32); // (l, a, b)
type Category = usize;

#[derive(Debug)]
struct PlaneParameters {
    t: Translation2d,
    /* ES params */
    e: Semiaxes,
    phi: Rotation,
    b_e: Slope,

    /* chromatic categories: red, orange, brown, yellow, green, blue, purple, pink */
    categories: [Option<(Rotation2d, Slope2d)>; 8]
}

const PLANES: [PlaneParameters; 6] = [
    PlaneParameters {
        t: (0.42, 0.25),
        e: (5.89, 7.47),
        phi: 2.32,
        b_e: 9.84,
        categories: [
            Some(((-2.24, -56.55), (0.90, 1.72))),
            None,
            Some(((33.45, 14.56), (1.72, 0.84))),
            None,
            Some(((104.56, 134.59), (0.84, 1.95))),
            Some(((224.59, -147.15), (1.95, 1.01))),
            Some(((-57.15, -92.24), (1.01, 0.90))),
            None
        ]
    },
    PlaneParameters {
        t: (0.23, 0.66),
        e: (6.46, 7.87),
        phi: 17.59,
        b_e: 6.03,
        categories: [
            Some(((2.21, -48.81), (0.52, 5.00))),
            None,
            Some(((41.19, 6.87), (5.00, 0.69))),
            None,
            Some(((96.87, 120.46), (0.69, 0.96))),
            Some(((210.46, -148.48), (0.96, 0.92))),
            Some(((-58.48, -105.72), (0.92, 1.10))),
            Some(((-15.72, -87.79), (1.10, 0.52)))
        ]
    },
    PlaneParameters {
        t: (-0.12, 0.52),
        e: (5.38, 6.98),
        phi: 19.58,
        b_e: 6.81,
        categories: [
            Some(((13.57, -45.54), (1.00, 0.57))),
            Some(((44.45, -28.76), (0.57, 0.52))),
            Some(((61.24, 6.65), (0.52, 0.84))),
            None,
            Some(((96.65, 109.38), (0.84, 0.60))),
            Some(((199.38, -148.24), (0.60, 0.80))),
            Some(((-58.24, -112.63), (0.80, 0.62))),
            Some(((-22.63, -76.43), (0.62, 1.00)))
        ]
    },
    PlaneParameters {
        t: (-0.47, 1.02),
        e: (5.99, 7.51),
        phi: 23.92,
        b_e: 7.76,
        categories: [
            Some(((26.70, -56.88), (0.91, 0.76))),
            Some(((33.12, -9.90), (0.76, 0.48))),
            None,
            Some(((80.10, 5.63), (0.48, 0.73))),
            Some(((95.63, 108.14), (0.73, 0.64))),
            Some(((198.14, -148.59), (0.64, 0.76))),
            Some(((-58.59, -123.68), (0.76, 5.00))),
            Some(((-33.68, -63.30), (5.00, 0.91)))
        ]
    },
    PlaneParameters {
        t: (-0.57, 1.16),
        e: (5.37, 6.90),
        phi: 24.75,
        b_e: 100.00,
        categories: [
            None,
            Some(((25.75, -15.85), (2.00, 0.84))),
            None,
            Some(((74.15, 12.27), (0.84, 0.86))),
            Some(((102.27, 98.57), (0.86, 0.74))),
            Some(((188.57, -150.83), (0.74, 0.47))),
            Some(((-60.83, -122.55), (0.47, 1.74))),
            Some(((-32.55, -64.25), (1.74, 2.00)))
        ]
    },
    PlaneParameters {
        t: (-1.26, 1.81),
        e: (6.04, 7.39),
        phi: -1.19,
        b_e: 100.00,
        categories: [
            None,
            Some(((25.74, -17.56), (1.03, 0.79))),
            None,
            Some(((72.44, 16.24), (0.79, 0.96))),
            Some(((106.24, 100.05), (0.96, 0.90))),
            Some(((190.05, -149.43), (0.90, 0.60))),
            Some(((-59.43, -122.37), (0.60, 1.93))),
            Some(((-32.37, -64.26), (1.93, 1.03)))
        ]
    }
];

pub fn membership(k: Category, (l, a, b): ColorSample) -> f32 {
    let plane = {
        if l <= 31.0 { 0 }
        else if l <= 41.0 { 1 }
        else if l <= 51.0 { 2 }
        else if l <= 66.0 { 3 }
        else if l <= 76.0 { 4 }
        else { 5 }
    };
    let PlaneParameters { t, e, phi, b_e, categories } = PLANES[plane];
    
    match categories[k] {
        Some((a2d, b2d)) =>
            triple_sigmoid_with_elliptical_center((a, b), t, a2d, b2d, e, phi, b_e),
        None =>
            0.0
    }
}

fn translation_matrix((t_a, t_b): Translation2d) -> Matrix3x3 {
    Matrix3x3::new(1.0, 0.0, -t_a,
                   0.0, 1.0, -t_b,
                   0.0, 0.0,  1.0)
}

fn rotation_matrix(a: Rotation) -> Matrix3x3 {
    Matrix3x3::new( a.cos(), a.sin(), 0.0,
                   -a.sin(), a.cos(), 0.0,
                        0.0,     0.0, 1.0)
}

#[cfg(test)]
mod tests {
    use membership;

    #[test]
    fn it_works() {
        for k in 0..7 {
            println!("k = {}, m = {}", k, membership(k, (84.22, 15.75, 13.75)));
        }
    }
}

reimplemented.rb

# Same algorithm, implemented in Ruby,
# which has the same problems as the Rust version.

# Example:
# Here are the results for #00dd6d,
# a light, highly saturated green:
#
# orange: 0.9999998431800354,
# blue: 0.999999999941878,
# ...
# sum = 1.9999998431219137
#
# Bluish orange is not a color we could perceive, now is it?

require 'matrix'

# Based on https://www.osapublishing.org/josaa/abstract.cfm?uri=josaa-25-10-2582
# (pdf link is in the right corner)

AXIS_A = Matrix.row_vector [1.0, 0.0, 0.0]
AXIS_B = Matrix.row_vector [0.0, 1.0, 0.0]

Tup = Struct.new(:a, :b)

def scalar(single_el_matrix)
  single_el_matrix[0, 0]
end

# Eq. 7
def S(p, t_tup, α, β, axis)
  transform = axis.map { |e| e * -β } * translation_m(t_tup) * rotation_m(α)
  1.0 / (1.0 + Math.exp(scalar(transform * p)))
end

# Eq. 8
def translation_m(t_tup)
  Matrix[ [1.0, 0.0, -t_tup.a],
          [0.0, 1.0, -t_tup.b],
          [0.0, 0.0, 1.0] ]
end

# Eq. 8
def rotation_m(α)
  Matrix[ [Math.cos(α), Math.sin(α), 0.0],
          [-Math.sin(α), Math.cos(α), 0.0],
          [0.0, 0.0, 1.0] ]
end

# Eq. 9
def DS(p, t_tup, α_tup, β_tup)
  S(p, t_tup, α_tup.b, β_tup.b, AXIS_A) * S(p, t_tup, α_tup.a, β_tup.a, AXIS_B)
end

# Eq. 10
def ES(p, t_tup, e_tup, ϕ, βe)
  curve = (
    (scalar(AXIS_A * rotation_m(ϕ) * translation_m(t_tup) * p) / e_tup.a) ** 2 +
    (scalar(AXIS_B * rotation_m(ϕ) * translation_m(t_tup) * p) / e_tup.b) ** 2 - 1.0)
  
  1.0 / (1.0 + Math.exp(-βe * curve))
end

# Eq. 11
def TSE(p, t_tup, α_tup, β_tup, e_tup, ϕ, βe)
  DS(p, t_tup, α_tup, β_tup) * ES(p, t_tup, e_tup, ϕ, βe)
end

def membership(category, l, a, b)
  # p = (a, b, 1)^T
  p = Matrix.column_vector [a, b, 1.0]
  
  # See Parameter Estimation section of
  # https://hal.inria.fr/hal-00640930/document
  plane = if l <= 31
      0
    elsif l <= 41
      1
    elsif l <= 51
      2
    elsif l <= 66
      3
    elsif l <= 76
      4
    else
      5
    end
  
  params = CHROMATICITY_PLANES[plane]
  
  return 0.0 if params == nil
  
  TSE(p, params[:t],
         params[:categories][category][:α],
         params[:categories][category][:β],
         params[:e],
         params[:ϕ],
         params[:βe])
end

CHROMATICITY_PLANES = [
  {},
  {},
  {},
  {},
  {},
  {t: Tup.new(-1.26, 1.81),
   e: Tup.new(6.04, 7.39),
   ϕ: -1.19,
   βe: 100.0,
   categories: {
     orange: { α: Tup.new(25.74, -17.56),
               β: Tup.new(1.03, 0.79) },
     yellow: { α: Tup.new(72.44, 16.24),
               β: Tup.new(0.79, 0.96) },
     green:  { α: Tup.new(106.24, 100.05),
               β: Tup.new(0.96, 0.90) },
     blue:   { α: Tup.new(190.05, -149.43),
               β: Tup.new(0.96, 0.90) },
     purple: { α: Tup.new(-59.43, -122.37),
               β: Tup.new(0.60, 1.93) },
     pink:   { α: Tup.new(-32.37, -64.26),
               β: Tup.new(1.93, 1.03) }
  }}
]

# A random color with L component > 76 in CIELAB colorspace
# (I haven't inserted the rest of the planes here, but my Rust
#  version had them and didn't work either)
cl, ca, cb = [77, -75, 41]

# Memberships for categories not present in the plane is assumed to be 0 
categories = CHROMATICITY_PLANES.last[:categories].keys
memberships = categories.map { |k| membership(k, cl, ca, cb) }

memberships.zip(categories).each { |m, k| puts "#{k}: #{m}" }

# This doesn't fulfill the unity-sum constraint.
# Sum of all memberships should be equal to 0,
# and as I didn't include achromatic categories,
# it could even be less than that, but not more.
sum = memberships.inject(0) {|sum, m| sum + m }