File Coverage

blib/lib/ICC/Support/nMIX.pm

Criterion	Covered	Total	%
statement	17	344	4.9
branch	1	142	0.7
condition	0	72	0.0
subroutine	5	29	17.2
pod	1	10	10.0
total	24	597	4.0

package ICC::Support::nMIX;

use strict;
use Carp;

our $VERSION = 0.31;

# revised 2016-05-17
#
# Copyright © 2004-2020 by William B. Birkett

# inherit from Shared
use parent qw(ICC::Shared);

# enable static variables
use feature 'state';

# create new nMIX object
# parameter may be a hash -or- an ICC::Support::Chart object
# columns is an optional column slice for the chart clut data
# hash may contain references to clut array or delta array
# hash keys are: ('array', 'delta')
# parameter: ()
# parameter: (ref_to_attribute_hash)
# parameter: (chart_object, [columns])
# returns: (ref_to_object)
sub new {

	# get object class
	my $class = shift();

	# create empty nMIX object
	my $self = [
		{},    # object header
		[],    # clut
		[],    # delta
		[],    # clut exp
		undef, # clut exp cache (for Lapack)
	];

	# if there are parameters
	if (@_) {

		# if one parameter, a hash reference
		if (@_ == 1 && ref($_[0]) eq 'HASH') {

			# make new nMIX object from attribute hash
			_new_from_hash($self, @_);

		# if one or two parameters, first parameter an ICC::Support::Chart object
		} elsif ((@_ == 1 || @_ == 2) && UNIVERSAL::isa($_[0], 'ICC::Support::Chart')) {

			# make new nMIX object from chart
			_new_from_chart($self, @_);

		} else {

			# error
			croak('parameter must be a hash reference or ICC::Support::Chart object');

		}

	}

	# bless object
	bless($self, $class);

	# return object reference
	return($self);

}

# get/set reference to header hash
# parameters: ([ref_to_new_hash])
# returns: (ref_to_hash)
sub header {

	# get object reference
	my $self = shift();

	# if there are parameters
	if (@_) {

		# if one parameter, a hash reference
		if (@_ == 1 && ref($_[0]) eq 'HASH') {

			# set header to copy of hash
			$self->[0] = {%{$_[0]}};

		} else {

			# error
			croak('parameter must be a hash reference');

		}

	}

	# return reference
	return($self->[0]);

}

# get/set reference to clut array
# parameters: ([ref_to_new_array])
# returns: (ref_to_array)
sub array {

	# get object reference
	my $self = shift();

	# if there are parameters
	if (@_) {

		# if one parameter, a 2-D array reference
		if (@_ == 1 && ref($_[0]) eq 'ARRAY' && @{$_[0]} == grep {ref() eq 'ARRAY'} @{$_[0]}) {

			# copy array to object and bless
			$self->[1] = bless(Storable::dclone($_[0]), 'Math::Matrix');

		# if one parameter, a Math::Matrix object
		} elsif (@_ == 1 && UNIVERSAL::isa($_[0], 'Math::Matrix')) {

			# copy matrix to object
			$self->[1] = Storable::dclone($_[0]);

		} else {

			# error
			croak('clut must be a 2-D array reference or Math::Matrix object');

		}

		# for each corner point
		for my $i (0 .. $#{$self->[1]}) {

			# for each spectral value
			for my $j (0 .. $#{$self->[1][$i]}) {

				# set value to zero, if negative
				$self->[1][$i][$j] = 0 if ($self->[1][$i][$j] < 0);

			}

		}

		# update arrays
		_update_clut_exp($self);

	}

	# return reference
	return($self->[1]);

}

# get/set reference to delta array
# scalar value parameter fills array with that value
# parameters: ([ref_to_new_array -or- scalar_value])
# returns: (ref_to_array)
sub delta {

	# get object reference
	my $self = shift();

	# if there are parameters
	if (@_) {

		# if one parameter, a 1-D array reference
		if (@_ == 1 && ref($_[0]) eq 'ARRAY' && @{$_[0]} == grep {Scalar::Util::looks_like_number($_)} @{$_[0]}) {

			# make delta array
			$self->[2] = [@{$_[0]}];

		# if one parameter, a scalar number
		} elsif (@_ == 1 && Scalar::Util::looks_like_number($_[0])) {

			# if array is defined
			if (defined($self->[1])) {

				# make delta array
				$self->[2] = [($_[0]) x @{$self->[1][0]}];

			} else {

				# error
				croak ('array must be defined when specifying delta as a scalar');

			}

		} else {

			# error
			croak('\'delta\' parameter must be a scalar or a 1-D array reference');

		}

		# update arrays
		_update_clut_exp($self);

	}

	# return reference
	return($self->[2]);

}

