#Original Function Location: /themis/lib/dav_lib/public_library/latest/image_stretch.dvrc
/themis/lib/dav_lib/library/image_stretch.dvrc
Source Code for Function: "get_color()" in "image_stretch.dvrc" (Public)
image_stretch_version=1.03
# dcs - taken from /u/gorelick/dcs.dv
# do_dcs - taken from /u/knowicki/dvrc/tools.dvrc
# do_rgb - taken from /u/knowicki/dvrc/tools.dvrc
# do_rdcs - written to keep in line with other "do_" tools
# do_hsv - written to keep in line with other "do_" tools
# sstretch_old - taken from /u/gorelick/.dvrc
# nstretch - taken from /u/gorelick/.dvrc
# stretch - taken from /u/knowicki/dvrc/stretch.dvrc
# cleandcs - wrapper for thm.cleandcs
# rdcs - taken from /u/cedwards/rstretch2.dvrc
# rsstretch - taken from /u/cedwards/rstretch2.dvrc
# rmnoise_pca - taken from /u/cedwards/christopher.dvrc
# rrmnoise_pca - taken from /u/cedwards/rstretch2.dvrc
# colorize - taken from /u/cedwards/colorize.dvrc
# get_color - taken from /u/cedwards/colorize.dvrc
# lstretch - Added by J. Bandfield 1/2008
define dcs(ignore,sample,variance) {
#c.edwards modified 6-24-10 to make the minimum stretched data value 1 and not 0
if ($ARGC != 1) {
printf("\ndcs()\n")
printf("performs a decorrelation stretch of data.\n")
printf("data will be centered around 127 with a sigma of +- 50.\n")
printf("requires either 1 or 3 bands of data to be stretched.\n\n")
printf("Syntax: dcs(input, [ignore], [sample], [variance])\n")
printf("example: pic=dcs(a, ignore=-32768, sample=b)\n")
printf("example: pic=dcs(cat(a[,,8],a[,,7],a[,,5],axis=z),ignore=0)\n")
printf(" 'input' can be any 1 or 3 band float/short/int array.\n")
printf(" 'ignore' is the null data value. Default = 0.\n")
printf(" 'sample' is an array of identical size to the input\n")
printf(" used to define the limits of the stretch.\n")
printf(" 'variance' is a number or set of 3 numbers defining\n")
printf(" the stretch per axis. (sigma) default = 50 \n\n")
printf("last updated Tue Jan 24 21:34:08 MST 2006\n")
printf("written by Noel Gorelick.\n")
return(null)
}
#
# this does a decorrelation stretch
# Data will be centered around 127, with a sigma of +/- 50.
#
# compute the symmetric covariance matrix and eigenvectors of it.
a = $1;
if (HasValue(sample)) {
if (HasValue(ignore)) {
c = covar(sample, axis=z, ignore=ignore);
#make the mask in this stupid way because of how ignores are handled
mask=byte(a[,,1]*0+1)
mask[where sum(a == ignore,axis=z)>=1 ]=0
} else {
c = covar(sample, axis=z);
mask=byte(a*0+1)
mask[where sum(a == 0,axis=z)>=1]=0
}
} else {
if (HasValue(ignore)) {
c = covar(a, axis=z, ignore=ignore);
mask=byte(a[,,1]*0+1)
mask[where sum(a == ignore,axis=z)>=1 ]=0
} else {
c = covar(a, axis=z);
mask=byte(a*0+1)
mask[where sum(a == 0,axis=z)>=1]=0
}
}
e = eigen(c);
# generate the scaling matrix with the scaling factors on the diagonal
# and compute the rotation matrix
if (HasValue(variance)) {
v = float(variance);
} else {
v = 50.0;
}
s = identity(dim(a)[3]) * v/sqrt(e[1]);
m = mxm(mxm(e[2:], s), translate(e[2:], x, y));
# We have to reshape the data to push it through mxm() because mxm is dumb.
# Force the input to be bip so we can use the forbidden resize() function.
d = dim(a);
a = bip(a);
resize(a, dim(a)[3], dim(a)[1]*dim(a)[2]);
# subtract the mean of each band, rotate and add back 127
if (HasValue(ignore)) {
avg = avg(a,y,ignore=ignore);
} else {
avg = avg(a,y);
}
a = a - avg;
a = mxm(a, m);
a = a + 127;
a = byte(bip(a));
# reshape the data back to its original size
resize(a, d[1], d[2], d[3]);
a[where mask==1 && a==0]=1
return(a);
}
define do_dcs(ignore,sample,variance) {
if($ARGC==0) {
printf("A nice shortcut to dcs THEMIS images with any three bands\n\n")
printf("Syntax: do_dcs(image,band1,band2,band3,ignore,sample,variance)\n")
printf(" 'ignore' null value of the data, default = -32768\n")
printf(" 'sample' a color sample array\n")
printf(" 'variance' specify color variance, default = 50\n")
printf("Examples:\n")
printf(" display(do_dcs(a,9,6,4))\n")
printf(" display(do_dcs(a,7,5,3,ignore=-32768))\n")
printf(" pic=do_dcs(data,variance=cat(80,50,50,axis=y))\n")
printf(" pic=do_dcs(data,ignore=0,sample=b)\n\n")
printf("if 3 bands aren't specified an 8-7-5 stretch is performed\n\n")
printf("last updated Tue Jan 24 21:34:08 MST 2006\n")
printf("written by Keith Nowicki\n")
return(null)
}
ignorev=-32768
var=cat(50,50,50,axis=y)
#$1 is the data
bands = dim($1)[3]
if(HasValue(variance)) var=variance
B1 = 8
B2 = 7
B3 = 5
if($ARGC == 4) {
B1 = $2
B2 = $3
B3 = $4
}
if(B1 > bands || B2 > bands || B3 > bands) {
printf("specified bands do not match data\n")
return(null)
}
printf("performing a %d-%d-%d stretch\n",B1,B2,B3)
if(HasValue(ignore)) ignorev=ignore
data=$1
if(HasValue(sample)==0) {
image=dcs(cat(data[,,B1],data[,,B2],data[,,B3],axis=z),ignore=ignorev,variance=var)
}
if(HasValue(sample)) {
image=dcs(cat(data[,,B1],data[,,B2],data[,,B3],axis=z),ignore=ignorev,sample=cat(sample[,,B1],sample[,,B2],sample[,,B3],axis=z),variance=var)
}
return(image)
}
define do_rgb(ignore,var) {
if($ARGC==0) {
printf("A nice shortcut to RGB images with any three bands\n\n")
printf("Syntax: do_rgb(image,band1,band2,band3,ignore,var)\n")
printf(" 'ignore' null value of the data, default = -32768\n")
printf(" 'var' variance value used in stretch algorithm, default = 40\n")
printf("Examples:\n")
printf(" display(do_rgb(a,9,6,4))\n")
printf(" display(do_rgb(a,7,5,3,ignore=-32768))\n")
printf("if 3 bands aren't specified an 1-2-3 stretch is performed\n")
printf("Calls sstretch() to do the stretching work\n\n")
printf("last updated Fri Oct 20 10:29:30 MST 2006\n")
printf("written by Keith Nowicki\n")
return(null)
}
ignorev=-32768;
bands = dim($1)[3]
B1 = 1;
B2 = 2;
B3 = 3;
variance = 40
if($ARGC == 4) {
B1 = $2
B2 = $3
B3 = $4
}
if(B1 > bands || B2 > bands || B3 > bands) {
printf("specified bands do not match data\n")
return(null)
}
if(HasValue(ignore)) ignorev=ignore;
if(HasValue(var)) variance=var;
data=$1;
data_pic = sstretch(data, var=variance, ignore=ignorev)
image=cat(data_pic[,,B1], data_pic[,,B2], data_pic[,,B3],axis=z);
return(image);
}
define do_rdcs(sample,ignore,variance,ysize,xsize,verbose) {
if($ARGC==0) {
printf("A nice shortcut to rdcs THEMIS images with any three bands\n\n")
printf("Syntax: do_rdcs(image,band1,band2,band3,ignore,sample,variance)\n")
printf(" 'ignore' null value of the data, default = -32768\n")
printf(" 'sample' a color sample array\n")
printf(" 'variance' specify color variance, default = 50\n")
printf("Examples:\n")
printf(" display(do_rdcs(a,9,6,4))\n")
printf(" display(do_rdcs(a,7,5,3,ignore=-32768))\n")
printf(" pic=do_rdcs(data,variance=cat(80,50,50,axis=y))\n")
printf(" pic=do_rdcs(data,ignore=0,sample=b)\n\n")
printf("if 3 bands aren't specified an 8-7-5 stretch is performed\n\n")
printf("last updated Fri Jun 16 10:39:32 MST 2006\n")
printf("written by Keith Nowicki\n")
return(null)
}
ignorev=-32768
var=cat(50,50,50,axis=y)
ys=1000
xs=4000
verb=0
#$1 is the data
bands = dim($1)[3]
if(HasValue(variance)) var=variance
if(HasValue(xsize)) xs=xsize
if(HasValue(ysize)) ys=ysize
if(HasValue(verbose)) verb=verbose
B1 = 8
B2 = 7
B3 = 5
if($ARGC == 4) {
B1 = $2
B2 = $3
B3 = $4
}
if(B1 > bands || B2 > bands || B3 > bands) {
printf("specified bands do not match data\n")
return(null)
}
printf("performing a %d-%d-%d stretch\n",B1,B2,B3)
if(HasValue(ignore)) ignorev=ignore
data=$1
if(HasValue(sample)==0) {
image=rdcs(cat(data[,,B1],data[,,B2],data[,,B3],axis=z),ignore=ignorev,variance=var,ysize=ys,xsize=xs)
}
if(HasValue(sample)) {
image=rdcs(cat(data[,,B1],data[,,B2],data[,,B3],axis=z),ignore=ignorev,sample=cat(sample[,,B1],sample[,,B2],sample[,,B3],axis=z),variance=var,ysize=ys,xsize=xs)
}
return(image)
}
define do_hsv(ignore,colorize) {
if($ARGC<2) {
printf("Perform an HSV overlay using a color and greyscale image pair\n")
printf("$1 = color image or numeric data to use as Hue and Saturation color\n")
printf("\tNote: if the data is not byte and 3 bands you must choose colorize > 0\n")
printf("$2 = black and white replacement image for Value\n")
printf("colorize = perform the colorize option on $1, same colors as in colorize() (Default is 0)\n")
printf("ignore = ignore value for the colorize option (Default is -32768)\n")
printf("\n")
printf("c.edwards 9/21/10\n")
return(null)
}
if(HasValue(colorize)==0) {
colorize=0
}
if(HasValue(ignore)==0) {
ignore=-32768
}
color=$1
bw=$2
if(dim(color)[3]==1 && colorize==0) {
printf("\bThe input color image is not 3 bands\n")
printf("Please select a colorize option or input a 3 band image\n\n")
return(null)
}
if(sum(dim(color)[:2]==dim(bw)[:2])!=2) {
printf("\n$1 and $2 do not have the same X and Y dimensions\n")
printf("Please make sure input data are of the same size\n\b")
return(null)
}
if(dim(color)[3]==1 && colorize!=0) {
color=colorize(color,ignore=ignore,colors=colorize)
}
if(type(color) != "byte") {
printf("\nSupply $1 as a byte array or a numeric array with colorize set to > 0\n\n")
return(null)
}
if(type(bw) != "byte") {
printf("\nSupply $2 as a byte array\n\n")
return(null)
}
color=rgb2hsv(byte(color))
color[,,3]=bw/255.
color=byte(hsv2rgb(color,maxval=255))
return(color)
}
define sstretch_old(variance,ignore) {
if($ARGC==0) {
printf("sstretch_old - sigma stretch the data\n\n")
printf("Syntax: sstretch_old(image,ignore)\n")
printf(" 'ignore' null value of the data, default = -32768\n")
printf("Examples:\n")
printf(" display(sstretch_old(a))\n")
printf(" pic=(sstretch_old(a,ignore=0))\n")
printf("last updated Tue Jan 24 21:34:08 MST 2006\n")
printf("written by Noel Gorelick\n")
return(null)
}
if (HasValue(ignore) == 0) ignore=-32768;
if (HasValue(variance) == 0) {
variance = 40;
}
if (type($1) == "STRUCT") {
b = avg(($1).data, ignore=ignore, both=1);
c = (($1).data - float(b.avg))*(variance/float(b.stddev))+127;
} else {
b = avg($1, ignore=ignore, both=1);
c = ($1 - float(b.avg))*(variance/float(b.stddev))+127;
}
return(byte(c));
}
define nstretch() {
if($ARGC==0) {
printf("nstretch - normalize stretch the data\n\n")
printf("Syntax: nstretch(image)\n")
printf("Examples:\n")
printf(" display(nstretch(a))\n")
printf(" pic=(nstretch(a)\n")
printf("last updated Tue Jan 24 21:34:08 MST 2006\n")
printf("written by Noel Gorelick\n")
return(null)
}
a = $1
b = moment(a)
c = (((a - b[1]) / (b[2]-b[1])) * 255)
return(byte(c))
}
define stretch(out_low, out_high, full, kluge, ignore,region) {
if ($ARGC == 0) {
printf ("\nImage stretch routine\n\n")
printf ("usage: stretch(array, 'stretch_type', param1, [param2], [out_low = i], [out_high = i], [full=1])\n\n")
printf (" where\n")
printf (" stretch_type = auto\n")
printf (" param1 = auto_end_val\n")
printf (" stretch_type = linear\n")
printf (" param1 = z_min\n")
printf (" param2 = z_max\n")
printf (" stretch_type = sigma\n")
printf (" param1 = # of sigmas away from mean
printf (" stretch_type = minmax\n")
printf (" stretch_type = table\n")
printf (" param1 = table mapping of input values to output values\n")
printf (" [,out_low] is low dn in output image (default = 20)\n")
printf (" [,out_high] is high dn in output image (default = 250)\n")
printf (" [,full=1] returns image, along with min and max input values used in stretch in structure\n\n")
printf (" example: out = stretch(in, 'sigma', 3) \n")
printf (" out = stretch(in, 'auto', .01, out_low=10) \n")
printf (" out = stretch(in, 'sigma', 3 out_low=10, full=1) \n")
printf (" out = stretch(in, 'minmax') \n")
printf (" out = stretch(in, 'table', table)\n\n")
printf ("last updated Tue Jan 24 21:34:08 MST 2006\n")
printf ("written by Phil Christensen\n")
return(null)
}
# p. christensen 2/6/00
# k. nowicki 01/17/06
input = $1;
b = byte(input*0.);
size = dim(input);
x1 = 1;
x2 = size[1];
y1 = 1;
y2 = size[2];
if(HasValue(region)) {
x1 = region[1];
x2 = region[2];
y1 = region[3];
y2 = region[4];
printf("statistics region = %d %d %d %d\n", x1, x2, y1, y2);
}
# set low and high values in output image for cosmetic purpose
low_out = 10;
high_out = 250;
if ($2 == "minmax") {
if(HasValue(ignore)) {
m = moment(input, ignore=0);
} else {
m = moment(input);
}
b = (((input - m[1]) / (m[2]-m[1])) * 255);
printf("stretching %.1e to 0 %.1e to 255\n", m[1], m[2]);
return(byte(b));
}
if ($2 == "auto") {
# do this one the "dumb" way - loop on bands
nbands = dim(input)[3];
low_dn_array = create(1,1,dim(input)[3],format=float);
high_dn_array = create(1,1,dim(input)[3],format=float);
for (iband = 1 ; iband <= nbands ; iband+=1) {
# set values less than 0 to 0
# tmp = input[,,iband]
# tmp[where(input[,,iband]<=0.)] = 0.
# a = histogram(input[,,iband], steps=2560, compress=1)
# changed to include region - prc - 10/30/02
a = histogram(input[x1:x2, y1:y2, iband], steps=2560, compress=1);
len = dim(a)[2];
# zero out 0 and 255 values - i.e. remove peaks at start and end
a[2,1] = 0;
a[2,len] = 0;
if(HasValue(kluge)) {
printf(" doing the kluge\n");
# zero out 20 points at high and low ends of histogram before to auto ends
a[2,1:100] = 0;
a[2,len-100:len] = 0;
}
total = sum(a[2,]);
auto_end = $3;
end = total * $3;
if(HasValue(out_low)) {
low_out = out_low;
}
if(HasValue(out_high)) {
high_out = out_high;
}
# this code is for historical reasons - param_map uses this form
if ($ARGC >= 4) {
low_out = $4;
high_out = $5;
}
low_dn = low_out;
high_dn = high_out;
# find auto-end point on low end of histogram
cum = 0;
for (i = 1 ; i <= len ; i+=1) {
cum = cum + a[2,i, 1];
if (cum > end) {
# found point at low DN limit
low_dn = a[1,i, 1];
break;
}
}
# find auto-end point on high end of histogram
cum = 0;
for (i = 1 ; i <= len ; i+=1) {
counter = len+1-i;
# counter = 255-i
cum = cum + a[2,counter, 1];
if (cum > end) {
# found point at high DN limit
high_dn = a[1,counter, 1];
break;
}
}
printf("doing auto_end stretch: plane = %d auto_ends = %.4e low input value = %.2e high input value = %.2e low output dn = %d high output dn = %d\n",iband, auto_end, low_dn, high_dn, low_out, high_out);
slope = float(high_out - low_out) / (high_dn - low_dn);
printf ("slope = %.2f\n",slope);
b[,,iband] = (input[,,iband] - low_dn) * slope + low_out;
low_dn_array[,,iband] = low_dn;
high_dn_array[,,iband] = high_dn;
}
}
if ($2 == "sigma") {
printf("doing sigma stretch\n")
if(HasValue(ignore)) {
# changed to use region - prc - 10/30/02
tmp = avg(input[x1:x2, y1:y2], axis=xy, both=1, ignore=0);
} else {
tmp = avg(input[x1:x2, y1:y2], axis=xy, both=1);
}
ave = tmp.avg;
stddev = tmp.stddev;
if ($ARGC >= 4) {
low_out = $4;
high_out = $5;
}
low_dn = ave - (stddev * $3);
high_dn = ave + (stddev * $3);
slope = float(high_out - low_out) / (high_dn - low_dn);
b = (input - low_dn) * slope + low_out;
# printf("linear stretch: iband = %d low_dn = %.2f high_dn = %.2f slope = %.2f\n",iband, low_dn, high_dn, slope)
}
if ($2 == "linear") {
printf("doing linear stretch\n");
low_dn = $3;
high_dn = $4;
if(HasValue(out_low)) {
low_out = out_low;
}
if(HasValue(out_high)) {
high_out = out_high;
}
slope = float(high_out - low_out) / (high_dn - low_dn);
b = (input - low_dn) * slope + low_out;
printf("linear stretch: low_out = %d high_out = %d low_dn = %.2f high_dn = %.2f slope = %.2f\n", low_out, high_out, low_dn, high_dn, slope);
}
if ($2 == "table") {
printf("doing table mapping stretch\n");
table = $3;
tdimx = dim(table)[1];
tdimy = dim(table)[2];
if(tdimy <= 1 || tdimx != 2) {
printf("\nThe table must a 2xNx1 array, where N must be larger than 1\n");
printf("column 1 is the input value, column 2 is the output value\n");
return(NULL);
}
b = input*0.0;
#loop through the mapping sets and perform linear stretches
for(i=1;i
low_in = float(table[1,i]);
high_in = float(table[1,i+1]);
low_out = float(table[2,i]);
high_out = float(table[2,i+1]);
slope = float(high_out - low_out) / (high_in - low_in);
b[where (input >= low_in && input < high_in)] = (input - low_in) * slope + low_out;
}
}
if(HasValue(full)) {
# return array and min, max values
struc = struct();
struc.out = b;
if($2 == "auto") {
struc.min = low_dn_array;
struc.max = high_dn_array;
} else {
struc.min = low_dn;
struc.max = high_dn;
}
return(struc);
} else {
return(byte(b));
}
}
define cleandcs(opt) {
if($ARGC==0) {
printf("cleandcs() - removes unwanted saturated pixels in dcs images\n")
printf("Syntax: cleandcs(image,opt)\n")
printf(" 'image' a byte color picture resulting from a dcs operation\n")
printf(" 'opt' option for level of cleaning. Default = 2\n")
printf(" opt=1 removes all pixels with 255 in one band and 0 in others\n")
printf(" opt=2 removes all pixels with 255 in two bands and 0 in other\n\n")
printf("last updated Tue Jan 24 21:34:08 MST 2006\n")
printf("written by Keith Nowicki\n")
return(null)
}
op = 2
img = $1
if(HasValue(opt)) op = int(opt)
if(op == 1) img[where max(img,axis=z) == 255 && max(img,axis=z,ignore=255) == 0] = 0
if(op == 2) img[where min(img,axis=z) == 0 && min(img,axis=z,ignore=0) == 255] = 0
return(img)
}
define rdcs(data,sample,ignore,variance,ysize,xsize,verbose) {
#c.edwards modified 6-24-10 to make the minimum stretched data value 1 and not 0
if($ARGC==0 && HasValue(data)==0) {
printf("\nCompute a Running DCS Stretch\n")
printf("Calls dcs() as stretching function\n\n")
printf("$1=data (3-band or 1-band)\n")
printf("data=data option for more efficient memory usage\n")
printf("ignore=value to ignore (Default=-32768)\n")
printf("ysize=size of the chunks to stretch in y (Default=1000)\n")
printf("xsize=size of the chunks to stretch in x (Default=4000)\n")
printf("NOTE: both xsize and ysize must be even integers\n")
printf("variance=variance of the stretch (Default=50)\n")
printf("sample=sample image\n\n")
printf("NOTE: can also be used to change variance of each band\n")
printf("variance=cat(50,20,40,axis=y)\n\n")
printf("verbose = turn on print statements\n\n")
printf("c.edwards 3/7/06\n\n")
return(null)
}
#initialize input
if(HasValue(data)==0) data=$1
xs=4000
ys=1000
ign=-32768
var=50.0
verb=0
ydiff=xdiff=0
dim=dim(data)
if(HasValue(verbose)) verb=1
if(HasValue(ysize)) ys=int(ysize)
if(HasValue(xsize)) xs=int(xsize)
if(HasValue(ignore)) ign=ignore
if(HasValue(variance)) var=variance
if(ys%2 == 1) ys-=1
if(xs%2 == 1) xs-=1
#create output pic, ramps and indices
pic=byte(data*0)
mask=byte(data[,,1]*0+1)
mask[where sum(data == ignore,axis=z)>=1]=0
yrange=xrange=clone(0,2,1,1)
yramp=create(1,int(ys/2),1,start=0,step=(1/(ys/2.)),format=float)
xramp=create(int(xs/2),1,1,start=0,step=(1/(xs/2.)),format=float)
#set maximum loop values
maxx=dim[1]-int(xs/2)
maxy=dim[2]-int(ys/2)
if(dim[1]
if(dim[2]
#loop through x to get x-direction
for(i=1;i<=maxx;i+=int(xs/2)) {
xrange[1]=i;xrange[2]=i+xs-1
if(xrange[2]>=dim[1]) {
xdiff=xrange[2]-dim[1]
xrange[2]=dim[1]
}
tmppic=byte(clone(0,xrange[2]-i+1,dim[2],dim[3]))
#loop through y to get y-direction
for(j=1;j<=maxy;j+=int(ys/2)) {
yrange[1]=j;yrange[2]=j+ys-1
if(yrange[2]>=dim[2]) {
ydiff=yrange[2]-dim[2]
yrange[2]=dim[2]
}
if(verb==1) printf("x=%i:%i y=%i:%i\n",xrange[1],xrange[2],yrange[1],yrange[2])
#stretch the data
if(HasValue(sample)) {
tmp=dcs(data[xrange[1]:xrange[2],yrange[1]:yrange[2],],ignore=ign,variance=var,sample=sample[xrange[1]:xrange[2],yrange[1]:yrange[2],])
}else{
tmp=dcs(data[xrange[1]:xrange[2],yrange[1]:yrange[2],],ignore=ign,variance=var)
}
#if its the first y-chunk then insert else ramp it into the y-column
if(j==1) tmppic[,yrange[1]:yrange[2],]=tmp
if((j!=1 && ydiff>int(ys/2))) break;
if(j!=1) {
tmppic[,yrange[1]:yrange[2]-int(ys/2)+ydiff]=tmppic[,yrange[1]:yrange[2]-int(ys/2)+ydiff]*(1-yramp)
tmp[,:int(ys/2),]=tmp[,:int(ys/2),]*yramp
tmppic[,yrange[1]:yrange[2],]=tmppic[,yrange[1]:yrange[2],]+tmp
}
ydiff=0
}
#if its the first x-chunk then insert else ramp it into the full picutre
if(i==1) pic[xrange[1]:xrange[2],]=tmppic
if((i!=1 && xdiff>int(xs/2))) return(bip(byte(pic)))
if(i!=1) {
pic[xrange[1]:xrange[2]-int(xs/2)+xdiff,]=pic[xrange[1]:xrange[2]-int(xs/2)+xdiff,]*(1-xramp)
tmppic[:int(xs/2),]=tmppic[:int(xs/2),]*xramp
pic[xrange[1]:xrange[2],]=pic[xrange[1]:xrange[2],]+tmppic
}
}
pic[where mask==1 && pic==0]=1
return(bip(byte(pic)))
}
define rsstretch(data,ignore,variance,ysize,xsize,verbose) {
if($ARGC==0 && HasValue(data)==0) {
printf("\nCompute a Running SStretch\n")
printf("Calls sstretch() as stretching function\n\n")
printf("$1=data\n")
printf("data=data option for more efficient memory usage\n")
printf("ignore=value to ignore (Default=-32768)\n")
printf("ysize=size of the chunks to stretch in y (Default=1000)\n")
printf("xsize=size of the chunks to stretch in x (Default=4000)\n")
printf("NOTE: both xsize and ysize must be even integers\n")
printf("variance=variance of the stretch (Default=40)\n\n")
printf("verbose = turn on print statements\n\n")
printf("c.edwards 3/7/06\n\n")
return(null)
}
#initialize input
if(HasValue(data)==0) data=$1
xs=4000
ys=1000
ign=-32768
var=40
verb=0
ydiff=xdiff=0
dim=dim(data)
if(HasValue(verbose)) verb=1
if(HasValue(ysize)) ys=int(ysize)
if(HasValue(xsize)) xs=int(xsize)
if(HasValue(ignore)) ign=ignore
if(HasValue(variance)) var=variance
if(ys%2 == 1) ys-=1
if(xs%2 == 1) xs-=1
#create output pic, ramps and indices
pic=byte(data)
pic[]=0
yrange=xrange=clone(0,2,1,1)
yramp=create(1,int(ys/2),1,start=0,step=(1/(ys/2.)),format=float)
xramp=create(int(xs/2),1,1,start=0,step=(1/(xs/2.)),format=float)
#set maximum loop values
maxx=dim[1]-int(xs/2)
maxy=dim[2]-int(ys/2)
if(dim[1]
if(dim[2]
#loop through x to get x-direction
for(i=1;i<=maxx;i+=int(xs/2)) {
xrange[1]=i;xrange[2]=i+xs-1
if(xrange[2]>=dim[1]) {
xdiff=xrange[2]-dim[1]
xrange[2]=dim[1]
}
tmppic=byte(clone(byte(0),xrange[2]-i+1,dim[2],dim[3]))
#loop through y to get y-direction
for(j=1;j<=maxy;j+=int(ys/2)) {
yrange[1]=j;yrange[2]=j+ys-1
if(yrange[2]>=dim[2]) {
ydiff=yrange[2]-dim[2]
yrange[2]=dim[2]
}
if(verb==1) printf("x=%i:%i y=%i:%i\n",xrange[1],xrange[2],yrange[1],yrange[2])
#stretch the data
tmp=sstretch(data[xrange[1]:xrange[2],yrange[1]:yrange[2],],ignore=ign,variance=var)
#if its the first y-chunk then insert else ramp it into the y-column
if(j==1) tmppic[,yrange[1]:yrange[2],]=tmp
if((j!=1 && ydiff>int(ys/2))) break;
if(j!=1) {
tmppic[,yrange[1]:yrange[2]-int(ys/2)+ydiff]=tmppic[,yrange[1]:yrange[2]-int(ys/2)+ydiff]*(1-yramp)
tmp[,:int(ys/2),]=tmp[,:int(ys/2),]*yramp
tmppic[,yrange[1]:yrange[2],]=tmppic[,yrange[1]:yrange[2],]+tmp
}
ydiff=0
}
#if its the first x-chunk then insert else ramp it into the full picutre
if(i==1) pic[xrange[1]:xrange[2],]=tmppic
if((i!=1 && xdiff>int(xs/2))) {
pic[where data!=ign && pic==0]=1
return(bip(byte(pic)))
}
if(i!=1) {
pic[xrange[1]:xrange[2]-int(xs/2)+xdiff,]=pic[xrange[1]:xrange[2]-int(xs/2)+xdiff,]*(1-xramp)
tmppic[:int(xs/2),]=tmppic[:int(xs/2),]*xramp
pic[xrange[1]:xrange[2],]=pic[xrange[1]:xrange[2],]+tmppic
}
}
pic[where data!=ign && pic==0]=1
return(bip(byte(pic)))
}
define rmnoise_pca(null,filt,noise_sample,debug,b10,mult) {
if($ARGC==0) {
printf("\nRemove image noise with cse.rmnoise() and other methods\n\n")
printf("This script uses cse.rmnoise() and eigen vectors\n")
printf("to help clean the data. Compared to\n")
printf("cse.rmnoise() more color information is retained\n\n")
printf("$1 = 10 band THEMIS IR data (will accept less bands but specify b10 location)\n")
printf("null = null value of the image (Default is -32768)\n")
printf("filt = the filter length for cse.rmnoise() (Default is 7)\n")
printf("b10 = the loacation of band 10 set to 0 if no band 10 exists\n\n")
printf("mult = multiplier (Default=10000, best for radiance data)\n")
printf("noise_sample = the sample area for eigen2 to determine vectors\n")
printf("This is found automatically but to use this to submit \n")
printf("an alternate sample of the noise returned from thm.white_noise_remove2()\n")
printf("This should be fewer than 4000 total points or eigen2 may fail\n")
printf("The number of points is related to the time required to process\n")
printf("\nNOTE: This WILL CHANGE the data and should only be used for\n")
printf("\"pretty pictures\". Band 10 will be smoothed by convolution\n")
printf("\nc.edwards 3/7/06\n\n")
return(null)
}
verbose=0
data=$1
mul=10000
if(HasValue(mult)) mul=mult
# get arguement and data properties
filter=7
if(HasValue(filt)) {
filter=filt
}
nullval=-32768
if(HasValue(null)) {
nullval=null
}
x=dim(data)[1]
y=dim(data)[2]
z=dim(data)[3]
band10=z
if(HasValue(b10)) band10=b10
if(band10==0) band10=z+1
#no user defined sample
if(HasValue(noise_sample)==0) {
#determine sample size for x and y
count=1
xsample=x
while(xsample>10) {
xsample=x/count
count+=1
}
xsamplen=count
count=1
ysample=y
while(ysample>40) {
ysample=y/count
count+=1
}
ysamplen=count
}
#run rmnoise on the data
clean=thm.white_noise_remove2(data,ignore=nullval,filt=filter,b10=band10)
#separate the old noise array
noise=data-clean
#sample the data
if(HasValue(noise_sample)==0) {
noise_mult_sample=noise[::xsamplen,::ysamplen,1:band10-1]*float(mul)
}
#user specified sample
if(HasValue(noise_sample)) {
noise_mult_sample=noise_sample
}
#run eigen2 on the sample data
printf("\nDetermining Eigen Vectors\n")
e_struct=eigen2(noise_mult_sample,band10-2)
#run jpca on data and eigen_struct to apply vectors to data
printf("Calculating Principal Components\n")
noise_pca=pca2(noise[,,:band10-1]*10000,e_struct)
#set the first principal component equal to zero
noise_pca[,,1]=noise_pca[,,1]*0
#reconstruct the noise (error message of division by zero normal)
noise_recon=pca2_recon(noise_pca,e_struct,band10-2)
noise_recon=noise_recon/float(mul)
noise_pca=0
#add the correction back
new_clean=clean
clean=0
new_clean[,,:band10-1]=new_clean[,,:band10-1]+(noise[,,:band10-1]-noise_recon)
noise=noise_recon=0
#take care of color shift caused by rmnoise
color_shift=avg(data-new_clean,axis=xy,ignore=nullval)
new_clean[,,:band10-1]=new_clean[,,:band10-1]+color_shift[,,:band10-1]
verbose=3
#if there is debug option selected then return extra info
if(HasValue(debug)) {
out={}
out.new_clean=float(new_clean)
out.eigen_struct=e_struct
out.color_shift=color_shift
return(out)
} else {
return(float(new_clean))
}
}
define rrmnoise_pca(data,ignore,ysize,xsize,verbose,b10,mult) {
if($ARGC==0 && HasValue(data)==0) {
printf("\nCompute a Running rmnoise_pca\n")
printf("$1=data\n")
printf("data=data option for more efficient memory usage\n")
printf("ignore=value to ignore (Default=-32768)\n")
printf("ysize=size of the chunks (Default=1000)\n")
printf("xsize=size of the chunks (Default=1000)\n")
printf("b10=band 10 location (Default=10, 0=no band 10)\n")
printf("mult=multiplier value (Default=10000, good for radiance data)\n")
printf("CAUTION: SLOWER THAN HELL!!!\n")
printf("verbose = turn on print statements\n\n")
printf("c.edwards 3/7/06\n\n")
return(null)
}
#initialize input
if(HasValue(data)==0) data=$1
xs=1000
ys=1000
ign=-32768
var=40
verb=0
ydiff=xdiff=0
mul=10000
dim=dim(data)
band10=dim[3]
if(HasValue(verbose)) verb=1
if(HasValue(ysize)) ys=int(ysize)
if(HasValue(xsize)) xs=int(xsize)
if(HasValue(ignore)) ign=ignore
if(HasValue(variance)) var=variance
if(HasValue(b10)) band10=b10
if(HasValue(mult)) mul=mult
if(ys%2 == 1) ys-=1
if(xs%2 == 1) xs-=1
#create output pic, ramps and indices
pic=float(data*0)
bmask=byte(pic*0)
bmask[where data!=ign]=1
yrange=xrange=clone(0,2,1,1)
yramp=create(1,int(ys/2),1,start=0,step=(1/(ys/2.)),format=float)
xramp=create(int(xs/2),1,1,start=0,step=(1/(xs/2.)),format=float)
#set maximum loop values
maxx=dim[1]-int(xs/2)
maxy=dim[2]-int(ys/2)
if(dim[1]
if(dim[2]
#loop through x to get x-direction
for(i=1;i<=maxx;i+=int(xs/2)) {
xrange[1]=i;xrange[2]=i+xs-1
if(xrange[2]>=dim[1]) {
xdiff=xrange[2]-dim[1]
xrange[2]=dim[1]
}
tmppic=float(clone(0,xrange[2]-i+1,dim[2],dim[3]))
#loop through y to get y-direction
for(j=1;j<=maxy;j+=int(ys/2)) {
yrange[1]=j;yrange[2]=j+ys-1
if(yrange[2]>=dim[2]) {
ydiff=yrange[2]-dim[2]
yrange[2]=dim[2]
}
if(verb==1) printf("x=%i:%i y=%i:%i\n",xrange[1],xrange[2],yrange[1],yrange[2])
#stretch the data
tmp=rmnoise_pca(data[xrange[1]:xrange[2],yrange[1]:yrange[2],],null=ign,mult=mul,b10=band10)
#if its the first y-chunk then insert else ramp it into the y-column
if(j==1) tmppic[,yrange[1]:yrange[2],]=tmp
if((j!=1 && ydiff>int(ys/2))) break;
if(j!=1) {
tmppic[,yrange[1]:yrange[2]-int(ys/2)+ydiff]=tmppic[,yrange[1]:yrange[2]-int(ys/2)+ydiff]*(1-yramp)
tmp[,:int(ys/2),]=tmp[,:int(ys/2),]*yramp
tmppic[,yrange[1]:yrange[2],]=tmppic[,yrange[1]:yrange[2],]+tmp
}
ydiff=0
}
#if its the first x-chunk then insert else ramp it into the full picutre
if(i==1) pic[xrange[1]:xrange[2],]=tmppic
if((i!=1 && xdiff>int(xs/2))) {
pic[where bmask==0]=ign
return(float(pic))
}
if(i!=1) {
pic[xrange[1]:xrange[2]-int(xs/2)+xdiff,]=pic[xrange[1]:xrange[2]-int(xs/2)+xdiff,]*(1-xramp)
tmppic[:int(xs/2),]=tmppic[:int(xs/2),]*xramp
pic[xrange[1]:xrange[2],]=pic[xrange[1]:xrange[2],]+tmppic
}
}
pic[where bmask==0]=ign
return(float(pic))
}
define colorize(colors,values,type,ignore) {
#Added $DV_EX support
#Added $DV_SCRIPT_FILES support - 5-13-08
#added color_tables directory to $DV_SCRIPT_FILES
if($ARGC==0) {
printf("Colorize any numeric array\n\n")
printf("$1 = numeric array to be colorized\n")
printf("ignore = ignore value of the data\n")
printf("type = rgb or hsv colorization (Default = \"rgb\")\n")
printf("colors = number of colors to use 1 = 6 colors, 2 = 7 colors\n")
printf(" 3 = xv color scale, 4 = tes colormap, 5 = tes daily colormap (Default = 1)\n")
printf("values = values of the data for each color (Assigned linearly by default)\n")
printf("\nn.gorelick\n")
printf("Modified by c.edwards Aug 15, 2006\n")
return(null)
}
in = $1
if (HasValue(ignore) == 0) {
ignore = -32768
}
color=1
if(HasValue(colors)) color=colors; colors=0
if(color==1) {
colors = {}
colors += { get_color("Magenta")}
colors += { get_color("Blue")}
colors += { get_color("Lime")}
colors += { get_color("Yellow")}
colors += { get_color("Orange")}
colors += { get_color("Red") }
}
if (color==2) {
colors = {}
colors += { get_color("Magenta")}
colors += { get_color("Blue")}
colors += { get_color("Cyan")}
colors += { get_color("Lime")}
colors += { get_color("Yellow")}
colors += { get_color("Orange")}
colors += { get_color("Red") }
}
if(color==3) {
scalebar=read($DV_SCRIPT_FILES+"/color_tables/xvscalebar.ppm")
colors={}
for(i=1;i<=dim(scalebar)[1];i+=1) {
colors+={scalebar[i,1]}
}
values=create(dim(scalebar)[1],1,1,bsq,int,0,1)
}
if(color==4) {
scalebar=read($DV_SCRIPT_FILES+"/color_tables/colormap.ppm")
colors={}
for(i=1;i<=dim(scalebar)[1];i+=1) {
colors+={scalebar[i,1]}
}
values=create(dim(scalebar)[1],1,1,bsq,int,0,1)
}
if(color==5) {
scalebar=read($DV_SCRIPT_FILES+"/color_tables/colormap_daily.ppm")
colors={}
for(i=1;i<=dim(scalebar)[1];i+=1) {
colors+={scalebar[i,1]}
}
values=create(dim(scalebar)[1],1,1,bsq,int,0,1)
}
if (HasValue(values) == 0) {
values = create(length(colors),1,1)/float(length(colors)-1)
maxval = max(in,ignore=ignore)
minval = min(in,ignore=ignore)
values = values * (maxval-minval) + minval
for (i = 1 ; i <= length(colors) ; i+=1) {
printf("%3.3d,%3.3d,%3.3d = %f\n", colors[i][,,1], colors[i][,,2], colors[i][,,3], values[i]);
}
}
if (HasValue(type) == 0) type = "rgb"
c = colors[1]
for (i = 2 ; i <= length(colors) ; i+=1) {
c = cat(c, colors[i], x)
}
if (type == "hsv") c = rgb2hsv(c,maxval=255)
rgb = clone(float(in)*0, 1, 1, 3)
rgb[,,1] = interp(c[,,1], values, in)
rgb[,,2] = interp(c[,,2], values, in)
rgb[,,3] = interp(c[,,3], values, in)
rgb[where in < values[1]] = c[1]
rgb[where in > values[length(values)]] = c[length(colors)]
rgb[where in==ignore] = 0
if (type == "hsv") rgb = hsv2rgb(rgb, maxval=255)
rgb = bip(byte(rgb))
return(rgb)
}
define get_color(display) {
colors = {"AliceBlue", 240//248//255, \
"AntiqueWhite", 250//235//215, \
"Aqua", 0//255//255, \
"Aquamarine", 127//255//212, \
"Azure", 240//255//255, \
"Beige", 245//245//220, \
"Bisque", 255//228//196, \
"Black", 0//0//0, \
"BlanchedAlmond", 255//235//205, \
"Blue", 0//0//255, \
"BlueViolet", 138//43//226, \
"Brown", 165//42//42, \
"Burlywood", 222//184//135, \
"CadetBlue", 95//158//160, \
"Chartreuse", 127//255//0, \
"Chocolate", 210//105//30, \
"Coral", 255//127//80, \
"CornflowerBlue", 100//149//237, \
"Cornsilk", 255//248//220, \
"Cyan", 0//255//255, \
"DarkBlue", 0//0//139, \
"DarkCyan", 0//139//139, \
"DarkGoldenrod", 184//134//11, \
"DarkGray", 169//169//169, \
"DarkGreen", 0//100//0, \
"DarkKhaki", 189//183//107, \
"DarkMagenta", 139//0//139, \
"DarkOliveGreen", 85//107//47, \
"DarkOrange", 255//140//0, \
"DarkOrchid", 153//50//204, \
"DarkRed", 139//0//0, \
"DarkSalmon", 233//150//122, \
"DarkSeaGreen", 143//188//143, \
"DarkSlateBlue", 72//61//139, \
"DarkSlateGray", 47//79//79, \
"DarkTurquoise", 0//206//209, \
"DarkViolet", 148//0//211, \
"DeepPink", 255//20//147, \
"DeepSkyBlue", 0//191//255, \
"DimGray", 105//105//105, \
"DodgerBlue", 30//144//255, \
"Firebrick", 178//34//34, \
"FloralWhite", 255//250//240, \
"ForestGreen", 34//139//34, \
"Fuschia", 255//0//255, \
"Gainsboro", 220//220//220, \
"GhostWhite", 255//250//250, \
"Gold", 255//215//0, \
"Goldenrod", 218//165//32, \
"Gray", 128//128//128, \
"Green", 0//128//0, \
"GreenYellow", 173//255//47, \
"Honeydew", 240//255//240, \
"HotPink", 255//105//180, \
"Ivory", 255//255//240, \
"Khaki", 240//230//140, \
"Lavender", 230//230//250, \
"LavenderBlush", 255//240//245, \
"LawnGreen", 124//252//0, \
"LemonChiffon", 255//250//205, \
"LightBlue", 173//216//230, \
"LightCoral", 240//128//128, \
"LightCyan", 224//255//255, \
"LightGoldenrod", 238//221//130, \
"LightGoldenrodYellow", 250//250//210, \
"LightGray", 211//211//211, \
"LightGreen", 144//238//144, \
"LightPink", 255//182//193, \
"LightSalmon", 255//160//122, \
"LightSeaGreen", 32//178//170, \
"LightSkyBlue", 135//206//250, \
"LightSlateBlue", 132//112//255, \
"LightSlateGray", 119//136//153, \
"LightSteelBlue", 176//196//222, \
"LightYellow", 255//255//224, \
"Lime", 0//255//0, \
"LimeGreen", 50//205//50, \
"Linen", 250//240//230, \
"Magenta", 255//0//255, \
"Maroon", 128//0//0, \
"MediumAquamarine", 102//205//170, \
"MediumBlue", 0//0//205, \
"MediumOrchid", 186//85//211, \
"MediumPurple", 147//112//219, \
"MediumSeaGreen", 60//179//113, \
"MediumSlateBlue", 123//104//238, \
"MediumSpringGreen", 0//250//154, \
"MediumTurquoise", 72//209//204, \
"MediumVioletRed", 199//21//133, \
"MidnightBlue", 25//25//112, \
"MintCream", 245//255//250, \
"MistyRose", 255//228//225, \
"Moccasin", 255//228//181, \
"NavajoWhite", 255//222//173 }
colors+={"Navy", 0//0//128}
colors+={"OldLace", 253//245//230, \
"Olive", 128//128//0, \
"OliveDrab", 107//142//35, \
"Orange", 255//128//0, \
"OrangeRed", 255//69//0, \
"Orchid", 218//112//214, \
"PaleGoldenrod", 238//232//170, \
"PaleGreen", 152//251//152, \
"PaleTurquoise", 175//238//238, \
"PaleVioletRed", 219//112//147, \
"PapayaWhip", 255//239//213, \
"PeachPuff", 255//218//185, \
"Peru", 205//133//63, \
"Pink", 255//192//203, \
"Plum", 221//160//221, \
"PowderBlue", 176//224//230, \
"Purple", 128//0//128, \
"Red", 255//0//0, \
"RosyBrown", 188//143//143, \
"RoyalBlue", 65//105//225, \
"SaddleBrown", 139//69//19, \
"Salmon", 250//128//114, \
"SandyBrown", 244//164//96, \
"SeaGreen", 46//139//87, \
"Seashell", 255//245//238, \
"Sienna", 160//82//45, \
"Silver", 192//192//192, \
"SkyBlue", 135//206//235, \
"SlateBlue", 106//90//205, \
"SlateGray", 112//128//144, \
"Snow", 255//250//250, \
"SpringGreen", 0//255//127, \
"SteelBlue", 70//130//180, \
"Tan", 210//180//140, \
"Teal", 0//128//128, \
"Thistle", 216//191//216, \
"Tomato", 255//99//71, \
"Turquoise", 64//224//208, \
"Violet", 238//130//238, \
"VioletRed", 208//32//144, \
"Wheat", 245//222//179, \
"White", 255//255//255, \
"WhiteSmoke", 245//245//245, \
"Yellow", 255//255//0, \
"YellowGreen", 154//205//50}
if ($ARGC == 0) {
printf("\nGets the rgb color values from a color name\n")
printf("The colors included are many (if not all) of the standard web colors\n\n")
printf("n.gorelick\n")
printf("Modified by c.edwards Aug 15, 2006\n\n")
printf("Color options are shown below\n\n")
for (i = 1 ; i < length(colors) ; i+=2) {
printf("%s\n", colors[i]);
}
} else {
for (i = 1 ; i < length(colors) ; i+=2) {
if (colors[i] == $1) {
if(HasValue(display)) {
display(bip(byte(clone(translate(colors[i+1],x,z),100,100))))
}
return(translate(colors[i+1],x,z))
}
}
printf("Color not found: %s\n", $1)
return(null)
}
}
define lstretch(ignore,in_low,in_high,out_low,out_high) {
if ($ARGC == 0) {
printf("Do a linear stretch for a single band image.\n")
printf(" \n")
printf("$1 = input image\n")
printf("in_low = in data low value\n")
printf("in_high = in data high value\n")
printf("out_low = out data low value\n")
printf("out_high = out data high value\n")
printf("All ins and outs are optional. Output is byte.\n")
printf("ignore = ignore value (default is -32768)\n")
printf(" \n")
printf(" \n")
return(0)
}
data=$1
# Make sure we only have a 1 band image
data=data[,,1]
if (HasValue(ignore) == 0) {
ignorev=-32768
} else {
ignorev=ignore
}
if (HasValue(in_low) == 0) {
in_lowv=min(data, ignore=ignorev)
} else {
in_lowv=float(in_low)
}
if (HasValue(in_high) == 0) {
in_highv=max(data, ignore=ignorev)
} else {
in_highv=float(in_high)
}
if (HasValue(out_low) == 0) {
out_lowv=float(0)
} else {
out_lowv=float(out_low)
}
if (HasValue(out_high) == 0) {
out_highv=float(255)
} else {
out_highv=float(out_high)
}
slope=float((out_highv-out_lowv))/float((in_highv-in_lowv))
intercept=out_lowv-in_lowv*slope
out=slope*data+intercept
out[where data!=ignorev && out<=out_lowv]=out_lowv
out[where (data==ignorev)]=0
return(byte(out))
}