function newModel = midToEdgeDistanceModel(oldPts,connectedPts,modelType)
%   Function generates a 2D kpar kper model from a 3D pointset.
%
%   oldPts - 3D pointset (from for example, a 3D AAMToolbox project or 3D
%   growth model).
%   connectedPts - is a csv file which contains correspondence for the vertex numbers for edge
%   and midvein points.
%   modelType - 'simple' or 'angles'. Simple model calculates Euclidean
%   distance from midvein points to corresponding edge pts (as designated by
%   the connectedPts csv file) and plots them at 90 degrees. The angles model does
%   the same but instead of plotting pts at 90 degree, plots them at a best
%   fit angle calculated from the original 3D model.
%   newModel - the generated 2D model.

if nargin<3
    modelType = 'simple';
end

%debug plot shows point plots for debugging
debugPlot = 0;

%% reshape the data into a useful format
data = reshape(oldPts,3,size(oldPts,1)/3)';

if(debugPlot)
    plot3(data(:,1),data(:,2),data(:,3),'rx');
    for j = 1:size(data,1)
        text(data(j,1),data(j,2),data(j,3),num2str(j))
    end
    axis equal
end

%% calculate the midvein distances
midVeinIdxs = connectedPts(:,2,:);
for i = 1:size(connectedPts,1)-1
    midVeinDistances(i,:) = sqrt(sum(abs(data(midVeinIdxs(i),:) - data(midVeinIdxs(i+1),:)).^2));
end

%% calculate the left edge to midvein distances
leftEdgeIdxs = connectedPts(:,1,:);
for i = 1:size(connectedPts,1)-1
    leftEdgeDistances(i,:) = sqrt(sum(abs(data(midVeinIdxs(i),:) - data(leftEdgeIdxs(i),:)).^2));
end

%% calculate the right edge to midvein distances
rightEdgeIdxs = connectedPts(:,3,:);
for i = 1:size(connectedPts,1)-1
    rightEdgeDistances(i,:) = sqrt(sum(abs(data(midVeinIdxs(i),:) - data(rightEdgeIdxs(i),:)).^2));
end


%% SIMPLE MODEL
%replot pts
if(modelType == 'simple')
    newModelDataM = [0,0];
    for i = 1:size(midVeinDistances,1)
        newModelDataM(i+1,2) = newModelDataM(i,2)+ midVeinDistances(i,:);
    end
    newModelDataL = [leftEdgeDistances,newModelDataM(1:end-1,2)];
    newModelDataR = [-rightEdgeDistances,newModelDataM(1:end-1,2)];
    
    newModel = [newModelDataL;newModelDataR;newModelDataM];
    
    if(debugPlot)
        figure;
        plot(newModel(:,1),newModel(:,2),'rx');
        for j = 1:size(newModel,1)
            text(newModel(j,1),newModel(j,2),num2str(j))
        end
        axis equal
    end
end

%% ANGLES MODEL
if(modelType == 'angles')
    for i = 1:size(midVeinIdxs,1)-1
        if(i == 1)
            firstPtFlag = 1;
        else
            firstPtFlag = 0;
        end
        %find the pts
        a = data(midVeinIdxs(i),:);
        b = data(midVeinIdxs(i+1),:);
        if(~firstPtFlag)
            d = data(midVeinIdxs(i-1),:);
        end
        c = data(rightEdgeIdxs(i),:);
        e = data(leftEdgeIdxs(i),:);
        %move vectors to origin
        a = a-a;
        b = b-a;
        c = c-a;
        e = e-a;
        if(~firstPtFlag)
            d = d-a;
        end
        
        %calculate angles between vectors
        if(~firstPtFlag)
            %dac = acos(atan2(norm(cross(d,c)),dot(d,c)))*180/pi;
            %ead = acos(atan2(norm(cross(e,d)),dot(e,d)))*180/pi;
            dac = acos(dot(d,c)/(norm(d)*norm(c)))*180/pi;
            ead = acos(dot(e,d)/(norm(e)*norm(d)))*180/pi;
            cab = acos(dot(c,b)/(norm(c)*norm(b)))*180/pi;
            bae = acos(dot(b,e)/(norm(b)*norm(e)))*180/pi;
        else
            dac = nan;
            ead = nan;
            cab = 90;
            bae = 90;
            %cab = acos(atan2(norm(cross(c,b)),dot(c,b)))*180/pi;
            %bae = acos(atan2(norm(cross(b,e)),dot(b,e)))*180/pi;
        end
        
        
        
        % normalise so dac + cab = 180 and ead + bae = 180
        if(~firstPtFlag)
            daccabNorm = (180 - dac - cab)/2;
            eadbaeNorm = (180 - ead - bae)/2;
        else
            daccabNorm = 0;
            eadbaeNorm = 0;
        end
        dacAngles(i) = dac + daccabNorm;
        cabAngles(i) = cab + daccabNorm;
        eadAngles(i) = ead + eadbaeNorm;
        baeAngles(i) = bae + eadbaeNorm;
    end
    
    newModelDataM = [0,0];
    for i = 1:size(midVeinDistances,1)
        newModelDataM(i+1,2) = newModelDataM(i,2)+ midVeinDistances(i,:);
    end
    
    
    for i = 1:size(rightEdgeDistances,1)
        if(cabAngles(i) < 90)
            x1 = newModelDataM(i,1) - cos((pi/2)-(cabAngles(i) * (pi/180))) * rightEdgeDistances(i);
            y1 = newModelDataM(i,2) - sin((pi/2)-(cabAngles(i) * (pi/180))) * rightEdgeDistances(i);
        elseif(cabAngles(i) > 90)
            y1 = newModelDataM(i,2) - cos(pi-(cabAngles(i) * (pi/180))) * rightEdgeDistances(i);
            x1 = newModelDataM(i,1) - sin(pi-(cabAngles(i) * (pi/180))) * rightEdgeDistances(i);
        else
            y1 = newModelDataM(i,2);
            x1 = newModelDataM(i,1) - rightEdgeDistances(i);
        end
        newModelDataR(i,:) = [x1,y1];
    end
    
    for i = 1:size(leftEdgeDistances,1)
        if(cabAngles(i) < 90)
            x1 = newModelDataM(i,1) + cos((pi/2)-(cabAngles(i) * (pi/180))) * leftEdgeDistances(i);
            y1 = newModelDataM(i,2) + sin((pi/2)-(cabAngles(i) * (pi/180))) * leftEdgeDistances(i);
        elseif(cabAngles(i) > 90)
            y1 = newModelDataM(i,2) + cos(pi-(cabAngles(i) * (pi/180))) * leftEdgeDistances(i);
            x1 = newModelDataM(i,1) + sin(pi-(cabAngles(i) * (pi/180))) * leftEdgeDistances(i);
        else
            y1 = newModelDataM(i,2);
            x1 = newModelDataM(i,1) + leftEdgeDistances(i);
        end
        newModelDataL(i,:) = [x1,y1];
    end
    
    newModel = [newModelDataL;newModelDataR;newModelDataM];
    
    if(debugPlot)
        figure;
        plot(newModelDataM(:,1),newModelDataM(:,2),'rx');
        hold on
        plot(newModelDataR(:,1),newModelDataR(:,2),'bx');
        plot(newModelDataL(:,1),newModelDataL(:,2),'bx');
        for j = 1:size(newModel,1)
            text(newModel(j,1),newModel(j,2),num2str(j))
        end
        axis equal
    end
    
    
    
end




end