function [ct,h] = imgct( imgfile, ctfile, width, levels, N )
%
% plots level sets (contour lines) of discrete values in 'ctfile'
% on top of image in 'imgfile'. contour levels are given by the 'levels' vector.
% to plot a single contour line, set, e.g., levels = [0 0] (default)
% contours are computed in the (Nx + 1) x (Ny + 1) regular grid which spams
% a domain whose length is determined by the x and y values in ctfile.
%

I = imread( imgfile );          % display image
[W,H] = size(I);
imshow( flipud(I) );
axis xy;

hold on;

A = load(ctfile);
A = unique(A, 'rows');          % remove duplicates and sort

if nargin < 3                   % delta function support width
  width = 3;
end

if nargin < 5
  Nx = sqrt(size(A,1));         % assume Nx = Ny
  Ny = Nx;
else
  Nx = N(1) + 1;
  Ny = N(2) + 1;
end

if nargin < 4
  levels = [0 0];
end

X = reshape( A(:,1), Ny, Nx );
Y = reshape( A(:,2), Ny, Nx );
Z = reshape( A(:,3), Ny, Nx );

% scaling
dx =  X(Ny,Nx) - X(1,1);
dy =  Y(Ny,Nx) - Y(1,1);

X = W * (X - X(1,1)) ./ dx;     % scale to fit in mesh box
Y = H * (Y - Y(1,1)) ./ dy;

%
% plot contour on top of image
%

maxval = 255;
levels = levels / maxval;
style = '-r';                   % solid, red contour lines

[ct,h] = contour( X, Y, Z, levels, style );

% if i want to plot levels of the heaviside instead...
%mu = 0.99;                      % confidence interval
%h = 1.0/(W - 1);                % mesh size
%alpha = log((1 + mu)/(1 - mu)) / (width * h);
%Z = 0.5 * ( 1.0 + tanh( alpha * Z ));  % the heaviside function
%[ct1,h1] = contour( X, Y, Z, [0.5 0.5], style );

hold on;
set( gca, 'position', [0 0 1 1] );
set( gcf, 'position', [20 650 W H] );

hold off;