Повторите вектор в 4 измерениях, не используя repmat

Question

Я должен повторить вектор в 4 измерениях, и я использовал для этого функцию repmat, но это занимает много времени. Итак, как я могу реализовать это без использования repmat?

Вот пример кода, который я использую.

A = rand(1, 240); %The vector I want to replicate

B = repmat(A, 240, 1, 180, 3); %Here I get a matrix of size 240 240 180 3

Вот функция, в которой я вызываю repmat. Функция в основном реализует три вложенные суммы:

function gr = g_transformation3(f, c, delta)
%f --> image to be transformed
%c --> control points matrix
%delta --> spacing in pixels, between the control points

[X, Y, Z] = size(f);
%[n, m] = size(c);

x0 = 0:(X - 1);
y0 = 0:(Y - 1);
z0 = 0:(Z - 1);

x1 = delta:(X + delta - 1);
y1 = delta:(Y + delta - 1);
z1 = delta:(Z + delta - 1);

x2 = (2*delta):(X + 2*delta - 1);
y2 = (2*delta):(Y + 2*delta - 1);
z2 = (2*delta):(Z + 2*delta - 1);

x3 = (3*delta):(X + 3*delta - 1);
y3 = (3*delta):(Y + 3*delta - 1);
z3 = (3*delta):(Z + 3*delta - 1);

i0 = floor(x0./delta) + 1;
j0 = floor(y0./delta) + 1;
k0 = floor(z0./delta) + 1;

i1 = floor(x1./delta) + 1;
j1 = floor(y1./delta) + 1;
k1 = floor(z1./delta) + 1;

i2 = floor(x2./delta) + 1;
j2 = floor(y2./delta) + 1;
k2 = floor(z2./delta) + 1;

i3 = floor(x3./delta) + 1;
j3 = floor(y3./delta) + 1;
k3 = floor(z3./delta) + 1;

u = x0./delta - floor(x0./delta);
v = y0./delta - floor(y0./delta);
w = z0./delta - floor(z0./delta);

Bu = zeros(4, numel(u));
Bv = zeros(4, numel(v));
Bw = zeros(4, numel(w));

%Computing the vector of B-splines
Bu(1, :) = ((1 - u).^3)/6;
Bu(2, :) = (3*u.^3 - 6*u.^2 + 4)/6;
Bu(3, :) = (-3*u.^3 + 3*u.^2 + 3*u + 1)/6;
Bu(4, :) = (u.^3)/6;

Bv(1, :) = ((1 - v).^3)/6;
Bv(2, :) = (3*v.^3 - 6*v.^2 + 4)/6;
Bv(3, :) = (-3*v.^3 + 3*v.^2 + 3*v + 1)/6;
Bv(4, :) = (v.^3)/6;

Bw(1, :) = ((1 - w).^3)/6;
Bw(2, :) = (3*w.^3 - 6*w.^2 + 4)/6;
Bw(3, :) = (-3*w.^3 + 3*w.^2 + 3*w + 1)/6;
Bw(4, :) = (w.^3)/6;

T00 = repmat(Bu(1, :)', 1, Y, Z, 3).*(repmat(Bv(1, :), X, 1, Z, 3).*c(i0, j0, k0, :) + ...
                                  repmat(Bv(2, :), X, 1, Z, 3).*c(i0, j1, k0, :) + ... 
                                  repmat(Bv(3, :), X, 1, Z, 3).*c(i0, j2, k0, :) + ...
                                  repmat(Bv(4, :), X, 1, Z, 3).*c(i0, j3, k0, :));

T01 = repmat(Bu(2, :)', 1, Y, Z, 3).*(repmat(Bv(1, :), X, 1, Z, 3).*c(i1, j0, k0, :) + ...
                                  repmat(Bv(2, :), X, 1, Z, 3).*c(i1, j1, k0, :) + ... 
                                  repmat(Bv(3, :), X, 1, Z, 3).*c(i1, j2, k0, :) + ... 
                                  repmat(Bv(4, :), X, 1, Z, 3).*c(i1, j3, k0, :));

T02 = repmat(Bu(3, :)', 1, Y, Z, 3).*(repmat(Bv(1, :), X, 1, Z, 3).*c(i2, j0, k0, :) + ... 
                                  repmat(Bv(2, :), X, 1, Z, 3).*c(i2, j1, k0, :) + ... 
                                  repmat(Bv(3, :), X, 1, Z, 3).*c(i2, j2, k0, :) + ...
                                  repmat(Bv(4, :), X, 1, Z, 3).*c(i2, j3, k0, :));

