function [y, m, b]=medmedline(x)
%function [y,m,b]=medmedline(x)
%  This function implements the median-median line approximation
%to the n x 2 data points in x.  The program will put x in the
%right format, so it can be input as n x 2 or 2 x n.  The program
%will also sort the data, so the domain need not be sorted already.
%
% Output:  y = n x 3 array, first column is domain, second
%              column are the values on the median-median line
%              and the third column is the set of residuals.
%          m = slope of the median-median line
%          b = intercept of the median-median line
%
%This method is most appropriate when the data has a few large
% errors, and the errors are in both domain and range.

%Get the matrix x in the desired format: Tall and ordered
[r,c]=size(x);
if r<c
    x=x';
end
[vals,idx]=sort(x(:,1));
x=x(idx,:);
[n,c]=size(x);  %n is the number of data points

%Separate into three groups
clear r c idx vals
idxgps=floor(n/3);
remainder=rem(n,3);
idx=idxgps*[1,1,1];  %holds the number of pts in each group

switch remainder
    case(1)
        idx(2)=idx(2)+1;
    case(2)
        idx(1)=idx(1)+1;
        idx(3)=idx(3)+1;
end
%Error check
if sum(idx)~=n
  fprintf('Error in dividing groups\n');
  y=0; m=0; b=0;
  return
end
numidx=cumsum(idx);
group1=x(1:numidx(1),:);
group2=x(numidx(1)+1:numidx(2),:);
group3=x(numidx(2)+1:end,:);

a=median(group1);
b=median(group2);
c=median(group3);

%Get the equation of the line
m=(c(2)-a(2))/(c(1)-a(1));  %slope
b1=-m*a(1)+a(2);            %intercept before shift

  %Determine the shift
  y2=m*b(1)+b1;
  shift1=(b(2)-y2)/3;
b=b1+shift1;  

%Output the median median line data in y
y1=m*x(:,1)+b;
resids=x(:,2)-y1;
y=[x(:,1), y1, resids];