ssd算法 matlab,Shear-Warp算法Matlab代码解读

Render函数

function render_image = render(volume,options)

此间省略部分。

%% Calculate the Shear and Warp

Matrices

if(ndims(data.Volume)==2)

sizes=[size(data.Volume) 1];

else

sizes=size(data.Volume);

end

求出Mshear，Mwarp及其逆矩阵，c是六种视角面(立方体只能看到三个面，六种可能性)。

[data.Mshear,data.Mwarp2D,data.c]=makeShearWarpMatrix(data.Mview,sizes);

data.Mwarp2Dinv=inv(double(data.Mwarp2D));

data.Mshearinv=inv(data.Mshear);

%% Store Volume sizes

data.Iin_sizex=size(data.Volume,1); data.Iin_sizey=size(data.Volume,2); data.Iin_sizez=size(data.Volume,3);

%% Create

Shear (intimidate) buffer

为什么乘以

？

data.Ibuffer_sizex=ceil(1.7321*max(size(data.Volume))+1);

data.Ibuffer_sizey=data.Ibuffer_sizex;

switch data.RenderType

case {'mip'}

data.Ibuffer=zeros([data.Ibuffer_sizex data.Ibuffer_sizey])+data.imin;

case {'bw'}

data.Ibuffer=zeros([data.Ibuffer_sizex data.Ibuffer_sizey]);

otherwise

data.Ibuffer=zeros([data.Ibuffer_sizex data.Ibuffer_sizey 3]);

end

这个的效果是减小了不透明度，为什么？

%% Adjust alpha table by voxel length because of

rotation and volume size

lengthcor=sqrt(1+data.Mshearinv(1,3)^2+data.Mshearinv(2,3)^2)*mean(size(data.Volume))/100;

data.AlphaTable=1 - (1-data.AlphaTable).^(1/lengthcor);

data.AlphaTable(data.AlphaTable<0)=0; data.AlphaTable(data.AlphaTable>1)=1;

%% Shading type -> Phong

values

switch lower(data.ShadingMaterial)

case {'shiny'}

data.material=[0.7, 0.6, 0.9, 15];

case {'dull'}

data.material=[0.7, 0.8, 0.0, 10];

case {'metal'}

data.material=[0.7, 0.3, 1.0, 20];

otherwise

data.material=[0.7, 0.6, 0.9, 20];

end

%% Normalize Light and Viewer

vectors

data.LightVector=[data.LightVector(:);0]; data.LightVector=data.LightVector./sqrt(sum(data.LightVector(1:3).^2));

data.ViewerVector=[data.ViewerVector(:);0]; data.ViewerVector=data.ViewerVector./sqrt(sum(data.ViewerVector(1:3).^2));

%% Shear Rendering

data = shear(data);

data = warp(data);

render_image = data.Iout;

makeShearWarpMatrix

function [Mshear,Mwarp2D,c]=makeShearWarpMatrix(Mview,sizes)

% Function MAKESHEARWARPMATRIX splits a View Matrix

in to

% a shear matrix and warp matrix, for efficient 3D

volume rendering.

%

%

[Mshear,Mwarp2D,c]=makeShearWarpMatrix(Mview,sizes)

%

% inputs,

% Mview: The 4x4 viewing matrix

% sizes: The sizes of the volume which will be

rendered

%

% outputs,

% Mshear: The shear matrix

% Mwarp2D: The warp matrix

% c: The principal viewing axis

1..6

%

% example,

%

% Mview=makeViewMatrix([45 45 0],[0.5 0.5 0.5],[0 0

0]);

% sizes=[512 512];

%

[Mshear,Mwarp2D,c]=makeShearWarpMatrix(Mview,sizes)

%

% Function is written by D.Kroon University of Twente

(October 2008)

% Find the principal viewing axis

Vo=[Mview(1,2)*Mview(2,3) - Mview(2,2)*Mview(1,3);

Mview(2,1)*Mview(1,3) - Mview(1,1)*Mview(2,3);

Mview(1,1)*Mview(2,2) - Mview(2,1)*Mview(1,2)];

[maxv,c]=max(abs(Vo));

% Choose the corresponding Permutation matrix

P

switch(c)

case 1, %yzx

