#===============================================
# Color
#===============================================
package Color;
use Exporter;
@ISA = qw(Exporter);

#公開したいサブルーチン
#@EXPORT = qw(aa );

use strict;

use Utils;
use JFile;
use CSV;
use GraphData;
use Sci::GeneralFileFormat;
use Sci::Science;

my $DataPath =  "D:\\Programs\\Perl\\Color\\data";
my $ColormapGIF               = Deps::MakePath($DataPath, "CIE-ColorMap.jpeg", 0);
my $ColorMatchingFunctionPath = Deps::MakePath($DataPath, "ColorFunction.txt", 0);
my $StandardLightPath         = Deps::MakePath($DataPath, "CIEStandardLight.txt", 0);
my $DayLightPath              = Deps::MakePath($DataPath, "DayLightS.txt", 0);
my $HalogenLampPath           = Deps::MakePath($DataPath, "HalogenLamp.lsc", 0);
my $SolarSpectraPath          = Deps::MakePath($DataPath, "SolarSpectra.txt", 0);

sub Spectrum                  { return shift->{Spectrum}; }
sub ColorMatchingFunctionType { return shift->{ColorMatchingFunctionType}; }
sub ColorMatchingFunction     { return shift->{ColorMatchingFunction}; }
sub LightSource               { return shift->{LightSource}; }

my @ImageCoordinatetoXYZMapping =
	( ( 34  , 367) => (0,   0),
	  (368.5, 367) => (0.8, 0),
	  ( 33  ,  30) => (0, 0.8)
	);

my $h  = Sci::h(); #6.62E-34;
my $c  = Sci::c(); #3.00E+08;
my $e  = Sci::e(); #1.602E-19;
my $pi = Sci::pi(); #3.141593;

sub ColorMapImagePath
{
	return $ColormapGIF;
}

sub GetDefinedxy
{
#$ColorMap: NTSC, HDTV, PrimaryColor
	my ($ColorMap) = @_;
	if($ColorMap =~ /ntsc/i) {
#R, G, B, W
		return ( (0.67, 0.3), (0.21, 0.71), (0.14, 0.08), (0.310, 0.316) );
	}
	elsif($ColorMap =~ /hdtv/i) {
		return ( (0.640, 0.330), (0.300, 0.600), (0.150, 0.060), (0.3127, 0.3290) );
	}
	elsif($ColorMap =~ /PrimaryColor/i) {
		return ( (0.735, 0.265), (0.274, 0.717), (0.167, 0.0090), (-1, -1) );
	}
	return undef;
}

sub XYZtoxyz
{
	my ($X, $Y, $Z) = @_;
	my $t = $X + $Y + $Z;
	return ($X/$t, $Y/$t, $Z/$t);
}

sub RGBtorgb
{
	my ($R, $G, $B) = @_;
	my $t = $R + $G + $B;
	return ($R/$t, $G/$t, $B/$t);
}

sub XYZtoRGB
{
	my ($X, $Y, $Z) = @_;
	my $R =  0.41844 * $X - 0.15866 * $Y - 0.08283 * $Z;
	my $G = -0.09117 * $X + 0.25242 * $Y + 0.01570 * $Z;
	my $B =  0.00092 * $X - 0.00255 * $Y + 0.17858 * $Z;
	return ($R, $G, $B);
}

sub ImageCoordinatetoXYZMapping
{
	my ($x, $y) = @_;
	my $x0  = $ImageCoordinatetoXYZMapping[0];
	my $y0  = $ImageCoordinatetoXYZMapping[1];
	my $vx0 = $ImageCoordinatetoXYZMapping[2];
	my $vy0 = $ImageCoordinatetoXYZMapping[3];
	my $x1  = $ImageCoordinatetoXYZMapping[4];
	my $y1  = $ImageCoordinatetoXYZMapping[5];
	my $vx1 = $ImageCoordinatetoXYZMapping[6];
	my $vy1 = $ImageCoordinatetoXYZMapping[7];
	my $x2  = $ImageCoordinatetoXYZMapping[8];
	my $y2  = $ImageCoordinatetoXYZMapping[9];
	my $vx2 = $ImageCoordinatetoXYZMapping[10];
	my $vy2 = $ImageCoordinatetoXYZMapping[11];
	my $dx = $x1 - $x0;
	my $dy = $y2 - $y0;
	my $vx = $vx0 + ($vx1-$vx0) / $dx * ($x - $x0);
	my $vy = $vy0 + ($vy2-$vy0) / $dy * ($y - $y0);
	return ($vx, $vy);
}

