File Coverage

blib/lib/Color/Model/Munsell/Util.pm

Criterion	Covered	Total	%
statement	190	230	82.6
branch	60	122	49.1
condition	16	33	48.4
subroutine	28	29	96.5
pod	7	7	100.0
total	301	421	71.5

# =============================================================================
package Color::Model::Munsell::Util;
# -----------------------------------------------------------------------------
$Color::Model::Munsell::Util::VERSION = '0.03';
# -----------------------------------------------------------------------------
use warnings;
use strict;

=head1 NAME

Color::Model::Munsell::Util - Utility functions for Color::Model::Munsell

=head1 SYNOPSIS

    use Color::Model::Munsell;
    use Color::Model::Munsell::Util;
    use Color::Model::RGB;

    my $m = Color::Model::Munsell->new("5R 4.5/14");
    printf("Munsell: %s = RGB: #%s\n", $m, Munsell2RGB($m));

=head1 DESCRIPTION

C gives some utility functions for color
conversion from Munsell to CIE xyY, XYZ or RGB, etc.

=cut

# =============================================================================
use Carp qw();
use List::Util qw(first);
use vars qw(@ISA @EXPORT @EXPORT_OK %EXPORT_TAGS $VERSION);
use base qw(Exporter);
@EXPORT    = qw(Munsell2RGB);
@EXPORT_OK = qw(
    huedegree calc_Yc
    Munsell2xyY Munsell2XYZ Munsell2XYZD65 Munsell2rgb
);

use Math::VectorReal;
use Color::Model::Munsell qw(@hue_order %hue_number);
use Color::Model::RGB;

my $_debug = 0;
sub _debug { warn join(" ",@_); }

# =============================================================================
my @_tableAll = ();
my %_tableC   = ();
my $_tableC_loaded = 0;


# -----------------------------------------------------------------------------
sub _table
{
    my $v = shift;
    $v = sprintf('%d',$v);
    _load_table_C();
    if ( $v && $_tableC{$v} ){
        # return same value data sorted by chroma
        return sort { $a->[2] <=> $b->[2] } @{$_tableC{$v}};
    } else {
        return wantarray? @_tableAll: \@_tableAll;
    }
}

sub _load_table_C
{
    unless ( $_tableC_loaded ){
        my $d;
        # append neutral color for convenience
        foreach my $v ( 0.2, 0.4, 0.6, 0.8, 1..10 ){
            $d = [ 'N', $v, 0, _neutralxyY($v) ];
            push @_tableAll, $d;
            push @{$_tableC{$v}}, $d;
        }

        while (){
            tr/\r\n//d;
            s/^ +//g;
            next unless $_;
            next if /^#/;
            my @d = split(/\s+/);
            $d[5] /= 100;                   # rescale Y to (0,1)
            $d = [ @d ];
            push @_tableAll, $d;
            push @{$_tableC{$d[1]}}, $d;    # key is value
        }
        $_tableC_loaded++;
    }
}


# -----------------------------------------------------------------------------

=head1 EXPORT SUBROUTINES

Only subroutine Munsell2RGB is exported by defalut.

=head2 huedegree()

    $n = huedegree( "5R" );

Return degree number of Hue; considered 10.0RP is 0, 10R to be 10, 10YR 20,
..., and ends 9.9RP as 99.9.
If bad formatted hue given, this returns undef.

=cut

sub huedegree
{
    if ( defined($_[0]) && $_[0] =~ /^(\d{1,2}|\d{1,2}\.\d)(R|YR|Y|GY|G|BG|B|PB|P|RP)$/ ){
        my ($n,$c) = ($1,$2);
        if ( $n <= 10 ){
            return $hue_number{$2} * 10 + $1;
        }
    }
    return undef;
}


# =============================================================================
# get xyY with linear interpolation
# -----------------------------------------------------------------------------
my @_Wcxy = (0.310061, 0.316150);

=head2 Munsell2xyY()

Munsell2xyY() returns an array of CIE x, y and Y which are calculated with
linear interpolation from Munsell-xyY table.

    use Color::Model::Munsell::Util qw(Munsell2xyY);

    $m = Color::Model::Munsell->new("2.5G 5.5/10");
    ($x, $y, $Y) = Munsell2xyY($m);

This Munsell-xyY table is from MCSL, R.I.T. which condition is using illuminant
C and the CIE 1931 2 degree observer.

=cut

