Skip to content
Snippets Groups Projects
Commit e2a2d609 authored by Laura Nervo's avatar Laura Nervo Committed by Nicola Vigano
Browse files

Boundaries : commenting and formatting


Signed-off-by: default avatarLaura Nervo <laura.nervo@esrf.fr>

git-svn-id: https://svn.code.sf.net/p/dct/code/trunk@581 4c865b51-4357-4376-afb4-474e03ccb993
parent bf1509e8
No related branches found
No related tags found
No related merge requests found
...@@ -2,12 +2,12 @@ ...@@ -2,12 +2,12 @@
function boundaries_structure = gtMakeBoundariesStructure(boundaries, grains) function boundaries_structure = gtMakeBoundariesStructure(boundaries, grains)
% from the boundaries and grains volumes, start to build the % from the boundaries and grains volume, start to build the
% boundaries_structure. % boundaries_structure.
nbounds=max(boundaries(:)); nbounds=max(boundaries(:));
for i=1:max(boundaries(:)) for i=1:nbounds
list=find(boundaries(:,1)==i); list=find(boundaries(:,1)==i);
list=[list;find(boundaries(:,2)==i)]; list=[list;find(boundaries(:,2)==i)];
boundaries_structure(i).count=length(list); boundaries_structure(i).count=length(list);
...@@ -19,12 +19,10 @@ for i=1:max(boundaries(:)) ...@@ -19,12 +19,10 @@ for i=1:max(boundaries(:))
disp('one voxel only') disp('one voxel only')
boundaries_structure(i).grain1=list(1); boundaries_structure(i).grain1=list(1);
boundaries_structure(i).grain2=[]; boundaries_structure(i).grain2=[];
else
keyboard
end end
if mod(i, 100)==0 if mod(i, 100)==0
disp(sprintf('done %d percent', round(100*i/nbounds))) fprintf('done %d percent\n', round(100*i/nbounds))
end end
end end
......
...@@ -41,12 +41,11 @@ load('parameters.mat'); ...@@ -41,12 +41,11 @@ load('parameters.mat');
spacegroup=parameters.cryst.spacegroup; spacegroup=parameters.cryst.spacegroup;
if spacegroup==225 || spacegroup==229 if spacegroup==225 || spacegroup==229
sym = gtGetCubicSymOp; sym = gtGetCubicSymOp();
elseif spacegroup==663 || spacegroup==194 elseif spacegroup==663 || spacegroup==194
sym= gtGetHexagonalSymOp_sab; sym= gtGetHexagonalSymOp_sab();
elseif spacegroup==167 elseif spacegroup==167
% sym= gtGetHexagonalSymOp_sab; sym=gtGetAluminaSymOp_sab();
sym=gtGetAluminaSymOp_sab;
disp('alumina sym ops testing...') disp('alumina sym ops testing...')
else else
disp('spacegroup not supported! cant do symmetry operations!') disp('spacegroup not supported! cant do symmetry operations!')
...@@ -117,6 +116,9 @@ if 1 % if using boundaries_structure ...@@ -117,6 +116,9 @@ if 1 % if using boundaries_structure
weights=[weights; weights]; weights=[weights; weights];
output=poles; output=poles;
output2=weights; output2=weights;
save poles poles
save weights weights
end end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
...@@ -162,3 +164,4 @@ if 0 ...@@ -162,3 +164,4 @@ if 0
count count
end end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end % end function
\ No newline at end of file
function [vol_out, boundaries, boundaries_count]=gt_find_grain_boundaries(vol) function [vol_out, boundaries, boundaries_count] = gt_find_grain_boundaries(vol)
% [vol_out, boundaries, boundaries_count]=gt_find_grain_boundaries(vol) % GT_FIND_GRAIN_BOUNDARIES Starts writing the boundaries_structure.mat file
% --------------------------------------------------------------------- % [vol_out, boundaries, boundaries_count] = gt_find_grain_boundaries(vol)
% Starts writing the boundaries_structure.mat file % -----------------------------------------------------------------------
% now need to read the grain boundaries from the stack with dilated grains % now need to read the grain boundaries from the stack with dilated grains
% Saves boundaries_structure as tmp.mat file % Saves boundaries_structure as tmp.mat file
% %
...@@ -15,7 +15,8 @@ function [vol_out, boundaries, boundaries_count]=gt_find_grain_boundaries(vol) ...@@ -15,7 +15,8 @@ function [vol_out, boundaries, boundaries_count]=gt_find_grain_boundaries(vol)
% %
% this actually shares functionality with gtMakeBoundariesStructure, which % this actually shares functionality with gtMakeBoundariesStructure, which
% also suppresses single voxel "boundaries". Should merge... % also suppresses single voxel "boundaries". Should merge...
%
% Version 002 by LNervo
boundaries=[]; boundaries=[];
boundaries_count=0; boundaries_count=0;
...@@ -24,8 +25,8 @@ boundary_number=1; ...@@ -24,8 +25,8 @@ boundary_number=1;
[size_x, size_y, size_z]=size(vol); [size_x, size_y, size_z]=size(vol);
vol_out=zeros(size(vol)); vol_out=zeros(size(vol));
%find the voxels that need to be studied % find the voxels that need to be studied
%dilates colours % dilates colours
vol_dilated = imdilate(vol, ones(3,3,3)); vol_dilated = imdilate(vol, ones(3,3,3));
vol_eroded = imerode(vol,ones(3,3,3)); vol_eroded = imerode(vol,ones(3,3,3));
dif = vol_dilated-vol_eroded; dif = vol_dilated-vol_eroded;
...@@ -34,10 +35,9 @@ todo_voxels = find(dif~=0); ...@@ -34,10 +35,9 @@ todo_voxels = find(dif~=0);
for i=1:length(todo_voxels) for i=1:length(todo_voxels)
[x,y,z]=ind2sub(size(vol), todo_voxels(i)); [x,y,z]=ind2sub(size(vol), todo_voxels(i));
%look at voxel neighbourhood % look at voxel neighbourhood
xstart=x-1; xstart=x-1;
xend=x+1; xend=x+1;
ystart=y-1; ystart=y-1;
...@@ -65,50 +65,52 @@ for i=1:length(todo_voxels) ...@@ -65,50 +65,52 @@ for i=1:length(todo_voxels)
neighbours=vol(xstart:xend, ystart:yend, zstart:zend); neighbours=vol(xstart:xend, ystart:yend, zstart:zend);
neighbours=neighbours(:); neighbours=neighbours(:);
%neighbours(find(neighbours==0))=[];%don't consider boundaries to zeros % neighbours(find(neighbours==0))=[]; % don't consider boundaries to zeros
voxel_colour=vol(x,y,z); voxel_colour=vol(x,y,z);
%find any different coloured neighbours % find any different coloured neighbours
neighbours=neighbours(find(neighbours~=voxel_colour)); neighbours=neighbours(find(neighbours~=voxel_colour));
if ~isempty(neighbours) if ~isempty(neighbours)
neighbour_colour=min(neighbours);%lowest non-zero neighbour colour neighbour_colour=min(neighbours); % lowest non-zero neighbour colour
%has this combination of voxel_colour and neighbour_colour already been % has this combination of voxel_colour and neighbour_colour already been
%treated? look for rows in boundaries % treated? look for rows in boundaries
dummy1=find(boundaries(:)==voxel_colour); dummy1=find(boundaries(:)==voxel_colour);
[dummy1,a]=ind2sub(size(boundaries),dummy1); [dummy1,a]=ind2sub(size(boundaries),dummy1);
dummy2=find(boundaries(:)==neighbour_colour); dummy2=find(boundaries(:)==neighbour_colour);
[dummy2,a]=ind2sub(size(boundaries),dummy2); [dummy2,a]=ind2sub(size(boundaries),dummy2);
%any common rows? % any common rows?
boundary_number=dummy1(find(ismember(dummy1,dummy2))); boundary_number=dummy1(find(ismember(dummy1,dummy2)));
if ~isempty(boundary_number) %a previously encountered boundary if ~isempty(boundary_number) % a previously encountered boundary
vol_out(x,y,z)=boundary_number; vol_out(x,y,z)=boundary_number;
boundaries_count(boundary_number)=boundaries_count(boundary_number)+1; boundaries_count(boundary_number)=boundaries_count(boundary_number)+1;
%add to boundaries_structure % add to boundaries_structure
boundaries_structure(boundary_number).count=boundaries_structure(boundary_number).count+1; boundaries_structure(boundary_number).count=boundaries_structure(boundary_number).count+1;
boundary_number=[]; boundary_number=[];
else %a new boundary else % a new boundary
boundary_number=size(boundaries,1)+1; boundary_number=size(boundaries,1)+1;
disp(['now on boundary' num2str(boundary_number)]) disp(['now on boundary: ' num2str(boundary_number)])
vol_out(x,y,z)=boundary_number; vol_out(x,y,z)=boundary_number;
boundaries(boundary_number,:)=[voxel_colour, neighbour_colour]; boundaries(boundary_number,:)=[voxel_colour, neighbour_colour];
boundaries_count(boundary_number)=1; boundaries_count(boundary_number)=1;
%add to boundaries_structure % add to boundaries_structure
boundaries_structure(boundary_number).grain1=voxel_colour; boundaries_structure(boundary_number).grain1=voxel_colour;
boundaries_structure(boundary_number).grain2=neighbour_colour; boundaries_structure(boundary_number).grain2=neighbour_colour;
boundaries_structure(boundary_number).count=1; boundaries_structure(boundary_number).count=1;
boundary_number=[]; boundary_number=[];
end end
end%if a boundary voxel end % if a boundary voxel
end %loop over todo voxels
save tmp boundaries_structure end % loop over todo voxels
save('tmp.mat','boundaries_structure');
end % end of function
% try % try
% if max(vol_out(:))<2^16 % if max(vol_out(:))<2^16
% edf_write(vol_out, '5_reconstruction/boundaries_volume.edf', 'uint16') % edf_write(vol_out, '5_reconstruction/boundaries_volume.edf', 'uint16')
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment