

/*********************					*
**************************				**
COMMON RENDERING ***************************
**************************				**
*********************					*/







/******************************************************************
function: Project 2D to left
******************************************************************/
image  MSA112_Render2DtoLeft(image data, number outputFormat)
{
number w,h
data.getSize(w,h)

image uneven,split,half 
if(outputFormat==0) uneven :=  MSA112_Make2DImageOfRightType(data.ImageGetDataType(),"",1,h)
else 			 	uneven :=  MSA112_Make2DImageOfRightType(2,"",1,h)	//real

number s = w,w1,w2

//initial split
if(outputFormat==0) split := data.imageclone()
else			 split = data	//real

while(s>1)
	{
	w1 = trunc(s/2)
//result("->"+w1)
	w2=2*w1
	half := split[0,0,h,w1]
	half += split[0,w1,h,w2]
	if(w2<s) uneven += split[0,w2,h,s]
	split := half
	s=w1
	}
				
split += uneven
		
return split
}	


/******************************************************************
function: Project 2D to bottom
******************************************************************/
image  MSA112_Render2DtoBottom(image data, number outputFormat)
{
number w,h
data.getSize(w,h)

image uneven,split,half
if(outputFormat==0) uneven :=  MSA112_Make2DImageOfRightType(data.ImageGetDataType(),"",w,1)
else				uneven :=  MSA112_Make2DImageOfRightType(2,"",w,1)	//real

number s = h,h1,h2

//initial split
if(outputFormat==0) split := data.imageclone()
else			 	split = data	//real


while(s>1)
	{
	h1 = trunc(s/2)
//result("->"+w1)
	h2=2*h1
	half := split[0,0,h1,w]
	half += split[h1,0,h2,w]
	if(h2<s) uneven += split[h2,0,s,w]
	split := half
	s=h1
	}
				
split += uneven
		
return split
}
/*
number height = 400
height = height**2
image in:= realimage("",4,2000,height)
in= icol
number t = getOstickCount()
image spec 
spec = in.MSA112_Render2DtoBottom()
	//spec = in.rebin(1,height)
result("\n"+CalcOSsecondsBetween(t,getOstickCount()))
*/	


/******************************************************************
function: Project 3D to spectrum
******************************************************************/
image  MSA112_Render3DtoPixel(image data)
{
return data.MSA112_Render3DtoLeft(1).MSA112_Render2DtoBottom(1)
}

/*
number width = 300
number depth = 2000  
image cube := realimage("",4,width,width,depth)
cube = iplane*icol //+icol

number t = getOstickCount()
image spec = cube.MSA112_Render3DtoPixel()
result("\n"+CalcOSsecondsBetween(t,getOstickCount()))
//spec.showimage()
*/


/******************************************************************
function: compact form of rendering commands , no output format, always real
******************************************************************/
image  MSA112_Render2DtoLeft(image data) return MSA112_Render2DtoLeft(data,1)
image  MSA112_Render2DtoBottom(image data) return MSA112_Render2DtoBottom(data,1)
image  MSA112_Render3DtoLeft(image data) return MSA112_Render3DtoLeft(data,1)
image  MSA112_Render3DtoFront(image data) return MSA112_Render3DtoFront(data,1)
image  MSA112_Render3DtoLeft3(image data) return MSA112_Render3DtoLeft3(data,1)
image  MSA112_Render3DtoBottom3(image data) return MSA112_Render3DtoBottom3(data,1)













/*********************					*
**************************				**
MATRIX OPERATIONS***************************
**************************				**
*********************					*/

	
	

/******************************************************************
function: Convert column image into rectangular image
******************************************************************/
image MSA112_ConvertColumnTo2D(image in,number width,number height)
{
return Slice2(in,0,0,0,		1,width,1,		1,height,width)	

}


/******************************************************************
function: Convert rectangular image into  column image
******************************************************************/
image MSA112_Convert2DtoColumn(image in)
{
number width,height
in.getSize(width,height)

object stream = NewStreamFromBuffer(0)
ImageWriteImageDataToStream(in,stream, 0)

//image out :=realImage("",4,1,width*height)
image out :=  MSA112_Make2DImageOfRightType(in.ImageGetDataType(),"",1,width*height)

stream.StreamSetPos(0,0)
ImageReadImageDataFromStream(out,stream, 0) 

return out
}

/******************************************************************
function: Convert 3D image into 2D image
******************************************************************/
image MSA112_Convert3Dto2D(image in)
{
number width,height,depth
in.get3Dsize(width,height,depth)

image out :=  MSA112_Make2DImageOfRightType(in.ImageGetDataType(),"",width*height,depth)

object stream = NewStreamFromBuffer(0)
ImageWriteImageDataToStream(in,stream, 0)
	
stream.StreamSetPos(0,0)
ImageReadImageDataFromStream(out,stream, 0) 
	
return out.MatrixTranspose()	
}

/******************************************************************
function: Convert dataCube into 2D matrix
******************************************************************/
image MSA112_ConvertCubeToMatrix(image DataCube,number Estart,number Eend)
{
//output data matrix is always REAL 4-Byte !!!!!
number width = DataCube.ImageGetDimensionSize(0)
number height = DataCube.ImageGetDimensionSize(1)
number depth = Eend-Estart

string DataName = DataCube.getName()

image dataMatrix :=realImage(DataName+" MSA",4,depth,width*height)
//always real matrix !!

dataMatrix = MSA112_Convert3Dto2D(DataCube[0,0,Estart,width,height,Eend])

return dataMatrix
}
	
/****************** simplified version without energy range **************/
image MSA112_ConvertCubeToMatrix(image DataCube)
{
number Depth =  DataCube.ImageGetDimensionSize(2)

return MSA112_ConvertCubeToMatrix(DataCube,0,Depth)
}


/******************************************************************
function: Convert 2D matrix into data cube
******************************************************************/
image MSA112_Convert2Dto3D(image in,number width,number height)
{
number Depth = in.ImageGetDimensionSize(0)

return Slice3(in,0,0,0,  1,width,1,  1,height,width, 0,depth,1)
}




/******************************************************************
Function Normalize vector
******************************************************************/	
/************************  calculate norm *********************/
number MSA112_Norm(image Vec) 
{
number width,height,norm
Vec.getSize(width,height)

if(width==1) norm = MatrixMultiply(MatrixTranspose(Vec),Vec).sum()
if(height==1) norm = MatrixMultiply(Vec,MatrixTranspose(Vec)).sum()
if((width!=1)&&(width!=1)) {result("\nVECTOR IS NEITHER ROW NOR COLUMN !!");return 1;}

norm = sqrt(norm)

return norm
}


/************** return also the norm***********************/	
image MSA112_NormalizeVector(image vector,number &Norm)
	{
	Norm = MSA112_Norm(vector)
	vector /= Norm
	
	return Vector
	}

/************** simple return ***********************/	
image MSA112_NormalizeVector(image vector)
	{
	number norm
	vector = MSA112_NormalizeVector(vector,norm)
	
	return vector
	}


/******************************************************************
Function Normalize matrix
******************************************************************/	
/************************  calculate norm *********************/
image MSA112_ColumnNorm(image Matrix) 
{
number width,height
Matrix.getSize(width,height)

image Squares = imageclone(Matrix)**2

image Norm = MSA112_Render2DtoBottom(Squares)
Norm = sqrt(Norm)

return Norm
}

/************** return also the columnNorm ***********************/	
image MSA112_NormalizeMatrix(image Matrix,image &ColumnNorm)
	{		
	ColumnNorm = MSA112_ColumnNorm(Matrix)
	Matrix /= MSA112_ExpandImageVertically(ColumnNorm,Matrix)
	
	return Matrix
	}

/************** simple return ***********************/	
image MSA112_NormalizeMatrix(image Matrix)
	{
	image ColumnNorm
	
	return MSA112_NormalizeMatrix(Matrix,ColumnNorm)
	}

/************** with L1 (not L2 norm) norm ***********************/
image MSA112_NormalizeMatrixL1(image matrix,image &matrixNorm,image &columnNorm)
	{
	columnNorm = MSA112_Render2DtoBottom(matrix)
	matrixNorm = MSA112_ExpandImageVertically(columnNorm,matrix)

	image normalizedMatrix =matrix/matrixNorm

	return normalizedMatrix
	}

/************** simple return ***********************/
image MSA112_NormalizeMatrixL1(image matrix)
	{
	image MatrixNorm,columnNorm
	
	return matrix.MSA112_NormalizeMatrixL1(MatrixNorm,columnNorm)
	}
	
	
/******************************************************************
function: Generate Random vector
******************************************************************/
image  MSA112_RandomVector(number dim)
{
image vec := realimage("",4,1,dim)

//for(number i=0;i<dim;i++) vec.setPixel(0,i,random()-0.5)
vec = random() -0.5

return MSA112_NormalizeVector(vec)
}

/******************************************************************
function: Generate Random matrix
******************************************************************/
image  MSA112_RandomMatrix(number dimX,number dimY)
{
image matrix := realimage("",4,1,dimX,dimY)

matrix = random() -0.5

return MSA112_NormalizeMatrix(matrix)
}


	
/******************************************************************
Function Find angular difference between vectors
******************************************************************/	
number MSA112_AngularDifference(image v1,image v2)
	{
	number norm1 = MatrixMultiply(MatrixTranspose(v1),v1).sum()
	image v1Normed = v1/sqrt(norm1)
	
	number norm2 =MatrixMultiply(MatrixTranspose(v2),v2).sum()
	image v2Normed = v2/sqrt(norm2)

	number proj = MatrixMultiply(MatrixTranspose(v1Normed),v2Normed).sum()
	number dev = acos(proj)*1000
//result("\n"+norm1+" "+norm2+" " +dev+" mrad")

	return dev
	}


/**************** 1D Gaussian smoothing **********************/
image MSA112_Smooth(image in,number X,number Y,number sigma)
	{
	//X=0,1 disable/enable X, Y=1 disable/enable Y
	// gaussian sigma in pixels,number of averaged images 1+2*N	
	//sigma 1: N = 1, images 3
	//sigma 2: N = 3, images 7
	//sigma 3: N = 4, images 9
	//sigma 4: N = 6, images 13
	//sigma 5: N = 7, images 15
	//sigma 6: N = 9, images 19
	//sigma 7: N = 10, images 21
		
	number minTail = 0.1
	number N = (sqrt(-log(minTail))*sigma).trunc()
	image out = in.imageclone()
	number Norm = 1
	for(number i=1;i<=N;i++)
		{
		number weight = exp(-(i/sigma)**2)

		out += in[icol-i*X,irow-i*Y]*weight
		out += in[icol+i*X,irow+i*Y]*weight
		Norm += 2*weight
		}
		
	return out/Norm
	}



image MSA112_der23(image in,number X,number Y)
	{
	//X=0,1 disable/enable X, Y=1 disable/enable Y
	//points three/five points used in estimation

	image Near = imageclone(in)
	Near = in[icol-X,irow-Y]
	Near += in[icol+X,irow+Y]

	return (Near-2*in)
	}

image MSA112_der25(image in,number X,number Y)
	{
	//X=0,1 disable/enable X, Y=1 disable/enable Y
	//points three/five points used in estimation

	image Near = imageclone(in)
	Near = in[icol-X,irow-Y]
	Near += in[icol+X,irow+Y]

	image Next = imageclone(in)
	Next = in[icol-2*X,irow-2*Y]
	Next += in[icol+2*X,irow+2*Y]
	
	return (16*Near-Next-30*in)/12
	}







/*********************					*
**************************				**
EXTRACT AVERAGE ***************************
**************************				**
*********************					*/

	
	
	
	
	
	
	
	
	

	
	
/******************************************************************
function: Extract meanImage from data cube (calibrated)
******************************************************************/
image MSA112_MakeMeanImage (image DataCube)
	{
	string Name = DataCube.ImageGetName()

	number depth = DataCube.ImageGetDimensionSize(2)
	
	image meanImage = DataCube.MSA112_Render3DtoFront()/depth	
	meanImage.ImageSetName(Name+" average image")

	number originX,originY,scaleX,scaleY
	string unitsX,unitsY
	
	ImageGetDimensionCalibration(DataCube,0,originX,scaleX,unitsX,1)
	ImageGetDimensionCalibration(DataCube,1,originY,scaleY,unitsY,1)

	ImageSetDimensionCalibration(meanImage,0,originX,scaleX,unitsX,1)
	ImageSetDimensionCalibration(meanImage,1,originY,scaleY,unitsY,1)
	
	return meanImage
	}


/******************************************************************
function: Extract averageSpectrum from data cube (calibrated)
******************************************************************/
image MSA112_MakeAverageSpectrum (image DataCube)
	{
	string Name = DataCube.ImageGetName()
	
	number width = DataCube.ImageGetDimensionSize(0)
	number height = DataCube.ImageGetDimensionSize(1)

	image avSpectrum = DataCube.MSA112_Render3DtoPixel()
	AvSpectrum /= (width*height)
	avSpectrum.ImageSetName(Name+" average spec")
	
	number originE,disp
	string units
	
	ImageGetDimensionCalibration(DataCube,2,originE,disp,units,1)

	ImageSetDimensionCalibration(avSpectrum,0,originE,disp,units,1)
	
	return avSpectrum
	}



/******************************************************************
function: Extract meanImage from matrix (column, UNcalibrated)
******************************************************************/
image MSA112_TakeMeanImage(image dataMatrix)
{
number Depth = dataMatrix.ImageGetDimensionSize(0)
image AvImage = MSA112_Render2DtoLeft(dataMatrix)
AvImage /= Depth

return AvImage
}


/******************************************************************
function: Extract averageSpectrum from matrix (row, UNcalibrated)
******************************************************************/
image MSA112_TakeAverageSpectrum(image dataMatrix)
{
number kN = dataMatrix.ImageGetDimensionSize(1)
image AvSpec = MSA112_Render2DtoBottom(dataMatrix)
AvSpec /= kN

return AvSpec
}


/******************************************************************
function: Reduce matrix by averageSpectrum
******************************************************************/
image MSA112_ReduceAverageSpectrum(image dataMatrix)
{
image AvSpec = MSA112_TakeAverageSpectrum(dataMatrix)

dataMatrix -= MSA112_ExpandImageVertically(AvSpec,dataMatrix)

return dataMatrix
}	







/*********************					*
**************************				**
CLUSTERS OPERATIONS ***********************
**************************				**
*********************					*/

	
	
	



	
/******************************************************************
Function Find ini and final clusters indexes
******************************************************************/
void MSA112_ClustersRange(number ClusterS,number ClusterN,number &Cluster0,number &ClusterF)
{
//if clusterN == 0 	return	Cluster0=0 , ClusterF=0
//if clusterN != 0	check if clusterSelected is valid and return Cluster0=ClusterS , ClusterF=ClusterS
//if ClusterS=(-1) makes package treatment 	Cluster0=1 , ClusterF=ClusterN
if(clusterN==0) 
	{
	if(ClusterS>0) if(!OkCancelDialog("there are no clusters in data.\nContinue without clustering ?")) throw("")
	Cluster0=0;//no matter what is clusterS
	ClusterF=0;
	}	
else 
	{
	if(ClusterS<(-1)) {OkDialog("Enter the valid cluster number !!");throw("");}
	if(ClusterS==(-1)) {Cluster0=1;ClusterF=ClusterN;} //package treatment
	if(ClusterS==0) {Cluster0=0;ClusterF=0;} //no clustering	
	if(ClusterS>=1) {Cluster0=ClusterS;ClusterF=ClusterS;}	//specific clusterS defined
	}
}


	
/******************************************************************
Function Take Cluster
******************************************************************/
/**************  reduces the dataMatrix, return kPoints - number of points in cluster **********************/
image MSA112_TakeCluster(image &dataMatrix,image cluster,number currCluster,number &kPoints)
{
number channels,Points
dataMatrix.getSize(channels,Points)

//extract kPoints - number points in cluater
image mask = imageclone(cluster)
if(currCluster>0) mask = tert(cluster==currCluster,1,0)
else 
	{
	kPoints = Points

	return dataMatrix
	}
kPoints = sum(mask)
if(kPoints==0) return NULL
//mask.matrixTranspose().showimage()

//prepare indexShifting image
//(increments as irow where mask==0, stays const where mask==1)
image Shifting = imageclone(mask)
Shifting = 0
//temp1 = tert(mask==0,offset(temp1,0,-1)+1,offset(temp1,0,-1))
Shifting = tert(mask==0,Shifting[0,irow-1]+1,Shifting[0,irow-1])
//Shifting.matrixTranspose().showimage()

//compress indexShiftImage
//(remove the areas of irow-increment)
image ShiftingComp = imageclone(mask)
ShiftingComp = 0
Index(ShiftingComp,0,irow-Shifting+1)= Shifting
//ShiftingComp.matrixTranspose().showimage()

//cut of dummy area
image maskRed := integerimage("reduced mask I",2,0,1,kPoints)
maskRed = ShiftingComp[0,0,kPoints,1]

//convert inot 2D real mask
image maskRed1D := realimage("reduced mask R",4,1,kPoints)	
maskRed1D = maskRed
image maskRed2D := realimage("reduced mask real",4,channels,kPoints)
maskRed2D = maskRed1D[0,irow]

//form reduced data
image dataMatrixRed := realimage("reduced mask",4,channels,kPoints)
dataMatrixRed = dataMatrix[icol,irow+maskRed2D]

return dataMatrixRed
}


/**************  additionally makes centering and return Average spectrum ******************/
image MSA112_TakeCluster(image dataMatrix,image cluster,number currCluster,number MakeCent,image &Center,number &kN)
{
image dataMatrixRed := MSA112_TakeCluster(dataMatrix,cluster,currCluster,kN)

try
	{
	if(dataMatrixRed.ImageIsValid()) Center = MSA112_TakeAverageSpectrum(dataMatrixRed)
	else //no points in cluster
		{
		number depth =dataMatrix.ImageGetDimensionSize(0)
		result(" size "+dataMatrix.ImageGetDimensionSize(0))
		center :=realimage("",4,depth,1) //zero center

		return NULL
		}
	
	if(MakeCent==1)
		{
		dataMatrixRed -= MSA112_ExpandImageVertically(Center,dataMatrixRed)
		}
	}
catch
	{
	OKdialog("TakeCluster: no enough memory !!")
	
	return NULL
	}
	
return dataMatrixRed
}


/***************** simplified command without returning Average Spectrum and the number of points kN  ******************************/
image MSA112_TakeCluster(image dataMatrix,number MakeCent,image cluster,number currCluster)
{
//dummy
number kPoints
image Center

return MSA112_TakeCluster(dataMatrix,cluster,currCluster,MakeCent,Center,kPoints)
}












/*********************					*
**************************				**
GET/SET info in Tags *******************
**************************				**
*********************					*/

	
	
	
	
	
	
	
	
	
	
			/****ALL NUMBERS UNsigned *****************/	
			/****  "S"	short (2-byte ******************/	
			/****  "I"	integer (4-byte ******************/
			/****  "L"	long (8-byte ******************/

void MSA112_setNumberInTags(TagGroup TargetTags,string name,string type,number Value)
	{
	if(type=="S") TargetTags.TagGroupSetTagAsShort(name,Value)
	if(type=="I") TargetTags.TagGroupSetTagAsUInt32(name,Value)
	if(type=="L") TargetTags.TagGroupSetTagAsLong(name,Value)
	if(type=="F") TargetTags.TagGroupSetTagAsFloat(name,Value) 	 	
	}

	
number MSA112_getNumberFromTags(TagGroup TargetTags,string name,string type)
	{
	number Value
	if(type=="S") TargetTags.TagGroupGetTagAsShort(name,Value)
	if(type=="I") TargetTags.TagGroupGetTagAsUInt32(name,Value)
	if(type=="L") TargetTags.TagGroupGetTagAsLong(name,Value)	
	if(type=="F") TargetTags.TagGroupGetTagAsFloat(name,Value) 
	
	return Value
	}

void MSA112_setStringInTags(TagGroup TargetTags,string name,string Value)
	{
	TargetTags.TagGroupSetTagAsString(name,Value)
	}

string MSA112_getStringFromTags(TagGroup TargetTags,string name)
	{
	string Value
	TargetTags.TagGroupGetTagAsString(name,Value)
	
	return Value
	}
	
TagGroup MSA112_GetOrCreateItemInTagList(TagGroup TargetTags,string path, number index)
	{
	TagGroup TargetList
	if(TargetTags.TagGroupGetTagAsTagGroup(path,TargetList))
		{	//taggroup exists
		if(TargetList.TagGroupIsList())
			{	//targetList is list
	//result("\ntaggroup is list")
			TagGroup TargetItem
			if(TargetList.TagGroupGetIndexedTagAsTagGroup(index,TargetItem))
				{
	//result("\nsuch item exists")
				return TargetItem
				}
			else	//create items up to the target index 
				{
				number ListSize = TagGroupCountTags(TargetList)
	//result("\nadd "+(ListSize-index+1)+" items")		
				for(number i=listSize;i<=index;i++) TargetList.TagGroupAddTagGroupAtEnd(NewTagGroup())
				TargetList.TagGroupGetIndexedTagAsTagGroup(index,TargetItem)
				
				return TargetItem
				}	
			}	
		}
		//if pass here then TargetList doesnot exist or is not a list
		
		//create TargetList
		TargetList =NewTagList()
		TargetTags.TagGroupsetTagAsTagGroup(path,TargetList)
		
		//create items up to the target index
		for(number i=0;i<=index;i++) TargetList.TagGroupAddTagGroupAtEnd(NewTagGroup())
		
		//return it
		TagGroup TargetItem
		TargetList.TagGroupGetIndexedTagAsTagGroup(index,TargetItem)
		
		return TargetItem
	}
		
void MSA112_setNumberInList(TagGroup TargetTags,string path,number index,string name,string type,number Value)
	{			//creates or extends list if nessesary
	TagGroup TargetItem = MSA112_GetOrCreateItemInTagList(TargetTags,path, index)
	if(type=="S") TargetItem.TagGroupSetTagAsShort(name,Value)
	if(type=="I") TargetItem.TagGroupSetTagAsUInt32(name,Value)
	if(type=="L") TargetItem.TagGroupSetTagAsLong(name,Value) 
	if(type=="F") TargetItem.TagGroupSetTagAsFloat(name,Value) 
	}

number MSA112_getNumberfromList(TagGroup TargetTags,string path,number index,string name,string type)
	{
/*	
	TagGroup TargetItem
	if(!TargetTags.TagGroupGetIndexedTagAsTagGroup(index,TargetItem))
		{
		result("\nWarning: no such index in target list")
		return 0
		}
TargetTags.TagGroupOpenBrowserWindow(0)
*/		
	number Value
	if(type=="S") TargetTags.TagGroupGetTagAsShort(path+":["+index+"]:"+name,Value)
	if(type=="I") TargetTags.TagGroupGetTagAsUInt32(path+":["+index+"]:"+name,Value)
	if(type=="L") TargetTags.TagGroupGetTagAsLong(path+":["+index+"]:"+name,Value) 
	if(type=="F") TargetTags.TagGroupGetTagAsFloat(path+":["+index+"]:"+name,Value) 
//result("\n"+path+":["+index+"]:"+name+"#")	

	return Value
	}

number MSA112_CheckNullImage(number width,number height,number type)
	{
	if(width==0) {result("\nWARNING !! cannot read image width from tags;");return 0;}
	if(height==0) {result("\nWARNING !! cannot read image height from tags;");return 0;}
	if(type==0) {result("\nWARNING !! cannot read image type from tags;");return 0;}
	
	return 1
	}
	
void MSA112_setImageinTags(TagGroup TargetTags,string name,image im)
	{	
	number width,height
	im.getSize(width,height)
	number type =ImageGetDataType(im)
	
	TargetTags.TagGroupSetTagAsLong(name+":width",width)
	TargetTags.TagGroupSetTagAsLong(name+":height",height)
	TargetTags.TagGroupSetTagAsLong(name+":type",type)
	
	TargetTags.TagGroupSetTagAsArray(name+":array",im)
	}

image MSA112_getImagefromTags(TagGroup TargetTags,string name)
	{
	number width,height,type
	TargetTags.TagGroupGetTagAsLong(name+":width",width)
	TargetTags.TagGroupGetTagAsLong(name+":height",height)
	TargetTags.TagGroupGetTagAsLong(name+":type",type)
	
	if(!MSA112_CheckNullImage(width,height,type)) return NULL
	
	image im :=MSA112_Make2DImageOfRightType(type,"",width,height)
	if(!TargetTags.TagGroupGetTagAsArray(name+":array",im)) return NULL
	
	return im
	}
	
void MSA112_setImageinList(TagGroup TargetTags,string path,number index,string name,image im)
	{	//creates or extends list if nessesary
	TagGroup TargetItem = MSA112_GetOrCreateItemInTagList(TargetTags,path, index)
	
	number width,height
	im.getSize(width,height)
	number type =ImageGetDataType(im)
	
	TargetItem.TagGroupSetTagAsLong(name+":width",width)
	TargetItem.TagGroupSetTagAsLong(name+":height",height)	
	TargetItem.TagGroupSetTagAsLong(name+":type",type)	
	TargetItem.TagGroupSetTagAsArray(name+":array",im)
	}

image MSA112_getImagefromList(TagGroup TargetTags,string path,number index,string name)
	{
	number width,height,type
	TargetTags.TagGroupGetTagAsLong(path+":["+index+"]:"+name+":width",width)
	TargetTags.TagGroupGetTagAsLong(path+":["+index+"]:"+name+":height",height)
	TargetTags.TagGroupGetTagAsLong(path+":["+index+"]:"+name+":type",type)

	if(!MSA112_CheckNullImage(width,height,type)) return NULL
	
	image im :=MSA112_Make2DImageOfRightType(type,"",width,height)
	TargetTags.TagGroupGetTagAsArray(path+":["+index+"]:"+name+":array",im)
	
	return im
	}
	















