-
Notifications
You must be signed in to change notification settings - Fork 5
/
mesh2model.m
107 lines (98 loc) · 4.56 KB
/
mesh2model.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
function phanimg = mesh2model(tissueprop,linv,xpts,ypts,zpts)
% Create 3D phantom with respiratory motion
%
% INPUTS:
% tissueprop [struct] -> tissue property generated by tissueproperty.m
% linv [nt x 3] -> 3D respiratory motion curves generated by genresp.m
% xpts [nx x 1] -> x resolution
% ypts [ny x 1] -> y resolution
% zpts [nz x 1] -> z resolution
%
% OUTPUT:
% phanimg [nx, ny, nz, nt] -> 3D phantom with respiratory motion
%
% -----------------------------------------------------------------------------------------
% Realistic 4D abdominal phantom for magnetic resonance imaging
% Wei-Ching Lo
% wxl317@case.edu
% Case Western Reserve University
% April 2018
% -----------------------------------------------------------------------------------------
% mtx = numel(xpts);
npar = numel(zpts);
nt = size(linv,1);
% scale and shift 3D model according to non-rigid motion
OUTPUTgridvalue = zeros(length(xpts),length(ypts),npar,nt);
for onum = 1:length(tissueprop)
selstl = fullfile('stlmodels_release',[tissueprop(onum).name '.stl']);
if strcmp(tissueprop(onum).name,'fat')
fat = onum;
elseif strcmp(tissueprop(onum).name,'skin')
skin = onum;
else
disp(['Start deformation and voxelisation: ' tissueprop(onum).name]);
origmesh = READ_stl(selstl);
organidx = false(length(xpts),length(ypts),npar,nt);
for scalev = 1:size(linv,1)
meshXYZ = origmesh;
% simple non-rigid transformation for abdominal phantom
meshZind = squeeze(meshXYZ(:,3,:) > 0);
meshZtemp = squeeze(meshXYZ(:,3,:));
meshZtemp(meshZind) = meshZtemp(meshZind)*linv(scalev,1); % scale along Z (SI)
meshXYZ(:,3,:) = meshZtemp+(linv(scalev,1)-1)*30; % shift along Z (SI)
meshYind = squeeze(meshXYZ(:,2,:) > min(ypts)+20);
meshYtemp = squeeze(meshXYZ(:,2,:))-(min(ypts)+20);
meshYtemp(meshYind) = meshYtemp(meshYind)./linv(scalev,2)+(linv(scalev,2)-1)*50; % scale along Y (AP)
meshXYZ(:,2,:) = meshYtemp+(min(ypts)+20);
meshXind = squeeze(meshXYZ(:,1,:) > min(xpts) & meshXYZ(:,1,:) < max(xpts));
meshXtemp = squeeze(meshXYZ(:,1,:));
meshXtemp(meshXind) = meshXtemp(meshXind)./linv(scalev,3); % scale along X (LR)
meshXYZ(:,1,:) = meshXtemp;
% voxelization
organidx(:,:,:,scalev) = VOXELISE(xpts,ypts,zpts,meshXYZ);
end
OUTPUTgridvalue(organidx) = onum;
end
end
% Define ellipse center for fat and skin mask
[~, xfirst] = max(sum(OUTPUTgridvalue,2) ~= 0,[],1);
minx = squeeze(xfirst);
[~, fillstart] = max(minx == 1,[],1);
xfill = find(fillstart > 1);
for i = 1:numel(xfill)
minx(fillstart(xfill(i)):end,xfill(i)) = minx(fillstart(xfill(i))-1,xfill(i));
end
[~, xlast] = max(flip(sum(OUTPUTgridvalue,2) ~= 0,1),[],1);
maxx = length(xpts)-squeeze(xlast)+1;
[~, fillstart] = max(maxx == length(xpts),[],1);
xfill = find(fillstart > 1);
for i = 1:numel(xfill)
maxx(fillstart(xfill(i)):end,xfill(i)) = maxx(fillstart(xfill(i))-1,xfill(i));
end
[~, yfirst] = max(sum(OUTPUTgridvalue,1) ~= 0,[],2);
miny = squeeze(yfirst);
[~, fillstart] = max(miny == 1,[],1);
yfill = find(fillstart > 1);
for i = 1:numel(yfill)
miny(fillstart(yfill(i)):end,yfill(i)) = miny(fillstart(yfill(i))-1,yfill(i));
end
[~, ylast] = max(flip(sum(OUTPUTgridvalue,1) ~= 0,2),[],2);
maxy = length(ypts)-squeeze(ylast)+1;
[~, fillstart] = max(maxy == length(ypts),[],1);
yfill = find(fillstart > 1);
for i = 1:numel(yfill)
maxy(fillstart(yfill(i)):end,yfill(i)) = maxy(fillstart(yfill(i))-1,yfill(i));
end
for b = 1:numel(zpts)*nt
cx = (minx(b)+maxx(b))/2; % Ellipse center point
cy = (miny(b)+maxy(b))/2;
xr_sq = ((maxy(b)-miny(b))/2+50/(256/length(ypts)*2))^2; % Ellipse radii squared
yr_sq = ((maxx(b)-minx(b))/2+25/(256/length(xpts)*2))^2;
[X, Y] = meshgrid(1:length(ypts),1:length(xpts));
ellipse_mask = (yr_sq * (X - cy -6+(256/length(xpts))) .^ 2 + xr_sq * (Y - cx) .^ 2 <= (xr_sq*yr_sq));
ellipse_skinmask = ((yr_sq*0.95) * (X - cy -6+(256/length(ypts))) .^ 2 + (xr_sq*0.95) * (Y - cx) .^ 2 <= (xr_sq*yr_sq));
fatmask = (ellipse_mask-logical(OUTPUTgridvalue(:,:,b)) == 1);
skinmask = (ellipse_skinmask-(fatmask+logical(OUTPUTgridvalue(:,:,b))) == 1);
OUTPUTgridvalue(:,:,b) = OUTPUTgridvalue(:,:,b)+fatmask*fat+skinmask*skin;
end
phanimg = permute(flip(flip(OUTPUTgridvalue,2),1),[2 1 3 4]);