T03 = repmat(Bu(4, :)', 1, Y, Z, 3).*(repmat(Bv(1, :), X, 1, Z, 3).*c(i3, j0, k0, :) + ...
                                  repmat(Bv(2, :), X, 1, Z, 3).*c(i3, j1, k0, :) + ... 
                                  repmat(Bv(3, :), X, 1, Z, 3).*c(i3, j2, k0, :) + ...
                                  repmat(Bv(4, :), X, 1, Z, 3).*c(i3, j3, k0, :));

matr = reshape(Bw(1, :), [1 1 Z]);
T0 = repmat(matr, X, X, 1, 3).*(T00 + T01 + T02 + T03);

clear T00 T01 T02 T03;

T10 = repmat(Bu(1, :)', 1, Y, Z, 3).*(repmat(Bv(1, :), X, 1, Z, 3).*c(i0, j0, k1, :) + ...
                                  repmat(Bv(2, :), X, 1, Z, 3).*c(i0, j1, k1, :) + ... 
                                  repmat(Bv(3, :), X, 1, Z, 3).*c(i0, j2, k1, :) + ...
                                  repmat(Bv(4, :), X, 1, Z, 3).*c(i0, j3, k1, :));

T11 = repmat(Bu(2, :)', 1, Y, Z, 3).*(repmat(Bv(1, :), X, 1, Z, 3).*c(i1, j0, k1, :) + ...
                                  repmat(Bv(2, :), X, 1, Z, 3).*c(i1, j1, k1, :) + ... 
                                  repmat(Bv(3, :), X, 1, Z, 3).*c(i1, j2, k1, :) + ... 
                                  repmat(Bv(4, :), X, 1, Z, 3).*c(i1, j3, k1, :));

T12 = repmat(Bu(3, :)', 1, Y, Z, 3).*(repmat(Bv(1, :), X, 1, Z, 3).*c(i2, j0, k1, :) + ... 
                                  repmat(Bv(2, :), X, 1, Z, 3).*c(i2, j1, k1, :) + ... 
                                  repmat(Bv(3, :), X, 1, Z, 3).*c(i2, j2, k1, :) + ...
                                  repmat(Bv(4, :), X, 1, Z, 3).*c(i2, j3, k1, :));

T13 = repmat(Bu(4, :)', 1, Y, Z, 3).*(repmat(Bv(1, :), X, 1, Z, 3).*c(i3, j0, k1, :) + ...
                                  repmat(Bv(2, :), X, 1, Z, 3).*c(i3, j1, k1, :) + ... 
                                  repmat(Bv(3, :), X, 1, Z, 3).*c(i3, j2, k1, :) + ...
                                  repmat(Bv(4, :), X, 1, Z, 3).*c(i3, j3, k1, :));

matr = reshape(Bw(2, :), [1 1 Z]);
T1 = repmat(matr, X, X, 1, 3).*(T10 + T11 + T12 + T13);

clear T10 T11 T12 T13;

T20 = repmat(Bu(1, :)', 1, Y, Z, 3).*(repmat(Bv(1, :), X, 1, Z, 3).*c(i0, j0, k2, :) + ...
                                  repmat(Bv(2, :), X, 1, Z, 3).*c(i0, j1, k2, :) + ... 
                                  repmat(Bv(3, :), X, 1, Z, 3).*c(i0, j2, k2, :) + ...
                                  repmat(Bv(4, :), X, 1, Z, 3).*c(i0, j3, k2, :));

T21 = repmat(Bu(2, :)', 1, Y, Z, 3).*(repmat(Bv(1, :), X, 1, Z, 3).*c(i1, j0, k2, :) + ...
                                  repmat(Bv(2, :), X, 1, Z, 3).*c(i1, j1, k2, :) + ... 
                                  repmat(Bv(3, :), X, 1, Z, 3).*c(i1, j2, k2, :) + ... 
                                  repmat(Bv(4, :), X, 1, Z, 3).*c(i1, j3, k2, :));

T22 = repmat(Bu(3, :)', 1, Y, Z, 3).*(repmat(Bv(1, :), X, 1, Z, 3).*c(i2, j0, k2, :) + ... 
                                  repmat(Bv(2, :), X, 1, Z, 3).*c(i2, j1, k2, :) + ... 
                                  repmat(Bv(3, :), X, 1, Z, 3).*c(i2, j2, k2, :) + ...
                                  repmat(Bv(4, :), X, 1, Z, 3).*c(i2, j3, k2, :));

T23 = repmat(Bu(4, :)', 1, Y, Z, 3).*(repmat(Bv(1, :), X, 1, Z, 3).*c(i3, j0, k2, :) + ...
                                  repmat(Bv(2, :), X, 1, Z, 3).*c(i3, j1, k2, :) + ... 
                                  repmat(Bv(3, :), X, 1, Z, 3).*c(i3, j2, k2, :) + ...
                                  repmat(Bv(4, :), X, 1, Z, 3).*c(i3, j3, k2, :));                              

