L*a*b*色空間上で色差を求める CIE DE 2000（ perl ）

cie_de_2000.pl

# 参考：http://www.brucelindbloom.com/index.html?Eqn_DeltaE_CIE2000.html

=comment
 $labx = {
   L => value,
   a => value,
   b => value,
 };
=cut

sub cie_de_2000 {
    my ($lab1, $lab2) = @_;
    my ($L1, $a1, $b1, $L2, $a2, $b2) = (
        $lab1->{L}, $lab1->{a}, $lab1->{b},
        $lab2->{L}, $lab2->{a}, $lab2->{b},
        );

    my $Ldash_ = ( $L1 + $L2 ) / 2;
    my $C1 = sqrt( $a1 ** 2 + $b1 ** 2 );
    my $C2 = sqrt( $a2 ** 2 + $b2 ** 2 );
    my $C_ = ( $C1 + $C2 ) / 2;
    my $G = ( 1 - sqrt( $C_ ** 7 / ($C_ ** 7 + 25 ** 7 ) ) ) / 2;
    my $a1dash = $a1 * ( 1 + $G );
    my $a2dash = $a2 * ( 1 + $G );
    my $C1dash = sqrt( $a1dash ** 2 + $b1 ** 2);
    my $C2dash = sqrt( $a2dash ** 2 + $b2 ** 2);
    my $Cdash_ = ( $C1dash + $C2dash )  / 2;

    my $h1dash = rad2deg( atan2( $b1, $a1dash ) );
    if ( 0 > $h1dash ) {
        $h1dash += 360;
    }
    my $h2dash = rad2deg( atan2( $b2 ,$a2dash ) );
    if ( 0 > $h2dash ) {
        $h2dash += 360;
    }
    my $Hdash_;
    if ( abs( $h1dash - $h2dash ) > 180 ) {
        $Hdash_ = ( $h1dash + $h2dash + 360) / 2;
    }
    else {
        $Hdash_ = ( $h1dash + $h2dash ) / 2;
    }

    my $T = 1
        - 0.17 * cos( $Hdash_ - 30 )
        + 0.24 * cos( 2 * $Hdash_)
        + 0.32 * cos( 3 * $Hdash_ + 6 )
        - 0.20 * cos( 4 * $Hdash_ - 63 )
        ;

    my $dhdash;
    if ( abs( $h2dash - $h1dash) <= 180 ) {
        $dhdash = $h2dash - $h1dash;
    }
    elsif ( ( abs( $h2dash - $h1dash) > 180 ) && ( $h2dash <= $h1dash ) ) {
        $dhdash = $h2dash - $h1dash + 360;
    }
    else {
        $dhdash = $h2dash - $h1dash - 360;
    }

    my $dLdash = $L2 - $L1;
    my $dCdash = $C2dash - $C1dash;
    my $dHdash = 2 * sqrt($C1dash * $C2dash) * sin( deg2rad( $dhdash / 2 ) );

    my $S_L = 1 + ( ( 0.015 * ($Ldash_ - 50) ** 2) / sqrt(20 + ($Ldash_ - 50) ** 2) );
    my $S_C = 1 + 0.045 * $Cdash_;
    my $S_H = 1 + 0.015 * $Cdash_ * $T;

    my $dtheta = 30 * exp( -( ($Hdash_ - 275) / 25) ** 2);

    my $R_C = 2 * sqrt( $Cdash_ ** 7 / ($Cdash_ ** 7 + 25 ** 7));
    my $R_T = -$R_C * sin( deg2rad( 2 * $dtheta ) );

    my $K_L = 1;
    my $K_C = 1;
    my $K_H = 1;

    my $dE = sqrt(
          ($dLdash / ($K_L * $S_L) ) ** 2
        + ($dCdash / ($K_C * $S_C) ) ** 2
        + ($dHdash / ($K_H * $S_H) ) ** 2
        + $R_T * ($dCdash / ($K_C * $S_C) ) * ($dHdash / ($K_H * $S_H) )
        );
    return $dE;
}