PDL::Transform::Color - Convert between color systems

gistfile1.pl

=head1 NAME

PDL::Transform::Color - Convert between color systems

=head1 SYNOPSIS

 # load an image (in this case the cartographic demo)
 use PDL::Transform::Cartography; $rgb = earth_image();

 # convert the image to CMYK
 $cmyk = $rgb->apply(t_cmyk());

=head1 DESCRIPTION

PDL::Transform::Color bundles several important color system
transformations.  Unlike most transforms, PDL::Transform::Color
transforms are optimized to work with color images.  Color images
typically have many pixels in the spatial direction and 3 or 4 pixels
in the color-index direction, but are sometimes stored with the
color index in either the 0 or 2 dim.  The PDL::Transform::Color
transformations should examine the first few dims of their input 
and, if only one suitably sized dim is found, they should transform
that dim rather than the 0 dim.

Conversions are relative to the RGB system, so (e.g.) C<t_cmyk>
converts RGB to CMYK and C<!t_cmyk> converts CMYK to RGB.

All values are promoted to floating point before conversion, to
avoid quantization problems associated with fixed-point/integer
arithmetic in the general case.

At present, the color representations all utilize simple linear 
theory -- no provision is made (other than gamma-encoding) for 
nonlinearities in visual perception, nor for the various issues that
arise with the particular spectral response from individual pigments.

=head1 STANDARD OPTIONS

Several options are standard.  They are:

=over 3 

=item * 'gamma' (or 'g') - gamma of the RGB space (default 1.0)

This is the gamma correction factor used to get physical values from
the RGB values in the RGB space.  Conversion is performed in a gamma=1
space -- i.e. if you specify a gamma to the forward transform, the
input RGB values are assumed to be gamma encoded, and are decoded to
linear physical values before processing.

=item * 'max' - max in-gamut value of the RGB space (default 1.0 or 255)

Some of the conversions (notably CMYK) require a range of the RGB
space to define their gamut.  The minimum is always taken to be 0.

=back


=head1 AUTHOR

Copyright 2012, Craig DeForest (deforest@boulder.swri.edu).  
This module may be modified and distributed under the same
terms as PDL itself.  The module comes with NO WARRANTY.

=head1 FUNCTIONS

This module defines and exports transform constructors ('t_<foo>') only.

=cut

use PDL::Transform;
package PDL::Transform::Color;
use PDL::Core ':Internal'; # Load 'topdl'

@ISA = ('Exporter','PDL::Transform');
$VERSION = "0.2";

BEGIN {
  use Exporter();
  @EXPORT_OK = qw(t_cmyk t_hsv);
  @EXPORT = @EXPORT_OK;
  %EXPORT_TAGS = (Func=>[@EXPORT_OK]);
}

use PDL;
use PDL::Transform;
use PDL::NiceSlice;

use Carp;

##############################
# Steal _opt from PDL::Transform.
*PDL::Transform::Color::_opt = \&PDL::Transform::_opt;

##############################
# Enable our own stringifier 
use overload '""' => \&_strval;
sub _strval {
  my($me) = shift;
  $me->stringify();
}

use strict;

use PDL::Constants;

sub _new { new('PDL::Transform::Color',@_); } # not exported
sub new {
    my $class = shift;
    my $opt = shift;

    my $me = PDL::Transform::new($class);

    $me->{name} = "generic color transform";
    $me->{idim} = 0;
    $me->{odim} = 0;
    $me->{itype} = ['red','green','blue'];
    $me->{iunit} = ['brightness','brightness','brightness'];
    $me->{func} = \*PDL::Transform::_identity;
    $me->{inv} = \*PDL::Transform::_identity;
    $me->{params} = {};

    # Parse standard options
    $me->{options}->{gamma} = _opt($opt,["gamma","g"]);
    $me->{options}->{max}   = _opt($opt,["max","m"]);

    bless $me,$class;
}

# Find the correct active and put it in front; return the dim number where it goes
# Also, promote integer types to float for internal work...
sub _rectify {
    my $input = shift;

    if($input->type !~ m/(double|float)/) {
	$input = float $input;
    }

    my $pos = shift;
    if(defined($pos)) {
	if($pos) {
	    return ($input->mv($pos,0),$pos);
	} else {
	    return ($input,undef);
	}
    }
    my $dims = pdl($input->dims);
    my $okdims = which($dims <= 5);
    if($okdims->nelem == 0) {
	die "PDL::Transform::Color: couldn't find an appropriate dim (size <= 5) for color\n  vectors. Specify position explicitly in the transform constructor.";
    }
    $pos = $okdims->at(0);
    return ($input->mv($pos,0),$pos);
}
    

