% Code for rendering three dimensional results generated by the gpca% algorithm.  Users should call plot_subspacesfunction 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 = [];