// $BACKGROUND$

number redon, greenon, blueon
number RL, RH, RC, GL, GH, GC, BL, BH, BC, SingleColorNormalization  // wavelengths for color bands in stack


/////////////////////				USER SET PARAMETERS: 				///////////////////////

number Bandwidth = 10		// <== BANDWIDTH

// if you would like unequal bandwidths, simply set low and high wavelengths manually

// RED Band:
RC = 560				// <== RED band CENTER WAVELENGTH
RL = RC - Bandwidth/2	// <== RED band Start wavelength
RH = RC + Bandwidth/2	// <== RED band End wavelength

// GREEN Band:
GC = 450				// <== GREEN band CENTER WAVELENGTH
GL = GC - Bandwidth/2	// <== GREEN band Start wavelength
GH = GC + Bandwidth/2	// <== GREEN band End wavelength

// BLUE Band:
BC = 370				// <== BLUE band CENTER WAVELENGTH
BL = BC - Bandwidth/2	// <== BLUE band Start wavelength
BH = BC + Bandwidth/2	// <== BLUE band End wavelength



// QUANTITATIVE RESULT: ////////////////////////////////////////////////////////////////////////
// "single color normalization" (quantitative):
// ON (TRUE) = 1; OFF (FALSE) =0

SingleColorNormalization = 0	// <== 0: Normalize each band to itself; 1: Normalize each band to band of maximum integrated intensity

redon = 1 		// Red band ON =1  OFF = 0
greenon = 1		// Green band ON =1  OFF = 0
blueon = 1		// Blue band ON =1  OFF = 0
// SCRIPT ///////////////////////////////////////////////////////////////////////////////////
// var init
image dataset := GetFrontImage()
number wavelength 
Number Nx = ImageGetDimensionSize(dataset, 0)
Number Ny = ImageGetDimensionSize(dataset, 1)
Number Nz = ImageGetDimensionSize(dataset, 2)

image NewData := RGBimage("Colorband Image",4, Nx, Ny, 1)
image ColorData := RGBimage("Colorband Image Stack",4, Nx, Ny, 3)
image rawimage := Realimage("redframe",4, Nx, Ny)
image redimage := Realimage("redframe",4, Nx, Ny)
image greenimage := Realimage("greenframe",4, Nx, Ny)
image blueimage := Realimage("blueframe",4, Nx, Ny)
image newimage := RGBimage("singleframe",4, Nx, Ny, 1)
image newblueimage := RGBimage("singleblueframe",4, Nx, Ny, 1)
image newgreenimage := RGBimage("singlegreenframe",4, Nx, Ny, 1)
image newredimage := RGBimage("singleredframe",4, Nx, Ny, 1)
number Scale = dataset.ImageGetDimensionScale(2)      // spectrum scale
number Origin = dataset.ImageGetDimensionOrigin(2)    // spectrum origin
string Unit = dataset.ImageGetDimensionUnitString(2)  // spectrum units
number Scalex = dataset.ImageGetDimensionScale(1)     // spatial scale
number Originx = dataset.ImageGetDimensionOrigin(1)   // spatial origin
string Unitx = dataset.ImageGetDimensionUnitString(1) // spatial unit

number sliceTowavelength(number plane, image data)
{
	wavelength = Scale*plane + Origin
	return wavelength
}
number j
for (j=0; j<Nz; j++)
{
	wavelength = sliceToWavelength(j,dataset)

	If (wavelength > BL && wavelength < BH)
	{// compile blue image
	rawimage = Slice2(dataset, 0, 0, j, 0, Nx, 1, 1, Ny, 1)
	blueimage = blueon*(blueimage + rawimage)
	}
	
	If (wavelength > GL && wavelength < GH)
	{// compile green image
	rawimage = Slice2(dataset, 0, 0, j, 0, Nx, 1, 1, Ny, 1)
	greenimage = greenon*(greenimage + rawimage)
	}
	
	If (wavelength > RL && wavelength < RH) 
	{// compile red image
	rawimage = Slice2(dataset, 0, 0, j, 0, Nx, 1, 1, Ny, 1)
	redimage = redon*(redimage + rawimage)
	}
}

if (SingleColorNormalization == 1)
{ //procedure for normalization of each band to the band with max intensity
number maxx
if (max(blueimage) > max(greenimage) && max(blueimage) > max(redimage))
{
maxx = max(blueimage)
}
if (max(greenimage) > max(blueimage) && max(greenimage) > max(redimage))
{
maxx = max(greenimage)
}
if (max(redimage) > max(greenimage) && max(redimage) > max(blueimage))
{
maxx = max(redimage)
}

blueimage = 255*blueimage/maxx 
newblueimage = RGB(0*blueimage,0*blueimage,1*blueimage)
greenimage = 255*greenimage/maxx 
newgreenimage = RGB(0*greenimage,1*greenimage,0*greenimage)
redimage = 255*redimage/maxx 
newredimage = RGB(1*redimage,0*redimage,0*redimage)
}

if (SingleColorNormalization == 0)
{ //procedure for normalization of each band to itself
blueimage = 255*blueimage/max(blueimage)
newblueimage = RGB(0*blueimage,0*blueimage,1*blueimage)
greenimage = 255*greenimage/max(greenimage)
newgreenimage = RGB(0*greenimage,1*greenimage,0*greenimage)
redimage = 255*redimage/max(redimage)
newredimage = RGB(1*redimage,0*redimage,0*redimage)
}

NewData = RGB(redimage,greenimage,blueimage) // generate composite RGB image

{// input original spatial calibrations
ImageSetDimensionScale(NewData,0,Scalex)
ImageSetDimensionOrigin(NewData,0,Originx)
ImageSetDimensionUnitString(NewData,0,Unitx)
ImageSetDimensionScale(NewData,1,Scalex)
ImageSetDimensionOrigin(NewData,1,Originx)
ImageSetDimensionUnitString(NewData,1,Unitx)
}
NewData.showImage() // show composite RGB image

{// create 3D image stack [R, G, B]
Slice2(ColorData, 0, 0, 0, 0, Nx, 1, 1, Ny, 1) = newredimage
Slice2(ColorData, 0, 0, 1, 0, Nx, 1, 1, Ny, 1) = newgreenimage
Slice2(ColorData, 0, 0, 2, 0, Nx, 1, 1, Ny, 1) = newblueimage
}

{// input original spatial calibrations
ImageSetDimensionScale(ColorData,0,Scalex)
ImageSetDimensionOrigin(ColorData,0,Originx)
ImageSetDimensionUnitString(ColorData,0,Unitx)

ImageSetDimensionScale(ColorData,1,Scalex)
ImageSetDimensionOrigin(ColorData,1,Originx)
ImageSetDimensionUnitString(ColorData,1,Unitx)
}

ColorData.showImage() // show 3D image stack