📄 avw_shrinkwrap.html
字号:
0189 <span class="comment">% Is the mean distance reasonable?</span>0190 <span class="keyword">if</span> mean(R) < 0.5,0191 error(<span class="string">'...mean distance < 0.5 voxel!\n'</span>);0192 <span class="keyword">end</span>0193 0194 <span class="comment">% MDifVal is the mean of the absolute difference</span>0195 <span class="comment">% between the vertex intensity and the fit intensity</span>0196 MDifVal = abs(intensityAtRMean(2) - fit.val);0197 0198 <span class="comment">% Is the mean difference within the tolerance range?</span>0199 <span class="keyword">if</span> MDifVal < fit.tol,0200 fprintf(<span class="string">'...mean intensity difference < tolerance (%5.2f +/- %5.2f)\n'</span>,<span class="keyword">...</span>0201 fit.val,fit.tol);0202 <span class="keyword">break</span>;0203 <span class="keyword">else</span>0204 fprintf(<span class="string">'...mean intensity difference > tolerance (%5.2f +/- %5.2f)\n'</span>,<span class="keyword">...</span>0205 fit.val,fit.tol);0206 <span class="keyword">end</span>0207 0208 <span class="comment">% How much has the intensity values changed?</span>0209 <span class="keyword">if</span> (i > 1) & intensityAtRMean(2) > 0,0210 <span class="keyword">if</span> intensityAtRMean(2) - intensityAtRMean(1) < fit.change,0211 fprintf(<span class="string">'...no significant intensity change (< %5.2f) in this iteration\n'</span>,<span class="keyword">...</span>0212 fit.change);0213 Fminima = Fminima + 1;0214 <span class="keyword">if</span> Fminima >= 5,0215 fprintf(<span class="string">'...no significant intensity change in last 5 iterations\n'</span>);0216 <span class="keyword">break</span>;0217 <span class="keyword">end</span>0218 <span class="keyword">else</span>0219 Fminima = 0;0220 <span class="keyword">end</span>0221 <span class="keyword">end</span>0222 0223 <span class="comment">% Ensure that iterations begin when MDifVal is truly initialised</span>0224 <span class="keyword">if</span> isnan(MDifVal),0225 i = 1;0226 <span class="keyword">else</span>,0227 i = i + 1;0228 <span class="keyword">end</span>0229 0230 <span class="keyword">end</span>0231 0232 FV = mesh_smooth(FV,origin,fit.vattr);0233 0234 <span class="comment">% Remove large edges matrix from FV</span>0235 Edges = FV.edge;0236 FV = struct(<span class="string">'vertices'</span>,FV.vertices,<span class="string">'faces'</span>,FV.faces);0237 0238 <span class="comment">% Now center the output vertices at 0,0,0 by subtracting</span>0239 <span class="comment">% the volume centroid</span>0240 FV.vertices(:,1) = FV.vertices(:,1) - origin(1);0241 FV.vertices(:,2) = FV.vertices(:,2) - origin(2);0242 FV.vertices(:,3) = FV.vertices(:,3) - origin(3);0243 0244 <span class="comment">% Scale the vertices by the pixdim values, after</span>0245 <span class="comment">% converting them from mm to meters</span>0246 xpixdim = double(avw.hdr.dime.pixdim(2)) / 1000;0247 ypixdim = double(avw.hdr.dime.pixdim(3)) / 1000;0248 zpixdim = double(avw.hdr.dime.pixdim(4)) / 1000;0249 0250 FV.vertices(:,1) = FV.vertices(:,1) .* xpixdim;0251 FV.vertices(:,2) = FV.vertices(:,2) .* ypixdim;0252 FV.vertices(:,3) = FV.vertices(:,3) .* zpixdim;0253 0254 0255 t=toc; fprintf(<span class="string">'...done (%5.2f sec).\n\n'</span>,t);0256 0257 <span class="keyword">return</span>0258 0259 0260 0261 0262 0263 0264 0265 0266 0267 0268 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>0269 <span class="comment">% Subfunctions...</span>0270 0271 0272 0273 0274 0275 0276 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>0277 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>0278 <a name="_sub1" href="#_subfunctions" class="code">function [fit] = estimate_scalp(FV,vol,origin,fit),</a>0279 0280 xo = origin(1); yo = origin(2); zo = origin(3);0281 0282 Nvert = size(FV.vertices,1);0283 0284 <span class="comment">% Estimate the scalp intensity, using 10% of vertices</span>0285 N = round(0.05 * Nvert);0286 scalp_intensity = zeros(N,1);0287 0288 fprintf(<span class="string">'...estimating scalp intensity from %d vertices\n'</span>, N);0289 fprintf(<span class="string">'...starting scalp intensity = %d\n'</span>, fit.val);0290 0291 indices = round(rand(1,N) * Nvert);0292 0293 i = 0;0294 <span class="keyword">for</span> v = indices,0295 0296 x = FV.vertices(v,1);0297 y = FV.vertices(v,2);0298 z = FV.vertices(v,3);0299 r = sqrt( (x-xo).^2 + (y-yo).^2 + (z-zo).^2 );0300 l = (x-xo)/r; <span class="comment">% cos alpha</span>0301 m = (y-yo)/r; <span class="comment">% cos beta</span>0302 n = (z-zo)/r; <span class="comment">% cos gamma</span>0303 0304 <span class="comment">% find discrete points from zero to the vertex</span>0305 POINTS = 250;0306 radial_distances = linspace(0,r,POINTS);0307 0308 L = repmat(l,1,POINTS);0309 M = repmat(m,1,POINTS);0310 N = repmat(n,1,POINTS);0311 0312 XI = (L .* radial_distances) + xo;0313 YI = (M .* radial_distances) + yo;0314 ZI = (N .* radial_distances) + zo;0315 0316 <span class="comment">% interpolate volume values at these points</span>0317 <span class="comment">% ( not sure why have to swap XI,YI here )</span>0318 VI = interp3(vol,YI,XI,ZI,<span class="string">'*nearest'</span>);0319 0320 <span class="comment">% find location in VI where intensity gradient is steep</span>0321 half_max = 0.5 * max(VI);0322 index_maxima = find(VI > half_max);0323 <span class="comment">% use the largest index value to locate the maxima intensity that lie</span>0324 <span class="comment">% furthest toward the outside of the head</span>0325 index_max = index_maxima(end);0326 0327 i = i + 1;0328 scalp_intensity(i,1) = VI(index_max);0329 0330 plot(radial_distances,VI)0331 0332 <span class="keyword">end</span>0333 0334 fit.val = mean(scalp_intensity);0335 fit.tol = std(scalp_intensity);0336 0337 fprintf(<span class="string">'...estimated scalp intensity = %g\n'</span>, fit.val);0338 fprintf(<span class="string">'...estimated tolerance intensity = %g\n'</span>, fit.tol);0339 0340 <span class="keyword">return</span>0341 0342 0343 0344 0345 0346 0347 0348 0349 0350 0351 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>0352 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>0353 <a name="_sub2" href="#_subfunctions" class="code">function [FV, intensityAtR, R] = locate_val(FV,vol,origin,fit),</a>0354 0355 xo = origin(1); yo = origin(2); zo = origin(3);0356 0357 Nvert = size(FV.vertices,1);0358 progress = round(.1 * Nvert);0359 0360 <span class="comment">% Initialise difference intensity & distance arrays</span>0361 intensityAtR = zeros(Nvert,1);0362 R = intensityAtR;0363 0364 <span class="comment">% Find distance and direction cosines for line from</span>0365 <span class="comment">% origin to all vertices</span>0366 XV = FV.vertices(:,1);0367 YV = FV.vertices(:,2);0368 ZV = FV.vertices(:,3);0369 RV = sqrt( (XV-xo).^2 + (YV-yo).^2 + (ZV-zo).^2 );0370 LV = (XV-xo)./RV; <span class="comment">% cos alpha</span>0371 MV = (YV-yo)./RV; <span class="comment">% cos beta</span>0372 NV = (ZV-zo)./RV; <span class="comment">% cos gamma</span>0373 0374 <span class="comment">% Check for binary volume data, if empty, binary</span>0375 binvol = find(vol > 1);0376 0377 <span class="comment">% Locate each vertex at a given fit value</span>0378 tic0379 <span class="keyword">for</span> v = 1:Nvert,0380 0381 <span class="keyword">if</span> v > progress,0382 fprintf(<span class="string">'...interp3 processed %4d of %4d vertices'</span>,progress,Nvert);0383 t = toc; fprintf(<span class="string">' (%5.2f sec)\n'</span>,t);0384 progress = progress + progress;0385 <span class="keyword">end</span>0386 0387 <span class="comment">% Find direction cosines for line from origin to vertex</span>0388 x = XV(v);0389 y = YV(v);0390 z = ZV(v);0391 d = RV(v);0392 l = LV(v); <span class="comment">% cos alpha</span>0393 m = MV(v); <span class="comment">% cos beta</span>0394 n = NV(v); <span class="comment">% cos gamma</span>0395 0396 <span class="comment">% find discrete points between origin</span>0397 <span class="comment">% and vertex + 20% of vertex distance</span>0398 POINTS = 250;0399 0400 Rarray = linspace(0,(d + .2 * d),POINTS);0401 0402 L = repmat(l,1,POINTS);0403 M = repmat(m,1,POINTS);0404 N = repmat(n,1,POINTS);0405 0406 XI = (L .* Rarray) + xo;0407 YI = (M .* Rarray) + yo;0408 ZI = (N .* Rarray) + zo;0409 0410 <span class="comment">% interpolate volume values at these points</span>0411 <span class="comment">% ( not sure why have to swap XI,YI here )</span>0412 VI = interp3(vol,YI,XI,ZI,<span class="string">'*linear'</span>);0413 0414 <span class="comment">% do we have a binary volume (no values > 1)</span>0415 <span class="keyword">if</span> isempty(binvol),0416 maxindex = max(find(VI>0));0417 <span class="keyword">if</span> maxindex,0418 R(v) = Rarray(maxindex);0419 <span class="keyword">end</span>0420 <span class="keyword">else</span>0421 <span class="comment">% find the finite values of VI</span>0422 index = max(find(VI(isfinite(VI))));0423 <span class="keyword">if</span> index,0424 0425 <span class="comment">% Find nearest volume value to the required fit value</span>0426 nearest = max(find(VI >= fit.val));0427 0428 <span class="comment">%[ nearest, value ] = NearestArrayPoint( VI, fit.val );</span>0429 0430 <span class="comment">% Check this nearest index against a differential</span>0431 <span class="comment">% negative peak value</span>0432 <span class="comment">%diffVI = diff(VI);</span>0433 <span class="comment">%if max(VI) > 1,</span>0434 <span class="comment">% diffindex = find(diffVI < -20);</span>0435 <span class="comment">%else</span>0436 <span class="comment">% probably a binary volume</span>0437 <span class="comment">% diffindex = find(diffVI < 0);</span>0438 <span class="comment">%end</span>0439 0440 <span class="comment">% now set d</span>0441 <span class="keyword">if</span> nearest,0442 R(v) = Rarray(nearest);0443 <span class="keyword">end</span>0444 <span class="keyword">end</span>0445 <span class="keyword">end</span>0446 0447 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>0448 <span class="comment">% Constrain relocation by fit.vattr,</span>0449 <span class="comment">% some % of distance from neighbours</span>0450 0451 vi = find(FV.edge(v,:)); <span class="comment">% the neighbours' indices</span>0452 X = FV.vertices(vi,1); <span class="comment">% the neighbours' vertices</span>0453 Y = FV.vertices(vi,2);0454 Z = FV.vertices(vi,3);0455 0456 <span class="comment">% Find neighbour distances</span>0457 RN = sqrt( (X-xo).^2 + (Y-yo).^2 + (Z-zo).^2 );0458 <span class="comment">% Find mean distance of neighbours</span>0459 neighbour_distances_mean = mean(RN);0460 0461 minattr = fit.vattr;0462 maxattr = 1 + fit.vattr;0463 0464 <span class="keyword">if</span> R(v) < (minattr * neighbour_distances_mean),0465 R(v) = minattr * neighbour_distances_mean;0466 <span class="keyword">end</span>0467 <span class="keyword">if</span> R(v) > (maxattr * neighbour_distances_mean),0468 R(v) = maxattr * neighbour_distances_mean;⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -