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;

Matlab: indexing by two bounds

OP.


Q: Suppose I have an array
v = 1:20
and two index bound vectors
i1 = [3 6 11 15];
i2 = [5 9 12 19];


I want to index the elements between i1(k) and i2(k)
so v(index) is [4 7 8 16 17 18]


A1:

index1 = cell2mat(arrayfun(@(x,y)(x+1:y-1),i1,i2,'uni',false));

A2:

s = (i2-i1) > 1;
x = zeros(size(v));
x(i1(s)+1) = 1;
x(i2(s)) = -1;

index2 = logical(cumsum(x));

Matlab: ismember in every row

Q:
Suppose I have

x = [1 2 3 4 5 6 7 8 9; 11 12 13 14 15 16 17 18 19]

I want to find out which row contains [3 7 9]

A1:

% Generate some data;

y = randi(50,1000,9);
t = [3 7 9];
% ensure that at least the first row has our target;
y(1:2,:) = x;

% Engine;

% This method assumes that the target is ascendingly sorted.
t1 = [3 3.1 7 7.1 9]; % 3.1 and 7.1 is arbitrary;
b = histc(y,t1,2);
idx1 = find(all(b(:,[1 3 5]),2));

A2:

% Engine;
idx2 = find(all(any(bsxfun(@eq,y,shiftdim(t,-1)),2),3));

Windows: Atbroker.exe failure workaround.

OP:
Atbroker.exe - Application Error in Vista


Quote:

"a workaround for now is to do CTRL-ALT-END and choose Logoff."

Matlab: grouping values according to ID

See also this post.

Q:

If we have n x 2 matrix x:

x = [11 19; 11 29; 10 34; 7 189; 6 123; 11 18; 21 182; 3 171; 6 145; 7 111]

Where the first column specifies the group ID and we want to break up the 2nd column into a cell array according to the grouping ID.

The desired output is:

c = {[171]; [145; 123]; [189; 111]; [34]; [18; 29; 19]; [182]}

A:

[u,dump,g] = unique(x(:,1))
c = accumarray(g,x(:,2),[length(u),1],@(x){x})

Matlab: handling nan for findstr

In many situations handling NaNs can be tricky if we want to write vectorized code.

For example, if I have vector

x = [1 2 3 nan 5 6 nan nan 9 10 11 12 nan 14]

If I want to locate 2 consequtive NaNs, I may use a combination of find, isnan, and diff.

However, what if I want to locate a subvector that also contains NaN? For instance, how to find

y = [3 nan 5]

in x?

If y does not contain NaNs, we know we can use findstr.

The straitforward solution is to convert NaNs to some value that doesn't appear in y.

For example,

sudoNan = max(x)+1;
x(isnan(x)) = sudoNan;
y(isnan(y)) = sudoNan;
findstr(x,y)

An even easier method is:

xc = typecast(x,'uint64');
yc = typecast(y,'uint64');
findstr(xc,yc)

Matlab: grouping and add up values according to ID

OP here.

Q:
Suppose I have n x 3 matrix x, I want a piece of code that checks the first and second columns and if they are equal then it adds up the correponding values in the third column.

x = [
1 2 3; 1 2 5; 1 2 3;...
1 3 1; 1 3 1; 2 1 1;...
2 1 1; 2 1 2; 2 1 2]

The desired answer should be

y = [
1 2 11; 1 3 2; 2 1 6]

A:

[val, dump, idx] = unique(x(:,1:2),'rows')
y = [val, accumarray(idx,x(:,3))]

Matlab: fast way to insert zeros into a vector

OP here.

Q:

Suppose I have vector

x = [1 2 3 4 5 6 7]

How do I insert zeros at every 3rd position to obtain

y = [1 2 0 3 4 0 5 6 0 7]

A1:

y = zeros(1,floor(length(x)*1.5));
z = 1:length(x);
y(z+ceil(z/2)-1) = x;

Matlab: indexing matrix of unknown number of dimensions

OP here.

Q:

If I know a multidimensional matrix x has 4 dims, I can index the first "row" of x using

x(1,:,:,:)

but how do I do this if I don't know ndims beforehand?

A:

c(1:ndims(x)-1)={':'};
x(1,c{:})

Matlab: counting zeros neighbouring each element

Q:

I have matrix x, I want to count for each element the number of zeros in the neighbouring 8 element.

For example, if

x = [
1 0 0
0 1 1
0 1 0];

I want to obtain count

c = [
2 2 1
2 5 3
1 3 0];

A:

c = conv2([1 1 1; 1 0 1; 1 1 1],double(~x),'same');

OR

c = conv2([1 1 1],[1 1 1],double(~x),'same')-~x;

Matlab: optimization

OP here.

Q:
Is there a way to optimize this for loop?

f = 0;
for i = 1:n
for j = i+1:n
z = y(i) - y(j) ;
f = f + z*z ;
end
end
f = 2*f;

A1:

yy = meshgrid(y);
f = sum(sum((yy - yy.').^2));

A2:
Assuming y is column vector:

f = 2*((n+1)*sum(y.^2)-2*cumsum(y')*y);

A little bit explanation:

Suppose

y = [y1 y2 y3 y4]

f =
(y1-y2).^2 + (y1-y3).^2 + (y1-y4).^2 +
(y2-y3).^2 + (y2-y4).^2 +
(y3-y4).^2
=
y1y1+y2y2-2*y1y2 + y1y1+y3y3-2*y1y3 + y1y1+y4y4-2*y1y4 +
y2y2+y3y3-2*y2y3 + y2y2+y4y4-2*y2y4 +
y3y3+y4y4-2*y3y4
=
y1y1 + y1y1 + y1y1 +
y2y2 + y2y2 + y2y2 +
y3y3 + y3y3 + y3y3 +
y4y4 + y4y4 + y4y4 - 2*(
y1y2 + y1y3 + y1y4 + y2y3 + y2y4 + y3y4
)
=
3*sum(y.^2) - 2*(
(y1).*y2
(y1+y2).*y3
(y1+y2+y3).*y4
)