scan_2_noise

Outline
Recenter isocenter
Refine isocenter
Create structuring elements
Noise measurements
External dependencies

Outline

This is a multistage process. First the scan is cropped to the area around the carina. The second step is to recenter the isocenter on that crop, to ensure a specific height above the carinal ridge, which is then refined in x and y to center it in the trachea. This is used to crop the image down to the final subsegment.

The next step is to create the structuring elements of the correct size. This will create a circular SE with a 1 cm² area, a spherical SE with the same radius (≈0.75 mm³), and a spherical SE that can be used to erode the outer part of the segmentation to remove edge artifacts.

An example CT file can be downloaded here (permalink). This zip file contains a text file with a citation for the data source, as well as a .mat file with the image, the relevant dicominfo fields, and a suggested isocenter.

function result=scan_2_noise(IM,info,isoc,HeightFromCarina,window)
%input:
% IM: CT image in CT numbers
% info: dicominfo with these fields:
%            - PixelSpacing (0028,0030)
%            - SpacingBetweenSlices (0018,0088)
%            - RescaleSlope (0028,1053)
%            - RescaleIntercept (0028,1052)
% isoc: initial seed, i.e. a coordinate in the trachea, not too far from carina
% height_from_ridge: the height in mm from the carinal ridge, rounded to the nearest slice
% window: [optional] a vector with the window of segIM (by default [-1500 100] i.e. W1600L-700)
%
%ouput:
% result: struct with these fields:
%            - ROI: 27x1 vector with noise measurements (negative if invalid)
%            - VOI: 27x1 vector with noise measurements (negative if invalid)
%            - seg: scalar with noise measurement
%            - segIM: RGB image with the scan and segmentations (window=[-1500 100] ie W1600L-700)

if nargin<=4,window=[-1500 100];end

PixelSize=[info.PixelSpacing info.SpacingBetweenSlices];

%use larger subsegment to recenter the isocenter
x=75;
dim1=(-x:x)+isoc(1);
dim2=(-x:x)+isoc(2);
dim3=0;
while any(dim3<=0) || any(dim3>size(IM,3))%account for scans with only few (thick) slices
    dim3=(-x:x)+isoc(3);x=x-1;
end
IM=IM(dim1,dim2,dim3);

%find the recentered isoc
if nargin<4
    HeightFromCarina=15;%in mm
end
isoc=recenter_isoc(IM,ceil(size(IM)/2),PixelSize,HeightFromCarina);
isoc=refine_isoc(IM,isoc,PixelSize);

%use updated isocenter for final subsegment
dim1=(-30:30)+isoc(1);
dim2=(-30:30)+isoc(2);
dim3=(-30:30)+isoc(3);
IM=IM(dim1,dim2,dim3);


%get appropriate SE size
[ROI,VOI,SE_reg]=load_SE(PixelSize);
%get the rescale and intercept from the dicominfo
RI=[info.RescaleSlope info.RescaleIntercept];

%measure SD
result=struct;
[result.ROI,result.VOI,result.seg,result.segIM]=...
    getSD(IM,ROI,VOI,SE_reg,ceil(size(IM)/2),RI,window);
end

Recenter isocenter

This function updates the isocenter to be 1.5cm above the carinal ridge. The height of the carinal ridge is defined as the first slice where the ridge is substantial enough to split the the trachea segmentation in two.

To do this it uses region growing to determine a rough trachea mask. The difference threshold for the region growing is automatically determined (it is the lowest threshold starting from 15 where the top slice contains a mask voxel). This mask is then cleaned up with a morphological close operation with a 2 voxel ball.

The skeletonization of the trachea is not a true skeletonization, as that is fairly complex to implement efficiently. The result would probably also require a lot of manual pruning. So instead an object tracking strategy is used. A binary labeling is used on each slice. At some point there will be two areas in the mask that overlap with pixels marked in the slice above it. This must be at the location of the carina, except for double tube tracheas, which don't exist in the current dataset.

Note that the threshold determination takes the longest time of this process. Each incorrect guess for the threshold will require an additional 10 seconds (assuming the total time is about 20 seconds if the threshold does not have to be adjusted).

function [new_isoc,mask,trachea,maxDiff]=recenter_isoc(IM,isoc,px_sz,H2Carina)
%IM is still uint16 here, which doesn't really matter, because the below steps will only rely on
%absolute and relative differences. The third and fourth inputs are the pixel size and the required
%height to the carina (in the same units).

%Repeat the top slice a few times to get a more consistent result with the region growing, imclose
%and skeletonization.
extra_slices=20;
top=IM(:,:,1);top=repmat(top,1,1,extra_slices);
isoc(3)=isoc(3)+extra_slices;
IM=cat(3,top,IM);

%Do a first order approximate segmentation of the trachea. A leak is not likely for such a small
%scan section, and coverage should be good, especially after an imclose operation.
maxDiff=find_maxDiff(IM,isoc);%account for high noise and/or bad seed
maxDiff=ceil(maxDiff*1.25);%don't use absolute minimum value
mask=RegGrow(IM,maxDiff,isoc,'silent',true);
mask=noToolbox__IPT('imclose',mask,2);%2 pixel radius ball

%Remove the extra slice to return to the native size.
mask(:,:,1:extra_slices)=[];

%do a de facto skeletonization
prev=noToolbox__IPT('bwlabel',mask(:,:,1))==1;%initialize with first slice
trachea=false(size(mask));full_trachea=trachea;
conn_min=2*ndims(prev);%minimal connectivity for bwlabel
conn_max=(3^ndims(prev)-1);%maximal connectivity for bwlabel
for slice=2:size(mask,3)
    %store result to array, so it can be used to find the centroid
    trachea(:,:,slice-1)=prev;

    current_slice=noToolbox__IPT('bwlabel',mask(:,:,slice),conn_min);
    x=current_slice(prev);
    bins=0:(max(x)+1);
    count=histc(x,bins-0.5); %#ok<HISTC>
    count(1)=[];bins(1)=[];%remove the 0-label (i.e. non-trachea)

    %orignally: if nnz(count)>2
    [count,order]=sort(count);
    if count(end-1) >= count(end)/3
        split_found=slice;
        break
    else
        split_found=0;
        %find the ID of the trachea (probably 1)
        x=bins(order(end));
    end
    prev=current_slice==x;

    %store a more generous label for determination of centerline (using max
    %connectivity)
    current_slice=noToolbox__IPT('bwlabel',mask(:,:,slice),conn_max);
    IDs_of_trachea=unique(current_slice(prev));
    IDs_of_trachea(IDs_of_trachea==0)=[];
    full_trachea(:,:,slice)=current_slice==IDs_of_trachea(1);
    for n=2:numel(IDs_of_trachea)
        full_trachea(:,:,slice)=current_slice==IDs_of_trachea(n);
    end