sub _getHueRange
{
    my ($h1, $h2);
    my ($hueStep, $hueCol);
    if ( $_[0] =~ /^([0-9\.]+)([A-Z]+)$/ ){
        ($hueStep, $hueCol) = ($1, $2);
    }
    if ( $hueStep < 2.5 ){
        my $pre_hC = ($hueCol eq 'R')? 'RP': $hue_order[$hue_number{$hueCol}-1];
        ($h1, $h2) = ( "10$pre_hC", "2.5$hueCol" );
    } elsif ( $hueStep < 5.0 ){
        ($h1, $h2) = ( "2.5$hueCol", "5$hueCol" );
    } elsif ( $hueStep < 7.5 ){
        ($h1, $h2) = ( "5$hueCol", "7.5$hueCol" );
    } else {
        ($h1, $h2) = ( "7.5$hueCol", "10$hueCol" );
    }
    return ($h1, $h2);
}

sub _neutralxyY
{
    my $value = shift;
    return ( @_Wcxy, calc_Yc($value)/100 );  # White Point of illuminant C
}

sub _getVectorOnValuePlane
{
    my ($h0,$v0,$c0) = @_;

    _debug(" * Cheking: $h0 $v0 / $c0\n") if $_debug;

    my ($h1,$h2) = _getHueRange($h0);
      _debug(" - - hue range = [ $h1, $h2 ]\n") if $_debug;

    my @luptable = _table($v0);
    my $eqcheck = first { $_->[0] eq $h0 && $_->[2] eq $c0 } @luptable;
    return vector(@{$eqcheck}[3,4,5]) if ( $eqcheck );  # found just same data!
    my @h1 = grep { $_->[0] eq 'N' or $_->[0] eq $h1 } @luptable;
    my @h2 = grep { $_->[0] eq 'N' or $_->[0] eq $h2 } @luptable;

    # h1,h2 ... hue line ( h1 < h2)
    # c1,c2 ... chroma line (c1 < c2 )
    #
    #  [ on some value-plane ]
    #          c1                 c2
    #          |     P          C/
    #   h2 ----+-----+----------/-- -
    #         B|     |         /     ) 1-l
    #          |     *Z       /     -
    #          |   / |       /        ` l
    #          | /   |      /         /
    #   h1 ----+-----+-----/---     -
    #         A      Q     D
    #          `  n  ^ 1-n '
    #
    #  (  = vector x )
    #   = (1-n) + n
    #   = n
    #   = m + (1-m)
    #       = (m-mn) + mn + (n-mn)
    #       = m + n + mn(--)
    #       = m + n + mn(-)

    my ($dA,$dB,$dC,$dD);   # lookup data
    my ($OA,$OB,$OC,$OD);   # vectors
    my ($c1,$c2);           # chroma values of two points

    # - on line of h1
    $dD = first { $_->[2] >= $c0 } @h1;
    unless ( $dD ){
        # not found, maybe target's chroma is larger than max chroma
        # of lookup table
        return undef;
    }
    $OD = vector( @{$dD}[3,4,5] );
    ($c1,$c2) = ($dD->[2]-2, $dD->[2]);

    $dA = first { $_->[2] == $c1 } @h1;
    $OA = vector( @{$dA}[3,4,5] );

    if ( $_debug ){
        _debug(" - - - $h1> OA = ($dA->[0] $dA->[1]/ $dA->[2]) \n");
        _debug(" - - - $h1> OD = ($dD->[0] $dD->[1]/ $dD->[2]) \n");
    }

    # - on line of h2
    $dC = first { $_->[2] >= $c0 } @h2;
    unless ( $dC ){
        # not found on lookup table
        return undef;
    }
    $OC = vector( @{$dC}[3,4,5] );
    ($c1,$c2) = ($dC->[2]-2, $dC->[2]);

    $dB = first { $_->[2] == $c1 } @h2;
    $OB = vector( @{$dB}[3,4,5] );

    if ( $_debug ){
        _debug(" - - - $h2> OB = ($dB->[0] $dB->[1]/ $dB->[2]) \n");
        _debug(" - - - $h2> OC = ($dC->[0] $dC->[1]/ $dC->[2]) \n");
    }

    # calculate ratio
    my $n = ($c0 - $c1) / ($c2 - $c1);              # ratio between chromas
    my ($h0n, $h1n, $h2n) = map { huedegree($_) } ($h0, $h1, $h2);
    my $l = ($h0n - $h1n) / ($h2n - $h1n);    # ratio between hues
    if ( $_debug ){
        _debug(" - - - Ratio chroma; ($c0 - $1) / ($c2 - $c1) = $n\n");
        _debug(" - - - Ratio hue;    ($h0n - $h1n) / ($h2n - $h1n) = $l\n");
    }

    my $AB = $OB - $OA;
    my $AD = $OD - $OA;
    my $BC = $OC - $OB;

    # return value is a vector  with Math::RealVector;
    my $AZ = $AB*$l + $AD*$n+ ($BC-$AD)*$l*$n;
    my $OZ = $OA + $AZ;

    _debug(map { sprintf("%s = (%.4f, %.4f, %.4f)\n", $_, eval"\$$_->array()") } qw(OA OB OC OD OZ)) if $_debug;
    return $OZ;
}

