what's the problem for my head pose estimate?

[xunkaixin]

this is my code to get the matrix 15*3, but the rank is always 13, and the resolution seems not correct, can you give me some advise?

%get_samples
[points, angle] = get_samples();

%regression
X=points\angle;

function [samples_points, samples_angle] = get_samples()
samples_points = [];
samples_angle = [];

model_3d = [ 6.775290,-0.730814,-0.012799;
            0.000000,-7.415691,4.070434;
            -6.775290,-0.730814,-0.012799;
            3.6973, 7.8784, 7.0262;
            -3.6973, 7.8784,7.0262;
            5.311432,5.485328,3.987654;
            3.550622,6.185143,5.712299;
            1.789930,5.393625,4.413414;
            3.530191,4.981603,4.937805;
            -1.789930,5.393625,4.413414;
            -3.550622,6.185143,5.712299;
            -5.311432,5.485328,3.987654;
            -3.530191,4.981603,4.937805;
            0.000000,0.348381,6.763430;
            1.930245,0.424351,5.914376;
            0.000000,0.000000,5.914376;
            -1.930245,0.424351,5.914376;
            2.774015,-2.080775,5.048531;
            -2.774015,-2.080775,5.048531;
            0.000000,-1.646444,6.704956;
            0.000000,-3.116408,6.097667 ];

        %get 7 points
        idx = 1+[5,7,9,11,15,17,18];
        model_3d = model_3d(idx,:);
        model_3d(:,1) = -model_3d(:,1);
        npts = size(idx,2);
        
   for irand=1:100
       shape = model_3d;
        %眼睛、鼻子、嘴巴随机变化
        eye_pos = randperm(2,1)-1;
        eye_dx = -1+2*rand(1,1);
        nose_dy = -0.4+1.6*rand(1,1);
        mouth_dx = -0.5+1.2*rand(1,1);
        depth_pos = randperm(3,1);
        dz = -0.5 + rand(1,1);
        if eye_pos>0
            shape(1,1) = shape(1,1) + eye_dx;
            shape(4,1) = shape(4,1) - eye_dx;
        else
            shape(2,1) = shape(2,1,1) + eye_dx;
            shape(3,1) = shape(3,1,1) - eye_dx;
        end
        shape(5,2) = shape(5,2) + nose_dy;
        shape(6,1) = shape(6,1) + mouth_dx;
        shape(7,1) = shape(7,1) - mouth_dx;
        
        if depth_pos==1
            for s=1:4
                shape(s,3) = shape(s,3) + dz;
            end
        elseif depth_pos==2
            shape(5,3) = shape(5,3) + dz;
        elseif depth_pos==3  
            shape(6,3) = shape(6,3) + dz;
            shape(7,3) = shape(7,3) + dz;
        end
        
        
        %随机角度
        for ang_rand=1:100
            ypr = randperm(3,1);
            if ypr==1
                yaw = -60 + 120*rand(1,1);
                pitch = 0;
                roll = 0;
            elseif ypr==2
                yaw = 0;
                pitch = -30 + 60*rand(1,1);
                roll = 0;
            elseif ypr==3
                yaw = 0;
                pitch = 0;
                roll = -30 + 60*rand(1,1);
            end
            
            %rotate
            roted_shape = [];
            for s = 1:size(shape,1)
                roted_pt = rotate(shape(s,:), yaw/180.0*pi, pitch/180.0*pi, roll/180.0*pi);
                roted_shape = [roted_shape;roted_pt];
            end

            %{
            figure(5);         
            hold on;
            for ii=1:npts
                plot(roted_shape(ii,1),roted_shape(ii,2),'*r');
                %pause;
            end
            hold off;
            %}
            
            %project
            pro_shape = roted_shape(:,1:2);

            %normalize
            avg_x = mean(pro_shape(:,1));
            avg_y = mean(pro_shape(:,2));
            dist_y = max(pro_shape(:,2))-min(pro_shape(:,2));
            nor_shape = [(pro_shape(:,1)-avg_x)/dist_y (pro_shape(:,2)-avg_y)/dist_y];
            nor_shape = [reshape(nor_shape,1,npts*2) 1];

            %add
            samples_points = [samples_points; nor_shape];
            samples_angle = [samples_angle;[yaw, pitch, roll]];
        end
  end

end

function [out] = rotate(M, yaw, pitch, roll)
% Calculate the rotation matrix
yaw_matrix = [ cos(yaw) 0 -sin(yaw);
0 1 0;
sin(yaw) 0 cos(yaw)];

pitch_matrix   = [1 0          0;
                0 cos(pitch) -sin(pitch);
                0 sin(pitch)  cos(pitch)];

roll_matrix  = [cos(roll) -sin(roll) 0;
                sin(roll)  cos(roll) 0;
                0         0        1];
rot_matrix = roll_matrix*pitch_matrix*yaw_matrix; 

out = rot_matrix*M';
out = out';

end