pow <- function(x, y){

	result <- x^y
	

}
RGBtoLAB <- function( Rin,  Gin,  Bin) 
{	
	# useful information at
	#http:#www.brucelindbloom.com
	  

	eps = 216.0/24389.0
	k = 24389.0/27.0
	   
	Xr = 0.964221 # reference white D50
	Yr = 1.0
	Zr = 0.825211
	
	# RGB to XYZ
	r = Rin/255.0 #R 0..1
	g = Gin/255.0 #G 0..1
	b = Bin/255.0 #B 0..1
	
	# assuming sRGB (D65)
	if (r <= 0.04045){
		r = r/12.0
	} else{
		r = pow((r+0.055)/1.055,2.4)
	}
	
	if (g <= 0.04045){
		g = g/12.0
	} else {
		g = pow((g+0.055)/1.055,2.4)
	}
	
	if (b <= 0.04045){
		b = b/12.0
	} else {
		b = pow((b+0.055)/1.055,2.4)
	}
	
	X =  0.436052025*r     + 0.385081593*g + 0.143087414 *b
	Y =  0.222491598*r     + 0.71688606 *g + 0.060621486 *b
	Z =  0.013929122*r     + 0.097097002*g + 0.71418547  *b
						
	# XYZ to Lab
	xr = X/Xr
	yr = Y/Yr
	zr = Z/Zr


			
	if ( xr > eps ){
		fx =  pow(xr, 1/3.0)
	} else{
		fx = ((k * xr + 16.0) / 116.0)
	}
	 
	if ( yr > eps ) {
		fy =  pow(yr, 1/3.0)
	} else {
		fy = ((k * yr + 16.0) / 116.0)
	}
	
	if ( zr > eps ) {
		fz =  pow(zr, 1/3.0)
	} else {
		fz = ((k * zr + 16.0) / 116.0)
	}
	
	Ls = ( 116 * fy ) - 16
	as = 500*(fx-fy)
	bs = 200*(fy-fz)
	
	# round up
	lab<-c(0, 0, 0)
	lab[0] = ceiling(Ls )
	lab[1] = ceiling (as )
	lab[2] = ceiling(bs )    
	
	cat(Rin, Gin, Bin, " converts to ", lab, "\n")
	
	return(lab)
} 

ColorDistance <- function(lab1, lab2)
	{
		 L1 = lab1[1]
		 a1 = lab1[2]
		 b1 = lab1[3]

		 L2 = lab2[1]
		 a2 = lab2[2]
		 b2 = lab2[3]
						
		Lbarprime = (L1 + L2)/2
		C1 = sqrt(pow(a1,2) + pow(b1,2))
		C2 = sqrt(pow(a2,2) + pow(b2,2))
		Cbar = (C1 + C2)/2
		
		G = (1- sqrt(pow(Cbar, 7)/(pow(Cbar,7)+pow(25,7)) ) )/2
		
		a1prime = a1*(1+G)
		a2prime = a2*(1+G)
	
		C1prime = sqrt(pow(a1prime,2) + pow(b1,2) )
		C2prime = sqrt(pow(a2prime,2) + pow(b2,2) )
		Cbarprime = (C1prime + C2prime)/2
		
		h1prime = 180*atan2(b1, a1prime)/pi
		if(h1prime<0){
			h1prime = h1prime + 360
		}
		h2prime = 180*atan2(b2, a2prime)/pi
		if(h2prime<0){
			h2prime = h2prime + 360		
			}
		
		if(abs(h1prime-h2prime)>180){
			Hbarprime = (h1prime +h2prime+360)/2
			} else{
			Hbarprime = (h1prime+h2prime)/2
		}
		
		T = 1 - 0.17*cos(pi*(Hbarprime-30)/180) + 0.24*cos(pi*(2*Hbarprime)/180) +0.32*cos(pi*(3*Hbarprime+6)/180) - 0.2*cos(pi*(4*Hbarprime-63)/180)
				
		
		if(abs(h2prime-h1prime)<=180){
			dhprime = h2prime-h1prime
			} else{ 
			if (h2prime<=h1prime){
				dhprime = h2prime-h1prime +360
				} else{
				dhprime = h2prime-h1prime - 360
				}
		}

		dLprime = L2 - L1
		dCprime = C2prime - C1prime
		dHprime  = 2* sqrt(C1prime*C2prime)*sin(dhprime/2*pi/180)
		Sl = 1 + 0.015*pow(Lbarprime-50,2)/sqrt(20+pow(Lbarprime-50,2))
		Sc = 1 + 0.045*Cbarprime
		Sh = 1 + 0.015*Cbarprime*T
		
		dtheta = 30*exp(-pow( (Hbarprime-275)/25,2))
		
		Rc = sqrt(pow(Cbarprime,7)/(pow(Cbarprime,7) + pow(25, 7)))
		
		Rt = -2*Rc * sin(2 * dtheta * pi/180) # from Lindbloom
		Kl = 1
		Kc = 1
		Kh = 1
		
		dE = sqrt( pow(dLprime/(Kl*Sl), 2) +pow(dCprime/(Kc*Sc), 2) + pow(dHprime/(Kh*Sh), 2)+Rt*dCprime/(Kc*Sc)*dHprime/(Kh*Sh) )

		return(dE)
		}