# get number of input channels
# returns: (number)
sub cin {

	# get object reference
	my $self = shift();

	# return
	return(int(log(@{$self->[1]})/log(2) + 1E-12));

}

# get number of output channels
# returns: (number)
sub cout {

	# get object reference
	my $self = shift();

	# return
	return(scalar(@{$self->[1][0]}));

}

# transform data
# supported input types:
# parameters: (list, [hash])
# parameters: (vector, [hash])
# parameters: (matrix, [hash])
# parameters: (Math::Matrix_object, [hash])
# parameters: (structure, [hash])
# returns: (same_type_as_input)
sub transform {

	# set hash value (0 or 1)
	my $h = ref($_[-1]) eq 'HASH' ? 1 : 0;

	# if input a 'Math::Matrix' object
	if (@_ == $h + 2 && UNIVERSAL::isa($_[1], 'Math::Matrix')) {

		# call matrix transform
		&_trans2;

	# if input an array reference
	} elsif (@_ == $h + 2 && ref($_[1]) eq 'ARRAY') {

		# if array contains numbers (vector)
		if (! ref($_[1][0]) && @{$_[1]} == grep {Scalar::Util::looks_like_number($_)} @{$_[1]}) {

			# call vector transform
			&_trans1;

		# if array contains vectors (2-D array)
		} elsif (ref($_[1][0]) eq 'ARRAY' && @{$_[1]} == grep {ref($_) eq 'ARRAY' && Scalar::Util::looks_like_number($_->[0])} @{$_[1]}) {

			# call matrix transform
			&_trans2;

		} else {

			# call structure transform
			&_trans3;

		}

	# if input a list (of numbers)
	} elsif (@_ == $h + 1 + grep {Scalar::Util::looks_like_number($_)} @_) {

		# call list transform
		&_trans0;

	} else {

		# error
		croak('invalid transform input');

	}

}

# compute Jacobian matrix
# nominal input range is (0 - 1)
# hash key 'ubox' enables unit box extrapolation
# clipped outputs are extrapolated using Jacobian
# parameters: (input_vector, [hash])
# returns: (Jacobian_matrix, [output_vector])
sub jacobian {

	# get parameters
	my ($self, $in, $hash) = @_;

	# local variables
	my ($ext, $out, $jac, $jac_bc, $d, $s, $dx);

	# check if ICC::Support::Lapack module is loaded
	state $lapack = defined($INC{'ICC/Support/Lapack.pm'});

	# if extrapolation required (any input outside the unit box)
	if ($hash->{'ubox'} && grep {$_ < 0.0 || $_ > 1.0} @{$in}) {

		# compute intersection with unit box
		($ext, $in) = _intersect($in);

	}

	# if ICC::Support::Lapack module is loaded
	if ($lapack) {

		# compute output using Lapack module
		$out = ICC::Support::Lapack::nMIX_vec_trans($in, $self->[4], $self->[2]);

		# compute Jacobian using Lapack module
		$jac = ICC::Support::Lapack::nMIX_jacobian($in, $self->[4], $self->[2], $out);

		# bless Jacobian as Math::Matrix object
		bless($jac, 'Math::Matrix');

	} else {

		# compute output values
		$out = _trans1($self, $in);

		# compute the barycentric jacobian
		$jac_bc = _barycentric_jacobian($in);

		# compute Jacobian (before exponentiation)
		$jac = $self->[3]->transpose * $jac_bc;

		# for each row (output)
		for my $i (0 .. $#{$jac}) {

			# get delta value
			$d = $self->[2][$i];

			# skip if delta is one
			next if ($d == 1);

			# get output value
			$s = $out->[$i];

			# compute exponentiation adjustment
			$dx = $s ? $d ? abs($s)**(1 - $d)/$d : $s : 0;

			# for each column (input)
			for my $j (0 .. $#{$jac->[0]}) {

				# adjust for exponentiation
				$jac->[$i][$j] *= $dx;

			}

		}

	}

	# if output values wanted
	if (wantarray) {

		# if output extrapolated
		if (defined($ext)) {

			# for each output
			for my $i (0 .. $#{$out}) {

				# add delta value
				$out->[$i] += ICC::Shared::dotProduct($jac->[$i], $ext);

			}

		}

		# return Jacobian and output vector
		return($jac, $out);

	} else {

		# return Jacobian only
		return($jac);

	}

}