sub Munsell2xyY
{
    my $m = shift;
    unless ( defined($m) && ref($m) eq 'Color::Model::Munsell' ){
        Carp::croak("Color::Model::Munsell object is not given");
    }
    _debug("<< Target: $m >>\n") if $_debug;
    my $v  = $m->value;

    if ( $m->isneutral ){
        my @ret = _neutralxyY($v);
        _debug(sprintf("xyY = (%.6f, %.6f, %.6f)\n",@ret)) if $_debug;
        return @ret;
    }

    _load_table_C();

    my ($LZ, $UZ, $Z);
    # get lower value-plane
    my $lv = ($v<1.0)? (int($v*10/2)*2/10): int($v);
    _debug(" - nearest lower value of $m is $lv\n") if $_debug;
    if ( $lv == 0 ){
        $LZ = vector( _neutralxyY($v) );
    } else {
        $LZ = _getVectorOnValuePlane( $m->hue, $lv, $m->chroma );
        unless ( defined($LZ) ){
            Carp::croak("$m is out of calculatable color space.");
        }
    }
    # get upper value-plane
    my $uv = ($v<1.0)? (int($v*10/2)*2/10+0.2): int($v+1);
    _debug(" - nearest upper value of $m is $uv\n") if $_debug;
    $UZ = _getVectorOnValuePlane( $m->hue, $uv, $m->chroma );
    unless ( defined($UZ) ){
        Carp::croak("$m is out of calculatable color space.");
    }

    _debug(" - value range = [ $lv, $uv ]\n") if $_debug;

    my $s  = ($v - $lv) / ($uv - $lv ); # ratio between values;
    _debug(" - value ratio = ($v - $lv) / ($uv - $lv) = $s\n") if $_debug;

    $Z = $UZ*$s + $LZ*(1-$s);
    _debug(sprintf("xyY = (%.6f, %.6f, %.6f)\n",$Z->array())) if $_debug;

    return $Z->array();   # ( xc, yc , Yc );
}
# -----------------------------------------------------------------------------


# -----------------------------------------------------------------------------
# multiply matrix * vector
# -----------------------------------------------------------------------------
sub _mmult
{
    my ($m,$v) = @_;
    return (
        $$m[0][0]*$$v[0] + $$m[0][1]*$$v[1] + $$m[0][2]*$$v[2],
        $$m[1][0]*$$v[0] + $$m[1][1]*$$v[1] + $$m[1][2]*$$v[2],
        $$m[2][0]*$$v[0] + $$m[2][1]*$$v[1] + $$m[2][2]*$$v[2]
    );
}


# =============================================================================
# get XYZ from Munsell with xyY
# -----------------------------------------------------------------------------
my %_matCAdaptD65 = (
    'XYZ' => [
        [ 0.9691356, 0.0000000, 0.0000000],
        [ 0.0000000, 1.0000000, 0.0000000],
        [ 0.0000000, 0.0000000, 0.9209267]],
    'Bradford' => [
        [ 0.9904476,-0.0071683,-0.0116156],
        [-0.0123712, 1.0155950,-0.0029282],
        [-0.0035635, 0.0067697, 0.9181569]],
    'vonKries' => [
        [ 0.9972812,-0.0093756,-0.0154171],
        [-0.0010298, 1.0007636, 0.0002084],
        [ 0.0000000, 0.0000000, 0.9209267]]
);

=head2 Munsell2XYZ()

=head2 Munsell2XYZD65()

    Munsell2XYZ( $m )
    Munsell2XYZD65( $m [, "ChromaticAdaptType" ] )

Munsell2XYZ() and Munsell2XYZD65() returns an array of CIE X,Y and Z.

    $m = Color::Model::Munsell->new("7.5B 6/10");
    printf("%s -> CIE XYZ (%.f, %.f, %.f)", $m, Munsell2XYZ($m));

    # perform cromatic adaptation from C to D65
    printf("%s -> CIE XYZ (%.f, %.f, %.f) via Chromatic Adaptation", $m, Munsell2XYZD65($m));