/*********************					*
**************************				**
GET/SET IN MSA IMAGE **********************
**************************				**
*********************					*/

	
	
	
	
	
	
	
	
	




/******************************************************************
function: Get matrix tags
******************************************************************/
TagGroup MSA112_GetMatrixTags(image MSAimage,number full)
	{ 	
	TagGroup MatrixTags = MSAimage.ImageGetTagGroup()
	string type =MatrixTags.MSA112_getStringFromTags(MSA112_Label()+":type")
//result("\n"+type+"#")
	if((type!="MSA matrix")&&(type!="MSA compressed")) 
		{
		OkDialog("It must be an MSA container in front !!")
		return NULL
		}
		
	if((full==1)&&(type=="MSA compressed"))
		{
		OkDialog("This operation requires full (not compressed) MSA matrix in front !!")
		return NULL
		}
	
	return MatrixTags
	}

	

	
/******************************************************************
function: Get/Set cluster mask 
******************************************************************/

image MSA112_GetCluster(TagGroup MatrixTags)
	{
	string MSALabel = MSA112_Label()
	
	number Points = MatrixTags.MSA112_getNumberFromTags(MSAlabel+":pixels","L")
	if(Points==0) {OKdialog("number of pixels not stored");throw("");}
	
	image Cluster:= integerImage("cluster",2,0,1,Points)
	MatrixTags.TagGroupGetTagAsArray(MSAlabel+":Cluster Mask",cluster)
	//if didnot exist return image with 0s
	
	return Cluster
	}


/********************  from MSAimage ********************************/
image MSA112_GetCluster(image MSAimage)
	{
	TagGroup MatrixTags =  MSA112_GetMatrixTags(MSAimage,0)
	
	return MSA112_GetCluster(MatrixTags)
	}
	

/********* Set cluster mask (used in clustering) ******************/
void MSA112_SetCluster(image MSAimage,image cluster)
{
string MSAlabel = MSA112_Label()
TagGroup MatrixTags = MSA112_GetMatrixTags(MSAimage,1)

//get previous number of clusters
//number clusterNold =MatrixTags.MSA112_getNumberFromTags(MSAlabel+":clusters number","S")	
TagGroup ClusterList = MatrixTags.TagGroupGetOrCreateTagList(MSAlabel+":clusters")
number clusterNold = TagGroupCountTags(ClusterList) -1
	
//as clusters are rearranged, previous info about their size and PCA must be erased
try
	{
	for(number clust =clusterNold;clust>=1 ; clust--)
		{ 	
		ClusterList.TagGroupDeleteTagWithIndex(clust)		
		}
	}
catch {break;}

//get new number of clusters	
number clusterN = max(cluster)

//set list
if(clusterNold ==(-1)) ClusterList.TagGroupCreateGroupTagAtEnd()	//virgin, create zero-cluster
for(number i = 1;i<=clusterN;i++) ClusterList.TagGroupCreateGroupTagAtEnd()

//set cluster number
MatrixTags.MSA112_setNumberInTags(MSAlabel+":clusters number","S",clusterN)

//set mask
MatrixTags.TagGroupSetTagAsArray(MSAlabel+":Cluster Mask",cluster)
}



/******************************************************************
functions: Dimensions and Calibration
******************************************************************/

/*************** ALL DIMENSIONS IN LONG (8byte) ***************/
void MSA112_SetNumberPoints(TagGroup MatrixTags,number Points)
	{
	string MSALabel = MSA112_Label()
	
	MatrixTags.MSA112_setNumberInTags(MSAlabel+":pixels","L",Points)
	}
	
number MSA112_GetNumberPoints(TagGroup MatrixTags)
	{
	string MSALabel = MSA112_Label()
	
	number Points =MatrixTags.MSA112_getNumberFromTags(MSAlabel+":pixels","L")
	
	return Points
	}
	
number MSA112_GetNumberPoints(image MSAimage)
	{
	TagGroup MatrixTags =  MSA112_GetMatrixTags(MSAimage,0)
	
	return MSA112_GetNumberPoints(MatrixTags)
	}


/******************* same in clusters ********************/
void MSA112_SetNumberPoints(TagGroup MatrixTags,number currCluster, number Points)
	{
	string 	clusterLabel = MSA112_GetClusterLabel(currCluster)
	
	MatrixTags.MSA112_setNumberInTags(ClusterLabel+":pixels","L",Points)	
	}

number MSA112_GetNumberPoints(TagGroup MatrixTags,number currCluster)
	{
	string 	clusterLabel = MSA112_GetClusterLabel(currCluster)
	
	number Points =MatrixTags.MSA112_getNumberFromTags(ClusterLabel+":pixels","L")
	
	//reserve path for older version (with "points")
	if(Points==0) 
		{
		Points =MatrixTags.MSA112_getNumberFromTags(ClusterLabel+":points","L")
		if(Points!=0) result("\nWARNING: pixel number stored in older version")
		}
		
	return Points
	}

	
	
void MSA112_GetSIWidthHeight(TagGroup MatrixTags,number &width,number &height)
	{
	string MSALabel = MSA112_Label()
	
	width  =MatrixTags.MSA112_getNumberFromTags(MSAlabel+":width","L")
	height =MatrixTags.MSA112_getNumberFromTags(MSAlabel+":height","L")
	}

void MSA112_GetSIWidthHeight(image MSAimage,number &width,number &height)
	{
	TagGroup MatrixTags =  MSA112_GetMatrixTags(MSAimage,0)
	
	MSA112_GetSIWidthHeight(MatrixTags,width,height)
	}

void MSA112_GetSISizes(tagGroup MatrixTags,number &Points,number &Depth,number &width,number &height)
	{
	string MSALabel = MSA112_Label()

	Points =MatrixTags.MSA112_getNumberFromTags(MSAlabel+":pixels","L")
	Depth  =MatrixTags.MSA112_getNumberFromTags(MSAlabel+":energy depth","L")
	width  =MatrixTags.MSA112_getNumberFromTags(MSAlabel+":width","L")
	height =MatrixTags.MSA112_getNumberFromTags(MSAlabel+":height","L")
	}

void MSA112_GetSISizes(image MSAimage,number &Points,number &Depth,number &width,number &height)
	{
	TagGroup MatrixTags =  MSA112_GetMatrixTags(MSAimage,0)
	
	MSA112_GetSISizes(MatrixTags,Points,Depth,width,height)
	}
		

void MSA112_SetCalibration(tagGroup MatrixTags,number scale,number origin, string units,string axis)
	{
	string MSALabel = MSA112_Label()

	if(Axis=="E")	//for compatibility with previous versions
		{
		MatrixTags.MSA112_SetStringInTags(MSAlabel+":energy units",units)
		MatrixTags.MSA112_SetNumberInTags(MSAlabel+":energy dispersion","F",scale)
		MatrixTags.MSA112_SetNumberInTags(MSAlabel+":energy origin","L",origin)		
		}
	else	
		{
		MatrixTags.MSA112_SetStringInTags(MSAlabel+":units "+axis,units)
		MatrixTags.MSA112_SetNumberInTags(MSAlabel+":scale "+axis,"F",scale)
		MatrixTags.MSA112_SetNumberInTags(MSAlabel+":origin "+axis,"L",origin)
		}
	}

void MSA112_SetCalibration(image MSAimage,number scale,number origin, string units,string axis)
	{
	TagGroup MatrixTags =  MSA112_GetMatrixTags(MSAimage,0)
	
	MSA112_SetCalibration(MatrixTags,scale,origin,units,axis)
	}
	
void MSA112_GetCalibration(tagGroup MatrixTags,number &scale,number &origin, string &units,string axis)
	{
	string MSALabel = MSA112_Label()

	if(Axis=="E")	//for compatibility with previous versions
		{
		units  =MatrixTags.MSA112_getStringFromTags(MSAlabel+":energy units")
		scale  =MatrixTags.MSA112_getNumberFromTags(MSAlabel+":energy dispersion","F")
		origin =MatrixTags.MSA112_getNumberFromTags(MSAlabel+":energy origin","L")
		}
	else	
		{
		units  =MatrixTags.MSA112_getStringFromTags(MSAlabel+":units "+axis)
		scale  =MatrixTags.MSA112_getNumberFromTags(MSAlabel+":scale "+axis,"F")
		origin =MatrixTags.MSA112_getNumberFromTags(MSAlabel+":origin "+axis,"L")		
		}
	}

void MSA112_GetCalibration(image MSAimage,number &scale,number &origin, string &units,string axis)
	{
	TagGroup MatrixTags =  MSA112_GetMatrixTags(MSAimage,0)
	
	MSA112_GetCalibration(MatrixTags,scale,origin,units,axis)
	}





	
/******************************************************************
functions: Average Spectrum (always real single-row calibrated image
******************************************************************/

/**************  GLOBAL average spectrum **********************/
image MSA112_GetAverageSpectrum(TagGroup MatrixTags)
	{
	string MSALabel = MSA112_Label()
	
	number dispersion,OriginE
	string units	
	MSA112_GetCalibration(MatrixTags,dispersion,originE,units,"E")	

	number Depth =MatrixTags.MSA112_GetNumberFromTags(MSALabel+":energy depth","L")
	
	image AverSpec := realimage("average spectrum",4,depth,1)
	MatrixTags.TagGroupGetTagAsArray(MSALabel+":Average Spectrum",AverSpec)	
	
	AverSpec.ImageSetDimensionCalibration(0,originE,dispersion,units,1)

	return AverSpec
	}


/********************  from MSAimage ********************************/
image MSA112_GetAverageSpectrum(image MSAimage)
	{
	TagGroup MatrixTags =  MSA112_GetMatrixTags(MSAimage,0)
	
	return MSA112_GetAverageSpectrum(MatrixTags)
	}

/*********************** SET ***********************************/
void MSA112_SetAverageSpectrum(TagGroup MatrixTags,image AverSpec)
	{
	string MSALabel = MSA112_Label()
	
	MatrixTags.TagGroupSetTagAsArray(MSALabel+":Average Spectrum",AverSpec)	
	}
	
/********************* in MSA image *****************/	
void MSA112_SetAverageSpectrum(image MSAimage,image AverSpec)
	{												//GLOBAL
	TagGroup MatrixTags =  MSA112_GetMatrixTags(MSAimage,0)
	
	MSA112_SetAverageSpectrum(MatrixTags,AverSpec)
	}
	
	

/**************  CLUSTERED average spectrum **********************/
image MSA112_GetAverageSpectrum(TagGroup MatrixTags,number currCluster)
	{
	string MSALabel = MSA112_Label()
	string 	clusterLabel = MSA112_GetClusterLabel(currCluster)	
	
	number dispersion,OriginE
	string units	
	MSA112_GetCalibration(MatrixTags,dispersion,originE,units,"E")	

	number Depth =MatrixTags.MSA112_GetNumberFromTags(MSALabel+":energy depth","L")	

	image AverSpec := realimage("average spectrum ("+currCluster+")",4,depth,1)
	MatrixTags.TagGroupGetTagAsArray(ClusterLabel+":Average Spectrum",AverSpec)	//LOCAL !!
	
	AverSpec.ImageSetDimensionCalibration(0,originE,dispersion,units,1)
		
	return AverSpec
	}
	
/********************** SET *********************/
void MSA112_SetAverageSpectrum(TagGroup MatrixTags,image AverSpec,number currCluster)
	{													//LOCAL !!
	string MSALabel = MSA112_Label()
	string 	clusterLabel = MSA112_GetClusterLabel(currCluster)	
	
	MatrixTags.TagGroupSetTagAsArray(ClusterLabel+":Average Spectrum",AverSpec)	
	}

/**************  CLUSTERED METHOD-DEPENDENT average spectrum (probably shifted) ************/
image MSA112_GetAverageSpectrum(TagGroup MatrixTags,number currCluster,string Method)
	{
	string MSALabel = MSA112_Label()
	string 	clusterLabel = MSA112_GetClusterLabel(currCluster)

	number Depth =MatrixTags.MSA112_getNumberFromTags(MSAlabel+":energy depth","L")		
	string units =MatrixTags.MSA112_getStringFromTags(MSAlabel+":energy units")
	number dispersion =MatrixTags.MSA112_getNumberFromTags(MSAlabel+":energy dispersion","F")
	number OriginE =MatrixTags.MSA112_getNumberFromTags(MSAlabel+":energy origin","L")	
	

	image AverSpec := realimage("average spectrum "+Method+" ("+currCluster+")",4,Depth,1)
	if(!MatrixTags.TagGroupGetTagAsArray(ClusterLabel+":"+Method+":Average Spectrum",AverSpec) ) //LOCAL METHOD-dependent
		{
		MatrixTags.TagGroupGetTagAsArray(ClusterLabel+":Average Spectrum",AverSpec)
		}
		
	return AverSpec
	}
	
/******************* Set **************************/
void MSA112_SetAverageSpectrum(TagGroup MatrixTags,image AverSpec,number currCluster,string Method)
	{
	string MSALabel = MSA112_Label()
	string 	clusterLabel = MSA112_GetClusterLabel(currCluster)	
	
	MatrixTags.TagGroupSetTagAsArray(ClusterLabel+":"+Method+":Average Spectrum",AverSpec) //LOCAL METHOD-dependent
	}	
	
	
	
	
/******************************************************************
functions: Average Image
******************************************************************/

/**************  GLOBAL average image (single column)**********************/
image MSA112_GetAverageImage(TagGroup MatrixTags)
	{
	string MSALabel = MSA112_Label()
	
	number Points =MatrixTags.MSA112_getNumberFromTags(MSAlabel+":pixels","L")

	image AvIm := realimage("Average image",4,1,Points)
	MatrixTags.TagGroupGetTagAsArray(MSALabel+":Average Image",AvIm)	

	return AvIm
	}


/********************  from MSAimage ********************************/
image MSA112_GetAverageImage(image MSAimage)
	{
	TagGroup MatrixTags =  MSA112_GetMatrixTags(MSAimage,0)
	
	return MSA112_GetAverageImage(MatrixTags)
	}
	
	
/**************  GLOBAL average image 2D calibrated **********************/
image MSA112_GetAverageImage2D(TagGroup MatrixTags)
	{
	string MSALabel = MSA112_Label()

	image AvIm = MSA112_GetAverageImage(MatrixTags)		

	number width,height
	MatrixTags.MSA112_GetSIWidthHeight(width,height)
	
	number ScaleX,ScaleY,originX,originY
	string unitsX,unitsY
	MatrixTags.MSA112_GetCalibration(ScaleX,originX,unitsX,"X")
	MatrixTags.MSA112_GetCalibration(ScaleY,originY,unitsY,"Y")		
	
	image AvIm2D = MSA112_ConvertColumnTo2D(AvIm,width,height)
	
	ImageSetDimensionCalibration(AvIm2D,0,originX,scaleX,unitsX,1)
	ImageSetDimensionCalibration(AvIm2D,1,originY,scaleY,unitsY,1)

	return AvIm2D
	}


/********************  from MSAimage ********************************/
image MSA112_GetAverageImage2D(image MSAimage)
	{
	TagGroup MatrixTags =  MSA112_GetMatrixTags(MSAimage,0)
	
	return MSA112_GetAverageImage2D(MatrixTags)
	}
	

/******************Same for SET ********************************/
void MSA112_SetAverageImage(TagGroup MatrixTags,image AvIm)
	{
	string MSALabel = MSA112_Label()

	MatrixTags.TagGroupSetTagAsArray(MSALabel+":Average Image",AvIm)	
	}

/******************* in MSA image *******************/
void MSA112_SetAverageImage(image MSAimage,image AvIm)
	{
	TagGroup MatrixTags =  MSA112_GetMatrixTags(MSAimage,0)
	
	MSA112_SetAverageImage(MatrixTags,AvIm)
	}
	
	
	
/******************************************************************
Functions: number of components
******************************************************************/

number MSA112_ReadNumberComponents(tagGroup MatrixTags,number currCluster,string Readmethod)
{
string 	clusterLabel = MSA112_GetClusterLabel(currCluster)
number CompN =MatrixTags.MSA112_getNumberFromTags(ClusterLabel+":"+readMethod+":number components","S")

return CompN
}

/************************* from MSA image *****************************/
number MSA112_ReadNumberComponents(image MSAimage,number currCluster,string Readmethod)
{
TagGroup MatrixTags =  MSA112_GetMatrixTags(MSAimage,0)

return MSA112_ReadNumberComponents(MatrixTags,currCluster,Readmethod)
}

/************************* from MSA image, PCA only***********************/
number MSA112_ReadNumberComponents(image MSAimage,number currCluster)
{
string Readmethod = "PCA"

return MSA112_ReadNumberComponents(MSAimage,currCluster,Readmethod)
}

/************************* from MSA image, 0-cluster only***********************/
number MSA112_ReadNumberComponents(image MSAimage,string Readmethod)
{
number currCluster = 0

return MSA112_ReadNumberComponents(MSAimage,currCluster,Readmethod)
}

/************************* from MSA image, PCA only, 0-cluster only*****************/
number MSA112_ReadNumberComponents(image MSAimage)
{
string Readmethod = "PCA"
number currCluster = 0

return MSA112_ReadNumberComponents(MSAimage,currCluster,Readmethod)
}





	
	
	
/******************************************************************
Functions: Loadings		columns, always real images UNcalibrated
******************************************************************/
image MSA112_GetLoadings(TagGroup MatrixTags,string Method,number currCluster)
	{
	string 	clusterLabel = MSA112_GetClusterLabel(currCluster)
	
	image Pall 
	Pall :=MatrixTags.MSA112_getImageFromTags(ClusterLabel+":"+Method+":Loadings")

	if(!Pall.ImageIsValid())	//reserve path for old versions
		{
result(" reserve path for older version");		
		string MSALabel = MSA112_Label()
		number Depth	
		MatrixTags.TagGroupGetTagAsLong(MSAlabel+":energy depth",Depth)		//GLOBAL !!
		number CompN = MSA112_ReadNumberComponents(MatrixTags,currCluster,Method)
		
		Pall := realimage("Pvec all",4,CompN,depth)	//GLOBAL !!

		if(!MatrixTags.TagGroupGetTagAsArray(ClusterLabel+":"+Method+":Loadings",Pall)) Pall :=NULL
		}
	
	//Pall.matrixTranspose().showimage()	
	
	return Pall
	}
	

/********************  from MSAimage *************************************/
image MSA112_GetLoadings(image MSAimage,string Method,number currCluster)
	{
	TagGroup MatrixTags =  MSA112_GetMatrixTags(MSAimage,0)
	if(!MatrixTags.TagGroupIsValid()) return NULL
	
	return MSA112_GetLoadings(MatrixTags,Method,currCluster)
	}


/********************  from MSAimage, PCA only***************************************/
image MSA112_GetLoadings(image MSAimage,number currCluster)
	{
	string Method = "PCA"
	
	return MSA112_GetLoadings(MSAimage,Method,currCluster)
	}
	
	
/********************  from MSAimage, 0-cluster only ********************************/
image MSA112_GetLoadings(image MSAimage,string Method)
	{
	number currCluster = 0
	
	return MSA112_GetLoadings(MSAimage,Method,currCluster)
	}

	
/********************  from MSAimage, PCA-only, 0-cluster only ***********************/
image MSA112_GetLoadings(image MSAimage)
	{
	string Method = "PCA"
	number currCluster = 0
	
	return MSA112_GetLoadings(MSAimage,Method,currCluster)
	}
	
/********************  same for SET *************************************/
void MSA112_SetLoadings(TagGroup MatrixTags,image Pall,string Method,number currCluster)
	{
	string 	clusterLabel = MSA112_GetClusterLabel(currCluster)

	MatrixTags.MSA112_setImageinTags(ClusterLabel+":"+Method+":Loadings",Pall)
	}

/********************  in MSAimage *************************************/
void MSA112_SetLoadings(image MSAimage,image Pall,string Method,number currCluster)
	{
	TagGroup MatrixTags =  MSA112_GetMatrixTags(MSAimage,0)
	if(MatrixTags.TagGroupIsValid())
		{
		MSA112_SetLoadings(MSAimage,Pall,Method,currCluster)
		}
	}

/********************  in MSAimage, PCA only***************************************/
void MSA112_SetLoadings(image MSAimage,image Pall,number currCluster)
	{
	string Method = "PCA"
	MSA112_SetLoadings(MSAimage,Pall,Method,currCluster)
	}

/******************** in MSAimage, 0-cluster only ********************************/
void MSA112_SetLoadings(image MSAimage,image Pall,string Method)
	{
	number currCluster = 0
	MSA112_SetLoadings(MSAimage,Pall,Method,currCluster)
	}

/********************  in MSAimage, PCA-only, 0-cluster only ***********************/
void MSA112_SetLoadings(image MSAimage,image Pall)
	{
	string Method = "PCA"
	number currCluster = 0
	MSA112_SetLoadings(MSAimage,Pall,Method,currCluster)
	}
	
	
/********************  ONE component, in MSAimage, 0-cluster only ***********************/
image MSA112_GetLoading(image MSAimage, number comp,string Method)
{
image loadings = MSA112_GetLoadings(getFrontImage(),Method)

number components,channels
loadings.getSize(components,channels)

if(comp<1) return NULL
if(comp>components) return NULL

return loadings[0,comp-1,channels,comp].MatrixTranspose()
}

/********************  ONE component, in MSAimage,PCA-only, 0-cluster only ***********************/
image MSA112_GetLoading(image MSAimage, number comp)
{
string Method = "PCA"

return MSA112_GetLoading(MSAimage,comp,Method)
}



	
/******************************************************************
Functions: Scores		columns, always real images UNcalibrated
******************************************************************/

image MSA112_GetScores(TagGroup MatrixTags,string Method,number currCluster)
	{
	string 	clusterLabel = MSA112_GetClusterLabel(currCluster)

	image Tall
	Tall :=MatrixTags.MSA112_getImageFromTags(ClusterLabel+":"+Method+":Scores")
	
	if(!Tall.ImageIsValid())	//reserve path for old versions
		{
result(" reserve path for older version");
		number kPoints
		MatrixTags.TagGroupGetTagAsLong(ClusterLabel+":points",kPoints)
		
		number CompN = MatrixTags.MSA112_ReadNumberComponents(currCluster,Method)
		
		Tall:= realimage("T-vectors",4,CompN,kPoints)
		if(!MatrixTags.TagGroupGetTagAsArray(ClusterLabel+":"+Method+":Scores",Tall)) Tall :=NULL		
		}
		
	return Tall
	}


/********************  from MSAimage *************************************/
image MSA112_GetScores(image MSAimage,string Method,number currCluster)
	{
	TagGroup MatrixTags =  MSA112_GetMatrixTags(MSAimage,0)
	if(!MatrixTags.TagGroupIsValid()) return NULL
	
	return MSA112_GetScores(MatrixTags,Method,currCluster)
	}

	
/********************  from MSAimage, PCA only***************************************/
image MSA112_GetScores(image MSAimage,number currCluster)
	{
	string Method = "PCA"
	
	return MSA112_GetScores(MSAimage,Method,currCluster)
	}


/********************  from MSAimage, 0-cluster only ********************************/
image MSA112_GetScores(image MSAimage,string Method)
	{
	number currCluster = 0
	
	return MSA112_GetScores(MSAimage,Method,currCluster)
	}
	
	
/********************  from MSAimage, PCA-only,0th-cluster only ********************************/

image MSA112_GetScores(image MSAimage)
	{
	string Method = "PCA"
	number currCluster = 0
	
	return MSA112_GetScores(MSAimage,Method,currCluster)
	}
	

/********************  from MSAimage one score only ***********************/

image MSA112_GetScore(image MSAimage,number comp,string Method)
	{
	image Scores =MSA112_GetScores(MSAimage,Method)
	
	number components,points
	Scores.getSize(components,points)

	if(comp<1) return NULL
	if(comp>components) return NULL
	
	return Scores[0,comp-1,points,comp].MatrixTranspose()
	}

image MSA112_GetScore(image MSAimage,number comp)
	{
	string Method ="PCA"

	return MSA112_GetScore(MSAimage,comp,Method)
	}
	
image MSA112_GetScore2D(image MSAimage,number comp)
	{
	TagGroup MatrixTags = MSA112_GetMatrixTags(MSAimage,0)
	string MSAlabel = MSA112_Label()

	number Points =MatrixTags.MSA112_getNumberFromTags(MSAlabel+":pixels","L")
	number width =MatrixTags.MSA112_getNumberFromTags(MSAlabel+":width","L")
	number height =MatrixTags.MSA112_getNumberFromTags(MSAlabel+":height","L")

	image Scores = MSAimage.MSA112_GetScores()
		
	return Scores[0,comp-1,Points,comp].MSA112_ConvertColumnTo2D(width,height)
	}	


	
	
/********************  same for SET *************************************/
void MSA112_SetScores(TagGroup MatrixTags,image Tall,string Method,number currCluster)
	{
	string 	clusterLabel = MSA112_GetClusterLabel(currCluster)

	MatrixTags.MSA112_setImageInTags(ClusterLabel+":"+Method+":Scores",Tall)
	}

/********************  in MSAimage *************************************/	
void MSA112_SetScores(image MSAimage,image Tall,string Method,number currCluster)
	{
	string 	clusterLabel = MSA112_GetClusterLabel(currCluster)

	TagGroup MatrixTags =  MSA112_GetMatrixTags(MSAimage,0)
	if(MatrixTags.TagGroupIsValid())
		{
		MSA112_SetScores(MatrixTags,Tall,Method,currCluster)
		}
	}

/********************  in MSAimage, PCA only***************************************/
void MSA112_SetScores(image MSAimage,image Tall,number currCluster)
	{
	string Method = "PCA"
	
	MSA112_SetScores(MSAimage,Tall,Method,currCluster)
	}

/********************  in MSAimage, 0-cluster only ********************************/
void MSA112_SetScores(image MSAimage,image Tall,string Method)
	{
	number currCluster = 0
	
	MSA112_SetScores(MSAimage,Tall,Method,currCluster)
	}

