--- ray/src/px/falsecolor.pl	2011/01/04 17:00:02	2.2
+++ ray/src/px/falsecolor.pl	2014/04/15 21:34:31	2.10
@@ -1,113 +1,129 @@
 #!/usr/bin/perl -w
-# RCSid $Id: falsecolor.pl,v 2.2 2011/01/04 17:00:02 greg Exp $
+# RCSid $Id: falsecolor.pl,v 2.10 2014/04/15 21:34:31 greg Exp $
 
+use warnings;
 use strict;
 use File::Temp qw/ tempdir /;
 use POSIX qw/ floor /;
 
-my $mult = 179.0;    # Multiplier. Default W/sr/m2 -> cd/m2
-my $label = 'cd/m2';    # Units shown in legend
-my $scale = 1000;    # Top of the scale
-my $decades = 0;    # Default is linear mapping
-my $redv = 'def_red(v)';    # Mapping function for R,G,B
-my $grnv = 'def_grn(v)';
-my $bluv = 'def_blu(v)';
-my $ndivs = 8;    # Number of lines in legend
+my @palettes = ('def', 'spec', 'pm3d', 'hot', 'eco');
+
+my $mult = 179.0;              # Multiplier. Default W/sr/m2 -> cd/m2
+my $label = 'cd/m2';           # Units shown in legend
+my $scale = 1000;              # Top of the scale
+my $decades = 0;               # Default is linear mapping
+my $pal = 'def';               # Palette
+my $redv = "${pal}_red(v)";    # Mapping functions for R,G,B
+my $grnv = "${pal}_grn(v)";
+my $bluv = "${pal}_blu(v)";
+my $ndivs = 8;                 # Number of lines in legend
 my $picture = '-';
 my $cpict = '';
-my $legwidth = 100;   # Legend width and height
+my $legwidth = 100;            # Legend width and height
 my $legheight = 200;
-my $docont = '';    # Contours
-my $loff = 0;    # Offset for drop-shadow
-my $doextrem = 0;    # Don't mark extrema
+my $docont = '';               # Contours: -cl and -cb
+my $doposter = 0;              # Posterization: -cp
+my $loff = 0;                  # Offset to align with values
+my $doextrem = 0;              # Don't mark extrema
 my $needfile = 0;