P=[0 1

0 0; 0 0

1 0; 1 0

0 0; 0 0

0 1;];

case 2, %

zxy

P=[0 0

1 0; 1 0

0 0; 0 1

0 0; 0 0

0 1;];

case 3, %

xyz

P=[1 0

0 0; 0 1

0 0; 0 0

1 0; 0 0

0 1;];

end

% Compute the

permuted view matrix from Mview and P

矩阵除法在MATLAB中，有两种矩阵除法运算：\和/，分别表示左除和右除。如果A矩阵是非奇异方阵，则A\B和B/A运算可以实现。A\B等效于A的逆左乘B矩阵，也就是inv(A)*B，而B/A等效于A矩阵的逆右乘B矩阵，也就是B*inv(A)。对于含有标量的运算，两种除法运算的结果相同。对于矩阵来说，左除和右除表示两种不同的除数矩阵和被除数矩阵的关系，一般A\B≠B/A。

Mview_p=Mview/P;

% 180 degrees rotate detection

if(Mview_p(3,3)<0), c=c+3; end

% Compute the shear coeficients from the permuted

view matrix

Si = (Mview_p(2,2)* Mview_p(1,3) - Mview_p(1,2)* Mview_p(2,3)) / (Mview_p(1,1)* Mview_p(2,2) - Mview_p(2,1)* Mview_p(1,2));

Sj = (Mview_p(1,1)* Mview_p(2,3) - Mview_p(2,1)* Mview_p(1,3)) / (Mview_p(1,1)* Mview_p(2,2) - Mview_p(2,1)* Mview_p(1,2));

% Compute the translation between the orgins of

standard object coordinates

% and intermdiate image coordinates

if((c==1)||(c==4)), kmax=sizes(1)-1; end

if((c==2)||(c==5)), kmax=sizes(2)-1; end

if((c==3)||(c==6)), kmax=sizes(3)-1; end

if ((Si>=0)&&(Sj>=0)),

Ti = 0; Tj = 0; end

if ((Si>=0)&&(Sj<0)),

Ti = 0; Tj = -Sj*kmax; end

if ((Si<0)&&(Sj>=0)),

Ti = -Si*kmax; Tj = 0; end

if ((Si<0)&&(Sj<0)),

Ti = -Si*kmax; Tj = -Sj*kmax; end

% Compute the shear matrix

Mshear=[1 0

Si Ti;

0 1

Sj Tj;

0 0

1 0;

0 0

0 1];

% Compute the 2Dwarp matrix

Mwarp2D=[Mview_p(1,1) Mview_p(1,2) (Mview_p(1,4)-Ti*Mview_p(1,1)-Tj*Mview_p(1,2));

Mview_p(2,1) Mview_p(2,2) (Mview_p(2,4)-Ti*Mview_p(2,1)-Tj*Mview_p(2,2));

0 0

1 ];

Shear

%% Shearwarp functions

function data=shear(data)

switch (data.c)

case 1

for z=0:(data.Iin_sizex-1);

% Offset calculation

xd=(-data.Ibuffer_sizex/2)+data.Mshearinv(1,3)*(z-data.Iin_sizex/2)+data.Iin_sizey/2;

yd=(-data.Ibuffer_sizey/2)+data.Mshearinv(2,3)*(z-data.Iin_sizex/2)+data.Iin_sizez/2;

xdfloor=floor(xd);

ydfloor=floor(yd);

�lculate the coordinates on which a image slice starts

and

%ends in the temporary shear image (buffer)

pystart=-ydfloor;

if(pystart<0),

pystart=0;

end

pyend=data.Iin_sizez-ydfloor;

if(pyend>data.Ibuffer_sizey),

pyend=data.Ibuffer_sizey;

end

pxstart=-xdfloor;

if(pxstart<0),

pxstart=0;

end

pxend=data.Iin_sizey-xdfloor;

if(pxend>data.Ibuffer_sizex),

pxend=data.Ibuffer_sizex;

end

data.py=(pystart+1:pyend-1);

data.px=(pxstart+1:pxend-1);

if(isempty(data.px)),

data.px=pxstart+1;

end

if(isempty(data.py)),

data.py=pystart+1;

end

% Determine x and y coordinates of pixel(s) which will be come

current pixel

yBas=data.py+ydfloor;