Munsell2XYZ() simply calculates XYZ from xyY. And Munsell2XYZD65() calculates
them with chromatic adaptatation to illuminant D65. Adaptation type must be
"XYZ", "vonKries", "Bradford" or "None". If Omitted, "Bradford" is used.
Specifying "None" is same as calling Munsell2XYZ().

=cut

sub Munsell2XYZ
{
    my $m = shift;                          # Color::Model::Munsell object
    unless ( defined($m) && ref($m) eq 'Color::Model::Munsell' ){
        Carp::croak("Munsell2XYZ() needs Color::Model::Munsell object.");
    }

    my ($x,$y,$Y) = Munsell2xyY($m);
    my $X = ($x * $Y) / $y;                 # X = (Y/y) * x
    my $Z = ( (1 - $x - $y) * $Y ) / $y;    # Z = ( (1-x-y)/y ) * Y
    return ($X, $Y, $Z);
}

sub Munsell2XYZD65
{
    my $m = shift;                          # Color::Model::Munsell object
    unless ( defined($m) && ref($m) eq 'Color::Model::Munsell' ){
        Carp::croak("Munsell2XYZD65() needs Color::Model::Munsell object.");
    }

    my $atype = shift || 'Bradford';        # Chromatic Adaptation type
    unless ( $atype =~ /^(None|XYZ|Bradford|vonKries)$/ ){
        Carp::croak(qq(Chromatic Adaptation must be "XYZ", "Bradford" or "vonKries"));
    }
      _debug(" - Chromatic Adaptation = $atype\n") if $_debug;

    my @XYZ;
    if ( $atype ne 'None' ){
        @XYZ = _mmult($_matCAdaptD65{$atype}, [ Munsell2XYZ($m) ]);
    } else {
        @XYZ = Munsell2XYZ($m);
    }
    _debug(sprintf("XYZ = (%.6f, %.6f, %.6F)\n",@XYZ)) if $_debug;
    return @XYZ;
}


# =============================================================================
# get RGB from XYZ
# -----------------------------------------------------------------------------
my %_matXYZ2rgb = (
    'AdobeRGB' => [     # D65
        [ 2.0413690, -0.5649464, -0.3446944],
        [-0.9692660,  1.8760108,  0.0415560],
        [ 0.0134474, -0.1183897,  1.0154096]],
    'AppleRGB' => [     # D65
        [ 2.9515373, -1.2894116, -0.4738445],
        [-1.0851093,  1.9908566,  0.0372026],
        [ 0.0854934, -0.2694964,  1.0912975]],
    'PAL' => [          # D65
        [ 3.0628971, -1.3931791, -0.4757517],
        [-0.9692660,  1.8760108,  0.0415560],
        [ 0.0678775, -0.2288548,  1.0693490]],
    'sRGB' => [         # D65
        [ 3.2404542, -1.5371385, -0.4985314],
        [-0.9692660,  1.8760108,  0.0415560],
        [ 0.0556434, -0.2040259,  1.0572252]],
    'NTSC' => [         # C
        [ 1.9099961, -0.5324542, -0.2882091],
        [-0.9846663,  1.9991710, -0.0283082],
        [ 0.0583056, -0.1183781,  0.8975535]],
);


=head2 Munsell2rgb()

=head2 Munsell2RGB()

    Munsell2rgb( $m [, "RGBModel" [, "ChromaticAdaptType" ]] )
    Munsell2RGB( $m [, "RGBModel" [, "ChromaticAdaptType" [, $gamma ]]] )

Munsell2rgb() returns an array of R, G and B values which calculated from
XYZ with transformation matrix.
And Munsell2RGB() returns RGB values with C object which applied gamma value.

    $m = Color::Model::Munsell->new("7PB 2.5/3");
    printf("%s -> RGB %s of sRGB", $m, Munsell2RGB($m));
    printf("%s -> RGB %s of AdobeRGB(1998)", $m, Munsell2RGB($m,"AdobeRGB"), 2.2);

RGBModel must be "sRGB", "AdobeRGB" that means Adobe RGB(1998), "AppleRGB"
or "NTSC". A gamma value will be used when RGB model is not sRGB.

=cut

