src/px/falsecolor.pl

#!/usr/bin/perl -w
# RCSid $Id: falsecolor.pl,v 2.2 2011/01/04 17:00:02 greg Exp $

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 $picture = '-';
my $cpict = '';
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 $needfile = 0;

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;
        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;

        # Switches
    } elsif ("$ARGV[0]" eq '-cl') {    # Contour lines
        $docont = 'a';
        $loff = 12;
    } elsif ("$ARGV[0]" eq '-cb') {    # Contour bands
        $docont = 'b';
        $loff = 13;
    } elsif ("$ARGV[0]" eq '-e') {
        $doextrem = 1;
        $needfile = 1;

        # Oops! Illegal option
    } else {
        die("bad option \"$ARGV[0]\"\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.
    $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;
}

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*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;

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

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

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

# Contour lines or bands
if ($docont ne '') {
    $pc0args .= " -e 'in=iscont$docont'";
}

if ($cpict eq '') {
    $pc1args .= " -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)'";
}

# 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
    my $sheight = floor($legheight / $ndivs + 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";
    }
    $cmd = "echo '$text' | psign -s -.15 -cf 1 1 1 -cb 0 0 0";
    $cmd .= " -h $sheight > $slabpic";
    system $cmd;
} else {
        # 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.
    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 = "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";

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 $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;

        # Add extrema labels to command line
    $cmd .= " $minvpic $minpos $maxvpic $maxpos";
}

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

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