function H = clusterPlot(X, Y, data, clusters, colmapfun)

hold on; 

% Plot data
H = imagesc(X, Y, data);
lims = get(gca, 'clim');

% Plot cluster
[meshX, meshY] = meshgrid(X, Y);

numCluster = length(clusters);

if nargin < 5
    colmapfun = @hsv;
end

colMap = colmapfun(numCluster);

for iCluster = 1:numCluster
    contour(meshX, meshY, clusters{iCluster}, 1, 'color', colMap(iCluster, :), 'linewidth', 1);
end

set(gca, 'clim', lims);

end


