Manual lung segmentation takes about 10 minutes and it requires a certain skill to get the same high-quality result as with automatic segmentation. Automatic segmentation takes about 15 seconds.

I assumed that without a neural network it would be possible to get an accuracy of no more than 70%. I also assumed, that morphological operations are only the preparation of an image for more complex algorithms. But as a result of processing of those, although few, 40 samples of tomographic data on hand, the algorithm segmented the lungs without errors. Moreover, after testing in the first five cases, the algorithm didn’t change significantly and correctly worked on the other 35 studies without changing the settings.

Also, neural networks have a disadvantage — for their training we need hundreds of training samples of lungs, which need to be marked up manually.

Content

Structure of the respiratory system
Haunsfield scale
Mathematical morphology
Lee algorithm and RLE compression
Threshold transform of the base volume
External air removal
Segmentation of the maximum volume object
Closing the vessels inside the lungs
Respiratory organs segmentation algorithm
Algorithm implementation in the MATLAB
Conclusion

Structure of the respiratory system

The respiratory system includes the airways and lungs. The airways are divided into the upper and lower airways. The point of separation between the lower and upper respiratory tract is the point of intersection of the esophagus and airways. All the organs that are above the larynx are the upper tract and the larynx and organs below are the lower tract.

List of the respiratory organs:

The nasal cavity is an air-filled cavity between the nose and the pharynx.
The pharynx is a channel through which food and air travels.
The trachea is a tube connected to the larynx and bronchi.
The larynx is responsible for the voice forming. It is at the level of the cervical vertebrae C4-C6.
The bronchi are the respiratory channels, the main part of which is located inside the lungs.
The lungs are the main respiratory organ.

Haunsfield scale

Godfrey Hounsfield was an English electrical engineer who, together with American theorist Allan Cormack, developed computed tomography, for which he shared the Nobel Prize in 1979.

The Hounsfield scale — is a quantitative scale for describing radiodensity, which is measured in Hounsfield units, denoted as HU.

Radiodensity is calculated on the basis of the radiation attenuation coefficient. Attenuation coefficient characterizes how easily a volume of material can be penetrated by the radiation.

Radiodensity is calculated by the formula:

${μ_{X}-μ_{water} \over μ_{water}-μ_{air}} \times 1000$

where $μ_X, μ_{water}, μ_{air}$ are linear attenuation coefficients for the measured material, water and air.

Radiodensity can be negative because zero radiodensity corresponds to water. This means that all materials having a lower attenuation coefficient than water attenuation coefficient (for example, lung tissue, air) will have a negative radiodensity.

The approximate radiodensities for various tissues are listed below:

Air: -1000 HU.
Respiratory organs: from -950 to -300 HU.
Blood (without contrasting the vessels): from 0 to 100 HU.
Bones: 100 to 1000 HU.

Links to Wikipedia: Hounsfield scale, Godfrey Hounsfield, attenuation coefficient.

Mathematical morphology

Morphological operations are the main tools which were used in the algorithm.
In the field of computer vision, morphological operations are called a group of algorithms for the transformation of the shape of objects. Most often, morphological operations are applied to binarized images where units correspond to object voxels, and zeros to empty.
The main morphological operations include:

Morphological dilation is the adding new voxels to all edge voxels of the objects. In other words, algorithm scans all the voxels corresponding to the edge of the object and enlarge them using a kernel with a given form (sphere, cube, cross etc). The kernel has its own radius (for the sphere), or the width of the side (for the cube). This operation is often used to merge multiple nearby objects into a single object.

Morphological erosion is the cut-off of all voxels lying on the border (boundary) of the objects. This operation is reverse to dilation. This operation is useful for removing noise which has the form of many small interconnected objects. However, this method of noise removal should be only used if the object of interest has a thickness much greater than the radius of erosion. Otherwise, the useful information will be lost.

Morphological closing is a dilation followed by an erosion. It is used to close holes inside objects and to merge nearby objects.

Morphological opening is an erosion followed by a dilation. It is used to remove small noise objects and to separate the objects into several objects.

Lee algorithm and RLE compression