end

if split_found~=0%not failed
    %Move up # cm and find the centroid of that slice's mask.
    z=round(split_found-H2Carina/px_sz(3));
    %fit a 3D line to estimate centroid
    %fit method from https://www.mathworks.com/matlabcentral/answers/78363#answer_88082
    x_fit_str='r(1)+(slice-r(3))/V(3,1)*V(1,1)';
    y_fit_str='r(2)+(slice-r(3))/V(3,1)*V(2,1)';
    try
        x_fit=eval(['@(slice,r,V) ' x_fit_str]);
        y_fit=eval(['@(slice,r,V) ' y_fit_str]);
    catch
        x_fit=inline(x_fit_str,'slice','r','V'); %#ok<DINLN>
        y_fit=inline(y_fit_str,'slice','r','V'); %#ok<DINLN>
    end
    upper_trachea=full_trachea;upper_trachea(:,:,z+1:end)=0;
    [xyz_x,xyz_y,xyz_z]=findND(upper_trachea);
    xyz=[xyz_x,xyz_y,xyz_z];
    r=mean(xyz);
    xyz=minus(xyz,repmat(r,size(xyz,1),1)); %equivalent to xyz=bsxfun(@minus,xyz,r);
    [ignore,ignore,V]=svd(xyz,0); %#ok<ASGLU>
    new_isoc=round([x_fit(z,r,V)   y_fit(z,r,V)   z]);

    % %This code works well in most cases, but causes a major shift if one bronchus is very
    % %horizontal. To avoid this, this code was removed. The refine_isoc function should be more
    % %succesful than this recentering anyway.
    % slice=full_trachea(:,:,z);
    % [x,y]=find(slice);
    % x=mean(x);y=mean(y);
    % %Treat this as an initial guess. Determine the middle along the specific centerlines.
    % isoc_guess=round([x y z]);
    %
    % new_isoc=isoc_guess;it=0;
    % while it==0 || (it<10 && sum(new_isoc-isoc_guess)>=2)
    %     it=it+1;isoc_guess=new_isoc;
    %     x=isoc_guess(1);y=isoc_guess(2);
    %     [x_,y_]=deal(x,y);
    %     x=mean(find(slice(:,round(y_))));
    %     y=mean(find(slice(round(x_),:)));
    %     new_isoc=round([x y z]);
    % end
else
    error('no carina found in image')
end
end
function maxDiff=find_maxDiff(IM,isoc)
%find lowest maxDiff value where the trachea mask reaches the top slice
maxDiff=15;%don't trust values lower than this
mask=[false false;false false];%initialize to start loop
while ~any(any(mask(:,:,1)))
    maxDiff=maxDiff+1;
    mask=RegGrow(IM,maxDiff,isoc,'silent',true);
end
end

Refine isocenter

This function refines the (x,y)-position of the isocenter. It will take into account three different isocenters:

The unrefined isocenter.
The weighted center of gravity.
The weighted center point (the point furthest from the edge)

The weights of the latter two are the number of voxels in the VOI for each slice. The latter two are only considered if they are sufficiently close to the unrefined isocenter, and are sufficiently close to each other.

The three isocenters are averaged, disregarding rejected isocenters. If only the original unrefined isocenter remains this will trigger a warning message. Such a situation may require manually confirming the isocenter.

function isoc=refine_isoc(IM,isoc,PixelSize)
%get the subimage
dim1=(-30:30)+isoc(1);
dim2=(-30:30)+isoc(2);
dim3=(-30:30)+isoc(3);
IM=IM(dim1,dim2,dim3);
centroid=ceil(size(IM)/2);

maxDiff=200;%overgrow to make sure all air is included
seg=RegGrow(IM,maxDiff,centroid,'nHood',26,'waitbar',false);
[ignore,VOI,SE_reg]=load_SE(PixelSize); %#ok<ASGLU>
seg=noToolbox__IPT('imerode',seg,SE_reg);%erode by 2mm to remove the edge

weights=zeros(size(IM,3),1);
d=struct;d.x=weights;d.y=weights;d.k=1;%short-hand for the deltas
d1=d;d2=d;
for z=centroid(3)+unique(VOI.z)'
    slice=seg(:,:,z);
    weights(z)=sum(VOI.z==(z-centroid(3)));

    %center of gravity
    [x,y,ignore_val]=find(slice); %#ok<ASGLU>
    d1.x(z)=mean(x)-centroid(1);
    d1.y(z)=mean(y)-centroid(2);

    %longest distance to edge
    dist_to_edge=noToolbox__IPT('bwdist',~slice);
    dist_to_edge=round(dist_to_edge*4)/4;%round to 0.25
    % %this code only returns the first pixel, and is therefore biased to the top left
    % [ign,ind]=max(dist_to_edge(:));[x,y]=ind2sub(size(IM(:,:,1)),ind); %#ok<ASGLU>
    [x,y]=find(dist_to_edge==max(dist_to_edge(:)));
    [ign,ind]=min(abs(x-mean(x)));x=x(ind); %#ok<ASGLU> Find the location closest to the middle,
    [ign,ind]=min(abs(y-mean(y)));y=y(ind); %#ok<ASGLU> which prevents 2 maxima causing issues.
    d2.x(z)=x-centroid(1);d2.y(z)=y-centroid(2);
end

