function [im, LH,PH,LD, PD,ang] = rotateImage(im, LH,PH,LD, PD)
    v = [1 0];
    u = [PH(1) - LH(1), PH(2)-LH(2)];
    ang = angleBVect(u,v);
    if(PH(2) <= LH(2))
        ang = -ang;
    end
    %Orange Points rotation
    
    R = [cosd(ang) -sind(ang); sind(ang) cosd(ang)];
    [H, W,~] = size(im);
    center1X = W/2;
    center1Y = H/2;
    LH(1) = LH(1) - center1X;
    LH(2) = center1Y - LH(2);
   

    LD(1) = LD(1) - center1X;
    LD(2) = center1Y - LD(2);

    PH(1) = PH(1) - center1X;
    PH(2) = center1Y - PH(2);

    PD(1) = PD(1) - center1X;
    PD(2) = center1Y - PD(2);

    LH = rot90(LH,-1);
    LD = rot90(LD,-1);
    PH = rot90(PH,-1);
    PD = rot90(PD,-1);

    LH2 = R*LH;
    LD2 = R*LD;
    PH2 = R*PH;
    PD2 = R*PD;

    LH(1) = LH2(1) + center1X;
    LH(2) = center1Y-LH2(2);

    LD(1) = LD2(1) + center1X;
    LD(2) = center1Y-LD2(2);

    PH(1) = PH2(1) + center1X;
    PH(2) = center1Y-PH2(2);

    PD(1) = PD2(1) + center1X;
    PD(2) = center1Y-PD2(2);

    LH = rot90(LH);
    LD = rot90(LD);
    PH = rot90(PH);
    PD = rot90(PD);

    im = imrotate(im,ang);

    [H, W,~] = size(im);
    center2X = W/2;
    center2Y = H/2;

    dx = center2X-center1X;
    dy = center2Y - center1Y;

    LH(1) = LH(1)+dx;
    LH(2) = LH(2)+dy;

    LD(1) = LD(1)+dx;
    LD(2) = LD(2)+dy;

    PH(1) = PH(1)+dx;
    PH(2) = PH(2)+dy;

    PD(1) = PD(1)+dx;
    PD(2) = PD(2)+dy;

end