// Cosine simialrity and diffraction similarity for 3D stack diffractions
// Input(3D): img3D
// Output   : Cosine_similarity, diffraction similarity 
// by Koji KIMOTO

// Cosine similarity
image img_CosSim(image img3D)
{
string nameX = "Cosine similarity_" + getname(img3D)  
number byte=4
number sizeX, sizeY, sizeZ			// Size  Get3DSize(img, sizeX, sizeY, sizeZ) 
	sizeX =ImageGetDimensionSize(img3D, 0)
	sizeY =ImageGetDimensionSize(img3D, 1)
	sizeZ =ImageGetDimensionSize(img3D, 2)
image tempA :=realimage("temp",byte, sizeX, sizeY)				
image tempB :=realimage("temp",byte, sizeX, sizeY)				
image imgX  :=realimage(nameX, byte, sizeZ, sizeZ)
 imgX= -9999
number i,j 
number numX  
 for(i=0;i<sizeZ;i++)
 {
	tempA = img3D[icol, irow, i]
	for(j=i;j<sizeZ;j++)
		{
		tempB = img3D[icol, irow, j]
			number Amp1 = sqrt(sum(tempA**2))
			number Amp2 = sqrt(sum(tempB**2))
			number ab = sum(tempA * tempB)
			numX = ab/(Amp1 * Amp2)
		SetPixel(imgX, i,j,numX)
	 // SetPixel(imgX, j,i,numX)
		}
 }
SetUnitString(imgX, "Slice")
  {
	TagGroup SourceTags =ImageGetTagGroup(img3D)
	TagGroup TargetTags =ImageGetTagGroup(imgX )
	TagGroupCopyTagsFrom(TargetTags, SourceTags)
  }
  return imgX
}

// R-phi transformation
image R_Phi(image img)
{
number sizeX,sizeY
GetSize(img, sizeX,sizeY)
number halfminor = min(sizeX, sizeY)/2
number centerX = sizeX/2
number centerY = sizeY/2 
number sampling = halfminor*4
image band
Band := RealImage( "R-Phi projection like by Mercator", 4, halfMinor, sampling )
number k = 2 * pi() / sampling
Band = warp(img, icol*sin(irow*k) + centerx, icol*cos(irow*k) + centery )
return Band
}

// Diffraction similarity
image img_difSim(image img3D)
{
string nameCC= "Diffraction similarity_" + getname(img3D)
number byte=4
number sizeX, sizeY, sizeZ			 
	sizeX =ImageGetDimensionSize(img3D, 0)
	sizeY =ImageGetDimensionSize(img3D, 1)
	sizeZ =ImageGetDimensionSize(img3D, 2)
image tempA :=realimage("temp",byte, sizeX, sizeY)				
image tempB :=realimage("temp",byte, sizeX, sizeY)				
//image tempC :=realimage("temp",byte, sizeX, sizeY)				
image imgCC :=realimage(nameCC,byte, sizeZ, sizeZ)
 imgCC=-9999 
number i,j 
number CC  
number CCmaxX, CCmaxY, CCmax
 for(i=0;i<sizeZ;i++)
 {
	tempA = img3D[icol, irow, i]
	image tempA_RPhi := R_Phi(tempA)
	//number sizeR, sizePhi
	//getsize(tempA_RPhi, sizeR, sizePhi)
	for(j=i;j<sizeZ;j++)
		{
		tempB = img3D[icol, irow, j]
		image tempB_RPhi := R_Phi(tempB)
	
		image tempC := crosscorrelation(tempA_RPhi,tempB_RPhi)
		tempC = tert(icol==iwidth/2, tempC, 0)
		CCmax = max(tempC, CCmaxX, CCmaxY)
		
		SetPixel(imgCC, i,j,CCmax)// or SetPixel(imgCC, i,j,CC)
		//SetPixel(imgCC, i,j,CCmax)
	
		}
 }
 
 setunitstring(imgCC, "Slice")
	{
	TagGroup SourceTags =ImageGetTagGroup(img3D)
	TagGroup TargetTags =ImageGetTagGroup(imgCC)
	TagGroupCopyTagsFrom(TargetTags, SourceTags)
	}

  return imgCC
}

// Psedudo color: (white) Blue - Gray - Red
void SetCustomLUT(image img)
{
 imageDisplay disp = img.ImageGetImageDisplay( 0 )
 // Contrast limit (display range)
 number low = -1 //min(img)
 number high = 1  //max(img)
    disp.ImageDisplaySetDoAutoSurvey( 0 ) 
	disp.ImageDisplaySetContrastLimits( low, high )
	ImageDisplaySetCaptionOn( disp, 1 )   
 // Brightness and contrast
	number bright   = 0.5			// 50%
	number contrast = 0.5			// 50 %
	disp.ImageDisplaySetContrastParameters( bright, contrast )
 // Apply ad-hoc myLUT
	image myLUT:=  RGBimage("(white<-1)-Blue-gray-red",4, 256, 1)
	//showimage(myLUT)
	number g = 30
	myLUT[0,0,1,128]   = RGB((255-g)*(icol/128), (255-g)*(icol/128), 255-g*(icol/128) )
	myLUT[0,128,1,256] = RGB((255)-g*((1-icol)/128), (255-g)*(1-(icol/128)), (255-g)*(1-(icol/128)) )
	myLUT[0,0,1,1]     = RGB(255,255,255)
 disp.ImageDisplaySetInputColorTable( myLUT ) 
}

void AlignWindowX(image img1, image img2)
{
 number offsetX=22
	number winsizeX, winsizeY, winposX, winposY
	getwindowsize(img1,winsizeX,winsizeY)
	getwindowposition(img1, winposX, winposY)
	Setwindowsize(img2,winsizeX, winsizeY)
	Setwindowposition(img2, winposX+winsizeX+OffsetX, winposY)
}


// main

image img3D:=getfrontimage()

// Cosine similarity 
	image img1 := img_CosSim(img3D)
	showimage(img1)
	SetCustomLut(img1)
	AlignWindowX(img3D, img1)
	
// Diffraciton similarity 
	image img2 := img_DifSim(img3D)
	showimage(img2)
	SetCustomLut(img2)
	AlignWindowX(img1, img2)