xBas=data.px+xdfloor;

switch (data.ShearInterp)

case {'bilinear'}

xBas1=xBas+1;

xBas1(end)=xBas1(end)-1;

yBas1=yBas+1;

yBas1(end)=yBas1(end)-1;

% Linear interpolation constants (percentages)

内插时的距离计算

xCom=xd-floor(xd);

yCom=yd-floor(yd);

perc=[(1-xCom)*(1-yCom)

(1-xCom)*yCom

xCom*(1-yCom)

xCom*yCom];

if(isempty(data.VolumeX))

% Get the intensities

亮度，squeeze Remove singleton dimensions.为什么不是相邻切片间内插？

if(data.Iin_sizez>1)

slice=double(squeeze(data.Volume(z+1,:,:)));

else

slice=double(data.Volume(z+1,:))';

end

intensity_xyz1=slice(xBas,

yBas);

intensity_xyz2=slice(xBas,

yBas1);

intensity_xyz3=slice(xBas1,

yBas);

intensity_xyz4=slice(xBas1,

yBas1);

else

slice=double(data.VolumeX(:,

:,z+1));

intensity_xyz1=slice(xBas,

yBas);

intensity_xyz2=slice(xBas,

yBas1);

intensity_xyz3=slice(xBas1,

yBas);

intensity_xyz4=slice(xBas1,

yBas1);

end

% Calculate the interpolated intensity

data.intensity_loc=(intensity_xyz1*perc(1)+intensity_xyz2*perc(2)+intensity_xyz3*perc(3)+intensity_xyz4*perc(4));

otherwise

if(isempty(data.VolumeX))

data.intensity_loc=double(squeeze(data.Volume(z+1,xBas,

yBas)));

else

data.intensity_loc=double(data.VolumeX(xBas,

yBas,z+1));

end

end

% Update the shear image buffer

switch (data.RenderType)

case {'mip'}

data=updatebuffer_MIP(data);

case {'color'}

data=updatebuffer_COLOR(data);

case {'bw'}

data=updatebuffer_BW(data);

case {'shaded'}

data=returnnormal(z+1,xBas,

yBas,data);

data=updatebuffer_SHADED(data);

end

end

case其他六种情况

updatebuffer_BW

function data=updatebuffer_BW(data)

% Calculate index in alpha transparency look up

table

if(data.imin~=0), data.intensity_loc=data.intensity_loc-data.imin; end

if(data.imaxmin~=1), data.intensity_loc=data.intensity_loc./data.imaxmin; end

betaA=(length(data.AlphaTable)-1);

通过亮度查出不透明度，越亮说明越不透明？

indexAlpha=round(data.intensity_loc*betaA)+1;

% calculate current alphaimage

alphaimage=data.AlphaTable(indexAlpha);

% 2D volume fix because alphaimage becomes a row instead of

column

if(data.Iin_sizez==1), alphaimage=reshape(alphaimage,size(data.Ibuffer(data.px,data.py))); end

alphaimage_inv=(1-alphaimage);

% Update the current pixel in the shear image buffer

Over操作，

data.Ibuffer(data.px,data.py)=alphaimage_inv.*data.Ibuffer(data.px,data.py)+alphaimage.*data.intensity_loc;

Warp

function data=warp(data)

% This function warp, will warp the shear rendered buffer

image

% Make Affine matrix

M=zeros(3,3);

M(1,1)=data.Mwarp2Dinv(1,1);

M(2,1)=data.Mwarp2Dinv(2,1);

M(1,2)=data.Mwarp2Dinv(1,2);

M(2,2)=data.Mwarp2Dinv(2,2);

M(1,3)=data.Mwarp2Dinv(1,3)+data.Mshearinv(1,4);

M(2,3)=data.Mwarp2Dinv(2,3)+data.Mshearinv(2,4);

% Perform the affine transformation

switch(data.WarpInterp)

case 'nearest', wi=5;

case 'bicubic', wi=3;

case 'bilinear', wi=1;

otherwise, wi=1;

end

data.Iout=affine_transform_2d_double(data.Ibuffer,M,wi,data.ImageSize);

解释：Affine transformation function (Rotation, Translation, Resize)

This function transforms a volume with a 3x3 transformation matrix