sub Munsell2rgb
{
    my $m = shift;                      # Color::Model::Munsell object
    unless ( defined($m) && ref($m) eq 'Color::Model::Munsell' ){
        Carp::croak("Munsell2rgb() needs Color::Model::Munsell object.");
    }
    my $rgbtype = shift || 'sRGB';      # RGB type
    my $atype   = shift || 'Bradford';  # Chromatic Adaptation type
    unless ( $rgbtype =~ /^(AdobeRGB|AppleRGB|PAL|sRGB|NTSC)$/ ){
        Carp::croak(qq(RGB type must be "AdobeRGB", "AppleRGB", "PAL", "NTSC" or "sRGB"));
    }
      _debug(" - RGB = $rgbtype\n") if $_debug;

    my @XYZ;
    if ( $rgbtype eq 'NTSC' ){
        @XYZ = Munsell2XYZ($m);
    } else {
        @XYZ = Munsell2XYZD65($m, $atype);
    }

    my @rgb = _mmult($_matXYZ2rgb{$rgbtype}, \@XYZ);
    _debug(sprintf("rgb = (%.6f, %.6f, %.6f)\n",@rgb)) if $_debug;
    return @rgb;
}

sub Munsell2RGB
{
    my $m = shift;                      # Color::Model::Munsell object
    unless ( defined($m) && ref($m) eq 'Color::Model::Munsell' ){
        Carp::croak("Munsell2RGB() needs Color::Model::Munsell object.");
    }
    my $rgbtype = shift || 'sRGB';      # RGB type
    my $atype   = shift || 'Bradford';  # Chromatic Adaptation type
    my $gamma   = shift || 2.2;         # Gamma value (uses when not sRGB)

    my @RGB;
    if ( $rgbtype eq 'sRGB' ){
        @RGB = map {
                $_ = $_>0.0031308?
                    (1.055 * ($_ ** (1/2.4)) - 0.055):
                    12.92 * $_;
                $_ = $_<0 ? 0: ( $_>1? 1: $_);
        } Munsell2rgb($m, 'sRGB', $atype);
    } else {
        @RGB = map {
                $_ = $_ ** (1/$gamma);
                $_ = $_<0 ? 0: ( $_>1? 1: $_);
        } Munsell2rgb($m, $rgbtype, $atype);
    }
    _debug(sprintf("RGB = (%.6f, %.6f, %.6f)\n",@RGB)) if $_debug;
    return Color::Model::RGB->new(@RGB);
}


# =============================================================================
# Calculation Y of illuminant C
# -----------------------------------------------------------------------------

=head2 calc_Yc()

    calc_Yc( $m )

calc_Yc() calculates Y value of illuminant C and 2 degree observer from Munsell
value with approximate calculation. Argument must be a I object or Munsell value.

    $m = Color::Model::Munsell->new("5R 4/14")
    printf "%.4f", calc_Yc( $m );
    printf "%.4f", calc_Yc( 7.5 );

=cut

use Scalar::Util qw(looks_like_number);

sub calc_Yc
{
#   Yc = 1.1913*V - 0.22532*V^2 + 0.23351*V^3 - 0.020483*V^4 + 0.00081936*V^5
#   (ref. JIS Z 8721(1993) )

    my $m = shift;
    my $v = undef;
    if ( defined($m) ){
        if ( ref($m) eq 'Color::Model::Munsell' ){
            $v = $m->value;
        } elsif ( looks_like_number($m) ){
            $v = $m;
        }
    }
    unless ( defined($v) ){
        Carp::croak("Bad argument. calc_Yc() needs a number");
    }
    my $v2 = $v*$v;
    my $v3 = $v*$v2;
    my $v4 = $v2*$v2;
    my $v5 = $v2*$v3;
    return 1.1913*$v - 0.22532*$v2 + 0.23351*$v3 - 0.020483*$v4 + 0.00081936*$v5;
}


=head1 AUTHOR

Takahiro Onodera, C<<  >>

=head1 BUGS

Please report any bugs or feature requests to C, or through
the web interface at L.  I will be notified, and then you'll
automatically be notified of progress on your bug as I make changes.

=head1 SEE ALSO

L, L

=head1 REFERENCES

Munsell Color Science Laboratory, R.I.T - L - Munsell-xyY data are from this site.

BruceLindbloom.com - L - Chromatic adaptation matrixes, transformation matrixes and important knowledge are from this site.

Japanese Industrial Standards(JIS) JIS Z 8721(1993)

=head1 LICENSE AND COPYRIGHT

Copyright 2010 Takahiro Onodera.

This program is free software; you can redistribute it and/or modify it
under the terms of either: the GNU General Public License as published
by the Free Software Foundation; or the Artistic License.

=cut

1; # End of Color::Model::Munsell::Util
__DATA__