% Random Walk on Stanford Bunny Point Cloud
%
% Loads the Stanford Bunny mesh, extracts vertices as a point cloud,
% subsamples to a manageable size, builds a k-nearest-neighbor graph,
% and simulates a random walk on that graph.
%
% The bunny has non-uniform point density (ears are sparser than body),
% so the stationary distribution is NOT uniform — denser regions get
% visited more often. This contrasts with the uniform sphere example.
%
% Animation shows:
%   1. 3D bunny with walker position and trail
%   2. Visit frequency heatmap on the bunny
%   3. Histogram of visit counts
%
% Requires: bunny.ply in the same directory (Stanford Bunny PLY file)

function rw_bunny()
    %% =====================================================================
    %  1. Load Stanford Bunny point cloud from PLY file
    %  =====================================================================
    plyFile = fullfile(fileparts(mfilename('fullpath')), 'bunny.ply');
    if ~exist(plyFile, 'file')
        error('bunny.ply not found. Place the Stanford Bunny PLY file in: %s', ...
              fileparts(mfilename('fullpath')));
    end

    allPoints = readPLY(plyFile);
    fprintf('Loaded %d vertices from bunny.ply\n', size(allPoints, 1));

    % Subsample to ~1500 points for reasonable computation
    Ntarget = 1500;
    if size(allPoints, 1) > Ntarget
        idx = randperm(size(allPoints, 1), Ntarget);
        points = allPoints(idx, :);
    else
        points = allPoints;
    end
    N = size(points, 1);

    % Center and scale for nice display
    points = points - mean(points);
    scale = max(max(points) - min(points));
    points = points / scale * 2;    % fit roughly in [-1, 1]

    fprintf('Using %d points (subsampled)\n', N);

    %% =====================================================================
    %  2. Build k-nearest-neighbor graph
    %  =====================================================================
    k = 12;
    D = pdist2(points, points);         % N x N distance matrix
    [~, sortIdx] = sort(D, 2);
    neighbors = sortIdx(:, 2:k+1);      % N x k (exclude self)

    %% =====================================================================
    %  3. Set up figure
    %  =====================================================================
    fig = figure('Name','Random Walk on Stanford Bunny','Color','w',...
                 'Position',[80 60 1400 550]);

    % Compute axis limits from data
    pad = 0.15;
    xLim = [min(points(:,1))-pad, max(points(:,1))+pad];
    yLim = [min(points(:,2))-pad, max(points(:,2))+pad];
    zLim = [min(points(:,3))-pad, max(points(:,3))+pad];

    % Left panel (large): 3D walker view
    ax1 = subplot(1,2,1,'Parent',fig);
    hold(ax1,'on');
    scatter3(ax1, points(:,1), points(:,2), points(:,3), ...
             8, [0.7 0.7 0.7], 'filled', 'MarkerFaceAlpha', 0.25);
    hWalker = scatter3(ax1, NaN, NaN, NaN, 200, [1 0.84 0], ...
                       'filled', 'MarkerEdgeColor','k', 'LineWidth', 1.5);
    hTrail = plot3(ax1, NaN, NaN, NaN, '-', 'Color', [1 0.5 0 0.6], ...
                   'LineWidth', 2);
    axis(ax1, 'equal');
    set(ax1, 'XLim', xLim, 'YLim', yLim, 'ZLim', zLim);
    xlabel(ax1,'x'); ylabel(ax1,'y'); zlabel(ax1,'z');
    title(ax1, 'Random Walk on Stanford Bunny', 'FontSize', 14);
    view(ax1, 30, 20);
    set(ax1, 'XGrid','on','YGrid','on','ZGrid','on');

    % Top-right: visit frequency heatmap
    ax2 = subplot(2,2,2,'Parent',fig);
    hold(ax2,'on');
    hHeat = scatter3(ax2, points(:,1), points(:,2), points(:,3), ...
                     15, zeros(N,1), 'filled');
    colormap(ax2, hot);
    cb = colorbar(ax2);
    ylabel(cb, 'Visit count');
    axis(ax2, 'equal');
    set(ax2, 'XLim', xLim, 'YLim', yLim, 'ZLim', zLim);
    title(ax2, 'Visit Frequency', 'FontSize', 12);
    view(ax2, 30, 20);

    % Bottom-right: histogram of visit counts
    ax3 = subplot(2,2,4,'Parent',fig);
    hold(ax3,'on');
    title(ax3, 'Visit Count Distribution', 'FontSize', 12);
    xlabel(ax3, 'Visit count'); ylabel(ax3, 'Number of points');
    set(ax3, 'XGrid','on','YGrid','on');

    % Status text
    iterText = annotation(fig, 'textbox', [0.25 0.01 0.5 0.04], ...
        'String', 'Step: 0', 'FontSize', 12, 'EdgeColor', 'none', ...
        'HorizontalAlignment', 'center', 'FontWeight', 'bold');

    %% =====================================================================
    %  4. Simulate random walk
    %  =====================================================================
    nSteps = 8000;
    updateEvery = 25;
    trailLen = 40;

    % Initialize
    state = randi(N);
    visitCount = zeros(N, 1);
    trailBuf = nan(trailLen, 3);    % circular buffer for trail

    for t = 1:nSteps
        % Record visit
        visitCount(state) = visitCount(state) + 1;

        % Update trail buffer
        bufIdx = mod(t-1, trailLen) + 1;
        trailBuf(bufIdx, :) = points(state, :);

        % Move to random neighbor
        nbrs = neighbors(state, :);
        state = nbrs(randi(k));

        % Animate
        if mod(t, updateEvery) == 0 || t == 1
            % Walker position
            set(hWalker, 'XData', points(state,1), ...
                         'YData', points(state,2), ...
                         'ZData', points(state,3));

            % Trail: reorder circular buffer for plotting
            trailOrder = [bufIdx+1:trailLen, 1:bufIdx];
            trailPts = trailBuf(trailOrder, :);
            validRows = ~isnan(trailPts(:,1));
            trailPts = trailPts(validRows, :);
            set(hTrail, 'XData', trailPts(:,1), ...
                        'YData', trailPts(:,2), ...
                        'ZData', trailPts(:,3));

            % Visit frequency heatmap
            set(hHeat, 'CData', visitCount);
            if max(visitCount) > 0
                caxis(ax2, [0 max(visitCount)]);
            end

            % Histogram
            cla(ax3);
            hold(ax3, 'on');
            if max(visitCount) > 0
                edges = 0:max(1, ceil(max(visitCount)/25)):max(visitCount)+1;
                histogram(ax3, visitCount, edges, ...
                          'FaceColor', [0.3 0.6 0.9], 'EdgeColor', 'w');
                expectedMean = t / N;
                yl = ylim(ax3);
                plot(ax3, [expectedMean expectedMean], yl, 'r--', ...
                     'LineWidth', 2);
                legend(ax3, {'Visits', sprintf('Expected mean = %.1f', expectedMean)}, ...
                       'Location', 'northeast');
            end
            xlabel(ax3, 'Visit count'); ylabel(ax3, 'Number of points');
            title(ax3, 'Visit Count Distribution', 'FontSize', 12);
            set(ax3, 'XGrid','on','YGrid','on');

            % Rotate 3D views slowly
            az = 30 + t * 0.02;
            view(ax1, az, 20);
            view(ax2, az, 20);

            % Status
            cv = 0;
            if mean(visitCount) > 0
                cv = std(visitCount) / mean(visitCount);
            end
            set(iterText, 'String', sprintf(...
                'Step: %d / %d  |  Mean visits: %.1f  |  CV: %.3f', ...
                t, nSteps, mean(visitCount), cv));

            drawnow;
            pause(0.02);
        end
    end

    %% =====================================================================
    %  5. Save figure
    %  =====================================================================
    saveas(fig, fullfile(fileparts(mfilename('fullpath')), 'rw_bunny.png'));

    %% =====================================================================
    %  6. Console summary
    %  =====================================================================
    fprintf('\n=============================================\n');
    fprintf('  Random Walk on Stanford Bunny — Summary\n');
    fprintf('=============================================\n');
    fprintf('Points: %d (subsampled from %d),  k-NN: k=%d\n', ...
            N, size(allPoints,1), k);
    fprintf('Total steps: %d\n', nSteps);
    fprintf('---------------------------------------------\n');
    fprintf('Mean visits per point:  %.2f  (expected: %.2f)\n', ...
            mean(visitCount), nSteps/N);
    fprintf('Std of visit counts:    %.2f\n', std(visitCount));
    fprintf('Min / Max visits:       %d / %d\n', ...
            min(visitCount), max(visitCount));
    fprintf('Coeff of variation:     %.4f\n', std(visitCount)/mean(visitCount));
    fprintf('Points never visited:   %d / %d\n', sum(visitCount==0), N);
    fprintf('=============================================\n\n');