# Compute color similarity and contrast for given RGB values

ColorBackground<-c(255, 255, 255) # White



Names <- c("fun", "csv", "link", "probability", "map", "stan", "bayesian", "statistics", "distribution", "ulam", "quap", "rethinking", "prior", "posterior")

numPixels <- c(3520, 2772, 3080, 12291, 4410, 3783, 9522, 8446, 10280, 4830, 4872, 11165, 5151, 9945)

# First word cloud color assignments

ColorTarget<-c(99, 185, 192)  # Cyan
ColorDistractor1 <-c(224, 168, 83)  # Orange
ColorDistractor2 <-c(85, 135, 77)  # Green
ColorDistractor3 <-c(209, 48, 53)  # Red

Colors<-list(ColorTarget, ColorDistractor1, ColorDistractor1, ColorTarget, ColorDistractor1, ColorDistractor3, ColorTarget, ColorDistractor3, ColorTarget, ColorTarget, ColorDistractor2, ColorDistractor2, ColorDistractor2, ColorDistractor3 )


# Second word cloud
ColorTarget<-c(112, 241, 187)  # Light green
ColorDistractor1 <-c(242, 242, 242)  # Light gray
ColorDistractor2 <-c(51, 28, 170)  # purple

Colors<-list(ColorTarget, ColorTarget, ColorDistractor1, ColorTarget, ColorDistractor1, ColorDistractor2, ColorTarget, ColorDistractor2, ColorTarget, ColorTarget, ColorDistractor1, ColorDistractor1, ColorDistractor1, ColorDistractor2 )

# Third word cloud
ColorTarget<-c(225, 131, 121)  # Reddish
ColorDistractor1 <-c(224, 222, 86 )  # yellow
ColorDistractor2 <-c(114, 200, 229)  # light blue
ColorDistractor3 <-c(40, 80, 154)  # dark blue

Colors<-list(ColorTarget, ColorDistractor1, ColorDistractor1, ColorTarget, ColorDistractor1, ColorDistractor3, ColorTarget, ColorDistractor3, ColorTarget, ColorTarget, ColorDistractor2, ColorDistractor2, ColorDistractor2, ColorDistractor3 )

# Convert colors to Lab space
labColorsL<-c()
labColorsA<-c()
labColorsB<-c()
for(rgb in Colors){
	lab<-RGBtoLAB(rgb[1], rgb[2], rgb[3])
	labColorsL<-c(labColorsL, lab[1])
	labColorsA<-c(labColorsA, lab[2])
	labColorsB<-c(labColorsB, lab[3])

}

cat("Got LAB colors\n")

# Compute color distances (weighted by size)
sum=0
maxDistance= -1
for(item in c(1:length(Names))){
	thisDistance = ColorDistance(c(labColorsL[1], labColorsA[1], labColorsB[1]), c(labColorsL[item], labColorsA[item], labColorsB[item]))
	sum = sum + numPixels[item]*thisDistance

	if(thisDistance>maxDistance){maxDistance=thisDistance}

}

cat("Color distance between target and cloud is ", sum/maxDistance, "\n")

# Compute contrast with white
maxWdistance= -1
targetContrast= -1
for(item in c(1:length(Names))){
	thisDistance = ColorDistance(c(1, 0, 0), c(labColorsL[item], labColorsA[item], labColorsB[item]))

	if(item==1){targetContrast=thisDistance}
	if(thisDistance>maxWdistance){maxWdistance=thisDistance}

}
cat("Contrast value is ", targetContrast/maxWdistance, "\n")