To select objects in binarized voxel volume, the Lie algorithm is used. Originally this algorithm was invented to find the shortest way out of the labyrinth. We use it to select the object and move it from one three-dimensional array of voxels to another. The basic idea of the algorithm consists in parallel movement in all possible directions from the starting point. For a three-dimensional case, it is possible to move in 26 or 6 directions from one voxel.

To optimize the performance, the run-length encoding algorithm was applied. The main idea of the algorithm is that the sequence of ones and zeros is replaced by a digit equal to the number of elements in the sequence. For example, the string “00110111” can be replaced by: “2; 2; 1; 3”. This reduces the number of memory accesses.

Wikipedia links: Lee algorithm, RLE algorithm.

Threshold transform of the base volume

Using the tomograph there were obtained the data about an radiodensity at each point of the scanned space. Air voxels have an radiodensity in the range from -1100 to -900 HU, and respiratory organs voxels have radiodensity from -900 to -300 HU. Therefore, we can remove all unnecessary voxels with radiodensity greater than -300 HU. As a result, we obtain a binarized voxel volume containing only the respiratory organs and air.

External air removal

To segment the internal air of the body, we delete all objects that are adjacent to the corners of the 3D volume. Thus, we remove the external air.

However, not in all cases the air inside the tomograph table will be removed, because it may not have a connection with the corners of the 3D volume.

Therefore, we will scan not only the corners, but also all the voxels lying on any of the edge planes of the f. But, as we can see in the image below, the lungs were also removed. It turns out that the trachea also had a connection with the upper plane of the 3D volume.

We'll have to exclude the top plane from the scan area. There are also studies in which the lungs were not fully captured and the lower plane is connected to the lungs.
So if you wish, you can exclude the lower plane too.

But this method only affects chest studies. In the case of the capture of the full volume of the body, the connection of internal and external air through the nasal cavity will appear Therefore, it is necessary to apply morphological erosion to separate the internal and external air.

After applying erosion, we can return to the previously obtained method of segmentation of external air which is based on connectivity of the external air with side planes of the 3D-volume.

After the external air segmentation, we could immediately subtract the external air from the whole volume of the air and lungs and to get the internal air of the body and lungs. But there is one problem. After erosion, some of the information about the external air was lost. To restore it, we apply dilatation of external air.

Next, we subtract the external air from the whole air and respiratory organs and get the internal air and respiratory organs.

Segmentation of the maximum volume object

Next, we segment the respiratory organs as the maximum object in volume. The respiratory organs have no connection with the air inside the gastrointestinal tract.

It is worth noting that the correct choice of the radiodensity threshold at the initial step is important. Otherwise, in some cases, there may be no connection between the two lungs as a result of low resolution. For example, if we assume that voxels of the respiratory organs have radiodensity from -500 HU and less, then in the case below, segmentation of the respiratory organs as the largest object in volume will lead to an error, because there is no connection between the two lungs. Therefore, the threshold should be increased to -300 HU.

Closing the vessels inside the lungs

To capture the vessels inside the lungs we will use morphological closing, that is, dilation followed by erosion with the same radius. Radiodensity of vessels (without contrasting) is about 0..100 HU.

In the image we can see that large blood channels are not closed. But this is not necessary. The purpose of this operation was to destroy the many small holes inside the lungs to simplify lung segmentation process in the next steps.

Respiratory organs segmentation algorithm

As a result, we obtain the following respiratory organs segmentation algorithm:

Threshold transform of the base volume with the threshold “< -300 HU”.
Morphological erosion with the 3 mm radius for the separation of external and internal air.
Segmentation of external air based on external air connectivity with the boundary side planes of the 3D volume.
Morphological dilation of the external air to restore the information lost after the erosion.
Subtraction of the external air from the whole air and respiratory organs to obtain the internal air and respiratory organs.
Segmentation of the maximum volume object.
Morphological closing of the vessels inside the lungs.

Algorithm implementation in MATLAB

Function "getRespiratoryOrgans"

% Returns the whole respiratory organs volume (lung and airway volume)
% without separating of the left and right lung.