sub xyToImageCoordinate
{
	my ($vx, $vy) = @_;
	my $x0  = $ImageCoordinatetoXYZMapping[0];
	my $y0  = $ImageCoordinatetoXYZMapping[1];
	my $vx0 = $ImageCoordinatetoXYZMapping[2];
	my $vy0 = $ImageCoordinatetoXYZMapping[3];
	my $x1  = $ImageCoordinatetoXYZMapping[4];
	my $y1  = $ImageCoordinatetoXYZMapping[5];
	my $vx1 = $ImageCoordinatetoXYZMapping[6];
	my $vy1 = $ImageCoordinatetoXYZMapping[7];
	my $x2  = $ImageCoordinatetoXYZMapping[8];
	my $y2  = $ImageCoordinatetoXYZMapping[9];
	my $vx2 = $ImageCoordinatetoXYZMapping[10];
	my $vy2 = $ImageCoordinatetoXYZMapping[11];
	my $dvx = $vx1 - $vx0;
	my $dvy = $vy2 - $vy0;
	my $x = $x0 + ($x1-$x0) / $dvx * ($vx - $vx0);
	my $y = $y0 + ($y2-$y0) / $dvy * ($vy - $vy0);
	return ($x, $y);
}

sub new
{
	my ($module) = @_;
	my $this = {};
	bless $this;

	return $this;
}

sub DESTROY
{
	my $this = shift;
}

sub CalculateColor
{
	my ($this) = @_;
#print "a: ", $this->Spectrum(), "\n";
	my $Spectrum              = $this->Spectrum()->GetGraphData(0);
#print "a: ", $this->LightSource(), "\n";
	my $LightSource           = $this->LightSource()->GetGraphData(0);
#print "a: ", $this->ColorMatchingFunctionType(), "\n";
	my $ColorMatchingFunction = $this->ColorMatchingFunction()->GetGraphData(0);
#my $out = new JFile;
#$out->Open("d:/a.csv", "w");
	my ($X, $Y, $Z, $Y0) = (0,0,0,0);
	my $nData = $Spectrum->nData();
#print "nData: $nData\n";
	my ($wlmin, $wlmax) = $Spectrum->GetXMinMax();
	my $wlstep = 0.5;
	my $nl = int(($wlmax - $wlmin) / $wlstep + 1.001);
	my ($xmin, $xmax) = $ColorMatchingFunction->GetXMinMax();
	for(my $i = 0 ; $i < $nl-1 ; $i++) {
		my $x0  = $wlmin + $wlstep * $i;
		next if($x0 < $xmin or $xmax < $x0);
my $Debug = 0;
$Debug = 1 if($i % 10 == 0);
$Debug = 0;
		my $RT0 = $Spectrum->YVal(0, $x0, $Debug);
		my $l00 = $LightSource->YVal(0, $x0);
		my $S00 = $ColorMatchingFunction->YVal(0, $x0);
		my $S10 = $ColorMatchingFunction->YVal(1, $x0);
		my $S20 = $ColorMatchingFunction->YVal(2, $x0);
		my $RTl0 = $RT0 * $l00;
printf "$i: [%6.3f][%6.3f*%6.3f=%6.3f] (%4.3f,%4.3f,%4.3f)\n",
	$x0, $RT0, $l00, $RTl0, $S00, $S10, $S20 if($i % 10 == 0);
#$out->print("$i,$x0,$RT0,$l00,$RTl0,$S00,$S10,$S20\n");
		my $x1  = $wlmin + $wlstep * ($i+1);
		next if($x1 < $xmin or $xmax < $x1);
		my $RT1 = $Spectrum->YVal(0, $x1);
		my $l01 = $LightSource->YVal(0, $x1);
		my $S01 = $ColorMatchingFunction->YVal(0, $x1);
		my $S11 = $ColorMatchingFunction->YVal(1, $x1);
		my $S21 = $ColorMatchingFunction->YVal(2, $x1);
		my $RTl1 = $RT1 * $l01;
#print "    1: [$x1][$RT1][$l01][$S01][$S11][$S21]\n";

		my $dx = abs($x1 - $x0) / 2.0;
		$X  += ($RTl0*$S00 + $RTl1*$S01) * $dx;
		$Y  += ($RTl0*$S10 + $RTl1*$S11) * $dx;
		$Z  += ($RTl0*$S20 + $RTl1*$S21) * $dx;
		$Y0 += ($l00*$S10 + $l01*$S11) * $dx;
	}
#$out->Close();
	$Y0 = 100.0 / $Y0;
	return 	($X*$Y0, $Y*$Y0, $Z*$Y0);
}

