/******************************************************************
function: MakeImageOfRightType
******************************************************************/
image MSA112_Make2DImageOfRightType(number type,string Name,number Width,number Height)
{
image out

if(type==14)out := binaryimage(Name,1,Width,Height)

if(type==2) out := realimage(Name,4,Width,Height)
if(type==12)out := realimage(Name,8,Width,Height)

if(type==3) out := compleximage(Name,8,Width,Height)
if(type==13)out := compleximage(Name,16,Width,Height)

if(type==9) out := integerimage(Name,1,1,Width,Height)
if(type==6) out := integerimage(Name,1,0,Width,Height)
if(type==1) out := integerimage(Name,2,1,Width,Height)
if(type==10)out := integerimage(Name,2,0,Width,Height)
if(type==7) out := integerimage(Name,4,1,Width,Height)
if(type==11)out := integerimage(Name,4,0,Width,Height)
if(type==39)out := integerimage(Name,8,1,Width,Height)
if(type==40)out := integerimage(Name,8,0,Width,Height)

return out
}

//same but with the reference-type image
image MSA112_Make2DImageOfRightType(image template,string Name,number Width,number Height)
{
number type =ImageGetDataType(template)

return MSA112_Make2DImageOfRightType(type,Name,Width,Height)
}


image  MSA112_Make3DImageOfRightType(number type,string Name,number Width,number Height,number Depth)
{
image out

if(type==14)out := binaryimage(Name,1,Width,Height,Depth)

if(type==2) out := realimage(Name,4,Width,Height,Depth)
if(type==12)out := realimage(Name,8,Width,Height,Depth)

if(type==3) out := compleximage(Name,8,Width,Height,Depth)
if(type==13)out := compleximage(Name,16,Width,Height,Depth)

if(type==9) out := integerimage(Name,1,1,Width,Height,Depth)
if(type==6) out := integerimage(Name,1,0,Width,Height,Depth)
if(type==1) out := integerimage(Name,2,1,Width,Height,Depth)
if(type==10)out := integerimage(Name,2,0,Width,Height,Depth)
if(type==7) out := integerimage(Name,4,1,Width,Height,Depth)
if(type==11)out := integerimage(Name,4,0,Width,Height,Depth)
if(type==39)out := integerimage(Name,8,1,Width,Height,Depth)
if(type==40)out := integerimage(Name,8,0,Width,Height,Depth)

return out
}

//same but with the reference-type image
image MSA112_Make3DImageOfRightType(image template,string Name,number Width,number Height,number Depth)
{
number type =ImageGetDataType(template)

return MSA112_Make3DImageOfRightType(type,Name,Width,Height,Depth)
}


	
/******************************************************************
function: Cliping for Gaussian Filtering
******************************************************************/	
number MSA112_CastleClip(number &sigma,number minTail)
	{
	number NumbCells = (sqrt(-log(minTail))*sigma).trunc()
	
	if(NumbCells>=4) {NumbCells=3; sigma = 4/sqrt(-log(minTail));} //max 3cells on vertical
	
	return NumbCells
	}

number MSA112_BishopClip(number sigma,number minTail)
	{
	//No upper limit needed because it should be done when CastleClip
	return (sqrt(-log(minTail))*sigma/sqrt(2)).trunc()
	}	

//		sigma >= 0.660	castle 1 bishop 0
//		sigma >= 0.932	castle 1 bishop 1
//		sigma >= 1.319	castle 2 bishop 1
//		sigma >= 1.864	castle 2 bishop 2
//		sigma >= 1.978	castle 3 bishop 2





/***********************************************************************
functions One ALS step
***********************************************************************/

/************ set lambda on all diagonals *********************/
image MSA112_MakeDiagonalMatrix(number k,number lambda)
	{
	image unity := realimage("",4,k,k)

	for(number i=0; i<k; i++) {unity.setPixel(i,i,lambda);}
	
	return unity
	}

/************ set numbers from array on  diagonals *********************/
image MSA112_MakeDiagonalMatrix(image lambda)
	{
	number k = lambda.ImageGetDimensionSize(0)
	image unity := realimage("",4,k,k)

	for(number i=0; i<k; i++) {unity.setPixel(i,i,lambda.getPixel(i,0));}
	
	return unity
	}
	
image MSA112_ALStepBase(image D,image T,image unity)
	{
	image M1 = (matrixMultiply(MatrixTranspose(T),T) +unity)
	
	image M2 = MatrixInverse(M1)
		
	image M3 = matrixMultiply(M2,MatrixTranspose(T))

//	number k,n,n1,m
//	M3.getsize(n,k)
//	D.getSize(m,n1)
//result(" k="+k+" n="+n+" n1="+n1+" m="+m)
	image PT = matrixMultiply(M3,D)

	return MatrixTranspose(PT)
	}
	
image MSA112_ALStepBase(image D,image T,number lambda)
	{
	number k =T.ImageGetDimensionSize(1)
	
	image unity :=MSA112_MakeDiagonalMatrix(k,lambda)

	return MSA112_ALStepBase(D,T,unity)
	}
	