Skip to content

Commit

Permalink
Repair AxonSeg
Browse files Browse the repository at this point in the history
  • Loading branch information
tanguyduval committed Feb 15, 2018
1 parent d1aa6ad commit f204b56
Show file tree
Hide file tree
Showing 10 changed files with 46 additions and 52 deletions.
4 changes: 2 additions & 2 deletions AxonSeg.m
Expand Up @@ -15,10 +15,10 @@
% ----------------------------------------------------------------------------------------------------
% Example:
% Segment and generate a SegParameters.mat file:
% AxonSeg test_image_OM_crop.tif
% AxonSeg test_image_OM_2.tif
%
% Segment a new image using the same SegParameters:
% AxonSeg test_image_2.tif SegParameters.mat -nogui
% AxonSeg test_image_OM_2.tif SegParameters.mat -nogui
%
% Segment myelin using a corrected axon mask:
% AxonSeg({'RawImage.tif','Seg_mask_axon_corr.tif'}, 'SegParameters.mat', '-nogui')
Expand Down
2 changes: 1 addition & 1 deletion code/Blockprocessing/as_listcell2axonlist.m
Expand Up @@ -27,7 +27,7 @@
% list of cell 2 one list only
axonlist(length(listcell_seg))=listcell_seg{end};
[axonlist.seg]=deal(listcell_seg{:});
axonlist=cat(2,axonlist.seg); axonlist=cat(2,axonlist.seg); axonlist=axonlist(logical([axonlist.myelinAera])); % weird but works (fast).. keep everything
axonlist=cat(2,axonlist.seg); axonlist=cat(2,axonlist.seg); axonlist=axonlist(logical([axonlist.gRatio])); % weird but works (fast).. keep everything

% remove axons that have been segmented twice
centroids=cat(1,axonlist.Centroid);
Expand Down
5 changes: 1 addition & 4 deletions code/Seg_myelin/as_myelinseg2axonlist.m
Expand Up @@ -10,16 +10,13 @@
for iaxon=size(seg,3):-1:1
[axonlist(iaxon).data(:,1), axonlist(iaxon).data(:,2)]=find(seg(:,:,iaxon));
axonlist(iaxon).data = single(axonlist(iaxon).data);
axonlist(iaxon).myelinAera=size(axonlist(iaxon).data,1);
axonlist(iaxon).axonID=repmat(iaxon,[axonlist(iaxon).myelinAera,1]);
tmp=regionprops(seg(:,:,iaxon),'Centroid');
if isempty(tmp), axonlist(iaxon).Centroid=[];
else
axonlist(iaxon).Centroid=tmp(1).Centroid([2 1]);
axonlist(iaxon).Centroid=tmp(1).Centroid([2 1]);
end
for istat=1:length(sfields)
tmp=stats.(sfields{istat});
axonlist(iaxon).(sfields{istat})=tmp(iaxon);
end
end

72 changes: 36 additions & 36 deletions code/Seg_myelin/as_stats.m
@@ -1,45 +1,45 @@
function AXstat =as_stats(bwseg,PixelSize,myelin,metric)
% statis=as_stats(bwseg,pixelSize,myelin)
function AXstat =as_stats(myelinseg,pixelSize,myelin)
% AXstat=as_stats(myelinseg,pixelSize)
if ~exist('myelin','var'), myelin=1; end


bwseg=reshape2D(~~bwseg,1);
Npix = length(bwseg(:));
cc = bwconncomp(bwseg, 8);
if myelin
myelinseg = bwseg;
bwseg = logical(imfill(myelinseg,'holes') - myelinseg);
end

stats = regionprops(bwseg,{'Eccentricity', 'EquivDiameter','Area'});
myelinseg=reshape2D(myelinseg,1);
cc = bwconncomp(myelinseg, 8);
prop =regionprops(cc, {'Area', 'FilledArea'});

% convert units pixel to meters
EquivDiameter = [stats.EquivDiameter]*PixelSize;
% Creation of a struct var for the stats
AXstat=struct('axonEquivDiameter',0,'gRatio',0,'myelinThickness',0);

AXstat.EquivDiameter01 = sum(EquivDiameter>0 & EquivDiameter<1);
AXstat.EquivDiameter14 = sum(EquivDiameter>1 & EquivDiameter<4);
AXstat.EquivDiameter48 = sum(EquivDiameter>4 & EquivDiameter<8);
AXstat.EquivDiameter812 = sum(EquivDiameter>8 & EquivDiameter<12);
AXstat.EquivDiameter1217 = sum(EquivDiameter>12 & EquivDiameter<17);
AXstat.EquivDiameter = mean(EquivDiameter);
AXstat.EquivDiameter_W = sum(EquivDiameter.^3)/sum(EquivDiameter.^2);
AXstat.EquivDiameter_std = std(EquivDiameter);
AXstat.Eccentricity = max([stats.Eccentricity]);
AXstat.AVF = sum([stats.Area])/Npix;
AXstat.Naxons = length(stats);

