1function [varargout] = rotate3D(rotation, varargin)
2%rotate3D rotate all objects (gradients) in the block by a rotation matrix
3%
4% [...] = rotate3D(rotation, obj <, obj> ...);
5%
6% Rotates the corresponding gradinet object(s) by the provided rotation
7% represented either as a 3x3 rotation matrix or a unit quaternion with
8% the scalar component at the first position, or a polar rotatiion
9% packaged into a 2D vector [phi, theta] expressed in radians or a single
10% scalar corresponding to a rotation about Z (e.g. phi) in radians; Keep
11% in mind that for a single Z axis rotation mr.rotate() can be faster.
12% Non-gradient objects are not affected.
13%
14% The objects that contain the field 'id' but are not modified by this
15% function are passed through without modification. For the gradients
16% that are modified, the field 'id' is removed automatically.
17%
18% Optional parameter list may include the keyword 'system' followed by a
19% system limits struct. The system can only be provided in the beginning
20% or at the ent of the list of optional parameters.
21%
22% Returns either a cell-array of objects if one return parameter is
23% provided or an explicit list of objects if multiple parameters are
24% given. Can be used directly as a parameter of seq.addBlock().
25%
26% See also mr.rotate, Sequence.addBlock
27%
28% Maxim Zaitsev <maxim.zaitsev@uniklinik-freiburg.de>
30% parse out the optional parameter 'system', which can only be at the beginning
31% or in the end of the optional parameters
32system=[];
33if ischar(varargin{1}) && strcmp(lower(varargin{1}),'system')
34 if ~isstruct(varargin{2}) || ~isfield(varargin{2},'gradRasterTime')
35 error('Error parsing input parameters, keyword ''system'' is not followed by a valid system struct');
36 end
37 system=varargin{2};
38 varargin=varargin(3:end);
39elseif length(varargin)>1 && ischar(varargin{end-1}) && strcmp(lower(varargin{end-1}),'system')
40 if ~isstruct(varargin{end}) || ~isfield(varargin{end},'gradRasterTime')
41 error('Error parsing input parameters, keyword ''system'' is not followed by a valid system struct');
42 end
43 system=varargin{end};
44 varargin=varargin(1:end-2);
45end
47% detect rotation matrix or quaternion formats
48if size(rotation)==[3 3]
49 rotMat=rotation;
50elseif length(rotation)==4
51 rotMat=mr.aux.quat.toRotMat(rotation);
52elseif length(rotation)==2
53 % phi & theta
54 phi = rotation(1);
55 theta = rotation(2);
56 assert( ( phi >= -pi ) && ( phi < 2*pi) , 'makeRotation:invalidTheta',...
57 'rotation angle phi (%.2f) is invalid. should be within [-pi,2*pi] radians',phi);
58 assert( ( theta >= -pi ) && ( theta <= pi) , 'makeRotation:invalidPhi',...
59 'rotation angle theta (%.2f) is invalid. should be within [-pi,pi] radians',theta);
60 q1=[cos(theta/2) 0 sin(theta/2) 0]; % y axis
61 q2=[cos(phi/2) 0 0 sin(phi/2)]; % z axis
62 rotQuaternion=mr.aux.quat.multiply(q2,q1); % ok, looks like the order is right
63 rotMat=mr.aux.quat.toRotMat(rotQuaternion);
64elseif numel(rotation)==1
65 % phi
66 phi = rotation;
67 assert( ( phi >= -pi ) && ( phi < 2*pi) , 'makeRotation:invalidTheta',...
68 'rotation angle phi (%.2f) is invalid. should be within [-pi,2*pi] radians',phi);
69 q=[cos(phi/2) 0 0 sin(phi/2)]; % z axis
70 rotMat=mr.aux.quat.toRotMat(q);
71else
72 error('The parameter ''rotation'' must either bi a 3x3 matrix or a quaternion');
73end
75% make this function accept ready-made blocks
76if isstruct(varargin{1}) && isfield(varargin{1}, 'rf')
77 varargin=mr.block2events(varargin);
78end
79% we need this to allow for nested mr.rotate() calls
80if 1==length(varargin) && iscell(varargin{1})
81 va=varargin{1};
82else
83 va=varargin;
84end
86% first create indexes of the objects to be bypassed or rotated
87ibypass=[];
88grads3_in=cell(1,3);
89axes={'x', 'y', 'z'};
91for i=1:length(va)
92 event = va{i};
93 if isempty(event)
94 continue;
95 end
96 if isnumeric(event) || ...
97 ((~strcmp(event.type,'grad') && ...
98 ~strcmp(event.type,'trap')))
99 ibypass=[ibypass i];
100 else
101 iAxis=find(strcmp(event.channel,axes));
102 if ~isempty(grads3_in{iAxis})
103 error('More than one gradient on the same axis %s provided', event.channel);
104 end
105 if isfield(event,'id')
106 grads3_in{iAxis}=rmfield(event,'id');
107 else
108 grads3_in{iAxis}=event;
109 end
110 end
111end
113max_mag=0; % measure of the relevant amplitude
114for i=1:3
115 if ~isempty(grads3_in{i})
116 max_mag=max(max_mag, getGradAbsMag(grads3_in{i}));
117 end
118end
119fthresh=1e-6;
120thresh=fthresh*max_mag;
122grads_out={};
123for j=1:3
124 grad_out_curr=[];
125 for i=1:3
126 if isempty(grads3_in{i}) || ...
127 abs(rotMat(j,i))<fthresh
128 continue;
129 end
130 g=mr.scaleGrad(grads3_in{i},rotMat(j,i));
131 g.channel=axes{j};
132 if isempty(grad_out_curr)
133 grad_out_curr=g;
134 else
135 if isempty(system)
136 grad_out_curr=mr.addGradients({grad_out_curr,g});
137 else
138 grad_out_curr=mr.addGradients({grad_out_curr,g},system);
139 end
140 end
141 end
142 % only output non-zero-amplitude gradients
143 if ~isempty(grad_out_curr) && getGradAbsMag(grad_out_curr) >= thresh
144 grads_out{end+1}=grad_out_curr;
145 end
146end
148% export
149bypass=va(ibypass);
150out={bypass{:},grads_out{:}};
152nout = nargout;
153varargout = cell(1,nout);
154if nout==1
155 varargout{1}=out;
156else
157 nr=min(nout,length(out));
158 if nout<length(out)
159 warning('insufficient number of return parameters, some rotated gradient components might go lost');
160 end
161 for k=1:nr
162 varargout{k}=out{k};
163 end
164end
166end
169function [out] = getGradAbsMag(grad)
170 if strcmp(grad.type,'trap')
171 out=abs(grad.amplitude);
172 else
173 out=max(abs(grad.waveform));
174 end
175end