function lines = proxSample(pic)

	debug = false;

	[xs, ys] = extractData(pic);
	data = [xs' ys'];
	lines = []
	
	
	%[Y,I]=sort(data(1,:));
	%data=data(:,I); %use the column indices from sort() to sort all columns of A.
	
	%scatter(data(:,1), data(:,2))
	xMin = min(xs);
	xMax = max(xs);
	yMax = max(ys);
	yMin = min(ys);
	numPoints = length(xs)
	
	gridSize = 10
	numRows= ceil(xMax/gridSize) +1
	numCols= ceil(yMax/gridSize) +1
	bins = cell( numRows, numCols);
	
	%Fill the bins
	for i=1:numPoints
		xIndex = ceil(xs(i)/gridSize);
		yIndex = ceil(ys(i)/gridSize);
		bins{xIndex, yIndex} = [bins{xIndex, yIndex}, i];
	end
	
	
	
	radius = 10;
	width = 1;
	threshold = 250;
	
	inliers = [];
	tracks = cell(0);
	numTracks = 0;
	
	for rowIn=1:numRows
		for colIn=1:numCols
			rowIn, colIn
			points = bins{rowIn,colIn};
			if length(points) == 0
				continue
			end
			
			i = sampleBin(bins, rowIn, colIn);
			x1 = xs(i);
			y1 = ys(i);
			%neighbors = arrayfun(@(idx) norm(data(idx,:) - [x1,y1]), i:size(data,1)) <= radius;
			%temp = data(i:end,:);
			%neighbors = temp(neighbors,:);
		
			xIndex = ceil(x1/gridSize);
			yIndex = ceil(y1/gridSize);
		
			neighbors = [];
			neighbors = [neighbors, sampleBin(bins, xIndex, yIndex)];
			neighbors = [neighbors, sampleBin(bins, xIndex+1, yIndex)];
			neighbors = [neighbors, sampleBin(bins, xIndex, yIndex+1)];
			neighbors = [neighbors, sampleBin(bins, max(1, xIndex-1), yIndex+1)];
			neighbors = [neighbors, sampleBin(bins, xIndex+1, yIndex+1)];
		
			numNeighbors = length(neighbors);
			for j=1:numNeighbors
				index = neighbors(j);
				x2 = xs(index);
				y2 = ys(index);
			
				%dis = dist(x1,y1, x2, y2)
				m = (y1 - y2)/(x1 - x2);		
				if abs(m) < 1,
					continue
				end
			
				yInt = y1 - m*x1;
			
				[track, indices] = getPoints(xs, ys, m, yInt, width);
				numInliers = length(indices);
				inliers = [inliers, numInliers];
			
				if numInliers >= threshold,
					numTracks = numTracks +1;
					tracks{numTracks,1} = track;
					tracks{numTracks,2} = [m, yInt];
					
				end
			end
		end
	end
	
	for k=1:(size(tracks))(1),
		disp('Tracks(k):');
		disp(k)
		%tracks{k}
		track = tracks{k,1};
		xTrack = track(:,1);
		yTrack = track(:,2);
		scatter(xTrack, yTrack);
		hold on
		line = tracks{k,2};
		disp('Slope, b, length')
		m = line(1)
		b = line(2)
		length(track)
		xRange = findEndPoints(xMin, xMax, yMin, yMax, m, b);
		plot(xRange, m*xRange + b)
		
		lines = [lines; [m,b, xRange(1), xRange(end), yMin, yMax]]; 
		
	end
	disp('Number of tracks: ')
	disp(length(tracks))
	
	sort(inliers)
		
end

function index = sampleBin(bins, i, j)
	points = bins{i,j};
	if length(points) == 0,
		index = [];
	else
		index = points(randsample(length(points), 1));
	end
end

function xRange = findEndPoints(xMin, xMax, yMin, yMax, m, b)
	x1 = (yMin - b)/m;
	x2 = (yMax - b)/m;

	xLow = min(x1, x2);
	xHig = max(x1, x2);
	xRange = [max(xLow,xMin), min(xHig, xMax)];
end

	
function dis = dist(x1,y1, x2, y2)
	dis = norm([x1- x2, y1 - y2]);

end