% V = base volume with radiodensity data in Hounsfield units.
% cr = radius of vessels morphological closing.
% ci = iteration count of vessels morphological closing (f.e. 3-times
% make dilation and after that 3-times make erosion.

function RO = getRespiratoryOrgans(V,cr,ci)
% Threshold transform of the base volume with the threshold level "<-300 HU".
AL=~imbinarize(V,-300);
% Morphological erosion with the 3 mm radius for the separation of external
% and internal air.
SE=strel('sphere',3);
EAL=imerode(AL,SE);
% Segmentation of external air based on external air connectivity with the
% boundary side planes of the 3D volume.
EA=getExternalAir(EAL);
% Morphological dilation of the external air to restore the information lost
% after the erosion.
DEA=EA;
for i=1:4
    DEA=imdilate(DEA,SE);
    DEA=DEA&AL;
end
% Subtraction of the external air from the whole air and respiratory organs to
% obtain the internal air and respiratory organs.
IAL=AL-DEA;
% Segmentation of the maximum volume object.
RO=getMaxObject(IAL);
Morphological closing of the vessels inside the lungs.
RO=closeVoxelVolume(RO,3,2);

Function "getExternalAir"

% Returns the volume which is connected with the edge surfaces of
% the voxel scene (except the top surface, because lungs can have
% a connection with the top surface).

% EAL = eroded lung-and-air binarized volume.

function EA = getExternalAir(EAL)
% The “bwlabeln” function segments objects: voxels of a one object
% equates   to 1,  and of another one  - to 2, etc.
V=bwlabeln(EAL);
% We request such characteristics of the objects as bounding box and
% list of all object voxels.
R=regionprops3(V,'BoundingBox','VoxelList');
n=height(R);
% Create a 3D matrix for storing external air voxels.
s=size(EAL);
EA=zeros(s,'logical');
% Scan all the objects to find objects belonging to the external air.
for i=1:n
    % Define the minimum and maximum x and y coordinates of the object.
    x0=R(i,1).BoundingBox(1);
    y0=R(i,1).BoundingBox(2);
    x1=x0+R(i,1).BoundingBox(4);
    y1=y0+R(i,1).BoundingBox(5);
   % If the edge voxels of the object are in contact with the side planes of
   % the 3D volume, then copy all the voxels of this object into the matrix EA.
    if (x0 < 1 || x1 > s(1)-1 || y0 < 1 || y1 > s(2)-1)
        % Convert the object voxel coordinates data to the matrix type:
        % [[x1 y1 z1] [x2 y2 z3] ... [xn yn zn]].
        mat=cell2mat(R(i,2).VoxelList);
        ms=size(mat);
        % Fill the matrix containing the voxels of the external air.
        for j=1:ms(1)
            x=mat(j,2);
            y=mat(j,1);
            z=mat(j,3);
            EA(x,y,z)=1;
        end
    end
end

Function "getMaxObject"

% Returns the largest object in the voxel volume V.

% O = the voxels of the largest object.
% m = the volume of the largest object.

function [O,m] = getMaxObject(V)
% Segment the objects. 
V=bwlabeln(V);
% Reqiest the information about the objects volume and objects voxel coordinates.
R=regionprops3(V,'Volume','VoxelList');
% Determine the index of the maximum volume object.
v=R(:,1).Volume;
[m,i]=max(v);
% Create the 3D matrix for storing the largest object voxels.
s=size(V);
O=zeros(s,'logical');
% Copy the largest object voxels to the new matrix.
mat=cell2mat(R(i,2).VoxelList);
ms=size(mat);
for j=1:ms(1)
    x=mat(j,2);
    y=mat(j,1);
    z=mat(j,3);
    O(x,y,z)=1;
end

The source code can be downloaded by the reference.

Conclusion

The following articles are going to be:

trachea and bronchi segmentation;
lung segmentation;
lung lobes segmentation.

The following algorithms will be considered:

distance transform;
nearest neighbor transform, also known as feature transform;
calculation of the eigenvalues of the Hessian matrix for segmentation of flat 3D objects;
watershed segmentation.

Automatic respiratory organ segmentation