end

%% =========================================================================
%  Read vertices from an ASCII PLY file
%  =========================================================================
function pts = readPLY(filename)
    fid = fopen(filename, 'r');
    if fid == -1
        error('Cannot open file: %s', filename);
    end

    nVerts = 0;
    nProps = 0;

    % Parse header
    while true
        line = fgetl(fid);
        if startsWith(line, 'element vertex')
            nVerts = sscanf(line, 'element vertex %d');
        elseif startsWith(line, 'property') && nVerts > 0 && ...
               ~startsWith(line, 'property list')
            nProps = nProps + 1;
        elseif startsWith(line, 'element face')
            % Stop counting vertex properties
            break;
        elseif strcmp(line, 'end_header')
            break;
        end
    end

    % If we broke out at 'element face', keep reading to end_header
    if ~strcmp(line, 'end_header')
        while true
            line = fgetl(fid);
            if strcmp(line, 'end_header')
                break;
            end
        end
    end

    % Read vertex data (first 3 columns are x, y, z)
    data = zeros(nVerts, nProps);
    for i = 1:nVerts
        line = fgetl(fid);
        vals = sscanf(line, '%f');
        data(i, 1:length(vals)) = vals';
    end
    fclose(fid);

    pts = data(:, 1:3);    % extract x, y, z
end