d0.k=1;
d1.x=sum(d1.x.*weights)/sum(weights);
d1.y=sum(d1.y.*weights)/sum(weights);
d2.x=sum(d2.x.*weights)/sum(weights);
d2.y=sum(d2.y.*weights)/sum(weights);
d_1_2=sqrt((d1.x-d2.x).^2+(d1.y-d2.y).^2);
d_0_1=sqrt((0-d1.x).^2+(0-d1.y).^2);
d_0_2=sqrt((0-d2.x).^2+(0-d2.y).^2);
d1.k=determine_k(d_0_1,d_1_2);
d2.k=determine_k(d_0_2,d_1_2);

d.x=(d1.x+d2.x)/(d0.k+d1.k+d2.k);
d.y=(d1.y+d2.y)/(d0.k+d1.k+d2.k);
new_isoc_delta=[d.x d.y 0];

if sum(d1.k+d2.k)==0
    warning('Refinement failed, this script will continue with the unrefined isocenter.');
else
    isoc=isoc+round(new_isoc_delta);
end
end
function k=determine_k(d_0_1,d_1_2)
%Determine weight depending on relative distances. This way the isoc will only move long distances
%if the two methods agree.
limit=[5 10];
if d_0_1<limit(1) %close to old isoc, assume this is correct
    k=1;
elseif d_0_1<limit(2) %medium distance, only allow if the two new ones are close to eachother
    if d_1_2<limit(1)
        k=1;
    else
        k=0;
    end
else %long distance, assume this is an error
    k=0;
end
end

Create structuring elements

The ROI and VOI SEs created in this function are returned as coordinate offsets.

For improved performance on muliple scans it will store the results in a persistent variable, which will facilitate quick loading. Pixel sizes are rounded to 4 digits when comparing to a previous input.

function [ROI,VOI,SE_reg]=load_SE(PixelSize)
%get appropriate SE size
%reg is an SE that shaves off 2 mm

persistent SE
if ~isa(SE,'struct')
    SE=struct('PixelSize',{{}},'ROI',{{}},'VOI',{{}},'SE_reg',{{}});
end

hash_str=sprintf('%.4f||',PixelSize);
L=ismember(SE.PixelSize,{hash_str});
if any(L)
    ROI=SE.ROI(L);VOI=SE.VOI(L);SE_reg=SE.SE_reg(L);
    %remove cell layer
    ROI=ROI{1};VOI=VOI{1};SE_reg=SE_reg{1};
    return
end
%the SE doesn't exist yet for this size voxel

%make a 1cm2 ROI and a 0.75cm3 VOI
%r should be in the same units as PixelSize
%1 cm2 is 100 mm2, A=pi*r^2 => r=sqrt(A/pi)
r=sqrt(100/pi);

radius_in_pixels=ceil(r/min(PixelSize));
r2=radius_in_pixels;

[X,Y,Z]=ndgrid(...
    (-r2:r2)*PixelSize(1),...
    (-r2:r2)*PixelSize(2),...
    (-r2:r2)*PixelSize(3));
VOI=sqrt(X.^2+Y.^2+Z.^2)<=r;
[x,y,z,val]=findND(VOI); %#ok<ASGLU>
VOI=struct;
[VOI.x,VOI.y,VOI.z]=deal(x-r2-1,y-r2-1,z-r2-1);

SE_reg=sqrt(X.^2+Y.^2+Z.^2)<=2;

[X,Y,Z]=ndgrid(...
    (-r2:r2)*PixelSize(1),...
    (-r2:r2)*PixelSize(2),...
    0);
ROI=sqrt(X.^2+Y.^2+Z.^2)<=r;
[x,y,z,val]=findND(ROI); %#ok<ASGLU>
ROI=struct;
[ROI.x,ROI.y,ROI.z]=deal(x-r2-1,y-r2-1,zeros(size(x)));

%save for later scans
SE.PixelSize(numel(SE.PixelSize)+1)={hash_str};
SE.ROI(numel(SE.ROI)+1)={ROI};
SE.VOI(numel(SE.VOI)+1)={VOI};
SE.SE_reg(numel(SE.SE_reg)+1)={SE_reg};
end

Noise measurements

In this function the actual noise measurement is performed.

Note that the noise is returned as a negative number if the measurement is outside the larger trachea segmentation.

The image this returns is a montage for the position around the isocenter.

function [o1,o2,o3,o4]=getSD(IM,ROI,VOI,SE_reg,cur_isoc,RI,window)
%seg data is from the eroded region growing
%other data is negative if outside of the non-eroded segmentation
%
%syntax (output names shortened to fit on one line):
%   [result_ROI,result_VOI,result_seg,segIM]=...
%       getSD(IM,ROI,VOI,SE_reg,cur_isoc,RI,window)

persistent shifts n_shift last_msg
if isempty(shifts)
    dims=3;
    kernel=logical(ones(3*ones(1,dims))); %#ok<LOGL>
    kernel((numel(kernel)+1)/2)=false;%put centroid last (skip first)
    c=cell(dims,1);
    [c{:}]=ind2sub(size(kernel),find(kernel(:)));
    for n=1:numel(c)
        c{n}=c{n}-(size(kernel,n)+1)/2;
    end
    shifts=[c{:}];
    shifts(end+1,1)=0;%add centroid back
    n_shift=size(shifts,1);
end

%use region growing on the image to segment the trachea
maxDiff=200;%overgrow to make sure all air is included
seg=RegGrow(IM,maxDiff,cur_isoc,'nHood',26,'waitbar',false);
seg_small=noToolbox__IPT('imerode',seg,SE_reg);%erode by 2mm to remove edge effects
%apply rescale and intercept:
data1=RI(1)*double(IM(seg_small))+RI(2);
result_seg=std(data1);

