function color = Colorize(img,cue)
	npixels = size(img,1)*size(img,2);
	entries = ones(npixels*9,3);
	entries(:,2) = 50;
	entries(:,3) = 0;
	
	% Construct Hessian
	for y=1:size(img,1)
		for x=1:size(img,2)
			index = (y-1)*size(img,2) + (x-1);
			entries(index*9+1:index*9+9,:) = NeighborhoodInfo(img,x,y,entries(index*9+1:index*9+9,:));
		end
	end
	hessian = sparse(entries(:,1),entries(:,2),entries(:,3),npixels,npixels);
	hessian = (hessian')*hessian;
	
	% Linear coeficients are zero
	linear = zeros(npixels,1);
	
	% Construct constraint equations
	entries = ones(npixels,4);
	entries(:,2) = 50;
	entries(:,3:4) = 0;
	numConstraints = 0;
	for y=1:size(img,1)
		for x=1:size(img,2)
			if(cue(y,x,2) > 0.1)
				numConstraints = numConstraints+1;
				index = (y-1)*size(img,2) + x;
				%entries(index,3:4) = cue(y,x,1:2);	% hsv
				entries(numConstraints,3:4) = cue(y,x,2:3);	% yiq
				entries(numConstraints,1) = numConstraints;
				entries(numConstraints,2) = index;
			end
		end
	end
	constraint_matrix = sparse(entries(1:numConstraints,1),entries(1:numConstraints,2),1,numConstraints,npixels);
	
	options = optimset('TolPCG', 0.00001);
	u = quadprog(hessian,linear,[],[],constraint_matrix,entries(1:numConstraints,3),[],[],[],options);
	v = quadprog(hessian,linear,[],[],constraint_matrix,entries(1:numConstraints,4),[],[],[],options);
	
	color = zeros(size(img,1),size(img,2),3);
	%color(:,:,3) = img(:,:,3);	% hsv
	color(:,:,1) = img(:,:,1);	% yiq
	for y=1:size(img,1)
		for x=1:size(img,2)
			index = (y-1)*size(img,2) + x;
			%color(y,x,1) = u(index);	% hsv
			%color(y,x,2) = v(index);
			color(y,x,2) = u(index);	% yiq
			color(y,x,3) = v(index);
		end
	end
end

function w = NeighborhoodInfo(img,x,y,w)
	% Create neighborhood array
	N = zeros(9,3);
	nsize = 0;
	for yn = max(1,y-1):min(size(img,1),y+1)
		for xn = max(1,x-1):min(size(img,2),x+1)
			if(yn == y && xn == x)
				continue;
			end
			nsize = nsize + 1;
			%val = img(yn,xn,3);	% hsv
			val = img(yn,xn,1);		% yiq
			N(nsize,:) = [yn,xn,val];
		end
	end
	N = N(1:nsize,:);
	mu = mean(N(:,3));
	variance = var(N(:,3));
	
	% Calculate weights
	%val = img(y,x,3)-mu;	% hsv
	val = img(y,x,1)-mu;	% yiq
	for i=1:size(N,1)
		if(variance == 0)
			w(i,3) = 1;
		else
			val2 = N(i,3) - mu;
			w(i,3) = 1 + (1/variance)*(val*val2);
		end
		w(i,2) = (N(i,1)-1)*size(img,2) + N(i,2);
	end
	w(:,3) = -w(:,3)/sum(w(:,3));
	
	% Add in affinity for self-loop
	w(9,3) = 1;
	w(9,1:2) = (y-1)*size(img,2) + x;
	
	% Index of the current pixel
	w(1:size(N,1),1) = w(9,2);
end