if myelin
stats = regionprops(myelinseg,{'EquivDiameter'});
myelinEquivDiameter = [stats.EquivDiameter]*PixelSize;
AXstat.MVF = sum(myelinseg)/Npix;
AXstat.myelinEquivDiameter = median(myelinEquivDiameter);
try
AXstat.gRatio = single(EquivDiameter ./ myelinEquivDiameter);
catch
warning('Some myelin have no holes')
if ~isempty(prop)
% Area gives myelin area & FilledArea gives (myelin+axon) area
area = [[prop.Area]' [prop.FilledArea]'-[prop.Area]'];
% Calculate myelin area by using pixel size
if myelin
myelinArea = single(area(:, 1).*pixelSize^2);
end
% Calculate axon area by using pixel size
if myelin
axonArea = single(area(:, 2).*pixelSize^2);
else
axonArea = single(area(:, 1).*pixelSize^2);
end
% Myelin equiv. diameter = sqrt(4*TotalArea/pi)
if myelin
myelinEquivDiameter = single(sqrt(4*(myelinArea + axonArea)/pi));
end
% Axon equiv. diameter = sqrt(4*AxonArea/pi)
AXstat.axonEquivDiameter = single(sqrt(4*axonArea/pi));
% Calculate gRatio
if myelin
AXstat.gRatio = single(AXstat.axonEquivDiameter ./ myelinEquivDiameter);
% Myelin thickness = (MyelinDiam - AxonDiam)/2
AXstat.myelinThickness = single((myelinEquivDiameter - AXstat.axonEquivDiameter) / 2);
end

% if nargout>0
% MorphoStat.circularity=;
% end
end

if exist('metric','var')
AXstat = AXstat.(metric);
end

7 changes: 2 additions & 5 deletions code/Stats/as_stats_fields.m
@@ -1,6 +1,3 @@
function fields=as_stats_fields
function fieldsval=as_stats_fields

fields={...
'gRatio'...
'axonEquivDiameter'...
'myelinThickness'};
fieldsval=fields(as_stats([1 1; 1 1],1));
2 changes: 1 addition & 1 deletion code/axonlist/as_display_label.m
Expand Up @@ -35,7 +35,7 @@
for i=Naxon:-1:1
if ~mod(i,1000), disp(i); end
if size(axonlist(i).data,1)>5
index=round(double(axonlist(i).data)+repmat(axonlist(i).Centroid,[size(axonlist(i).data,1),1]));
index=round(double(axonlist(i).data));

% If 'axon' display type is specified, find axon index instead of
% myelin index
Expand Down
4 changes: 2 additions & 2 deletions code/batch_multiple_segmentations.m
Expand Up @@ -44,14 +44,14 @@
% imshow(display_1);

maxdiam=ceil(prctile(cat(1,axonlist.axonEquivDiameter),99));
imwrite(display_1,[output 'axonEquivDiameter_(axons)_0µm_' num2str(maxdiam) 'µm.jpg'])
imwrite(display_1,[output 'axonEquivDiameter_(axons)_0um_' num2str(maxdiam) 'um.jpg'])


%%

% reduce size of axonlist - eliminate all unecessary fields

fields_to_remove = {'data','axonID'}; % erasing these 2 fields reduces axonlist size by 90% 1 MB -> 113 KB for example
fields_to_remove = {'data'}; % erasing these 2 fields reduces axonlist size by 90% 1 MB -> 113 KB for example
axonlist_reduced=rmfield(axonlist,fields_to_remove);
% save('axonlist_reduced.mat','axonlist_reduced','-v7.3');

Expand Down
Binary file removed code/data/test_image_OM.tif
Binary file not shown.
Binary file added code/data/test_image_OM_2.tif
Binary file not shown.
2 changes: 1 addition & 1 deletion code/imageprocessing/as_improc_blockwising.m
Expand Up @@ -24,7 +24,7 @@
end
end

parfor(block=1:n_blocki*n_blockj,parforArg)
for (block=1:n_blocki*n_blockj)
j=1+(blocksize-overlap)*(mod(block-1,n_blockj));
i=1+(blocksize-overlap)*(block-mod(block-1,n_blockj)-1)/n_blockj;

Expand Down

0 comments on commit f204b56

Please sign in to comment.