NIRS_ICA_V1.0/ica_algorithms/tentative/STICA_skew_demo/STICA_skew_demo/STICA_UTIL/jjrotate.m
function jjrotate_surf_dots(h,azel,alpha,origin,xd,yd,zd)
%JJROTATE Rotate objects about specified origin and direction.
%   JJROTATE(H,[THETA PHI],ALPHA) rotates the objects with handles H
%   through angle ALPHA about an axis described by the 2-element
%   direction vector [THETA PHI] (spherical coordinates).  
%   All the angles are in degrees.  The handles in H must be children
%   of the same axes.
%
%   THETA is the angle in the xy plane counterclockwise from the
%   positive x axis.  PHI is the elevation of the direction vector
%   from the xy plane (see also SPH2CART).  Positive ALPHA is defined
%   as the righthand-rule angle about the direction vector as it
%   extends from the origin.
%
%   ROTATE(H,[X Y Z],ALPHA) rotates the objects about the direction
%   vector [X Y Z] (cartesian coordinates). The direction vector
%   is the vector from the center of the plot box to (X,Y,Z).
%
%   ROTATE(...,ORIGIN) uses the the point ORIGIN = [x0,y0,y0] as
%   the center of rotation instead of the center of the plot box.
%
%   See also SPH2CART, CART2SPH.

% Determine the default origin (center of plot box).
if nargin < 4
  ax = get(h(1),'parent');
  if isempty(ax) | ax==0,
    error('H must contain axes children only.');
  end
  origin = sum([get(ax,'xlim')' get(ax,'ylim')' get(ax,'zlim')'])/2;
end

origin

% find unit vector for axis of rotation
if prod(size(azel)) == 2 % theta, phi
    phi = pi*azel(1)/180;
    theta = pi*azel(2)/180;
    u = [cos(theta)*cos(phi); cos(theta)*sin(phi); sin(theta)];
elseif prod(size(azel)) == 3 % direction vector
    u = azel(:)/norm(azel);
    theta = asin(u(3));
    phi = atan2(u(2),u(1))-pi/2;
end

u'

alph = alpha*pi/180;
cosa = cos(alph);
sina = sin(alph);
vera = 1 - cosa;
x = u(1);
y = u(2);
z = u(3);
rot = [cosa+x^2*vera x*y*vera-z*sina x*z*vera+y*sina; ...
       x*y*vera+z*sina cosa+y^2*vera y*z*vera-x*sina; ...
       x*z*vera-y*sina y*z*vera+x*sina cosa+z^2*vera]';

for i=1:prod(size(h)),
  t = get(h(i),'type');
  skip = 0;
 % strcmp(t,'surface')==TRUE
  if strcmp(t,'surface') | strcmp(t,'patch') | strcmp(t,'line')
    x = get(h(i),'xdata')-origin(1);
    y = get(h(i),'ydata')-origin(2);
    z = get(h(i),'zdata')-origin(3);
    if isempty(z)
       z = -origin(3)*ones(size(y));
    end
    [m,n] = size(z);
    if prod(size(x)) < m*n % NOT TRUE
      [x,y] = meshgrid(x,y);
    end
  elseif strcmp(t,'text')
    p = get(h(i),'position');
    x = p(1); y = p(2); z = p(3);
  elseif strcmp(t,'image')
    x = get(h(i),'xdata');
    y = get(h(i),'zdata');
    z = zeros(size(x));
  else
    skip = 1;
  end

  if ~skip,
    [m,n] = size(x);
    newxyz = [x(:), y(:), z(:)];
    % size(newxyz) % == 2601x3
    newxyz = newxyz*rot;
    newx = origin(1) + reshape(newxyz(:,1),m,n);
    newy = origin(2) + reshape(newxyz(:,2),m,n);
    newz = origin(3) + reshape(newxyz(:,3),m,n);
  end

  if strcmp(t,'surface') | strcmp(t,'patch') | strcmp(t,'line')
    set(h(i),'xdata',newx,'ydata',newy,'zdata',newz);
  elseif strcmp(t,'text')
    set(h(i),'position',[newx newy newz])
  elseif strcmp(t,'image')
    set(h(i),'xdata',newx,'ydata',newy)
  end
end