# compute parametric Jacobian matrix
# parameters: (input_vector)
# returns: (parametric_jacobian_matrix)
sub parajac {

	# get parameters
	my ($self, $in) = @_;

	# return Jacobian matrix
	return(Math::Matrix->diagonal(_parametric($self, $in)));

}

# print object contents to string
# format is an array structure
# parameter: ([format])
# returns: (string)
sub sdump {

	# get parameters
	my ($self, $p) = @_;

	# local variables
	my ($s, $fmt);

	# resolve parameter to an array reference
	$p = defined($p) ? ref($p) eq 'ARRAY' ? $p : [$p] : [];

	# get format string
	$fmt = defined($p->[0]) && ! ref($p->[0]) ? $p->[0] : 'undef';

	# set string to object ID
	$s = sprintf("'%s' object, (0x%x)\n", ref($self), $self);

	# return
	return($s);

}

# transform list
# parameters: (object_reference, list, [hash])
# returns: (list)
sub _trans0 {

	# local variables
	my ($self, $hash);
	my ($coef, $sum, $product, @out);

	# check if ICC::Support::Lapack module is loaded
	state $lapack = defined($INC{'ICC/Support/Lapack.pm'});

	# get object reference
	$self = shift();

	# get optional hash
	$hash = pop() if (ref($_[-1]) eq 'HASH');

	# if ICC::Support::Lapack module is loaded
	if ($lapack) {

		# compute output using Lapack module
		@out = @{ICC::Support::Lapack::nMIX_vec_trans(\@_, $self->[4], $self->[2])};

	} else {

		# compute output using '_trans1'
		@out = @{_trans1($self, \@_)};

	}

	# return
	return(@out);

}

# transform vector
# parameters: (object_reference, vector, [hash])
# returns: (vector)
sub _trans1 {

	# get parameters
	my ($self, $in, $hash) = @_;

	# local variables
	my ($coef, $sum, $out);

	# check if ICC::Support::Lapack module is loaded
	state $lapack = defined($INC{'ICC/Support/Lapack.pm'});

	# if ICC::Support::Lapack module is loaded
	if ($lapack) {

		# compute output using Lapack module
		$out = ICC::Support::Lapack::nMIX_vec_trans($in, $self->[4], $self->[2]);

	} else {

		# compute barycentric coefficients
		$coef = _barycentric($in);

		# for each output value
		for my $i (0 .. $#{$self->[3][0]}) {

			# initialize sum
			$sum = 0;

			# for each coefficient
			for my $j (0 .. $#{$coef}) {

				# add product to sum
				$sum += $self->[3][$j][$i] * $coef->[$j] if ($coef->[$j]);

			}

			# save result
			$out->[$i] = _pow1p($sum, $self->[2][$i]);

		}

	}

	# return
	return($out);

}

# transform matrix (2-D array -or- Math::Matrix object)
# parameters: (object_reference, matrix, [hash])
# returns: (matrix)
sub _trans2 {

	# get parameters
	my ($self, $in, $hash) = @_;

	# local variables
	my ($coef, $sum, $product, $mean, $ratio, $out);

	# check if ICC::Support::Lapack module is loaded
	state $lapack = defined($INC{'ICC/Support/Lapack.pm'});

	# if ICC::Support::Lapack module is loaded
	if ($lapack) {

		# compute output using Lapack module
		$out = ICC::Support::Lapack::nMIX_mat_trans($in, $self->[4], $self->[2]);

	} else {

		# for each input vector
		for my $k (0 .. $#{$in}) {

			# compute barycentric coefficients
			$coef = _barycentric($in->[$k]);

			# for each output value
			for my $i (0 .. $#{$self->[3][0]}) {

				# initialize sum
				$sum = 0;

				# for each coefficient
				for my $j (0 .. $#{$coef}) {

					# add product to sum
					$sum += $self->[3][$j][$i] * $coef->[$j] if ($coef->[$j]);

				}

				# save result
				$out->[$k][$i] = _pow1p($sum, $self->[2][$i]);

			}

		}

	}

	# return
	return(UNIVERSAL::isa($in, 'Math::Matrix') ? bless($out, 'Math::Matrix') : $out);

}

