function [controlPoints, contour] = gibbsSampler(U,controlPoints,contour,h,Temp,iCP,nSamples)

% Subfunction of the Simulated Annealing program.
% Updates a control point by picking a new value at random
% from the local conditional probability distribution.


% Initialize auxiliary variables
nOrder = 3;
dev = 10;

% Generate sample set
n = floor(nSamples/2);
x = -dev:(dev/n):dev;
L = length(x);
newVals = x + controlPoints(iCP);

% Compute the energy of the contour given each new value
% of the control point
for i = 1:L,
    H(i) = computeH(U, controlPoints, h, iCP, newVals(i));
end

% Convert the energy array to probability weights
H = H - min(H);
H = exp( -H / Temp );
H = H ./ sum(H);

% Choose one of the candidate locations based on the
% probability of each one of them
rVal = rand;
k = 1;
[sortedH, index] = sort(H);
probArray = cumsum(sortedH);
while ( rVal > probArray(k) )
    k = k + 1;
end

% Update coefficient array
controlPoints(iCP) = newVals(index(k));

% Update contours array
nControlPoints = length(controlPoints);
nCBegin = max(1, iCP - nOrder);
nCEnd = min(iCP + nOrder, nControlPoints);
subCP = controlPoints(nCBegin:nCEnd);
subContour = formBspline(subCP,h);
ntBegin = ( nCBegin - 1 ) * h + 1;
ntEnd = ntBegin + length(subContour) - 1;
contour(ntBegin:ntEnd) = subContour;