function xp = rok_camera_project_points (X, camera, scaling)
    % xp = ROK_CAMERA_PROJECT_POINTS (X, camera, scale)
    %
    % Projects 3-D point into camera's image plane, using the camera model
    % used by Camera Calibration Toolbox for Matlab.
    %
    % Input:
    %  - X = 3xM vector of points' world coordinates
    %  - camera = camera parameters structure (note: it is assumed that the
    %    alpha (skew) is 0!)
    %  - scaling = optional scaling parameters, for projection on resized images
    %
    % Output:
    %  - xp = 2xM vector of coordinates in image plane
    %
    % Cpyright (C) 2011, Rok Mandeljc <rok.mandeljc@fe.uni-lj.si>
    
    if nargin < 3,
        scaling = 1;
    end
    
    % Copy parameters, for easier access
    R = camera.R;
    T = camera.T;
    f = camera.f;
    c = camera.c;
    k = camera.k;
    %alpha = camera.alpha;
    
    n = size(X, 2); % Number of points

    % Transform into camera's coordinate system
    X = R*X + repmat(T, [1 n]);

    % Prune the points that fall behind the camera (Z <= 0)
    idx = X(3,:) <= 0;
    X(:, idx) = NaN;

    % Projective transformation
    inv_Z = 1./X(3,:);
    x = (X(1:2,:) .* (ones(2,1) * inv_Z)) ;

    % Add distortion
    r2 = x(1,:).^2 + x(2,:).^2;
    r4 = r2.^2;

    % Radial distortion
    cdist = 1 + k(1) * r2 + k(2) * r4;
    xd1 = x .* (ones(2,1)*cdist);

    % Tangential distortion:
    a1 = 2*x(1,:).*x(2,:);
    a2 = r2 + 2*x(1,:).^2;
    a3 = r2 + 2*x(2,:).^2;

    delta_x = [ k(3)*a1 + k(4)*a2;
                k(3)*a3 + k(4)*a1 ];

    xd2 = xd1 + delta_x;

    % Pixel coordinates:
    if length(f) > 1
        xp = xd2 .* (f * ones(1,n))  +  c*ones(1,n);
    else
        xp = f * xd2 + c*ones(1,n);
    end
    
    % Scaling
    xp = scaling*xp;
end