/********************  in MSAimage, PCA-only,0th-cluster only ***********************/
void MSA112_SetScores(image MSAimage,image Tall)
	{
	string Method = "PCA"
	number currCluster = 0
	
	MSA112_SetScores(MSAimage,Tall,Method,currCluster)
	}




/******************************************************************
Functions: EndMembers		rows, always real images UNcalibrated
******************************************************************/
image MSA112_GetEMSpectra(TagGroup MatrixTags,number currCluster)
	{
	string 	clusterLabel = MSA112_GetClusterLabel(currCluster)
	
	image EMSSpectra
	string Emembering ="EndMembering"
	EMSSpectra :=MatrixTags.MSA112_getImageFromTags(ClusterLabel+":"+Emembering+":EMSpectra")

	
	//Pall.matrixTranspose().showimage()	
	
	return EMSSpectra
	}
	
/********************  from MSAimage *************************************/
image MSA112_GetEMSpectra(image MSAimage,number currCluster)
	{
	TagGroup MatrixTags =  MSA112_GetMatrixTags(MSAimage,0)
	if(!MatrixTags.TagGroupIsValid()) return NULL
	
	return MSA112_GetEMSpectra(MatrixTags,currCluster)
	}

/********************  from MSAimage, 0-cluster only ********************************/
image MSA112_GetEMSpectra(image MSAimage)
	{
	number currCluster = 0
	
	return MSA112_GetEMSpectra(MSAimage,currCluster)
	}

/******************************************************************
Functions: EndMembers Abundancies	columns, always real images UNcalibrated
******************************************************************/
image MSA112_GetEMAbundancies(TagGroup MatrixTags,number currCluster)
	{
	string 	clusterLabel = MSA112_GetClusterLabel(currCluster)
	
	image Abundancies
	string Emembering ="EndMembering"
	Abundancies :=MatrixTags.MSA112_getImageFromTags(ClusterLabel+":"+Emembering+":Abundancies")

	
	//Pall.matrixTranspose().showimage()	
	
	return Abundancies
	}
	
/********************  from MSAimage *************************************/
image MSA112_GetEMAbundancies(image MSAimage,number currCluster)
	{
	TagGroup MatrixTags =  MSA112_GetMatrixTags(MSAimage,0)
	if(!MatrixTags.TagGroupIsValid()) return NULL
	
	return MSA112_GetEMAbundancies(MatrixTags,currCluster)
	}

/********************  from MSAimage, 0-cluster only ********************************/
image MSA112_GetEMAbundancies(image MSAimage)
	{
	number currCluster = 0
	
	return MSA112_GetEMAbundancies(MSAimage,currCluster)
	}
	
	


/******************************************************************
Functions: EigenValues
******************************************************************/


/* *****************caclulate variance *************** */
number MSA112_Variance(image scoreP)
	{
	image score = scoreP.imageclone()
	
	score -=mean(score)

	number kPoints =score.ImageGetDimensionSize(1)	
		
	if(kPoints >1)	//column vector
		{
		return sum(MatrixMultiply(MatrixTranspose(score),score))/kPoints
		}
	if(kPoints ==1)	//row vector
		{ 
		kPoints =score.ImageGetDimensionSize(0)
		return sum(MatrixMultiply(score,MatrixTranspose(score)))/kPoints	
		} 
	}
	
	
/********* Calculate Kurtosis **********************************/
number MSA112_Kurtosis(image wTx2P)
{ //entrance must be SQUARED score
image wTx2 =wTx2P.imageclone()
wTx2 /= mean(wTx2)

return (mean(wTx2**2)-3)		
}


/********* Calculate EigenValues******************************/
image MSA112_ExtractEigenValues(TagGroup MatrixTags,number CurrCluster,string Method,string ValName)
	{
	image Tall :=MSA112_GetScores(MatrixTags,Method,currCluster)
	
	number CompN,kPoints
	Tall.getSize(CompN,kPoints)

	image Values := realimage("",4,CompN,1)
	for(number i=0;i<CompN;i++)
		{
		image Tone = Tall[0,i,kPoints,i+1]
		number Value
		if(ValName=="Variance") Value = MSA112_Variance(Tone)
		if(ValName=="Kurtosis") Value = MSA112_Kurtosis(Tone**2)

		Values.setPixel(i,0,Value)
		}
		
	return Values
	}
	
/********************  calculate eigenValues (PCA only, variance only *********************/
image MSA112_ExtractEigenValues(image MSAimage)
	{
	TagGroup MatrixTags =  MSA112_GetMatrixTags(MSAimage,0)
	string Method = "PCA"
	string ValName = "Variance"
	number currClust = 0
	
	return MSA112_ExtractEigenValues(MatrixTags,currClust,Method,ValName)
	}



/********************  get stored eigen values or calculate new ones *******************/
image MSA112_GetEigenValues(TagGroup MatrixTags,string ValName,number currCluster, string Method)
	{
	string 	clusterLabel = MSA112_GetClusterLabel(currCluster)
	image Values
	Values :=MatrixTags.MSA112_getImageFromTags(ClusterLabel+":"+Method+":"+ValName)
	
	if(!Values.ImageIsValid())	//reserve path for old versions
		{
result(" reserve path for older version");		
		number CompN = MatrixTags.MSA112_ReadNumberComponents(currCluster,Method)	
		
		Values := realimage("",4,CompN,1)	

		if(!MatrixTags.TagGroupGetTagAsArray(ClusterLabel+":"+Method+":"+ValName,Values))
			{
			//was not store even in the old format
			Values =MSA112_ExtractEigenValues(MatrixTags,CurrCluster,Method,ValName) 
			}
		}
		
	return Values
	}


/********************  from MSAimage ********************************/
image MSA112_GetEigenValues(image MSAimage,string ValName,number currCluster,string Method)
	{
	TagGroup MatrixTags = MSA112_GetMatrixTags(MSAimage,0)
	if(!MatrixTags.TagGroupIsValid()) return NULL
	
	return MSA112_GetEigenValues(MatrixTags,ValName,currCluster,Method)
	}


/********************  from MSAimage, PCA only********************************/
image MSA112_GetEigenValues(image MSAimage,string ValName,number currCluster)
	{
	string Method = "PCA"
	
	return MSA112_GetEigenValues(MSAimage,ValName,currCluster,Method)
	}


/********************  from MSAimage, 0-cluster only ********************************/
image MSA112_GetEigenValues(image MSAimage,string ValName,string Method)
	{
	number currCluster = 0
	
	return MSA112_GetEigenValues(MSAimage,ValName,currCluster,Method)
	}


/********************  from MSAimage, PCA-only, 0-cluster only ********************************/
image MSA112_GetEigenValues(image MSAimage,string ValName)
	{
	string Method = "PCA"
	number currCluster = 0
	
	return MSA112_GetEigenValues(MSAimage,ValName,currCluster,Method)
	}


/********************  from MSAimage, PCA-only, 0-cluster only, variance only *********************/
image MSA112_GetEigenValues(image MSAimage)
	{
	string ValName = "Variance"
	string Method = "PCA"
	number currCluster = 0
	
	return MSA112_GetEigenValues(MSAimage,ValName,currCluster,Method)
	}
	
	
/**************************Same for set ***********************/

void MSA112_SetEigenValues(TagGroup MatrixTags,image Values,string ValName,number currCluster, string Method)
	{
	number CompN = MatrixTags.MSA112_ReadNumberComponents(currCluster,Method)	
		
	string 	clusterLabel = MSA112_GetClusterLabel(currCluster)
	MatrixTags.MSA112_setImageinTags(ClusterLabel+":"+Method+":"+ValName,Values)
	}
	

/********************  in MSAimage ********************************/
void MSA112_SetEigenValues(image MSAimage,image Values,string ValName,number currCluster, string Method)
	{
	TagGroup MatrixTags =  MSA112_GetMatrixTags(MSAimage,0)
	
	MSA112_SetEigenValues(MatrixTags,Values,ValName,currCluster,Method)
	}
	
/********************  in MSAimage, PCA only********************************/
void MSA112_SetEigenValues(image MSAimage,image Values,string ValName,number currCluster)
	{
	string Method ="PCA"
	
	MSA112_SetEigenValues(MSAimage,Values,ValName,currCluster,Method)
	}
	
/********************  in MSAimage, 0-cluster only ********************************/
void MSA112_SetEigenValues(image MSAimage,image Values,string ValName,string Method)
	{
	number currCluster =0
	
	MSA112_SetEigenValues(MSAimage,Values,ValName,currCluster,Method)
	}
	
/********************  in MSAimage, PCA-only, 0-cluster only ********************************/
void MSA112_SetEigenValues(image MSAimage,image Values,string ValName)
	{
	number currCluster =0
	string Method ="PCA"
	
	MSA112_SetEigenValues(MSAimage,Values,ValName,currCluster,Method)
	}
	
/********************  in MSAimage, PCA-only, 0-cluster only, variance only *********************/
void MSA112_SetEigenValues(image MSAimage,image Values)
	{
	number currCluster =0
	string Method ="PCA"
	string ValName ="Variance"
	
	MSA112_SetEigenValues(MSAimage,Values,ValName,currCluster,Method)
	}
	
	
	
/******************************************************************
functions: MSA matrix 	(always real image)
******************************************************************/
image MSA112_GetMSAMatrix(TagGroup MatrixTags)
	{
	string MSALabel = MSA112_Label()
	
	number Points =MatrixTags.MSA112_getNumberFromTags(MSAlabel+":pixels","L")
	number Depth =MatrixTags.MSA112_getNumberFromTags(MSAlabel+":energy depth","L")
	if(Points==0) {OKdialog("number of pixels not stored");throw("");}
	if(Depth==0) {OKdialog("number of channels not stored");throw("");}
	
	image DataMatrix := realimage("",4,Depth,Points)	
	MatrixTags.TagGroupGetTagAsArray(MSAlabel+":Matrix",DataMatrix)
	
	return DataMatrix
	}

/********************  from MSAimage ********************************/
image MSA112_GetMSAMatrix(image MSAimage)
	{
	TagGroup MatrixTags =  MSA112_GetMatrixTags(MSAimage,1)	//full MSAimage needed
	
	return MSA112_GetMSAMatrix(MatrixTags)
	}
	
	
/********************* set MSA matrix ********************/
void MSA112_SetMSAMatrix(TagGroup MatrixTags,image DataMatrix)
	{
	string MSALabel = MSA112_Label()
	
	MatrixTags.TagGroupSetTagAsArray(MSAlabel+":Matrix",DataMatrix)
	}

/****************same in MSA image *****************/
void MSA112_SetMSAMatrix(image MSAimage,image DataMatrix)
	{
	TagGroup MatrixTags =  MSA112_GetMatrixTags(MSAimage,1)	//full MSAimage needed
	
	MSA112_SetMSAMatrix(MatrixTags,DataMatrix)
	}
	

	
/**************** cluster dependent ********************************/
image MSA112_getMSAMatrix(TagGroup MatrixTags,number currCluster)
	{
	//this may use compressed MSA, thus PCA in certain cluster needed
	//for the moment only one cluster reconstucted
	string MSAlabel = MSA112_Label()
	
	string type =MatrixTags.MSA112_getStringFromTags(MSA112_Label()+":type")
	if(type=="MSA compressed") if(!OKCancelDialog("Original data are not available\nThe PCA compressed data will be used for ALS")) exit(0)

	number Points =MatrixTags.MSA112_getNumberFromTags(MSAlabel+":pixels","L")
	number depth =MatrixTags.MSA112_getNumberFromTags(MSAlabel+":energy depth","L")	
	image Cluster:= integerImage("cluster",2,0,1,Points)
	MatrixTags.TagGroupGetTagAsArray(MSAlabel+":Cluster Mask",cluster)

	image DataMatrix,dataMatrixRed
	
	if(type=="MSA matrix")
		{
		DataMatrix = MSA112_GetMSAMatrix(MatrixTags)
		
		if(currCluster==0) dataMatrixRed := dataMatrix		
		else 
			{
			number kPoints
			dataMatrixRed = MSA112_TakeCluster(dataMatrix,cluster,currCluster,kPoints)
			}
		}
		
	if(type=="MSA compressed")	
		{	
		image Tall:= MSA112_GetScores(MatrixTags,"PCA",currCluster)
		image Pall:= MSA112_GetLoadings(MatrixTags,"PCA",currCluster)

		dataMatrixRed = matrixMultiply(Tall,MatrixTranspose(Pall))
		image AvSpecWeighted := MSA112_GetAverageSpectrum(MatrixTags,CurrCluster)	
		dataMatrixRed += MSA112_ExpandImageVertically(AvSpecWeighted,dataMatrixRed)
		//dataMatrixRed.showimage()
		}
		
	return dataMatrixRed
	}
	
	





/******************************************************************
Functions Check components
******************************************************************/
number MSA112_CheckValidityOfComponents(number First,number Second,number CompN,number CurrCluster)
	{
	if(First<=0||Second<=0) {OKdialog("component must have index greater then zero !!");return 0;} 
	if((First>CompN)||(Second>CompN)) {OKdialog("there are only "+CompN+" components\n calculated in cluster "+CurrCluster);return 0;}

	return 1 
	}

number MSA112_CheckNumberOfComponents(number CurrCluster,number &CompN,TagGroup MatrixTags,string Method)
	{
	//return 1 if there are PCA components in requested  cluster
	//return 0 if something misfit
	//return (-1) if there are many cluster, but requested 0 (forget about clustering?)
	
	//check whether input cluster makes sense
	if(Currcluster<(-1)) {OKdialog("Enter the valid cluster number !!");return 0;}
	if(Currcluster==(-1)) {OKdialog("Batch treatment is not available for this operation !!");return 0;} 

	//check if there is requested Method in cluster	
	string 	clusterLabel = MSA112_GetClusterLabel(currCluster)
	CompN =MatrixTags.MSA112_getNumberFromTags(ClusterLabel+":"+Method+":number components","S") 

	if(CompN==0) {OKdialog("there are no "+Method+" components in cluster "+currCluster+" !!");return 0;} 

	//ensures that user ignores clustering
	number clusterN =MatrixTags.MSA112_getNumberFromTags(MSA112_Label()+":clusters number","S") 

	if(clusterN!=0&&Currcluster==0) return (-1)
			
	return 1
	}



// this is needed in display loads
image MSA112_GetComponents(TagGroup MatrixTags,string ReadMethod,number currCluster)
	{
	string 	clusterLabel = MSA112_GetClusterLabel(currCluster)

	number CompN = MSA112_ReadNumberComponents(MatrixTags,currCluster,ReadMethod)
	image Components :=MatrixTags.MSA112_getImageFromTags(ClusterLabel+":"+ReadMethod+":components")

	if(!Components.ImageIsValid())	Components :=NULL
	
	return Components
	}


	
/******************************************************************
Functions related to component list
******************************************************************/


/************** make component list **********************/
TagGroup MSA112_MakeComponentList(image Allcomponent,TagGroup &ComponentList,number Cluster0,number ClusterF,number kDim,string action)
	{	
	TagGroup DialogList =  DLGCreateDialog(Action).DLGTableLayout(ClusterF-Cluster0+1,1,0)	

	for(number CurrCluster = cluster0;CurrCluster <=ClusterF;CurrCluster++)
		{
		TagGroup ClusterBox = DLGCreateBox("cluster "+CurrCluster)
		for(number j = 1;j<=kDim;j++)
			{
			TagGroup Check
			number Value = Allcomponent.getPixel(CurrCluster,j)
			//0 - not existing/disabled , 1 - PCA

			if(value!=0)
				{			
				string Caption = Action+" "
				if(Value==1) Caption +="PCA" 
				
				Caption += (" "+j)				
				Check= DLGCreateCheckBox(Caption,value).DLGanchor("west")
				
				}
			else
				{
				Check = DLGCreateLabel("")
				}
			ComponentList.TagGroupAddTagGroupAtEnd(Check)
			//components lists needs for reading the choice
			ClusterBox.DLGAddElement(Check)
			}
			
		DialogList.DLGAddElement(ClusterBox)
		}
	
	return DialogList
	}	



/************** check if enough components to rotate **********************/
number MSA112_CheckComponents(number CompN)
{
if(CompN<=0)
	{
	OKdialog("No components found !!")
	return 0;
	}
if(CompN<2) 
	{
	OKdialog("At least 2 components should exist for this operation !!")
	return 0;
	}
else return 1
}


/************** confirm components **********************/
image MSA112_ConfirmComponents(number Cluster0,number ClusterF,number kDim,TagGroup MatrixTags,string MSAlabel,string action)
{
//find max number of components among clusters
number mMax
for(number Clust=Cluster0;Clust<=ClusterF;Clust++)
	{
	string ClusterLabel = MSA112_GetClusterLabel(Clust)

	number CompN = MSA112_ReadNumberComponents(MatrixTags,Clust,"PCA")
	if(CompN==0) OKDialog("PCA decomposition in cluster "+Clust+ " was not performed !!")
	if(CompN>mMax) mMax =CompN
result("\ncluster "+Clust+" m "+CompN+" mMax "+mMax)
	}
kDim = min(kDim,mMax)	//max extracted number of components or input number

image Allcomponents:=integerimage("",1,0,ClusterF+1,kDim+1) 

//read exsisting info about components
for(number Clust=Cluster0;Clust<=ClusterF;Clust++)
	{
	string ClusterLabel = MSA112_GetClusterLabel(Clust)

	number CompN = MSA112_ReadNumberComponents(MatrixTags,Clust,"PCA")
	image Components := integerimage("",1,0,1,CompN+1)
	Components[1,0,CompN+1,1] = 1 //all PCA components are default to be rotated

	if(kDim<CompN) AllComponents[0,Clust,kDim+1,Clust+1]  = Components[0,0,kDim+1,1]	//only kDim requested components)
	if(kDim>=CompN) AllComponents[0,Clust,CompN+1,Clust+1] = Components[0,0,CompN+1,1]	//only CompN available components
	}
	
TagGroup ComponentList = newTagList()
object CompListFrame = Alloc(UIframe).init(MSA112_MakeComponentList(AllComponents,ComponentList,Cluster0,ClusterF,kDim,action))			
if(CompListFrame.Pose()) 
	{ 
	for(number CurrCluster = cluster0;CurrCluster <= ClusterF;CurrCluster++)
		{
		for(number j = 1;j<=kDim;j++)
			{
			TagGroup Catched 
			ComponentList.TagGroupGetIndexedTagAsTagGroup((CurrCluster-Cluster0)*kDim+j-1,Catched)
			number Flag
			if(catched.DLGgetValue(Flag))	//if component not exist then always "0"
				{
				if(Flag==0)
					{	
					AllComponents.setPixel(CurrCluster,j,0) 	//not use/not rotate
					} 
				}
			}
		}
	}
else {result("\ncanceled");throw("canceled");}

return Allcomponents
}


image MSA112_GetComponentList(TagGroup MatrixTags)
{
string MSAlabel =MSA112_Label()
image AllComponents :=MatrixTags.MSA112_getImageFromTags(MSAlabel+":Components List")

if(!AllComponents.ImageIsValid())	//reserve path for old versions
	{
	result(" reserve path for older version");
	number Cwidth,Cheight
	MatrixTags.TagGroupGetTagAsShort(MSAlabel+":components width",Cwidth)
	MatrixTags.TagGroupGetTagAsShort(MSAlabel+":components height",Cheight)

	if((Cwidth>=1)&&(Cheight>=1)) AllComponents :=integerimage("",1,0,Cwidth,Cheight)
	else return NULL
	
	if(!MatrixTags.TagGroupGetTagAsArray(MSAlabel+":Components list",AllComponents)) 
		{
		AllComponents :=NULL	
		}
	}
return AllComponents
}



/******************************************************************
Function Flip components
******************************************************************/
void MSA112_FlipComponent(TagGroup MatrixTags, number Comp,string Method,number currCluster)
	{
	image Pall = MSA112_GetLoadings(MatrixTags,Method,currCluster)
	image Tall = MSA112_GetScores(MatrixTags,Method,currCluster)

	number Points,depth,CompN
	Pall.getSize(CompN,depth)
	Tall.getSize(CompN,Points)

	if(Comp<1) throw("Component not valid !!")
	if(Comp>CompN) throw("Requested component not extracted !!")
		
	Pall[0,Comp-1,depth,Comp] *= (-1) 
	Tall[0,Comp-1,Points,Comp] *= (-1) 

	MSA112_SetLoadings(MatrixTags,Pall,Method,currCluster)
	MSA112_SetScores(MatrixTags,Tall,Method,currCluster)
	}

void MSA112_FlipComponent(image MSAimage, number Comp,string Method,number currCluster)
	{
	TagGroup MatrixTags =  MSA112_GetMatrixTags(MSAimage,0)
	if(MatrixTags.TagGroupIsValid())
		{
		MSA112_FlipComponent(MatrixTags,Comp,Method,currCluster)
		}
	}

void MSA112_FlipComponent(image MSAimage, number Comp,number currCluster)
	{
	string Method = "PCA"
	
	MSA112_FlipComponent(MSAimage,Comp,Method,currCluster)
	}
	
void MSA112_FlipComponent(image MSAimage, number Comp,string Method)
	{
	number currCluster = 0
	
	MSA112_FlipComponent(MSAimage,Comp,Method,currCluster)
	}

void MSA112_FlipComponent(image MSAimage, number Comp)
	{
	string Method = "PCA"
	number currCluster = 0
	
	MSA112_FlipComponent(MSAimage,Comp,Method,currCluster)
	}

















/*********************					*
**************************				**
for Advanced         **********************
**************************				**
*********************					*/












/******************************************************************
Functions: basic lowCut command
******************************************************************/
image MSA112_LBfilter(image Data,image AvSpecIn,number bond)
{
//AvSpecIn must have same dimension as Data
image filtered = imageclone(Data)
filtered = tert(Data<AvSpecIn*bond,0,Data)

return filtered
}



/******************************************************************
Functions: Spikes code
******************************************************************/
/**************** encode spikes code *************/
number MSA112_Encode(number OutTreat,number OutCorr,number shiftCor,number Aver)
{
//OutTreat = 0: do not treat outliers		
//OutTreat = 1: treat score outliers
//OutTreat = 2: treat loading outliers
//OutCorr = 1 : remove by shifting pixels
//OutCorr = 2 : remove by shifting channels
//ShiftCorr : shifting neighbourhood to correct a spike 1...9 (default 1)
//Aver : Gaussian smooth in extracting smooth loading 2...99(default 7)

number SpikesCode
if(OutTreat==0) SpikesCode=0	//first digit is ON/OFF swich
else SpikesCode = Aver*1000+ShiftCor*100+OutCorr*10+OutTreat

return SpikesCode
}

/**************** decode spikes code *************/
number MSA112_Decode(number &Code, number Coding)
{
//it is decoded sequentially, from right to left, each time code is reduced
number decoded = Code -(trunc(Code/Coding))*Coding
Code = trunc(Code/Coding)

return decoded
}









/*********************					*
**************************				**
Filtering         *************************
**************************				**
*********************					*/

	
	
	
	
	
	
	
	

	
	
/******************************************************************
Function Transfer info from image to image
******************************************************************/
void MSA112_TransferCalibration(image cube,image &final,number binX,number binY,number clipE)	
	{
	number origin,originX,originY,disp,scaleX,scaleY
	string units,unitsX,unitsY
	ImageGetDimensionCalibration(cube,0,originX,scaleX,unitsX,1)
	ImageGetDimensionCalibration(cube,1,originY,scaleY,unitsY,1)
	ImageGetDimensionCalibration(cube,2,origin,disp,units,1)

	scaleX *= binX
	originX /= binX
	scaleY *= binY
	originY /= binY
	
	origin -= clipE
	
	ImageSetDimensionCalibration(final,0,originX,scaleX,unitsX,1)
	ImageSetDimensionCalibration(final,1,originY,scaleY,unitsY,1)
	ImageSetDimensionCalibration(final,2,origin,disp,units,1)
	}

TagGroup MSA112_TransferImageTags(image in,image &out)	
	{
	TagGroup InTags = in.ImageGetTagGroup()
	TagGroup OutTags = out.ImageGetTagGroup()
	TagGroupCopyTagsFrom(OutTags,InTags)
	
	return OutTags
	}
	
	
/******************************************************************
Function Clipping
******************************************************************/
image MSA112_Xclipping(image cube,number clip)
	{
	number width,height,depth
	cube.get3Dsize(width,height,depth)
	
	return cube[clip,clip,0,width-clip,height-clip,Depth]		
	}
	
image MSA112_Eclipping(image cube,number clip)
	{
	number width,height,depth
	cube.get3Dsize(width,height,depth)
	
	return cube[0,0,clip,width,height,Depth-clip]
	}	

image MSA112_binning(image cube, number binX, number binY)
	{
	number running,steps,stepped
	object frm
	
	image binned = MSA112_binning(frm,cube,binX,binY,running,steps,stepped)
			
	MSA112_TransferCalibration(cube,binned,binX,binY,0)
	TagGroup BinnedTags = MSA112_TransferImageTags(cube,binned)
	
	return binned
	}

image MSA112_Xfiltering(image cube, number sigma,number minTail)
	{
	number running,steps,stepped
	object frm
	
	image filtered = frm.MSA112_Xfiltering(cube,sigma,minTail,running,steps,stepped)
	
	MSA112_TransferCalibration(cube,filtered,1,1,0)
	TagGroup FinalTags = MSA112_TransferImageTags(cube,filtered)
	
	return filtered
	}

image MSA112_Efiltering(image cube, number sigma,number minTail)
	{
	number running,steps,stepped
	object frm
	
	image filtered = frm.MSA112_Efiltering(cube,sigma,minTail,running,steps,stepped)
	
	MSA112_TransferCalibration(cube,filtered,1,1,0)
	TagGroup FinalTags = MSA112_TransferImageTags(cube,filtered)
		
	return filtered	
	}