sub ReadSpectrum
{
#$Observation: R or T
	my ($this, $path, $Observation) = @_;
	return undef if(!defined $path or $path eq '');
	my $Obs = substr($Observation, 0, 1);
	my $CSV = new CSV;
	if(!$CSV->Read($path)) {
		return undef;
	}
	my $LabelArray = $CSV->LabelArray;
	my $DataArray = $CSV->DataArray;
	my ($ObsIndex, $WLIndex);
	for(my $i = 0 ; $i < @$LabelArray ; $i++) {
print "$i: $LabelArray->[$i]\n";
		if($LabelArray->[$i] =~ /^$Obs/i) {
			$ObsIndex = $i;
		}
		if($LabelArray->[$i] =~ /^wl/i) {
			$WLIndex = $i;
		}
	}
print "index = ($WLIndex, $ObsIndex)\n";

	my $GraphArray = new GraphDataArray;
	my $Data0 = new GraphData;
	$Data0->SetTitle("Spectrum ($Observation)");

	$Data0->{'x0'}      = $DataArray->[$WLIndex];
	$Data0->{'x0_Name'} = "lambda / nm";
	$Data0->{'y0'}      = $DataArray->[$ObsIndex];
	$Data0->{'y0_Name'} = "Spectrum($Observation)";
	$Data0->CalMinMax();
	$GraphArray->AddGraphData($Data0);

	return $this->{Spectrum} = $GraphArray;
}