# transform structure
# parameters: (object_reference, structure, [hash])
# returns: (structure)
sub _trans3 {

	# get parameters
	my ($self, $in, $hash) = @_;

	# local variables
	my ($out);

	# transform the array structure
	_crawl($self, $in, $out = [], $hash);

	# return
	return($out);

}

# recursive transform
# array structure is traversed until scalar arrays are found and transformed
# parameters: (ref_to_object, input_array_reference, output_array_reference, hash)
sub _crawl {

	# get parameters
	my ($self, $in, $out, $hash) = @_;

	# if input is a vector (reference to a scalar array)
	if (@{$in} == grep {! ref()} @{$in}) {

		# transform input vector and copy to output
		@{$out} = @{_trans1($self, $in, $hash)};

	} else {

		# for each input element
		for my $i (0 .. $#{$in}) {

			# if an array reference
			if (ref($in->[$i]) eq 'ARRAY') {

				# transform next level
				_crawl($self, $in->[$i], $out->[$i] = [], $hash);

			} else {

				# error
				croak('invalid transform input');

			}

		}

	}

}

# update exponentiated CLUT
# parameter: (ref_to_object)
sub _update_clut_exp {

	# get object reference
	my $self = shift();

	# check if ICC::Support::Lapack module is loaded
	state $lapack = defined($INC{'ICC/Support/Lapack.pm'});

	# if clut and delta are defined and not null
	if (defined($self->[1]) && @{$self->[1]} && defined($self->[2]) && @{$self->[2]}) {

		# if ICC::Support::Lapack module is loaded
		if ($lapack) {

			# compute mixture array
			$self->[3] = ICC::Support::Lapack::nMIX_power($self->[1], $self->[2]);

			# create cached mixture array
			$self->[4] = ICC::Support::Lapack::cache_2D($self->[3]);

		} else {

			# for each row
			for my $i (0 .. $#{$self->[1]}) {

				# for each column
				for my $j (0 .. $#{$self->[1][0]}) {

					# exponentiate corner point value
					$self->[3][$i][$j] = _powm1($self->[1][$i][$j], $self->[2][$j]);

				}

			}

		}

	}

	# bless as Math::Matrix object
	bless($self->[3], 'Math::Matrix');

}

# compute barycentric coefficients
# parameter: (input_vector)
# returns: (coefficient_array)
sub _barycentric {

	# get parameter
	my $dev = shift();

	# local variables
	my ($devc, $coef);

	# compute complement values
	$devc = [map {1 - $_} @{$dev}];

	# initialize coefficient array
	$coef = [(1.0) x 2**@{$dev}];

	# for each coefficient
	for my $i (0 .. $#{$coef}) {

		# for each device value
		for my $j (0 .. $#{$dev}) {

			# if $j-th bit set
			if ($i >> $j & 1) {

				# multiply by device value
				$coef->[$i] *= $dev->[$j];

			} else {

				# multiply by (1 - device value)
				$coef->[$i] *= $devc->[$j];

			}

		}

	}

	# return
	return($coef);

}

# compute barycentric Jacobian matrix
# parameter: (input_vector)
# returns: (Jacobian_matrix)
sub _barycentric_jacobian {

	# get parameter
	my $dev = shift();

	# local variables
	my ($devc, $rows, $jac);

	# compute complement values
	$devc = [map {1 - $_} @{$dev}];

	# compute matrix rows
	$rows = 2**@{$dev};

	# for each matrix row
	for my $i (0 .. $rows - 1) {

		# initialize row
		$jac->[$i] = [(1.0) x @{$dev}];

		# for each matrix column
		for my $j (0 .. $#{$dev}) {

			# for each device value
			for my $k (0 .. $#{$dev}) {

				# if $k-th bit set
				if ($i >> $k & 1) {

					# multiply by device value -or- 1 (skip)
					$jac->[$i][$j] *= $dev->[$k] if ($j != $k);

				} else {

					# multiply by (1 - device value) -or- -1
					$jac->[$i][$j] *= ($j != $k) ? $devc->[$k] : -1;

				}

			}

		}

	}

	# return
	return(bless($jac, 'Math::Matrix'));

}