/*********************					*
**************************				**
SCATTER PLOTS *****************************
**************************				**
*********************					*/

	
	
	
	
	
	
	
	
	
	
	
	
	
/******************************************************************
functions: Cluster Colors
******************************************************************/
rgbnumber  MSA112_ClusterColor(number currCluster)
{
rgbnumber col
if(currCluster==0) col = rgb(0,0,0)			//black
if(currCluster==1) col = rgb(255,0,0)		//red
if(currCluster==2) col = rgb(0,255,0)		//green
if(currCluster==3) col = rgb(0,0,255)		//blue
if(currCluster==4) col = rgb(255,0,255) 	//pink
if(currCluster==5) col = rgb(255,255,0) 	//yellow
if(currCluster==6) col = rgb(0,255,255) 	//lightblue

if(currCluster>6) 
	{
	//color cyclically repeated
	currCluster = mod(currCluster-1,6) +1
	
	return MSA112_ClusterColor(currCluster)
	}

return col
}

number  MSA112_GetClusterNo(string Color)
{
number  CurrCluster
if(Color=="red") CurrCluster = 1
if(Color=="green") CurrCluster = 2
if(Color=="blue") CurrCluster = 3
if(Color=="pink") CurrCluster = 4
if(Color=="yellow") CurrCluster = 5
if(Color=="lightblue") CurrCluster = 6

return CurrCluster
}

string  MSA112_ColorString(number currCluster)
{
string Color 
if(currCluster==1) Color = "red"
if(currCluster==2) Color = "green"
if(currCluster==3) Color = "blue"
if(currCluster==4) Color = "pink"
if(currCluster==5) Color = "yellow"
if(currCluster==6) Color = "lightblue"


if(currCluster>6) 
	{
	//color cyclically repeated
	currCluster = mod(currCluster-1,6) +1
	
	return MSA112_ColorString(currCluster)
	}
	
return Color
}

//Currently not is use
/*
string  MSA112_ClusterColorList(number Subclusters)
{
string ClusterList ="choose cluster to correct:\n"
for(number i=1;i<=Subclusters;i++) ClusterList += (MSA112_ColorString(i)+"\n")

return ClusterList
}
*/


/******************************************************************
Functions:  plot Scatter Plots 
******************************************************************/
/**************paint one cluster **********************/
void  MSA112_PaintCluster(rgbimage Scatter,image score1,image score2, number currCluster)
	{	
	rgbnumber col = MSA112_ClusterColor(currCluster)
	
	Index(Scatter,score1,score2) = col	
	}


/************** Repaint Scatter Plot ************/
void  MSA112_RepaintCluster(rgbimage Scatter,image cluster2D,image cXini,image cYini,number currCluster)
	{
	//reduce to current cluster
	image cXCol := imageclone(cXini)
	image cYCol := imageclone(cYini)
	cXCol = tert(cluster2D==currCluster,cXCol,0)
	cYCol = tert(cluster2D==currCluster,cYCol,0)
	
	rgbnumber col = MSA112_ClusterColor(currCluster)
		
	Index(Scatter,cXCol,cYCol) = col
	Scatter.SetPixel(0,0,rgb(256,256,256)) 	//restore zero pixel in white
	}
	

	
/**************Find range for scatter plot**********************/	

void MSA112_SetPlainMaximumScatterPlot(number PlainMax) 
	{
	string parName ="PlainMax"
	setPersistentNumberNote(MSA112_MSApath()+":"+parName,PlainMax)
	}

number MSA112P_GetPlainMaximumScatterPlot()
	{
	string parName ="PlainMax"
	number PlainMax
	getPersistentNumberNote(MSA112_MSApath()+":"+parName,PlainMax)	
	
	return PlainMax
	}
	
void MSA112_FindMinMaxPlot(image score,number &mini,number &maxi)
	{
	number stDev = 5	//5 std deviations range
	number aver = mean(score)
	number sigma = sqrt(MSA112_Variance(score))

	number minV = min(score)
	number maxV = max(score)
	number minS = aver - stDev*sigma
	number maxS = aver + stDev*sigma
//result("\nminV "+minV+" minS "+minS+" maxV "+maxV+" maxS "+maxS)
	
	number PlainMax =MSA112P_GetPlainMaximumScatterPlot()
	
	if(PlainMax)
		{
		mini = minV
		maxi = maxV	
		}
	else
		{
		//choose either min or max or +- 5 standard deviations if extreme outliers		
		mini = max(minV,minS)
		maxi = min(maxV,maxS)
		}
	}

/**************core plotting function**********************/

void MSA112_PlotAxis(image scatter,number zero,string Axis, number red, number green, number blue)
	{	//CURRENTLY NOT IN USE
	number width,height
	scatter.getSize(width,height)
	
	imagedisplay disp = ImageGetImageDisplay(scatter,0)
	number AxisID 
	if(Axis =="X") AxisID =CreateLineAnnotation(scatter,zero,0,zero,width)
	if(Axis =="Y") AxisID =CreateLineAnnotation(scatter,0,zero,height,zero)
	
	//set line properties
	component Axis =disp.ComponentGetChildbyID(AxisID)
	Axis.ComponentSetDrawingMode(2)	//no background
	Axis.ComponentSetForegroundColor(red,green,blue)	//color
	Axis.ComponentSetMovable(0)	//not movable
	Axis.ComponentSetResizable(0)	//not resizable
	}
	
	
image  MSA112_ScatterPlot2components(image &score1,image &score2,image cluster,number F,number S,number currCluster,number sparsity)
{
//curCluster = -1	plot all clusters differently colored		kPoints = Points
//curCluster = 0	plot all clusters black						kPoints = Points
//curCluster > 0	plot one specific cluster black		Tall must correspond kPoints = Tall height 

string MSAlabel = MSA112_Label()

number Points = cluster.ImageGetDimensionSize(1)

//evaluate scatter range
number Mi1,Ma1,Mi2,Ma2
MSA112_FindMinMaxPlot(score1,Mi1,Ma1)
MSA112_FindMinMaxPlot(score2,Mi2,Ma2)

//find scatter width  and height
number disp = sqrt((Ma1-Mi1)*(Ma2-Mi2)/Points/sparsity)	//sparsity defines scale
number width = trunc((Ma1-Mi1)/disp)+2
number height = trunc((Ma2-Mi2)/disp)+2

//rescale scores to scatter range (in pxl) 
score1 = (score1-Mi1)/disp
score2 = (score2-Mi2)/disp
score1 = round(score1)
score2 = round(score2)

//make scatter plot and calibrate it 
image Scatter := rgbimage("",4,width,height)
Scatter = rgb(256,256,256)
number zero0 = -Mi1/disp
number zero1 = -Mi2/disp
//result("\nzeroF "+zero0+" zeroS "+zero1)
scatter.ImagesetDimensionCalibration(0,zero0,disp,"",1)
scatter.ImagesetDimensionCalibration(1,zero1,disp,"",1)

//plot axis and axis labels
if(zero0>=0&&zero0<width) scatter[0,zero0,height,zero0+1] = rgba(0,0,0,0) 
if(zero1>=0&&zero1<height) scatter[zero1,0,zero1+1,width] = rgba(0,0,0,0)
Scatter.CreateTextAnnotation (zero1,0.9*width,F+"")
Scatter.CreateTextAnnotation (0.9*height,zero0,S+"")
	
//set pixels
if(currCluster==(-1))
	{
	//cluster-reduced scores
	image score1C = imageclone(score1)
	image score2C = imageclone(score2)
	
	number clusterN = max(cluster)
	for(number cCluster=0;cCluster<=clusterN;cCluster++)
		{
		//reduce to current cluster
		score1C = tert(cluster==cCluster,score1,0)
		score2C = tert(cluster==cCluster,score2,0)
		
		//plot pixels
		MSA112_PaintCluster(Scatter,score1C,score2C,cCluster)			//specific color
		}
	}
else
	{
	//plot pixels
	MSA112_PaintCluster(Scatter,score1,score2,0)	//all black
	}	
	
//cleanup borders (paint to white)
Scatter[0,0,1,width] = rgb(256,256,256)		//top border
Scatter[0,0,height,1] = rgb(256,256,256)		//left border
Scatter[height-1,0,height,width] = rgb(256,256,256)	//bottom border
Scatter[0,width-1,height,width] = rgb(256,256,256)	//right border
	
return Scatter
}


taggroup MSA112_getScatterTags(image scatter)
	{
	TagGroup ImTags = Scatter.ImageGetTagGroup()	
	TagGroup ScatterTags 
	try
		{
		//check if this is MSA-type image, otherwise error
		ImTags.TagGroupGetTagAsTagGroup(MSA112_label(),ScatterTags)
		
		//check if it is scatterplot
		string type =ScatterTags.MSA112_getStringFromTags("type")
		if(type!="MSA scatterPlot") throw("")
		}
	catch 
		{
		OKdialog("It must be an scatterPlot in front !!");
		exit(0);
		}
	
	return ScatterTags
	}
	

/**************store global ID of MSA matrix in scatterplot **********************/	
void  MSA112_StoreGlobalID(image in,image out,string path,string type)
	{
	object ID = in.ImageGetUniqueID()
	number ID1,ID2,ID3,ID4
	ID.GetValues(ID1,ID2,ID3,ID4)
	//result("\n"+ID1+" "+ID2+" "+ID3+" "+ID4)

	TagGroup Tags = out.ImageGetTagGroup()
	Tags.MSA112_setNumberInTags(Path+":"+type+" ID:1","L",ID1)
	Tags.MSA112_setNumberInTags(Path+":"+type+" ID:2","L",ID2)
	Tags.MSA112_setNumberInTags(Path+":"+type+" ID:3","L",ID3)
	Tags.MSA112_setNumberInTags(Path+":"+type+" ID:4","L",ID4)
	}


/**************find MSAimage from its global ID stored in scatterplot *************/
image  MSA112_OpenImageFromGlobalID(taggroup Tags,string path,string type)
	{
	number ID1,ID2,ID3,ID4
	ID1 =Tags.MSA112_getNumberFromTags(path+":"+type+" ID:1","L")
	ID2 =Tags.MSA112_getNumberFromTags(path+":"+type+" ID:2","L")
	ID3 =Tags.MSA112_getNumberFromTags(path+":"+type+" ID:3","L")
	ID4 =Tags.MSA112_getNumberFromTags(path+":"+type+" ID:4","L")
//result("\n"+ID1+" "+ID2+" "+ID3+" "+ID4)

	object ID =alloc(imageid)
	ID.init(ID1,ID2,ID3,ID4)
//result("\n"+ID.IsOpen())
	image MSAimage 
	if(ID.CanOpen()) 
		{
		MSAimage := ID.Open()
		return MSAimage
		}
	else throw(type+" is not opened")
	}
	
	
	
	
	
/************** head command to make scatter plot **********************/
image MSA112_PlotXY(TagGroup MatrixTags,number CurrCluster,string Method,number First,number Second,number sparsity)
	{
	string MSALabel = MSA112_Label()
	
	//read cluster
	image Cluster:= MSA112_GetCluster(MatrixTags)
	
	number width,height
	MSA112_GetSIWidthHeight(MatrixTags,width,height)

	//read and check the number of components
	string 	clusterLabel = MSA112_GetClusterLabel(currCluster)
	number CompN = MatrixTags.MSA112_ReadNumberComponents(currCluster,Method)
	if(CompN==0) {OKdialog("No PCA components in cluster "+CurrCluster);throw("");}
	if(!MSA112_CheckValidityOfComponents(First,Second,CompN,CurrCluster)) return NULL

	//read the number of points in cluster
	number kPoints = MatrixTags.MSA112_GetNumberPoints(CurrCluster)
//result("\ncluster "+CurrCluster+" ; "+m+ " components"+ " kPoints "+kPoints)	

	//read scores
	image Tall:= MatrixTags.MSA112_getScores(Method,currCluster)

	//extract two scores
	image score1 = Tall[0,First-1,kPoints,First]
	image score2 = Tall[0,Second-1,kPoints,Second]

	//make scatterplot

	rgbimage scatter := MSA112_ScatterPlot2components(score1,score2,cluster,First,Second,currCluster,sparsity)
													//scores are converted into scatter plot ccordinates
	//convert scores to integer
	image score1I := integerimage("",2,0,1,kPoints)
	image score2I := integerimage("",2,0,1,kPoints)
	score1I = score1
	score2I = score2

	//prepare mask images
	image cluster2D := MSA112_ConvertColumnTo2D(cluster,width,height)

	image CoordX,CoordY
	if(currCluster<=0)
		{
		CoordX := MSA112_ConvertColumnTo2D(score1I,width,height)
		CoordY := MSA112_ConvertColumnTo2D(score2I,width,height)
		}
//CoordY.showimage()
	number clusterN = max(cluster)
	
	//set tags
	TagGroup ScatterTags = Scatter.ImageGetTagGroup()
	ScatterTags.MSA112_setStringinTags(MSAlabel+":type","MSA scatterPlot")
	ScatterTags.MSA112_setNumberInTags(MSAlabel+":first","S",First) 
	ScatterTags.MSA112_setNumberInTags(MSAlabel+":second","S",Second)  
	ScatterTags.MSA112_setNumberInTags(MSAlabel+":clusters number","S",clusterN) 
	ScatterTags.MSA112_setNumberInTags(MSAlabel+":width","L",width) 
	ScatterTags.MSA112_setNumberInTags(MSAlabel+":height","L",height) 
	ScatterTags.MSA112_setImageinTags(MSAlabel+":mask",cluster2D)

	if(currCluster<=0)
		{
	ScatterTags.MSA112_setImageinTags(MSAlabel+":Coordinates-X",CoordX)
	ScatterTags.MSA112_setImageinTags(MSAlabel+":Coordinates-Y",CoordY)
		}

	return scatter
	}


string MSA112_TruncateMSAname(string Dataname,number ToCut)
	{
	return left(DataName,len(DataName)-ToCut)
	}
	
/********************  from MSAimage *******************************/
image MSA112_PlotXY(image MSAimage,number CurrCluster,string Method,number First,number Second,number sparsity)
	{
	TagGroup MatrixTags =  MSA112_GetMatrixTags(MSAimage,0)
	
	image scatter :=MSA112_PlotXY(MatrixTags,currCluster,Method,First,Second,sparsity)
 
	string DataName = MSAimage.getName().MSA112_TruncateMSAname(4)	//remove " MSA"
	scatter.setName(DataName+" SP "+First+"-"+Second+"("+currCluster+") "+Method)	

	string MSAlabel = MSA112_Label() 
	MSA112_StoreGlobalID(MSAimage,Scatter,MSAlabel,"MSA Matrix")
		
	return scatter
	}
	
/********************  from MSAimage , 0th-cluster only********************************/
image MSA112_PlotXY(image MSAimage,number First,number Second,string Method, number sparsity)
	{
	number currClust = 0
	
	return MSA112_PlotXY(MSAimage,currClust,Method,First,Second,sparsity)
	}

/********************  from MSAimage , PCA,0th-cluster only********************************/
image MSA112_PlotXY(image MSAimage,number First,number Second,number sparsity)
	{
	string Method = "PCA"
	number currClust = 0
	
	return MSA112_PlotXY(MSAimage,currClust,Method,First,Second,sparsity)
	}
	
	












/*********************					*
**************************				**
DISPLAY LOADS *****************************
**************************				**
*********************					*/

	
	
	
	
	
	
	
	
	
	


void  MSA112_DisplayLoads(image AverSpec,image Pall,number Comp0,number CompF,string Name,string Method,image Components)
{
number Edepth,compN
Pall.getSize(compN,Edepth)

if(Comp0<1) Throw("component must have index greater then zero !!")
if(CompF<Comp0) Throw("First component must have index not greater then that of the Last component !!")
if(CompF>compN) CompF = compN

image ClippedComponents
number defined
if(ImageIsValid(Components))if(Components.ImageGetDimensionSize(1)>CompF) defined = 1 
if(defined ==1) 
	{
	ClippedComponents = Components[0,0,CompF+1,1]
	}
else 
	{
	ClippedComponents :=integerimage("",1,0,1,CompF+1)
	ClippedComponents  = 1
	}
	
imagedisplay imdisp
AverSpec.setName(Name)
AverSpec.showimage()
imdisp = ImageGetImageDisplay(AverSpec,0)
ImageDisplaySetSliceLabelByID(imdisp,ImageDisplayGetSliceIDByIndex(imdisp,0),"average")

number originE,scaleE
string unitE
AverSpec.ImageGetDimensionCalibration(0,originE,scaleE,unitE,1)

for(number i=Comp0;i<=CompF;i++)
	{
	string currMethod = Method
	if(!ClippedComponents.getPixel(0,i))  currMethod = "PCA" //override unrotated component
	
	image spec = MatrixTranspose(Pall[0,i-1,Edepth,i])
	spec.ImageSetDimensionCalibration(0,originE,scaleE,unitE,1)
	
	imdisp.imageDisplayAddImage(spec,currMethod+" "+i)
	}
	
imdisp.LinePlotImageDisplaySetLegendShown(1)
}	


/*************************** simplified , no clipped components *********************/
void  MSA112_DisplayLoads(image AverSpec,image Pall,number Comp0,number CompF,string Name,string Method)
{
image Components

MSA112_DisplayLoads(AverSpec,Pall,Comp0,CompF,Name,Method,Components)
}	

/*************************** simplified , PCA only *********************/
void  MSA112_DisplayLoads(image AverSpec,image Pall,number Comp0,number CompF,string Name)
{
image Components
string Method= "PCA"

MSA112_DisplayLoads(AverSpec,Pall,Comp0,CompF,Name,Method,Components)
}
	
	
	
imagedisplay MSA112_PlotLine(image imag,number level,string Label)
	{
	imagedisplay imdisp = ImageGetImageDisplay(imag,0)
	imdisp.ImageDisplaySetSliceLabelByID(ImageDisplayGetSliceIDByIndex(imdisp,0),"aniso")
		
	image ImLevel :=imag.ImageClone()
	ImLevel =level
	
	imdisp.imageDisplayAddImage(ImLevel,Label)
	}	


imagedisplay MSA112_PlotCut(image imag,number comp,string Label)
	{
	imagedisplay imdisp = ImageGetImageDisplay(imag,0)
	
	ROI cut =NewRoi()
	cut.ROISetRange(0,comp)
	
	cut.ROISetVolatile(0)
	cut.ROISetColor(1,0,0)
	cut.ROISetLabel(Label)
	
	imdisp.ImageDisplayAddROI(cut)
	
	return imdisp
	}








/*********************					*
**************************				**
DISPLAY EIGENVALUES ***********************
**************************				**
*********************					*/

	
	
	
	
	
	
	
	
	
	
void MSA112_DisplayEigenValues(image Kurtosis,string Name,string Method,image Components)
	{
	//this is cast for Kurtosis but can be variance as well
	image KurtPCA = Kurtosis.imageclone()
	image KurtICA = Kurtosis.imageclone()
	KurtPCA=0
	KurtICA=0	
	
	number CompN = Kurtosis.ImageGetDimensionSize(0)
	if(!components.ImageIsValid()) components := integerImage("",1,0,1,CompN+1)
	
	for(number i=0;i<CompN;i++)
		{	
		if(Components.getPixel(0,i+1)) KurtICA.setPixel(i,0,Kurtosis.getPixel(i,0))
		else KurtPCA.setPixel(i,0,Kurtosis.getPixel(i,0))
		}
		
	KurtPCA.showimage()
	KurtPCA.setName(Name)		
	KurtPCA.ImagesetDimensionCalibration(0,-0.5,1,"components",1) //such as line plot: 1st comp ==1, 2nd==2,...

	imagedisplay imdisp = ImageGetImageDisplay(KurtPCA,0)
	imdisp.ImageDisplaySetSliceLabelByID(ImageDisplayGetSliceIDByIndex(imdisp,0),"PCA")	
	imdisp.LinePlotImageDisplaySetSliceDrawingStyle(0,1)		//draw as line
	imdisp.LinePlotImageDisplaySetSliceComponentColor(0,0,0,0,0)//set line color black
	imdisp.MSA112_SetLineThickness(0,2)
	
	if(sum(Components)>0)
		{
		//this is only for kurtosis of rotated loadings
		imdisp.imageDisplayAddImage(KurtICA,Method)
		imdisp.LinePlotImageDisplaySetSliceDrawingStyle(1,2)	//draw as filled
		imdisp.LinePlotImageDisplaySetSliceComponentColor(1,1,0,0,1)
		}
	}
	















/*********************					*
**************************				**
Make Container ****************************
**************************				**
**********************					*/














/******************************************************************
Function Calculate Sparsity
******************************************************************/
number MSA112_Sparsity(image Cube)
{
number width,height,depth
Cube.Get3dSize(width,height,depth)

image mask = cube.imageclone()
mask = tert(cube==0,0,1)

number Sparsity = sum(mask)/width/height/depth

return Sparsity
}





/******************************************************************
Function MakeMatrix
******************************************************************/

image MSA112_SelectEnergyRange(image dataCube,number &Estart,number &Eend)
	{
	//THIS FUNCTION WORKS IN THREAD (BACKGROUND ONLY !!!
	number Depth = DataCube.ImageGetDimensionSize(2)
	
	Estart = 0
	Eend = depth/2
	
	//select energy range
	image Overview := MSA112_MakeAverageSpectrum(DataCube)
	Overview.showimage()
	imageDisplay disp = Overview.ImageGetImageDisplay(0)
	roi Range  = NewRoi()
	Range.ROISetRange(0,Depth/2) 
	disp.ImageDisplayAddRoi(Range)
	
	if(MSA112_SetFeatureInImage("Select the energy range by moving ROI.\nClick Cancel if the whole range is required.","Treat selected"))
		{
		Range.ROIGetRange(Estart,Eend)
		}
	else Eend = Depth
	
	disp.ImageDisplayDeleteROI(Range)

	image AvSpec = OverView[0,Estart,1,Eend]
	
	OverView.deleteimage()

	//return average spectrum	
	return AvSpec
	}
	

image MSA112_MakeMSAimage(image dataCube,image &AvImage, number Estart,number Eend)
	{												//AvSpec and AverImage returned because they are costly
	//extact all parameters from dataCube
	string DataName = DataCube.getName()
	number width,height,Depth 
	DataCube.get3Dsize(width,height,Depth )
	number Points = width*height
	number ScaleX,ScaleY,ScaleE,ScaleI
	number OriginX,OriginY,OriginE,OriginI
	string unitX,unitY,unitE,unitI
	DataCube.ImageGetDimensionCalibration(0,OriginX,ScaleX,unitX,1)
	DataCube.ImageGetDimensionCalibration(1,OriginY,ScaleY,unitY,1)
	DataCube.ImageGetDimensionCalibration(2,OriginE,ScaleE,unitE,1)
	ScaleI = DataCube.ImageGetIntensityScale()
	OriginI = DataCube.ImageGetIntensityOrigin()
	unitI = DataCube.ImageGetIntensityUnitString()
//result("\nSI data "+width+"x"+height+"="+N+"pixels, energy depth "+Depth+ " scale "+ScaleX+" dispersion "+ScaleE+unitE)	
	
	//make MSA image
	image MeanImage := MSA112_MakeMeanImage(DataCube)	
	image MSAimage := MeanImage.ImageClone()
	MSAimage.setName(DataName+" MSA")
	
	//make average image
	AvImage = MeanImage.MSA112_Convert2DtoColumn()

	//find pivot channel or quadrants
	number pivot = depth/2 - Estart
	if(pivot<0) pivot =0			//pivot is left of the selected range
	if(pivot>(Eend-Estart)) pivot ==0 //pivot is right of the selected range
	
	//correct depth and originE according selected energy range
	Depth = Eend - Estart 
	OriginE -= Estart
	
	//make matrix tags
	TagGroup IniTags = DataCube.ImageGetTagGroup()
	TagGroup MatrixTags = MSAimage.ImageGetTagGroup()

	string MSALAbel = MSA112_Label()
	MatrixTags.MSA112_setStringInTags(MSAlabel+":type","MSA matrix")
	MatrixTags.MSA112_setNumberInTags(MSAlabel+":width","L",width)
	MatrixTags.MSA112_setNumberInTags(MSAlabel+":height","L",height)
	MatrixTags.MSA112_setNumberInTags(MSAlabel+":pixels","L",Points)
	MatrixTags.MSA112_setNumberInTags(MSAlabel+":energy depth","L",Depth)	
	MatrixTags.MSA112_setNumberInTags(MSAlabel+":energy pivot","L",pivot)	
	MatrixTags.MSA112_setStringInTags(MSAlabel+":energy units",unitE)
	MatrixTags.MSA112_setNumberInTags(MSAlabel+":energy dispersion","F",scaleE)
	MatrixTags.MSA112_setNumberInTags(MSAlabel+":energy origin","L",OriginE)
	MatrixTags.MSA112_setNumberInTags(MSAlabel+":scale X","F",ScaleX)
	MatrixTags.MSA112_setNumberInTags(MSAlabel+":origin X","L",originX)
	MatrixTags.MSA112_setStringInTags(MSAlabel+":units X",unitX)
	MatrixTags.MSA112_setNumberInTags(MSAlabel+":scale Y","F",ScaleY)
	MatrixTags.MSA112_setNumberInTags(MSAlabel+":origin Y","L",originY)
	MatrixTags.MSA112_setStringInTags(MSAlabel+":units Y",unitY)
	MatrixTags.MSA112_setNumberInTags(MSAlabel+":scale Int","F",ScaleI)
	MatrixTags.MSA112_setNumberInTags(MSAlabel+":origin Int","L",originI)
	MatrixTags.MSA112_setStringInTags(MSAlabel+":units Int",unitI)
	
	TagGroupCopyTagsFrom(MatrixTags,IniTags)
	
	return MSAimage	
	}















