Matlab: operate on matched dimension with different shaped matrices

OP.

Q: Suppose I have

a = [1 2 3; 7 8 9; 6 5 4];

b = [7 5; 1 3; 9 4];

I want calculate the product between every column of a and b. Is there a vectorized way of this loop:

p = zeros(3,6);

c = 0;

for k1 = 1:3

for k2 = 1:2

c = c+1;

p(:,c) = a(:,k1).*b(:,k2);

end

end

A:

[m,n] = size(a);

p = bsxfun(@times,reshape(a,[m 1 n]),b);

p = reshape(p,n,[]);

Note:

permute(a,[1 3 2]) has the same effect as reshape(a,[m 1 n]) , but reshape is much faster.

Geometry in 3D: Orthonormal Cartesian basis transformation

Transforming the expression of a vector from one orthonormal Cartesian basis frame cf1 to another orthonormal Cartesian basis cf2, given cf2In1, the expression of basis cf2 in cf1;

Orthonormal transformation is no other than a rotation around the origin (plus possibly some flip of axes);

Suppose

cf2In1 == [c1;c2;c3] == [c1x,c1y,c1z;c2x,c2y,c2z;c3x,c3y,c3z]; % c1,c2,c3 are the basis vectors;

orthonormality requires that

cf2In1 * cf2In1 .' == eye(3) % Identity; because dot(c1,c1) == 1; dot(c1,c2) == 0; etc;

it can be shown, via some straightforward linear algebra, that coordinates expressed in cf1 transformed into cf2 is just

vIn2 == vIn1*tr1To2 == vIn1*(cf2In1.'); % if orthonormality holds, by doing dot product with each axis of cf2In1 we get the transformation here; notice the transposition that makes the dot product appropriate (assume vIn1 being 1 x 3 vector);

On the other hand, if we have coordinates valued in cf2 and want its expression in cf1:

% figure out the inverse transformation;

vIn1 == vIn2*tr2To1 == vIn2*inv(cf2In1.');

however since cf2In1*cf2In1.' == eye(3) we have inv(cf2In1.') == cf2In1;

thus the transformation is simply cf2In1:

vIn1 == vIn2*cf2In1; % tr2To1 == cf2In1;

Some numeric examples:

un = @(x)(x./norm(x)); % make norm 1;

tocol = @(x)(x(:));

% suppose n1 = [1 0 0] and cf1 = eye(3), the default Cartesian frame;

n2 = un(randn(1,3)); % x-direction of cf2In1;

% get cf2In1, using qr;

[cf2In1,r] = qr(tocol(n2)); % r is 3 x 1 with (eigen) values [1;0;0]; cf2In1 are orthonormal basis with first row equal to n2; each row forms one axis of the basis;

% now suppose we have some vector vIn2 that can be expressed as [2 3 4] in cf2; we want to know vIn1, its expression in cf1; intuitively, we can multiply the ratio to each axis in cf2, then find the total projection on coordinate frame 1; which can be expressed as:

vIn2 = [2 3 4]

sum([2*cf2In1(1,:);3*cf2In1(2,:);4*cf2In1(3,:)],1)

% or equivalently

vIn1 = vIn2*cf2In1

% check if this is indeed the answer; find the expression of vIn1 in cf2;

vIn2 = vIn1*(cf2In1.')

Suppose the scenario that v is defined in cf2; we need to express v in cf3In1;

suppose we have cf3In1

[cf3In1,r] = qr(randn(1,3)); % r(1) is also the L2 norm of the random vector;

first express v in cf1;

vIn1 == vIn2*cf2In1;

next express vIn1 in cf3;

vIn3 == vIn1*(cf3In1.');

thus

vIn3 == vIn2*cf2In1*(cf3In1.');

tr2To3 == cf2In1*cf3In1.';

tr2To3 is essentially a rotation from cf2In1 to cf3In1

we established a direct way to get the rotation from cf2In1 to cf3In1:

tr2To3 == cf2In1*cf3In1.';

or from a vector vA to another vector vB all in default coordinate frame 1:

[cfA,~] = qr(vA(:));

[cfB,~] = qr(vB(:));

trAToB = cfA*cfB.';

Geometry in 3D: projection

A vector a projected onto another vector b effectively gives a vector with same orientation as b and with magnitude norm(a).*cos(span(a,b));

define

unitize = @(a)(a./norm(a));

which gives the unit vector along a;

The projection can be equivalently written in the original coordinate frame as

projection(a,b) == dot(a,b)./norm(b).*unitize(b) == dot(a,b).*b./(norm(b)^2);

Recall that

norm(b)^2 == torow(b)*tocol(b);

dot(a,b) == torow(a)*tocol(b);

suppose a, b are 1 x 3 row vectors

We have

projection(a,b) == (a*b.').*b./(b*b.') == a*(b.'*b)./(b*b.') == a*pm;
pm == (b.'*b)./(b*b.');

pm is called the projection matrix of b;

Geometry in 3D: dot, cross, norm

First we define some utility functions;

torow = @(x)(x(:).');

tocol = @(x)(x(:));

norm (L2, Euclidean) of a vector is a scalar defined as

norm(a) == sqrt(sum(a.^2)) == sqrt(torow(a)*tocol(a)) == sqrt(dot(a,a));

dot product of two vectors is a scalar,

dot(a,b) == sum(a.*b) == torow(a)*tocol(b);

cross product of two vectors is a vector orthogonal to both, following the right hand rule:

cross(a,b) == ([a(2)*b(3)-a(3)*b(2),...

a(3)*b(1)-a(1)*b(3),...

a(1)*b(2)-a(2)*b(1)]); % each row corresponds to an axis;

Right hand rule:

(Image from Wikipedia)

Some very commonly used properties:

dot(a,b) == norm(a).*norm(b).*cos(span(a,b));
span(a,b) == acos(dot(a,b)./norm(a)./norm(b));

where span() is the angle spanned between two vectors;

norm(cross(a,b)) == norm(a).*norm(b).*sin(span(a,b));

span(a,b) == atan2(norm(cross(a,b)),dot(a,b));

This formula of span(a,b) is more robust than the acos version at near pi/2; the convention of such angle is measured by the shortest great circle path between two vectors;

Matlab: column circshift by each row

Q:
Suppose T is a triangular mesh, f is the face list, v is the vertex list;

For some vertex k, I have the info of its neighboring faces va, which is an index into f.

fva = f(va,:) is length(va) by 3 matrix about vertices involved in va, and k is present in every row of fva.

I want to re-arrange columns in fva, row by row, so that vertex k is always in the first column, while the rotation order of vertices (clockwise,cc) is retained for each row.

% Set up the data;
load seamount;
v = [x,y,z];
f = delaunay(x,y);
k = 10;
va = find(any(f == 10,2));
fva = f(va,:);

A1: straightforward for-loop solution;

for c = 1:length(va)
if fva(c,2) == k, fva(c,:) = fva(c,[2 3 1]); end
if fva(c,3) == k, fva(c,:) = fva(c,[3 1 2]); end
end

A2:

fvaT = fva.';

[ii,~] = find(fvaT == 10);

cOrder = [1 2 3; 2 3 1; 3 1 2];

fva = fvaT(bsxfun(@plus,cOrder(:,ii),(0:3:numel(ii)*3-3))).';

Generally A1 is faster than A2;