matr = reshape(Bw(3, :), [1 1 Z]);
T2 = repmat(matr, X, X, 1, 3).*(T20 + T21 + T22 + T23);

clear T20 T21 T22 T23;

T30 = repmat(Bu(1, :)', 1, Y, Z, 3).*(repmat(Bv(1, :), X, 1, Z, 3).*c(i0, j0, k3, :) + ...
                                  repmat(Bv(2, :), X, 1, Z, 3).*c(i0, j1, k3, :) + ... 
                                  repmat(Bv(3, :), X, 1, Z, 3).*c(i0, j2, k3, :) + ...
                                  repmat(Bv(4, :), X, 1, Z, 3).*c(i0, j3, k3, :));

T31 = repmat(Bu(2, :)', 1, Y, Z, 3).*(repmat(Bv(1, :), X, 1, Z, 3).*c(i1, j0, k3, :) + ...
                                  repmat(Bv(2, :), X, 1, Z, 3).*c(i1, j1, k3, :) + ... 
                                  repmat(Bv(3, :), X, 1, Z, 3).*c(i1, j2, k3, :) + ... 
                                  repmat(Bv(4, :), X, 1, Z, 3).*c(i1, j3, k3, :));

T32 = repmat(Bu(3, :)', 1, Y, Z, 3).*(repmat(Bv(1, :), X, 1, Z, 3).*c(i2, j0, k3, :) + ... 
                                  repmat(Bv(2, :), X, 1, Z, 3).*c(i2, j1, k3, :) + ... 
                                  repmat(Bv(3, :), X, 1, Z, 3).*c(i2, j2, k3, :) + ...
                                  repmat(Bv(4, :), X, 1, Z, 3).*c(i2, j3, k3, :));

T33 = repmat(Bu(4, :)', 1, Y, Z, 3).*(repmat(Bv(1, :), X, 1, Z, 3).*c(i3, j0, k3, :) + ...
                                  repmat(Bv(2, :), X, 1, Z, 3).*c(i3, j1, k3, :) + ... 
                                  repmat(Bv(3, :), X, 1, Z, 3).*c(i3, j2, k3, :) + ...
                                  repmat(Bv(4, :), X, 1, Z, 3).*c(i3, j3, k3, :));

matr = reshape(Bw(4, :), [1 1 Z]);
T3 = repmat(matr, X, X, 1, 3).*(T30 + T31 + T32 + T33);

clear T30 T31 T32 T33;

gr = T0 + T1 + T2 + T3;

Спасибо.

Answer 1

Если у вас есть MATLAB R2016b или более поздняя версия, вы можете реализовать все эти операции напрямую, без каких-либо вызовов repmat:

T00 = repmat(Bu(1, :)', 1, Y, Z, 3).*(repmat(Bv(1, :), X, 1, Z, 3).*c(i0, j0, k0, :) + ...
                                      repmat(Bv(2, :), X, 1, Z, 3).*c(i0, j1, k0, :) + ... 
                                      repmat(Bv(3, :), X, 1, Z, 3).*c(i0, j2, k0, :) + ...
                                      repmat(Bv(4, :), X, 1, Z, 3).*c(i0, j3, k0, :));

совпадает с

T00 = Bu(1, :).' .* (Bv(1, :) .* c(i0, j0, k0, :) + ...
                     Bv(2, :) .* c(i0, j1, k0, :) + ... 
                     Bv(3, :) .* c(i0, j2, k0, :) + ...
                     Bv(4, :) .* c(i0, j3, k0, :));

Если у вас более старая версия MATLAB, вы можете использовать bsxfun для реализации каждого из этих операторов. Хотя это будет выглядеть намного страшнее.

---

Обратите внимание, что я заменил ' на .'. Первый - это сложное сопряженное транспонирование, которое, вероятно, не то, что вы намереваетесь использовать. Чтобы изменить ориентацию вектора, вы всегда должны предпочитать .'. Для вещественных матриц они одинаковы, но привыкание к правильному оператору предотвратит множество трудно обнаруживаемых ошибок, если вы когда-нибудь будете работать со сложными данными в будущем.

Answer 2

Как насчет прямой поэлементной работы? С последними версиями Matlab он делает все размеры для вас:

B = repmat(A, 240, 1, 180, 3);  % Elapsed time is 0.264423 seconds.
C = A.*ones(240,1,180,3);       % Elapsed time is 0.027216 seconds.

Отказ от ответственности : время было измерено с помощью R2019a на ноутбуке с tic и toc, поэтому вы, вероятно, получите разные значения и даже другое соотношение. Однако я думаю, что разница во времени достаточно значительна, чтобы попробовать.