result_VOI=zeros(n_shift,1);
for n=1:n_shift
    x=VOI.x+cur_isoc(1)+shifts(n,1);
    y=VOI.y+cur_isoc(2)+shifts(n,2);
    z=VOI.z+cur_isoc(3)+shifts(n,3);
    remove_ind=...
        (x<1 | y<1 |z<1) | ...
        (x>size(IM,1) | y>size(IM,2) |z>size(IM,3));
    x(remove_ind)=[];y(remove_ind)=[];z(remove_ind)=[];
    inds=sub2ind(size(IM),x,y,z);
    VOI_seg=zeros(size(IM));VOI_seg(inds)=1;VOI_seg=logical(VOI_seg);
    %apply rescale and intercept:
    data2=RI(1)*double(IM(inds))+RI(2);
    result_VOI(n)=std(data2);
    %set to negative if outside of segmentation
    if any(~seg(inds))
        result_VOI(n)=-std(data2);
    end
    %if any(~seg_small(inds))
    if any(~seg(inds))
        str=sprintf('warning: %d/%d of %d outside of seg',...
            sum(~seg_small(inds)),sum(~seg(inds)),numel(inds));
        if ~strcmp(last_msg,str)
            fprintf('%s\n',str)
            last_msg=str;
        end
    end
end

result_ROI=zeros(n_shift,1);
for n=1:n_shift
    x=ROI.x+cur_isoc(1)+shifts(n,1);
    y=ROI.y+cur_isoc(2)+shifts(n,2);
    z=ROI.z+cur_isoc(3)+shifts(n,3);
    remove_ind=...
        (x<1 | y<1 |z<1) | ...
        (x>size(IM,1) | y>size(IM,2) |z>size(IM,3));
    x(remove_ind)=[];y(remove_ind)=[];z(remove_ind)=[];
    inds=sub2ind(size(IM),x,y,z);
    ROI_seg=zeros(size(IM));ROI_seg(inds)=1;ROI_seg=logical(ROI_seg);
    %apply rescale and intercept:
    data3=RI(1)*double(IM(inds))+RI(2);
    result_ROI(n)=std(data3);
    %set to negative if outside of segmentation
    if any(~seg(inds))
        result_ROI(n)=-std(data3);
    end
    %if any(~seg_small(inds))
    if any(~seg(inds))
        str=sprintf('warning: %d/%d of %d outside of seg',...
            sum(~seg_small(inds)),sum(~seg(inds)),numel(inds));
        if ~strcmp(last_msg,str)
            fprintf('%s\n',str)
            last_msg=str;
        end
    end
end

%create montage to confirm correct segmentation

%Create the image with the scan as a base, marking big_seg with green, seg with yellow, VOI_seg
%with blue, and ROI_seg with red.
I=repmat(RI(1)*double(IM)+RI(2),1,1,1,3);
I(I<window(1))=window(1);
I(I>window(2))=window(2);
I=I-window(1);I=I/(window(2)-window(1));
[sz1,sz2,sz3]=size(IM);
green =repmat(cat(4,  0,1.0,  0),sz1,sz2,sz3,1);
yellow=repmat(cat(4,1.0,1.0,  0),sz1,sz2,sz3,1);
red   =repmat(cat(4,1.0,  0,  0),sz1,sz2,sz3,1);
blue  =repmat(cat(4,  0,0.4,1.0),sz1,sz2,sz3,1);%'Brandeis blue'

L=repmat(seg,1,1,1,3);
I(L)=green(L);
L=repmat(seg_small,1,1,1,3);
I(L)=yellow(L);
L=repmat(VOI_seg,1,1,1,3);
I(L)=blue(L);
L=repmat(ROI_seg,1,1,1,3);
I(L)=red(L);

%call custom implementation for silent operation
segIM=noToolbox__IPT('montage',permute(I,[1 2 4 3]),'NoFigure',true);

%keep syntax short
[o1,o2,o3,o4]=deal(result_ROI,result_VOI,result_seg,segIM);
end

External dependencies

The code above needs external depencies to work. Below are all the functions require to allow them to work.

