src/px/falsecolor.pl

#!/usr/bin/perl -w
# 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 @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 $matchscale = 1;            # Adjust top of scale to match given label
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 $legheight = 200;
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) {
    $_ = 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].*/) {
            $matchscale = 0;
            $needfile = 1;
       } else {
            $matchscale = 1;
        }
    } 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 (m/-cl/) {         # Contour lines
        $docont = 'a';
        $loff = 0.48;
    } elsif (m/-cb/) {         # Contour bands
        $docont = 'b';
        $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
    } else {
        die("bad option \"$_\"\n");
    }
}

if (($legwidth <= 20) || ($legheight <= 40)) {
    # Legend is too small to be legible. Don't bother doing one.
    $legwidth = 0;
    $legheight = 0;
    $loff = 0;
    $matchscale = 0;
}

# 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.
    $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 @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;
} elsif ($matchscale) {
    # Adjust scale so legend reflects -s setting
    if ($decades > 0) {
        $scale *= 10**($decades/(2.*$ndivs));
    } else {
        $scale *= $ndivs/($ndivs - .5);
    }
}

my $pc0 = "$td/pc0.cal";
open(FHpc0, ">$pc0");
print FHpc0 <<EndOfPC0;
PI : 3.14159265358979323846 ;
scale : $scale ;
mult : $mult ;
ndivs : $ndivs ;
gamma : 2.2;

or(a,b) : if(a,a,b);
EPS : 1e-7;
neq(a,b) : if(a-b-EPS,1,b-a-EPS);
btwn(a,x,b) : if(a-x,-1,b-x);
clip(x) : if(x-1,1,if(x,x,0));
frac(x) : x - floor(x);
boundary(a,b) : neq(floor(ndivs*a+.5),floor(ndivs*b+.5));

spec_red(x) = 1.6*x - .6;
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*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));

def_redp(i):select(i,0.18848,0.05468174,
0.00103547,8.311144e-08,7.449763e-06,0.0004390987,0.001367254,
0.003076,0.01376382,0.06170773,0.1739422,0.2881156,0.3299725,
0.3552663,0.372552,0.3921184,0.4363976,0.6102754,0.7757267,
0.9087369,1,1,0.9863);
def_red(x):interp_arr(x/0.0454545+1,def_redp);

def_grnp(i):select(i,0.0009766,2.35501e-05,
0.0008966244,0.0264977,0.1256843,0.2865799,0.4247083,0.4739468,
0.4402732,0.3671876,0.2629843,0.1725325,0.1206819,0.07316644,
0.03761026,0.01612362,0.004773749,6.830967e-06,0.00803605,
0.1008085,0.3106831,0.6447838,0.9707);
def_grn(x):interp_arr(x/0.0454545+1,def_grnp);

def_blup(i):select(i,0.2666,0.3638662,0.4770437,
0.5131397,0.5363797,0.5193677,0.4085123,0.1702815,0.05314236,
0.05194055,0.08564082,0.09881395,0.08324373,0.06072902,
0.0391076,0.02315354,0.01284458,0.005184709,0.001691774,
2.432735e-05,1.212949e-05,0.006659406,0.02539);
def_blu(x):interp_arr(x/0.0454545+1,def_blup);

isconta = if(btwn(0,v,1),or(boundary(vleft,vright),boundary(vabove,vbelow)),-1);
iscontb = if(btwn(0,v,1),btwn(.4,frac(ndivs*v),.6),-1);

ra = 0;
ga = 0;
ba = 0;

in = 1;

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");
print FHpc1 <<EndOfPC1;
norm : mult/scale/le(1);

v = map(li(1)*norm);

vleft = map(li(1,-1,0)*norm);
vright = map(li(1,1,0)*norm);
vabove = map(li(1,0,1)*norm);
vbelow = map(li(1,0,-1)*norm);

map(x) = x;

ra = ri(nfiles);
ga = gi(nfiles);
ba = bi(nfiles);
EndOfPC1
close FHpc1;

