FlyOrient

FLYORIENT
   Usage:
       orientation = FlyOrient(subset_frame, threshold)

 This function takes in a subset frame around the fly (calculated by
 FindFly) and discards the 3D data by placing points where the pixel intensity
 is larger than the user chosen threshold.  Then, Principal Components
 Analysis (by way fot the pca1 function) is performed on the resulting 
 scatter plot to find the direction of maximum variance --- this direction
 is taken to be the fly's (ambiguous) orientation.  
 orientation is a vector consisting of two angles (complements) that comprise 
 the body axis. The first element is an angle in the upper half plane; the
 second element is an angle in the lower half plane.

0001 function orientation  = FlyOrient(subset_frame, threshold)
0002 
0003 %FLYORIENT
0004 %   Usage:
0005 %       orientation = FlyOrient(subset_frame, threshold)
0006 %
0007 % This function takes in a subset frame around the fly (calculated by
0008 % FindFly) and discards the 3D data by placing points where the pixel intensity
0009 % is larger than the user chosen threshold.  Then, Principal Components
0010 % Analysis (by way fot the pca1 function) is performed on the resulting
0011 % scatter plot to find the direction of maximum variance --- this direction
0012 % is taken to be the fly's (ambiguous) orientation.
0013 % orientation is a vector consisting of two angles (complements) that comprise
0014 % the body axis. The first element is an angle in the upper half plane; the
0015 % second element is an angle in the lower half plane.
0016 
0017 % Written by Dan Valente
0018 % 11 October 2006
0019 
0020 %Normalize frame data by pixel of maximum intensity
0021 subset_frame = subset_frame/max(max(subset_frame));
0022 
0023 % Put dots where fly is and do PCA on reduced data set
0024 [rows, cols] = find(subset_frame >= threshold);
0025 rows = length(subset_frame(:,1))-rows+1;
0026 x = [cols';rows'];
0027 [xnew, PC, V, data] = pca1(x);
0028 
0029 % Find orientation vectors (two, mirrored across diagonal), and group into
0030 % upper half and lower half planes.
0031 a1 = PC(1,1);
0032 b1 = PC(2,1);
0033 a2 = -PC(1,1);
0034 b2 = -PC(2,1);
0035 if (b1 >= 0 );
0036     orientUHP = atan2(b1,a1);
0037     orientLHP = atan2(b2,a2);
0038 elseif (b2 >=0);
0039     orientUHP = atan2(b2,a2);
0040     orientLHP = atan2(b1,a1);
0041 else
0042 end
0043 
0044 % The vector we will return
0045 orientation = [orientUHP orientLHP];
0046 
0047 return;