Объемная проекция в оттенках серого 3d на плоскость 2D

Question

У меня есть объемная шкала серого, соответствующая данным УЗИ.В Matlab этот трехмерный том является просто трехмерной матрицей MxNxP .Интересующая меня структура ориентирована не вдоль оси z , а вдоль уже известной локальной системы координат ( x'y'z ').До этого момента у меня есть что-то вроде рисунка, показанного ниже, с изображением оригинальной ( xyz ) и локальной системы координат ( x'y'z '):

Я хочу получить двумерную проекцию этого объема (т. Е. Изображение) через конкретную плоскость в локальной системе координат, скажем, в z '= z0 .Как я могу это сделать?

Если бы объем был ориентирован вдоль оси z , эта проекция могла бы быть легко достигнута.т. е. если объем в Matlab равен V , то:

projection = sum(V,3);

, таким образом, проекция может быть вычислена как сумма по 3-му измерению массива.Однако с изменением ориентации проблема становится более сложной.

Я смотрел на радоновое преобразование (2D, которое применяется только к двумерным изображениям, а не на объемы), а также рассматривал ортографические проекции, но вя не знаю, что делать!

Спасибо за любой совет!

Answer 1

Новая попытка решения:

После обучения http://blogs.mathworks.com/steve/2006/08/17/spatial-transformations-three-dimensional-rotation/ и внесения небольших изменений я могу кое-что вам помочь. Имейте в виду, у меня мало или нет опыта работы с объемными данными в MATLAB, поэтому реализация довольно хакерская.

В приведенном ниже коде я использую tformarray (), чтобы вращать структуру в пространстве. Сначала данные центрируются, а затем поворачиваются с использованием вращениемmat3D для создания пространственного преобразования, прежде чем данные возвращаются в исходное положение.

Поскольку я никогда раньше не использовал tformarray, я управлял точками данных, выходящими за пределы определенной области после поворота, просто заполняя матрицу данных (NxMxP) нулями вокруг. Если кто-нибудь знает лучший способ, пожалуйста, сообщите нам:)

код:

    %Synthetic dataset, 25x50x25
blob = flow();

%Pad to allow for rotations in space. Bad solution, 
%something better might be possible to better understanding
%of tformarray()
blob = padarray(blob,size(blob));

f1 = figure(1);clf;
s1=subplot(1,2,1);
p = patch(isosurface(blob,1));
set(p, 'FaceColor', 'red', 'EdgeColor', 'none');
daspect([1 1 1]);
view([1 1 1])
camlight
lighting gouraud

%Calculate center
blob_center = (size(blob) + 1) / 2;

%Translate to origin transformation
T1 = [1 0 0 0
    0 1 0 0
    0 0 1 0
    -blob_center 1];

%Rotation around [0 0 1]
rot = -pi/3;
Rot = rotationmat3D(rot,[0 1 1]);
T2 = [ 1  0  0   0
       0  1  0   0
       0  0  1   0
       0  0  0   1];
T2(1:3,1:3) = Rot;   

%Translation back
T3 = [1 0 0 0
    0 1 0 0
    0 0 1 0
    blob_center 1];

%Total transform
T = T1 * T2 * T3;

%See http://blogs.mathworks.com/steve/2006/08/17/spatial-transformations-three-dimensional-rotation/
tform = maketform('affine', T);
R = makeresampler('linear', 'fill');
TDIMS_A = [1 2 3];
TDIMS_B = [1 2 3];
TSIZE_B = size(blob);
TMAP_B = [];
F = 0;
blob2 = ...
tformarray(blob, tform, R, TDIMS_A, TDIMS_B, TSIZE_B, TMAP_B, F);

s2=subplot(1,2,2);
p2 = patch(isosurface(blob2,1));
set(p2, 'FaceColor', 'red', 'EdgeColor', 'none');
daspect([1 1 1]);
view([1 1 1])
camlight
lighting gouraud

Произвольная визуализация, приведенная ниже, просто для подтверждения того, что данные повернуты, как и ожидалось, на графике замкнутой поверхности, когда данные приняли значение «1». С blob2 вы должны знать, что сможете проецировать, используя простые суммы.

figure(2)
subplot(1,2,1);imagesc(sum(blob,3));
subplot(1,2,2);imagesc(sum(blob2,3));

Answer 2

Предполагая, что у вас есть доступ к базису координат R = [x 'y' z '], и что эти векторы ортонормированы, вы можете просто извлечь представление в этом базисе, умножив ваши данные на матрицу 3x3 R, где x ', y', z '- векторы столбцов.

С данными, хранящимися в D (Nx3), вы можете получить представление с помощью R, умножив на него: Dmarked = D * R;

и теперь D = Dmarked * inv (R), так что возвращаться туда и обратно очень трудно.

Следующий код может помочь увидеть трансформацию. Здесь я создаю синтетический набор данных, поворачиваю его, а затем поворачиваю обратно. В этом случае сумма (DR (:, 3)) будет вашей суммой по z '

%#Create synthetic dataset
N1 = 250;
r1 = 1;
dr1 = 0.1;
dz1 = 0;
mu1 = [0;0];
Sigma1 = eye(2);
theta1 = 0 + (2*pi).*rand(N1,1);
rRand1 = normrnd(r1,dr1,1,N1);
rZ1 = rand(N1,1)*dz1+1;
D = [([rZ1*0 rZ1*0] + repmat(rRand1',1,2)).*[sin(theta1) cos(theta1)] rZ1];

%Create roation matrix
rot = pi/8;
R = rotationmat3D(rot,[0 1 0]);
% R =       0.9239          0       0.3827
%           0               1.0000  0
%           -0.3827         0       0.9239

Rinv = inv(R);

%Rotate data
DR = D*R;

%#Visaulize data
f1 = figure(1);clf
subplot(1,3,1);
plot3(DR(:,1),DR(:,2),DR(:,3),'.');title('Your data')
subplot(1,3,2);
plot3(DR*Rinv(:,1),DR*Rinv(:,2),DR*Rinv(:,3),'.r');
view([0.5 0.5 0.2]);title('Representation using your [xmarked ymarked zmarked]');
subplot(1,3,3);
plot3(D(:,1),D(:,2),D(:,3),'.');
view([0.5 0.5 0.2]);title('Original data before rotation');

Answer 3

Если у вас есть два нормализованных 3x1 вектора x2 и y2, соответствующих вашей локальной системе координат (x 'и y').

Затем для позиции P, его локальная координата будет xP=P'x2 и yP=P'*y2.

Так что вы можете попытаться спроецировать свой объем, используя accumarray:

[x y z]=ndgrid(1:M,1:N,1:P);
posP=[x(:) y(:) z(:)];
xP=round(posP*x2);
yP=round(posP*y2);
xP=xP+min(xP(:))+1;
yP=yP+min(yP(:))+1;
V2=accumarray([xP(:),yP(:)],V(:));

Если вы предоставите свои данные, я протестируюэто.