#Original Function Location: /themis/lib/dav_lib/library/dshadow.dvrc Source Code for Function: "kt2()" in "dshadow.dvrc" (Beta)

define dcss(ignore,sample,variance) {
if (\$ARGC != 1) {
printf("\ndcs()\n")
printf("performs a decorrelation stretch of data.\n")
printf("data will be centered around 0 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 = -32768.\n")
printf("  'sample'   is an array of identical size to the input\n")
printf("             used to set 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 = 12850 \n\n")
printf("last updated Sun Jan 31 17:48:00 PST 2010\n")
printf("based on dcs() by Noel Gorelick - modified for 16-bit use by Ben Lincoln.\n")
return(null)
}

#
# this does a decorrelation stretch using 16-bit data
# Data will be centered around 0, with a sigma of +/- 12850.
#

# compute the symmetric covariance matrix and eigenvectors of it.
a = \$1;

if (HasValue(sample)) {
if (HasValue(ignore)) {
c = covar(sample, axis=z, ignore=ignore);
} else {
c = covar(sample, axis=z, ignore=-32768);
}
} else {
if (HasValue(ignore)) {
c = covar(a, axis=z, ignore=ignore);
} else {
c = covar(a, axis=z, ignore=-32768);
}
}
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 = 12850.0;
}

s = identity(dim(a)) * v/sqrt(e);
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), dim(a)*dim(a));

# subtract the mean of each band and rotate
if (HasValue(ignore)) {
avg = avg(a,y,ignore=ignore);
} else {
avg = avg(a,y,ignore=-32768);
}
a = a - avg;
a = mxm(a, m);
a = short(bip(a));

# reshape the data back to its original size
resize(a, d, d, d);
return(a);
}