/*********************					*
**************************				**
PCA  decomposition  ***********************
**************************				**
*********************					*/

	
	
	
	
	
	
	
	
	
	
	

	
/******************************************************************
ensure cluster Taglist
******************************************************************/

TagGroup  MSA112_EnsureClusterList(tagGroup MatrixTags)
{
string MSAlabel = MSA112_Label()

//set tag about clusters	
TagGroup ClusterList = MatrixTags.TagGroupGetOrCreateTagList(MSAlabel+":clusters")
number Points
if(!MatrixTags.TagGroupGetTagAsLong(MSAlabel+":clusters:[0]"+":pixels",Points))
	{	
	//[0]-cluster must exist consisting number of points Points 
	ClusterList.TagGroupCreateGroupTagAtEnd()
	MatrixTags.TagGroupSetTagAsLong(MSAlabel+":clusters:[0]"+":pixels",Points)
	}
return ClusterList	
}

/******************************************************************
functions prepare PCA standard images
******************************************************************/
image MSA112_PCAimageEigens(number maxIt,number deviaSpec,number converg)
{
image eigens := realimage("",4,1,maxIt+1)
if(converg!=0) eigens.setPixel(0,0,deviaSpec/100) //0-th pixel for deviaSpec (interpreted in %) 

return eigens
}

image MSA112_PCAimageDevia(number maxIt,number mincos)
{
image devia := realimage("",4,1,maxIt+1)
devia.setPixel(0,0,minCos)		//0-th pixel for specified minimal cos of angle

return devia
}

image MSA112_PCAimageSpikes(number maxIt,number spikesCode)
{
image spikes := integerimage("",2,0,1,maxIt+1)
spikes.setPixel(0,0,SpikesCode)	//0-th pixel keeps a code for all spikes parameters

return spikes
}

image MSA112_PCAimagesAll(image PCAImageOne,number Components)
{	//makes many-column image fro one-column one
number type = PCAImageOne.ImageGetDataType()
number Height = PCAImageOne.ImageGetDimensionSize(1)

return MSA112_Make2DImageOfRightType(type,"",Components,Height)
}



/******************************************************************
Function Resort Descending
******************************************************************/
image MSA112_ResortDescending(image input, image Score)
{ //Score: {1,3,2,5,4} -> {last,next,then,largest,high}
image ScoreCopy = Score.ImageClone()
image output = input.imageClone()

number compN,height
input.getSize(compN,height)
for(number i = 0;i<compN;i++)
	{
	number largest,dummy
	ScoreCopy.max(largest,dummy)
	output[0,i,height,i+1] = input[0,largest,height,largest+1]
	ScoreCopy.setPixel(largest,0,0)
	}
	
return output
}


/******************************************************************
Function Reduce matrix on one component
******************************************************************/
image MSA112_ReduceDataMatrix(image dataMatrix,image Pvec,image Tvec)
{
return  (DataMatrix - matrixMultiply(Tvec,MatrixTranspose(Pvec)) )
}


/******************************************************************
Functions Store PCA tags
******************************************************************/
//basic parameters
void MSA112_StorePCAtags(TagGroup MatrixTags,number currCluster,number Points,number maxItt,number deviaSpec)
	{
	string ClusterLabel = MSA112_GetClusterLabel(currCluster)
	
	MatrixTags.MSA112_setNumberInTags(ClusterLabel+":PCA:number components","S",Points)
	MatrixTags.MSA112_setNumberInTags(ClusterLabel+":PCA:max iterations","I",maxItt)
	MatrixTags.MSA112_setNumberInTags(ClusterLabel+":PCA:required precision","F",deviaSpec)
	}

//spikes info	
	void MSA112_StorePCAtags(TagGroup MatrixTags,number currCluster,number OutTreat,number OutCorr,number OutThres, image spikesXY)
	{
	string ClusterLabel = MSA112_GetClusterLabel(currCluster)

	MatrixTags.MSA112_setNumberInTags(ClusterLabel+":PCA:Outlier find","S",OutTreat)
	MatrixTags.MSA112_setNumberInTags(ClusterLabel+":PCA:Outlier remove","S",OutCorr)
	MatrixTags.MSA112_setNumberInTags(ClusterLabel+":PCA:Outlier threshold","F",OutThres)
	
	MatrixTags.TagGroupDeleteTagWithLabel(ClusterLabel+":PCA:Outliers coord")
	if(spikesXY.ImageIsValid())
		{
		number spikesN = spikesXY.ImageGetDimensionSize(1)
	MatrixTags.MSA112_setNumberInTags(ClusterLabel+":PCA:outliers number","L",SpikesN)
	MatrixTags.MSA112_setImageinTags(ClusterLabel+":PCA:Outliers coord",spikesXY)
		}	
	}

	void MSA112_StorePCAtags(TagGroup MatrixTags,string extra,number currCluster,image Pall,image Tall,image finEVs)
	{
	string ClusterLabel = MSA112_GetClusterLabel(currCluster)

	//to ensure that tags from previous versions notinterfere
	MatrixTags.TagGroupDeleteTagWithLabel(ClusterLabel+":PCA:"+extra+"Loadings")
	MatrixTags.TagGroupDeleteTagWithLabel(ClusterLabel+":PCA:"+extra+"Scores")
	MatrixTags.TagGroupDeleteTagWithLabel(ClusterLabel+":PCA:"+extra+"Variance")
	
	MatrixTags.MSA112_setImageinTags(ClusterLabel+":PCA:"+extra+"Loadings",Pall)
	MatrixTags.MSA112_setImageinTags(ClusterLabel+":PCA:"+extra+"Scores",Tall)
	MatrixTags.MSA112_setImageinTags(ClusterLabel+":PCA:"+extra+"Variance",finEVs)	
	}

	void MSA112_StorePCAtags(TagGroup MatrixTags,number currCluster,image Pall,image Tall,image finEVs)
	{
	MSA112_StorePCAtags(MatrixTags,"",currCluster,Pall,Tall,finEVs)	
	}
	
	void MSA112_StorePCAtags(TagGroup MatrixTags,number currCluster,image EVsAll,image DeviaAll,image spikesAll,image Itterations)
	{
	string ClusterLabel = MSA112_GetClusterLabel(currCluster)

	//to ensure that tags from previous versions notinterfere
	MatrixTags.TagGroupDeleteTagWithLabel(ClusterLabel+":PCA:Valuesflow")
	MatrixTags.TagGroupDeleteTagWithLabel(ClusterLabel+":PCA:Precisions")
	MatrixTags.TagGroupDeleteTagWithLabel(ClusterLabel+":PCA:Outliers")
	MatrixTags.TagGroupDeleteTagWithLabel(ClusterLabel+":PCA:Iterations")
	
	MatrixTags.MSA112_setImageInTags(ClusterLabel+":PCA:Valuesflow",EVsAll)
	MatrixTags.MSA112_setImageInTags(ClusterLabel+":PCA:Precisions",DeviaAll)
	MatrixTags.MSA112_setImageInTags(ClusterLabel+":PCA:Outliers",spikesAll)
	MatrixTags.MSA112_setImageInTags(ClusterLabel+":PCA:Iterations",Itterations)
	}

	void MSA112_StorePCAtags(TagGroup MatrixTags,number currCluster,number uncentered,number points,image AvSpec,number quadrant)
	{
	string ClusterLabel = MSA112_GetClusterLabel(currCluster)
		
	MatrixTags.MSA112_setNumberInTags(ClusterLabel+":pixels","L",points)
	MatrixTags.TagGroupSetTagAsArray(ClusterLabel+":Average Spectrum",AvSpec)
	MatrixTags.MSA112_setNumberInTags(ClusterLabel+":uncentered","S",uncentered)
	MatrixTags.MSA112_setNumberInTags(ClusterLabel+":quadrant","S",quadrant)	
	}
	


















/*********************					*
**************************				**
Reconstruction ****************************
**************************				**
*********************					*/

	
	
	
	
	
	
	
	
	
	
	
	


/******************************************************************
Function: Project clustered score to all pixels
******************************************************************/
image MSA112_ProjectClusteredScore(TagGroup MatrixTags,image Tall,number currCluster)
	{			//extends cluster-specific score (Tall) to full-size score (RTall)
	string MSALabel = MSA112_Label()
	
	number Points = MSA112_GetNumberPoints(MatrixTags)

	number compN = Tall.ImageGetDimensionSize(0)
	
	image Cluster:= MSA112_GetCluster(MatrixTags) 
	
	image RTall:= realimage("T-vectors",4,compN,Points)	//full Points not kPoints
	number ii
	for(number i=0;i<Points;i++)	//pixel-by-pixel loop
		{
		if(cluster.getPixel(0,i)==currCluster) 
			{
			try{RTall[i,0,i+1,compN] = Tall [ii,0,ii+1,compN];}
			catch{break;}
			ii++
			}
		}
	
	return RTall				
	}


/******************************************************************
Function: Set out-cluster pixel to standard value
******************************************************************/
image MSA112_SetOutClusterPixels(TagGroup MatrixTags,image RTall,number currCLuster)
	{
	string MSALabel = MSA112_Label()
	
	number Points = MSA112_GetNumberPoints(MatrixTags)
	
	number compN = RTall.ImageGetDimensionSize(0)
	
	image Cluster:= MSA112_GetCluster(MatrixTags)
		
	number OutClusterValue = MSA112_OutClusterValue()
	
	if(currCluster!=(-1))	//set pixels out of given cluster
		{
		for(number i=0;i<Points;i++)	//pixel-by-pixel loop
			{
			if(cluster.getPixel(0,i)!=currCluster) 
				{
				RTall[i,0,i+1,compN] = OutClusterValue
				}
			}
		}
	else					//set pixels out of any cluster
		{
		for(number i=0;i<Points;i++)
			{
			if(cluster.getPixel(0,i)==0) 
				{
					RTall[i,0,i+1,compN] = OutClusterValue
				}
			}
		}	
	
	return RTall				
	}
	
	

	

	
/******************************************************************
Function: roughly correct spectrum weighting for previous normalization on Average image
******************************************************************/
 image MSA112_CorrectSpectrumForAverageImage(image spectrum, TagGroup MatrixTags,number currCluster)
	{
	//data in cluster was divided on sqrt(AvImage)	
	//all spectra from that data were divided too
	//to restore the spectrum : multiply on mean(AvImage) (this is however rough)
	
	string MSALabel = MSA112_Label()
	
	number WeightSpecOnly =MatrixTags.MSA112_getNumberFromTags(MSAlabel+":weighted spectra only","S")
	
	if(!WeightSpecOnly)		
		{					
		number ImScaling				
		image AvIm = MSA112_GetAverageImage(MatrixTags)
			
		if(currCluster == 0)
			{
			//mean of AverImage
			ImScaling = mean(AvIm)
			}
		else
			{
			//mean of AverImage on currCluster only
			image cluster = MSA112_GetCluster(MatrixTags)		
			cluster = tert(cluster==currCluster, 1,0)
			AvIm *= cluster
					
			number kPoints = MSA112_GetNumberPoints(MatrixTags,currCluster)
			ImScaling = sum(AvIm)/kPoints				
			}
			
		spectrum *= sqrt(ImScaling)	
		
		return 	spectrum
		}
		
	else return spectrum	//in this case AvImage = 1	
	}		
	
	
	
/******************************************************************
Function: Restore apectrum after weighting 
******************************************************************/
image MSA112_RestoreWeightedSpectrum(image Spectrum,TagGroup MatrixTags,number currCluster)
	{	//(cluster dependent because account avImage in cluster)
	string MSALabel = MSA112_Label()
	
	image OutSpectrum := imageclone(Spectrum)
	
	number weighting =MatrixTags.MSA112_getNumberFromTags(MSAlabel+":weighted","S")

	if(Weighting)
		{
		image GlobalAvSpectrum = MSA112_GetAverageSpectrum(MatrixTags)
		
		GlobalAvSpectrum = sqrt(GlobalAvSpectrum)		
		GlobalAvSpectrum  = MSA112_CorrectSpectrumForAverageImage(GlobalAvSpectrum,MatrixTags,currCluster)
		
		OutSpectrum = Spectrum * GlobalAvSpectrum[icol,0] //spectrum may be multiple in Y
		}	
	
	return OutSpectrum
	}



/******************************************************************
Function re-Weight apectrum 
******************************************************************/
image MSA112_ReWeighSpectrum(image Spectrum,TagGroup MatrixTags,number currCluster)
	{
	string MSALabel = MSA112_Label()
	
	image OutSpectrum = imageclone(Spectrum)
	
	number weighting =MatrixTags.MSA112_getNumberFromTags(MSAlabel+":weighted","S")

	if(Weighting)
		{
		image GlobalAvSpectrum = MSA112_GetAverageSpectrum(MatrixTags)
		
		GlobalAvSpectrum = sqrt(GlobalAvSpectrum)
		GlobalAvSpectrum  = MSA112_CorrectSpectrumForAverageImage(GlobalAvSpectrum,MatrixTags,currCluster)
				
		OutSpectrum = Spectrum / GlobalAvSpectrum[icol,0] 
		//spectrum may be multiple in Y
		}	
	
	return OutSpectrum
	}
	

	




void MSA112_CalibrateRecCube(tagGroup MatrixTags,image &recCube)
	{
	//read calibration
	number disp,originE
	number ScaleX,ScaleY,ScaleI,originX,originY,originI
	string unitE,unitX,unitY,unitI
	MatrixTags.MSA112_GetCalibration(disp,originE,unitE,"E")
	MatrixTags.MSA112_GetCalibration(ScaleX,originX,unitX,"X")
	MatrixTags.MSA112_GetCalibration(ScaleY,originY,unitY,"Y")
	MatrixTags.MSA112_GetCalibration(ScaleI,originI,unitI,"Int")

	//calibrate and copy previous info
	if(ScaleX==0) ScaleX = 1
	if(ScaleY==0) ScaleY = 1
	if(disp==0) disp = 1
	RecCube.ImageSetDimensionCalibration(0,OriginX,ScaleX,unitX,1)
	RecCube.ImageSetDimensionCalibration(1,OriginY,ScaleY,unitY,1)
	RecCube.ImageSetDimensionCalibration(2,OriginE,disp,unitE,1)
	RecCube.ImageSetIntensityScale(ScaleI)
	RecCube.ImageSetIntensityOrigin(OriginI)
	RecCube.ImageSetIntensityUnitString(unitI)
		
	}
	
void MSA112_Weighting(image &dataMatrix,image AvSpec, image AvImage)
	{
	//prepare square root
	image AvSpecSqrt = imageclone(AvSpec)
	image AvImageSqrt = imageclone(AvImage)
	
	AvSpecSqrt = tert(AvSpec>0,sqrt(AvSpec),sqrt(-AvSpec))
	AvImageSqrt = tert(AvImage>0,sqrt(AvImage),sqrt(-AvImage))
	
	try{
		//prepare weight for spectrum
		image ToWeight = MatrixMultiply(AvImageSqrt,AvSpecSqrt)
			
		//weighting
		dataMatrix = tert(ToWeight!=0,dataMatrix/ToWeight,0)
		}
	catch{
		OKdialog("Not enough memory !!")
		throw("")
		}
	}

void MSA112_DeWeighting(image &dataMatrix,image AvSpec, image AvImage)
	{
	//prepare square root
	image AvSpecSqrt = imageclone(AvSpec)
	image AvImageSqrt = imageclone(AvImage)
	
	AvSpecSqrt = tert(AvSpec>0,sqrt(AvSpec),sqrt(-AvSpec))
	AvImageSqrt = tert(AvImage>0,sqrt(AvImage),sqrt(-AvImage))
	
	//prepare weight for spectrum
	image ToWeight = MatrixMultiply(AvImageSqrt,AvSpecSqrt)
		
	//deweighting
	dataMatrix = dataMatrix*ToWeight
	}
	
void MSA112_AttachImages(image &MSAimage,image dataMatrix,image AvSpec,image AvImage,number weighting,number WeightSpecOnly)
	{
	TagGroup MatrixTags = MSAimage.ImageGetTagGroup()
	
	string MSALAbel = MSA112_Label()
	
	//these arrays are set without extar info (width,height,type)
	MatrixTags.TagGroupSetTagAsArray(MSAlabel+":Average Spectrum",AvSpec)
	MatrixTags.TagGroupSetTagAsArray(MSAlabel+":Average Image",AvImage)	
	MatrixTags.TagGroupSetTagAsArray(MSAlabel+":Matrix",dataMatrix)
	
	MatrixTags.MSA112_SetNumberInTags(MSAlabel+":weighted","S",weighting)
	MatrixTags.MSA112_SetNumberInTags(MSAlabel+":weighted spectra only","S",WeightSpecOnly)
	}


image MSA112_CalculateRecMatrix(TagGroup MatrixTags,string MSAlabel)
	{
	//max available components are used
	image Tall:= MSA112_GetScores(MatrixTags,"PCA",0) 
	image Pall:= MSA112_GetLoadings(MatrixTags,"PCA",0) 

	//check if not too many components
	number maxComp =50
	number comp =Pall.ImageGetDimensionSize(0)
	if(comp >maxComp) if(!OKcanceldialog("There are "+comp+" principal components extracted\ndo you want to reconstruct with such a huge number ?")) throw("")
		
	image AvSpec := MSA112_GetAverageSpectrum(MatrixTags,0,"PCA") 

	image RecMatrix = matrixMultiply(Tall,MatrixTranspose(Pall))
	
	RecMatrix += MSA112_ExpandImageVertically(AvSpec,RecMatrix) 

	return RecMatrix
	}

image MSA112_CalculateResidialMatrix(TagGroup MatrixTags,string MSAlabel,image RecMatrix)
	{	
	image IniMatrix := MSA112_GetMSAMatrix(MatrixTags) 
//IniMatrix.showimage()
	
	image ResidialMatrix = IniMatrix - RecMatrix 
		
	return ResidialMatrix
	}


//from variance residual matrix
number MSA112_CalculateNoiseVar_Malinowski(image VarSpec, number PrimaryComponents)
	{
	number MeanErrorVar = mean(VarSpec)

	number Edepth = VarSpec.ImageGetDimensionSize(0)

	//"real error" estimated from secondary eigenvalues, but renormalizes for whole matrix
	return	MeanErrorVar*Edepth/(Edepth-PrimaryComponents)
	}


//from MatrixTags
number MSA112_CalculateNoiseVarianceMalinowski(taggroup MatrixTags)	
	{	
	image RecMatrix := MSA112_CalculateRecMatrix(MatrixTags,MSA112_label())

	image ResidialMatrix = MSA112_CalculateResidialMatrix(MatrixTags,MSA112_label(),RecMatrix)

	ResidialMatrix = ResidialMatrix**2
	
	//noise spectrum
	image VarSpec = MSA112_Render2DtoBottom(ResidialMatrix)	
	number Points = RecMatrix.imagegetDimensionSize(1)
	VarSpec /=Points
		
	
	string 	clusterLabel = MSA112_GetClusterLabel(0)
	number CompN =MatrixTags.MSA112_getNumberFromTags(ClusterLabel+":"+"PCA"+":number components","S")
	
	return MSA112_CalculateNoiseVar_Malinowski(VarSpec,CompN)
	}










/*********************					*
**************************				**
functions for EELS  ***********************
**************************				**
*********************					*/

	
	
	
	
	
	
	
	
	
	
		



/*************** get channel from eV ****************/
number MSA112_eVtoCh(number eV,number origin,number disp)
	{
	return origin + eV/disp 
	}
	
	
/*************** Linear LS interpolation 3D ****************/
void MSA112_LinearLSInterpolation3D(image &A, image &B,image X,image Y)
	{
	number H,W,D,H2,W2,D2
	X.get3DSize(W,H,D)
	Y.get3DSize(W2,H2,D2)
	if((W!=W2)||(H!=H2)||(D!=D2)) Throw("different X,Y arrays in LS fitting !!")

	image X2 = X**2
	image XY=X*Y
	 
	image sumX = MSA112_Render3DtoFront(X)
	image sumY = MSA112_Render3DtoFront(Y)
	image sumX2 = MSA112_Render3DtoFront(X2)
	image sumXY = MSA112_Render3DtoFront(XY)

	image Znam = D*sumX2-sumX**2
	A = (sumY*sumX2-sumX*sumXY)/Znam
	B = (D*sumXY-sumX*SumY)/Znam
	}


/******* Filter out singularities************************/
image MSA112_FilterSingularities(image in)
{
return tert(in>-infinity(),in,0)
}


/******* EELS background *****************************/	
image MSA112_FindEELSBackground(image in,number startEv,number endEv)
	{
	number origin,disp
	string units
	ImageGetDimensionCalibration(in,2,origin,disp,units,1)
	number start = MSA112_eVtoCh(startEv,origin,disp)
	number end = MSA112_eVtoCh(endEv,origin,disp)

	number H,W,D
	in.get3DSize(W,H,D)
	image Y = imageclone(in[0,0,start,W,H,end])
	Y = log(Y)
	image X = imageclone(Y)
	X = log((iplane-origin+Start)*Disp)

	image A,R
	MSA112_LinearLSInterpolation3D(A,R ,X,Y)
	//R.showimage()
	//A.showimage()
	image bckgrd = imageclone(in)
	bckgrd = A[icol,irow,0]+R[icol,irow,0]*log((iplane-origin)*disp)
	bckgrd = exp(bckgrd)
	bckgrd[0,0,0,W,H,start] = 0
	ImageSetDimensionCalibration(bckgrd,2,origin,disp,units,1)

	return bckgrd 
	}


/**************  extract EELS signal ******************************/
image ExtractEELSsignal112(image Spectrum,number PreStartEv,number PreEndEv,number PostStartEv,number PostEndEv)
	{
	number H,W,D,H2,W2,D2
	Spectrum.get3DSize(W,H,D)
	number origin,originX,originY,disp,scaleX,scaleY
	string units,unitsX,unitsY
	ImageGetDimensionCalibration(Spectrum,0,originX,scaleX,unitsX,1)
	ImageGetDimensionCalibration(Spectrum,1,originY,scaleY,unitsY,1)
	ImageGetDimensionCalibration(Spectrum,2,origin,disp,units,1)
	
	number PostStart = MSA112_eVtoCh(PostStartEv,origin,disp)
	number PostEnd   = MSA112_eVtoCh(PostEndEv,origin,disp)

	image bckgrd:= MSA112_FindEELSBackground(Spectrum,PreStartEv,PreEndEv)
	image signal3D = Spectrum[0,0,PostStart,W,H,PostEnd] - bckgrd[0,0,PostStart,W,H,PostEnd]

	//bckgrd.showimage()
	image Signal = MSA112_Render3DtoFront(signal3D)
	Signal = MSA112_FilterSingularities(Signal)
	
	ImageSetDimensionCalibration(Signal,0,originX,scaleX,unitsX,1)
	ImageSetDimensionCalibration(Signal,1,originY,scaleY,unitsY,1)
	
	return Signal
	}


/************ Find average count in SI ***********************/
number MSA112_SIaverageCount(image in)
	{
	string Name = in.getName()
	number w,h,d
	in.get3Dsize(w,h,d)
	image summed = in.MSA112_Render3DtoFront()
	//result("\n"+summed.mean()+" "+w+" "+h+" "+d)
	number AvCount = summed.mean()/d
result("\n"+name+" count per channel "+AvCount)

	return AvCount
	}












/*********************					*
**************************				**
ROTATION **********************************
**************************				**
*********************					*/













/******************************************************************
Function Gram-Schmidt Orthogonalization
******************************************************************/
/*************** one couple of axis*********************/
image MSA112_GramSchmidtStep(image CurrVec,image OrtSet)
	{
	//OrtSet(k,n) 	set of k-vectors already orthogonalized
	//CurrVec 		(k+1)the vector to be orthogonalized to that set 
	
	number n,k
	OrtSet.getSize(k,n)

	image CorrWeights  = MatrixMultiply(MatrixTranspose(CurrVec),OrtSet)
	
	image Proj2D = OrtSet.imageclone()
	Proj2D = CorrWeights[icol,0]
	Proj2D *= OrtSet
	image Proj = Proj2D.MSA112_Render2DtoLeft()

	CurrVec -= Proj
	
	return CurrVec
	}

image MSA112_GramSchmidt(image Set)
	{
	number n,m
	Set.getSize(m,n)
	image OrtSet = Set.imageclone()
	
	OrtSet[0,0,n,1] = Set[0,0,n,1].MSA112_NormalizeVector()
	
	for(number k=1;k<m;k++)
		{
		image CurrVec = Set[0,k,n,k+1]		
		OrtSet[0,k,n,k+1] = MSA112_GramSchmidtStep(CurrVec,OrtSet[0,0,n,k]).MSA112_NormalizeVector()
		}
		
	return OrtSet	
	}	
	