sub ReadLightSource
{
#$LightSource: Solar(All|Direct), HalogenLamp, DayLight, StdLight(A,D65,C,etc)
#              BlackBody
	my ($this, $LightSource, $T) = @_;

	my $FileName = '';
	if($LightSource =~ /Solar/i) {
		$FileName = $SolarSpectraPath;
	}
	elsif($LightSource =~ /HalogenLamp/i) {
		$FileName = $HalogenLampPath;
	}
	elsif($LightSource =~ /DayLight/i) {
		$FileName = $DayLightPath;
	}
	elsif($LightSource =~ /StdLight/i) {
		$FileName = $StandardLightPath;
	}

	my $GraphArray = new GraphDataArray;
	my $Data0 = new GraphData;
	$Data0->SetTitle($LightSource);
	my @X;
	my @Y1;
	my @Y2;
	if($LightSource =~ /BlackBody/i) {
		my $l0 = 200.0;
		my $l1 = 3100.0;
		my $lstep = 20.0;
		my $nl = int( ($l1 - $l0) / $lstep + 1.001);
		for(my $i = 0 ; $i < $nl ; $i++) {
			my $l = $l0 + $lstep * $i;
			my $nu = $c / $l / 0.000000001;
			my $hv = $h * $nu / $e;
			my $y = $nu**5 / (exp($hv/($T/300.0*0.026))-1) * 8 * $pi * $h / $c**4;
print "$l\t$nu\t$hv\t$y\n";
			$X[$i]  = $l;
			$Y1[$i] = $y;
		}

		$Data0->{'x0'}      = \@X;
		$Data0->{'x0_Name'} = "lambda / nm";
		$Data0->{'y0'}      = \@Y1;
		$Data0->{'y0_Name'} = "BlackBody($T K)";
	}
	elsif($LightSource =~ /Monochromatic/i) {
		my $l0 = 200.0;
		my $l1 = 800.0;
		my $lstep = 0.5;
		my $nl = int( ($l1 - $l0) / $lstep + 1.001);
		for(my $i = 0 ; $i < $nl ; $i++) {
			my $l = $l0 + $lstep * $i;
			$X[$i] = $l;
			$Y1[$i] = 0.0;
		}
		my $idx = ($T - $l0) / $lstep;
print "idx=$idx T=$T\n";
		for(my $i = $idx - 1 ; $i <= $idx + 1 ; $i++) {
			next if($idx < 0 or $idx >= $nl);
			$Y1[$i] = 1.0;
		}

		$Data0->{'x0'}      = \@X;
		$Data0->{'x0_Name'} = "lambda / nm";
		$Data0->{'y0'}      = \@Y1;
		$Data0->{'y0_Name'} = "Monochromatic($T nm)";
	}
	elsif($LightSource =~ /DayLight/i) {
		my $in = new JFile;
		if(!$in->Open($FileName, "r")) {
			return undef;
		}
		my $FirstLine = $in->ReadLine();
		my $SecondLine = $in->ReadLine();

		my $count = 0;
		my (@lambda, @S0, @S1, @S2);
		while(!$in->eof()) {
			my $line = $in->ReadLine();
			last if($in->eof());
			my @a = Utils::Split("\\s+", $line);
print "$a[0]\t$a[1]\t$a[2]\t$a[3]\n";
			$lambda[$count] = $a[0];
			$S0[$count]     = $a[1];
			$S1[$count]     = $a[2];
			$S2[$count]     = $a[3];
			$count++;
		}
		$in->Close();

		my $xD;
		if($T >= 4000.0 and $T <= 7000.0) {
			$xD = -4.6070e9 / $T**3 + 2.9678e6 / $T**2 + 0.09911e3 / $T + 0.244063;
		}
		else {
			$xD = -2.0064e9 / $T**3 + 1.9081e6 / $T**2 + 0.24748e3 / $T + 0.237040;
		}
		my $yD = -3.000*$xD**2 + 2.870 * $xD - 0.275;
		my $M1 = (-1.3515 - 1.7703 * $xD + 5.9114 * $yD) /
				(0.0241 + 0.2562 * $xD - 0.7341 * $yD);
		my $M2 = (0.0300 - 31.442 * $xD + 30.0717 * $yD) /
				(0.0241 + 0.2562 * $xD - 0.7341 * $yD);
print "(xD,yD) = ($xD, $yD)\n";
print "(M1,M2) = ($M1, $M2)\n";

		for(my $i = 0 ; $i < $count ; $i++) {
			$X[$i]  = $lambda[$i];
			$Y1[$i] = $S0[$i] + $M1*$S1[$i] + $M2*$S2[$i];
		}
		$Data0->{'x0'}      = \@X;
		$Data0->{'x0_Name'} = "lambda / nm";
		$Data0->{'y0'}      = \@Y1;
		$Data0->{'y0_Name'} = "DayLight($T K)";
		$Data0->{'x1'}      = \@X;
		$Data0->{'x1_Name'} = "lambda / nm";
		$Data0->{'y1'}      = \@S0;
		$Data0->{'y1_Name'} = "DayLight($T K)";
		$Data0->{'x2'}      = \@X;
		$Data0->{'x2_Name'} = "lambda / nm";
		$Data0->{'y2'}      = \@S1;
		$Data0->{'y2_Name'} = "DayLight($T K)";
		$Data0->{'x3'}      = \@X;
		$Data0->{'x3_Name'} = "lambda / nm";
		$Data0->{'y3'}      = \@S2;
		$Data0->{'y3_Name'} = "DayLight($T K)";
	}
	elsif($LightSource =~ /StdLight/i) {
		my $in = new JFile;
		if(!$in->Open($FileName, "r")) {
			return undef;
		}
		my $FirstLine = $in->ReadLine();
		my $SecondLine = $in->ReadLine();

		my $count = 0;
		while(!$in->eof()) {
			my $line = $in->ReadLine();
			last if($in->eof());
			my @a = Utils::Split("\\s+", $line);
			my $y;
			if($LightSource =~ /StdLight.*A/i) {
				$y = $a[1];
			}
			elsif($LightSource =~ /StdLight.*D65/i) {
				$y = $a[2];
			}
			elsif($LightSource =~ /StdLight.*C/i) {
				$y = $a[3];
			}
			elsif($LightSource =~ /StdLight.*D50/i) {
				$y = $a[4];
			}
			elsif($LightSource =~ /StdLight.*D55/i) {
				$y = $a[5];
			}
			elsif($LightSource =~ /StdLight.*D75/i) {
				$y = $a[6];
			}
			elsif($LightSource =~ /StdLight.*B/i) {
				$y = $a[7];
			}
print "$a[0]\t$y\n";
			next if($y eq '-');
			$X[$count]  = $a[0];
			$Y1[$count] = $y;
			$count++;
		}
		$in->Close();

		$Data0->{'x0'}      = \@X;
		$Data0->{'x0_Name'} = "lambda / nm";
		$Data0->{'y0'}      = \@Y1;
		$Data0->{'y0_Name'} = "Solar(Total)";
	}
	else {
		my $in = new JFile;
		if(!$in->Open($FileName, "r")) {
			return undef;
		}
		my $FirstLine = $in->ReadLine();

		my $count = 0;
		while(!$in->eof()) {
			my $line = $in->ReadLine();
			last if($in->eof());
			my @a = Utils::Split("\\s+", $line);
			if($LightSource =~ /HalogenLamp/i) {
				$a[0] *= 1000.0;
				$a[1] = $a[2] / 1000.0;
			}
print "$a[0]\t$a[1]\t$a[2]\n";
			$X[$count] = $a[0];
			$Y1[$count]  = $a[1];
			$Y2[$count] = $a[2];
			$count++;
		}
		$in->Close();

		$Data0->{'x0'}      = \@X;
		$Data0->{'x0_Name'} = "lambda / nm";
		if($LightSource =~ /HalogenLamp/i) {
			$Data0->{'y0'}      = \@Y1;
			$Data0->{'y0_Name'} = "HalogenLamp";
		}
		elsif($LightSource =~ /Solar.*(All)/) {
			$Data0->{'y0'}      = \@Y1;
			$Data0->{'y0_Name'} = "Solar(Total)";
		}
		elsif($LightSource =~ /Solar.*(Direct)/) {
			$Data0->{'y0'}      = \@Y2;
			$Data0->{'y0_Name'} = "Solar(Direct)";
		}
		else {
			$Data0->{'x1'}      = \@X;
			$Data0->{'x1_Name'} = "lambda / nm";
			$Data0->{'y0'}      = \@Y1;
			$Data0->{'y0_Name'} = "Solar(Total)";
			$Data0->{'y1'}      = \@Y2;
			$Data0->{'y1_Name'} = "Solar(Direct)";
		}
	}

	$Data0->CalMinMax();
	$GraphArray->AddGraphData($Data0);
	return $this->{LightSource} = $GraphArray;
}

