【0048】
リスト1
1 % 複数のCATスキャンを読み込み、それらのスキャンを積み重ねることによって3D画像を形成する
2 % 複数のスキャンをカラーで結合し、3D画像を2Dプロット上にレンダリングする
3 % イニシャル画像をその他から減算して、望ましくないものの一部を除去する
4 % 例えば肋骨などを返す
5
6 clear
7 close all
8
9 % コントラスト設定
10 bias = -1050;
11 gain = 0.3;
12
13 % 減算する第1の画像の量を選択する(0-1)
14 subtraction_fraction = 0;
15
16 ROTXY = 1;
17 ROTXZ = 2;
18 ROTYZ = 3;
19
20 % このラン(run)での回転の方向を選択する
21 rot_axis = ROTXY;
22
23 % 3D画像を2Dに投影するのに使用される統計に関する重みを選択する
24 c1 = 0;
25 c2 = 1;
26 c3 = 0;
27 c4 = 0;
28
29 % メジアンフィルタ(3次元フィルタ)の幅を選択する
30 NMED = 3;
31
32 % スキャンを読み込む
33
34 load sample_patient_1_a3_no_filt.mat
35 scan1 = a3;
36 load sample_patient_2_a3_no_filt.mat
37 scan2 = a3;
38 load sample_patient_3_a3_no_filt.mat
39 scan3 = a3;
40 load sample_patient_4_a3_no_filt.mat
41 scan4 = a3;
42
43 scan21=scan2-subtraction_fraction*scan1;
44 scan31=scan3-subtraction_fraction*scan1;
45 scan41=scan4-subtraction_fraction*scan1;
46
47 [Nx,Ny,Nz] = size(scan21);
48
49 % バイアスし、フィルタリングし、そして、色に割り当てる
50 red = medfilt3(scan21+bias,NMED);
51 grn = medfilt3(scan31+bias,NMED);
52 blu = medfilt3(scan41+bias,NMED);
53
54 % ゼロ未満のものをクリッピング
55 red = max(red,0);
56 grn = max(grn,0);
57 blu = max(blu,0);
58
59 % アレイを初期化
60 red1=zeros(size(red),'single');
61 grn1=zeros(size(grn),'single');
62 blu1=zeros(size(blu),'single');
63
64 [Nx,Ny,Nz] = size(red);
65
66 if rot_axis == ROTXZ
67 front_lighting = ones(Nx,1) * (1:Nz).^2;
68 elseif rot_axis == ROTYZ
69 front_lighting = ones(Ny,1) * (1:Nz).^2;
70 elseif rot_axis == ROTXY
71 front_lighting = [zeros(1,120) (Nx-120-50:-1:1) zeros(1,50)]' * ones(1,Ny);
72 front_lighting = front_lighting.^4 ./ Nx/Nx/Nx;
73 end
74
75 ANG_LIMIT = 90;
76
77 ANG_STEP = round( ANG_LIMIT/45);
78
79 ang_array = [(0:-ANG_STEP:-ANG_LIMIT),...
80 (-ANG_LIMIT+1:ANG_STEP:ANG_LIMIT),(ANG_LIMIT-1:-ANG_STEP:1)];
81
82 img_count = 0;
83
84 % 多数の観察角度にわたってループさせて、3D画像が以下のように観察され得るようにする
85 % 多数の2D画像にてスクリーン上に投影される
86
87 for ang = ang_array
88
89 for n = 1 : Ny*(rot_axis==ROTXZ) + Nx*(rot_axis==ROTYZ) + Nz*(rot_axis==ROTXY)
90
91 if rot_axis == ROTXZ
92 red_temp = squeeze( red(:,n,:) );
93 grn_temp = squeeze( grn(:,n,:) );
94 blu_temp = squeeze( blu(:,n,:) );
95 elseif rot_axis == ROTYZ
96 red_temp = squeeze( red(n,:,:) );
97 grn_temp = squeeze( grn(n,:,:) );
98 blu_temp = squeeze( blu(n,:,:) );
99 elseif rot_axis == ROTXY
100 red_temp = squeeze( red(:,:,n) );
101 grn_temp = squeeze( grn(:,:,n) );
102 blu_temp = squeeze( blu(:,:,n) );
103 end
104
105 % 観察のために3つの軸のいずれかでスクリーン上で画像を回転させる
106 red_tempr = imrotate(red_temp,ang,'bilinear','crop');
107 grn_tempr = imrotate(grn_temp,ang,'bilinear','crop');
108 blu_tempr = imrotate(blu_temp,ang,'bilinear','crop');
109
110 if rot_axis == ROTXZ
111 red1(:,n,:) = red_tempr .* front_lighting;
112 grn1(:,n,:) = grn_tempr .* front_lighting;
113 blu1(:,n,:) = blu_tempr .* front_lighting;
114 elseif rot_axis == ROTYZ
115 red1(n,:,:) = red_tempr .* front_lighting;
116 grn1(n,:,:) = grn_tempr .* front_lighting;
117 blu1(n,:,:) = blu_tempr .* front_lighting;
118 elseif rot_axis == ROTXY
119 red1(:,:,n) = red_tempr .* front_lighting;
120 grn1(:,:,n) = grn_tempr .* front_lighting;
121 blu1(:,:,n) = blu_tempr .* front_lighting;
122 end
123 end
124
125 if rot_axis == ROTXZ || rot_axis == ROTYZ
126
127 % 統計を用いて画像を3Dから2Dへと投影
128 red_im =
squeeze(c1*mean(red1,3)+c2*max(red1,[],3)+c3*std(red1,[],3));
129 grn_im =
squeeze(c1*mean(grn1,3)+c2*max(grn1,[],3)+c3*std(grn1,[],3));
130 blu_im =
squeeze(c1*mean(blu1,3)+c2*max(blu1,[],3)+c3*std(blu1,[],3));
131
132 col_im(:,:,1) = red_im;
133 col_im(:,:,2) = grn_im;
134 col_im(:,:,3) = blu_im;
135
136 if rot_axis == ROTXZ
137 image(max(0,min(1,gain*col_im/mean(red_im(:)))))
138 else
139 image(max(0,min(1,gain*col_im/mean(red_im(:)))))
140 end
141
142 elseif rot_axis == ROTXY
143
144 front_lighting1 = zeros(size(red1));
145 for z = 1 : Nz
146 front_lighting1(:,:,z) = front_lighting;
147 end
148
149 % 統計を用いて画像を3Dから2Dへと投影
150 red_im =
squeeze(c1*mean(red1,1)+c2*max(red1,[],1)+c3*std(red1,[],1))';
151 grn_im =
squeeze(c1*mean(grn1,1)+c2*max(grn1,[],1)+c3*std(grn1,[],1))';
152 blu_im =
squeeze(c1*mean(blu1,1)+c2*max(blu1,[],1)+c3*std(blu1,[],1))';
153 end
154
155 % 画像をMatlab RGBアレイに入れる 1=red, 2=grn, 3=blu in 3rd dimension
156 % imresizeを用いてスクリーン上で拡大する
157 col_im(:,:,1) = imresize(red_im,850/size(red_im,1),'bicubic');
158 col_im(:,:,2) = imresize(grn_im,850/size(red_im,1),'bicubic');
159 col_im(:,:,3) = imresize(blu_im,850/size(red_im,1),'bicubic');
160
161 image(max(0,min(1,gain*col_im/mean(red_im(:)))))
162 truesize
163 drawnow
164
165 img_count = img_count + 1;
166
167 Mov(img_count) = getframe;
168
169 end