plot_subspaces4.m

一个用EM算法的源程序

% Code for rendering three dimensional results generated by the gpca
% algorithm.  Users should call plot_subspaces
function pl = plot_subspaces4(X, labels, groupBases, basisDimensions)

GRID_LINE_COUNT = 15;

% Generate a list, existList, of labels that have had points assigned to them.
highestLabel = max(labels);
existList = zeros(highestLabel, 1);  % existList will have at most highestLabel elements.
groupCount = 0;
for labelIndex = 1:highestLabel,
    if (sum(labels == labelIndex) > 0) % If the group has at least one element in it,
        groupCount = groupCount + 1;
        existList(groupCount) = labelIndex;
    end;
end;
existList = existList(1:groupCount);

markers = ['.ox+*svp^h<>']; markerCount = length(markers); markerIndex = 1;
colors = ['rbkmgcy']; colorCount = length(colors); colorIndex = 1;

% Iterate through every group that contains at least one element.
for groupIndex = 1:groupCount,
        
    basis = groupBases{existList(groupIndex)}(:,1:basisDimensions(existList(groupIndex))); 
    groupData = X(:, labels == existList(groupIndex));
    groupSampleCount = sum(labels == existList(groupIndex));
    center = mean(groupData,2);
    
    if (size(basis,2) == 2)
        % No lines yet!
        
        % Plot the one dimensional subspace (line)
%         x = basis(1); 
%         y = basis(2); 
%         z = basis(3);
%         
%         if(CENTER_DATA),
%             projected_data = (groupData' - ones(size(groupData , 2) , 1) * center') * basis;
%             upper_bound = max(projected_data) * basis + center;
%             lower_bound = min(projected_data) * basis + center;
%         else,
%             projected_data = (groupData') * basis;
%             upper_bound = max(projected_data) * basis;
%             lower_bound = min(projected_data) * basis;
%         end
% 
%         l = line([upper_bound(1) lower_bound(1)], [upper_bound(2) lower_bound(2)], [upper_bound(3) lower_bound(3)]);
%         set(l, 'Color', colors(colorIndex));
%         set(l, 'LineWidth', 3);
        
    elseif (size(basis,2) == 3)
        % Plot the plane as a surface object
        
        % Compute a set of basis vectors that align with the plot axes....
        % (hack)
%         row_d(1) = norm(basis(1,:));
%         row_d(2) = norm(basis(2,:));
%         row_d(3) = norm(basis(3,:));
%         [ignore sort_index] = sort(row_d);
%         sort_index = fliplr(sort_index);
%         basis(sort_index,:) = rref(basis(sort_index,:)')';
        
        % Compute Bounds for the data.
        
            homogeneousValue = groupData(end,1)
            homogeneousProjectedData = basis*basis'*groupData;
            rescaledData = homogeneousProjectedData * diag(sparse(homogeneousValue ./ homogeneousProjectedData(end,:)));
            center = mean(rescaledData, 2);
            
            centeredData = rescaledData - repmat(center, 1, groupSampleCount);
            [u s v] = svd(centeredData(1:3,:)',0);
            basis = v(:,1:2);
            projectedData = basis' * centeredData(1:3,:);
            
            max_pd = max(projectedData, [], 2);
            min_pd = min(projectedData, [], 2);

        ub1 = max_pd(1);
        ub2 = max_pd(2);
        lb1 = min_pd(1);
        lb2 = min_pd(2);
        tempArray = zeros(GRID_LINE_COUNT, GRID_LINE_COUNT, 3);  % Pages of this array are X, Y, and Z values.
        AIndex = 1;
        for A = linspace(lb1, ub1, GRID_LINE_COUNT),
            BIndex = 1;
            for B = linspace(lb2, ub2, GRID_LINE_COUNT),
                tempArray(AIndex, BIndex, :) = A * basis(:,1) + B * basis(:,2) + center(1:3);
                BIndex = BIndex + 1;
            end
            AIndex = AIndex + 1;
        end
        
        % Define vectors for the various group colors.
        % colors = 'rbkmgcy'
        groupColors = [1 0 0; 0 1 0; 0 0 0; 1 0 1; 0 1 0; 0 1 1; 1 1 0];
        %CData = groupColors(colorIndex,:);
        
%        h = surf(tempArray(:,:,1), tempArray(:,:,2), tempArray(:,:,3));
        h = surf(tempArray(:,:,1), tempArray(:,:,2), tempArray(:,:,3));
        set(h,'FaceColor',colors(colorIndex),'EdgeColor','black')
        set(h, 'EdgeAlpha', .2);
        alpha .2;
        
%         for i = 0:0.1:1
%             %% Draw the horizontal lines
%             l1 = line([ub1(1)+i*ub2(1)+(1-i)*lb2(1) lb1(1)+i*ub2(1)+(1-i)*lb2(1)], ... 
%                 [ub1(2)+i*ub2(2)+(1-i)*lb2(2) lb1(2)+i*ub2(2)+(1-i)*lb2(2)], ...
%                 [ub1(3)+i*ub2(3)+(1-i)*lb2(3) lb1(3)+i*ub2(3)+(1-i)*lb2(3)]);
%             set(l1, 'Color', colors(colorIndex));
%             %set(l1, 'LineWidth', 2.5);    
%             %% Draw the vertical lines
%             l2 = line([i*ub1(1)+(1-i)*lb1(1)+ub2(1) i*ub1(1)+(1-i)*lb1(1)+lb2(1)], ...
%                 [i*ub1(2)+(1-i)*lb1(2)+ub2(2) i*ub1(2)+(1-i)*lb1(2)+lb2(2)], ...
%                 [i*ub1(3)+(1-i)*lb1(3)+ub2(3) i*ub1(3)+(1-i)*lb1(3)+lb2(3)]);
%             set(l2, 'Color', colors(colorIndex));
%             %set(l2, 'LineWidth', 2.5);    
%         end
    end
    
    % Advance to the next color and marker.
    markerIndex = markerIndex + 1; if (markerIndex > markerCount) markerIndex = 1; end;
    colorIndex = colorIndex+1; if (colorIndex > colorCount) colorIndex = 1; end;
end;

xlabel('axis 1')
ylabel('axis 2')
zlabel('axis 3')


pl = [];