function a = crossdiv(varargin)
%CROSSDIV Vector cross division.
% This function exploits the BAXFUN engine (MATLAB Central, file #23084),
% which enables multiple cross divisions and array expansion (AX).
%
% Orthogonal cross division (OCD):
% Aorth = CROSSDIV(C, B) is equivalent to Aorth = OCD(C, B)
% Aorth = CROSSDIV(C, B, ...) is equivalent to Aorth = OCD(C, B, ...)
%
% When B and C are 3-element vectors (e.g. 31, 13, or 113 arrays):
%
% A = CROSSDIV(TYPE, C, B, EXTRA), where C = A x B, returns a
% definite cross division C / B of type TYPE, with additional input
% argument EXTRA.
%
% More generally, when B and C are arrays of any size containing one or
% more 3-element vectors:
%
% A = CROSSDIV(TYPE, C, B, EXTRA) is equivalent to
% A = CROSSDIV(TYPE, C, B, EXTRA, IDC, IDB), where
% IDA and IDB are the first dimensions of A and B whose length is 3.
%
% A = CROSSDIV(TYPE, C, B, EXTRA, DIM) is equivalent to
% A = CROSSDIV(TYPE, C, B, EXTRA, IDC, IDB), where IDC = IDB = DIM.
%
% A = CROSSDIV(TYPE, C, B, EXTRA, IDC, IDB), where C = A x B, returns
% the definite cross divisions between the vectors contained in C
% along dimension IDC and those contained in B along dimension IDB.
% These vectors must have 3 elements. C and B are viewed as "block
% arrays". IDC and IDB are referred to as their "internal dimensions"
% (IDs). For instance, a 362 array may be viewed as an array
% containing twelve 3-el. blocks. In this case, its size is denoted
% by (3)62, and its ID is 1. Since AX is enabled, C and B may have
% different size, and IDC may not coincide with IDB (see MULTIPROD).
%
% TYPE A string containing one specifier for the type of
% additional input (EXTRA), and optionally one specifier
% for the type of missing input, for minimum-input CDs.
% White spaces are ignored. TYPE is case insensitive.
% Specifiers for the type of EXTRA:
% '+s' Cross division coefficient.
% '+theta' (or '+t') Angle between A and B (radians).
% '+Ax' (or '+x') Known component of vector A.
% '+Ay' (or '+y')
% '+Az' (or '+z')
% Specifiers for the type of missing input:
% '' No missing input.
% '-Bx' Unknown component of vector B.
% '-By'
% '-Bz'
% '-Cx' (or '-x') Unknown component of vector C.
% '-Cy' (or '-y')
% '-Cz' (or '-z')
% A, B, C Array Block size Internal dimension
% ---------------------------------------------------------
% C 3 (1-D) IDC
% B 3 (1-D) IDB
% A 3 (1-D) MAX(IDC, IDB)
% ---------------------------------------------------------
% For minimum-input CDs, an arbitrary value (e.g. NaN) must
% be assigned to the missing component of each vector
% contained in B or C. If SIZE(B) == SIZE(C) and
% IDB == IDC, then SIZE(A) = SIZE(B) = SIZE(C).
% EXTRA Additional input value(s) (see TYPE). EXTRA must be
% either a scalar or an array of the same length as A along
% all dimensions except IDA. SIZE(EXTRA, IDA) must be 1.
% IDC, IDB The "internal dimensions" along which the vectors are
% contained in A and B, respectively.
%
% Two main kinds of cross divisions (CDs) exist: indefinite (ICDs) and
% definite (DCDs). DCDs allow to write in a short vectorial form
% equations such as those used to compute the point of application of a
% force based on dynamometric data (see examples). The definition of the
% DCD is based on that of the ICD. The symbol of the ICD is //. As well
% as indefinite integrals yield infinite antiderivatives, ICDs yield
% infinite anticrossproducts (antiCPs). There are infinite vectors X
% meeting the condition X x B = C, all laying on a plane normal to C. If
% they include A, they are called the "antiCPs of A x B, where A x B =
% C. These antiCPs are determined by C // B, defined as follows:
% C // B = X = Aorth + Xpar (1)
% Aorth = B x C / (B * B) (Orthogonal term)
% Xpar = S * UNIT(B) (Parallel term)
% Each antiCP X shares with A its component Aorth = Xorth, orthogonal to
% both B and C, and differs from the others only by Xpar, its component
% parallel to B. Since Aorth is unique, it can be determined as shown
% above. Notice that Aorth is one of the antiCPs (Aorth x B = C).
% Although Xpar is nonunique, it can be partially determined. Being it
% parallel to B, it is equal to S * UNIT(B), where S is a scaling factor,
% called "cross division coefficient". Thus, the only fully undetermined
% value in equation 1 is S. Luckily, in some cases S is known or it can
% be determined otherwise, and this allows for the definition of four
% kinds of DCDs (see below). In the definition of an ICD, S has the same
% role as the arbitrary integration constant in an indefinite integral.
%
% Orthogonal cross division (OCD)
% The OCD is the simplest form of DCD. It is determined by setting
% S = 0, and is equal to the unique orthogonal term of an ICD:
% C /o B = Aorth = B x C / (B * B) (2)
% When the orthogonal antiCP Aorth happens to coincide with A,
% C /o B = A (2b)
%
% Basic DCD
% When the value Sa for S which uniquely identifies A in eq. 1 is
% known, A can be fully determined as follows:
% C /s B = Aorth + Apar = C /o B + Sa * UNIT(B) = A (3)
%
% Computed-coefficient DCDs
% In equation 3, Sa can be computed as a function of:
% - THETA (angle between A and B) or
% - a component of A along an axis non-orthogonal to B.
%
% Minimum-input DCDs
% In any of the above mentioned DCDs, one of the components of B or C
% can be omitted, provided it is relative to an axis non-orthogonal
% to C or B. Since B and C are by definition orthogonal, the missing
% component can be computed by solving the equation B * C = 0.
%
% Warning:
% By definition, in an ICD (C // B) the divisor (B) must be the 2nd
% operand of a cross product (A x B = C). When you use as divisor the
% 1st operand (A) of the cross product, rather than the 2nd, the ICD
% does not yield the expected set of antiCPs (X including B).
% Since A x B = -B x A, it yields the opposite set (X including B):
% C // A = - X = - Borth - S * UNIT(A) (4)
% Being derived by the ICD, all the DCDs share the same limitation.
% E.g.,
% C /o A = - Borth (5)
%
% Reference:
% de Leva, P (2008). Anticrossproducts and cross divisions.
% Journal of Biomechanics, 8, 1790-1800.
%
% Examples:
% If the moment of a force F is R x F = M:
% 1) R = CROSSDIV('+t', M, F, THETA) [Computed-coeff. DCD]
% 2) R = CROSSDIV('+z', M, F, Rz) [Computed-coeff. DCD]
% 3) R = CROSSDIV('+z-z', M, F, Rz) [Minimum-input DCD]
% 4) F = - CROSSDIV('+t', M, R, THETA) [Equation 4]
% See also the examples in the help text of function OCD.
% NOTE 1 - When Rz is used as input argument, axis z must be
% non-orthogonal to F (i.e. Fz must be non-null).
% NOTE 2 - The third instruction can be used to compute the point of
% application of a force based on the output of a force
% plate (typically, z is vertical, Rz known, Mz unknown).
%
% If C and B are both ................ (3)62 arrays of vectors,
% and EXTRA is either a scalar or .. a (1)62 array
% C = CROSSDIV(TYPE, C, B, EXTRA) is a (3)62 array of vectors.
%
% A single vector C is divided by thirty vectors contained in B:
% If C is .......................... a (3)1 vector,
% B is .......................... a 56(3) array of 30 vectors,
% and EXTRA is either a scalar or ... a 56(1) array
% C = CROSSDIV(TYPE,C,B,EXTRA,1,3) is a 56(3) array of 30 vectors.
%
% See also OCD, DOT2, CROSS2, OUTER, MAGN, UNIT, PROJECTION, REJECTION,
% TESTXDIV.
% $ Version: 3.0 $
% CODE by: Paolo de Leva (IUSM, Rome, IT) 2006 Jul 20
% optimized by: Code author 2009 Feb 12
% COMMENTS by: Code author 2009 Feb 26
% OUTPUT tested by: Code author 2009 Feb 26
% -------------------------------------------------------------------------
% SIMPLEST CASE: A simple OCD (no missing value)
if ~ischar(varargin{1})
a = ocd(varargin{:});
return
end
% Checking and adjusting input arguments
error( nargchk(4, 6, nargin) );
[type, c, b, extra] = varargin{1:4};
type = type( type~=' ' ); % Removing white spaces
type = lower(type); % Lower case
sizeC = size(c);
sizeB = size(b);
sizeE = size(extra);
% Dividing TYPE into two substrings
findplus = strfind(type, '+'); % Mandatory specifier
findminus = strfind(type, '-'); % Optional specifier
if length(findplus)~=1 || length(findminus)>1
error('CROSSDIV:InvalidType', ['Invalid TYPE: ', type]);
elseif isempty( findminus )
extratype = type;
missingtype = '';
elseif findplus > 1
missingtype = type( 1 : findplus-1 );
extratype = type( findplus : end );
else % findminus > 1
extratype = type( 1 : findminus-1 );
missingtype = type( findminus : end );
end
% Setting flags specifying type of EXTRA
switch extratype
case {'+ax', '+x'}, extraflag = 'x'; extraaxis = 1;
case {'+ay', '+y'}, extraflag = 'y'; extraaxis = 2;
case {'+az', '+z'}, extraflag = 'z'; extraaxis = 3;
case '+s' , extraflag = 's'; extraaxis = 0;
case {'+theta', '+t'}, extraflag = 't'; extraaxis = 0;
otherwise, error('CROSSDIV:InvalidType', ['Invalid TYPE: ', type]);
end
% Setting flag specifying type of missing value
switch missingtype
case '', missingaxis = 0;
case '-bx', missingaxis = 1;
case '-by', missingaxis = 2;
case '-bz', missingaxis = 3;
case {'-cx', '-x'}, missingaxis = 4;
case {'-cy', '-y'}, missingaxis = 5;
case {'-cz', '-z'}, missingaxis = 6;
otherwise, error('CROSSDIV:InvalidType', ['Invalid TYPE: ', type]);
end
% Setting IDC and/or IDB
switch nargin
case 4
idC = find(size(c)==3, 1, 'first'); % First dim. of length 3
idB = find(size(b)==3, 1, 'first');
if isempty(idC) || isempty(idB)
error('CROSSDIV:InvalidSize',...
'C and B must have at least one dimension of length 3.');
end
case 5
idC = varargin{5};
idB = idC;
case 6
[idC, idB] = varargin{[5 6]};
end
% Checking block size
if (sizeC(idC)~=3) || (sizeB(idB)~=3),
error('CROSSDIV:InvalidBlockSize',...
['C and B must be of length 3 in the dimensions\n'...
'in which the cross division is taken.'])
end
% Dimension shift (first step of array expansion)
% NOTE: The BAXFUN engine is implemented here without calling
% BAXFUN, to avoid repeated shift of B and C
diff = idC - idB;
if diff == 0
id = idB; % ID = IDA = IDB = IDC
elseif diff > 0
id = idC;
sizeB = [ones(1, diff) sizeB];
b = reshape(b, sizeB);
% if sizeB(end)==1, sizeB = sizeB(end-1); end % when B is 31
else % diff < 0
id = idB;
sizeC = [ones(1,-diff) sizeC];
c = reshape(c, sizeC);
% if sizeC(end)==1, sizeC = sizeC(end-1); end % when C is 31
end
% Computing size of A
NdimsB = length(sizeB);
NdimsC = length(sizeC);
NdimsA = max(NdimsB, NdimsC);
NsinglB = NdimsA - NdimsB; % Number of added trailing singletons
NsinglC = NdimsA - NdimsC;
adjsizeB = [sizeB ones(1,NsinglB)];
adjsizeC = [sizeC ones(1,NsinglC)];
sizeA = max(adjsizeB, adjsizeC) .* (adjsizeB>0 & adjsizeC>0);
% Checking size of EXTRA
fullsizeE = sizeA;
fullsizeE(id) = 1;
if ~isequal(sizeE, [1 1])
NsinglE = NdimsA - length(sizeE);
adjsizeE = [sizeE ones(1,NsinglE)];
if ~isequal(adjsizeE, fullsizeE)
error('CROSSDIV:InvalidSizeExtra', 'Invalid size of EXTRA.');
end
end
% 1 - Minimum-input CDs
if missingaxis
notmissing = [1 2 3 1 2 3 1 2];
notmissing1 = notmissing(missingaxis + 1);
notmissing2 = notmissing(missingaxis + 2);
% Vectorized indices
% NOTE: If By or Cy is missing, it is important (in step 1a) that
% the two NOTMISSING axes are [3 1] rather than [1 3]
idxA = ivector(sizeA);
idxB = ivector(sizeB);
idxC = ivector(sizeC);
% 1a - Simplified minimum-input CD
if missingaxis > 3 && (missingaxis-3) == extraaxis
a = zeros(sizeA);
% Known component of A
idxA{id} = extraaxis;
a(idxA{:}) = extra + zeros(fullsizeE);
% A1 = EXTRA*B1 - C2
idxA{id} = notmissing1;
idxB{id} = notmissing1;
idxC{id} = notmissing2;
temp = bsxfun( @times, extra, b(idxB{:}) );
a(idxA{:}) = bsxfun( @minus, temp, c(idxC{:}) );
% A2 = EXTRA*B2 + C1)
idxA{id} = notmissing2;
idxB{id} = notmissing2;
idxC{id} = notmissing1;
temp = bsxfun( @times, extra, b(idxB{:}) );
a(idxA{:}) = bsxfun( @plus, temp, c(idxC{:}) );
% Unknown components (A1,A2) = (A1,A2) / B0
idxA{id} = [notmissing1 notmissing2];
idxB{id} = extraaxis;
a(idxA{:}) = bsxfun( @rdivide, a(idxA{:}), b(idxB{:}) );
return
% 1b - Computing the missing component of B or C by solving B*C = 0
else
idxB{id} = notmissing1;
idxC{id} = notmissing1;
temp = bsxfun( @times, b(idxB{:}), c(idxC{:}) ); % B1 * C1
idxB{id} = notmissing2;
idxC{id} = notmissing2;
temp = temp + ...
bsxfun( @times, b(idxB{:}), c(idxC{:}) ); % B2 * C2
% B0 = - (B1*C1 + B2*C2) / C0
if missingaxis < 4
if ~isequal(sizeA, sizeB)
b = bsxfun(@plus, b, zeros(sizeA)); % Expanding B
idxB = idxA;
end
idxB{id} = missingaxis;
idxC{id} = missingaxis;
b(idxB{:}) = - bsxfun( @rdivide, temp, c(idxC{:}) );
% C0 = - (B1*C1 + B2*C2) / B0
else % missingaxis > 3
if ~isequal(sizeA, sizeC)
c = bsxfun(@plus, c, zeros(sizeA)); % Expanding C
idxC = idxA;
end
idxB{id} = missingaxis - 3;
idxC{id} = missingaxis - 3;
c(idxC{:}) = - bsxfun( @rdivide, temp, b(idxB{:}) );
end
end
end
clear temp
% 2 - Orthogonal cross division (A = Aorth = C /o B)
a = ocd(c, b, id);
% 3 - Computed-coefficient DCD
switch extraflag
case 't' % Angle-input CD (Sa = MAGN(Aorth) / TAN(EXTRA))
SimplifiedBasic = false;
if isequal(sizeE, [1 1]) && mod(extra, pi) == 0.5*pi
return; % OCD (A = Aorth)
end
sa = bsxfun( @rdivide, magn(a, id), tan(extra) );
case {'x','y','z'} % Component-input CD (SA0 = (EXTRA-Aorth0) / B0)
SimplifiedBasic = true;
if ~missingaxis
idxA = ivector(sizeA);
idxB = ivector(sizeB);
end
idxA{id} = extraaxis;
idxB{id} = extraaxis;
sa0 = extra - a(idxA{:});
sa0 = bsxfun( @rdivide, sa0, b(idxB{:}) );
case 's' % (Basic DCD)
SimplifiedBasic = false;
sa = extra;
end
% 4 - Basic DCD (A = Aorth + Apar = Aorth + Sa * UNIT(B))
if SimplifiedBasic
Apar = bsxfun( @times, sa0, b );
a = a + Apar;
a(idxA{:}) = extra; % (removes rounding errors)
else
Apar = bsxfun( @times, sa, unit(b, id) );
a = a + Apar;
end
%--------------------------------------------------------------------------
function indices = ivector(sizeA)
%IVECTOR Vectorizing the indices of an array.
Ndims = length(sizeA);
indices = cell(1,Ndims); % preallocating
for d = 1 : Ndims
indices{d} = 1:sizeA(d);
end