+my $showpal = 0;               # Show availabel colour palettes
 
 while ($#ARGV >= 0) {
-	# Options with qualifiers
-    if ("$ARGV[0]" eq '-lw') {    # Legend width
-        $legwidth = $ARGV[1];
-        shift @ARGV;
-    } elsif ("$ARGV[0]" eq '-lh') {    # Legend height
-        $legheight = $ARGV[1];
-        shift @ARGV;
-    } elsif ("$ARGV[0]" eq '-m') {    # Multiplier
-        $mult = $ARGV[1];
-        shift @ARGV;
-    } elsif ("$ARGV[0]" eq '-s') {    # Scale
-        $scale = $ARGV[1];
-        shift @ARGV;
+    $_ = shift;
+    # Options with qualifiers
+    if (m/-lw/) {              # Legend width
+        $legwidth = shift;
+    } elsif (m/-lh/) {         # Legend height
+        $legheight = shift;
+    } elsif (m/-m/) {          # Multiplier
+        $mult = shift;
+    } elsif (m/-spec/) {
+        die("depricated option '-spec'. Please use '-pal spec' instead.");
+    } elsif (m/-s/) {          # Scale
+        $scale = shift;
         if ($scale =~ m/[aA].*/) {
             $needfile = 1;
         }
-    } elsif ("$ARGV[0]" eq '-l') {    # Label
-        $label = $ARGV[1];
-        shift @ARGV;
-    } elsif ("$ARGV[0]" eq '-log') {    # Logarithmic mapping
-        $decades = $ARGV[1];
-        shift @ARGV;
-    } elsif ("$ARGV[0]" eq '-r') {
-        $redv = $ARGV[1];
-        shift @ARGV;
-    } elsif ("$ARGV[0]" eq '-g') {
-        $grnv = $ARGV[1];
-        shift @ARGV;
-    } elsif ("$ARGV[0]" eq '-b') {
-        $bluv = $ARGV[1];
-        shift @ARGV;
-    } elsif ("$ARGV[0]" eq '-pal') {
-	$redv = "$ARGV[1]_red(v)";
-	$grnv = "$ARGV[1]_grn(v)";
-	$bluv = "$ARGV[1]_blu(v)";
-	shift @ARGV;
-    } elsif ("$ARGV[0]" eq '-i') {    # Image for intensity mapping
-        $picture = $ARGV[1];
-        shift @ARGV;
-    } elsif ("$ARGV[0]" eq '-p') {    # Image for background
-        $cpict = $ARGV[1];
-        shift @ARGV;
-    } elsif ("$ARGV[0]" eq '-ip' || "$ARGV[0]" eq '-pi') {
-        $picture = $ARGV[1];
-        $cpict = $ARGV[1];
-        shift @ARGV;
-    } elsif ("$ARGV[0]" eq '-n') {    # Number of contour lines
-        $ndivs = $ARGV[1];
-        shift @ARGV;
+    } elsif (m/-l$/) {         # Label
+        $label = shift;
+    } elsif (m/-log/) {        # Logarithmic mapping
+        $decades = shift;
+    } elsif (m/-r/) {          # Custom palette functions for R,G,B
+        $redv = shift;
+    } elsif (m/-g/) {
+        $grnv = shift;
+    } elsif (m/-b/) {
+        $bluv = shift;
+    } elsif (m/-pal$/) {        # Color palette
+        $pal = shift;
+        if (! grep $_ eq $pal, @palettes) {
+            die("invalid palette '$pal'.\n");
+        }
+        $redv = "${pal}_red(v)";
+        $grnv = "${pal}_grn(v)";
+        $bluv = "${pal}_blu(v)";
+    } elsif (m/-i$/) {          # Image for intensity mapping
+        $picture = shift;
+    } elsif (m/-p$/) {         # Image for background
+        $cpict = shift;
+    } elsif (m/-ip/ || m/-pi/) {
+        $picture = shift;
+        $cpict = $picture;
+    } elsif (m/-n/) {          # Number of contour lines
+        $ndivs = shift;
 
-	# Switches
-    } elsif ("$ARGV[0]" eq '-cl') {    # Contour lines
+    # Switches
+    } elsif (m/-cl/) {         # Contour lines
         $docont = 'a';
-        $loff = 12;
-    } elsif ("$ARGV[0]" eq '-cb') {    # Contour bands
+        $loff = 0.48;
+    } elsif (m/-cb/) {         # Contour bands
         $docont = 'b';
-        $loff = 13;
-    } elsif ("$ARGV[0]" eq '-e') {
+        $loff = 0.52;
+    } elsif (m/-cp/) {              # Posterize
+        $doposter = 1;
+    } elsif (m/-palettes/) {        # Show all available palettes
+        $scale   = 45824;           # 256 * 179
+        $showpal = 1;
+    } elsif (m/-e/) {
         $doextrem = 1;
         $needfile = 1;
 
-	# Oops! Illegal option
+    # Oops! Illegal option
     } else {
-        die("bad option \"$ARGV[0]\"\n");
+        die("bad option \"$_\"\n");
     }
-    shift @ARGV;
 }
 
 # Temporary directory. Will be removed upon successful program exit.
 my $td = tempdir( CLEANUP => 1 );
 
 if ($needfile == 1 && $picture eq '-') {
-    $picture = $td . '/' . $picture;
+    # Pretend that $td/stdin.rad is the actual filename.
+    $picture = "$td/stdin.hdr";
     open(FHpic, ">$picture") or
             die("Unable to write to file $picture\n");
+    # Input is from STDIN: Capture to file.
+    while (<>) {
+        print FHpic;
+    }
     close(FHpic);
+
+    if ($cpict eq '-') {
+        $cpict = "$td/stdin.hdr";
+    }
 }
 
 # Find a good scale for auto mode.
 if ($scale =~ m/[aA].*/) {
-    my $LogLmax = `phisto $picture| tail -2| sed -n '1s/[0-9]*\$//p'`;
+    my @histo = split(/\s/, `phisto $picture| tail -2`);
+    # e.g. $ phisto tests/richmond.hdr| tail -2
+    # 3.91267	14
+    # 3.94282	6
+    my $LogLmax = $histo[0];
     $scale = $mult / 179 * 10**$LogLmax;
 }
 
-my $pc0 = $td . '/pc0.cal';
+my $pc0 = "$td/pc0.cal";
 open(FHpc0, ">$pc0");
 print FHpc0 <<EndOfPC0;
 PI : 3.14159265358979323846 ;
@@ -129,13 +145,17 @@ spec_grn(x) = if(x-.375, 1.6-1.6*x, 8/3*x);
 spec_blu(x) = 1 - 8/3*x;
 
 pm3d_red(x) = sqrt(x) ^ gamma;
-pm3d_grn(x) = x*x*x ^ gamma;
-pm3d_blu(x) = clip(sin(2*PI*x)) ^ gamma;
+pm3d_grn(x) = (x*x*x) ^ gamma;
+pm3d_blu(x) = clip(sin(2*PI*clip(x))) ^ gamma;
 
 hot_red(x) = clip(3*x) ^ gamma;
 hot_grn(x) = clip(3*x - 1) ^ gamma;
 hot_blu(x) = clip(3*x - 2) ^ gamma;
 
+eco_red(x) = clip(2*x) ^ gamma;
+eco_grn(x) = clip(2*(x-0.5)) ^ gamma;
+eco_blu(x) = clip(2*(0.5-x)) ^ gamma;
+
 interp_arr2(i,x,f):(i+1-x)*f(i)+(x-i)*f(i+1);
 interp_arr(x,f):if(x-1,if(f(0)-x,interp_arr2(floor(x),x,f),f(f(0))),f(1));
 
@@ -173,8 +193,9 @@ ro = if(in,clip($redv),ra);
 go = if(in,clip($grnv),ga);
 bo = if(in,clip($bluv),ba);
 EndOfPC0
+close FHpc0;
 
-my $pc1 = $td . '/pc1.cal';
+my $pc1 = "$td/pc1.cal";
 open(FHpc1, ">$pc1");
 print FHpc1 <<EndOfPC1;
 norm : mult/scale/le(1);
@@ -192,62 +213,94 @@ ra = ri(nfiles);
 ga = gi(nfiles);
 ba = bi(nfiles);
 EndOfPC1
+close FHpc1;
 
 my $pc0args = "-f $pc0";
 my $pc1args = "-f $pc1";
 
-# Contour lines or bands
+if ($showpal == 1) {
+    my $pc = "pcompos -a 1";
+    foreach my $pal (@palettes) {
+        my $fcimg = "$td/$pal.hdr";
+        my $lbimg = "$td/${pal}_label.hdr";
+        system "psign -cb 0 0 0 -cf 1 1 1 -h 20 $pal > $lbimg";
+
+        my $cmd = qq[pcomb $pc0args -e "v=x/256"];
+        $cmd .= qq[ -e "ro=clip(${pal}_red(v));go=clip(${pal}_grn(v));bo=clip(${pal}_blu(v))"];
+        $cmd .= qq[ -x 256 -y 30 > $fcimg];
+        system "$cmd";
+        $pc .= " $fcimg $lbimg";
+    }
+    system "$pc";
+    exit 0;
+}
+
+# Contours
 if ($docont ne '') {
-    $pc0args .= " -e 'in=iscont$docont'";
+    # -cl -> $docont = a
+    # -cb -> $docont = b
+    $pc0args .= qq[ -e "in=iscont$docont"];
+} elsif ($doposter == 1) {
+    # -cp -> $doposter = 1
+    $pc0args .= qq[ -e "ro=${pal}_red(seg(v));go=${pal}_grn(seg(v));bo=${pal}_blu(seg(v))"];
+    $pc0args .= q[ -e "seg(x)=(floor(v*ndivs)+.5)/ndivs"];
 }
 
 if ($cpict eq '') {
-    $pc1args .= " -e 'ra=0;ga=0;ba=0'";
+    $pc1args .= qq[ -e "ra=0;ga=0;ba=0"];
 } elsif ($cpict eq $picture) {
     $cpict = '';
 }
 
 # Logarithmic mapping
 if ($decades > 0) {
-    $pc1args .= " -e 'map(x)=if(x-10^-$decades,log10(x)/$decades+1,0)'";
+    $pc1args .= qq[ -e "map(x)=if(x-10^-$decades,log10(x)/$decades+1,0)"];
 }
 
 # Colours in the legend
-my $scolpic = $td . '/scol.hdr';
+my $scolpic = "$td/scol.hdr";
+
 # Labels in the legend
-my $slabpic = $td . '/slab.hdr';
+my $slabpic = "$td/slab.hdr";
 my $cmd;
 
-if ($legwidth > 20 && $legheight > 40) {
-	# Legend: Create the background colours
-    $cmd = "pcomb $pc0args -e 'v=(y+.5)/yres;vleft=v;vright=v'";
-    $cmd .= " -e 'vbelow=(y-.5)/yres;vabove=(y+1.5)/yres'";
-    $cmd .= " -x $legwidth -y $legheight > $scolpic";
-    system $cmd;
-
-	# Legend: Create the text labels
+if (($legwidth > 20) && ($legheight > 40)) {
+    # Legend: Create the text labels
     my $sheight = floor($legheight / $ndivs + 0.5);
+    $legheight = $sheight * $ndivs;
+    $loff = floor($loff * $sheight + 0.5);
     my $text = "$label";
-    my $value;
-    my $i;
-    for ($i=0; $i<$ndivs; $i++) {
+    for (my $i=0; $i<$ndivs; $i++) {
         my $imap = ($ndivs - 0.5 - $i) / $ndivs;
+        my $value = $scale;
         if ($decades > 0) {
-            $value = $scale * 10**(($imap - 1) * $decades);
+            $value *= 10**(($imap - 1) * $decades);
         } else {
-            $value = $scale * $imap;
+            $value *= $imap;
         }
+
         # Have no more than 3 decimal places
         $value =~ s/(\.[0-9]{3})[0-9]*/$1/;
         $text .= "\n$value";
     }
-    $cmd = "echo '$text' | psign -s -.15 -cf 1 1 1 -cb 0 0 0";
-    $cmd .= " -h $sheight > $slabpic";
+    open PSIGN, "| psign -s -.15 -cf 1 1 1 -cb 0 0 0 -h $sheight > $slabpic";
+    print PSIGN "$text\n";
+    close PSIGN;
+
+    # Legend: Create the background colours
+    $cmd = qq[pcomb $pc0args];
+    $cmd .= q[ -e "v=(y+.5)/yres;vleft=v;vright=v"];
+    $cmd .= q[ -e "vbelow=(y-.5)/yres;vabove=(y+1.5)/yres"];
+    $cmd .= qq[ -x $legwidth -y $legheight > $scolpic];
     system $cmd;
 } else {
-	# Legend is too small to be legible. Don't bother doing one.
+    # Legend is too small to be legible. Don't bother doing one.
     $legwidth = 0;
     $legheight = 0;
+    $loff = 0;
+
+    # Create dummy colour scale and legend labels so we don't
+    # need to change the final command line.
     open(FHscolpic, ">$scolpic");
     print FHscolpic "\n-Y 1 +X 1\naaa\n";
     close(FHscolpic);
@@ -256,15 +309,21 @@ if ($legwidth > 20 && $legheight > 40) {
     close(FHslabpic);
 }
 
+# Legend: Invert the text labels (for dropshadow)
+my $slabinvpic = "$td/slabinv.hdr";
+$cmd = qq[pcomb -e "lo=1-gi(1)" $slabpic > $slabinvpic];
+system $cmd;
+
 my $loff1 = $loff - 1;
+
 # Command line without extrema
-$cmd = "pcomb $pc0args $pc1args $picture $cpict";
-$cmd .= "| pcompos $scolpic 0 0 +t .1";
-$cmd .= " \"\!pcomb -e 'lo=1-gi(1)' $slabpic\" 2 $loff1";
-$cmd .= " -t .5 $slabpic 0 $loff - $legwidth 0";
+$cmd = qq[pcomb $pc0args $pc1args "$picture"];
+$cmd .= qq[ "$cpict"] if ($cpict);
+$cmd .= qq[ | pcompos $scolpic 0 0 +t .1 $slabinvpic 2 $loff1];
+$cmd .= qq[ -t .5 $slabpic 0 $loff - $legwidth 0];
 
 if ($doextrem == 1) {
-	# Get min/max image luminance
+    # Get min/max image luminance
     my $cmd1 = 'pextrem -o ' . $picture;
     my $retval = `$cmd1`;
     # e.g.
@@ -273,35 +332,25 @@ if ($doextrem == 1) {
     # 211 202 1.292969e+00 1.308594e+00 1.300781e+00
 
     my @extrema = split(/\s/, $retval);
-    my $lxmin = $extrema[0] + $legwidth;
-    my $ymin = $extrema[1];
-    my $rmin = $extrema[2];
-    my $gmin = $extrema[3];
-    my $bmin = $extrema[4];
-    my $lxmax = $extrema[5] + $legwidth;
-    my $ymax = $extrema[6];
-    my $rmax = $extrema[7];
-    my $gmax = $extrema[8];
-    my $bmax = $extrema[9];
+    my ($lxmin, $ymin, $rmin, $gmin, $bmin, $lxmax, $ymax, $rmax, $gmax, $bmax) = @extrema;
+    $lxmin += $legwidth;
+    $lxmax += $legwidth;
 
-	# Weighted average of R,G,B
+    # Weighted average of R,G,B
     my $minpos = "$lxmin $ymin";
     my $minval = ($rmin * .27 + $gmin * .67 + $bmin * .06) * $mult;
     $minval =~ s/(\.[0-9]{3})[0-9]*/$1/;
-    my $maxpos = "$lxmax $ymax";
     my $maxval = ($rmax * .27 + $gmax * .67 + $bmax * .06) * $mult;
     $maxval =~ s/(\.[0-9]{3})[0-9]*/$1/;
 
-	# Create the labels for min/max intensity
-    my $minvpic = $td . '/minv.hdr';
-    $cmd1 = "psign -s -.15 -a 2 -h 16 $minval > $minvpic";
-    system $cmd1;
-    my $maxvpic = $td . '/maxv.hdr';
-    $cmd1 = "psign -s -.15 -a 2 -h 16 $maxval > $maxvpic";
-    system $cmd1;
+    # Create the labels for min/max intensity
+    my $minvpic = "$td/minv.hdr";
+    system "psign -s -.15 -a 2 -h 16 $minval > $minvpic";
+    my $maxvpic = "$td/maxv.hdr";
+    system "psign -s -.15 -a 2 -h 16 $maxval > $maxvpic";
 
-	# Add extrema labels to command line
-    $cmd .= " $minvpic $minpos $maxvpic $maxpos";
+    # Add extrema labels to command line
+    $cmd .= qq[ $minvpic $minpos $maxvpic $lxmax $ymax];
 }
 
 # Process image and combine with legend