=head2 t_gamma - expand/decode encoded data to physical (gamma=1, max=1).

=for usage 

  $im_phys = $im->apply(t_gamma(gamma=>2.2));

=for ref

t_gamma is mostly used internally to handle gamma conversion of RGB
values before other transforms are applied.  It is automatically used
by the other transforms if you include a C<gamma> option to the
constructor.

=cut

sub t_gamma {
    my $opt = shift;

    my $me = _new($opt);

    $me->{name} = "Gamma decoding and scaling";
    $me->{func} = sub {
	my($d,$o) = @_;
	my $out;
	if($d->type !~ m/(double|float)/) {
	    $out = float $d;
	} else {
	    $out = $d->copy;
	}
	
	my $max = $o->{max};
	if(!defined($max)) {
	    # guess max based on type of input
	    if($d->type =~ m/(byte|short|ushort|long)/) {
		$max = 255;
	    } else {
		$max = 1.0;
	    }
	}

	$out /= $max;
	$out **= $o->{gamma} // 1;
	
	return $out;
    };
    $me->{inv} = sub {
	my($d,$o) = @_;
	my $out = double $d;
	die "t_gamma inverse: gamma encoding - can't encode with gamma=0!" unless($o->{gamma}//1);
	$out **= (1.0/( $o->{gamma} // 1) );
	$out *= $o->{max} // 1;
	return $out;
    };

    return $me;
}

=head2 t_cmyk - convert RGB to CMYK (or vice versa)

=for usage

 $cmyk = $im->apply(t_cmyk);

=for ref 

C<t_cmyk> converts to four-color separation subtractive process
values, maximizing black ink at the expense of the cyan, magenta, and
yellow channels.  Standard options (notably C<max> and C<gamma>) are
accepted, but the CMYK representation is always scaled 0-1, with a
gamma of unity.

Linear subtractive conversion is used -- thus the CMYK values
represent corrected halftone fraction with idealized subtractive
pigments that are exactly conjugate to the RGB colors.

Like most of the color conversions, C<t_cmyk> doesn't necessarily work
in the 0 dim of the input -- it attempts to find the color dim in one
of the first three dims of the input PDL.  That is because some image
manipulation code puts the colors in the 0 dim and some in the 2 dim.

=cut

sub t_cmyk {
    my $opt = shift;
    my $me = _new($opt);
    
    $me->{name} = "CMYK conversion";
    
    # Function and inverse work in 0-1 linear physical space (t_gamma composition fixes scaling)
    $me->{func} = sub {
	my($d,$o) = @_;

	my ($d2,$where) = _rectify($d,$o->{pos});  

	# Expand the color dim by one in the output (to make room for K in CMYK)
	my @dims = $d->dims;
	$dims[$where]++;

	# Generate the output to match the expanded input dims, and make a working
	# link into it with the active dim at 0
	my $out = PDL->new_from_specification($d2->type, @dims);
	my $oo = ($where) ? $out->mv($where,0) : $out;

	# Convert RGB->CMY and copy any extra information over
	$oo->(0:2) .= 1 - $d2->(0:2);
	if($oo->dim(0) > 4) {
	    $oo->(4:-1) .= $d2->(2:-1);
	}

	# Find the K channel
	$oo->((3)) .= $oo->(0:2)->minimum;
	$oo->(0:2) -= $oo->(3);

	# The $oo stuff flowed to $out; return that to preserve shape
	return $out;	  
    };

    $me->{inv} = sub {
	my($d,$o) = @_;

	my($d2, $where) = _rectify($d,$o->{pos});
	
	my @dims = $d->dims;
	if($dims[$where] < 4)  {
	    die "t_cmyk inverse: color dim has size ".$dims[0].", too small for cmyk (4 needed)\n";
	}
	$dims[$where]--;
	my $out = PDL->new_from_specification($d2->type, @dims);
	my $oo = ($where) ? $out->mv($where,0) : $out;
	
	# Convert CMY to RGB
	$oo->(0:2) .= 1 - $d2->(0:2);
	
	# Correct RGB downward for the K portion.
	$oo->(0:2) -= $d2->(3);
		
	return $out;
    };

    $me->{otype} = ['cyan','magenta','yellow','black'];
    $me->{ounit} = ['ink fraction','ink fraction','ink fraction','ink fraction'];
    return t_compose(t_gamma($me->{options}), $me);

}

##############################
# _t_hs_ivl - handle hsi, hsv, or hsl -- which differ only in their treatment
# of brightness.  Hue is normalized 0..1
#
# Definitions from en.wikipedia.org...
sub _hue_from_rgb {
    my $rgb = shift;
    my $m = shift;
    my $C = shift() - $m;

    # Find index of max...
    my $maxdex = $rgb->qsorti->(-1); # index (R,G,B) of maximum component
    my $dexes = ($maxdex + pdl(0,1,2)) % 3;
    my $wonky = $rgb->(:,*3)->index($dexes); # R,G,B permuted 
    my $offset = $maxdex->((0)) * 2;
    my $H = ( ($wonky->(1:2) * pdl(1,-1))->sumover ) / ($C + ($C==0)) + $offset ;
    $H += 6 * ($H<0);
    $H /= 6;
    return $H;
}

=head2 t_hsv - convert RGB to HSV (or vice versa)

=for usage

 $hsv = $im->apply(t_hsv);

=for ref

HSV is Hue/Saturation/Value.

=cut
sub t_hsv {
    my $opt = shift;
    
    my $me = _new($opt);
    $me->{name} = "HSV conversion";
    
    $me->{func} = sub {
	my ($d,$o) = @_;
	my ($d2, $where) = _rectify($d, $o->{pos});
	
	# Generate output to match the input dims, and make a working link
	# to it with the active dim at 0
	my $out = PDL->new_from_specification( $d2->type, $d->dims );
	my $oo = ($where) ? $out->mv($where,0) : $out;
	
	my $m = $d2->minimum;
	my $M = $d2->maximum;
	$oo->((0)) .= _hue_from_rgb($d2, $m, $M);  # H

	# Find min/max and chroma
	$oo->((2)) .= $M;                          # V
	
	$oo->((1)) .= ($M-$m) / ($M + ($M==0));    # S

	# Copy ancillary info...
	if($d2->dim(0) > 3) {
	    $oo->(3:-1) .= $d2->(3:-1);
	}
	return $out;
    };

    $me->{inv} = sub {
	my ($d,$o) = @_;
	my ($d2, $where) = _rectify($d, $o->{pos});
	
	# Generate output to match the input dims, and make a working link
	# to it with the active dim at 0
	my $out = PDL->new_from_specification( $d2->type, $d->dims );
	my $oo = ($where) ? $out->mv($where,0) : $out;

	my $Hp = $d2->((0)) * 6;
	$Hp -= 6 * ($Hp/6)->floor;

	my $wonky = PDL->new_from_specification($d2->type, $d2->dims );
	$wonky->((2)) .= 0;
	$wonky->((0)) .= $d2->((1)) * $d2->((2));      # C = S * V
	$wonky->((1)) .= $wonky->((0)) * (1 - ($Hp - 2 * (($Hp/2)->floor) - 1)->abs ); # X

	my $lookup = pdl([0,1,2],[1,0,2],[2,0,1],[2,1,0],[1,2,0],[0,2,1]);
	my $dex = ($Hp->floor);
	my $l2 = $lookup->mv(0,-1)->slice(":". (",*1"x$dex->ndims))->index($dex)->mv(-1,0);

	$oo->(0:2) .= $wonky->(:,*1)->index(  $l2  );  # R'G'B' 
	$oo += $d2->(2) - $wonky->(0);                 # Add in unsaturated part

	# Copy ancillary info...
	if($d2->dim(0) > 3) {
	    $oo->(3:-1) .= $d2->(3:-1);
	}
	return $out;
    };

    $me->{otype} = ['Hue','Saturation','Value'];
    $me->{ounit} = ['degrees','scaled','brightness'];
    return t_compose(t_gamma($me->{options}),$me);
}
	

=head2 t_hsv - convert RGB to HSV (or vice versa)

=head2 t_hsl - convert RGB to HSL (or vice versa)

=head2 t_yuv - convert RGB to YUV (or vice versa)




=head2 t_cielab - convert RGB to CIELAB 1976 (or vice versa)

=head2 t_ciexyz - convert RGB to CIE 1931 XYZ (or vice versa)
1;