# find unit box intersection
# with line from input to box-center
# parameters: (input_vector)
# returns: (extrapolation_vector, intersection_vector)
sub _intersect {

	# get input values
	my ($in) = shift();

	# local variables
	my (@cin, $dmax, $ubox, $ext);

	# compute input to box-center difference
	@cin = map {$_ - 0.5} @{$in};

	# initialize
	$dmax = 0;

	# for each difference
	for (@cin) {

		# if larger absolute value
		if (abs($_) > $dmax) {

			# new max difference
			$dmax = abs($_);

		}

	}

	# multiply max difference by 2
	$dmax *= 2;

	# compute intersection vector (on surface of unit box)
	$ubox = [map {$_/$dmax + 0.5} @cin];

	# compute extrapolation vector (as Math::Matrix object)
	$ext = [map {$in->[$_] - $ubox->[$_]} (0 .. $#{$in})];

	# return
	return($ext, $ubox);

}

# compute parametric partial derivatives
# parameters: (object_reference, input_vector)
# returns: (partial_derivative_vector)
sub _parametric {

	# get parameters
	my ($self, $in) = @_;

	# local variables
	my ($bc, $cp, $d, $s1, $sum1, $sum2, $pj, $dk, $pd, $r);

	# calculate barycentric coordinates
	$bc = _barycentric($in);

	# get corner point matrix
	$cp = $self->[1];

	# for each delta/output value
	for my $i (0 .. $#{$self->[2]}) {

		# get delta value
		$d = $self->[2][$i];

		# if delta non-zeroish
		if (abs($d) >= 1E-5) {

			# initialize sums
			$sum1 = $sum2 = 0;

			# for each corner point
			for my $j (0 .. $#{$cp}) {

				# accumulate sums
				$sum1 += $s1 = $bc->[$j] * $cp->[$j][$i]**$d;
				$sum2 += $s1 * log($cp->[$j][$i]);

			}

			# compute partial derivative
			$pj->[$i] = $sum1**(1/$d) * ($sum2/$sum1 - log($sum1)/$d)/$d;

		} else {

			# for delta +/- 1E-5
			for my $k (0 .. 1) {

				# set delta
				$dk = $k ? -1E-5 : 1E-5;

				# initialize sums
				$sum1 = $sum2 = 0;

				# for each corner point
				for my $j (0 .. $#{$cp}) {

					# accumulate sums
					$sum1 += $s1 = $bc->[$j] * $cp->[$j][$i]**$dk;
					$sum2 += $s1 * log($cp->[$j][$i]);

				}

				# compute partial derivative
				$pd->[$k] = $sum1**(1/$dk) * ($sum2/$sum1 - log($sum1)/$dk)/$dk;

			}

			# compute interpolation ratio
			$r = 0.5 - $d/2E-5;

			# interpolate partial derivative
			$pj->[$i] = $r * $pd->[0] + (1 - $r) * $pd->[1];

		}

	}

	# return array of partial derivatives
	return($pj);

}

# exponentiate corner point
# uses expm1 function for small exponents
# parameters: (base, exponent)
# returns: (base^exponent)
sub _powm1 {

	# get parameters
	my ($base, $exp) = @_;

	if ($exp == 0.0) {

		if ($base > 0.0) {

			return(log($base))

		} else {

			return(-(DBL_MAX))

		}

	} elsif ($exp < 1.0) {

		if ($base > 0.0) {

			return(ICC::Support::Lapack::expm1(log($base) * $exp))

		} else {

			return(-1.0)

		}

	} else {

		if ($base > 0.0) {

			return(POSIX::pow($base, $exp))

		} else {

			return(0)

		}

	}

}

# exponentiate barycentric sum
# uses log1p function for small exponents
# parameters: (base, exponent)
# returns: (base^(1/exponent))
sub _pow1p {

	# get parameters
	my ($base, $exp) = @_;

	if ($exp == 0.0) {

		return(exp($base))

	} elsif ($exp < 1.0) {

		if ($base > -1.0) {

			return(exp(ICC::Support::Lapack::log1p($base)/$exp))

		} else {

			return(0.0)

		}

	} else {

		if ($base > 0.0) {

			return(POSIX::pow($base, 1.0/$exp))

		} else {

			return(0.0)

		}

	}

}

# make new nMIX object from attribute hash
# hash may contain pointers to clut array or delta array
# hash keys are: ('array', 'delta')
# object elements not specified in the hash are unchanged
# parameters: (ref_to_object, ref_to_attribute_hash)
sub _new_from_hash {

	# get parameters
	my ($self, $hash) = @_;

	# local variable
	my ($value);

	# if 'array' attribute
	if (defined($value = $hash->{'array'})) {

		# if reference to a 2-D array
		if (ref($value) eq 'ARRAY' && @{$value} == grep {ref() eq 'ARRAY'} @{$value}) {

			# copy array to object and bless
			$self->[1] = bless(Storable::dclone($value), 'Math::Matrix');

		# if reference to a Math::Matrix object
		} elsif (UNIVERSAL::isa($value, 'Math::Matrix')) {

			# copy matrix to object
			$self->[1] = Storable::dclone($value);

		} else {

			# wrong data type
			croak('\'array\' must be a 2-D array reference or Math::Matrix object');

		}

		# for each corner point
		for my $i (0 .. $#{$self->[1]}) {

			# for each spectral value
			for my $j (0 .. $#{$self->[1][$i]}) {

				# set value to zero, if negative
				$self->[1][$i][$j] = 0 if ($self->[1][$i][$j] < 0);

			}

		}

	}

	# if 'delta' attribute
	if (defined($value = $hash->{'delta'})) {

		# if reference to a 1-D array (vector)
		if (ref($value) eq 'ARRAY' && @{$value} == grep {Scalar::Util::looks_like_number($_)} @{$value}) {

			# make delta array
			$self->[2] = [@{$value}];

		# if scalar number
		} elsif (Scalar::Util::looks_like_number($value)) {

			# if array is defined
			if (defined($self->[1])) {

				# make delta array
				$self->[2] = [($value) x @{$self->[1][0]}];

			} else {

				# error
				croak ('array must be defined when specifying delta as a scalar');

			}

		} else {

			# wrong data type
			croak('\'delta\' must be a scalar or a 1-D array reference');

		}

	}

	# update arrays
	_update_clut_exp($self);

}

# make new nMIX object from ICC::Support::Chart object
# chart must contain device values and, spectral or XYZ values
# copies corner points into clut array, warns if corner point(s) missing
# parameters: (ref_to_object, ref_to_chart, [columns])
sub _new_from_chart {

	# get parameters
	my ($self, $chart, $cols) = @_;

	# local variables
	my ($dev, $fmt, $devc, $cs);

	# verify chart has device values
	($dev = $chart->device()) or croak('chart must have device values');

	# verify clut column slice is defined
	(defined($cols) || ($cols = $chart->spectral() || $chart->xyz())) or croak('clut data slice is undefined');

	# get format keys
	$fmt = $chart->fmt_keys($cols);

	# verify spectral or XYZ data
	((@{$fmt} == grep {m/^(?:(.*)\|)?(?:nm|SPECTRAL_NM_|SPECTRAL_NM|SPECTRAL_|NM_|R_)\d{3}$/} @{$fmt}) || (3 == grep {m/^(?:(.*)\|)?XYZ_[XYZ]$/} @{$fmt})) || warn('clut data neither spectral nor XYZ');

	# for each corner point
	for my $i (0 .. 2**@{$dev} - 1) {

		# for each device channel
		for my $j (0 .. $#{$dev}) {

			# get device value
			$devc->[$j] = $i >> $j & 1;

		}

		# get corner point samples
		($cs = $chart->ramp(sub {@{$devc} == grep {$devc->[$_] == $_[$_]} (0 .. $#{$devc})})) or croak("no chart samples for corner point [@{$devc}]\n");

		# save clut vector
		$self->[1][$i] = $chart->slice($chart->add_avg($cs), $cols)->[0];

		# discard avg sample
		pop(@{$chart->array()});

		# for each spectral value
		for my $j (0 .. $#{$self->[1][$i]}) {

			# if value is negative
			if ($self->[1][$i][$j] < 0) {

				# set value to zero
				$self->[1][$i][$j] = 0;

				# print warning
				print "clut value [$i][$j] was negative, set to 0\n";

			}

		}

	}

	# bless clut
	bless($self->[1], 'Math::Matrix');

}

1;