function [A,bdry] = imgseg3(filename,support)
% function [A,bdry] = imgseg3(filename,support)
% 
% computes the matrix related to the image segmentation problem
% for the image fn
% it identifies all the bdry nodes in bdry
%
% support is the size: minimum is 1

sig = 0.1;

if (ischar(filename)),
  img = double(imread(filename));
else
  img = filename;
end

  
sz = size(img);
xsz = sz(1)
ysz = sz(2)

imgv = double(reshape(img,xsz*ysz,3));

ai = [];
aj = [];
av = [];


for xi = 0:support,
  for yi = 0:(support-xi),
    
    if ((xi ~= 0) || (yi ~= 0)),
      tm = zeros(xsz,ysz);
      tm1 = tm;
      tm2 = tm;
      tm1(1:end-xi, 1:end-yi) = 1;
      tm2(1+xi:end, 1+yi:end) = 1;
      
      ind1 = find(tm1(:));
      ind2 = find(tm2(:));
      difv = imgv(ind1,:) - imgv(ind2,:);
      dif = sum(difv.^2,2);
      
      ai = [ai;ind1];
      aj = [aj;ind2];
      av = [av;dif];
      
      if (xi*yi ~= 0)
	tm = zeros(xsz,ysz);
	tm1 = tm;
	tm2 = tm;
	tm1(1+xi:end, 1:end-yi) = 1;
	tm2(1:end-xi, 1+yi:end) = 1;
	
	ind1 = find(tm1(:));
	ind2 = find(tm2(:));
	difv = imgv(ind1,:) - imgv(ind2,:);
	dif = sum(difv.^2,2);
	
	ai = [ai;ind1];
	aj = [aj;ind2];
	av = [av;dif];
	
      end
    end
  end
end

av = av/max(av);
av = exp(-av/sig^2);


n = xsz * ysz;

A = sparse(ai,aj,av, n ,n);
A = A + A';

tm = zeros(xsz,ysz);
tm(1,:) = 1;
tm(xsz,:) = 1;
tm(:,1) = 1;
tm(:,ysz) = 1;
bdry = (tm(:) ==1);