% Used with permission.
% Original license:
% CC by-nc-sa 4.0 ( https://creativecommons.org/licenses/by-nc-sa/4.0 )
% H.J. Wisselink
%
% included functions (internal dependencies not listed):
%  For minify (version 1.0.1) see https://www.mathworks.com/matlabcentral/fileexchange/84717
%  noToolbox__IPT (version 0.2 [unpublished alpha])
%  findND (version 1.2.1)
%   For the normal version see https://www.mathworks.com/matlabcentral/fileexchange/64383
%  RegGrow (version 1.1.1)
%   For the normal version see https://www.mathworks.com/matlabcentral/fileexchange/72944


function varargout=findND(v000,varargin),...
if ~(isnumeric(v000) || islogical(v000)) || numel(v000)==0,error('HJW:findND:FirstInput',...
'Expected first input (X) to be a non-empty numeric or logical array.'),end,switch ...
nargin,case 1,v001='first';v002=inf;case 2,v001='first';v002=varargin{1};if ~(isnumeric(v002) ...
|| islogical(v002)) || numel(v002)~=1 || any(v002<0),error('HJW:findND:SecondInput',...
'Expected second input (K) to be a positive numeric or logical scalar.'),...
end,case 3,v002=varargin{1};if ~(isnumeric(v002) ...
|| islogical(v002)) || numel(v002)~=1 || any(v002<0),error('HJW:findND:SecondInput',...
'Expected second input (K) to be a positive numeric or logical scalar.'),...
end,v001=varargin{2};if ~isa(v001,...
'char') || ~( strcmpi(v001,'first') || strcmpi(v001,'last')),error('HJW:findND:ThirdInput',...
'Third input must be either ''first'' or ''last''.'),end,v001=lower(v001);otherwise,...
error('HJW:findND:InputNumber','Incorrect number of inputs.'),end,v003=length(size(v000));
if nargout>1 && nargout<v003,error('HJW:findND:Output','Incorrect number of output arguments.'),...
end,persistent v004,if isempty(v004),v004=f09('<',7,'Octave',...
'<',3);end,varargout=cell(nargout,1);if v004,if nargout>v003,[v005,v006,v007]=find(v000(:));
if length(v005)>v002,if strcmp(v001,'first'),v005=v005(1:v002);v007=v007(1:v002);
else,v005=v005((end-v002+1):end);v007=v007((end-v002+1):end);end,end,[varargout{1:(end-1)}] ...
= ind2sub(size(v000),v005);varargout{end}=v007;else,v005=find(v000);if length(v005)>v002,...
if strcmp(v001,'first'),v005=v005(1:v002);else,v005=v005((end-v002):end);end,end,[varargout{:}] ...
= ind2sub(size(v000),v005);end,else,if nargout>v003,[v005,v006,v007]=find(v000(:),...
v002,v001);[varargout{1:(end-1)}] = ind2sub(size(v000),v005);varargout{end}=v007;
else,v005=find(v000,v002,v001);[varargout{:}] = ind2sub(size(v000),v005);end,end,end
function varargout=noToolbox__IPT(v000,varargin),persistent ...
v001 v002 v003,if isempty(v001),v001={ 'imopen',@f05;'imclose',@f02;'imdilate',@f03;'imerode',@f04;
'bwlabel',@f00;'bwlabeln',@f00;'bwdist',@f10;'montage',@f06};for v004=1:size(v001,...
1),v001{v004,3}=v001{v004,2};v001{v004,2}=func2str(v001{v004,2});end,v002={'help','-help','h',...
'-h','doc'};v003=func2str(@noToolbox__IPT);end,if nargin==0,clc,help(mfilename),return,end,if ...
isa(v000,'string'),v000=char(v000);elseif ~isa(v000,'char'),error('HJW:noToolbox:InvalidInput',...
'The first parameter must be either a function name or ''help''.'),...
end,if ismember({lower(v000)},v002),if nargin==1,error('HJW:noToolbox:InvalidInput',...
['The functionname parameter must be the ','name of an implemented function.']),end,...
v005=find(ismember(v001(:,1),varargin(1)));if isempty(v005),error('HJW:noToolbox:InvalidInput',...
['The functionname parameter must be the ','name of an implemented function.']),...
end,if strcmpi(v000,'doc'),doc(sprintf('%s>%s',v003,v001{v005,...
2})),else,help(sprintf('%s>%s',v003,v001{v005,2})),end,return,end,v005=find(ismember(v001(:,...
1),{lower(v000)}));if isempty(v005),error('HJW:noToolbox:InvalidInput',...
['The functionname parameter must be the ','name of an implemented function.']),...
else,v006=v001{v005,3};end,varargout=cell(nargout,1);[varargout{:}]=v006(varargin{:});end
function [v000,v001] = RegGrow(v002,v003,v004,varargin),if ...
nargout>2,error('HJW:RegGrow:nargout','Incorrect number of output arguments.'),end,v005=cell(1,...
nargin);if nargin>=1,v005{1}=v002;end,if nargin>=2,v005{2}=v003;end,if nargin>=3,v005{3}=v004;
end,if nargin>=4,v005(4:end)=varargin;end,[v006,v007,v008]=f07(v005{:});if ~v006,rethrow(v008),...
else,[v002,v003,v004,v009,v010,v011]=deal( v007.I,v007.maxDiff,v007.seed,v007.kernel,...
v007.overflow,v007.waitbar);end,if v011.bar,v011.h=waitbar(0,'processing...');v011.count=0;
v012=0;end,v001=false;v000=logical(zeros(size(v002)));v004=num2cell(v004);v000(v004{:})=true;
v013=true;while v013,if v010.check,if sum(v000(:))>v010.volume,v001=true;break,end,end,if ...
v011.bar,v011.count=v011.count+1;if mod(v011.count,5)==0,try v012=mod(v012+1,3);v014=repmat('.',...
1,1+v012);v015=sprintf('processing%s',v014);waitbar(sum(v000(:))/v011.max,v011.h,v015),...
drawnow,catch,v008=lasterror;if strcmp(v008.identifier,'MATLAB:waitbar:InvalidSecondInput'),...
v011.bar=false;end,end,end,end,v016= ~v000 & convn(v000,v009,...
'same')>0.5;v017=abs(v002-mean(v002(v000)))>=v003;v016(v017)=false;if any(v016(:)),v000(v016)=true;
else,v013=false;end,end,if v011.bar,delete(v011.h),end,if nargout==0,clear v000,end,end
function v000=f00(v001,v002),if nargin>2 ...
|| nargin==0 || nargout~=1,error('HJW:bwlabel:Nargs','Incorrect number of inputs or outputs'),...
end,if ~isa(v001,'logical'),error('HJW:bwlabel:InvalidInput_IM','IM must be a logical array'),...
end,v003=ndims(v001);if nargin==1,v002=(3^v003-1);end,if isa(v002,'logical'),...
v004=v002;if any(size(v004)~=3) || ndims(v004)~=v003,error('HJW:bwlabel:InvalidInput_SE_conn',...
'SE must be a 3x...x3 logical array with the same number of dimensions as IM.'),...
end,else,if numel(v002)~=1 || round(v002)~=v002,error('HJW:bwlabel:InvalidInput_SE_conn',...
'The second input must be a valid connectivity or a 3x...x3 logical array.'),...
end,if v002==(2*v003),...
v005=cell(v003,1);[v005{:}]=ndgrid(-1:1);v005=cat(v003+1,v005{:});v005=sqrt(sum(v005.^2,...
v003+1));v004=v005<=1;elseif v002==(3^v003-1),v004=logical(ones(3*ones(1,v003)));
else,error('HJW:bwlabel:conn_not_supported',['The provided connectivity is not supported. ',...
char(10),'The value of conn must be either the minimal or maximal connectivity',...
char(10),'for the size of IM. You can provide a custom SE instead.']),...
end,end,v000=zeros(size(v001));v006=cell(v003,1);v007=0;v008=[];while any(v001(:)),...
[v006{:}]=findND(v001,1);v009=[v006{:}];v010=RegGrow(v001,0.5,v009,'kernel',v004,'waitbar',false);
v007=v007+1;v000(v010)=v007;v001(v010)=false;v008(v007)=sum(v010(:));end,[v011,v012]=sort(v008);
v012=(numel(v008)+1)-v012;v000=-v000;for v013=1:v007,v000(v000==-v013)=v012(v013);end,end
function ...
[v000,v001,v002,v003]=f01(v000,v001,v004),v002=false;v003=struct;if ~isa(v000,'logical'),...
v003.identifier=sprintf('HJW:%s:InvalidInput_IM',v004);v003.message='IM must be a logical array.';
return,end,if ~isa(v001,'logical'),v005=v001;if ~isnumeric(v005) ...
|| numel(v005)~=1 || v005<=0,v003.identifier=sprintf('HJW:%s:InvalidInput_SE_r',...
v004);v003.message=['The second input must be a ','scalar (radius), or logical (SE).'];
return,end,v006=ceil(v005);v007=ndims(v000);v008=cell(v007,1);[v008{:}]=ndgrid(-v006:v006);
v008=cat(v007+1,v008{:});v008=sqrt(sum(v008.^2,v007+1));v001=v008<=v005;
else,if ndims(v001)~=ndims(v000),v003.identifier=sprintf('HJW:%s:InvalidInput_SE',...
v004);v003.message='SE must have the same number of dimensions as IM.';return,end,end,v002=true;end
function v000=f02(v001,v002),if nargin~=2 || nargout~=1,error('HJW:imclose:Nargs',...
'Incorrect number of inputs or outputs'),end,[v001,v002,v003,v004]=f01(v001,v002,...
'imclose');if ~v003,rethrow(v004),end,persistent v005,if isempty(v005),v005=f09('>=','R2017a',...
'Octave','<',0);end,if v005,v006=ceil(size(v002)/2);v007=logical(zeros(size(v001)+2*v006));
v008=ndims(v001);switch v008,case 2,v007((1+v006(1)):(end-v006(1)),...
(1+v006(2)):(end-v006(2)))=v001;case 3,v007((1+v006(1)):(end-v006(1)),(1+v006(2)):(end-v006(2)),...
(1+v006(3)):(end-v006(3)))=v001;otherwise,v009=cell(v008,1);v010=v009;for ...
v011=1:v008,v010{v011}=logical([zeros(1,v006(v011)) ones(1,size(v001,v011)) zeros(1,v006(v011))]);
end,[v009{:}]=ndgrid(v010{:});v009=cat(v008+1,v009{:});v009=all(v009,v008+1);v007(v009)=v001;
end,v001=v007;end,v000=f04(f03(v001,v002),v002);if v005,switch ndims(v001),case 2,v000=v000( (1+...
v006(1)):(end-v006(1)),(1+v006(2)):(end-v006(2)));case 3,v000=v000( (1+v006(1)):(end-v006(1)),...
(1+v006(2)):(end-v006(2)),(1+v006(3)):(end-v006(3)));otherwise,v000=v000(v009);end,end,end
function v000=f03(v001,v002),if nargin~=2 || ...
nargout~=1,error('HJW:imdilate:Nargs','Incorrect number of inputs or outputs'),end,[v001,v002,...
v003,v004]=f01(v001,v002,'imdilate');if ~v003,rethrow(v004),end,v000=convn(v001,v002,'same')>0;end
function v000=f04(v001,v002),...
if nargin~=2 || nargout~=1,error('HJW:imerode:Nargs','Incorrect number of inputs or outputs'),...
end,[v001,v002,v003,v004]=f01(v001,v002,'imerode');if ~v003,rethrow(v004),...
end,v005=convn(v001,v002,'same');v006=convn(ones(size(v001)),v002,'same');v000= v005 >= v006;end
function v000=f05(v001,v002),if nargin~=2 ...
|| nargout>1,error('HJW:imopen:Nargs','Incorrect number of inputs or outputs'),end,[v001,v002,...
v003,v004]=f01(v001,v002,'imopen');if ~v003,rethrow(v004),end,v000=f03(f04(v001,v002),v002);end
function v000=f06(v000,...
varargin),if ~isequal(varargin,{'NoFigure',true}),error('HJW:montage:NotImplemented',...
'Only one syntax is supported, see documentation.'),end,v001=eval([class(v000) '(0)']);
if issparse(v000),v001=sparse(v001); end,if ~isreal(v000), v001=complex(v001);
end,v001=repmat(v001, size(v000,1),size(v000,2),size(v000,3) );v002=ceil(sqrt(size(v000,4)));
v003=ceil(size(v000,4)/v002);v000=mat2cell(v000,size(v000,1),size(v000,...
2),size(v000,3),ones(1,size(v000,4)));v000=squeeze(v000);if numel(v000)<(v002*v003),v000( ...
(numel(v000)+1) : (v002*v003))={v001};end,v000=reshape(v000,v002,v003)';v000=cell2mat(v000);end
function [v000,v001,v002]=f07(varargin),v000=false;v001=struct;v002=struct('identifier',...
'','message','');persistent v003,if isempty(v003),v004=exist('OCTAVE_VERSION',...
 'builtin');if v004,v005 = get(0,'DefaultImageCData');v005 = v005/max(v005(:));
else,v006 = pow2(get(0,'DefaultImageCData'),47);v005 = bitshift(v006,-37);v005 = fix(v005);
v005 = bitand(v005,31);v005 = v005/max(v005(:));end,v003.I=v005;v003.maxDiff=[];v003.seed=[];
v003.nHood=[];v003.kernel=[];v003.shifts=[];v003.overflow=[];v003.KeepShellVoxelsOnly=false;
v003.waitbar=[];end,if nargin==0,v001=f08(v003);v001.seed=[2 2];v000=true;return,...
end,v007=0;v008={'I','maxDiff','seed'};while numel(varargin)>0 && ~(isa(varargin{1},'char') ...
|| isa(varargin{1},'struct')),v007=v007+1;v001.(v008{v007})=varargin{1};varargin(1)=[];end,...
if numel(varargin)>0,if isa(varargin{1},'struct'),v009=v001;v001=varargin{1};v010=fieldnames(v009);
for v011=1:numel(v010),v001.(v010{v011})=v009.(v010{v011});end,else,try while numel(varargin)>0,...
v001.(varargin{1})=varargin{2};varargin(1:2)=[];end,catch,v002.identifier='HJW:RegGrow:ParseError';
v002.message='Parsing of the inputs failed.';return,end,end,end,v010=fieldnames(v001);[v012,...
v013]=sort(lower(v010));v010=v010(v013);for v014=1:numel(v010),v015=v010{v014};v016=v001.(v015);
v002.identifier=['HJW:RegGrow:incorrect_input_opt_' lower(v015)];switch lower(v015),case ...
'i',try v016=double(v016);catch,v002.message='The image input must be convertible to a double .';
return,end,v001.I=v016;case 'kernel',try if isfield(v001,'I'),v005=v001.I;
else,v005=v003.I;end,v016=logical(v016);if ndims(v016)~=ndims(v005) || ~all(mod(size(v016),...
2)==1),disp(v016(-1)),end,v017=v016;v018=(size(v017)+1)/2;v018=num2cell(v018);v017(v018{:})=false;
v016(v018{:})=true;if sum(v017(:))==0,warning('HJW:RegGrow:CentroidOnlyKernel',...
['Only the centroid of the kernel is marked as true.\n',...
'This will cause the result to only contain the seed position.']),...
end,catch,v002.message=['The kernel input must a logical with ndims(IM) dimensions.',...
char(10),'Each dimension must be an odd length.'];return,end,v001.kernel=v016;
case 'maxdiff',try v016=double(v016);catch,end,if ~isa(v016,'double') || numel(v016)~=1 || ...
v016<0 || isnan(v016),v002.message=['The maxDiff input must be a non-negative numeric scalar.',...
char(10),'It must also be convertible to a double.'];return,end,v001.maxDiff=v016;case 'nhood',...
if any(ismember(lower(v010),{'kernel','shifts'})),continue,end,if isfield(v001,'I'),v005=v001.I;
else,v005=v003.I;end,v019=['The nHood parameter must be either the maximal or minimal value, ',...
char(10),'so 4 or 8 for 2D images, 6 or 26 for 3D images, etc.',...
char(10),'This parameter must be set if the kernel parameter is not set.'];
try v020=ndims(v005);if ~ismember(v016,...
[3^v020-1 2*v020]),v002.message=v019;return,end,catch,v002.message=v019;return,end,v001.nHood=v016;
case 'overflow',if ~isempty(v016) && ( ~isnumeric(v016) || numel(v016)>1 || abs(v016-...
round(v016))>eps),v002.message='The overflow volume must be empty, or a numeric scalar integer.';
return,end,v001.overflow=v016;case 'seed',if isfield(v001,'I'),v005=v001.I;else,v005=v003.I;
end,v019='The seed must be a vector with a valid position.';try if numel(v016)~=ndims(v005),...
v002.message=v019;return,end,v021=v016(:)';v022=num2cell(v021);v023=v005(v022{:});
catch,v002.message=v019;return,end,v002.message='';v001.seed=v021;case {'waitbar',...
'silent'},[v024,v016]=f11(v016);if ~v024,v002.message='waitbar must be a logical scalar.';
return,end,if strcmpi(v015,'silent'),v016= ~v016;end,v001.waitbar=v016;case ...
'shifts',v001.shifts=v016;otherwise,v002.message=sprintf('Name,Value pair not recognized: %s',...
v015);v002.identifier='HJW:RegGrow:incorrect_input_NameValue';
return,end,end,v010=fieldnames(v003);for v014=1:numel(v010),if ...
~isfield(v001,v010(v014)),v001.(v010{v014})=v003.(v010{v014});end,end,v000=true;v001=f08(v001);end
function v000=f08(v000),if isempty(v000.seed),v001=ones(ndims(v000.I),1);
v000.seed=v001;end,if isempty(v000.maxDiff),v000.maxDiff=2/3*std(v000.I(:));end,...
if ~isempty(v000.shifts),v002=v000.shifts;v003=max(abs(v002(:)));v004=zeros( (1+2*v003)*ones(1,...
size(v002,2)) );v002=v002+v003+1;v002=mat2cell(v002,size(v002,1),ones(1,size(v002,2)));
v005=sub2ind(size(v004),v002{:});v004(v005)=1;v004=logical(v004);v000.kernel=v004;end,if ...
isempty(v000.kernel),v006=v000.nHood;if isempty(v006),v006=2*ndims(v000.I);end,v007=ndims(v000.I);
if v006==(2*v007),v008=cell(v007,1);[v008{:}]=ndgrid(-1:1);v008=cat(v007+1,v008{:});
v008=sqrt(sum(v008.^2,v007+1));v004=v008<=1;elseif v006==(3^v007-1),v004=logical(ones(3*ones(1,...
v007)));end,v000.kernel=v004;end,v009=v000.overflow;if isempty(v009) || isnan(v009) || ...
isinf(v009),v000.overflow=struct('check',false,'volume',inf);else,v000.overflow=struct('check',...
true,'volume',v009);end,if isempty(v000.waitbar),v000.waitbar= ndims(v000.I)~=2 ;
end,v000.waitbar=struct('bar',v000.waitbar,'max',min(numel(v000.I),v000.overflow.volume));end
function v000=f09(v001,...
v002,v003,v004,v005),persistent v006 v007 v008,if isempty(v006),v008=exist('OCTAVE_VERSION',...
 'builtin');v006=version;v009=strfind(v006,'.');if numel(v009)~=1,v006(v009(2):end)='';
v009=v009(1);end,v006=[str2double(v006(1:(v009-1))) str2double(v006((v009+1):end))];
v006=v006(1)+v006(2)/100;v006=round(100*v006);v007={ 'R13' 605;'R13SP1' 605;'R13SP2' 605;
'R14' 700;'R14SP1' 700;'R14SP2' 700;'R14SP3' 701;'R2006a' 702;'R2006b' 703;'R2007a' 704;
'R2007b' 705;'R2008a' 706;'R2008b' 707;'R2009a' 708;'R2009b' 709;'R2010a' 710;'R2010b' 711;
'R2011a' 712;'R2011b' 713;'R2012a' 714;'R2012b' 800;'R2013a' 801;'R2013b' 802;'R2014a' 803;
'R2014b' 804;'R2015a' 805;'R2015b' 806;'R2016a' 900;'R2016b' 901;'R2017a' 902;'R2017b' 903;
'R2018a' 904;'R2018b' 905;'R2019a' 906;'R2019b' 907;'R2020a' 908;'R2020b' 909};end,if v008,...
if nargin==2,warning('HJW:ifversion:NoOctaveTest',['No version test for Octave was provided.',...
char(10),'This function might return an unexpected outcome.']),...
if isnumeric(v002),v010=0.1*v002+0.9*fix(v002);v010=round(100*v010);
else,v011=ismember(v007(:,1),v002);if sum(v011)~=1,warning('HJW:ifversion:NotInDict',...
'The requested version is not in the hard-coded list.'),v000=NaN;return,else,...
v010=v007{v011,2};end,end,elseif nargin==4,[v001,v010]=deal(v003,v004);v010=0.1*v010+0.9*fix(v010);
v010=round(100*v010);else,[v001,v010]=deal(v004,v005);v010=0.1*v010+0.9*fix(v010);
v010=round(100*v010);end,else,if isnumeric(v002),v010=0.1*v002+0.9*fix(v002);v010=round(100*v010);
else,v011=ismember(v007(:,1),v002);if sum(v011)~=1,warning('HJW:ifversion:NotInDict',...
'The requested version is not in the hard-coded list.'),v000=NaN;return,else,v010=v007{v011,2};
end,end,end,switch v001,case '==', v000= v006 == v010;case '<' , v000= v006 < v010;
case '<=', v000= v006 <= v010;case '>' , v000= v006 > v010;case '>=', v000= v006 >= v010;end,end
function [v000,v001]=f10(v002,v003),if (nargin<1 || nargin>2) ...
|| (nargout<1 || nargout>2),error('HJW:bwdist:Nargs','Incorrect number of inputs or outputs'),...
end,if nargout==2,v001=[];error('HJW:bwdist:NotImplemented',...
'The second output is not implemented'),end,if ~isa(v002,'logical'),...
try if isa(v002,'char') || isa(v002,'string'),error('trigger error'),end,v002=logical(v002);
if ~isa(v002,'logical'),v002=gather(v002);end,catch,error('HJW:bwdist:InvalidIM',...
'The image input is not logical or numeric'),end,end,if nargin==1,v003='euclidean';
else,if ~(isa(v003,'char') || isa(v003,'string')),error('HJW:bwdist:InvalidMethod',...
'The distance method input is not a char'),end,v003=lower(char(v003));if ~ismember(v003,...
{'euclidean','chessboard','cityblock','quasi-euclidean'}),error('HJW:bwdist:InvalidMethod',...
['Distance method input not recognised, should one of the following:\n',...
'euclidean, chessboard, cityblock, quasi-euclidean']),...
end,end,v004=cell(1,ndims(v002));v005=1+ndims(v002);for ...
v006=1:numel(v004),v004{v006}=0:(size(v002,v006)-1);end,[v004{:}]=ndgrid(v004{:});v004=cat(v005,...
v004{:});switch v003,case 'euclidean',v004=sqrt(sum(v004.^2,v005));case 'chessboard',...
v004=max(abs(v004),[],v005);case 'cityblock',v004=sum(abs(v004),v005);case 'quasi-euclidean',...
error('HJW:bwdist:NotImplemented','The quasi-euclidean metric is not implemented.'),otherwise,...
v007={'If you''re seeing this, the code is in what';'I thought was an unreachable state.';
' ';'I could give you advice for what to do,';'but honestly, why should you trust me?';
'I clearly screwed this up. I''m writing a';'message that should never appear, yet';
'I know it will probably appear someday.';' ';'On a deep level, I know I''m not';
'up to this task. I''m so sorry.';' ';'(c) R. Munroe, xkcd.com/2200'};uiwait(errordlg(v007)),...
error('unreachable line'),end,v008=unique(v004(:));v009=true(3*ones(1,ndims(v002)));v010=v002 & ...
~f04(v002,v009);v011=4;v012=cell(1,ndims(v002));v004=cell(1,ndims(v002));for v006=1:numel(v004),...
v012{v006}=floor(-size(v002,v006)/v011):ceil(size(v002,v006)/v011);v004{v006}=(-size(v002,...
v006)):(size(v002,v006));end,[v004{:}]=ndgrid(v004{:});[v012{:}]=ndgrid(v012{:});
v004=cat(v005,v004{:});v012=cat(v005,v012{:});switch v003,case 'euclidean',...
v004=sqrt(sum(v004.^2,v005));v012=sqrt(sum(v012.^2,v005));v013=floor(min(size(v012)/2));
case 'chessboard',v004=max(abs(v004),[],v005);v012=max(abs(v012),[],v005);
v013=floor(min(size(v012)/2));case 'cityblock',v004=sum(abs(v004),v005);v012=sum(abs(v012),...
v005);v013=floor(min(size(v012)/2));case 'quasi-euclidean',error('HJW:bwdist:NotImplemented',...
'The quasi-euclidean metric is not implemented.'),otherwise,error('unreachable line'),...
end,v014=v012;v015=v002;v000=single(zeros(size(v002)));for v016=2:numel(v008),...
if ~any(~v015(:)),break,end,if v008(v016)>=v013,v014=v004;end,v017=v014<=v008(v016);
v017=convn(v010,v017,'same');v018= ~v015 & v017 ;v000(v018)=v008(v016);v015=v017 | v015;end,end
function [v000,v001]=f11(v001),...
persistent v002,if isempty(v002),v002={true,false;1,0;'on','off';'enable','disable';
'enabled','disabled'};try v002(end+1,:)=eval('{"on","off"}');catch,end,end,v000=true;try if ...
isa(v001,'char') || isa(v001,'string'),try v001=lower(v001);catch,end,end,for v003=1:size(v002,...
1),for v004=1:2,if isequal(v001,v002{v003,v004}),v001=v002{1,v004};return,end,end,end,...
if isa(v001,'matlab.lang.OnOffSwitchState'),v001=logical(v001);return,end,catch,end,v000=false;end