define rdcss(data,sample,ignore,variance,ysizepercent,xsizepercent,verbose) {

if(\$ARGC==0 && HasValue(data)==0) {
printf("\nCompute a Running DCS Stretch\n")
printf("Calls dcs16() 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("ysizepercent=size of the chunks to stretch in y (Default=1/4 of the X dimension of the data)\n")
printf("xsizepercent=size of the chunks to stretch in x (Default=1/4 of the Y dimension of the data)\n")
printf("NOTE: both xsize and ysize must be even integers\n")
printf("variance=variance of the stretch (Default=12850)\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("last updated Sun Jan 31 17:55:00 PST 2010\n")
printf("based on rdcs() by c.edwards 3/7/06 - modified for 16-bit use by Ben Lincoln.\n\n")
return(null)
}

#initialize input
if(HasValue(data)==0) data=\$1
ign=-32768
var=12850.0
verb=0
ydiff=xdiff=0
dim=dim(data)

xs=int(dim[1,0,0] * 0.25)
ys=int(dim[2,0,0] * 0.25)

if(HasValue(verbose)) verb=1
if(HasValue(ysizepercent)) ys=int(dim[1,0,0] * ysizepercent)
if(HasValue(xsizepercent)) xs=int(dim[2,0,0] * xsizepercent)
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=short(data*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-int(xs/2)
maxy=dim-int(ys/2)
if(dim   if(dim
#loop through x to get x-direction
for(i=1;i<=maxx;i+=int(xs/2)) {

xrange=i;xrange=i+xs-1
if(xrange>=dim) {
xdiff=xrange-dim
xrange=dim
}

tmppic=short(clone(0,xrange-i+1,dim,dim))

#loop through y to get y-direction
for(j=1;j<=maxy;j+=int(ys/2)) {

yrange=j;yrange=j+ys-1
if(yrange>=dim) {
ydiff=yrange-dim
yrange=dim
}

if(verb==1) printf("x=%i:%i y=%i:%i\n",xrange,xrange,yrange,yrange)

#stretch the data
if(HasValue(sample)) {
tmp=dcss(data[xrange:xrange,yrange:yrange,],ignore=ign,variance=var,sample=sample[xrange:xrange,yrange:yrange,])
}else{
tmp=dcss(data[xrange:xrange,yrange:yrange,],ignore=ign,variance=var)
}
#if its the first y-chunk then insert else ramp it into the y-column
if(j==1) tmppic[,yrange:yrange,]=tmp
if((j!=1 && ydiff>int(ys/2))) break;
if(j!=1) {
tmppic[,yrange:yrange-int(ys/2)+ydiff]=tmppic[,yrange:yrange-int(ys/2)+ydiff]*(1-yramp)
tmp[,:int(ys/2),]=tmp[,:int(ys/2),]*yramp
tmppic[,yrange:yrange,]=tmppic[,yrange:yrange,]+tmp
}
ydiff=0
}

#if its the first x-chunk then insert else ramp it into the full picutre
if(i==1) pic[xrange:xrange,]=tmppic
if((i!=1 && xdiff>int(xs/2))) return(bip(short(pic)))
if(i!=1) {
pic[xrange:xrange-int(xs/2)+xdiff,]=pic[xrange:xrange-int(xs/2)+xdiff,]*(1-xramp)
tmppic[:int(xs/2),]=tmppic[:int(xs/2),]*xramp
pic[xrange:xrange,]=pic[xrange:xrange,]+tmppic
}
}
return(bip(short(pic)))
}

define echoarraystats(arrayName){
if((\$ARGC != 1) && (\$ARGC != 2)) {
printf("\nError: no array was passed.\n")
return(null)
}

arrName = ""

if (HasValue(arrayName)==1) {
arrName = arrayName
}

arrmax = max(\$1, axis=z)
arrmax = max(arrmax, axis=y)
arrmax = max(arrmax, axis=x)

arrmin = min(\$1, axis=z)
arrmin = min(arrmin, axis=y)
arrmin = min(arrmin, axis=x)

arravg = min(\$1, axis=z)
arravg = min(arravg, axis=y)
arravg = min(arravg, axis=x)

if (HasValue(arrName)==1) {
printf("Array: %s - ", arrName)
}

printf("min: %f ", arrmin)
printf("max: %f", arrmax)
printf("avg: %f\n", arravg)
}

define dcsf(ignore,sample,variance) {
if (\$ARGC != 1) {
printf("\ndcsf()\n")
printf("performs a decorrelation stretch of data using floating point values.\n")
printf("data will be centered around 0.5 with a sigma of +- (50.0 / 255.0).\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=0, 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 set 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.0 / 255.0) \n\n")
printf("Last updated 2010--06-25.\n")
printf("Based on dcs by Noel Gorelick, modified by Ben Lincoln.\n")
printf("Magic numbers in this function all the floating-point equivalents of the byte values in Noel's function.\n")
return(null)
}

#
# this does a decorrelation stretch
# Data will be centered around 0.5, with a sigma of +/- (50.0 / 255.0).
#

# compute the symmetric covariance matrix and eigenvectors of it.
a = \$1;

if (HasValue(sample)) {
if (HasValue(ignore)) {
c = covar(sample, axis=z, ignore=ignore);
} else {
c = covar(sample, axis=z);
}
} else {
if (HasValue(ignore)) {
c = covar(a, axis=z, ignore=ignore);
} else {
c = covar(a, axis=z);
}
}
e = eigen(c);
c = 0;

# 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 / 255.0);
}

# wrapped the reference to e in abs() to avoid values involving imaginary
# numbers - their presence would cause the resulting image to be a black rectangle
# I believe their presence is related to a bug in DaVinci regarding very small negative
# floating-point values.
#s = identity(dim(a)) * v/sqrt(e);

# following line split up to optimize memory use
#s = identity(dim(a)) * v/sqrt(abs(e));
v1 = abs(e);
v1 = sqrt(v1);
v1 = v / v1;
s = dim(a);
s = identity(s);
s = s * v1;
v1 = 0;

# following line split up to optimize memory use
#m = mxm(mxm(e[2:], s), translate(e[2:], x, y));
m = mxm(e[2:], s);
m2 = translate(e[2:], x, y);
m = mxm(m, m2);
m2 = 0;
s = 0;
e = 0;

# 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), dim(a)*dim(a));

# subtract the mean of each band, rotate and add back 0.5
if (HasValue(ignore)) {
avg = avg(a,y,ignore=ignore);
} else {
avg = avg(a,y,ignore=0);
}
y = 0;
a = a - avg;
avg = 0;

a
m

# uncomment the following line and then comment out the one after that to perform the matrix multiplication using
# double-precision floating-point values - this will use a lot of memory.
#a = mxm(a, m);
a = mxm(a, m, format=float);

m = 0;
a = a + 0.5;
a = bip(a);

# reshape the data back to its original size
resize(a, d, d, d);
return(a);
}

define rdcsf(data,sample,ignore,variance,ysizepercent,xsizepercent,ysize,xsize,verbose) {

if(\$ARGC==0 && HasValue(data)==0) {
printf("\nCompute a Running DCS Stretch using single-precision floating point values.\n")
printf("Calls dcsf() 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=0)\n")
printf("ysizepercent=relative size of the chunks to stretch in y (Default=1/4 of the X dimension of the data)\n")
printf("xsizepercent=relative size of the chunks to stretch in x (Default=1/4 of the Y dimension of the data)\n")
printf("ysize=absolute size of the chunks to stretch in y (Default=none)\n")
printf("xsize=absolute size of the chunks to stretch in x (Default=none)\n")
printf("\nIMPORTANT: The xsize and ysize parameters are included to provide drop-in replacement compatibility ")
printf("with the original rdcs() function. If they are specified, then the values of the xsizepercent and ")
printf("ysizepercent parameters are ignored.\n")
printf("variance=variance of the stretch (Default=(50.0 / 255.0))\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("Last updated 2010-03-21.\n")
printf("Based on rdcs by c.edwards, modified by Ben Lincoln.\n")
printf("Magic numbers in this function are the floating-point equivalents of the byte values in the original function.\n")
return(null)
}

#initialize input
if(HasValue(data)==0) data=\$1
ign=0
var=(50.0 / 255.0)
verb=0
ydiff=xdiff=0
dim=dim(data)
xs=int(dim[1,0,0] * 0.25)
ys=int(dim[2,0,0] * 0.25)

if(HasValue(verbose)) verb=1

if(HasValue(ysizepercent)) ys=int(float(dim[1,0,0]) * ysizepercent)
if(HasValue(xsizepercent)) xs=int(float(dim[2,0,0]) * xsizepercent)

# if absolute tile sizes are specified, then ignore the relative sizes
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=bip(tofloat(data*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-int(xs/2)
maxy=dim-int(ys/2)
if(dim   if(dim
#loop through x to get x-direction
for(i=1;i<=maxx;i+=int(xs/2)) {

xrange=i;xrange=i+xs-1
if(xrange>=dim) {
xdiff=xrange-dim
xrange=dim
}

tmppic=float(clone(0,xrange-i+1,dim,dim))
#tmppic=bip(float(clone(0,xrange-i+1,dim,dim)))

#loop through y to get y-direction
for(j=1;j<=maxy;j+=int(ys/2)) {

yrange=j;yrange=j+ys-1
if(yrange>=dim) {
ydiff=yrange-dim
yrange=dim
}

#next line commented out for efficiency
#if(verb==1) printf("x=%i:%i y=%i:%i\n",xrange,xrange,yrange,yrange)

#stretch the data
if(HasValue(sample)) {
tmp=dcsf(data[xrange:xrange,yrange:yrange,],ignore=ign,variance=var,sample=sample[xrange:xrange,yrange:yrange,])
}else{
tmp=dcsf(data[xrange:xrange,yrange:yrange,],ignore=ign,variance=var)
}
#if its the first y-chunk then insert else ramp it into the y-column
if(j==1) tmppic[,yrange:yrange,]=tmp
if((j!=1 && ydiff>int(ys/2))) break;
if(j!=1) {
tmppic[,yrange:yrange-int(ys/2)+ydiff]=tmppic[,yrange:yrange-int(ys/2)+ydiff]*(1-yramp)
tmp[,:int(ys/2),]=tmp[,:int(ys/2),]*yramp
tmppic[,yrange:yrange,]=tmppic[,yrange:yrange,]+tmp
}
ydiff=0
}

#if its the first x-chunk then insert else ramp it into the full picture
if(i==1) pic[xrange:xrange,]=tmppic
if((i!=1 && xdiff>int(xs/2))) return(bip(pic))
if(i!=1) {
pic[xrange:xrange-int(xs/2)+xdiff,]=pic[xrange:xrange-int(xs/2)+xdiff,]*(1-xramp)
tmppic[:int(xs/2),]=tmppic[:int(xs/2),]*xramp
pic[xrange:xrange,]=pic[xrange:xrange,]+tmppic
}
}
return(bip(pic))
}

define dvrgb2yuv(r, g, b, minval, maxval, noscale, weightred, weightblue, nouvscale) {
if(\$ARGC == 1) {
rgb = float(\$1)
} else {
if ((HasValue(rgb) == 0) && (HasValue(r) == 0 || HasValue(g) == 0 || HasValue(b) == 0)) {
printf("\nConverts RGB colour data to YUV (Luminance/Chrominance).\n")
printf("\nAccepts either a single RGB object or three individual channels.\n")
printf("\nReturns values in a range of 0.0 - 1.0 unless the noscale=1 parameter is specified.\n")
printf("\nThe optional minval and maxval arguments allows scaling the input values.\n")
printf("\nFor example, to convert a 16-bit (short) RGB image to YUV:\n")
printf("\nyuv = rgb2yuv(rgb, minval=-32768, maxval=32768)\n")
printf("\nNote: Always returns an object of type Float.\n")
return(null)
} else {
rgb = cat(float(r), float(g), axis=z)
rgb = cat(rgb, float(b), axis=z)
}
}
if (HasValue(minval) == 0) {
minval = 0.0
}
if (HasValue(maxval) == 0) {
maxval = 1.0
}
if (HasValue(noscale) == 0) {
noscale = 0
}
if (HasValue(nouvscale) == 0) {
nouvscale = 0
}
if (HasValue(weightred) == 0) {
weightred = 0.299
}
if (HasValue(weightblue) == 0) {
weightblue = 0.114
}
weightgreen = 1.0 - weightred - weightblue

# not sure if these are calculated or magic
umax = 0.436
vmax = 0.615

if (noscale == 0) {
range = maxval - minval
if (minval != 0.0) {
rgb = rgb - minval
}
if (range != 1.0) {
rgb = rgb / range
}
}

r=rgb[0,0,1]
g=rgb[0,0,2]
b=rgb[0,0,3]

y = (r * weightred) + (g * weightgreen) + (b * weightblue)
u = umax * ((b - y) / (1 - weightblue))
v = vmax * ((r - y) / (1 - weightred))

if (nouvscale == 0) {
u = u + umax
u = u / (umax * 2)
v = v + vmax
v = v / (vmax * 2)
}

yuv = cat(y, u, axis=z)
yuv = cat(yuv, v, axis=z)
return(yuv)
}

define dvyuv2rgb(y, u, v, minval, maxval, noscale, weightred, weightblue, nouvscale) {
if(\$ARGC == 1) {
yuv = float(\$1)
} else {
if ((HasValue(yuv) == 0) && (HasValue(y) == 0 || HasValue(u) == 0 || HasValue(v) == 0)) {
printf("\nConverts YUV colour data to RGB.\n")
printf("\nAccepts either a single YUV object or three individual channels.\n")
printf("\nReturns values in a range of 0.0 - 1.0 unless the noscale=1 parameter is specified.\n")
printf("\nThe optional minval and maxval arguments allows scaling the output values.\n")
printf("\nFor example, to convert YUV values to a 16-bit (short) RGB image:\n")
printf("\nrgb = short(yuv2rgb(yuv, minval=-32768, maxval=32768))\n")
printf("\nNote: Always returns an object of type Float.\n")
return(null)
} else {
yuv = cat(float(y), float(u), axis=z)
yuv = cat(yuv, float(v), axis=z)
}
}
if (HasValue(minval) == 0) {
minval = 0.0
}
if (HasValue(maxval) == 0) {
maxval = 1.0
}
if (HasValue(noscale) == 0) {
noscale = 0
}
if (HasValue(nouvscale) == 0) {
nouvscale = 0
}
if (HasValue(weightred) == 0) {
weightred = 0.299
}
if (HasValue(weightblue) == 0) {
weightblue = 0.114
}
weightgreen = 1.0 - weightred - weightblue

# not sure if these are calculated or magic
umax = 0.436
vmax = 0.615

y = yuv[0,0,1]
u = yuv[0,0,2]
v = yuv[0,0,3]

if (nouvscale == 0) {
u = u * (umax * 2)
u = u - umax
v = v * (vmax * 2)
v = v - vmax
}

r = y + (v * ((1 - weightred) / vmax))
g = y - (u * ((weightblue * (1 - weightblue)) / (umax * weightgreen))) - (v * ((weightred * (1 - weightred)) / (vmax * weightgreen)))
b = y + (u * ((1 - weightblue) / umax))

rgb = cat(r, g, axis=z)
rgb = cat(rgb, b, axis=z)

if (noscale == 0) {
range = maxval - minval
if (range != 1.0) {
rgb = rgb * range
}
if (minval != 0.0) {
rgb = rgb + minval
}
}

return(rgb)
}

define blendcolourrgb() {
if(\$ARGC != 2) {
printf("\nPerforms an operation very similar to the 'Colour' blend mode in Photoshop.\n")
printf("\n(This is a wrapper for the blendcolourhsy() function).\n")
printf("\nRequires two NxNx3 objects representing images in an RGB colourspace,\n")
printf("or one NxNx1 object representing a greyscale luminance image and one NxNx3 \n")
printf("object representing an image in the RGB colourspace.\n")
printf("\nThe colour of object 2 is applied to the luminance of object 1.\n")
printf("\nExample: blended = colourblendrgb(baseimage, colourimage).\n\n")
return(null)
}

# comment out this next line to use the old version of this function
return(blendcolourpdfspec(\$1, \$2))

base = \$1
colour = \$2

if((dim(base)[3,1,1] != 1) && (dim(base)[3,1,1] != 3)) {
printf("\nError: the luminance object must have either one or three data planes.\n")
return(null)
}

if((dim(colour)[3,1,1] != 3)) {
printf("\nError: the colour object must have three data planes.\n")
return(null)
}

if(dim(base)[3,1,1] == 1) {
base = threeclone(base)
}
return hsy2rgb(blendcolourhsy(rgb2hsy(tofloat(base)),rgb2hsy(tofloat(colour))))
}

# Applies a transformation very similar to the "colour" blend mode in Photoshop
define blendcolourhsy() {
if(\$ARGC != 2) {
printf("\nPerforms an operation very similar to the 'Colour' blend mode in Photoshop.\n")
printf("\nRequires two NxNx3 objects representing images in an HSY colourspace.\n")
printf("\nThe colour of object 2 is applied to the luminance of object 1.\n")
printf("\nMore specifically, the resulting object consists of:\n")
printf("\n- The luminance of object 1\n")
printf("\n- The hue of object 2\n")
printf("\n- A saturation which is the inverse of the absolute value of object 1's\n")
printf("  luminosity channel's distance from object 2's luminosity channel, capped at the value of\n")
printf("  object 2's saturation.\n")
printf("\nExample: blended = colourblendhsy(baseimage, colourimage).\n\n")
return(null)
}

# comment out this next line to use the old version of this function
return(rgb2hsy(blendcolourpdfspec(hsy2rgb(\$1), hsy2rgb(\$2))))

base = \$1
colour = \$2

if((dim(base)[3,1,1] != 3) || (dim(colour)[3,1,1] != 3)) {
printf("\nError: both input objects must have three data planes (Hue, Saturation, and Luminosity).\n")
return(null)
}
if((dim(base)[1,1,1] != dim(colour)[1,1,1]) || (dim(base)[2,1,1] != dim(colour)[2,1,1])) {
printf("\nError: the input objects must be of equal dimensions.\n")
return(null)
}

y = base[0,0,3]
h = colour[0,0,1]
s = base[0,0,3]
cy = colour[0,0,3]
#s = abs(s - 0.5)
#s = s * 2.0

# note: uncomment the following line (remove the pound sign at the beginning) to use a simple HSY colour blend
# return(cat(cat(colour[0,0,1], colour[0,0,2], axis=z), base[0,0,3], axis=z))

# note 2: this will improve processing times, but will cause significantly different results in most cases

# delta from colour luminance
s = s - cy
spos = s
spos[0,0,1] = 0.0
spos = cat(spos, s, axis=z)
spos = max(spos, axis=z)

sneg = s
sneg[0,0,1] = 0.0
sneg = cat(sneg, s, axis=z)
sneg = abs(min(sneg, axis=z))

# rescale the delta values
sposceil = cy
sposceil[0,0,1] = 1.0
sposceil = sposceil - cy
spos = spos / sposceil

sneg = sneg / cy

cy = 0

s = cat(spos, sneg, axis=z)

spos = 0
sneg = 0

s = max(s, axis=z)

# clip saturation to a range of 0 - 1
s2 = s
s2[0,0,1] = 1.0
s = cat(s, s2, axis=z)
s = min(s, axis=z)
s2[0,0,1] = 0.0
s = cat(s, s2, axis=z)
s = max(s, axis=z)
# end clip
# invert the calculated saturation value
s2 = s
s2[0,0,1] = 1.0
s = s2 - s
s2 = 0
# end inversion
# cap the saturation at the value from the colour image
s = cat(s, colour[0,0,2], axis=z)
s = min(s, axis=z)
# end cap

hsy = cat(h, s, axis=z)
hsy = cat(hsy, y, axis=z)
return(hsy)
}

define blendcolourpdfspec() {
if(\$ARGC != 2) {
printf("\nPerforms an operation based on the 'Colour' blend mode as described in Adobe's PDF specification.\n")
printf("\nRequires two NxNx3 objects representing images in an RGB colourspace,\n")
printf("or one NxNx1 object representing a greyscale luminance image and one NxNx3 \n")
printf("object representing an image in the RGB colourspace.\n")
printf("\nThe colour of object 2 is applied to the luminance of object 1.\n")
printf("\nExample: blended = blendcolourpdfspec(baseimage, colourimage).\n\n")
return(null)
}

base = \$1
#colour = \$2

if((dim(base)[3,1,1] != 1) && (dim(base)[3,1,1] != 3)) {
printf("\nError: the luminance object must have either one or three data planes.\n")
return(null)
}

if((dim(\$2)[3,1,1] != 3)) {
printf("\nError: the colour object must have three data planes.\n")
return(null)
}

if(dim(base)[3,1,1] == 1) {
base = threeclone(base)
}

basey = rgb2hsy(tofloat(base))[0,0,3]
coloury = rgb2hsy(tofloat(\$2))[0,0,3]

base = 0

# following line expanded to reduce memory usage:
#result = \$2 + (basey - coloury)
result = \$2
result = result + basey
result = result - coloury

#colour = 0
coloury = 0

# uncomment the next two lines to perform a simple HSY luminance-channel swap
#result = clipdata(result, min = 0.0, max = 1.0)
#return(result)

colourmin = min(result, axis = z)
colourmax = max(result, axis = z)

resultmin = threeclone(basey)
resultmax = resultmin

# following line expanded to reduce memory usage:
#resultmin = resultmin + (((result - basey) * basey) / (basey - colourmin))
r1 = result - basey
r1 = r1 * basey
#r2 = basey - colourmin
#r1 = r1 / r2
#r2 = 0
r1 = r1 / (basey - colourmin)
resultmin = resultmin + r1
r1 = 0

one = basey
one[0,0,0] = 1.0

# following line expanded to reduce memory usage:
#resultmax = resultmax + (((result - basey) * (one - basey)) / (colourmax - basey))
r1 = result - basey
r2 = one - basey
r1 = r1 * r2
#r2 = colourmax - basey
#r1 = r1 / r2
#r2 = 0
r1 = r1 / (colourmax - basey)
resultmax = resultmax + r1
r1 = 0

basey = 0

# following line expanded to reduce memory usage:
#minmask = ceil(abs(clipdata(colourmin, min=-1.0, max=0.0, quiet=1)))
minmask = clipdata(colourmin, min=-1.0, max=0.0, quiet=1)

colourmin = 0

# following line expanded to reduce memory usage:
#maxmask = ceil(clipdata((colourmax - 1.0), min=0.0, max=1.0, quiet=1))
maxmask = colourmax - 1.0

colourmax = 0

#one[0,0,0] = 1.0

# following line expanded to reduce memory usage:

one = 0
result = result * noclipmask

# in the PDF specification, the max case takes precedence over the min case
# following line expanded to reduce memory usage:

# following line expanded to reduce memory usage:
#result = (result * noclipmask) + (resultmin * minmask) + (resultmax * maxmask)
# following line moved up so noclipmask can be discarded earlier
#result = result * noclipmask
r1 = resultmin * minmask
resultmin = 0
result = result + r1
r1 = resultmax * maxmask
resultmax = 0
result = result + r1
r1 = 0

return(result)
}

# Converts numeric values to double-precision floating-point in the range of 0.0 - 1.0
define todouble(minval, maxval) {
if(\$ARGC != 1) {
printf("\nConverts numeric values to double-precision floating-point in the range of 0.0 - 1.0 (or the range specified by the optional minval and maxval parameters).\n")
printf("\nWill automatically convert byte, short, and integer types.\n")
printf("\nFloat and double types require that the maximum and minimum values be specified, or no scaling will take place.\n\n")
return(null)
}

inval = \$1

foundtype = 0

intype = type(inval)

if (intype == "byte") {
foundtype = 1
if (HasValue(minval) == 0) {
minval = 0.0
}
if (HasValue(maxval) == 0) {
maxval = 255.0
}
}

if (intype == "short") {
foundtype = 1
if (HasValue(minval) == 0) {
#minval = -32767.0
minval = -32768.0
}
if (HasValue(maxval) == 0) {
maxval = 32767.0
}
}

if (intype == "int") {
foundtype = 1
if (HasValue(minval) == 0) {
#minval = -2147483647.0
minval = -2147483648.0
}
if (HasValue(maxval) == 0) {
maxval = 2147483647.0
}
}

if (foundtype == 0) {
if (HasValue(minval) == 0) {
minval = 0.0
}
if (HasValue(maxval) == 0) {
maxval = 1.0
}
if ((intype == "float") || (intype == "double")) {
foundtype = 1
}
}

if (foundtype == 0) {
printf("\nThe object type %s cannot be converted by this function. No conversion will take place.\n", intype)
return(inval)
}

if ((intype == "double") && (minval == 0.0) && (maxval == 1.0)) {
printf("\nThe output type is identical to the input type. No conversion will take place.\n")
return(inval)
}

if (intype == "double") {
doubled = inval
} else {
doubled = double(inval)
}

range = maxval - minval
if (minval != 0.0) {
doubled = doubled - minval
}
if (range != 1.0) {
doubled = doubled / range
}

return(doubled)
}

# a wrapper for todouble for use when only single-precision is required
# note: for improved processing speed, a modified copy of the todouble() function is used in the production version.
define tofloatdev(minval, maxval) {
if(\$ARGC != 1) {
printf("\nConverts numeric values to single-precision floating-point in the range of 0.0 - 1.0.\n")
printf("\nWill automatically convert byte, short, and integer types.\n")
printf("\nFloat and double types require that the maximum and minimum values be specified, or no scaling will take place.\n\n")
return(null)
}

inval = \$1

intype = type(inval)

if (intype == "float") {
printf("\nThe output type is identical to the input type. No conversion will take place.\n")
return(inval)
}

if ((HasValue(minval) == 1) && (HasValue(maxval) == 1)) {
minv = minval
maxv = maxval
return(float(todouble(\$1, minval = minv, maxval = maxvl)))
}
if (HasValue(maxval) == 1) {
maxv = maxval
return(float(todouble(\$1, maxval = maxval)))
}
if (HasValue(minval) == 1) {
minv = minval
return(float(todouble(\$1, minval = minval)))
}
return(float(todouble(\$1)))
}

# Converts numeric values to double-precision floating-point in the range of 0.0 - 1.0
define tofloat(minval, maxval) {
if(\$ARGC != 1) {
printf("\nConverts numeric values to single-precision floating-point in the range of 0.0 - 1.0.\n")
printf("\nWill automatically convert byte, short, and integer types.\n")
printf("\nThe optional minval and maxval parameters specify the minimum and maximum values to use for the input data. ")
printf("For integer input, if they are not specified, they default to the minimum and maximum that can be represented by that integer type. ")
printf("E.g, 0 to 255 for bytes, -32768 to 32767 for shorts.\n")
printf("\nFloat and double types require that the maximum and minimum values be specified, or no scaling will take place.\n\n")
return(null)
}

inval = \$1

foundtype = 0

intype = type(inval)

if (intype == "byte") {
foundtype = 1
if (HasValue(minval) == 0) {
minval = 0.0
}
if (HasValue(maxval) == 0) {
maxval = 255.0
}
}

if (intype == "short") {
foundtype = 1
if (HasValue(minval) == 0) {
#minval = -32767.0
minval = -32768.0
}
if (HasValue(maxval) == 0) {
maxval = 32767.0
}
}

if (intype == "int") {
foundtype = 1
if (HasValue(minval) == 0) {
#minval = -2147483647.0
minval = -2147483648.0
}
if (HasValue(maxval) == 0) {
maxval = 2147483647.0
}
}

if (foundtype == 0) {
if (HasValue(minval) == 0) {
minval = 0.0
}
if (HasValue(maxval) == 0) {
maxval = 1.0
}
if ((intype == "float") || (intype == "double")) {
foundtype = 1
}
}

if (foundtype == 0) {
printf("\nThe object type %s cannot be converted by this function. No conversion will take place.\n", intype)
return(inval)
}

if ((intype == "float") && (minval == 0.0) && (maxval == 1.0)) {
printf("\nThe output type is identical to the input type. No conversion will take place.\n")
return(inval)
}

if (intype == "float") {
floated = inval
} else {
floated = float(inval)
}

range = maxval - minval

if (minval != 0.0) {
floated = floated - minval
}
if (range != 1.0) {
floated = floated / range
}

return(floated)
}

# Converts floating-point values to integers
define tointeger(integertype, minval, maxval) {
if(\$ARGC != 1) {
printf("\nConverts floating-point values to integers (or integers to other types of integers). The range defaults to the maximum range for the specified integer type (0 - 255 for bytes, +/-32767 for shorts, etc.)\n")
printf("\nThis range can be overridden using the option minval and maxval parameters.\n\n")
return(null)
}

inval = \$1
intype = type(inval)

if (HasValue(integertype) == 0) {
printf("\nThe type of integer ('byte', 'short', or 'int') must be specified using the integertype parameter. No conversion will take place.\n")
return(inval)
}

if (intype == integertype) {
printf("\nThe output type is identical to the input type. No conversion will take place.\n")
return(inval)
}

if ((intype != "float") && (intype != "double")) {
#outval = todouble(inval)
outval = tofloat(inval)
} else {
# do nothing - probably no need to convert doubles to floats or vice-versa
#outval = double(inval)
#outval = float(inval)
outval = inval
}

foundtype = 0

if (integertype == "byte") {
foundtype = 1
if (HasValue(minval) == 0) {
minval = 0.0
}
if (HasValue(maxval) == 0) {
maxval = 255.0
}
if (minval < 0.0) {
printf("\nWarning: the minimum value specified (%d) is less than the value that can be represented by the 'byte' object type. Values lower than 0 will be truncated to 0.\n", minval)
}
if (maxval > 255.0) {
printf("\nWarning: the maximum value specified (%d) is greater than the value that can be represented by the 'byte' object type. Values greater than 255 will be truncated to 255.\n", maxval)
}
}

if (integertype == "short") {
foundtype = 1
if (HasValue(minval) == 0) {
#minval = -32767.0
minval = -32768.0
}
if (HasValue(maxval) == 0) {
maxval = 32767.0
}
# if minval < -32767.0
if (minval < (-32768.0)) {
printf("\nWarning: the minimum value specified (%d) is less than the value that can be represented by the 'short' object type. Values lower than -32767 will be truncated to -32767.\n", minval)
}
if (maxval > 32767.0) {
printf("\nWarning: the maximum value specified (%d) is greater than the value that can be represented by the 'short' object type. Values greater than 32767 will be truncated to 32767.\n", maxval)
}
}

if (integertype == "int") {
foundtype = 1
if (HasValue(minval) == 0) {
#minval = -2147483647.0
minval = -2147483648.0
}
if (HasValue(maxval) == 0) {
maxval = 2147483647.0
}
# if minval < 0.0
if (minval < (-2147483648.0)) {
printf("\nWarning: the minimum value specified (%d) is less than the value that can be represented by the 'int' object type. Values lower than -2147483647 will be truncated to 2147483647.\n", minval)
}
if (maxval > 2147483647.0) {
printf("\nWarning: the maximum value specified (%d) is greater than the value that can be represented by the 'int' object type. Values greater than 2147483647 will be truncated to 2147483647.\n", maxval)
}
}

if (foundtype == 0) {
printf("\nThis function cannot convert to the object type %s. No conversion will take place.\n", intype)
return(inval)
}

range = maxval - minval

if (range != 1.0) {
outval = outval * range
}
if (minval != 0.0) {
outval = outval + minval
}

if (integertype == "byte") {
outval = byte(outval)
}

if (integertype == "short") {
outval = short(outval)
}

if (integertype == "int") {
# Hack/workaround to avoid wraparound due to relatively low accuracy of single-precision floating-point numbers
# Build a mask of values in which the temporary (still floating-point) value is greater than where the output value
# (which has been converted to an integer) is converted back to a float, then substitute the maximum int value
# wherever that mask is true
outvaltemp = outval
outval = int(outval)
imask = int(normalize(outvaltemp - float(outval)))
outval = outval * imaskinv
outval = outval + imask
}

return(outval)
}

# a wrapper for tointeger that converts to 32-bit int
define toint(minval, maxval) {
if ((HasValue(minval) == 1) && (HasValue(maxval) == 1)) {
minv = minval
maxv = maxval
return(tointeger(\$1, integertype = "int", minval = minv, maxval = maxvl))
}
if (HasValue(maxval) == 1) {
maxv = maxval
return(tointeger(\$1, integertype = "int", maxval = maxvl))
}
if (HasValue(minval) == 1) {
minv = minval
return(tointeger(\$1, integertype = "int", minval = minv))
}
return(tointeger(\$1, integertype = "int"))

}

# a wrapper for tointeger that converts to 16-bit short
define toshort(minval, maxval) {
if ((HasValue(minval) == 1) && (HasValue(maxval) == 1)) {
minv = minval
maxv = maxval
return(tointeger(\$1, integertype = "short", minval = minv, maxval = maxvl))
}
if (HasValue(maxval) == 1) {
maxv = maxval
return(tointeger(\$1, integertype = "short", maxval = maxvl))
}
if (HasValue(minval) == 1) {
minv = minval
return(tointeger(\$1, integertype = "short", minval = minv))
}
return(tointeger(\$1, integertype = "short"))
}

# a wrapper for tointeger that converts to 8-bit byte
define tobyte(minval, maxval) {
if ((HasValue(minval) == 1) && (HasValue(maxval) == 1)) {
minv = minval
maxv = maxval
return(tointeger(\$1, integertype = "byte", minval = minv, maxval = maxvl))
}
if (HasValue(maxval) == 1) {
maxv = maxval
return(tointeger(\$1, integertype = "byte", maxval = maxvl))
}
if (HasValue(minval) == 1) {
minv = minval
return(tointeger(\$1, integertype = "byte", minval = minv))
}
return(tointeger(\$1, integertype = "byte"))
}

# concatenates three planes of the same dimensions and type into a single cube
# does not alter the object type
define threeplane() {
if(\$ARGC != 3) {
printf("\nConcatenates three planes of the same x/y dimensions into a single cube, but does not alter the data type.\n")
return(null)
}

p1 = \$1
p2 = \$2
p3 = \$3

p1type = type(p1)
p2type = type(p2)
p3type = type(p3)

p1org = org(p1)
p2org = org(p2)
p3org = org(p3)

p1dim = dim(p1)
p2dim = dim(p2)
p3dim = dim(p3)

# make sure each object has only a single plane
if (p1dim[3,1,1] != 1) {
printf("\nError: The first specified object has a z-dimension of %d. Only single planes are allowed.\n", p1dim[3,1,1])
return(null)
}
if (p2dim[3,1,1] != 1) {
printf("\nError: The second specified object has a z-dimension of %d. Only single planes are allowed.\n", p2dim[3,1,1])
return(null)
}
if (p3dim[3,1,1] != 1) {
printf("\nError: The third specified object has a z-dimension of %d. Only single planes are allowed.\n", p3dim[3,1,1])
return(null)
}

# verify that the dimensions match
if ((p1dim != p2dim) || (p1dim != p3dim) || (p2dim != p3dim)) {
printf("\nError: The specified objects do not have the same dimensions.\n")
printf("\tObject 1: x = %d, y = %d", p1dim[1,1,1], p1dim[2,1,1])
printf("\tObject 2: x = %d, y = %d", p2dim[1,1,1], p2dim[2,1,1])
printf("\tObject 3: x = %d, y = %d\n", p3dim[1,1,1], p3dim[2,1,1])
return(null)
}

# verify that the types match
if ((p1type != p2type) || (p1type != p3type) || (p2type != p3type)) {
printf("\nError: The specified objects are not of the same type.\n")
printf("\tObject 1: %s", p1type)
printf("\tObject 2: %s", p2type)
printf("\tObject 3: %s\n", p3type)
return(null)
}

# verify that the organizations match
if ((p1org != p2org) || (p1org != p3org) || (p2org != p3org)) {
printf("\nError: The specified objects are not of the same data organization.\n")
printf("\tObject 1: %s", p1org)
printf("\tObject 2: %s", p2org)
printf("\tObject 3: %s\n", p3org)
return(null)
}

outval = cat(p1, p2, axis = z)
outval = cat(outval, p3, axis = z)

return(outval)
}

# this is a workaround for an apparent bug in the DaVinci scripting syntax
# if the same object is passed multiple times as an argument to a function, only the first occurence is visible to the script code
define threeclone() {
if(\$ARGC != 1) {
printf("\nCreates a three-plane clone of a single plane.\n")
return(null)
}

plane1 = \$1
plane2 = \$1
plane3 = \$1

return (threeplane(plane1, plane2, plane3))
}

define momentcubesimple(){
if (\$ARGC != 1) {
printf("\nCreates a higher-order image cube whose contents are the result of passing the input image cube through the standard statistical equations, ")
printf("as well as the median and range values.\n")
printf("\n")
printf("\$1 = the input image cube\n")
printf("\n")
printf("Plane  1 = Minimum\n")
printf("Plane  2 = Maximum\n")
printf("Plane  3 = Average (Mean)\n")
printf("Plane  4 = Average Deviation\n")
printf("Plane  5 = Standard Deviation\n")
printf("Plane  6 = Variance\n")
printf("Plane  7 = Skewness\n")
printf("Plane  8 = Kurtosis\n")
printf("Plane  9 = Median\n")
printf("Plane 10 = Range\n")
printf("\n")
return(null)
}

sourcecube = \$1

xsize = dim(sourcecube)[1,1,1]
ysize = dim(sourcecube)[2,1,1]
zsize = dim(sourcecube)[3,1,1]

resultcube = clone(sourcecube[1,1,1], xsize, ysize, 10)

ycurrent = 1

# standard moments
while (ycurrent <= ysize) {
xcurrent = 1
while (xcurrent <= xsize) {
resultcube[xcurrent, ycurrent, 1:10] = resize(momentext(sourcecube[xcurrent, ycurrent, 0]), 1, 1, 10)
xcurrent++
}
ycurrent++
}

return(resultcube)
}

define momentcube(windowSize){
if (\$ARGC != 1) {
printf("\nCreates a higher-order image cube whose contents are the result of passing the input image cube through the standard statistical equations.\n")
printf("as well as several semi- and non-standard equations.\n")
printf("\n")
printf("\$1         = the input image cube\n")
printf("windowSize = x/y dimension of the window to use for calculating stats (default = 1)\n")
printf("\n")
printf("Plane  1 = Minimum\n")
printf("Plane  2 = Maximum\n")
printf("Plane  3 = Average (Mean)\n")
printf("Plane  4 = Average Deviation\n")
printf("Plane  5 = Standard Deviation\n")
printf("Plane  6 = Variance\n")
printf("Plane  7 = Skewness\n")
printf("Plane  8 = Kurtosis\n")
printf("Plane  9 = Median\n")
printf("Plane 10 = Range\n")
printf("\n")
return(null)
}

xsize = dim(\$1)[1,1,1]
ysize = dim(\$1)[2,1,1]
zsize = dim(\$1)[3,1,1]

winSize = 1

if(HasValue(windowSize) == 1) {
winSize = windowSize
}

if (winSize == 1) {
return(momentcubesimple(\$1))
}

if (winSize % 2 == 0) {
winSize = winSize + 1
}

resultcube = clone((\$1)[1,1,1], xsize, ysize, 10)

resultcube[0, 0, 0] = 0.0

# extended standard moments
printf("\nDebug: getting extended standard moments.\n")

yRange=xRange=clone(0,2,1,1)
xStart = 0
xEnd = 0
yStart = 0
yEnd = 0
xcurrent = 1
ycurrent = 1
wDiff = (winSize - 1) / 2
if (winSize == 1) {
while (ycurrent <= ysize) {
xcurrent = 1
while (xcurrent <= xsize) {
# the following line was expanded to reduce memory use:
#resultcube[xcurrent, ycurrent, 1:10] = resize(momentext(sourcecube[xcurrent, ycurrent, 0]), 1, 1, 10)
r1 = momentext((\$1)[xcurrent, ycurrent, 0])
resultcube[xcurrent, ycurrent, 1:10] = resize(r1, 1, 1, 10)
r1 = 0
xcurrent++
}
ycurrent++
}
} else {
while (ycurrent <= ysize) {
#yRange = max(cat(1, ycurrent - wDiff, axis = x), axis = x)
#yRange = min(cat(ysize, ycurrent + wDiff, axis = x), axis = x)
yRange = ycurrent - wDiff
if (yRange < 1) {
yRange = 1
}
yRange = ycurrent + wDiff
if (yRange > ysize) {
yRange = ysize
}

xcurrent = 1

while (xcurrent <= xsize) {

#xRange = max(cat(1, xcurrent - wDiff, axis = x), axis = x)
#xRange = min(cat(xsize, xcurrent + wDiff, axis = x), axis = x)
xRange = xcurrent - wDiff
if (xRange < 1) {
xRange = 1
}
xRange = xcurrent + wDiff
if (xRange > xsize) {
xRange = xsize
}

#yRange = max(cat(1, ycurrent - wDiff, axis = x), axis = x)
#yRange = min(cat(ysize, ycurrent + wDiff, axis = x), axis = x)

# the following line was expanded to reduce memory use:
#resultcube[xcurrent, ycurrent, 1:10] = resize(momentext(sourcecube[xRange:xRange, yRange:yRange, 0]), 1, 1, 10)
r1 = momentext((\$1)[xRange:xRange, yRange:yRange, 0])
resultcube[xcurrent, ycurrent, 1:10] = resize(r1, 1, 1, 10)
r1 = 0

xcurrent++
}
ycurrent++
}
}

sourcecube = 0
xRange = 0
yRange = 0

return(resultcube)
}

define dras(lbound, ubound) {
printf("\nData Reposition And Scale")
if (\$ARGC != 1) {
printf("\n\nA utility function for repositioning the minimum value in a set of input data and then (optionally) rescaling the data ")
printf("so that the maximum value is equal to the specified new upper bound.\n")
printf("\n")
printf("\$1    = the input array\n")
printf("lbound    = the new lower bound for the data (default = 0.0)\n")
printf("ubound    = the new upper bound for the data (default = 1.0)\n")
return(null)
}

newLBound = 0.0
newUBound = 1.0

if(HasValue(lbound) == 1) {
newLBound = lbound
}

if(HasValue(ubound) == 1) {
newUBound = ubound
}

inVal = float(\$1)

currentMin = min(inVal)
currentMax = max(inVal)

reposition = 1
rescale = 1

if (currentMin == newLBound) {
reposition = 0
}

if ((currentMax + currentMin) == newUBound) {
rescale = 0
}

if ((reposition == 0) && (rescale == 0)) {
printf(" - the input data already occupies the specified range. No change is necessary.")
return(inVal)
}

if ((reposition == 0) && (rescale == 1)) {
printf(" - scaling data to the new upper bound.")
return(inVal)
}

if ((reposition == 1) && (rescale == 0)) {
printf(" - repositioning data to the new lower bound.")
return(inVal)
}

if ((reposition == 1) && (rescale == 1)) {
printf(" - repositioning and scaling data to match the new upper and lower bounds.")
return(inVal)
}

outVal = (inVal)

norm = 1

if ((currentMin == 0.0) && (currentMax == 1.0)) {
norm = 0
}

if (norm == 1) {
outVal = normalize(outVal)
}

if (rescale == 1) {
dataRange = currentMax - currentMin
outputRange = newUBound - newLBound
scalingFactor = outputRange / dataRange
outVal = outVal * scalingFactor
}

if (reposition == 1) {
outVal = outVal + newLBound
}

return(outVal)
}

define anormalize(steps, backlash, threshold, threshrepeat, searchwidth){
if (\$ARGC != 1) {
printf("\nAggressively normalizes the specified array.\n")
printf("\n")
printf("\$1           = the input array\n")
printf("steps          = the number of steps between 0 and 1 to use when generating the histogram (default = 100)\n")
printf("backlash     = the number of steps to \"back off\" the new min and max values when rescaling (default = 1)\n")
printf("threshold    = the relative count that must be exceeded by a step in the histogram to trigger a response (default = 0.05)\n")
printf("threshrepeat = the number of steps that must exceed the threshold in order to trigger a response (default = 3)\n")
printf("searchwidth  = the number of steps after the first which exceeds the threshold that will be searched to find a repetition (default = 5)\n")
return(null)
}

stepCount = 100
bl = 1
thresh = 0.05
trepeat = 3
swidth = 5

if(HasValue(steps) == 1) {
stepCount = steps
}

if(HasValue(backlash) == 1) {
bl = backlash
}

if(HasValue(threshold) == 1) {
thresh = threshold
}

if(HasValue(threshrepeat) == 1) {
trepeat = threshrepeat
}

if(HasValue(searchwidth) == 1) {
swidth = searchwidth
}

printf("Performing an aggressive normalization.\n")

# first, perform a standard normalization to make the rest of the calculations easier
inVal = float(\$1)

inMin = min(inVal)
inMax = max(inVal)

if ((inMin != 0.0) || (inMax != 1.0)) {
inVal = inVal - inMin
inVal = inVal / max(inVal)
}

# get the new values to use
newVals = findAggMinMax(inVal, steps = stepCount, backlash = bl, threshold = thresh, threshrepeat = trepeat, searchwidth = swidth)
newMin = newVals[1,1,1]
newMax = newVals[2,1,1]

#printf("Debug: new minimum: %f.\n", newMin)
#printf("Debug: new maximum: %f.\n", newMax)

dNewMin = dim(newMin)
dNewMinX = dNewMin[1,1,1]
dNewMinY = dNewMin[2,1,1]
dNewMinZ = dNewMin[3,1,1]

#printf("Debug: newMin dimensions\n")
#printf("\t%f\n", dNewMinX)
#printf("\t%f\n", dNewMinY)
#printf("\t%f\n", dNewMinZ)

dNewMax = dim(newMax)
dNewMaxX = dNewMax[1,1,1]
dNewMaxY = dNewMax[2,1,1]
dNewMaxZ = dNewMax[3,1,1]

#printf("Debug: newMax dimensions\n")
#printf("\t%f\n", dNewMaxX)
#printf("\t%f\n", dNewMaxY)
#printf("\t%f\n", dNewMaxZ)

#printf("Debug: outVal = inVal - newMin\n")
if (newMin > 0.0) {
#printf("Debug: subtracting %f from the input data.\n", newMin)
outVal = inVal - newMin
} else {
if (newMin == 0.0) {
#printf("Debug: %f is equal to zero, so it will not be subtracted from the input data.\n", newMin)
} else {
#printf("Debug: %f is less than zero, so it will not be subtracted from the input data - this indicates a calculation error.\n", newMin)
}
outVal = inVal
}

#printf("Debug: newMax = newMax - newMin\n")
newMax = newMax - newMin

#printf("Debug: outVal = outVal / newMax\n")
if (newMax < 1.0) {
#printf("Debug: dividing the input data by %f.\n", newMax)
outVal = outVal / newMax
} else {
#printf("Debug: %f is greater than or equal to 1.0, so the input data will not be divided by it.\n", newMin)
outVal = inVal
}

# crop values outside 0.0 - 1.0

#printf("Debug: xMax = dim(outVal)[1,1,1]\n")
xMax = dim(outVal)[1,1,1]

#printf("Debug: yMax = dim(outVal)[2,1,1]\n")
yMax = dim(outVal)[2,1,1]

#printf("Debug: zMax = dim(outVal)[3,1,1]\n")
zMax = dim(outVal)[3,1,1]

xCurrent = 1
yCurrent = 1
zCurrent = 1

#printf("Debug: Final pass to clip values to 0.0 - 1.0.\n")

outVal = clipdata(outVal, min = 0.0, max = 1.0)

return(outVal)

}

define clipdata(min, max, quiet){
if (\$ARGC != 1) {
printf("\nClips the input data to the specified range.\n")
printf("\n")
printf("\$1    = the input data (required)\n")
printf("min    = clip to a minimum of this value (default = min(\$1))\n")
printf("max    = clip to a maximum of this value (default = max(\$1))\n")
printf("quiet    = do not display informational messages (default = 0)\n")
return(null)
}

willChange = 0
noChangeMin = 0
noChangeMax = 0
q = 0
if(HasValue(quiet) == 1) {
if (quiet == 1) {
q = 1
}
}

inVal = float(\$1)

cMax = max(inVal)
cMin = min(inVal)

if(q == 0) {
printf("Current minimum value: %f.\n", cMin)
printf("Current maximum value: %f.\n", cMax)
}

if(HasValue(min) == 1) {
if (min > cMin) {
cMin = min
willChange = 1
} else {
if(q == 0) {
printf("The input data already has a minimum value of %f, so its lower bound will not be clipped.\n", cMin)
}
noChangeMin = 1
}
}

if(HasValue(max) == 1) {
if (max < cMax) {
cMax = max
willChange = 1
} else {
if(q == 0) {
printf("The input data already has a maximum value of %f, so its upper bound will not be clipped.\n", cMax)
}
noChangeMax = 1
}
}

if (willChange == 0) {
if(q == 0) {
printf("Clipping values were not provided. No change will be made to the data.\n")
}
return(inVal)
}

if ((noChangeMin == 1) && (noChangeMax == 1)) {
if(q == 0) {
printf("The input data is already within the specified bounds. No change will be made to the data.\n")
}
return(inVal)
}

if(q == 0) {
printf("Clipping input values to a range of (%f) to (%f).\n", cMin, cMax)
}

outVal = inVal

zMax = dim(outVal)[3,1,1]

zCurrent = 1

zTemp = inVal[0,0,zCurrent]

while (zCurrent <= zMax) {
zTemp = inVal[0,0,zCurrent]
zTemp = cat(zTemp, zTemp, axis = z)
zTemp[0,0,2] = cMin
zTemp = max(zTemp, axis = z)

zTemp = cat(zTemp, zTemp, axis = z)
zTemp[0,0,2] = cMax
zTemp = min(zTemp, axis = z)
outVal[0,0,zCurrent] = zTemp

zCurrent++
}

zTemp = 0

return(outVal)
}

define ndi(){
if (\$ARGC != 2) {
printf("\nGenerates a Normalized Difference Index map based on the input data.\n")
printf("This function can be used to generate any normalized difference index that uses the form ")
printf("x = (a - b) / (a + b).\n")
printf("\n")
printf("\$1    = the value to use for the 'a' variable (required)\n")
printf("\$2      = the value to use for the 'b' variable (required)\n")
printf("\n")
printf("Note: values are in the range of -1.0 to 1.0.\n")

printf("\n")
printf("Examples:\n")
printf("\tNormalized Difference Vegetation Index (a = NIR, b = red)\n")
printf("\tGreen Normalized Difference Vegetation Index (a = NIR, b = green)\n")
printf("\tNormalized Burn Ratio (a = NIR, b = 2220nm)\n")
printf("\tND53 (a = 1625nm, b = red)\n")
printf("\tND54 (a = 1625nm, b = NIR)\n")
printf("\tND57 (a = 1625nm, b = 2220nm)\n")
printf("\tND32 (a = red, b = green)\n")
printf("\n")
printf("When calculating the NDVI, values less than 0 indicate no vegetation.\n")
printf("Values greater than 0 indicate the health/stress level of vegetation, with 1.0 representing dense, healthy vegetation.\n")
printf("\n")
printf("Some sources describe the NDVI calculation using the generic term 'visible light' as the 'b' variable. ")
printf("This can be accomplished by collapsing the red, green, and blue data to a single greyscale band ")
printf("(e.g. using the max(), min(), avg(), calculating the luminance, or other statistical methods).\n")
printf("\n")
return(null)
}

nir = \$1
vis = \$2

nirPlanes = dim(nir)[3,1,1]
visPlanes = dim(vis)[3,1,1]

printf("Generating a Normalized Difference Index image.\n")

if (nirPlanes != 1) {
printf("Error: the first dataset must have exactly one plane. The input value has %i planes.\n", nirPlanes)
return(null)
}

if (visPlanes != 1) {
printf("Error: the second dataset must have exactly one plane. The input value has %i planes.\n", visPlanes)
return(null)
}

outVal = (nir - vis) / (nir + vis)

return(outVal)
}

define evi(l, c1, c2, gain, clipMin, clipMax, wrapMax, wrapTh, smartWrap){
if (\$ARGC != 3) {
printf("\nGenerates an Enhanced Vegetation Index map based on the input data.\n")
printf("\nNote: requires near-infrared, red, and blue spectral band data.\n")
printf("\n")
printf("\$1        = the near-infrared band of the input image (required)\n")
printf("\$2          = the red band of the input image (required)\n")
printf("\$3          = the blue band of the input image (required)\n")
printf("l         = canopy background adjustment (default = 1.0)\n")
printf("c1         = canopy background adjustment (default = 6.0)\n")
printf("c2         = canopy background adjustment (default = 7.5)\n")
printf("gain         = canopy background adjustment (default = 1.0)\n")
printf("clipMin        = clip formula result to a minimum of this value (default = do not clip)\n")
printf("clipMax        = clip formula result to a maximum of this value (default = do not clip)\n")
printf("wrapMax        = wrap data with the maximum value back to the minimum value for a more natural appearance (default = 0)\n")
printf("wrapTh        = range from the actual maximum that will be included when using the wrapMax option (default = 0.0)\n")
printf("smartWrap    = attempts to automatically determine the best point at which to wrap values back to the minimum (default = 0)\n")
printf("\n")
printf("The basic formula and defaults used are as indicated by Professor Paul R. Baumann in \"Tropical Wet Realms Of Central Africa, Part II\".\n")
# note: (http://employees.oneonta.edu/baumanpr/geosat2/Central%20Africa%20II/Central%20Africa%20Part%20II.htm)
printf("\n")
printf("Note: the original formula used for MODIS non-normalized 12-bit integer data specifies a gain of 2.5.\n")
return(null)
}

nir = \$1
red = \$2
blue = \$3

printf("Generating an Enhanced Vegetation Index image.\n")

numPlanes = dim(nir)[3,1,1]
if (numPlanes != 1) {
printf("Error: the near-infrared dataset must have exactly one plane. The input value has %i planes.\n", numPlanes)
return(null)
}

numPlanes = dim(red)[3,1,1]
if (numPlanes != 1) {
printf("Error: the red dataset must have exactly one plane. The input value has %i planes.\n", numPlanes)
return(null)
}

numPlanes = dim(blue)[3,1,1]
if (numPlanes != 1) {
printf("Error: the blue dataset must have exactly one plane. The input value has %i planes.\n", numPlanes)
return(null)
}

lval = 1.0
c1val = 6.0
c2val = 7.5
g = 1.0
wMax = 0
wrapThreshold = 0.0
sWrap = 0

if(HasValue(l) == 1) {
lval = l
}

if(HasValue(c1) == 1) {
c1val = c1
}

if(HasValue(c2) == 1) {
c2val = c2
}

if(HasValue(gain) == 1) {
g = gain
}

if(HasValue(wrapMax) == 1) {
wMax = wrapMax
}

if(HasValue(wrapTh) == 1) {
wrapThreshold = wrapTh
}

if(HasValue(smartWrap) == 1) {
sWrap = smartWrap
}

# if the smart wrap option is enabled, force the wrapMax option to enabled as well
if(sWrap == 1) {
wMax = 1
}

# this is the actual formula
outVal = ((nir - red) / (nir + c1val * red - c2val * blue + lval))
if (g != 1.0) {
outVal = g * outVal
}

if((HasValue(clipMin) == 1) || (HasValue(clipMax) == 1)) {
cMin = -1.0
if (HasValue(clipMin) == 1) {
cMin = clipMin
} else {
cMin = min(outVal)
}

cMax = 1.0
if (HasValue(clipMax) == 1) {
cMax = clipMax
} else {
cMax = max(outVal)
}

outVal = clipdata(outVal, min = clipMin, max = clipMax)
}

if (wMax == 1) {
xMax = dim(outVal)[1,1,1]
yMax = dim(outVal)[2,1,1]
maxVal = max(outVal) - wrapThreshold
minVal = min(outVal)

if (sWrap == 1) {
maxVal = findAggMinMax(outVal, backlash = 5, threshold = 0.001)[2,1,1]
# that function returns values in a range of 0 - 1, so convert back to whatever range is being used here
cRange = maxVal - minVal
maxVal = (maxVal * cRange) - (cRange / 2)
}

xCurrent = 1
yCurrent = 1

while (yCurrent <= yMax) {
xCurrent = 1
while (xCurrent <= xMax) {
#printf("Debug: current position: x = %i, y = %i.\n", xCurrent, yCurrent)
if (outVal[xCurrent, yCurrent, 1] >= maxVal) {
outVal[xCurrent, yCurrent, 1] = minVal
}
xCurrent++
}
yCurrent++
}
}

return(outVal)

}

define arvi(){
if (\$ARGC != 3) {
printf("\nGenerates an Atmospherically-Resistant Vegetation Index map based on the input data.\n")
printf("\nNote: requires near-infrared, red, and blue spectral band data.\n")
printf("\n")
printf("\$1        = the near-infrared band of the input image (required)\n")
printf("\$2          = the red band of the input image (required)\n")
printf("\$3          = the blue band of the input image (required)\n")
printf("\n")
printf("Uses the formula as published in \"Above-Ground Biomass Estimation of Successional and Mature Forests Using TM Images in the Amazon Basin\", ")
printf("by Dengsheng Lu, Paul Mausel, Eduardo Brondizio, and Emilio Moran.\n")
# note: http://www.isprs.org/proceedings/XXXIV/part4/pdfpapers/059.pdf
return(null)
}

nir = \$1
red = \$2
blue = \$3

numPlanes = dim(nir)[3,1,1]
if (numPlanes != 1) {
printf("Error: the near-infrared dataset must have exactly one plane. The input value has %i planes.\n", numPlanes)
return(null)
}

numPlanes = dim(red)[3,1,1]
if (numPlanes != 1) {
printf("Error: the red dataset must have exactly one plane. The input value has %i planes.\n", numPlanes)
return(null)
}

numPlanes = dim(blue)[3,1,1]
if (numPlanes != 1) {
printf("Error: the blue dataset must have exactly one plane. The input value has %i planes.\n", numPlanes)
return(null)
}

printf("Generating an Atmospherically-Resistant Vegetation Index image.\n")

outVal = (nir - (red * 2.0) + blue) / (nir + (red * 2.0) + blue)

return(outVal)
}

define asvi(){
if (\$ARGC != 3) {
printf("\nGenerates an Atmospheric and Soil Vegetation Index map based on the input data.\n")
printf("\nNote: requires near-infrared, red, and blue spectral band data.\n")
printf("\n")
printf("\$1        = the near-infrared band of the input image (required)\n")
printf("\$2          = the red band of the input image (required)\n")
printf("\$3          = the blue band of the input image (required)\n")
printf("\n")
printf("Uses the formula as published in \"Above-Ground Biomass Estimation of Successional and Mature Forests Using TM Images in the Amazon Basin\", ")
printf("by Dengsheng Lu, Paul Mausel, Eduardo Brondizio, and Emilio Moran.\n")
# note: http://www.isprs.org/proceedings/XXXIV/part4/pdfpapers/059.pdf
return(null)
}

nir = \$1
red = \$2
blue = \$3

printf("Generating an Atmospheric and Soil Vegetation Index image.\n")

numPlanes = dim(nir)[3,1,1]
if (numPlanes != 1) {
printf("Error: the near-infrared dataset must have exactly one plane. The input value has %i planes.\n", numPlanes)
return(null)
}

numPlanes = dim(red)[3,1,1]
if (numPlanes != 1) {
printf("Error: the red dataset must have exactly one plane. The input value has %i planes.\n", numPlanes)
return(null)
}

numPlanes = dim(blue)[3,1,1]
if (numPlanes != 1) {
printf("Error: the blue dataset must have exactly one plane. The input value has %i planes.\n", numPlanes)
return(null)
}

outVal = nir + 0.5 - (0.5 * ( sqrt(pow((nir * 2.0 + 1), 2) - 8 * (nir - (red * 2.0) + blue))))

return(outVal)
}

define gemi(){
if (\$ARGC != 2) {
printf("\nGenerates a Global Environmental Monitoring Index map based on the input data.\n")
printf("\nNote: requires near-infrared and red spectral band data.\n")
printf("\n")
printf("\$1        = the near-infrared band of the input image (required)\n")
printf("\$2          = the red band of the input image (required)\n")
printf("\n")
printf("Uses the formula as published in \"Above-Ground Biomass Estimation of Successional and Mature Forests Using TM Images in the Amazon Basin\", ")
printf("by Dengsheng Lu, Paul Mausel, Eduardo Brondizio, and Emilio Moran.\n")
# note: http://www.isprs.org/proceedings/XXXIV/part4/pdfpapers/059.pdf
return(null)
}

nir = \$1
red = \$2

printf("Generating a Global Environmental Monitoring Index image.\n")

numPlanes = dim(nir)[3,1,1]
if (numPlanes != 1) {
printf("Error: the near-infrared dataset must have exactly one plane. The input value has %i planes.\n", numPlanes)
return(null)
}

numPlanes = dim(red)[3,1,1]
if (numPlanes != 1) {
printf("Error: the red dataset must have exactly one plane. The input value has %i planes.\n", numPlanes)
return(null)
}

# the formula as printed in the paper uses an unusual glyph to represent a sub-component of the equation which is used twice
# it turns out the glyph is the Greek letter "Xi"
squigglyLine = ( (pow(nir, 2) - pow(red, 2)) * 2.0 + (nir * 1.5) + (red * 0.5) ) / (nir + red + 0.5)

invert = nir
invert[0,0,0] = 1.0

outVal = squigglyLine * (invert - (squigglyLine * 0.25)) - ( (red - 0.125)  / (invert - red) )

return(outVal)
}

# note: according to the "Vegetation Detection for Mobile Robot Navigation" paper referenced below, the original
# (IE "non-modified") SAVI formula only differs from MSAVI in that it has a manually-specified variable called "L"
# which "is selected based on where isovegetation lines intersect the soil line", and that MSAVI and MSAVI2
# automatically calculate this value.
# Therefore, I have not included the original SAVI formula as a function.

define msavi(){
if (\$ARGC != 2) {
printf("\nGenerates a Modified Soil-Adjusted Vegetation Index map based on the input data.\n")
printf("\nNote: requires near-infrared and red spectral band data.\n")
printf("\n")
printf("\$1        = the near-infrared band of the input image (required)\n")
printf("\$2          = the red band of the input image (required)\n")
printf("\n")
printf("Uses the formula as published in \"Above-Ground Biomass Estimation of Successional and Mature Forests Using TM Images in the Amazon Basin\", ")
printf("by Dengsheng Lu, Paul Mausel, Eduardo Brondizio, and Emilio Moran.\n")
# note: http://www.isprs.org/proceedings/XXXIV/part4/pdfpapers/059.pdf
return(null)
}

nir = \$1
red = \$2

printf("Generating a Modified Soil-Adjusted Vegetation Index image.\n")

numPlanes = dim(nir)[3,1,1]
if (numPlanes != 1) {
printf("Error: the near-infrared dataset must have exactly one plane. The input value has %i planes.\n", numPlanes)
return(null)
}

numPlanes = dim(red)[3,1,1]
if (numPlanes != 1) {
printf("Error: the red dataset must have exactly one plane. The input value has %i planes.\n", numPlanes)
return(null)
}

outVal = nir + 0.5 - ((sqrt( pow((nir * 2.0 + 1), 2) - ((nir - (red * 2.0)) * 8.0) )) * 0.5)

return(outVal)
}

define msavi2(){
if (\$ARGC != 2) {
printf("\nGenerates a Modified Soil-Adjusted Vegetation Index 2 map based on the input data.\n")
printf("\nNote: requires near-infrared and red spectral band data.\n")
printf("\n")
printf("\$1    = the near-infrared band of the input image (required)\n")
printf("\$2    = the red band of the input image (required)\n")
printf("\n")
printf("Uses the formula as published in \"Vegetation Detection for Mobile Robot Navigation\", ")
printf("by David M. Bradley, Scott M. Thayer, Anthony Stentz, and Peter Rander.\n")
printf("\n")
printf("\nNote: Useful output is in the range of 0.0 - 1.0.\n")
return(null)
}

nir = \$1
red = \$2

printf("Generating a Modified Soil-Adjusted Vegetation Index 2 image.\n")

numPlanes = dim(nir)[3,1,1]
if (numPlanes != 1) {
printf("Error: the near-infrared dataset must have exactly one plane. The input value has %i planes.\n", numPlanes)
return(null)
}

numPlanes = dim(red)[3,1,1]
if (numPlanes != 1) {
printf("Error: the red dataset must have exactly one plane. The input value has %i planes.\n", numPlanes)
return(null)
}

outVal = (((nir + 1.0) * 2.0) - sqrt( pow(((nir * 2.0) + 1.0), 2) - ((nir - red) * 8.0))) / 2

return(outVal)
}

define avi(mwnir, mwred, mwgreen){
if (\$ARGC != 3) {
printf("\nGenerates an Angular Vegetation Index map based on the input data.\n")
printf("\nNote: requires near-infrared, red, and green spectral band data.\n")
printf("\n")
printf("\$1    = the near-infrared band of the input image (required)\n")
printf("\$2    = the red band of the input image (required)\n")
printf("\$3    = the green band of the input image (required)\n")
printf("mwnir    = the middle wavelength of the near-infrared band (default = 825)\n")
printf("mwred    = the middle wavelength of the red band (default = 660\n")
printf("mwgreen    = the middle wavelength of the green band (default = 565)\n")
printf("\n")
printf("Uses the formula as published in \"Bilko Mini-Lesson 8: Vegetation Indices\"\n")
# note: http://www.noc.soton.ac.uk/bilko/miniless/l8_100.php
printf("\n")
printf("Note: the units used for the optional middle wavelength values are unimportant as long as the same unit is used for all three values.\n")
return(null)
}

nir = \$1
red = \$2
green = \$3

printf("Generating an Angular Vegetation Index image.\n")

numPlanes = dim(nir)[3,1,1]
if (numPlanes != 1) {
printf("Error: the near-infrared dataset must have exactly one plane. The input value has %i planes.\n", numPlanes)
return(null)
}

numPlanes = dim(red)[3,1,1]
if (numPlanes != 1) {
printf("Error: the red dataset must have exactly one plane. The input value has %i planes.\n", numPlanes)
return(null)
}

numPlanes = dim(green)[3,1,1]
if (numPlanes != 1) {
printf("Error: the green dataset must have exactly one plane. The input value has %i planes.\n", numPlanes)
return(null)
}

mnir = 825.0
mred = 660.0
mgreen = 565.0

if(HasValue(mwnir) == 1) {
mnir = mwnir
}

if(HasValue(mwred) == 1) {
mred = mwred
}

if(HasValue(mwgreen) == 1) {
mgreen = mwgreen
}

recip = nir
recip[0,0,0] = 1.0

outVal = atan( ((mnir - mred)/mred) * (recip / (nir - red)) ) + atan( ((mred - mgreen)/mred) * (recip / (green - red)))
return(outVal)
}

# note: the calculations for the KT1 - KT6 images are all very simple. I have created functions for them in order
# to avoid having to memorize a bunch of coefficients or treat them as magic numbers elsewhere.

define kthelp(){
printf("\nNote: requires six spectral bands, which correspond to channels 1-5 and 7 of Landsat satellite data.\n")
printf("\n")
printf("\$1    = blue band data (TM1 / 450nm - 515nm) (required)\n")
printf("\$2      = green band data (TM2 / 525nm - 605nm) (required)\n")
printf("\$3      = red band data (TM3/ 630nm - 690nm) (required)\n")
printf("\$4      = near-infrared band data (TM4 / 750nm    - 900nm) (required)\n")
printf("\$5      = mid-infrared band 1 data (TM5 / 1500nm - 1750nm) (required)\n")
printf("\$6      = mid-infrared band 2 data (TM7 / 2090nm - 2350nm (required)\n")
printf("mode    = processing mode (default = 0)\n")
printf("      0 = LandSat 4 reflectance values\n")
#printf("      1 = LandSat 5 reflectance values\n")
printf("      2 = LandSat 7 reflectance values\n")
printf("      3 = LandSat 4 DN values\n")
printf("      4 = LandSat 5 DN values\n")
#printf("      5 = LandSat 7 DN values\n")
printf("\n")
printf("Note: TM6 (thermal/long-wave IR) data is not used in this calculation.\n")
printf("\n")
printf("Uses the formula as published in \"Above-Ground Biomass Estimation of Successional and Mature Forests Using TM Images in the Amazon Basin\", ")
printf("(Dengsheng Lu, Paul Mausel, Eduardo Brondizio, and Emilio Moran, 2002),\n")
# note: http://www.isprs.org/proceedings/XXXIV/part4/pdfpapers/059.pdf
printf("\"Tasseled Cap Transformation: Mississippi Coastal Corridor July 24, 2000\",")
printf("(Dr. Roger King, Dr. Charles O'Hara, 2001),\n")
# note: http://www.ncrste.msstate.edu/archive/publications/posters/ncrste_tg002-9.pdf
printf("and \"The Tasseled Cap Transformation in Remote Sensing\",")
printf("(Thayer Watkins, date unknown).\n")
# note: http://www.sjsu.edu/faculty/watkins/tassel.htm
}

define kt1(mode){
if (\$ARGC != 6) {
printf("\nGenerates a KT1 (Tasseled Cap - Brightness) image map based on the input data.\n")
kthelp()
return(null)
}

tm1 = \$1
tm2 = \$2
tm3 = \$3
tm4 = \$4
tm5 = \$5
tm7 = \$6

calcMode = 0

if(HasValue(mode) == 1) {
calcMode = mode
}

printf("Generating a KT1 (Tasseled Cap - Brightness) image.\n")

#outVal = (tm1 * 0.304) + (tm2 * 0.279) + (tm3 * 0.474) + (tm4 * 0.559) + (tm5 * 0.508) + (tm7 * 0.186)
#return(outVal)

# LandSat 4 reflectance values
if(calcMode == 0) {
outVal = (tm1 * 0.2043) + (tm2 * 0.4158) + (tm3 * 0.5524) + (tm4 * 0.5741) + (tm5 * 0.3124) + (tm7 * 0.2303)
return(outVal)
}

# LandSat 7 reflectance values
if(calcMode == 2) {
outVal = (tm1 * 0.3561) + (tm2 * 0.3972) + (tm3 * 0.3904) + (tm4 * 0.6966) + (tm5 * 0.2286) + (tm7 * 0.1596)
return(outVal)
}

# LandSat 4 DN values
if(calcMode == 3) {
outVal = (tm1 * 0.3037) + (tm2 * 0.2793) + (tm3 * 0.4743) + (tm4 * 0.5585) + (tm5 * 0.5082) + (tm7 * 0.1863)
return(outVal)
}

# LandSat 5 DN values
if(calcMode == 4) {
outVal = (tm1 * 0.2909) + (tm2 * 0.2493) + (tm3 * 0.4806) + (tm4 * 0.5568) + (tm5 * 0.4438) + (tm7 * 0.1706)
return(outVal)
}

printf("Error: invalid KT mode specified.\n")
return(null)
}

define kt2(mode){
if (\$ARGC != 6) {
printf("\nGenerates a KT2 (Tasseled Cap - Greenness) image map based on the input data.\n")
kthelp()
return(null)
}

tm1 = \$1
tm2 = \$2
tm3 = \$3
tm4 = \$4
tm5 = \$5
tm7 = \$6

calcMode = 0

if(HasValue(mode) == 1) {
calcMode = mode
}

printf("Generating a KT2 (Tasseled Cap - Greenness) image.\n")

#outVal = (tm1 * (-0.285)) - (tm2 * 0.244) - (tm3 * 0.544) + (tm4 * 0.704) + (tm5 * 0.084) - (tm7 * 0.18)
#return(outVal)

# LandSat 4 reflectance values
if(calcMode == 0) {
outVal = (tm1 * -0.1063) + (tm2 * -0.2819) + (tm3 * -0.4934) + (tm4 * 0.794) + (tm5 * -0.0002) + (tm7 * -0.1446)
return(outVal)
}

# LandSat 7 reflectance values
if(calcMode == 2) {
outVal = (tm1 * -0.3344) + (tm2 * -0.3544) + (tm3 * -0.4556) + (tm4 * 0.6966) + (tm5 * -0.0242) + (tm7 * -0.263)
return(outVal)
}

# LandSat 4 DN values
if(calcMode == 3) {
outVal = (tm1 * -0.2848) + (tm2 * -0.2435) + (tm3 * -0.5436) + (tm4 * 0.7243) + (tm5 * 0.084) + (tm7 * -0.18)
return(outVal)
}

# LandSat 5 DN values
if(calcMode == 4) {
outVal = (tm1 * -0.2728) + (tm2 * -0.2174) + (tm3 * -0.5508) + (tm4 * 0.722) + (tm5 * 0.0733) + (tm7 * -0.1648)
return(outVal)
}

printf("Error: invalid KT mode specified.\n")
return(null)
}

define kt3(mode){
if (\$ARGC != 6) {
printf("\nGenerates a KT3 (Tasseled Cap - Wetness) image map based on the input data.\n")
kthelp()
return(null)
}

tm1 = \$1
tm2 = \$2
tm3 = \$3
tm4 = \$4
tm5 = \$5
tm7 = \$6

calcMode = 0

if(HasValue(mode) == 1) {
calcMode = mode
}

printf("Generating a KT3 (Tasseled Cap - Wetness) image.\n")

#outVal = (tm1 * 0.151) + (tm2 * 0.197) + (tm3 * 0.328) + (tm4 * 0.341) - (tm5 * 0.711) - (tm7 * 0.457)
#return(outVal)

# LandSat 4 reflectance values
if(calcMode == 0) {
outVal = (tm1 * 0.0315) + (tm2 * 0.2021) + (tm3 * 0.3102) + (tm4 * 0.1594) + (tm5 * -0.6806) + (tm7 * -0.6109)
return(outVal)
}

# LandSat 7 reflectance values
if(calcMode == 2) {
outVal = (tm1 * 0.2626) + (tm2 * 0.2141) + (tm3 * 0.0926) + (tm4 * 0.0656) + (tm5 * -0.7629) + (tm7 * -0.5388)
return(outVal)
}

# LandSat 4 DN values
if(calcMode == 3) {
outVal = (tm1 * 0.1509) + (tm2 * 0.1973) + (tm3 * 0.3279) + (tm4 * 0.3406) + (tm5 * -0.7112) + (tm7 * -0.4572)
return(outVal)
}

# LandSat 5 DN values
if(calcMode == 4) {
outVal = (tm1 * 0.1446) + (tm2 * 0.1761) + (tm3 * 0.3322) + (tm4 * 0.3396) + (tm5 * -0.621) + (tm7 * -0.4186)
return(outVal)
}

printf("Error: invalid KT mode specified.\n")
return(null)
}

define kt4(mode){
if (\$ARGC != 6) {
printf("\nGenerates a KT4 (Tasseled Cap - Fourth Component) image map based on the input data.\n")
kthelp()
return(null)
}

tm1 = \$1
tm2 = \$2
tm3 = \$3
tm4 = \$4
tm5 = \$5
tm7 = \$6

calcMode = 0

if(HasValue(mode) == 1) {
calcMode = mode
}

printf("Generating a KT4 (Tasseled Cap - Fourth Component) image.\n")

#outVal = (tm1 * 0.151) + (tm2 * 0.197) + (tm3 * 0.328) + (tm4 * 0.341) - (tm5 * 0.711) - (tm7 * 0.457)
#return(outVal)

# LandSat 4 reflectance values
if(calcMode == 0) {
outVal = (tm1 * -0.2117) + (tm2 * -0.0284) + (tm3 * 0.1302) + (tm4 * -0.1007) + (tm5 * 0.6529) + (tm7 * -0.7078)
return(outVal)
}

# LandSat 7 reflectance values
if(calcMode == 2) {
outVal = (tm1 * 0.0805) + (tm2 * -0.0498) + (tm3 * 0.195) + (tm4 * -0.1327) + (tm5 * 0.5752) + (tm7 * -0.7775)
return(outVal)
}

# LandSat 4 DN values
if(calcMode == 3) {
outVal = (tm1 * 0.8832) + (tm2 * -0.0819) + (tm3 * -0.458) + (tm4 * -0.0032) + (tm5 * -0.0563) + (tm7 * 0.013)
return(outVal)
}

# LandSat 5 DN values
if(calcMode == 4) {
outVal = (tm1 * 0.8461) + (tm2 * -0.0731) + (tm3 * -0.464) + (tm4 * -0.0032) + (tm5 * -0.0492) + (tm7 * 0.0119)
return(outVal)
}

printf("Error: invalid KT mode specified.\n")
return(null)
}

define kt5(mode){
if (\$ARGC != 6) {
printf("\nGenerates a KT5 (Tasseled Cap - Fifth Component) image map based on the input data.\n")
kthelp()
return(null)
}

tm1 = \$1
tm2 = \$2
tm3 = \$3
tm4 = \$4
tm5 = \$5
tm7 = \$6

calcMode = 0

if(HasValue(mode) == 1) {
calcMode = mode
}

printf("Generating a KT5 (Tasseled Cap - Fifth Component) image.\n")

#outVal = (tm1 * 0.151) + (tm2 * 0.197) + (tm3 * 0.328) + (tm4 * 0.341) - (tm5 * 0.711) - (tm7 * 0.457)
#return(outVal)

# LandSat 4 reflectance values
if(calcMode == 0) {
outVal = (tm1 * -0.8669) + (tm2 * -0.1835) + (tm3 * 0.3856) + (tm4 * 0.0408) + (tm5 * -0.1132) + (tm7 * 0.2272)
return(outVal)
}

# LandSat 7 reflectance values
if(calcMode == 2) {
outVal = (tm1 * -0.7252) + (tm2 * -0.0202) + (tm3 * 0.6683) + (tm4 * 0.0631) + (tm5 * -0.1494) + (tm7 * -0.0274)
return(outVal)
}

# LandSat 4 DN values
if(calcMode == 3) {
outVal = (tm1 * 0.0573) + (tm2 * -0.026) + (tm3 * 0.0335) + (tm4 * -0.1943) + (tm5 * 0.4766) + (tm7 * -0.8545)
return(outVal)
}

# LandSat 5 DN values
if(calcMode == 4) {
outVal = (tm1 * 0.0549) + (tm2 * -0.0232) + (tm3 * 0.0339) + (tm4 * -0.1937) + (tm5 * 0.4162) + (tm7 * -0.7823)
return(outVal)
}

printf("Error: invalid KT mode specified.\n")
return(null)
}

define kt6(mode){
if (\$ARGC != 6) {
printf("\nGenerates a KT6 (Sixth Component) image map based on the input data.\n")
kthelp()
return(null)
}

tm1 = \$1
tm2 = \$2
tm3 = \$3
tm4 = \$4
tm5 = \$5
tm7 = \$6

calcMode = 0

if(HasValue(mode) == 1) {
calcMode = mode
}

printf("Generating a KT6 (Sixth Component) image.\n")

#outVal = (tm1 * 0.151) + (tm2 * 0.197) + (tm3 * 0.328) + (tm4 * 0.341) - (tm5 * 0.711) - (tm7 * 0.457)
#return(outVal)

# LandSat 4 reflectance values
if(calcMode == 0) {
outVal = (tm1 * 0.3677) + (tm2 * 0.82) + (tm3 * 0.4354) + (tm4 * 0.0518) + (tm5 * -0.0066) + (tm7 * -0.0104)
return(outVal)
}

# LandSat 7 reflectance values
if(calcMode == 2) {
outVal = (tm1 * 0.4) + (tm2 * -0.8172) + (tm3 * 0.3832) + (tm4 * 0.0602) + (tm5 * -0.1095) + (tm7 * 0.0985)
return(outVal)
}

# LandSat 4 DN values
if(calcMode == 3) {
outVal = (tm1 * 0.1238) + (tm2 * -0.9038) + (tm3 * 0.4041) + (tm4 * 0.0573) + (tm5 * -0.0261) + (tm7 * 0.024)
return(outVal)
}

# LandSat 5 DN values
if(calcMode == 4) {
outVal = (tm1 * 0.1186) + (tm2 * -0.8069) + (tm3 * 0.4094) + (tm4 * 0.0571) + (tm5 * -0.0228) + (tm7 * 0.22)
return(outVal)
}

printf("Error: invalid KT mode specified.\n")
return(null)
}

define slavi(){
if (\$ARGC != 3) {
printf("\nGenerates a Specific Leaf Area Vegetation Index image map based on the input data.\n")
printf("\nNote: requires red, near-infrared, and mid-infrared spectral band data.\n")
printf("\n")
printf("\$1      = mid-infrared data (TM5 / 1500nm - 1750nm or TM7 / 2090nm - 2350nm) (required)\n")
printf("\$2      = near-infrared band data (TM4 / 750nm    - 900nm) (required)\n")
printf("\$3    = red band data (required)\n")
printf("\n")
printf("Note: the traditional version of this calculation uses landsat band 7 (2090nm - 2350nm) for the mid-infrared data.\n")
printf("\n")
printf("Uses the formula as published in \"Landsat Image Processing: ES 351, ES 771, ES 775\", ")
printf("by James S. Aber.\n")
return(null)
}

mir = \$1
nir = \$2
red = \$3

printf("Generating a Specific Leaf Area Vegetation Index image.\n")

outVal = nir / (red + mir)

return(outVal)
}

define findAggMinMax(steps, backlash, threshold, threshrepeat, searchwidth){
if (\$ARGC != 1) {
printf("\nFinds new minimum and maximum values for an input set - for use in aggressive normalization.\n")
printf("\nReturns an array of size [2,1,1] containing the calculated minimum and maximum.\n")
printf("\n")
printf("\$1           = the input array\n")
printf("steps          = the number of steps between 0 and 1 to use when generating the histogram (default = 100)\n")
printf("backlash     = the number of steps to \"back off\" the new min and max values when rescaling (default = 1)\n")
printf("threshold    = the relative count that must be exceeded by a step in the histogram to trigger a response (default = 0.05)\n")
printf("threshrepeat = the number of steps that must exceed the threshold in order to trigger a response (default = 3)\n")
printf("searchwidth  = the number of steps after the first which exceeds the threshold that will be searched to find a repetition (default = 5)\n")
return(null)
}

stepCount = 100
bl = 1
thresh = 0.05
trepeat = 3
swidth = 5

if(HasValue(steps) == 1) {
stepCount = steps
}

if(HasValue(backlash) == 1) {
bl = backlash
}

if(HasValue(threshold) == 1) {
thresh = threshold
}

if(HasValue(threshrepeat) == 1) {
trepeat = threshrepeat
}

if(HasValue(searchwidth) == 1) {
swidth = searchwidth
}

printf("Attempting to find aggressive normalization parameters.\n")

# first, perform a standard normalization to make the rest of the calculations easier
inVal = float(\$1)

inMin = min(inVal)
inMax = max(inVal)

if ((inMin != 0.0) || (inMax != 1.0)) {
inVal = inVal - inMin
inVal = inVal / max(inVal)
}

stepSize = (1.0 / stepCount)

hist = histogram(inVal, start = 0.0, size = stepSize, steps = stepCount)

newMin = 0.0
newMax = 1.0

maxCount = max(hist[2,0,0])

tAbs = maxCount * thresh

currentStep = 1

foundInitial = 0
foundMatch = 0
tNewVal = newMin
windowIteration = 0
hitCount = 0

#printf("Debug: tAbs = %i\n", tAbs)

#printf("Debug: finding new minimum.\n")

while (currentStep <= stepCount) {
#printf("Debug: currentStep = %i\n", currentStep)
if (hist[2,currentStep,1] >= tAbs) {
if (foundInitial == 0) {
#printf("Debug: found initial hit at step %i\n", currentStep)
tNewVal = hist[1,currentStep,1]
foundInitial = 1
windowIteration = 0
hitCount = 1
} else {
#printf("Debug: found secondary hit at step %i\n", currentStep)
hitCount++
if (hitCount >= trepeat) {
foundMatch = 1
newMin = tNewVal
break
}
}
} else {
if (foundInitial == 1) {
windowIteration++
if (windowIteration > swidth) {
#printf("Debug: gave up on current hit at %i\n", currentStep)
foundInitial = 0
}
}
}
currentStep++
}

#printf("Debug: new minimum: %f.\n", newMin)
#printf("Debug: finding new maximum.\n")
currentStep = stepCount
foundInitial = 0
hitCount = 0
while (currentStep >= 1) {
#printf("Debug: currentStep = %i\n", currentStep)
if (hist[2,currentStep,1] >= tAbs) {
if (foundInitial == 0) {
#printf("Debug: found initial hit at step %i - value = %f\n", currentStep, hist[2,currentStep,1])
tNewVal = hist[1,currentStep,1]
foundInitial = 1
windowIteration = 0
hitCount = 1
} else {
#printf("Debug: found secondary hit at step %i - value = %f\n", currentStep, hist[2,currentStep,1])
hitCount++
if (hitCount >= trepeat) {
foundMatch = 1
newMax = tNewVal
break
}
}
} else {
if (foundInitial == 1) {
windowIteration++
if (windowIteration > swidth) {
#printf("Debug: gave up on current hit at %i\n", currentStep)
foundInitial = 0
}
}
}
currentStep--
}

#printf("Debug: new maximum: %f.\n", newMax)

if ((newMax <= newMin) || (newMax == 0.0) || (newMin == 1.0)) {
printf("\nError: Unable to obtain new minimum and maximum values using the specified search parameters. The result will only have standard normalization applied.\n")
return(inVal)
}

if (newMin > 0) {
if (newMin > (stepSize * bl)) {
newMin = newMin - (stepSize * bl)
} else {
newMin = 0.0
}
}
#if (newMin < 0.0) {
#    newMin = 0.0
#}
#apparently DaVinci is not always so great at determining if a value is less than zero
#if (newMin != abs(newMin)) {
#    newMin = 0.0
#}

if (newMax < 1.0) {
if ((newMax + (stepSize * bl)) < 1.0) {
newMax = newMax + (stepSize * bl)
} else {
newMax = 1.0
}
}
#newMax = newMax + (stepSize * bl)
#if (newMax > 1.0) {
#    newMax = 1.0
#}

printf("New minimum: %f.\n", newMin)
printf("New maximum: %f.\n", newMax)

outVal = clone(float(0.0),2,1,1)
outVal[1,1,1] = newMin
outVal[2,1,1] = newMax

return(outVal)
}

# because clock() errors out for me and I didn't need anything that complicated
# most of this code is copy/pasted from clock()
define getdatetime(){
istr = "Error obtaining date/time"
ios = getos()
if (equals(ios, "win") == 1) { # Windows
istr = syscall("echo %date% %time%")
} else if (equals(ios, "mac") == 1) { # Mac
istr = syscall("date \"+%Y-%m-%d %H:%M:%S\"")
} else { # Does this work for everything that isn't Windows or Mac? Or just Linux?
istr = syscall("date \"+%F %T.%N %z\"")
}
return(istr)
}

define echodatetime(){
dt = getdatetime()
printf("Local date/time: %s\n", dt[1:])
}

goodsyntax = 1
# minimum of four arguments (position 1, colour 1, position 2, colour 2)
if (\$ARGC < 4) {
goodsyntax = 0
}
# number of arguments must be even
if ((\$ARGC % 2) == 1) {
goodsyntax = 0
}
if (goodsyntax == 0) {
printf("\nBuilds a 1-dimensional colour gradient\n")
printf("\n")
printf("width        = the horizontal size of the gradient (default = 65536)\n")
printf("smoothwindow    = width in pixels to smooth gradient transitions (default = 1)\n")
printf("\$1        = the first colour position (integer, 1 <= n <= width)\n")
printf("\$2        = the first colour\n")
printf("\$3        = the second colour position (integer, 1 <= n <= width)\n")
printf("\$4        = the second colour\n")
printf("\n")
printf("Additional colours may be specified by adding more arguments. The number of arguments ")
printf("(not including the optional width and smoothwind parameters) must be even (because each colour requires a position).\n")
printf("\n")
printf("Note: the number of channels (greyscale, RGB, ARGB) is determined by the Z-dimension of the first colour.\n")
printf("\n")
return(null)
}

gWidth = 65536
if(HasValue(width) == 1) {
gWidth = width
}

smoothWin = 1
if(HasValue(smoothwindow) == 1) {
smoothWin = smoothwindow
}

colourCount = \$ARGC / 2
channelCount = dim(\$2)[3,1,1]

printf("Generating a gradient using %i colours and %i channels.\n", colourCount, channelCount)

#printf("Debug: colour count = %i.\n", colourCount)
#printf("Debug: channel count = %i.\n", channelCount)

# build an empty (black) gradient
gradient = clone(\$2, gWidth, 1, 1)

currentColourNum = 1

while (currentColourNum < colourCount) {
#printf("Debug: current colour number = %i.\n", currentColourNum)

ccIndex = (currentColourNum * 2) - 1
#printf("Debug: \$ARGV array index = %i.\n", ccIndex)

#printf("Debug: getting beginning position.\n")
beginPosition = \$ARGV[ccIndex]

#printf("Debug: getting ending position.\n")
endPosition = \$ARGV[ccIndex + 2]

#printf("Debug: getting beginning colour.\n")
beginColour = \$ARGV[ccIndex + 1]

#printf("Debug: getting ending colour.\n")
endColour = \$ARGV[ccIndex + 3]

#printf("Debug: done accessing \$ARGV.\n")

stepCount = (endPosition - beginPosition) + 1
#printf("Debug: step count = %i.\n", stepCount)

if (stepCount <= 0) {
printf("Error: colour number %i of the gradient has a position that is less than or equal to the next colour.\n", currentColourNum)
return(null)
}

#printf("Debug: start position for this range = %i.\n", beginPosition)
#printf("Debug: end position for this range = %i.\n", endPosition)
#printf("Debug: step count for this range = %i.\n", stepCount)

inc = (float(endColour - beginColour)) / float(stepCount)

#printf("Debug: beginning colour:\n")
#echoarraycontents(beginColour)
#printf("Debug: ending colour:\n")
#echoarraycontents(endColour)
#printf("Debug: colour increment per step:\n")
#echoarraycontents(inc)

stepNum = 1
basePosition = beginPosition - 1
while (stepNum <= stepCount) {
#channelNum = 1
#while (channelNum <= channelCount) {
#    channelVal = beginColour[1, 1, channelNum]
#    channelVal += (inc[1, 1, channelNum] * stepNum)
#    gradient[basePosition + stepNum, 1, channelNum] = channelVal
#    channelNum++
#}

channelVal = beginColour[1, 1, 0]
channelVal += (inc[1, 1, 0] * stepNum)
gradient[basePosition + stepNum, 1, 0] = channelVal

stepNum++
}

if (smoothWin > 1) {
# perform a sliding smooth operation on the gradient section that was just generated
#printf("Debug: smoothing transition %i.\n", currentColourNum)

#smoothSteps = ceil((endPosition - beginPosition) / smoothWin)
#printf("Debug: smoothing steps needed: %i.\n", smoothSteps)

stepNum = 1
basePosition = beginPosition - 1
while (stepNum <= stepCount) {
currentPosition = basePosition + stepNum

#sectionWidth = float(endPosition - beginPosition)
#perc = 1.0 - abs(((float(currentPosition) - float(beginPosition)) / sectionWidth) - 0.5)
#perc = abs(((float(currentPosition) - float(beginPosition)) / sectionWidth) - 0.5)
#acWinSize = ceil(float(smoothWin) * 0.5) * pow(perc, 2)

#perc = (abs((float(stepNum) / float(stepCount)) - 0.5)) * 2.0
perc = 1.0 - ((abs((float(stepNum) / float(stepCount)) - 0.5)) * 2.0)

#printf("Debug: perc: %f.\n", perc)

acWinSize = ceil((float(smoothWin) * 0.5) * pow(perc, 3)) - 1

if (acWinSize < 0) {
acWinSize = 0
}

winBegin = int(currentPosition - acWinSize)
if (winBegin < beginPosition) {
winBegin = beginPosition
}

winEnd = int(winBegin + acWinSize)
if (winBegin > endPosition) {
winEnd = endPosition
}

#printf("Debug: winBegin: %i.\n", winBegin)
#printf("Debug: winEnd: %i.\n", winEnd)

performSmooth = 1

# do not attempt a smoothing operation of this section if...
# ...there is nothing to smooth between
if (winBegin == winEnd) {
performSmooth = 0
}

# ...there was a calculation error
if (winBegin > winEnd) {
performSmooth = 0
#printf("Debug: error: winBegin > winEnd.\n")
}

if (performSmooth == 1)
{
channelNum = 1
while (channelNum <= channelCount) {
channelVal = avg(gradient[winBegin:winEnd, 1, channelNum])
sGradient[basePosition + stepNum, 1, channelNum] = channelVal
channelNum++
}
}

stepNum++
}

}

currentColourNum++
}

}

goodsyntax = 1
# minimum of three arguments (colour 1 (upper left), colour 2 (upper right), colour 3 (lower left))
if (\$ARGC < 3) {
goodsyntax = 0
}

if (goodsyntax == 0) {
printf("\nBuilds a simple 2-dimensional colour gradient\n")
printf("\n")
printf("width        = the horizontal size of the gradient (default = 512)\n")
printf("height        = the vertical size of the gradient (default = 512)\n")
printf("\$1        = the first colour (will appear in the upper-left corner of the image)\n")
printf("\$2        = the second colour (will appear in the upper-right corner of the image)\n")
printf("\$3        = the third colour (will appear in the lower-left corner of the image)\n")
printf("\n")
printf("Note: the number of channels (greyscale, RGB, ARGB) is determined by the Z-dimension of the first colour.\n")
printf("\n")
return(null)
}

hWidth = 512
if(HasValue(width) == 1) {
hWidth = width
}

hHeight = 512
if(HasValue(height) == 1) {
hHeight = height
}

colour1 = \$1
colour2 = \$2
colour3 = \$3

printf("Generating a 2d gradient.\n")

# build an empty (black) gradient
gradient = clone(colour1, hWidth, hHeight, 1)

# build the y-axis gradient first
yaGradient = buildgradient(1, colour1, hHeight, colour3, width=int(hHeight))

yPosition = 1
while (yPosition <= hHeight) {
yPosition++
}

}

define echoarraycontents(){
arrayDim = dim(\$1)

xDim = arrayDim[1,1,1]
yDim = arrayDim[2,1,1]
zDim = arrayDim[3,1,1]

printf("Array size: X = %i, Y = %i, Z = %i.\n", xDim, yDim, zDim)

zCurrent = 1
while (zCurrent <= zDim) {
yCurrent = 1
while (yCurrent <= yDim) {
xCurrent = 1
while (xCurrent <= xDim) {
currentVal = \$1
currentVal = currentVal[xCurrent, yCurrent, zCurrent]
printf("%i, %i, %i = %f\n", xCurrent, yCurrent, zCurrent, currentVal)
xCurrent++
}
yCurrent++
}
zCurrent++
}
}

define applypalette(min, max){
if (\$ARGC != 2) {
printf("\nApplies a palette (such as a gradient) to a greyscale image.\n")
printf("\n")
printf("\$1    = the greyscale image (required)\n")
printf("\$2    = the 1-dimensional palette image to apply (required)\n")
printf("min    = the value in the image to map to the leftmost colour of the palette (default = 0.0).\n")
printf("max    = the value in the image to map to the rightmost colour of the palette (default = 1.0).\n")
printf("\n")
return(null)
}

printf("Applying the specified palette to the specified image.\n")

gsImage = \$1
palette = \$2

minVal = float(0.0)
maxVal = float(1.0)

if(HasValue(min) == 1) {
minVal = float(min)
}
if(HasValue(max) == 1) {
maxVal = float(max)
}

range = maxVal - minVal

pDim = dim(palette)
pMax = pDim[1,1,1]
#paletteMax = float(pMax) - 1.0
paletteMax = float(pMax)
#paletteStep = float(1.0) / float(paletteMax)

iDim = dim(gsImage)
xSize = iDim[1,1,1]
ySize = iDim[2,1,1]
zSize = pDim[3,1,1]

outImage = clone(0.0, xSize, ySize, zSize)

clippedCountMax = 0
clippedCountMin = 0

#printf("Debug: range = %f\n", range)
#printf("Debug: paletteMax = %f\n", paletteMax)
#printf("Debug: paletteStep = %f\n", paletteStep)

yCurrent = 1
while (yCurrent <= ySize) {
xCurrent = 1
while (xCurrent <= xSize) {
currentVal = float(gsImage[xCurrent, yCurrent, 1])

mappedPosition = int(((currentVal - minVal) / range) * paletteMax)

if (mappedPosition > paletteMax) {
clippedCountMax++
mappedPosition = int(paletteMax)
}

mappedPosition = int(mappedPosition) + 1

if (mappedPosition > paletteMax) {
mappedPosition = int(paletteMax)
}

if (mappedPosition < 1) {
clippedCountMin++
mappedPosition = 1
}

outImage[xCurrent, yCurrent, 0] = palette[mappedPosition, 1, 0]

xCurrent++
}
yCurrent++
}

if (clippedCountMax > 0) {
printf("Warning: %i values exceeded the maximum allowed value and were clipped to the maximum value as a result.\n", clippedCountMax)
}
if (clippedCountMin > 0) {
printf("Warning: %i values were lower than the minimum allowed value and were clipped to the minimum value as a result.\n", clippedCountMin)
}

return(outImage)
}

define iiaz(){
if (\$ARGC != 1) {
printf("\nInvert Independently About Zero\n")
printf("\nInverts positive and negative values independently based on the current minimum and maximum values.\n")
printf("For example, if the input values consist of -1.0, -0.3, 0.25, 0.9, and 1.0, the output values will consist of 0.0, -0.7, 0.75, 0.1, and 0.0.\n")
printf("\n")
printf("\$1    = the input value (required)\n")
return(null)
}

inVal = \$1

# handle positive values
posVals = clipdata(inVal, min = 0.0)
mask[0, 0, 0] = max(posVals)
posVals = mask - posVals

# handle negative values
negVals = clipdata(inVal, max = 0.0)
mask[0, 0, 0] = min(negVals)
negVals = mask - negVals

outVal = inVal

zMax = dim(outVal)[3, 1, 1]
zCurrent = 1

while (zCurrent <= zMax) {
maskTemp[0, 0, 2] = 1.0

outVal[0, 0, zCurrent] = posVals[0, 0, zCurrent] * maskTemp

maskTemp[0, 0, 2] = 1.0

outVal[0, 0, zCurrent] = outVal[0, 0, zCurrent] + (negVals[0, 0, zCurrent] * maskTemp)

zCurrent++
}

return(outVal)
}

define naav(center, ubound, lbound, independent, anormalize, steps, backlash, threshold, threshrepeat, searchwidth){
if (\$ARGC != 1) {
printf("\nNormalize About Arbitrary Value.\n")
printf("\nNormalizes positive values up and negative values down in relation to an arbitrary center-point.\n")
printf("\n")
printf("\$1        = the input value (required)\n")
printf("center        = the value which determines the cutoff between positive and negative (default = 0.0).\n")
printf("ubound        = the upper boundary for positive values (default = 1.0).\n")
printf("lbound        = the lower boundary for negative values (default = -1.0).\n")
printf("independent    = normalizes positive and negative values using separate scales (default = 0).\n")
printf("anormalize    = perform an aggressive normalization on the data (default = 0).\n")
printf("\n")
printf("Important: enabling aggressive normalization implicitly enables independent normalization.\n")
printf("\n")
printf("Additional options available when anormalize is set to 1:\n")
printf("steps              = the number of steps between 0 and 1 to use when generating the histogram (default = 100)\n")
printf("backlash         = the number of steps to \"back off\" the new min and max values when rescaling (default = 1)\n")
printf("threshold        = the relative count that must be exceeded by a step in the histogram to trigger a response (default = 0.05)\n")
printf("threshrepeat     = the number of steps that must exceed the threshold in order to trigger a response (default = 3)\n")
printf("searchwidth      = the number of steps after the first which exceeds the threshold that will be searched to find a repetition (default = 5)\n")
printf("\n")
printf("Note: if independent is set to 0 (the default), then the lesser of the absolute values of ")
printf("ubound and lbound will generally take precedence (assuming they are equal distances from the center value).\n")
printf("\n")
return(null)
}

cen = 0.0
ub = 1.0
lb = -1.0
ind = 0
anorm = 0

if(HasValue(center) == 1) {
cen = center
}

if(HasValue(ubound) == 1) {
ub = ubound
}

if(HasValue(lbound) == 1) {
lb = lbound
}

if(HasValue(independent) == 1) {
ind = independent
}

if(HasValue(anormalize) == 1) {
anorm = anormalize
}

if (anorm == 1) {
ind = 1
}

stepCount = 100
bl = 1
thresh = 0.05
trepeat = 3
swidth = 5

if(HasValue(steps) == 1) {
stepCount = steps
}

if(HasValue(backlash) == 1) {
bl = backlash
}

if(HasValue(threshold) == 1) {
thresh = threshold
}

if(HasValue(threshrepeat) == 1) {
trepeat = threshrepeat
}

if(HasValue(searchwidth) == 1) {
swidth = searchwidth
}

printf("Normalizing the input value about %f, with an upper bound of %f and a lower bound of %f.\n", cen, ub, lb)

if (ind == 1) {
printf("Positive and negative values will be normalized independently.\n")
} else {
printf("Positive and negative values will be normalized using the same scaling factor.\n")
}

inVal = \$1 - cen

posVals = clipdata(inVal, min = 0.0, max = max(inVal))
negVals = abs(clipdata(inVal, min = min(inVal), max = 0.0))

if (anorm == 1) {
printf("Aggressive normalization is enabled.\n")
posVals = anormalize(posVals, steps = stepCount, backlash = bl, threshold = thresh, threshrepeat = trepeat, searchwidth = swidth)
negVals = anormalize(negVals, steps = stepCount, backlash = bl, threshold = thresh, threshrepeat = trepeat, searchwidth = swidth)
}

maxPos = max(posVals)
maxNeg = max(negVals)

lb = abs(lb)

posSF = ub / maxPos
negSF = lb / maxNeg

if (independent == 0) {
sf = posSF
if (posSF > negSF) {
sf = negSF
}
posSF = sf
negSF = sf
}

posVals = (posVals * posSF)
negVals = (negVals * negSF)

outVal = posVals + (negVals * (-1.0))
outVal = outVal + cen
return(outVal)
}

# this is really slow and should only be used for debugging
goodsyntax = 1
# minimum of four arguments (position 1, colour 1, position 2, colour 2)
if (\$ARGC < 1) {
goodsyntax = 0
}
# number of arguments must be even
if ((\$ARGC > 2) == 1) {
goodsyntax = 0
}
if (goodsyntax == 0) {
printf("\nReplaces imaginary, infinite, undefined, and other problematic values in an array.\n")
printf("\n")
printf("r        = the value to replace with (default = 0)\n")
printf("\$1        = the input data\n")
printf("\n")
return(null)
}

input = \$1

r = 0

if(HasValue(replaceValue) == 1) {
r = replaceValue
}

result = input

arrayDim = dim(input)

xDim = arrayDim[1,1,1]
yDim = arrayDim[2,1,1]
zDim = arrayDim[3,1,1]

replacedCount = 0

zCurrent = 1
while (zCurrent <= zDim) {
yCurrent = 1
while (yCurrent <= yDim) {
xCurrent = 1
while (xCurrent <= xDim) {
currentVal = input[xCurrent, yCurrent, zCurrent]
if (currentVal != currentVal) {
replacedCount++
currentVal = 0
}
result[xCurrent, yCurrent, zCurrent] = currentVal
xCurrent++
}
yCurrent++
}
zCurrent++
}

if (replacedCount > 0) {
printf("Replaced %i values.\n", replacedCount)
}

return(result)
}

# for use in emulating applications and equations that don't take negative numbers into account properly
# the "sc" is for "special case"
define powsc(){
if (\$ARGC != 2) {
printf("\nRaises values to the specified exponent, but scales the input values to the range 0.0 - 1.0 during the conversion.\n")
printf("(Values are expanded to their previous range after being raised to the specified exponent).\n")
printf("\nUsage is the same as the standard pow() function.\n")
}

input = \$1
exp = \$2

inmax = max(input)
inmin = min(input)

if ((inmax == 1.0) && (inmin == 0.0)) {
return(pow(input, exp))
}

orng = inmax - inmin

result = normalize(input)
result = pow(result, exp)
result = (result * orng) + inmin

return(result)
}

define findRedundantImageData(){
img1 = \$1
img2 = \$2

startMultiplier = 0.01
endMultiplier = 0.90
step = 0.01

m = startMultiplier

while (m <= endMultiplier) {
img1sep = clipdata((img1 - (img2 * m)), min = 0, max = 1, quiet=1) / (1.0 - m)
img2sep = clipdata((img2 - (img1 * m)), min = 0, max = 1, quiet=1) / (1.0 - m)

#img1overlap = (img1 - img1sep) / m
#img2overlap = (img2 - img2sep) / m

#img1overlap = normalize(clipdata((img1 - img2sep), min = 0, max = 1, quiet=1) / m)
#img2overlap = normalize(clipdata((img2 - img1sep), min = 0, max = 1, quiet=1) / m)

#img1overlap = clipdata((img1 - img2sep), min = 0, max = 1, quiet=1) / m
#img2overlap = clipdata((img2 - img1sep), min = 0, max = 1, quiet=1) / m

img1overlap = clipdata((img1 - img2sep), min = 0, max = 1, quiet=1) / m
img2overlap = clipdata((img2 - img1sep), min = 0, max = 1, quiet=1) / m

diff = img1overlap - img2overlap
absdiff = abs(diff)

diffAvg = avg(diff)
diffMed = med(diff)
absdiffAvg = avg(absdiff)
absdiffMed = med(absdiff)

printf("%f\t%f\t%f\t%f\t%f\n", m, diffAvg, diffMed, absdiffAvg, absdiffMed)
#printf("%f\t%f\t%f\n", m, diffAvg, absdiffAvg)

m += step
}

}

define divideNZ(default){
if (\$ARGC != 2) {
printf("\nPerforms a division, but returns zero (or another value) where a divide-by-zero error would occur.\n")
printf("\$1\t=\tnumerator\n")
printf("\$2\t=\tdenominator\n")
printf("default\t=\tdefault value\n")
printf("\ne.g. result = divideNZ(a, b).\n")
}

n = \$1
d = \$2

df = 0.0

if(HasValue(default) == 1) {
df = default
}

dmask = clipdata(abs(ceil(d)), min=0.0, max=1.0, quiet=1)
one[0, 0, 0] = 1.0

n1 = n * dmask
d1 = d * dmask

result = (n1 / (d1 + dmaskneg)) + (dmaskneg * df)
return(result)
}

define getMaxValForPS(inC1v, inC2v, inC3v){
result = inC1v

c2mask = clipdata(abs(ceil(inC2v)), min=0.0, max=1.0, quiet=1)
c2 = inC2v * c2mask

c3mask = clipdata(abs(ceil(inC3v)), min=0.0, max=1.0, quiet=1)
c3 = inC3v * c3mask

ratioC1ToC2 = divideNZ(inC1v, c2, default=1.0)
ratioC1ToC3 = divideNZ(inC1v, c3, default=1.0)

c2c3mask = clipdata(abs(ceil(inC2v - inC3v)), min=0.0, max=1.0, quiet=1)
one[0, 0, 0] = 1.0

ratioC1ToC2 = ratioC1ToC2 * c2c3mask
ratioC1ToC3 = ratioC1ToC3 * c3c2mask

r = cat(ratioC1ToC2, ratioC1ToC3, axis = z)
r = max(r, axis = z)

rpos = ceil(clipdata((r - 0.5), min = 0.0, quiet=1))

result = clipdata(inC1v + rpos, max = 1.0, quiet=1)

return(result)
}

define colourIsolationInternal(omode, inC1, inC2, inC3, c1c2Isolation, c1c3Isolation){
o1 = inC1
o2 = inC1
maxval = inC1
maxval[0,0,0] = 1.0
one = maxval

if (c1c2Isolation != 0.0){
if (omode == "ps"){
maxval = getMaxValForPS(inC1v=inC1, inC2v=inC2, inC3v=inC3)
}

if ((omode == "p") || (omode == "s")){
maxval = max(cat(inC1, inC2, axis = z), axis = z)
}
o1 = (inC1 - (inC2 * c1c2Isolation)) / (one - c1c2Isolation)
o1 = min(cat(o1, maxval, axis = z), axis = z)
}
if (c1c3Isolation != 0.0){
if (omode == "ps"){
maxval = getMaxValForPS(inC1v=inC1, inC3v=inC3, inC2v=inC2)
}

if ((omode == "p") || (omode == "s")){
maxval = max(cat(inC1, inC3, axis = z), axis = z)
}
o2 = (inC1 - (inC3 * c1c3Isolation)) / (one - c1c3Isolation)
o2 = min(cat(o2, maxval, axis = z), axis = z)
}

o1mask = clipdata(ceil(o1 - o2), min=0.0, max=1.0, quiet=1)
#o2mask = clipdata(ceil(o2 - o1), min=0.0, max=1.0, quiet=1)

if (omode == "ps"){
maxval = getMaxValForPS(inC1v=inC1, inC3v=inC3, inC2v=inC2)
}
if ((omode == "p") || (omode == "s")){
maxval = max(cat(inC1, inC3, axis = z), axis = z)
}
result1 = (o1 - (inC3 * c1c3Isolation)) / (one - c1c3Isolation)
result1 = min(cat(result1, maxval, axis = z), axis = z)

if (omode == "ps"){
maxval = getMaxValForPS(inC1v=inC1, inC2v=inC2, inC3v=inC3)
}
if ((omode == "p") || (omode == "s")){
maxval = max(cat(inC1, inC2, axis = z), axis = z)
}
result2 = (o2 - (inC2 * c1c2Isolation)) / (one - c1c2Isolation)
result2 = min(cat(result2, maxval, axis = z), axis = z)

result1 = result1 * o1mask
result2 = result2 * o2mask

result = result1 + result2

return(result)
}

define colourIsolation(mode, master, rg, gb, rb){
goodsyntax = 1

if (\$ARGC != 1) {
goodsyntax = 0
}

om = "u"

if(HasValue(mode) == 1) {
om = mode
}

if ((om != "p") && (om != "s") && (om != "ps") && (om != "u")) {
goodsyntax = 0
}

if (goodsyntax == 0) {
printf("\nMagnifies the difference between colours.\n\n")
printf("\$1\t= the input image\n")
printf("mode\t= operational mode, which is one of:\n")
printf("\t  \"p\" - isolate primary colours\n")
printf("\t  \"s\" - isolate secondary colours\n")
printf("\t  \"ps\" - isolate primary and secondary colours\n")
printf("\t  \"u\" - unrestricted\n")
printf("\t  (default: \"u\"\n")
printf("master\t= overall effect multiplier (> 0.0, < 1.0, default = 0.3)\n")
printf("rg\t= red/green isolation (0.0 - 1.0, default = 1.0)\n")
printf("gb\t= green/blue isolation (0.0 - 1.0, default = 0.7)\n")
printf("\n")
return(null)
}

inimg = \$1

mi = 0.3
rgi = 1.0
gbi = 0.7
rbi = 0.1

if(HasValue(master) == 1) {
mi = master
}

if(HasValue(rg) == 1) {
rgi = rg
}

if(HasValue(gb) == 1) {
gbi = gb
}

if(HasValue(rb) == 1) {
rbi = rb
}

rgi = rgi * mi
gbi = gbi * mi
rbi = rbi * mi

printf("Operating Mode: ")
if (om == "p") {
printf("Isolate Primary Colours\n")
}
if (om == "s") {
printf("Isolate Secondary Colours\n")
}
if (om == "ps") {
printf("Isolate Primary And Secondary Colours\n")
}
if (om == "u") {
printf("Unrestricted\n")
}

printf("Red/Green Isolation: %f\n", rgi)
printf("Green/Blue Isolation: %f\n", gbi)
printf("Red/Blue Isolation: %f\n", rbi)

inR = inimg[0, 0, 1]
inG = inimg[0, 0, 2]
inB = inimg[0, 0, 3]

outR = inR
outG = inG
outB = inB

if (om == "s") {
one = inR
one[0, 0, 0] = 1.0
inC = one - inR
inM = one - inG
inY = one - inB

outC = colourIsolationInternal(omode="s", inC1=inC, inC2=inM, inC3=inY, c1c2Isolation=rgi, c1c3Isolation=rbi)
outM = colourIsolationInternal(omode="s", inC1=inM, inC2=inC, inC3=inY, c1c2Isolation=rgi, c1c3Isolation=gbi)
outY = colourIsolationInternal(omode="s", inC1=inY, inC2=inM, inC3=inC, c1c2Isolation=gbi, c1c3Isolation=rbi)

outR = one - outC
outG = one - outM
outB = one - outY
# This next set of else ifs is a workaround for DaVinci's apparent inability to pass string variables to other functions
} else if (om == "p") {
outR = colourIsolationInternal(omode="p", inC1=inR, inC2=inG, inC3=inB, c1c2Isolation=rgi, c1c3Isolation=rbi)
outG = colourIsolationInternal(omode="p", inC1=inG, inC2=inR, inC3=inB, c1c2Isolation=rgi, c1c3Isolation=gbi)
outB = colourIsolationInternal(omode="p", inC1=inB, inC2=inG, inC3=inR, c1c2Isolation=gbi, c1c3Isolation=rbi)
} else if (om == "ps") {
outR = colourIsolationInternal(omode="ps", inC1=inR, inC2=inG, inC3=inB, c1c2Isolation=rgi, c1c3Isolation=rbi)
outG = colourIsolationInternal(omode="ps", inC1=inG, inC2=inR, inC3=inB, c1c2Isolation=rgi, c1c3Isolation=gbi)
outB = colourIsolationInternal(omode="ps", inC1=inB, inC2=inG, inC3=inR, c1c2Isolation=gbi, c1c3Isolation=rbi)
} else {
outR = colourIsolationInternal(omode="u", inC1=inR, inC2=inG, inC3=inB, c1c2Isolation=rgi, c1c3Isolation=rbi)
outG = colourIsolationInternal(omode="u", inC1=inG, inC2=inR, inC3=inB, c1c2Isolation=rgi, c1c3Isolation=gbi)
outB = colourIsolationInternal(omode="u", inC1=inB, inC2=inG, inC3=inR, c1c2Isolation=gbi, c1c3Isolation=rbi)
}

result = cat(cat(outR, outG, axis = z), outB, axis = z)

return(result)

}