/******************************************************************
Function Calculate one fast ICA vector
******************************************************************/		
number MSA112_FindICAcomponent(image dataMatrix,number p,image &Wall,image Kurt,image &devia,number COntrast,number deviaSpec, number maxItt)	
	{
	number N,m
	GetSize(dataMatrix,N,m)
		
	image Wvec := realimage("W-vector",4,1,m)
	image WvecOld := realimage("W-vector",4,1,m)		
	//WvecOld.setPixel(0,p,1)	// =PCA component
	WvecOld = MSA112_randomvector(m)

	image wTx  :=realimage("",4,N,1)	
	image wTx2 :=realimage("",4,N,1)
	image wTxE :=realimage("",4,N,1)
	image Der1 :=realimage("",4,N,1)
	image Der2 :=realimage("",4,N,1)
	image Der1_2D :=realimage("",4,N,m)
	number MeanDev2

	number dev=3
	number diverge = 0
	number Nitt=0
	while(dev>deviaSpec)
		{
		Nitt++
		wTx = MatrixMultiply(MatrixTranspose(WvecOld),dataMatrix)
		wTx2 = wTx**2
		//calculate previous Kurtosis 
		Kurt.setPixel(0,Nitt-1,MSA112_Kurtosis(wTx2))	

		//new itteration

		if(Contrast==1)
			{	
			Der1 = wTx2*wTx/3
			MeanDev2 = 1
			}		
		if(Contrast==2)
			{
			Der1 = tanh(wTx)
			Der2 = 1-Der1**2
			MeanDev2 = mean(der2)
			}		
		if(Contrast==3)
			{	
			wTxE  =exp(-wTx2/2)
			Der1 = wTx*wTxE 
			Der2 = (1-wTx2)*wTxE
			MeanDev2 = mean(der2)
			}
			
		Der1_2D = Der1[icol,0]	

		Wvec =	MeanDev2*WvecOld - MSA112_Render2DtoLeft(dataMatrix*Der1_2D)/N
		MSA112_GramSchmidtStep(Wvec,Wall).MSA112_NormalizeVector()			

		number Ccos = sum(MatrixMultiply(MatrixTranspose(Wvec),WvecOld))
		if(Ccos<0) Wvec = -Wvec
		dev = Acos(modulus(Ccos))
		devia.setPixel(0,Nitt,dev)	
		
		WvecOld = Wvec
		if(Nitt>=MaxItt)
			{
			diverge =1
			break
			}
		}	
		
	//final Kurtosis
	wTx = MatrixMultiply(MatrixTranspose(Wvec),dataMatrix)
	wTx2 =  wTx**2
	Kurt.setPixel(0,Nitt,MSA112_Kurtosis(wTx2))
		
	//final Wvec			
	if(diverge==1)
		{	
		Wvec=0	//trash component
		}
	Wall[0,p-1,m,p] = Wvec

	return Nitt
	}


/******************************************************************
Function rotate Decart
******************************************************************/	
void MSA112_rotate(number fi,image &score1,image &score2)
	{
	number ccos = cos(fi)
	number ssin = sin(fi)
//result("\n"+ccos+" "+ssin)
	
	image score11 = ccos*score1+ssin*score2
	score2 = -ssin*score1+ccos*score2
	score1 = score11
	}


/******************************************************************
Function Calculate 4th moment
******************************************************************/	
number MSA112_VRMCriteria(image Score)
	{
	number N = score.imagegetDimensionSize(1)
	
	image ScoreQ = Score**2		//^2	
	ScoreQ = ScoreQ**2			//^4
	image SQQ = ScoreQ.MSA112_Render2DtoBottom()
	SQQ /= N					//sum(^4)*N
	
	return sum(SQQ)
	}

	
image MSA112_VRMCriteria(image Score1,image Score2,number M,number sign)
	{
	number N = score1.imagegetDimensionSize(1)
	
	image ScoreA:= realimage("criteria",4,M,N)

	ScoreA = score1[0,irow]*cos(icol*90/M/57.2958) +sign*score2[0,irow]*sin(icol*90/M/57.2958)

	ScoreA = ScoreA**2		//^2
	ScoreA = ScoreA**2		//^4
	image SQQ = ScoreA.MSA112_Render2DtoBottom()
	SQQ /= N				//sum(^4)*N
	
	return SQQ
	}
	
	
/******************************************************************
Function Calculate VRM angle
******************************************************************/	
number MSA112_VRMAngle(image score1,image score2,number res )
	{
	number   M = 3.14159/2/res //range 0...90 grad
//result("\nM "+M)
	number N = score1.imagegetDimensionSIze(1)
	
	image norm = (score1**2+score2**2).sqrt()
	image scoreN1 = score1/norm
	image scoreN2 = score2/norm
	
	image Criteria =  MSA112_VRMCriteria(scoreN1,scoreN2,M,1)+MSA112_VRMCriteria(scoreN2,scoreN1,M,-1)

	//criteria.showimage()
	
	number Xmax,dummy
	Criteria.max(Xmax,dummy)

	return Xmax * res
	}

	
/******************************************************************
Function number of pair
******************************************************************/
number MSA112_combination(number k)
{
if(k>2) return (k-1+MSA112_combination(k-1))
else return 1
}




/******************************************************************
Class Loading Viewer
******************************************************************/
class MSA_LoadingViewer : object
{
number ListArr
number ScatterID,LoadID
number AxisID

image load1,load2

object init(object self, image scatter,number AxisIDP,image Load1P,image Load2P)
	{
	result("\nLoadingViewer initialized") 
	load1=load1P
	load2=load2P
	ScatterID = scatter.ImageGetID()
	AxisID = AxisIDP
	
	imagedisplay disp = scatter.ImageGetImageDisplay(0)

	image Load := imageclone(load1P)
	Load.SetName("current loading")
	LoadID = Load.ImageGetID()
	load.showimage()

	ListArr = disp.ImageDisplayAddEventListener(self,"component_changed,component_removed:ArrowChanged")
	
	return self
	}
		
void ArrowChanged(Object self, Number e_fl, ImageDisplay disp, Number c_fl, Number c_fl2, Component c )
	{
	if(e_fl==33554432) 
		{
		disp.ImageDisplayRemoveEventListener(ListArr)
		exit(0)
		}
		
	number top, left, bottom, right
	ScatterID.FindImageByID().getannotationrect(AxisID, top,left,bottom,right)
	number dx = right-left
	number dy = bottom-top
	number length = sqrt(dx**2+dy**2)
	dx /= length
	dy /= length
	image load := LoadID.FindImageByID()
	load = dx*load1 +dy*load2	
	}
		
~MSA_LoadingViewer(object self) 
	{
	result("\nLoadingViewer destroyed")
	}
}











/*********************					*
**************************				**
EndMembers ********************************
**************************				**
**********************					*/
















interface EndMembers
{
//Tableau function used in EndMember Starter
void newTableau(object ob);
void loadTableau(object ob);
void StoreTableau(object ob);
void CloseTableau(object self);
TagGroup GetTableauTags(object ob);
number NewScatterPlot(object ob,number First,number Second,number currCluster);
number GetNumberScatterPlots(object ob);
image GetScatterPlot(object ob,number scatt);
void UpdateAllOvals(object self);
image GetLoading(object ob,number member);
image GetFirstLoading(object ob);
image GetLastLoading(object ob);
number GetNumberLoadings(object ob);
void AddEndMember(object ob);
void OptimizeOneMember(object self,image scatter,number minPoints);
void FitALS(Object ob,number lambda);

// EndMember Starter function used in Tableau
void orderLoadings(object ob);
void EnablingButtons(object ob,number State);
void CloseWaiting(object ob,string buttonIdentifier);

//Tableau function used in EndMembersScatter 
void RemoveScatterPlot(object ob,number ScatterID);
}



/***********************************************************************
function Make EndMembers maps
***********************************************************************/
void MSA112_EndMemberMaps(image EndMembersAll, image Tall,TagGroup MatrixTags,number currCluster,string DataName,number MapsCubed)
	{
	string MSAlabel = MSA112_Label()
	string clusterLabel = MSA112_GetClusterLabel(currCluster)
	
	number Edepth = EndMembersAll.ImageGetDimensionSize(0)
	number endMembN = EndMembersAll.ImageGetDimensionSize(1)
	
	//read info
	number width,height
	MatrixTags.MSA112_GetSIWidthHeight(width,height)
	number scaleX,originX,scaleY,originY
	string unitX,unitY
	MatrixTags.MSA112_GetCalibration(scaleX,originX,unitX,"X")
	MatrixTags.MSA112_GetCalibration(scaleY,originY,unitY,"Y")

	//get averageSpec
	image AvSpecWeighted:= MSA112_GetAverageSpectrum(MatrixTags,CurrCluster)
	
	//restore weighting
	image AvSpec := MSA112_RestoreWeightedSpectrum(AvSpecWeighted,MatrixTags,CurrCluster)
	number scaleE,originE
	string unitE
	AvSpec.ImageGetDimensionCalibration(0,originE,scaleE,unitE,1)
		
	//display EndMembers
		AvSpec.showimage()
		AvSpec.setName(DataName+"("+currCluster+") Endmembers")
		imagedisplay imdisp = ImageGetImageDisplay(AvSpec,0)
		imdisp.LinePlotImageDisplaySetLegendShown(1)
		ImageDisplaySetSliceLabelByID(imdisp,ImageDisplayGetSliceIDByIndex(imdisp,0),"average")	
	
		for(number i=0;i<endMembN;i++)
			{
			image spec =EndMembersAll[i,0,i+1,Edepth]
			spec.ImageSetDimensionCalibration(0,originE,scaleE,unitE,1)
			
			imdisp.imageDisplayAddImage(spec,"EM"+(i+1))
			rgbnumber col =MSA112_ClusterColor(i+1)
			imdisp.LinePlotImageDisplaySetSliceComponentColor(i+1,0,red(col),green(col),blue(col))
			}


	//tune clustering
	image RTall
	if(currCluster!=0) 
		{
		RTall = MSA112_ProjectClusteredScore(MatrixTags,Tall,currCluster)
		}
	else RTall = Tall
	
	//make datacube
	image MapCube := realimage(DataName+"("+currCluster+") EM maps",4,width,height,endMembN)
	for(number j=0;j<height;j++)
		{
		MapCube[irow,j,icol] = RTall[j*width,0,(j+1)*width,endMembN]
		}
		
	//normalize	
	image norm := realimage("",4,width,height)
	norm[icol,irow,0] += MapCube
	image norm3D = imageclone(MapCube)
	norm3D = norm[icol,irow,0]
	
	//set out-off-cluater pixels
	number OutCluster = MSA112_OutClusterValue()
	MapCube = tert(modulus(norm3D)>0,MapCube/norm3D,OutCluster)
	
	if(MapsCubed==0) for(number j=0;j<endMembN;j++)
		{
		image Map = MapCube[0,0,j,width,height,j+1]
		
		Map.ImageSetDimensionCalibration(0,OriginX,ScaleX,unitX,1)	
		Map.ImageSetDimensionCalibration(1,OriginY,ScaleY,unitY,1)
		
		Map.ImageSetName(DataName+"("+currCluster+") EM map "+(j+1))
		Map.showimage()
		}
	else 
		{
		MapCube.ImageSetDimensionCalibration(0,OriginX,ScaleX,unitX,1)	
		MapCube.ImageSetDimensionCalibration(1,OriginY,ScaleY,unitY,1)
		
		MapCube.showimage()
		}
		
	}


void MSA112_ConstructMaps(image &EndMembers,image &Abundancies,TagGroup MatrixTags,image Tableau,image LoadAll,image AvSpec,number currCluster,number MSAID,number endMembN,number lambda)
	{	//return EndMembers(all Specs as columns) and maps Abundancies(all maps as columns)
			
	image MSAimage := getImageFromID(MSAID)
	TagGroup MatrixTags = MSA112_GetMatrixTags(MSAimage,0)
	string clusterLabel = MSA112_GetClusterLabel(currCluster)	
	
	//get dataMatrix
	image DataMatrix = MSA112_getMSAMatrix(MatrixTags,currCluster)
	number Depth = DataMAtrix.ImageGetDimensionSize(0)
//result("\npoints "+DataMatrix.ImageGetDimensionSize(1))
	
	//make EndMembers matrix from current endmember Loadings
	EndMembers = MatrixMultiply(MatrixTranspose(Tableau),LoadAll)[0,0,endMembN,depth]
		//[endMembN x depth] matrix
		
	//add average
	EndMembers += AvSpec[icol,0]
	
	//weight
	image EndMembersWeighted = MSA112_ReWeighSpectrum(EndMembers,MatrixTags,currCluster) 
	
	//convert to row-matrix
	image Sall = MatrixTranspose(EndMembersWeighted)			
//Sall.MatrixTranspose().showimage()

	//filter out peculiarities
	Sall = MSA112_FilterSingularities(Sall)
	
	//adjust Lambda
	lambda *= mean(Sall)
	
	//ALS on weighted set
	image UnityMatrix = MSA112_MakeDiagonalMatrix(endMembN,lambda)

	Abundancies = MSA112_ALStep(MatrixTranspose(DataMatrix),Sall,UnityMatrix)	
//Abundancies.MatrixTranspose().showimage()
	}
	

documentWindow MSA112_GetLoadingWindow(image Loading)
	{
	imagedocument doc =loading.ImageGetOrCreateImageDocument()

	if(MSA112_CheckImageWindow(doc))
		{
		return doc.ImageDocumentGetWindow()
		}
	else return NULL	
	}

void MSA112_GetLoadingPositionAndSize(image loading,number &Xpos,number &Ypos,number &width,number &height)	
	{
	documentWindow win = MSA112_GetLoadingWindow(Loading)

	win.WindowGetContentSize(width,height)
	win.WindowGetFramePosition(Xpos,Ypos)
	}

void MSA112_SetLoadingPositionAndSize(image loading,number Xpos,number Ypos,number width,number height)	
	{
	documentWindow win = MSA112_GetLoadingWindow(Loading)

	win.WindowSetContentSize(width,height)
	win.WindowSetFramePosition(Xpos,Ypos)
	}

	
void MSA112_storeLoadingPositionAndSize(image Loading,string path)
	{
	number X,Y,width,height
	Loading.MSA112_GetLoadingPositionAndSize(X,Y,width,height)
	
	TagGroup pers = GetPersistentTagGroup()
	TagGroup TargetTags
	pers.TagGroupGetTagAsTagGroup(MSA112_MSApath(),TargetTags)
	
	TargetTags.MSA112_SetNumberInTags(path+":X","I",X)
	TargetTags.MSA112_SetNumberInTags(path+":Y","I",Y)
	TargetTags.MSA112_SetNumberInTags(path+":width","I",width)
	TargetTags.MSA112_SetNumberInTags(path+":height","I",height)	
	}

void MSA112_restoreLoadingPositionAndSize(image Loading,string path)
	{
	TagGroup pers = GetPersistentTagGroup()
	TagGroup TargetTags
	pers.TagGroupGetTagAsTagGroup(MSA112_MSApath(),TargetTags)

	number X =TargetTags.MSA112_getNumberFromTags(path+":X","I")
	number Y =TargetTags.MSA112_getNumberFromTags(path+":Y","I")
	number width =TargetTags.MSA112_getNumberFromTags(path+":width","I")
	number height =TargetTags.MSA112_getNumberFromTags(path+":height","I")
	
	if((width !=0)&&(height !=0))
		{
		Loading.MSA112_SetLoadingPositionAndSize(X,Y,width,height)
		}
	}	


void MSA112_setOneEndMember(image scatter,image point,number rad,number colCode)
	{	
	//CURRENTLY NOT IN USE IN MAIN PACKAGE
	number width,height
	scatter.getSize(width,height)
	
	number scale,originX,originY
	string units
	scatter.ImageGetDimensionCalibration(0,originX,scale,units,1)
	scatter.ImageGetDimensionCalibration(1,originY,scale,units,1)
	
	string label =MSA112_Label()	//"temDM MSA"
	TagGroup ScatterTags = scatter.ImageGetTagGroup()

	number first=ScatterTags.MSA112_getNumberFromTags(label+":"+"first","S")
	number second =ScatterTags.MSA112_getNumberFromTags(label+":"+"second","S")
//result("\n"+first+"-"+second)	
	
	number X =point.getPixel(first-1,0)
	number Y =point.getPixel(second-1,0)

//result("\nrad "+rad+" "+member)	
	x =(x/scale+originX).round()
	y =(y/scale+originY).round()
//result(" x "+X +" y "+y)

	imagedisplay disp = ImageGetImageDisplay(scatter,0)
	number OvalID = CreateOvalAnnotation(scatter,Y-rad,X-rad,Y+rad,X+rad)
	component Oval = disp.ComponentGetChildbyID(OvalID)
	
	rgbnumber col = MSA112_ClusterColor(colCode)
	number red = red(col)
	number green = green(col)
	number blue = blue(col)
//result("\n"+member+" "+red+"-"+green+"-"+blue)	
	Oval.ComponentSetForegroundColor(red,green,blue)
	Oval.ComponentSetBackgroundColor(red,green,blue)
	Oval.ComponentSetFillMode(1)

	Oval.ComponentSetResizable(0)
	Oval.ComponentSetDeletable(1)
	}

void MSA112_setEndMembers(image scatter,image Endmembers,number upper,number rad,number colCode)
	{
	//CURRENTLY NOT IN USE IN MAIN PACKAGE
	number membN = Endmembers.ImageGetDimensionSize(1)
	number dim = Endmembers.ImageGetDimensionSize(0)-2
	number Limit = membN
	if(Upper!=0) Limit =Upper
	
	for(number i=0;i<Limit;i++)
		{	
		scatter.MSA112_setOneEndMember(EndMembers[i,0,i+1,dim],rad,colCode)
		}
	}
	
	

/***********************************************************************
Class Scatterplot with EndMembers
***********************************************************************/
	
class MSA_EndMembersScatter : MSA112_Frame
{
number TableauID
number ScatterID
number Comp1,Comp2

void init(object self, imagedocument scatterDoc,string Package,number Comp1P,number Comp2P, number TableauIDP)
	{
	TableauID = TableauIDP
	Comp1 = Comp1P
	Comp2 = Comp2P
	
	ScatterID = scatterDoc.ImageDocumentGetImage(0).ImageGetID()
	
	self.MSA112_initFrame(ScatterDoc,Package,"Scatter "+Comp1+"-"+Comp2)	
	}	

~MSA_EndMembersScatter(object self) 
	{
	try{
		object Tableau = GetScriptObjectFromID(TableauID)
		
		//remove Tags
		Tableau.RemoveScatterPlot(ScatterID)
		}
	catch{break;}	
	}
}


/***********************************************************************
Class Loading with EndMember
***********************************************************************/
//CURRENTLY NOT IN USE
class MSA_EndMemberLoad : MSA112_Frame
{
number member
number TableauID

void init(object self, imageDocument loadDoc,string Package,number memberP, number TableauIDP)
	{
	member = memberP
	TableauID = TableauIDP
	
	self.MSA112_initFrame(loadDoc,Package,"Load "+(member+1))	//Loads in List 0,1,2... Load -> 1,2,3...
	}
	
		
~MSA_EndMemberLoad(object self) 
	{
	try{
	//result("\nTableau ID " + TableauID)
		object Tableau = GetScriptObjectFromID(TableauID)
		//Tableau.StoreLoadRangeAndContrast(member)
		//Tableau.RemoveMember(member)
		}
	catch{break;}
	}
}


/***********************************************************************
Class Tableau showing Endmembers mixing
***********************************************************************/

