--- ray/src/px/falsecolor.pl 2011/03/23 19:33:29 2.3 +++ ray/src/px/falsecolor.pl 2013/04/04 02:59:20 2.9 @@ -1,122 +1,125 @@ #!/usr/bin/perl -w -# RCSid $Id: falsecolor.pl,v 2.3 2011/03/23 19:33:29 greg Exp $ +# RCSid $Id: falsecolor.pl,v 2.9 2013/04/04 02:59:20 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 '-') { - # Pretend that $td/stdin.rad is the actual filename. + # 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; - } + # Input is from STDIN: Capture to file. + while (<>) { + print FHpic; + } close(FHpic); - if ($cpict eq '-') { - $cpict = "$td/stdin.hdr"; - } + if ($cpict eq '-') { + $cpict = "$td/stdin.hdr"; + } } # Find a good scale for auto mode. if ($scale =~ m/[aA].*/) { - 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]; + 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; } @@ -142,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)); @@ -186,6 +193,7 @@ 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"; open(FHpc1, ">$pc1"); @@ -205,63 +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"; + # Labels in the legend 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 + # Legend: Create the text labels my $sheight = floor($legheight / $ndivs + 0.5); + $legheight = $sheight * $ndivs; + $loff = floor($loff * $sheight + 0.5); my $text = "$label"; for (my $i=0; $i<$ndivs; $i++) { my $imap = ($ndivs - 0.5 - $i) / $ndivs; - my $value = $scale; + my $value = $scale; if ($decades > 0) { $value *= 10**(($imap - 1) * $decades); } else { $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; - # Create dummy colour scale and legend labels so we don't - # need to change the final command line. + $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); @@ -272,17 +311,18 @@ if (($legwidth > 20) && ($legheight > 40)) { # Legend: Invert the text labels (for dropshadow) my $slabinvpic = "$td/slabinv.hdr"; -$cmd = "pcomb -e 'lo=1-gi(1)' $slabpic > $slabinvpic"; +$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 $slabinvpic 2 $loff1"; -$cmd .= " -t .5 $slabpic 0 $loff - $legwidth 0"; +$cmd = qq[pcomb $pc0args $pc1args $picture $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. @@ -292,27 +332,24 @@ if ($doextrem == 1) { my @extrema = split(/\s/, $retval); my ($lxmin, $ymin, $rmin, $gmin, $bmin, $lxmax, $ymax, $rmax, $gmax, $bmax) = @extrema; - $lxmin += $legwidth; - $lxmax += $legwidth; + $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 + # Create the labels for min/max intensity my $minvpic = "$td/minv.hdr"; - $cmd1 = "psign -s -.15 -a 2 -h 16 $minval > $minvpic"; - system $cmd1; + system "psign -s -.15 -a 2 -h 16 $minval > $minvpic"; my $maxvpic = "$td/maxv.hdr"; - $cmd1 = "psign -s -.15 -a 2 -h 16 $maxval > $maxvpic"; - system $cmd1; + 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