% mainhyper.m
m = 100;
rand ('seed', 10);
randn ('seed', 0);
P = rand (2, m) - 0.5;
P (3, :) = 0.5 * P (1, :) + 0.5 * P (2, :) + ...
           0.25 * randn (1, m);

[V, p] = hyper (P');

proj = p * ones (1, m) + ...
       V (:, 1:2) * V (:, 1:2)' * (P - p * ones (1, m));

corners = [1 -1 -1  1
           1  1 -1 -1];
corners = p * ones (1, 4) + V (:, 1:2) * corners;

straight = p * ones (1, 2) + V (:, 1) * [-1 1];

nframes = 36;  phi = 2 * pi / nframes;
rot = [cos(phi) -sin(phi) 0
       sin(phi)  cos(phi) 0
       0         0        1];
clf;  axis ('off');  colormap (gray);
frames = moviein (nframes);

for k = 1:nframes,
   clf;
   axis ([-1.25 1.25 -1.25 1.25 -1.25 1.25]);  axis ('off');
   view (viewmtx (0, 0, 0));  caxis ([-1.3 1.3]);

   % plot points behind the plane
   for j = 1:m,
      if P (2, j) > proj (2, j),
         line ([P(1,j) proj(1,j)], ...
               [P(2,j) proj(2,j)], ...
               [P(3,j) proj(3,j)], ...
               'Color', 'yellow', 'EraseMode', 'none');
      end
   end
   drawnow;

   % plot plane
   patch (corners (1, :), corners (2, :), corners (3, :), ...
          -corners (2, :), ...
          'EdgeColor', 'none', 'EraseMode', 'none');
   drawnow;

   % plot least squares line
   line (straight (1, :), straight (2, :), straight (3, :), ...
         'Color', 'red', 'EraseMode', 'none');
   drawnow;

   % plot points before the plane
   for j = 1:m,
      if P (2, j) <= proj (2, j),
         line ([P(1,j) proj(1,j)], ...
               [P(2,j) proj(2,j)], ...
               [P(3,j) proj(3,j)], ...
               'Color', 'yellow', 'EraseMode', 'none');
      end
   end
   drawnow;

   frames (:, k) = getframe;

   % rotate points
   P = rot * P;  proj = rot * proj;
   corners = rot * corners;  straight = rot * straight;
end;
clf;  axis ('off');  movie (frames, inf);