class MSA_EndMembersTableau //: MSA112_Frame
{
/*
EndMembers Tags:

EndMembers:		[0]:	Load ID		//ID of Memeber's Load image 
						low			//min contrast
						high		//max contrast
						left		//left displayed channel
						right		//right displayed channel
				[1] .....
						
ScatterPlots:	[0]:	Scatter ID 	//ID of scatterPLot image			
						Comp 1     	//first PCA component in scatter
						Comp 2	   	//second PCA component in scatter
						Listener ID	//display listener ID
						Members:	[0]:	Oval ID	//ID of the Oval showimng position of a member
									[1] .....

TableauTags.TagGroupOpenBrowserWindow(0)
									
									
Extra Tags in Stored Endmembers

number EndMembers
number ScatterPLots
number components

Tableau		//stored image with contribution of PCA loadings in EndMebers
*/


//TABLAYU DESTROED WITH CLOSING LAST SCATTER PLOT !!

//Setting parameters
string MSAlabel		
number sparsity
number radius 		//radius of Oval
	
number StarterID		//to call Starter 
number MSAID			//to call MSA image
number currCluster		//fixed cluster for given EMTableau instance, columns -endmembers spectra expressed in PCA loadings
number dim		//fixed parameters kept constant in the given instance

image DataMatrix		//full data matrix loaded only when neccesary
TagGroup TableauTags	//kept in memory and copied in Tableau image
image Tableau			//image with contribution of PCA loadings in EndMebers
image AvSpec,LoadAll	//AverageSpec (UNWEIGHTED) and PCA loadings (UNWEIGHTED) for quick access
number scattN,endMembN	//current number of scatterPLots and EndMembers

number cleaningUpper //temporary flag when deleting a load




/* *******  TableauTags ************* */

void CreateTableauTags(object self)
	{	
	TableauTags = NewTagGroup()
	TableauTags.TagGroupSetTagAsTagGroup("ScatterPlots",NewTagList())
	TableauTags.TagGroupSetTagAsTagGroup("EndMembers",NewTagList())
	
	//add it to Tableau
	TagGroup ImageTags = Tableau.ImageGetTagGroup()
	TagGroupSetTagAsTagGroup(ImageTags,MSA112_Label(),TableauTags)	
	}

		
TagGroup GetTableauTags(object self) return TableauTags




/* ******* Scatter plots ************* */
		
image GetScatterPlot(object self,number scatt)
	{		
	number scatterID = TableauTags.MSA112_getNumberfromList("ScatterPlots",scatt,"Scatter ID","L")
//result("\nscattID "+scatterID)
	try{	
		image scatter := getImageFromID(scatterID)
		return scatter
		}
	catch{break;}
	
	return NULL
	}
	
number GetNumberScatterPlots(object self) return scattN
	
number FindScatterIndex(object self,image scatter) 
	{
	number scatterID = scatter.ImageGetID()
	for (number i = 0;i<scattN;i++)
		{
		if(scatterID==TableauTags.MSA112_getNumberfromList("ScatterPlots",i,"Scatter ID","L")) return i
		}
	return (-1)	
	}


number CheckIfScatterExist(object self,number first,number second) 
	{
	for (number i = 0;i<scattN;i++)
		{
		number Comp1 = TableauTags.MSA112_getNumberfromList("ScatterPlots",i,"Comp 1","S")
		number Comp2 = TableauTags.MSA112_getNumberfromList("ScatterPlots",i,"Comp 2","S")
//result("\nscatter "+i+" "+Comp1+"-"+Comp2+ " new "+first+"-"+second)		
		number exist = ((comp1==first)&&(comp2==second)) || ((comp1==second)&&(comp2==first))
//result(" exist "+exist)		
		if(exist)	return 0
		}
		
	return 1
	}
	
void GetScatterCalibration(object self,image scatter,number &originX,number &originY,number &scale) 
	{
	string unit
	scatter.ImageGetDimensionCalibration(0,originX,scale,unit,1)
	scatter.ImageGetDimensionCalibration(1,originY,scale,unit,1)
	}

number GetOvalCoord(object self,number scale,number origin,number dx) 
	{
	return (origin + dx/scale)
	}

number GetOvalCoordCalibrated(object self,number scale,number origin,number x) 
	{
	return ((x-origin)* scale)
	}
	
number createOval(object self,imageDisplay ScatDisp,number member,number X,number Y) 
	{
	image scatter := ScatDisp.ImageDisplayGetImage()
	number OvalID = CreateOvalAnnotation(scatter,Y-radius,X-radius,Y+radius,X+radius)
//result("\novalID "+ovalID)	
	component Oval = ScatDisp.ComponentGetChildbyID(OvalID)
	
	rgbnumber col = MSA112_ClusterColor(member+1)
	number red = red(col)
	number green = green(col)
	number blue = blue(col)

	number colCodeB 
	//6 fullcolor	- 6 color with blackrand	- 6 color with red rand
	if(member<6) colCodeB =(member+1)
	if((member>=6)&&(member<12)) colCodeB =0
	if(member==12) colCodeB =3
	if((member>12)&&(member<18)) colCodeB =1

	rgbnumber colB = MSA112_ClusterColor(colCodeB)
	number redB = red(colB)
	number greenB = green(colB)
	number blueB = blue(colB)
	
	Oval.ComponentSetForegroundColor(red,green,blue)
	Oval.ComponentSetBackgroundColor(redB,greenB,blueB)
	Oval.ComponentSetFillMode(1)
	Oval.ComponentSetResizable(0)
	Oval.ComponentSetDeletable(0)
	
	//make mark on Endmember 
	tagGroup OvalTags = Oval.ComponentGetTagGroup()
	OvalTags.MSA112_setNumberInTags("EndMember","S",member)
	
	return OvalID
	}	

number SetOvalInScatter(object self,number scatt,number member) 
	{			
	image scatter := self.GetScatterPlot(scatt)	
	imagedisplay ScatDisp = scatter.ImageGetImageDisplay(0)	
	
	number first = TableauTags.MSA112_getNumberfromList("ScatterPlots",scatt,"Comp 1","S")	
	number second = TableauTags.MSA112_getNumberfromList("ScatterPlots",scatt,"Comp 2","S")

//result(" "+member+" "+endMembN+" "+first+"-"+second)
	number originX,originY,scale
	self.GetScatterCalibration(scatter,originX,originY,scale)
	number X = self.GetOvalCoord(scale,originX,Tableau.getPixel(member,first-1))
	number Y = self.GetOvalCoord(scale,originY,Tableau.getPixel(member,second-1))
//result(" x "+x+ " y "+y)

	//set position of existing Oval, creates oval if ot exists
	number ovalID = TableauTags.MSA112_getNumberfromList("ScatterPlots:["+scatt+"]:Members",member,"Oval ID","L")
	component Oval 
	try {Oval = ScatDisp.ComponentGetChildbyID(OvalID);} catch {break;}
	if(Oval.ComponentIsValid())
		{
		Oval = ScatDisp.ComponentGetChildbyID(OvalID)
//result("\novalID "+ovalID)		
		Oval.ComponentSetRect(Y-radius,X-radius,Y+radius,X+radius)
		}
	else
		{
		//create Oval	
		OvalID = self.CreateOval(ScatDisp,member,X,Y)

		//write oval ID in scatter tags
		TableauTags.MSA112_setNumberinList("ScatterPlots:["+scatt+"]:Members",member,"Oval ID","L",OvalID)		
		}
 		
	return OvalID
	}
		
void AddScatterPlot(object self,number CurrCluster,number first,number second) 
	{
	//add new Scatter Plot but without Ovals
	if(self.CheckIfScatterExist(first,second)) 
		{
//result("\nadd scatterplot "+first +"-"+second)
		image MSAimage:= getImageFromID(MSAID)
		TagGroup MatrixTags = MSA112_GetMatrixTags(MSAimage,0)

		//create scatterPLot
		string Method = "PCA"	
		rgbimage scatter := MSA112_PlotXY(MSAimage,currCluster,Method,first,second,sparsity)
		scatter.showimage()

		//add listener for changing Ovals
		imagedisplay ScatDisp = scatter.ImageGetImageDisplay(0)	
		number ListenerID = ScatDisp.ImageDisplayAddEventListener(self,"component_changed,component_removed:MemberChanged")
	
		//wrap into Frame112
		imageDocument Scatdoc = scatter.ImageGetOrCreateImageDocument()
		Scatdoc.ImageDocumentShow()
		alloc(MSA_EndMembersScatter).init(Scatdoc,MSA112_MSApath()+":EndMembers",first,second,self.ScriptObjectGetID())
result(" (ScatterPlot "+first+"-"+second+")")	
	
		//make scatter tags
		TableauTags.MSA112_setNumberinList("ScatterPlots",scattN,"Scatter ID","L",scatter.imageGetID())
		TableauTags.MSA112_setNumberinList("ScatterPlots",scattN,"Listener ID","L",ListenerID)
		TableauTags.MSA112_setNumberinList("ScatterPlots",scattN,"Comp 1","S",first)
		TableauTags.MSA112_setNumberinList("ScatterPlots",scattN,"Comp 2","S",second)
		
		//create empty member list
		TableauTags.TagGroupSetTagAsTagGroup("ScatterPlots:["+scattN+"]:Members",NewTagList())	
		
		scattN ++
		}
	}		


void LoadExistingEndMembers(object self,number scatt) 
	{
	for (number j = 0;j<endMembN;j++)
		{
//result("\nadd oval "+j)
		self.SetOvalInScatter(scatt,j)
		}
	}

	
number NewScatterPlot(object self,number First,number Second,number currCluster) 
	{
	if(!self.CheckIfScatterExist(First,Second)) return 0

	self.AddScatterPlot(CurrCluster,first,second) 	//new ScatterPlot
	self.LoadExistingEndMembers(scattN-1)			//new ovals
	
	return 1
//TableauTags.TagGroupOpenBrowserWindow(0)	
	}



/* ******* EndmembersLoadings ************* */
	
image GetLoading(object self, number member)
	{
	number LoadID = TableauTags.MSA112_getNumberfromList("EndMembers",member,"Load ID","L")	
	try{
//result("\ndeleted loading "+LoadID)	
		return getImageFromID(LoadID)
		}
	catch{break;}

	return NULL
	}

number FindLoadingIndex(object self,image Loading) 
	{
	number LoadingID = Loading.ImageGetID()

	for (number i = 0;i<endMembN;i++)
		{		
		if(LoadingID==TableauTags.MSA112_getNumberfromList("EndMembers",i,"Load ID","L")) return i
		}
	return (-1)	
	}
	
image GetFirstLoading(object self) return self.GetLoading(0)

image GetLastLoading(object self) return self.GetLoading(endMembN-1)	

number GetNumberLoadings(object self) return endMembN
	
image UpdateLoad(Object self,number member)
	{
	image contribution = Tableau[0,member,dim,member+1]
	
	return (AvSpec + MatrixMultiply(MatrixTranspose(contribution),LoadAll))
	}
	

//creates spectrum profile	
void CreateEndMember(object self,number member) 
	{	
	//create load
	image load := imageclone(AvSpec)
	load.setName("endMember "+(member+1))
	ImageDocument loadDoc = CreateImageDocument("New ImageDocument")
	ImageDisplay imDisp = loadDoc.ImageDocumentAddImageDisplay(load,3)
	rgbnumber col = MSA112_ClusterColor(member+1)
	imDisp.LinePlotImageDisplaySetSliceComponentColor(0,1,red(col),green(col),blue(col))
	loadDoc.ImageDocumentShow()
	
	//attach listener
	DocumentWindow LoadWin = MSA112_GetLoadingWindow(Load)	
	number Listener = LoadWin.WindowAddWindowListener(self,"window_closed:RemoveMember")
//result("\nListener "+Listener)
	//set range 
	number left = TableauTags.MSA112_getNumberfromList("EndMembers",member,"left","I")	
	number right = TableauTags.MSA112_getNumberfromList("EndMembers",member,"right","I")
//result("\ncontrast low "+low+" high "+high) 
	if((left==0)&&(right==0)) {;}	
	else	imDisp.LinePlotImageDisplaySetDisplayedChannels(left,right)
	
	//set contrast	
	number low  = TableauTags.MSA112_getNumberfromList("EndMembers",member,"low","F")
	number high	= TableauTags.MSA112_getNumberfromList("EndMembers",member,"high","F")		
	if((low==0)&&(high==0)) {;}
	else
		{
		imDisp.LinePlotImageDisplaySetDoAutoSurvey(0,0)
		imDisp.LinePlotImageDisplaySetContrastLimits(low,high)
		}
		
	//update Load
	load = self.UpdateLoad(endMembN)
	
//result("\nEndMember "+(member+1)+" created")	

	//make member tags
	TableauTags.MSA112_setNumberinList("EndMembers",member,"Load ID","L",load.imageGetID())	
	TableauTags.MSA112_setNumberinList("EndMembers",member,"Listener ID","L",Listener)		
	}
	
	
//creates spectrum and oval	
void AddEndMember(object self) 
	{	
	//create ovala in all scatterPLots
	for (number i = 0;i<scattN;i++)
		{
		self.SetOvalInScatter(i,endMembN)
		}	

	//create load profileb
	self.CreateEndMember(endMembN)
	
	endMembN ++
	
	//update member information in EM tool
	object Starter = GetScriptObjectFromID(StarterID)
	Starter.lookupElement("EndMembN").DLGTitle(": "+endMembN)
	}


void UpdateEndMembers(object self,number member,number first,number second) 
	{	
	//scroll over ScatterPLots
	for (number i = 0;i<scattN;i++)
		{
//result("\nupdate scatt "+i)
		number Comp1 = TableauTags.MSA112_getNumberfromList("ScatterPlots",i,"Comp 1","S")
		number Comp2 = TableauTags.MSA112_getNumberfromList("ScatterPlots",i,"Comp 2","S")

		//check if scatterplot consists of changed components
		number any = (comp1==first)||(comp2==first)||(comp1==second)||(comp2==second)
		number all = (comp1==first)&&(comp2==second)
		if(any&&(!all))	//changes in scatterplots where components are mentioned
						//EXCEPT the original scatterplot that triggered change 
			{
			self.SetOvalInScatter(i,member)
			}
		}	
	}



/* ******* initiation ************* */
number init(object self, number MSAIDP,number CurrClusterP,number StarterIDP,number sparsityP,number radiusP)
	{
	MSAlabel = MSA112_Label()	
	sparsity = sparsityP
	radius = radiusP
	MSAID =MSAIDP
	
	StarterID = StarterIDP
	currCluster = currClusterP
	
	//enable Tuning Box in starter
	StarterID.getScriptObjectFromID().EnablingButtons(1)
	
	//get MSA tags
	image MSAimage := getImageFromID(MSAID)
	string DataName = MSAImage.getName()
	TagGroup MatrixTags = MSA112_GetMatrixTags(MSAimage,0)
	//MSAID = MSAimage.ImageGetID()
	string clusterLabel = MSA112_GetClusterLabel(currCluster)	
	
	//read info from tags
	string ReadMethod = "PCA"
	dim = MSA112_ReadNumberComponents(MatrixTags,currCluster,ReadMethod)
	if(!MSA112_CheckComponents(dim)) throw("")
	image AvSpecWeighted := MSA112_GetAverageSpectrum(MatrixTags,CurrCluster)		
	image LoadAllWeighted = matrixTranspose(MSA112_GetLoadings(MatrixTags,ReadMethod,currCluster))
	
	//correct for weighting
	AvSpec := MSA112_RestoreWeightedSpectrum(AvSpecWeighted,MatrixTags,currCluster)
	LoadAll = MSA112_RestoreWeightedSpectrum(LoadAllWeighted,MatrixTags,currCluster)
		
	//make Tableau
	Tableau := realimage("PCA Fractions",4,2*dim,dim)
	
	//Make Tableau tags
	self.CreateTableauTags()	

	number selfID =self.ScriptObjectGetID()
result("\nTableau "+ selfID +" initialized")
	
	return selfID
	}
	

~MSA_EndMembersTableau(object self) 
	{
	//disable Tuning Box in starter
	object Starter = StarterID.getScriptObjectFromID()
	if(Starter.ScriptObjectIsValid()) Starter.EnablingButtons(0)
//result("\nMSA_EndMembersTableau destroyed")
	
result("\nTableau "+ self.ScriptObjectGetID()+" destroyed")	
	}
	
TagGroup StoredTags(object self)
	{	
	image MSAimage := MSAID.getImageFromID()

	string DataName = MSAimage.getName()
	TagGroup MatrixTags = MSA112_GetMatrixTags(MSAimage,0)
	string 	clusterLabel = MSA112_GetClusterLabel(currCluster)	
			
	TagGroup StoredEndMembers
	MatrixTags.TagGroupGetTagAsTagGroup(ClusterLabel+":EndMembering",StoredEndMembers)
	
	return StoredEndMembers
	}
	
void newTableau(object self)
	{	
	if(self.StoredTags().TagGroupIsValid())
		{
		if(!OKCanceldialog("There are already endmembers stored.\nDo you want to design from scratch ?")) throw("")
		}

	//make default scatterplot
	self.AddScatterPlot(CurrCluster,1,2)
	
	//add default two endMembers
	for(number j= 0;j<2;j++){self.AddEndMember();}	
	}
	
void NoStored(object self)
	{
	OKdialog("No stored endmembers !!")
	throw("")	
	}
	
void loadTableau(object self)
	{
	TagGroup StoredEndMembers = self.StoredTags()
	
	if(!StoredEndMembers.TagGroupIsValid()) self.NoStored()
	
	//load Tableau
	Tableau := StoredEndMembers.MSA112_getImageFromTags("Tableau")

	//take number of scatters and endMembers
	TagGroup Members,Scatters
	if(!StoredEndMembers.TagGroupGetTagAsTagGroup("EndMembers",Members)) self.NoStored()
	if(!StoredEndMembers.TagGroupGetTagAsTagGroup("ScatterPlots",Scatters)) self.NoStored()
	number storedScattN = Scatters.TagGroupCountTags()
	number storedEndMembN = Members.TagGroupCountTags()

	//copy to TableauTags
	TagGroupCopyTagsFrom(TableauTags,StoredEndMembers)	
	
	//make scatterPlots
	for (number i = 0;i<storedScattN;i++)
		{
		number comp1 =StoredEndMembers.MSA112_getNumberFromTags("ScatterPlots:["+i+"]:Comp 1","S") 
		number comp2 =StoredEndMembers.MSA112_getNumberFromTags("ScatterPlots:["+i+"]:Comp 2","S") 
		
		self.AddScatterPlot(CurrCluster,Comp1,Comp2)		
		}
	
	//make members
	for (number j = 0;j<storedEndMembN;j++)
		{
		self.AddEndMember()
		}
	}


image getTableau(object self)
	{
	return Tableau
	}
	
	
	
/* **************** REMOVING ***************/

//This remove ScatterPlot tags if the Scatter image was closed
void RemoveScatterPlot(Object self, number ScatterID)
	{  
//result("\nremove scatterID "+scatterID)
	number ScatterIndex = self.FindScatterIndex(getImageFromID(ScatterID))
//result(" index "+ScatterIndex)
	TagGroup ScatterPlots
	TableauTags.TagGroupGetTagAsTagGroup("ScatterPlots",ScatterPlots)
	ScatterPlots.TagGroupDeleteTagWithIndex(ScatterIndex)
	
	scattN --
	}


	
void RemoveOval(Object self) //only last member
	{
	//clean up Oval annotations in all scatterplots
	for (number i = 0;i<scattN;i++)
		{
		//find scatterplot
		image scatter 
		imageDisplay ScatDisp
		try	{
			scatter:= self.GetScatterPlot(i)
			ScatDisp = scatter.ImageGetImageDisplay(0)
			number OvalID = TableauTags.MSA112_getNumberfromList("ScatterPlots:["+i+"]:Members",EndMembN-1,"Oval ID","L")
			scatter.DeleteAnnotation(ovalID)
			
			TagGroup Ovals
//result("\ndelete ovalTags in "+i) 		
			TableauTags.TagGroupGetTagAsTagGroup("ScatterPlots:["+i+"]:Members",Ovals)		
			Ovals.TagGroupDeleteTagWithIndex(EndMembN-1)	
			}
		catch{break;}
		}
	}


void StoreLoadRangeAndContrast(object self,number member)	
	{
	try{
		image load := self.GetLoading(member)
		imageDisplay imdisp = load.ImageGetImageDisplay(0)
		
		number low,high,left,right
		imdisp.LinePlotImageDisplayGetContrastLimits(low,high)
		//result("\ncontrast low "+low+" high "+high) 	
		imdisp.LinePlotImageDisplayGetDisplayedChannels(left,right)
		//result("\nTableau ID " + TableauID)
			
		TableauTags.MSA112_setNumberinList("EndMembers",member,"low","F",low)	
		TableauTags.MSA112_setNumberinList("EndMembers",member,"high","F",high)	
		TableauTags.MSA112_setNumberinList("EndMembers",member,"left","I",left)	
		TableauTags.MSA112_setNumberinList("EndMembers",member,"right","I",right)	
		}
	catch{break;}
	}
		
void RemoveMember(Object self)	//only last member
	{
	//oval always deleted if loading closed
	try{self.RemoveOval();} catch{break;}
	
	//reduce and 
	number ListenerID = TableauTags.MSA112_getNumberFromList("EndMembers",EndMembN-1,"Listener ID","L")
	DocumentWindow LoadWin = MSA112_GetLoadingWindow(self.getLoading(EndMembN-1))
	LoadWin.WindowRemoveWindowListener(ListenerID)
	
	//rearrange spectra	
	try{
		image loading := self.GetLastLoading()
		loading.deleteImage()
		}
	catch{break;}
			
	EndMembN --	
	}

void RemoveMember(Object self, number member)	//this is for any member, also in the middle
	{	
	number endMembNOld = endMembN
		
	//clean member from Tableaue, rearrange Tableau
	number dim = Tableau.ImageGetDimensionSize(1)
	if(member<(endMembNOld-1)) //member in the middle
		{
		Tableau[0,member,dim,endMembNOld-1] = Tableau[0,member+1,dim,endMembNOld]
		}
	Tableau[0,endMembNOld-1,dim,endMembNOld] =0

	object Starter = GetScriptObjectFromID(StarterID)
	
	if(Starter.ScriptObjectIsValid())
		{
		//all upper members including the target member are deleted
		for(number j=(endMembNOld-1);j>=member;j--)	//if member =(EndMembNold-1) then only once
			{	
			self.StoreLoadRangeAndContrast(j)
			self.removeMember()
			}
//result(" removed")
		}
		
	//remove Tags of deleted member (the rest Tags are rearranged automatically)
	TagGroup MembersTags
	TableauTags.TagGroupGetTagAsTagGroup("EndMembers",MembersTags)
	MembersTags.TagGroupDeleteTagWithIndex(member)

	if(Starter.ScriptObjectIsValid())
		{	
		//restore upper members
		for(number j= (member+1);j<endMembNOld;j++) {self.AddEndMember();} 
						//if Last loading, member =(EndMembNold-1) then no adding

		//restore loadings positions
		if(endMembN>=1)
			{
			MSA112_restoreLoadingPositionAndSize(self.GetFirstLoading(),"EndMembers:first Loading")	
			
			if(endMembN>=2) //otherwise the last loading is the first one
				{
			MSA112_restoreLoadingPositionAndSize(self.GetLastLoading(),"EndMembers:last Loading")
				}
			}
		}	
	
	if(Starter.ScriptObjectIsValid())
		{
		Starter.lookupElement("EndMembN").DLGTitle(""+EndMembN)	
		if(endMembN>2)Starter.orderLoadings();
		}
		
//TableauTags.TagGroupOpenBrowserWindow(0)	
	}
	
	
//when loading window is manually closed
void RemoveMember(object self, number e_fl, DocumentWindow Win)
	{
	imagedocument doc =Win.ImageWindowGetImageDocument()
	image loading :=doc.ImageDocumentGetImage(0)
	number member = self.FindLoadingIndex(Loading)
//result("\nwin closed "+member)	
	self.RemoveMember(member)
	}

void CloseTableau(object self)
	{
	number scattNold =scattN
	number endMembNold =endMembN

	for (number i = (endMembNold-1);i>=0;i--)
		{
		self.removeMember()		
		}	
		
	for (number i = (scattNold-1);i>=0;i--)	//descending becouse TagList is reduced from end
		{
		number ListenerID = TableauTags.MSA112_getNumberfromList("ScatterPlots",i,"Listener ID","L")
//result("\nscatt "+i+" listener "+ListenerID)
		image scatter := self.GetScatterPlot(i)	
		imagedisplay ScatDisp = scatter.ImageGetImageDisplay(0)	
		ScatDisp.ImageDisplayRemoveEventListener(ListenerID)
//result(" removed")		
		scatter.deleteImage()
		}
		


	Tableau.deleteImage()	
	}


/* ********************* LIVE TUNE **********************/
	
void MemberChanged(Object self, Number e_fl, ImageDisplay disp, Number c_fl, Number c_fl2, Component c )
	{
	//retrieve member
	tagGroup OvalTags = c.ComponentGetTagGroup()
	number member
	if(OvalTags.TagGroupGetTagAsShort("EndMember",member))	//exceptional: TagGroupGetTagAsShort
		{
		//get scatterplot
		image scatter := disp.ImageDisplayGetImage()
		TagGroup ScatterTags = scatter.ImageGetTagGroup()	
		number first =ScatterTags.MSA112_getNumberFromTags(MSAlabel+":first","S")
		number second =ScatterTags.MSA112_getNumberFromTags(MSAlabel+":second","S")
		
		if(e_fl==33554432) 
			{
			//deleted Oval annotation		
			}
		else
			{	
			number scatterID = disp.ImageDisplayGetImage().imageGetID() 	
			if(scatterID == GetFrontImageID())	
				{//only react on moving ovals in frontmost image 
				
				//get changes in loading
				number top, left, bottom, right
				c.ComponentGetRect(top, left, bottom, right)
				number originX,originY,scale
				self.GetScatterCalibration(scatter,originX,originY,scale)
				number dX = self.GetOvalCoordCalibrated(scale,originX,(left+right)/2)
				number dY = self.GetOvalCoordCalibrated(scale,originY,(top+bottom)/2)
	//result("\n"+(left+right)/2+" "+(top+bottom)/2+" dx "+dx+ " dy "+dy)

				//update Tableau
				Tableau.setPixel(member,first-1,dX)
				Tableau.setPixel(member,second-1,dY)	
							
				//update loading				
				number LoadID = TableauTags.MSA112_getNumberfromList("EndMembers",member,"Load ID","L")
				image load := LoadID.FindImageByID()	
//result("\nmember "+member+" EndmemberID "+LoadID +" valid "+load.ImageIsValid())
				load = self.UpdateLoad(member)

				self.UpdateEndMembers(member,first,second)
				}
			}
		}	
	}	

	

void UpdateAllOvals(object self)
	{
	//update all Endmembers
	for (number i = 0;i<scattN;i++)
		{
//result("\nscatter "+i)		
		image scatter := self.GetScatterPlot(i)
		scatter.ImageGetOrCreateImageDocument().ImageDocumentGetWindow().WindowSelect()
		for (number j = 0;j<endMembN;j++)
			{
//result("   member "+j)
			self.SetOvalInScatter(i,j)
			}
		}
	}

void OptimizeOneMember(object self,image scatter,number minPoints)
	{
	//endmember is looling for in only one scatter plot
	number Scatt = self.FindScatterIndex(scatter)

	number Comp1 = TableauTags.MSA112_getNumberfromList("ScatterPlots",Scatt,"Comp 1","S")
	number Comp2 = TableauTags.MSA112_getNumberfromList("ScatterPlots",Scatt,"Comp 2","S")
	
	imagedisplay ScatDisp = scatter.ImageGetImageDisplay(0)
	component Oval
	number member
	for(member=0;member<endMembN;member++)
		{
		number ovalID = TableauTags.MSA112_getNumberfromList("ScatterPlots:["+scatt+"]:Members",member,"Oval ID","L")
		Oval = ScatDisp.ComponentGetChildbyID(OvalID)

		if(ComponentIsSelected(Oval)) break		//selected endmember is found
		if(member == (endMembN-1)) {OKdialog("select endMember marker !!");exit(0);}
		}		
//result("\nmember "+member+" comp "+Comp1+" "+Comp2)	
	
	//get position of the selected member
	number top, left, bottom, right
	Oval.ComponentGetRect(top, left, bottom, right)
	number cX =(left+right)/2
	number cY =(top+bottom)/2					

	//get PCA scores
	image MSAimage := getImageFromID(MSAID)
	TagGroup MatrixTags = MSA112_GetMatrixTags(MSAimage,0)
	string Method ="PCA"
	number currCluster =0
	image scores =MSA112_GetScores(MatrixTags,Method,currCluster)
	number points = scores.ImageGetDimensionSize(1)
	image score1 = scores[0,Comp1-1,points,Comp1]
	image score2 = scores[0,Comp2-1,points,Comp2]	
	image cluster = MSA112_GetCluster(MatrixTags)
//score1.matrixTranspose().showimage()
//score2.matrixTranspose().showimage()	
	image mask =cluster.imageclone()
	
	//get calibration
	number originX,originY,scale
	self.GetScatterCalibration(scatter,originX,originY,scale)

	//try to extend area around target pixel to get min number of points
	number delPix =1
	image redScores,Aver	//reduced scores and their average	
	for(number i=1;i<10;i++)	//maximal extending til +-10pixels in scatted plot
		{
		number Xmin = self.GetOvalCoordCalibrated(scale,originX,cX-delPix)
		number Ymin = self.GetOvalCoordCalibrated(scale,originY,cY-delPix)
		number Xmax = self.GetOvalCoordCalibrated(scale,originX,cX+delPix)
		number Ymax = self.GetOvalCoordCalibrated(scale,originY,cY+delPix)
//result("\nXmin "+Xmin+" Xmax "+Xmax+" Ymin "+Ymin+" Ymax "+Ymax)
		
		mask =tert((score1>Xmin)&&(score1<Xmax)&&(score2>Ymin)&&(score2<Ymax),1,0)
	//mask.matrixTranspose().showimage()	
		
		number kPoints
		if((max(mask)>0))
			{
			redScores := MSA112_TakeCluster(scores,mask,1,0,Aver,kPoints) //return reddataMatrix (NOT used)
			}
//result("\nkPoints "+kPoints)
		if(kPoints <minPoints) delPix++
		else break
		}
		
	//set the most dense scores	
	number dim = Tableau.ImageGetDimensionSize(1)	
	Tableau[0,member,dim,member+1] = Aver.MatrixTranspose()
	
	self.UpdateAllOvals()
	}
	
	
void FitALS(Object self,number lambda)
	{
	image MSAimage := getImageFromID(MSAID)
	TagGroup MatrixTags = MSA112_GetMatrixTags(MSAimage,0)

	//get dataMatrix
	if(!DataMatrix.ImageIsValid()) DataMatrix = MSA112_getMSAMatrix(MatrixTags,currCluster)

	//make EndMembersMatrix matrix (only active endMember counted)
	number depth = LoadAll.ImageGetDimensionSize(0)
	image EndMembersAll = MatrixMultiply(MatrixTranspose(Tableau),LoadAll)[0,0,endMembN,depth]
	//EndMembersAll.showimage()	
	
	//add average
	EndMembersAll += AvSpec[icol,0]
//EndMembersAll.showimage()

	//weighting
	number weighting = MatrixTags.MSA112_getnumberFromTags(MSAlabel+":weighted","S")
	if(weighting) EndMembersAll = MSA112_ReWeighSpectrum(EndMembersAll,MatrixTags,currCluster)
	
	//convert to row-matrix
	image Sall = MatrixTranspose(EndMembersAll)			
		
	Sall = MSA112_FilterSingularities(Sall)
		
	//adjust Lambda
	lambda *= mean(Sall)
//result("\nlambda "+lambda)

	//ALS on weighted set
	image UnityMatrix = MSA112_MakeDiagonalMatrix(endMembN,lambda) 
	image Mall = MSA112_ALStep(MatrixTranspose(DataMatrix),Sall,UnityMatrix)
	Sall = MSA112_ALStep(DataMatrix,Mall,UnityMatrix)
	//Sall.imageclone().showimage()	
		
	//convert back to column matrix
	EndMembersAll = MatrixTranspose(Sall)	

	//subtract weighted Average
	image AvSpecWeighted = MSA112_ReWeighSpectrum(AvSpec,MatrixTags,currCluster)
	EndMembersAll -= AvSpecWeighted[icol,0]
	
	//project to weighted PCA Loads
	image LoadAllWeighted 
	if(weighting) LoadAllWeighted = MSA112_ReWeighSpectrum(LoadAll,MatrixTags,currCluster)
	else LoadAllWeighted = LoadAll


	LoadAllWeighted = MSA112_FilterSingularities(LoadAllWeighted)
	EndMembersAll = MSA112_FilterSingularities(EndMembersAll)
	//update Tableau
	image correctedTableau = MatrixMultiply(LoadAllWeighted,MatrixTranspose(EndMembersAll))
//correctedTableau.imageclone().showimage()
	Tableau[0,0,dim,endMembN] = CorrectedTableau
		
	//update all Endmembers
	self.UpdateAllOvals()
		
	StarterID.getScriptObjectFromID().CloseWaiting("ALS")	
	}


	
void MakeMaps(Object self,image &EndMembers,image &Abundancies,number lambda)
	{	
	image MSAimage := getImageFromID(MSAID)
	TagGroup MatrixTags = MSA112_GetMatrixTags(MSAimage,0)

	//construct EndMembers maps
	MSA112_ConstructMaps(EndMembers,Abundancies,MatrixTags,Tableau,LoadAll,AvSpec,currCluster,MSAID,endMembN,lambda)
//Abundancies.showimage()	
//EndMembers.showimage()
			
	MSA112_setImageinTags(TableauTags,"EMSpectra",EndMembers)
	MSA112_setImageinTags(TableauTags,"Abundancies",Abundancies)
//TableauTags.TagGroupOpenBrowserWindow(0)		
	}

	
void StoreTableau(Object self)
	{
	string clusterLabel = MSA112_GetClusterLabel(currCluster)	

	image MSAimage := MSAID.getImageFromID()
	TagGroup MatrixTags = MSA112_GetMatrixTags(MSAimage,0)	
	
	//store range and contrast of loadings in TableauTags
	for (number j = 0;j<EndMembN;j++)
		{
		self.StoreLoadRangeAndContrast(j)
		}
	
	//copy endMembers tags
	TagGroup StoredEndMembers = NewTagGroup()
	TagGroupCopyTagsFrom(StoredEndMembers,TableauTags)

	//extend TableauTags with Tableau image
	StoredEndMembers.MSA112_setImageinTags("Tableau",Tableau)
	
	MatrixTags.TagGroupSetTagAsTagGroup(clusterLabel+":EndMembering",StoredEndMembers)
	}	

}