sub ReadColorMatchingFunctions
{
#$Type: XYZ, X10Y10Z10
	my ($this, $Type) = @_;
	my $FileName = $ColorMatchingFunctionPath;
	$Type = 'XYZ' if(!defined $Type);

	my $GraphArray = new GraphDataArray;
	my $Data0 = new GraphData;
	$Data0->SetTitle("Color Matching Functions");

	my (@X, @x, @y, @z, @x10, @y10, @z10);

	my $in = new JFile;
	if(!$in->Open($FileName, "r")) {
		return undef;
	}
	my $FirstLine = $in->ReadLine();
	my $count = 0;
	while(!$in->eof()) {
		my $line = $in->ReadLine();
		last if($in->eof());
		my @a = Utils::Split("\\s+", $line);
print "$a[0]\t$a[1]\t$a[2]\t$a[3]\t$a[4]\t$a[5]\t$a[6]\n";
		$X[$count]   = $a[0];
		$x[$count]   = $a[1];
		$y[$count]   = $a[2];
		$z[$count]   = $a[3];
		$x10[$count] = $a[4];
		$y10[$count] = $a[5];
		$z10[$count] = $a[6];
		$count++;
	}
	$in->Close();

	$Data0->{'x0'}      = \@X;
	$Data0->{'x0_Name'} = "lambda / nm";
	$Data0->{'x1'}      = \@X;
	$Data0->{'x1_Name'} = "lambda / nm";
	$Data0->{'x2'}      = \@X;
	$Data0->{'x2_Name'} = "lambda / nm";
	if($Type =~ /X10/i) {
		$Data0->{'y0'}      = \@x10;
		$Data0->{'y0_Name'} = "x10";
		$Data0->{'y1'}      = \@y10;
		$Data0->{'y1_Name'} = "y10";
		$Data0->{'y2'}      = \@z10;
		$Data0->{'y2_Name'} = "z10";
	}
	else {
		$Data0->{'y0'}      = \@x;
		$Data0->{'y0_Name'} = "x";
		$Data0->{'y1'}      = \@y;
		$Data0->{'y1_Name'} = "y";
		$Data0->{'y2'}      = \@z;
		$Data0->{'y2_Name'} = "z";
	}

	$Data0->CalMinMax();
	$GraphArray->AddGraphData($Data0);

	$this->{ColorMatchingFunctionType} = $Type;
	return $this->{ColorMatchingFunction} = $GraphArray;
}

1;