vb_test_morphology

 Copyright (C) 2011, ATR All Rights Reserved.
 License : New BSD License(see VBMEG_LICENSE.txt)

0001 %
0002 % Copyright (C) 2011, ATR All Rights Reserved.
0003 % License : New BSD License(see VBMEG_LICENSE.txt)
0004 clear all
0005 pack
0006 
0007 udir= [getenv('MATHOME') '/SBIdata-new/BEM/head/'];
0008 
0009 fname='100008d-bem.voxel.mat';
0010 %fname='100008d-bem.mask.mat';
0011 
0012 load([udir fname]);
0013 
0014 % データの間引き
0015 step   = 2;
0016 % マスク閾値
0017 vlevel = 0.5;
0018 % Morfology operation
0019 R = 5; 
0020 
0021 profile_on = 0;
0022 
0023 % Make mask image
0024 [N1,N2,N3] = size(voxel);
0025 voxel       = voxel( 1:step:N1, 1:step:N2, 1:step:N3 );
0026 [N1,N2,N3] = size(voxel);
0027 
0028 ix = find( voxel > vlevel);
0029 voxel = zeros(N1,N2,N3);
0030 voxel(ix) = 1;
0031 
0032 %z=30;
0033 %h=image(voxel(:,:,z),'CDataMapping','scaled');
0034 
0035 %return
0036 
0037 
0038 if profile_on == 0,
0039     tic
0040 elseif profile_on > 0,
0041     profile off
0042     profile on -detail builtin
0043     tic
0044 end;
0045 
0046 % test morphology
0047 
0048 %B   = vb_dilation_3d(voxel,R);
0049 %
0050 %vb_ptime(toc);tic
0051 %
0052 %B1 = vb_dilation_3d_int(voxel,R);
0053 %
0054 %vb_ptime(toc);tic
0055 %
0056 %B2 = dilation_3d_new(voxel,R);
0057 %
0058 %vb_ptime(toc);tic
0059 
0060 %B3 = dilation_3d_bak(voxel,R);
0061 
0062 B = vb_erosion_3d(voxel,R);
0063 vb_ptime(toc);tic
0064 B1 = erosion_3d_bak(voxel,R);
0065 vb_ptime(toc);tic
0066 %B2 = vb_erosion_3d_int(voxel,R);
0067 %vb_ptime(toc);tic
0068 
0069 B1   = double(B1);
0070 B2   = B1 - B;
0071 err1 = sum(abs(B(:) - B1(:)))
0072 %err2 = sum(abs(B(:) - B2(:)))
0073 %err3 = sum(abs(B(:) - B2(:)))
0074 
0075 % 3次元データの平滑化
0076 %B = smooth3(B,'gaussian',3);
0077 % マスク境界面の抽出
0078 %[V,F]  = vb_get_boundary_mesh(B, vlevel);
0079 % 等特性サーフェスの抽出
0080 %[F,V] = isosurface(B,Val);
0081 
0082 if profile_on==0,
0083     vb_ptime(toc)
0084 elseif profile_on==1,
0085     vb_ptime(toc)
0086     if vb_matlab_version >= 7
0087         profile off
0088         profsave(profile('info'),'profile_job')
0089     else
0090         profile report profile_job
0091     end
0092 elseif profile_on==2,
0093     vb_ptime(toc)
0094     profile viewer
0095 end
0096 
0097 return
0098 
0099 %%%%%%%%%% Plot surf %%%%%%%%%%
0100 mode  = 4;        % 1: 5-class, 2: 6-class,  3:'gray', 4:'jet'
0101 zindx = [100 , 100, 180];    % slice coordinate
0102 slice = {'x'; 'y'; 'z'};    % slice cut direction
0103 
0104 zindx = [130:10:160];    % slice coordinate
0105 slice = 'z';    % slice cut direction
0106 NX=2;
0107 NY=2;
0108 
0109 figure;
0110 
0111 nfig = length(zindx);
0112 
0113 for n=1:nfig
0114     subplot(NX,NY,n);
0115     vb_plot_3d_image(B, zindx(n), mode, slice);
0116     title([slice '-slice (' num2str(zindx(n),3) ')'])
0117 %    vb_plot_3d_image(B, zindx(n), mode, slice{n});
0118 %    title([slice{n} '-slice (' num2str(zindx(n),3) ')'])
0119 end;
0120 
0121 return
0122 %%%%%%%%%% Plot surf %%%%%%%%%%
0123 vangle   = [80 30];
0124 fcolor     = [0.9, 0.8, 0.8];
0125 egcolor  = 'k';
0126 %egcolor  = 'none';
0127 
0128 patch('Faces',F,'Vertices',V,'FaceColor',fcolor,'EdgeColor',egcolor);
0129 
0130 %camlight(30,80);
0131 %camlight(-60,0);
0132 %camlight(60,0);
0133 %camlight(180,0);
0134 lighting phong; 
0135 material dull;
0136 
0137 view(vangle);
0138 axis equal
0139 
0140 zoomrate=1.5;
0141 
0142 zoom(zoomrate);