This tutorial is automatically generated from the file lung/test/tutorials/TestSimpleImpedanceProblemTutorial.hpp at revision 196f6e705993/git_repo. Note that the code is given in full at the bottom of the page.

An example showing how to calculate transfer impedance of an airway tree using a simple impedance model

In this tutorial we demonstrate the use of SimpleImpedanceProblem to calculate transfer impedance on an airway tree model. We further demonstrate post-processing of the output using ImpedancePostProcessor to calculate a number of clinically relevant measures.

Note that SimpleImpedanceProblem uses Poiseuille formulas to calculate impedance rather than more accurate acoustic impedance equations. For the more accurate version see ImpedanceProblem.

The usual headers are included 

#include <cxxtest/TestSuite.h>
#include "TrianglesMeshReader.hpp"

SimpleImpedanceProblem does most of the work in calculating impedance. 

#include "SimpleImpedanceProblem.hpp"

ImpedancePostProcessor allows easy calculation of a number of clinically relevant measures. 

#include "ImpedancePostProcessor.hpp"

Define the test 

class TestSimpleImpedanceProblemTutorial : public CxxTest::TestSuite
{
public: // Tests should be public!

    void TestCalculateImpedance()
    {
        EXIT_IF_PARALLEL;

First, we load up a mesh containing the centre lines and radii of the a complete conducting airway tree. The mesh will typically have been developed using a combination of computed tomography (CT) image segmentation and algorithmic airway generation. 

        TetrahedralMesh<1,3> mesh;
        TrianglesMeshReader<1,3> mesh_reader("lung/test/data/TestSubject002");
        mesh.ConstructFromMeshReader(mesh_reader);

Note that the mesh defined above was developed using a CT scan taken at full inspiration. Impedance is more commonly recorded during tidal breathing. Here we use a simple scaling to bring the airway radii down into the tidal breathing range. 

        for (TetrahedralMesh<1,3>::NodeIterator node_iter = mesh.GetNodeIteratorBegin();
             node_iter != mesh.GetNodeIteratorEnd();
             ++node_iter)
        {
            node_iter->rGetNodeAttributes()[0] *= 0.7;
        }

Setup a SimpleImpedanceProblem and tell it that the given mesh is defined in millimetres 

        SimpleImpedanceProblem problem(mesh, 0u);
        problem.SetMeshInMilliMetres();

This vector lists the input frequencies at which to calculate impedance. They must be monotonically increasing. 

        std::vector<double> test_frequencies;
        test_frequencies.push_back(1.0);
        test_frequencies.push_back(2.0);
        test_frequencies.push_back(3.0);
        test_frequencies.push_back(5.0);
        test_frequencies.push_back(10.0);
        test_frequencies.push_back(20.0);
        test_frequencies.push_back(30.0);
        problem.SetFrequencies(test_frequencies);               //Set & get frequencies for coverage

The simple impedance model defines a linear spring at each terminal of the airway tree. This method allows us to set the elastance of the whole lung (in Pa/m3). This elastance is then evenly distributed over the terminals. 

        problem.SetElastance(5.8*98.0665*1e3);

Calculates the impedance at the given frequencies 

        problem.Solve();

Get the calculated impedances. The impedance at each frequency is made up of a real component (the resistance) and a complex component (the elastance). 

        std::vector<std::complex<double> > impedances = problem.rGetImpedances();

The impedances calculated above could at this stage be written to a file and plotted. Instead, we make use of ImpedancePostProcessor to calculate a number of common clinical summary statistics from the data. 

        ImpedancePostProcessor processor(test_frequencies, impedances);

        std::cout << "\n";
        std::cout << "R5       = " << real(impedances[3])*1e-6 << " kPa.s.L^-1" << std::endl;
        std::cout << "R20      = " << real(impedances[5])*1e-6 << " kPa.s.L^-1" << std::endl;
        std::cout << "R5 - R20 = " << processor.GetR5MinusR20()*1e-6 << " kPa.s.L^-1" << std::endl;
        std::cout << "Rrs      = " << processor.GetRrs()*1e-6 << " kPa.s.L^-1" << std::endl;
        std::cout << "\n";
        std::cout << "X5       = " << imag(impedances[5])*1e-6 << " kPa.s.L^-1" << std::endl;
        std::cout << "Fres     = " << processor.GetResonantFrequency() << " Hz " << std::endl;
        std::cout << "Ax       = " << processor.GetAx()*1e-6 << " kPa.L^-1" << std::endl;
        std::cout << "\n";
    }
};

Code

The full code is given below

File name TestSimpleImpedanceProblemTutorial.hpp

#include <cxxtest/TestSuite.h>
#include "TrianglesMeshReader.hpp"

#include "SimpleImpedanceProblem.hpp"

#include "ImpedancePostProcessor.hpp"

class TestSimpleImpedanceProblemTutorial : public CxxTest::TestSuite
{
public: // Tests should be public!

    void TestCalculateImpedance()
    {
        EXIT_IF_PARALLEL;

        TetrahedralMesh<1,3> mesh;
        TrianglesMeshReader<1,3> mesh_reader("lung/test/data/TestSubject002");
        mesh.ConstructFromMeshReader(mesh_reader);

        for (TetrahedralMesh<1,3>::NodeIterator node_iter = mesh.GetNodeIteratorBegin();
             node_iter != mesh.GetNodeIteratorEnd();
             ++node_iter)
        {
            node_iter->rGetNodeAttributes()[0] *= 0.7;
        }

        SimpleImpedanceProblem problem(mesh, 0u);
        problem.SetMeshInMilliMetres();

        std::vector<double> test_frequencies;
        test_frequencies.push_back(1.0);
        test_frequencies.push_back(2.0);
        test_frequencies.push_back(3.0);
        test_frequencies.push_back(5.0);
        test_frequencies.push_back(10.0);
        test_frequencies.push_back(20.0);
        test_frequencies.push_back(30.0);
        problem.SetFrequencies(test_frequencies);               //Set & get frequencies for coverage

        problem.SetElastance(5.8*98.0665*1e3);

        problem.Solve();

        std::vector<std::complex<double> > impedances = problem.rGetImpedances();

        ImpedancePostProcessor processor(test_frequencies, impedances);

        std::cout << "\n";
        std::cout << "R5       = " << real(impedances[3])*1e-6 << " kPa.s.L^-1" << std::endl;
        std::cout << "R20      = " << real(impedances[5])*1e-6 << " kPa.s.L^-1" << std::endl;
        std::cout << "R5 - R20 = " << processor.GetR5MinusR20()*1e-6 << " kPa.s.L^-1" << std::endl;
        std::cout << "Rrs      = " << processor.GetRrs()*1e-6 << " kPa.s.L^-1" << std::endl;
        std::cout << "\n";
        std::cout << "X5       = " << imag(impedances[5])*1e-6 << " kPa.s.L^-1" << std::endl;
        std::cout << "Fres     = " << processor.GetResonantFrequency() << " Hz " << std::endl;
        std::cout << "Ax       = " << processor.GetAx()*1e-6 << " kPa.L^-1" << std::endl;
        std::cout << "\n";
    }
};