my $pc0args = "-f $pc0";
my $pc1args = "-f $pc1";

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 '') {
    # -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 .= qq[ -e "ra=0;ga=0;ba=0"];
} elsif ($cpict eq $picture) {
    $cpict = '';
}

# Logarithmic mapping
if ($decades > 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 > 0) {
    # 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;
        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";
    }
    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 {
    # 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);
    open(FHslabpic, ">$slabpic");
    print FHslabpic "\n-Y 1 +X 1\naaa\n";
    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 = 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
    my $cmd1 = 'pextrem -o ' . $picture;
    my $retval = `$cmd1`;
    # e.g.
    # x   y   r            g            b
    # 193 207 3.070068e-02 3.118896e-02 1.995850e-02
    # 211 202 1.292969e+00 1.308594e+00 1.300781e+00

    my @extrema = split(/\s/, $retval);
    my ($lxmin, $ymin, $rmin, $gmin, $bmin, $lxmax, $ymax, $rmax, $gmax, $bmax) = @extrema;
    $lxmin += $legwidth;
    $lxmax += $legwidth;

    # 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 $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";
    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 .= qq[ $minvpic $minpos $maxvpic $lxmax $ymax];
}

# Process image and combine with legend
system "$cmd";

#EOF
Revision:	2.11
Committed:	Tue Sep 5 22:42:52 2017 UTC (8 years, 1 month ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 2.10:	+22 -8 lines
Log Message:	Altered falsecolor behavior to match given scale (-s) value in legend
#	User	Rev	Content
1	greg	2.1	#!/usr/bin/perl -w
2	greg	2.11	# RCSid $Id: falsecolor.pl,v 2.10 2014/04/15 21:34:31 greg Exp $
3	greg	2.1
4	greg	2.5	use warnings;
5	greg	2.1	use strict;
6			use File::Temp qw/ tempdir /;
7			use POSIX qw/ floor /;
8
9	greg	2.9	my @palettes = ('def', 'spec', 'pm3d', 'hot', 'eco');
10	greg	2.5
11			my $mult = 179.0; # Multiplier. Default W/sr/m2 -> cd/m2
12			my $label = 'cd/m2'; # Units shown in legend
13			my $scale = 1000; # Top of the scale
14	greg	2.11	my $matchscale = 1; # Adjust top of scale to match given label
15	greg	2.5	my $decades = 0; # Default is linear mapping
16			my $pal = 'def'; # Palette
17			my $redv = "${pal}_red(v)"; # Mapping functions for R,G,B
18			my $grnv = "${pal}_grn(v)";
19			my $bluv = "${pal}_blu(v)";
20			my $ndivs = 8; # Number of lines in legend
21	greg	2.1	my $picture = '-';
22			my $cpict = '';
23	greg	2.5	my $legwidth = 100; # Legend width and height
24	greg	2.1	my $legheight = 200;
25	greg	2.6	my $docont = ''; # Contours: -cl and -cb
26			my $doposter = 0; # Posterization: -cp
27	greg	2.5	my $loff = 0; # Offset to align with values
28			my $doextrem = 0; # Don't mark extrema
29	greg	2.1	my $needfile = 0;
30	greg	2.6	my $showpal = 0; # Show availabel colour palettes
31	greg	2.1
32			while ($#ARGV >= 0) {
33	greg	2.5	$_ = shift;
34			# Options with qualifiers
35			if (m/-lw/) { # Legend width
36			$legwidth = shift;
37			} elsif (m/-lh/) { # Legend height
38			$legheight = shift;
39			} elsif (m/-m/) { # Multiplier
40			$mult = shift;
41	greg	2.9	} elsif (m/-spec/) {
42			die("depricated option '-spec'. Please use '-pal spec' instead.");
43	greg	2.5	} elsif (m/-s/) { # Scale
44			$scale = shift;
45	greg	2.1	if ($scale =~ m/[aA].*/) {
46	greg	2.11	$matchscale = 0;
47	greg	2.1	$needfile = 1;
48	greg	2.11	} else {
49			$matchscale = 1;
50			}
51	greg	2.5	} elsif (m/-l$/) { # Label
52			$label = shift;
53			} elsif (m/-log/) { # Logarithmic mapping
54			$decades = shift;
55			} elsif (m/-r/) { # Custom palette functions for R,G,B
56			$redv = shift;
57			} elsif (m/-g/) {
58			$grnv = shift;
59			} elsif (m/-b/) {
60			$bluv = shift;
61			} elsif (m/-pal$/) { # Color palette
62			$pal = shift;
63			if (! grep $_ eq $pal, @palettes) {
64			die("invalid palette '$pal'.\n");
65			}
66			$redv = "${pal}_red(v)";
67			$grnv = "${pal}_grn(v)";
68			$bluv = "${pal}_blu(v)";
69			} elsif (m/-i$/) { # Image for intensity mapping
70			$picture = shift;
71			} elsif (m/-p$/) { # Image for background
72			$cpict = shift;
73			} elsif (m/-ip/ \|\| m/-pi/) {
74			$picture = shift;
75			$cpict = $picture;
76			} elsif (m/-n/) { # Number of contour lines
77			$ndivs = shift;
78	greg	2.1
79	greg	2.5	# Switches
80			} elsif (m/-cl/) { # Contour lines
81	greg	2.1	$docont = 'a';
82	greg	2.4	$loff = 0.48;
83	greg	2.5	} elsif (m/-cb/) { # Contour bands
84	greg	2.1	$docont = 'b';
85	greg	2.4	$loff = 0.52;
86	greg	2.6	} elsif (m/-cp/) { # Posterize
87			$doposter = 1;
88			} elsif (m/-palettes/) { # Show all available palettes
89			$scale = 45824; # 256 * 179
90			$showpal = 1;
91	greg	2.5	} elsif (m/-e/) {
92	greg	2.1	$doextrem = 1;
93			$needfile = 1;
94
95	greg	2.5	# Oops! Illegal option
96	greg	2.1	} else {
97	greg	2.5	die("bad option \"$_\"\n");
98	greg	2.1	}
99			}
100
101	greg	2.11	if (($legwidth <= 20) \|\| ($legheight <= 40)) {
102			# Legend is too small to be legible. Don't bother doing one.
103			$legwidth = 0;
104			$legheight = 0;
105			$loff = 0;
106			$matchscale = 0;
107			}
108
109	greg	2.1	# Temporary directory. Will be removed upon successful program exit.
110			my $td = tempdir( CLEANUP => 1 );
111
112			if ($needfile == 1 && $picture eq '-') {
113	greg	2.5	# Pretend that $td/stdin.rad is the actual filename.
114	greg	2.3	$picture = "$td/stdin.hdr";
115	greg	2.1	open(FHpic, ">$picture") or
116			die("Unable to write to file $picture\n");
117	greg	2.5	# Input is from STDIN: Capture to file.
118			while (<>) {
119			print FHpic;
120			}
121	greg	2.1	close(FHpic);
122	greg	2.3
123	greg	2.5	if ($cpict eq '-') {
124	greg	2.6	$cpict = "$td/stdin.hdr";
125	greg	2.5	}
126	greg	2.1	}
127
128			# Find a good scale for auto mode.
129			if ($scale =~ m/[aA].*/) {
130	greg	2.5	my @histo = split(/\s/, `phisto $picture\| tail -2`);
131			# e.g. $ phisto tests/richmond.hdr\| tail -2
132			# 3.91267 14
133			# 3.94282 6
134			my $LogLmax = $histo[0];
135	greg	2.1	$scale = $mult / 179 * 10**$LogLmax;
136	greg	2.11	} elsif ($matchscale) {
137			# Adjust scale so legend reflects -s setting
138			if ($decades > 0) {
139			$scale = 10($decades/(2.$ndivs));
140			} else {
141			$scale *= $ndivs/($ndivs - .5);
142			}
143	greg	2.1	}
144
145	greg	2.3	my $pc0 = "$td/pc0.cal";
146	greg	2.1	open(FHpc0, ">$pc0");
147			print FHpc0 <<EndOfPC0;
148			PI : 3.14159265358979323846 ;
149			scale : $scale ;
150			mult : $mult ;
151			ndivs : $ndivs ;
152			gamma : 2.2;
153
154			or(a,b) : if(a,a,b);
155			EPS : 1e-7;
156			neq(a,b) : if(a-b-EPS,1,b-a-EPS);
157			btwn(a,x,b) : if(a-x,-1,b-x);
158			clip(x) : if(x-1,1,if(x,x,0));
159			frac(x) : x - floor(x);
160			boundary(a,b) : neq(floor(ndivsa+.5),floor(ndivsb+.5));
161
162			spec_red(x) = 1.6*x - .6;
163			spec_grn(x) = if(x-.375, 1.6-1.6x, 8/3x);
164			spec_blu(x) = 1 - 8/3*x;
165
166			pm3d_red(x) = sqrt(x) ^ gamma;
167	greg	2.5	pm3d_grn(x) = (xxx) ^ gamma;
168	greg	2.6	pm3d_blu(x) = clip(sin(2PIclip(x))) ^ gamma;
169	greg	2.1
170			hot_red(x) = clip(3*x) ^ gamma;
171			hot_grn(x) = clip(3*x - 1) ^ gamma;
172			hot_blu(x) = clip(3*x - 2) ^ gamma;
173
174	greg	2.9	eco_red(x) = clip(2*x) ^ gamma;
175			eco_grn(x) = clip(2*(x-0.5)) ^ gamma;
176			eco_blu(x) = clip(2*(0.5-x)) ^ gamma;
177
178	greg	2.1	interp_arr2(i,x,f):(i+1-x)f(i)+(x-i)f(i+1);
179			interp_arr(x,f):if(x-1,if(f(0)-x,interp_arr2(floor(x),x,f),f(f(0))),f(1));
180
181			def_redp(i):select(i,0.18848,0.05468174,
182			0.00103547,8.311144e-08,7.449763e-06,0.0004390987,0.001367254,
183			0.003076,0.01376382,0.06170773,0.1739422,0.2881156,0.3299725,
184			0.3552663,0.372552,0.3921184,0.4363976,0.6102754,0.7757267,
185			0.9087369,1,1,0.9863);
186			def_red(x):interp_arr(x/0.0454545+1,def_redp);
187
188			def_grnp(i):select(i,0.0009766,2.35501e-05,
189			0.0008966244,0.0264977,0.1256843,0.2865799,0.4247083,0.4739468,
190			0.4402732,0.3671876,0.2629843,0.1725325,0.1206819,0.07316644,
191			0.03761026,0.01612362,0.004773749,6.830967e-06,0.00803605,
192			0.1008085,0.3106831,0.6447838,0.9707);
193			def_grn(x):interp_arr(x/0.0454545+1,def_grnp);
194
195			def_blup(i):select(i,0.2666,0.3638662,0.4770437,
196			0.5131397,0.5363797,0.5193677,0.4085123,0.1702815,0.05314236,
197			0.05194055,0.08564082,0.09881395,0.08324373,0.06072902,
198			0.0391076,0.02315354,0.01284458,0.005184709,0.001691774,
199			2.432735e-05,1.212949e-05,0.006659406,0.02539);
200			def_blu(x):interp_arr(x/0.0454545+1,def_blup);
201
202			isconta = if(btwn(0,v,1),or(boundary(vleft,vright),boundary(vabove,vbelow)),-1);
203			iscontb = if(btwn(0,v,1),btwn(.4,frac(ndivs*v),.6),-1);
204
205			ra = 0;
206			ga = 0;
207			ba = 0;
208
209			in = 1;
210
211			ro = if(in,clip($redv),ra);
212			go = if(in,clip($grnv),ga);
213			bo = if(in,clip($bluv),ba);
214			EndOfPC0
215	greg	2.6	close FHpc0;
216	greg	2.1
217	greg	2.3	my $pc1 = "$td/pc1.cal";
218	greg	2.1	open(FHpc1, ">$pc1");
219			print FHpc1 <<EndOfPC1;
220			norm : mult/scale/le(1);
221
222			v = map(li(1)*norm);
223
224			vleft = map(li(1,-1,0)*norm);
225			vright = map(li(1,1,0)*norm);
226			vabove = map(li(1,0,1)*norm);
227			vbelow = map(li(1,0,-1)*norm);
228
229			map(x) = x;
230
231			ra = ri(nfiles);
232			ga = gi(nfiles);
233			ba = bi(nfiles);
234			EndOfPC1
235	greg	2.6	close FHpc1;
236	greg	2.1
237			my $pc0args = "-f $pc0";
238			my $pc1args = "-f $pc1";
239
240	greg	2.6	if ($showpal == 1) {
241			my $pc = "pcompos -a 1";
242			foreach my $pal (@palettes) {
243			my $fcimg = "$td/$pal.hdr";
244			my $lbimg = "$td/${pal}_label.hdr";
245			system "psign -cb 0 0 0 -cf 1 1 1 -h 20 $pal > $lbimg";
246
247			my $cmd = qq[pcomb $pc0args -e "v=x/256"];
248			$cmd .= qq[ -e "ro=clip(${pal}_red(v));go=clip(${pal}_grn(v));bo=clip(${pal}_blu(v))"];
249			$cmd .= qq[ -x 256 -y 30 > $fcimg];
250			system "$cmd";
251			$pc .= " $fcimg $lbimg";
252			}
253			system "$pc";
254			exit 0;
255			}
256
257			# Contours
258	greg	2.1	if ($docont ne '') {
259	greg	2.6	# -cl -> $docont = a
260			# -cb -> $docont = b
261			$pc0args .= qq[ -e "in=iscont$docont"];
262			} elsif ($doposter == 1) {
263			# -cp -> $doposter = 1
264			$pc0args .= qq[ -e "ro=${pal}_red(seg(v));go=${pal}_grn(seg(v));bo=${pal}_blu(seg(v))"];
265			$pc0args .= q[ -e "seg(x)=(floor(v*ndivs)+.5)/ndivs"];
266	greg	2.1	}
267
268			if ($cpict eq '') {
269	greg	2.6	$pc1args .= qq[ -e "ra=0;ga=0;ba=0"];
270	greg	2.1	} elsif ($cpict eq $picture) {
271			$cpict = '';
272			}
273
274			# Logarithmic mapping
275			if ($decades > 0) {
276	greg	2.6	$pc1args .= qq[ -e "map(x)=if(x-10^-$decades,log10(x)/$decades+1,0)"];
277	greg	2.1	}
278
279			# Colours in the legend
280	greg	2.3	my $scolpic = "$td/scol.hdr";
281	greg	2.6
282	greg	2.1	# Labels in the legend
283	greg	2.3	my $slabpic = "$td/slab.hdr";
284	greg	2.1	my $cmd;
285
286	greg	2.11	if ($legwidth > 0) {
287	greg	2.5	# Legend: Create the text labels
288	greg	2.1	my $sheight = floor($legheight / $ndivs + 0.5);
289	greg	2.4	$legheight = $sheight * $ndivs;
290			$loff = floor($loff * $sheight + 0.5);
291	greg	2.1	my $text = "$label";
292	greg	2.3	for (my $i=0; $i<$ndivs; $i++) {
293	greg	2.1	my $imap = ($ndivs - 0.5 - $i) / $ndivs;
294	greg	2.5	my $value = $scale;
295	greg	2.1	if ($decades > 0) {
296	greg	2.3	$value = 10(($imap - 1) $decades);
297	greg	2.1	} else {
298	greg	2.3	$value *= $imap;
299	greg	2.1	}
300	greg	2.6
301	greg	2.1	# Have no more than 3 decimal places
302			$value =~ s/(\.[0-9]{3})[0-9]*/$1/;
303			$text .= "\n$value";
304			}
305	greg	2.8	open PSIGN, "\| psign -s -.15 -cf 1 1 1 -cb 0 0 0 -h $sheight > $slabpic";
306			print PSIGN "$text\n";
307			close PSIGN;
308	greg	2.4
309	greg	2.5	# Legend: Create the background colours
310	greg	2.6	$cmd = qq[pcomb $pc0args];
311			$cmd .= q[ -e "v=(y+.5)/yres;vleft=v;vright=v"];
312			$cmd .= q[ -e "vbelow=(y-.5)/yres;vabove=(y+1.5)/yres"];
313			$cmd .= qq[ -x $legwidth -y $legheight > $scolpic];
314	greg	2.4	system $cmd;
315	greg	2.1	} else {
316	greg	2.5	# Create dummy colour scale and legend labels so we don't
317			# need to change the final command line.
318	greg	2.1	open(FHscolpic, ">$scolpic");
319	greg	2.2	print FHscolpic "\n-Y 1 +X 1\naaa\n";
320	greg	2.1	close(FHscolpic);
321			open(FHslabpic, ">$slabpic");
322	greg	2.2	print FHslabpic "\n-Y 1 +X 1\naaa\n";
323	greg	2.1	close(FHslabpic);
324			}
325
326	greg	2.3	# Legend: Invert the text labels (for dropshadow)
327			my $slabinvpic = "$td/slabinv.hdr";
328	greg	2.6	$cmd = qq[pcomb -e "lo=1-gi(1)" $slabpic > $slabinvpic];
329	greg	2.3	system $cmd;
330
331	greg	2.1	my $loff1 = $loff - 1;
332	greg	2.6
333	greg	2.1	# Command line without extrema
334	greg	2.10	$cmd = qq[pcomb $pc0args $pc1args "$picture"];
335			$cmd .= qq[ "$cpict"] if ($cpict);
336	greg	2.6	$cmd .= qq[ \| pcompos $scolpic 0 0 +t .1 $slabinvpic 2 $loff1];
337			$cmd .= qq[ -t .5 $slabpic 0 $loff - $legwidth 0];
338	greg	2.1
339			if ($doextrem == 1) {
340	greg	2.5	# Get min/max image luminance
341	greg	2.1	my $cmd1 = 'pextrem -o ' . $picture;
342			my $retval = `$cmd1`;
343			# e.g.
344			# x y r g b
345			# 193 207 3.070068e-02 3.118896e-02 1.995850e-02
346			# 211 202 1.292969e+00 1.308594e+00 1.300781e+00
347
348			my @extrema = split(/\s/, $retval);
349	greg	2.3	my ($lxmin, $ymin, $rmin, $gmin, $bmin, $lxmax, $ymax, $rmax, $gmax, $bmax) = @extrema;
350	greg	2.5	$lxmin += $legwidth;
351			$lxmax += $legwidth;
352	greg	2.1
353	greg	2.5	# Weighted average of R,G,B
354	greg	2.1	my $minpos = "$lxmin $ymin";
355			my $minval = ($rmin * .27 + $gmin * .67 + $bmin * .06) * $mult;
356			$minval =~ s/(\.[0-9]{3})[0-9]*/$1/;
357			my $maxval = ($rmax * .27 + $gmax * .67 + $bmax * .06) * $mult;
358			$maxval =~ s/(\.[0-9]{3})[0-9]*/$1/;
359
360	greg	2.5	# Create the labels for min/max intensity
361	greg	2.3	my $minvpic = "$td/minv.hdr";
362	greg	2.6	system "psign -s -.15 -a 2 -h 16 $minval > $minvpic";
363	greg	2.3	my $maxvpic = "$td/maxv.hdr";
364	greg	2.6	system "psign -s -.15 -a 2 -h 16 $maxval > $maxvpic";
365	greg	2.1
366	greg	2.5	# Add extrema labels to command line
367	greg	2.6	$cmd .= qq[ $minvpic $minpos $maxvpic $lxmax $ymax];
368	greg	2.1	}
369
370			# Process image and combine with legend
371			system "$cmd";
372
373			#EOF