plot_subspaces3-1.m

一个用EM算法的源程序

% Code for rendering three dimensional results generated by the gpca
% algorithm.  Users should call plot_subspaces
function pl = plot_subspaces3(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));

    if (size(basis,2) == 1)
        % Plot the one dimensional subspace (line)
        projected_data = groupData' * basis;
        upper_bound = max(projected_data) * basis;
        lower_bound = min(projected_data) * basis;

        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) == 2)
        % 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,:)')';
        basis = normalize(basis);

        % Compute Bounds for the data.
        projected_data = (X') * basis;

        % The data bounds in the new basis.
        max_pd = max(projected_data, [], 1);
        min_pd = min(projected_data, [], 1);

        % Compute an array of values in three space in an order that can be
        % used for generating a surface plot.
        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);
                BIndex = BIndex + 1;
            end
            AIndex = AIndex + 1;
        end

        h = surf(tempArray(:,:,1), tempArray(:,:,2), tempArray(:,:,3));
        set(h,'FaceColor',colors(colorIndex),'EdgeColor','black')
        set(h, 'EdgeAlpha', .2);
        alpha .2;
        
    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 = [];