/***************************End of Endmembers classes ******************************/






/*******************************************
Plotting functions
*******************************************/



number MSA112_plotOneEndMember(image scatter,image point,number fill,number move,number resize, number delet,number rad,rgbnumber col)
	{	
	number width,height
	scatter.getSize(width,height)
	
	number scale,originX,originY
	string units
	scatter.ImageGetDimensionCalibration(0,originX,scale,units,1)
	scatter.ImageGetDimensionCalibration(1,originY,scale,units,1)
	
	string label =MSA112_Label()	//"temDM MSA"
	TagGroup ScatterTags = scatter.ImageGetTagGroup()

	number first=ScatterTags.MSA112_getNumberFromTags(label+":"+"first","S")
	number second =ScatterTags.MSA112_getNumberFromTags(label+":"+"second","S")
//result("\n"+first+"-"+second)	
	
	number dim =point.ImageGetDimensionSize(0)
	number X,Y
	if(first<=dim) X =point.getPixel(first-1,0)
	else X=0
	if(second<=dim) Y =point.getPixel(second-1,0)
	else Y=0
	
//result("\nrad "+rad+" "+member)	
	x =(x/scale+originX).round()
	y =(y/scale+originY).round()
//result(" x "+X +" y "+y)

	imagedisplay disp = ImageGetImageDisplay(scatter,0)

	number red = red(col)/256
	number green = green(col)/256
	number blue = blue(col)/256
	
	number OvalID = CreateOvalAnnotation(scatter,Y-rad,X-rad,Y+rad,X+rad)
	component Oval = disp.ComponentGetChildbyID(OvalID)
		
	//result("\n"+member+" "+red+"-"+green+"-"+blue)	
	Oval.ComponentSetForegroundColor(red,green,blue)
	if(move) Oval.ComponentSetBackgroundColor(red,green,blue) 	//make radius solid when GLOBAL background filled
	number fillMode												//otherwise always thin even if background filled
	if(fill) fillMode =1
	else fillMode =2
	Oval.ComponentSetFillMode(fillMode)	
	Oval.ComponentSetSelectable(move)	//move implies selectable
	Oval.ComponentSetMovable(move)
	Oval.ComponentSetResizable(resize)
		
	return OvalID
	}

number MSA112_plotOneEndMember(image scatter,image point,number fill,number move,number resize, number delet,number rad,number colCode)
	{
	rgbnumber col =MSA112_ClusterColor(ColCode)
	
	return MSA112_plotOneEndMember(scatter,point,fill,move,resize,delet,rad,col)
	}
	
image MSA112_plotEndMembers(image scatter,image Endmembers,number fill,number move,number rad,number colCode)
	{
	number membN = Endmembers.ImageGetDimensionSize(1)
	number dim = Endmembers.ImageGetDimensionSize(0)-2

	image OvalIDs := realimage("",4,1,membN)
	
	for(number i=0;i<membN;i++)
		{	
		rgbnumber col
		if(ColCode==(-1)) col =MSA112_ClusterColor(i+1)
		else 	col = MSA112_ClusterColor(colCode)
		
		number ID =scatter.MSA112_plotOneEndMember(EndMembers[i,0,i+1,dim],fill,move,0,0,rad,col)
		OvalIDs.setPixel(0,i,ID)
		}
		
	return OvalIDs	
	}
	
number MSA112_PlotRadius(image Scatter, number Rad,number move,number resize, number delet, number colCode)
	{
	image NullVec :=realimage("",4,100,1)
	NullVec =0//Rad
	
	number scaleScatter =Scatter.ImageGetDimensionScale(0)
	Rad /=scaleScatter 
			
		number fill =0
	number RadID =Scatter.MSA112_plotOneEndMember(NullVec,fill,move,resize,delet,Rad,MSA112_ClusterColor(colCode))
	
	return RadID
	}

void MSA112_ReplotRadius(image scatter,number Rad,number RadID)
	{
	imagedisplay disp = ImageGetImageDisplay(scatter,0)

	component Oval = disp.ComponentGetChildbyID(RadID)
	number top,left,bottom,right
	Oval.ComponentGetRect(top,left,bottom,right)
	number X = (left+right)/2
	number Y = (top+bottom)/2

	number scale =scatter.ImageGetDimensionScale(0)
	Rad  /=scale

	Oval.ComponentSetRect(Y-rad,X-rad,Y+rad,X+rad)
	}


void MSA112_removeOneEndMember(image scatter,number OvalID)

	{
	Scatter.deleteAnnotation(OvalID)
	}
		
void MSA112_removeEndMembers(image scatter,Image OvalIDs)
	{
	number HitsN =OvalIDs.ImageGetDimensionSize(1)
	for(number i=0;i<HitsN;i++)
		{
		number OvalID =OvalIDs.getPixel(0,i)

		Scatter.MSA112_removeOneEndMember(OvalID)
		}
	}
	







/*********************					*
**************************				**
Clustering ********************************
**************************				**
**********************					*/















void MSA112_unrestrictROIdisplay(roi targetROI)
{	
number roiID =targetROI.ROIGetID()
string NoRestrict ="GetROIFromID("+roiID+").ROISetRestrictToDisplay(0)"
try{ExecuteScriptString(NoRestrict);}
catch{break;} 
}

/******************************************************************
Class Ellipse selector
******************************************************************/

Class MSA_EllipseSelector : object
{
number    x1,y1,x2,y2,xp1,yp1,xp2,yp2
			//Dragon  along coord + perp coord
number xe1,ye1,xe2,ye2,eccent
			//Elipse long axis coord + eccentricity 
	
number width,height
number Vortexes,gap
image coordX,coordY,Mask //X,Y coordinates and current Mask
number scatterID,PerpID,AlongID,ellipseID

//void SetFrameName(object self,string NameP) Name = NameP
//string GetFrameName(object self) return Name

				
void DecartAlpha(object self,number x1,number y1,number x2,number y2,number &L,number &alpha)
{
L = sqrt((x2-x1)**2+(y2-y1)**2)
alpha =acos((x2-x1)/L)
if((y2-y1)>0) alpha *=(-1)
}


void FindPerpDragon(object self,number x1,number y1,number x2,number y2,number e,number g,number &xp1,number &yp1,number &xp2,number &yp2)
{
number L,alpha 
self.DecartAlpha(x1,y1,x2,y2,L,alpha)

number x0 =(x2+x1)/2
number y0 =(y2+y1)/2
xp1 = x0 + (L*e/2+g)*sin(alpha)
xp2 = x0 - (L*e/2+g)*sin(alpha)
yp1 = y0 + (L*e/2+g)*cos(alpha)
yp2 = y0 - (L*e/2+g)*cos(alpha)
}


void ExtractEccentFromDragon(object self,number x1,number y1,number x2,number y2,number &e,number g,number xp1,number yp1,number xp2,number yp2)
{
number Length,alpha 
self.DecartAlpha(x1,y1,x2,y2,Length,alpha)

number x0 =(x2+x1)/2
number y0 =(y2+y1)/2
number d1= sqrt((xp1-x0)**2+(yp1-y0)**2)-g
number d2= sqrt((xp2-x0)**2+(yp2-y0)**2)-g
e = (d1+d2)/Length
if(e>1) e=1
}


void ExtractEllipseFromDragon(object self,number x1,number y1,number x2,number y2,number g,number &xe1,number &ye1,number &xe2,number &ye2)
{
number Length = sqrt((x2-x1)**2+(y2-y1)**2)
number f=gap/Length

if(Length>2*gap){xe1 = x1+f*(x2-x1);xe2 = x2-f*(x2-x1);}
if(Length>2*gap){ye1 = y1+f*(y2-y1);ye2 = y2-f*(y2-y1);}
else {xe1 =0;xe2 = 0;ye1 =0;ye2 = 0;}
}

	
ROI DrawEllipseROI(object self,ROI ellipse,number x1,number y1,number x2,number y2,number e,number w,number h,number Vortexes)
{
number x,y,x0,y0
number L = sqrt((x2-x1)**2+(y2-y1)**2)	
if(Vortexes>L/4) Vortexes = round(L/4)
number alpha =acos((x2-x1)/L)
if((y2-y1)>0) alpha *=(-1)
//result("\nL "+L+" alpha "+alpha+" "+cos(alpha)+" "+sin(alpha))
	x = x1
	y = y1 
//result("\n"+0+" x "+x+" y "+y)
	ellipse.ROIAddVertex(x,y) //if(x>=0&&y>=0&&x<w&&y<h) 

for(number i=1;i<Vortexes;i++)
	{
	x0 = L/((Vortexes/i)**2+1)
	y0 = -e*Vortexes/i*x0 
	x = x1+cos(alpha)*x0+sin(alpha)*y0
	y = y1-sin(alpha)*x0+cos(alpha)*y0
//result("\n"+i+" x "+x+" y "+y)
	ellipse.ROIAddVertex(x,y)
	}

for(number i=Vortexes;i>=1;i--)
	{
	x0 = L/((i/Vortexes)**2+1)
	y0 = -i*e/Vortexes*x0 
	x = x1+cos(alpha)*x0+sin(alpha)*y0
	y = y1-sin(alpha)*x0+cos(alpha)*y0
//result("\n"+i+" x "+x+" y "+y)
	ellipse.ROIAddVertex(x,y)	
	}
for(number i=0;i<Vortexes;i++)
	{
	x0 = L/((i/Vortexes)**2+1)
	y0 = i*e/Vortexes*x0 
	x = x1+cos(alpha)*x0+sin(alpha)*y0
	y = y1-sin(alpha)*x0+cos(alpha)*y0
//result("\n"+i+" x "x+" y "+y)
	ellipse.ROIAddVertex(x,y)	
	}
	
for(number i=Vortexes;i>=1;i--)
	{	
	x0 = L/((Vortexes/i)**2+1)
	y0 = e*Vortexes/i*x0 
	x = x1+cos(alpha)*x0+sin(alpha)*y0
	y = y1-sin(alpha)*x0+cos(alpha)*y0
//result("\n"+i+" x "+x+" y "+y)	
	ellipse.ROIAddVertex(x,y)	
	}
	
ellipse.ROISetIsClosed(1)
ellipse.ROISetColor(1,0,0)
ellipse.MSA112_unrestrictROIdisplay()	
//ellipse.ROISetRestrictToDisplay(0)
ellipse.ROISetVolatile(0)
	//BUG
	//number opacity = 0.5
	//ellipse.ROISetDrawFilled(1)
	//ellipse.ROISetFillProperties(opacity,1,0,0)
	
return ellipse
}


void UpdateCluster(object self)
{
number e = sqrt(1-eccent**2)
number r = sqrt((xe2-xe1)**2+(ye2-ye1)**2)

//coordinates of f1,f2 ellipse focuses
number xf1 = (xe1+xe2)/2 + e*(xe2-xe1)/2
number xf2 = (xe1+xe2)/2 - e*(xe2-xe1)/2
number yf1 = (ye1+ye2)/2 + e*(ye2-ye1)/2
number yf2 = (ye1+ye2)/2 - e*(ye2-ye1)/2

Mask = 0
//sum of distances to both focuses less than ellipse long axis
mask = tert(sqrt((coordX-xf1)**2+(coordY-yf1)**2)+sqrt((coordX-xf2)**2+(coordY-yf2)**2)< r,1,0)
}


object init(object self,number Ellpoints,number gapP)
{
Vortexes = Ellpoints
gap = gapP

image scatter := GetFrontImage() 
getsize(scatter,width,height)
scatterID = scatter.ImageGetID()

TagGroup ScatterTags =MSA112_getScatterTags(scatter)	

coordX := MSA112_GetImageFromTags(ScatterTags,"Coordinates-X") 
coordY := MSA112_GetImageFromTags(ScatterTags,"Coordinates-Y") 
//coordX.showimage()
//coordY.showimage()
Mask = coordX.imageclone()
mask = 0

Mask.showimage()
mask.imageSetName("current mask")
mask.ImageGetImageDisplay(0).ImageDisplaySetContrastLimits(0,1)

	number minLength = gap - width/4 + 1
	if(minLength<0) minLength = 0
	x1=width/4 - minLength
	y1=height/2
	x2=3/4*width + minLength
	y2=height/2
	eccent =0.5		//default excentricity
		
roi Along = NewRoI()
Along.RoiSetLine(x1,y1,x2,y2)
roi Perp  = NewRoi()
self.FindPerpDragon(x1,y1,x2,y2,eccent,gap,xp1,yp1,xp2,yp2)
Perp.RoiSetLine(xp1,yp1,xp2,yp2)
Along.ROISetColor(1,0,0)
Perp.ROISetColor(1,0,0)
Along.ROISetVolatile(0)
Perp.ROISetVolatile(1)
Along.MSA112_unrestrictROIdisplay()
Perp.MSA112_unrestrictROIdisplay()

self.ExtractEllipseFromDragon(x1,y1,x2,y2,gap,xe1,ye1,xe2,ye2)
	//result("\n"+x1+" "+y1+" "+x2+" "+y2)
	//result("\n"+xp1+" "+yp1+" "+xp2+" "+yp2)
	//result("\n"+xe1+" "+ye1+" "+xe2+" "+ye2)

roi Ellipse = newROI()
Ellipse = self.DrawEllipseROI(Ellipse,xe1,ye1,xe2,ye2,eccent,width,height,Vortexes)
Ellipse.ROISetMoveable(0)	//not moved manually
Ellipse.ROISetVolatile(0)	

imageDisplay disp = scatterID.FindImageByID().ImageGetImageDisplay(0)
disp.ImageDisplayAddRoi(ellipse )
disp.ImageDisplayAddRoi(Perp)	
disp.ImageDisplayAddRoi(Along)
disp.ImageDisplaySetROISelected(Along,1)

PerpID =   Perp.RoiGetID()
AlongID = Along.RoiGetID()
EllipseID = Ellipse.RoiGetID()
	//Result("\nw "+w+" h "+h+" roiID "+PerpID+ " "+alongID+" "+EllipseID)

ConnectObject(PerpID,"changed","Connection-IDA",self,"ChangePerp")
ConnectObject(AlongID,"changed","Connection-IDB",self,"ChangeAlong")

self.UpdateCluster()

return self
}


~MSA_EllipseSelector(object self)
{
Result("\n ellipse destroyed")
}


void selectScatterPlot(object self)
	{
	scatterID.FindImageByID().ImageGetOrCreateImageDocument().ImageDocumentGetWindow().WindowSelect()
	}


void ChangePerp(object self, ROI Perp )
	{
	roi Along = AlongID.GetRoiFromID()
	Along.ROIGetLine(x1,y1,x2,y2)
	Perp.ROIGetLine(xp1,yp1,xp2,yp2)
	self.ExtractEccentFromDragon(x1,y1,x2,y2,eccent,gap,xp1,yp1,xp2,yp2) //extract eccent
//result("\neccent "+eccent)
	self.FindPerpDragon(x1,y1,x2,y2,eccent,gap,xp1,yp1,xp2,yp2) //move perp line after along line
	Perp.ROISetLine(xp1,yp1,xp2,yp2)
	
	imageDisplay disp = scatterID.FIndIMageByID().ImageGetImageDisplay(0)
	roi Ellipse = disp.ImageDisplayLookupROIByID(EllipseID)
	self.ExtractEllipseFromDragon(x1,y1,x2,y2,gap,xe1,ye1,xe2,ye2)
	if(Ellipse.ROIIsValid()) Ellipse.ROIClearVertices()
	else 
		{
		Ellipse = NewRoi();
		disp.ImageDisplayAddRoi(ellipse );
		EllipseID = Ellipse.RoiGetID();
		}
	Ellipse = self.DrawEllipseROI(Ellipse,xe1,ye1,xe2,ye2,eccent,width,height,Vortexes)
	disp.ImageDisplaySetROISelected(Ellipse,0)
	
	self.UpdateCluster()
	//Mask.showimage()
	}

	
void ChangeAlong(object self, ROI Along)
	{
	Along.ROIGetLine(x1,y1,x2,y2)
	roi Perp = PerpID.GetRoiFromID()
	self.FindPerpDragon(x1,y1,x2,y2 ,eccent,gap,xp1,yp1,xp2,yp2)
	Perp.ROISetLine(xp1,yp1,xp2,yp2)
	}


image getMask(object self) return mask


void getEllipseCoordinates(object self,number &x1P,number &y1P,number &x2P,number &y2P,number &eccentP)
	{
	x1P =x1
	y1P =y1	
	x2P =x2
	y2P =y2
	eccentP = eccent
	}
		

void StopSelection(object self)
	{
	image scatter := scatterID.FindImageByID()
	imageDisplay disp = scatter.ImageGetImageDisplay(0)
	roi Ellipse = disp.ImageDisplayLookupROIByID(EllipseID)
	if(Ellipse.ROIIsValid()) disp.ImageDisplayDeleteROI(Ellipse)
	roi Along = disp.ImageDisplayLookupROIByID(AlongID)
	if(Along.ROIIsValid()) disp.ImageDisplayDeleteROI(Along)
	roi Perp = disp.ImageDisplayLookupROIByID(PerpID)
	if(Perp.ROIIsValid()) disp.ImageDisplayDeleteROI(Perp)

	}

}



/***********************************************************************
Class Scatterplot for clustering
***********************************************************************/
	
class MSA_ClusteringScatter : MSA112_Frame
{
number ControlID
number ScatterID
number Comp1,Comp2

void init(object self, imagedocument scatterDoc,string Package,number Comp1P,number Comp2P, number ControlIDP)
	{
	ControlID = ControlIDP
	Comp1 = Comp1P
	Comp2 = Comp2P
	
	ScatterID = scatterDoc.ImageDocumentGetImage(0).ImageGetID()
	
	self.MSA112_initFrame(scatterDoc,Package,"Scatter "+Comp1+"-"+Comp2)

		//DocumentWindow Win = scatterDoc.ImageDocumentGetWindow()
		//number closeID = Win.WindowAddWindowListener(ControlID.GetScriptObjectFromID(),"window_closed:ScatterClosed")	
	}	

~MSA_ClusteringScatter(object self) 
	{
	try{
		object Control = GetScriptObjectFromID(ControlID)
		
		//remove Tags
		Control.RemoveScatterPlot(ScatterID)
		}
	catch{break;}	
	}
}


/***********************************************************************
Class ControlTags for clustering in scatter plots
***********************************************************************/

class MSA_ClusteringControl 
{
/*
ClusteringControl

						
ScatterPlots:	[0]:	Scatter ID 	//ID of scatterPLot image			
						Comp 1     	//first PCA component in scatter
						Comp 2	   	//second PCA component in scatter
						Listener ID	//display listener ID
						Members:	[0]:	Oval ID	//ID of the Oval showimng position of a member
									[1] .....

ControlTags.TagGroupOpenBrowserWindow(0)
									
*/


//Setting parameters
string MSAlabel		
number sparsity

number MSAID
number StarterID	//to call Starter 
TagGroup ControlTags
number scattN//current number of scatterPLots and EndMembers

image Ranges



/* *******  ControlTags ************* */

TagGroup GetControlTags(object self) return ControlTags

void CreateControlTags(object self)
	{	
	ControlTags = NewTagGroup()
	ControlTags.TagGroupSetTagAsTagGroup("ScatterPlots",NewTagList())

	//add it to Ranges
	TagGroup ImageTags = Ranges.ImageGetTagGroup()
	TagGroupSetTagAsTagGroup(ImageTags,MSA112_Label(),ControlTags)	
	}
			
		
image GetScatterPlot(object self,number scatt)
	{
	number scatterID = ControlTags.MSA112_getNumberfromList("ScatterPlots",scatt,"Scatter ID","L")
//result("\nscattID "+scatterID)
	try{	
		image scatter := getImageFromID(scatterID)
		return scatter
		}
	catch{break;}
	
	return NULL
	}

number FindScatterIndex(object self,image scatter) 
	{
	number scatterID = scatter.ImageGetID()
	for (number i = 0;i<scattN;i++)
		{
		if(scatterID==ControlTags.MSA112_getNumberfromList("ScatterPlots",i,"Scatter ID","L")) return i
		}
	return (-1)	
	}

number CheckIfScatterExist(object self,number first,number second) 
	{
	for (number i = 0;i<scattN;i++)
		{

		number Comp1 = ControlTags.MSA112_getNumberfromList("ScatterPlots",i,"Comp 1","S")
		number Comp2 = ControlTags.MSA112_getNumberfromList("ScatterPlots",i,"Comp 2","S")
//result("\nscatter "+i+" "+Comp1+"-"+Comp2+ " new "+first+"-"+second)		
		number exist = ((comp1==first)&&(comp2==second)) || ((comp1==second)&&(comp2==first))
//result(" exist "+exist)		
		if(exist)
			{
			OKdialog("Scatterplot of such components is already drawn !!") 
			return 0
			}
		}
		return 1
	}


void AddScatterPlot(object self,number first,number second) 
	{
	//add new scatter plot
//result("\nSP first "+first+" second "+second)	
	if(self.CheckIfScatterExist(first,second)) 
		{
//result("\nadd scatterplot "+first +"-"+second)
		image MSAimage:= getImageFromID(MSAID)	
		TagGroup MatrixTags = MSA112_GetMatrixTags(MSAimage,0)
		number CurrCluster =0
		
		//create scatterPLot	
		string Method = "PCA"		
		rgbimage scatter := MSA112_PlotXY(MSAimage,currCluster,Method,first,second,sparsity)
		scatter.showimage()

		//add listener for changing annotations
		imagedisplay ScatDisp = scatter.ImageGetImageDisplay(0)	
		number ListenerID = ScatDisp.ImageDisplayAddEventListener(self,"component_changed,component_removed:MemberChanged")
	
		//wrap into Frame112
		imageDocument Scatdoc = scatter.ImageGetOrCreateImageDocument()
		Scatdoc.ImageDocumentShow()
		
		DocumentWindow Win = Scatdoc.ImageDocumentGetWindow()
		number closeID = Win.WindowAddWindowListener(self,"window_closed:ScatterClosed")
		
		alloc(MSA_ClusteringScatter).init(Scatdoc,MSA112_MSApath()+":EndMembers",first,second,self.ScriptObjectGetID())
//result(" - ScatterPlot "+first+"-"+second)	
	
		
		//make scatter tags
		ControlTags.MSA112_setNumberinList("ScatterPlots",scattN,"Scatter ID","L",scatter.imageGetID())
		ControlTags.MSA112_setNumberinList("ScatterPlots",scattN,"Listener ID","L",closeID)
		ControlTags.MSA112_setNumberinList("ScatterPlots",scattN,"Comp 1","S",first)
		ControlTags.MSA112_setNumberinList("ScatterPlots",scattN,"Comp 2","S",second)
		
		scattN ++

ControlTags.TagGroupOpenBrowserWindow(0)		
		}
	}

void scatterClosed(object self, number e_fl, DocumentWindow Win)
    {
	result("\nclicked")
	}	

	
/* ******* initiation ************* */	
number init(object self, image MSAimage,number StarterIDP,number sparsityP)
	{
	MSAlabel = MSA112_Label()	
	sparsity = sparsityP

	StarterID = StarterIDP
	
	//enable Tuning Box in starter
	StarterID.getScriptObjectFromID().EnablingButtons(1)
	
	MSAID = MSAimage.ImageGetID()
	
	//get MSA tags
	string DataName = MSAImage.getName()
	TagGroup MatrixTags = MSA112_GetMatrixTags(MSAimage,0)
		
	
	//read info from tags
	string ReadMethod = "PCA"
	number currCluster =0
	number m = MSA112_ReadNumberComponents(MatrixTags,currCluster,ReadMethod)
	if(!MSA112_CheckComponents(m)) throw("")	
	
	//make Ranges
	Ranges := realimage("PCA Fractions",4,m,2)
//ranges.showimage()
	//Make Control tags
	self.CreateControlTags()
	
	result("\nClusteringControl "+self.ScriptObjectGetID()+" initialized")	
	
	
	return self.ScriptObjectGetID()
	}
	

~MSA_ClusteringControl(object self) 
	{
	//disable Tuning Box in starter
	object Starter = StarterID.getScriptObjectFromID()
	if(Starter.ScriptObjectIsValid()) Starter.EnablingButtons(0)
//result("\nMSA_EndMembersTableau destroyed")
	
result("\nClusteringControl "+ self.ScriptObjectGetID()+" destroyed")	
	}	
}