#!/usr/bin/perl

use strict;
use Image::Magick;

$| = 1;

my $T0 = 1.0/400;
my $A0 = 3.5;
my $T1 = 1.0/400;
my $A1 = 6.3;
my $T2 = 1.0/80;
my $A2 = 3.5;

my $E01 = ($T0/$T1)*(($A1*$A1)/($A0*$A0));
my $E02 = ($T0/$T2)*(($A2*$A2)/($A0*$A0));
#print "0 => 1: $E01; 0 => 2: $E02\n";
#$E01 = 3.1;
#$E01 = 2.0;
#$E02 = 0.5;

my $img0 = Image::Magick->new;
$img0->Read('000000.tif');
$img0->Resize(geometry => '80%', support => 1.1);
my $img1 = Image::Magick->new;
$img1->Read('000001.tif');
$img1->Resize(geometry => '80%', support => 1.1);
my $img2 = Image::Magick->new;
$img2->Read('000002.tif');
$img2->Resize(geometry => '80%', support => 1.1);
#my $img0 = Image::Magick->new;
#$img0->Read('../../2006/11/IMG_0934.JPG');
#my $img1 = Image::Magick->new;
#$img1->Read('../../2006/11/IMG_0935.JPG');
#my $img2 = Image::Magick->new;
#$img2->Read('../../2006/11/IMG_0936.JPG');

sub w {
    #32768/(16384+abs(32768-$_[0]));
    1;
}

## Create output
my $out = $img0->Clone();
$out->Set('depth',16);
my $w = $out->Get('width');
my $h = $out->Get('height');
print "${w}x$h\n";

if (0) {
## Crop output
my $lx = 0;
my $ty = 0;
my $rx = $w;
my $by = $h;
for (my $y = 0; $y < $h; ++$y) {
    for (my $x = 0; $x < $w; ++$x) {
        my ($r,$g,$b,$o) = split /,/, $out->Get("pixel[$x,$y]");
        next if $o != 65535;
        $lx = $x if $x > $lx;
        $ty = $y if $y > $ty;
        $rx = $x if $x < $rx;
        $by = $y if $y < $by;
    }
}
my $g = ($rx-$lx)."x".($by-$ty)."+$lx+$ty";
print "${w}x$h => $g\n";
$out->Crop(geometry => $g);
exit 0;
}

#my ($C,$R01,$G01,$B01,$R201,$G201,$B201,$R02,$G02,$B02,$R202,$G202,$B202) = 
#    (0,0,0,0,0,0,0,0,0,0,0,0,0);
for (my $y = 0; $y < $h; ++$y) {
    for (my $x = 0; $x < $w; ++$x) {
        my ($r0,$g0,$b0,$o0) = split /,/, $img0->Get("pixel[$x,$y]");
        my ($r1,$g1,$b1,$o1) = split /,/, $img1->Get("pixel[$x,$y]");
        my ($r2,$g2,$b2,$o2) = split /,/, $img2->Get("pixel[$x,$y]");
        if ($o1 == 65535 || $o2 == 65535) {
            $out->Set("pixel[$x,$y]",'0,0,0,65535');
            next;
        }
        #if ($r1 && $g1 && $b1 && $r2 && $g2 && $b2) {
        #    ++$C;
        #    my $r01 = $r0/$r1;
        #    my $g01 = $g0/$g1;
        #    my $b01 = $b0/$b1;
        #    my $r02 = $r0/$r2;
        #    my $g02 = $g0/$g2;
        #    my $b02 = $b0/$b2;
        #    $R01 += $r01; $R201 += $r01*$r01;
        #    $G01 += $g01; $G201 += $g01*$g01;
        #    $B01 += $b01; $B201 += $b01*$b01;
        #    $R02 += $r02; $R202 += $r02*$r02;
        #    $G02 += $g02; $G202 += $g02*$g02;
        #    $B02 += $b02; $B202 += $b02*$b02;
        #}
        #print "${x}x$y: ($r0,$g0,$b0,$o0); ($r1,$g1,$b1,$o1); ($r2,$g2,$b2,$o2)\n";
        #$r1 *= $E01;
        #$g1 *= $E01;
        #$b1 *= $E01;
        #if (abs($r0-$r1) > 128 || abs($g0-$g1) > 128 || abs($b0-$b1) > 128) {
        #    $out->Set("pixel[$x,$y]",'0,0,0,0');
        #    print "${x}x$y (1): ($r0,$g0,$b0); ($r1,$g1,$b1)\n";
        #    next;
        #}
        #$r2 *= $E02;
        #$g2 *= $E02;
        #$b2 *= $E02;
        #if (abs($r0-$r2) > 128 || abs($g0-$g2) > 128 || abs($b0-$b2) > 128) {
        #    $out->Set("pixel[$x,$y]",'65535,65535,65535,0');
        #    print "${x}x$y (1): ($r0,$g0,$b0); ($r2,$g2,$b2)\n";
        #    next;
        #}
        my $w0 = w($r0);
        my $w1 = w($r1);
        my $w2 = w($r2);
        my $r = int(($r0*$w0+$r1*$w1*$E01+$r2*$w2*$E02)/($w0+$w1+$w2));
        $r = 65535 if $r > 65535;
        $w0 = w($g0);
        $w1 = w($g1);
        $w2 = w($g2);
        my $g = int(($g0*$w0+$g1*$w1*$E01+$g2*$w2*$E02)/($w0+$w1+$w2));
        $g = 65535 if $g > 65535;
        $w0 = w($b0);
        $w1 = w($b1);
        $w2 = w($b2);
        my $b = int(($b0*$w0+$b1*$w1*$E01+$b2*$w2*$E02)/($w0+$w1+$w2));
        $b = 65535 if $b > 65535;
        $out->Set("pixel[$x,$y]","$r,$g,$b,0");
    }
}
#printf 
#    "R01=%f (%f); G01=%f (%f); B01=%f (%f)\nR02=%f (%f); G02=%f (%f); B02=%f (%f)\n",
#    $R01/$C, sqrt(($R201-$R01*$R01/$C)/$C), 
#    $G01/$C, sqrt(($G201-$G01*$G01/$C)/$C), 
#    $B01/$C, sqrt(($B201-$B01*$B01/$C)/$C), 
#    $R02/$C, sqrt(($R202-$R02*$R02/$C)/$C), 
#    $G02/$C, sqrt(($G202-$G02*$G02/$C)/$C), 
#    $B02/$C, sqrt(($B202-$B02*$B02/$C)/$C); 

$out->Write('00.tif');
