Table Of ContentASTOKESFLOW BOUNDARYINTEGRAL MEASUREMENTOF TUBULAR STRUCTURE
CROSSSECTIONSINTWO DIMENSIONS
MarcNiethammer,Eric Pichon,AllenTannenbaum PeterJ.Mucha
fmarcn,eric,[email protected] [email protected]
Georgia InstituteofTechnology Georgia InstituteofTechnology
SchoolofElectricalandComputerEngineering SchoolofMathematics
Atlanta,GA,30332-0250,USA Atlanta,GA,30332-0160,USA
ABSTRACT image. TosolvetheStokesequationweuseaboundaryin-
tegralformulation, thus reducing thedimensionalityofthe
In this paper we will develop a method to determine
problem by one dimension. There is no need to discretize
cross sections of arbitrary two-dimensional tubular struc-
theinterioroftheobject.Werestrictourselvestothetwodi-
tures, which are allowed to branch, by means of a Stokes
mensionalproblem(flowinaplane).However,themethod-
flowbasedboundaryintegralformulation. Themeasurefor
ology developed can be extended to three dimensions (the
the cross sections for a point on the boundary of a given
boundary integral formulation for the Stokes equation is
structure will be the path obtained by integrating perpen-
standard in three dimensions). Other fluid flow equations
dicularlytotheflowlinesfromonesideoftheboundaryto
have recently been used in computer graphics and image
the other. Special emphasis will be put on the behavior at
processing,notablyforimageinpainting/reconstruction[3].
branchingpoints,thebehavioratvortices,andthenecessary
Examplesforusingtheboundaryelementapproachincom-
boundaryconditions. Themethodcanbeextendedtothree
putervisioncanbefoundin[4].
dimensionalproblems.
Section 2 introduces the Stokes equations and
presents the boundary integral formulation as given
1. INTRODUCTION by Pozrikidis [5]. Section 3 deals with the discretization
of the boundary integral equations, and introduces the
Measuring cross sections of structures in a consistent,
formalism for solving the discretized equations. Section 4
meaningful way is important for many applications. In
describes the methodology for measuring cross sections
medicalimaging,anatomicalstructuresoftentimesexhibit
once a flow field has been obtained. Examples are given
complicated, convoluted shapes. Manual thickness mea-
in Section 5. The paper ends with a conclusion and
surementsbasedonimagedata(especiallyin threedimen-
suggestionsforfutureworkinSection6.
sions) then easily become error-prone. Jones et al. [1] use
Laplace’s equation to measure cortical thickness in three-
2. BOUNDARYINTEGRALFORMULATIONOF
dimensional images whose variations can be associated
THESTOKESEQUATIONS
with many pathologies: e.g. Alzheimer’s disease. Yezzi
and Prince [2] improve the speed of Jones’ computational
Pozrikidis [5] gives a good introduction to the boundary
methodandeliminatetheneedfortheexplicitcomputation
integral method for linearized viscous flow. Higdon [6]
oftrajectories. However,asin[1],theapproachisconfined
treatstwo-dimensionalStokesflowproblemsforarbitrarily
tostructuresthatcanbedescribedbytwosimplyconnected
shapeddomains.Zebetal.[7]presentanimplementationof
boundaries. Budding structures (see Section 5) pose se-
Higdon’sapproachandextendittohandlepressurebound-
rious problems for the aforementioned approaches. Note
ary conditions. This section aims at reviewing the basics
that, since the temperature at the two boundaries is fixed,
for the formulationofthe Stokesflowprobleminterms of
nosensiblethicknesswouldbeassignedtothecavityofthe
boundaryintegrals,following[5,7].
budding structure by a Laplace equation based method: a
solutionofLaplace’sequationwillnotexhibitvortices.
This paper will deal with structures that cannot be de- 2.1. TheStokesflowequations
scribed by two simply connected boundaries, e.g. struc-
Incompressible fluid flow problems where viscous forces
tures that branch (e.g. blood vessels), and/or exhibit bud-
dominatecanbedescribedbythedivergence-freecondition
ding/constriction of boundary part(s). Instead of defining
andtheStokesequation
thickness by means of the solution of Laplace’s equation
we use the Stokes equation (describing fluid flow at low r(cid:1)u=0; (cid:0)rP +(cid:22)r2u+(cid:26)b=0: (1)
Reynoldsnumbers),withfreeslipboundaryconditions.We
Hereuisthefluidvelocity,(cid:22)theviscosityofthefluid,P is
assume that the objectto be measured is givenas a digital
the pressure and b is the body force. r indicates the gra-
ThisworkwassupportedbyAFOSR,ARO,NSF,andHEL-MRI. dient, r(cid:1) represents the divergence operator, and r2 the
Report Documentation Page Form Approved
OMB No. 0704-0188
Public reporting burden for the collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and
maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information,
including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson Davis Highway, Suite 1204, Arlington
VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to a penalty for failing to comply with a collection of information if it
does not display a currently valid OMB control number.
1. REPORT DATE 3. DATES COVERED
2003 2. REPORT TYPE 00-00-2003 to 00-00-2003
4. TITLE AND SUBTITLE 5a. CONTRACT NUMBER
A Stokes Flow Boundary Integral Measurement of Tubular Structure
5b. GRANT NUMBER
Cross Sections in Two Dimensions
5c. PROGRAM ELEMENT NUMBER
6. AUTHOR(S) 5d. PROJECT NUMBER
5e. TASK NUMBER
5f. WORK UNIT NUMBER
7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 8. PERFORMING ORGANIZATION
Georgia Institute of Technology,School of Electrical and Computer REPORT NUMBER
Engineering,Atlanta,GA,30332
9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSOR/MONITOR’S ACRONYM(S)
11. SPONSOR/MONITOR’S REPORT
NUMBER(S)
12. DISTRIBUTION/AVAILABILITY STATEMENT
Approved for public release; distribution unlimited
13. SUPPLEMENTARY NOTES
14. ABSTRACT
15. SUBJECT TERMS
16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF 18. NUMBER 19a. NAME OF
ABSTRACT OF PAGES RESPONSIBLE PERSON
a. REPORT b. ABSTRACT c. THIS PAGE 4
unclassified unclassified unclassified
Standard Form 298 (Rev. 8-98)
Prescribed by ANSI Std Z39-18
Laplacian. In the sequel we incorporate the body force in 3. BOUNDARYDISCRETIZATION
the modifiedpressurePMOD = P (cid:0)(cid:26)b(cid:1)x. To avoidun-
necessary heavy notation P is understood to mean PMOD Following [7] we assume a piecewise linear parameteriza-
inwhatfollows. tion of the boundary. The boundary @D is divided into n
straight lines@D ;k = 1;2;:::n, wheref (x) andu (x)
k i i
areassumedconstantalongtheseelements.Withthesesim-
2.2. Boundaryintegralformulation
plifications, the governing boundary integral equations (3)
TheincompressiblesolutionoftheStokesequation become
2
(cid:0)rP +(cid:22)r u+g(cid:14)(x(cid:0)x0)=0; u1(x0)
where(cid:14)(x(cid:0)x0)isatwo-dimensionalDirac-delta-function (cid:17)(x0)(cid:18)u2(x0)(cid:19)=
icsengtievreendbayt1x0 andg isanarbitraryvector-valuedconstant, 1 n n (x ) T~11k(xk;x0) T~21k(xk;x0) u1(xk)
1 4(cid:25) Xk=1 k k (cid:18)T~12k(xk;x0) T~22k(xk;x0)(cid:19)(cid:18)u2(xk)(cid:19)
ui(x)= 4(cid:25)(cid:22)Gij(x;x0)gj: (2) (cid:0) 1 G~11(xk;x0) G~21(xk;x0) f1(xk) (4)
(2) is the solution for the Stokes flow caused by a two- (cid:22)(cid:18)G~12(xk;x0) G~22(xk;x0)(cid:19)(cid:18)f2(xk)(cid:19)
dimensional point force at x0. Gij are the Green’s func- wherethetwodimensionalStokesletintegralsaregivenby
tions,whichdescribetheinfluenceofaforceatpositionx0
on the solution at position x. Since the Stokes flow prob-
lem is linear, one can express the solution to an arbitrary G~11(xk;x0)=Z (cid:0)lnr+ x^r1x2^1 dl
forcedistributionbysuperpositionoftherespectiveGreen’s @Dk
funcTtihoensstorelusstiaosnssofcoiarteeadcwhipthoitnhtisfoflrcoew. maybewritten G~12(xk;x0)=Z x^r1x2^2 dl=G~21(xk;x0)
@Dk
(cid:27)ij(x)= 41(cid:25)Tijk(x;x0)gk; G~22(xk;x0)=Z (cid:0)lnr+ x^r2x2^2 dl; (5)
@Dk
wherethestresstensorT is
ijk
andthestresstensorintegralsby
Tijk(x;x0)=
(cid:0)(cid:14)ikpj(x;x0)+ @@Gxij(x;x0)+ @@Gxkj(x;x0): T~11k(xk;x0)nk(xk)=(cid:0)4x^ini(xk)Z@Dkx^r1x4^1 dl
k i
TheGreen’sfunctionsGij (alsocalledthetwo-dimensional T~12k(xk;x0)nk(xk)=(cid:0)4x^ini(xk)Z x^r1x4^2 dl=T~21k
Stokeslets)andthestresstensorTijk aregivenby[5] @Dk
Gij =(cid:0)(cid:14)ijlnr+ x^rix2^j; Tijk =(cid:0)4x^ix^rj4x^k: T~22k(xk;x0)nk(xk)=(cid:0)4x^ini(xk)Z@Dkx^r2x4^2 dl: (6)
Here (cid:14)ik is the Kronecker-delta-function,r = kx(cid:0)x0k2, dl is the differential along the boundary element. Equa-
andx^ =x(cid:0)x0. tions(5,6)canbeevaluatedanalytically; theydescribethe
Forapointx0 insideorontheboundarythesolutionto influence of a whole boundary segment on a point x0 lo-
theStokesequation(1)intermsofboundaryintegralequa- cated in or on the boundary of the domain D. The inte-
tionsis grands of Equations (5, 6) exhibit singularities. These are
removable (except for the case of points which lie exactly
(cid:17)(x0)uj(x0)=(cid:0) 1 fi(x)Gij(x;x0)dl(x)+ ontheendpointofaboundaryelement)2.
4(cid:25)(cid:22)Z Toevaluatethediscretizedboundaryintegrals(4)weuse
@D
1 thefollowingtwostepapproach(see[7]):
ui(x)Tijk(x;x0)nk(x)dl(x); (3)
4(cid:25) Z@D a) Evaluate the velocity boundary integrals for the n
forthevelocities,wheref (x) = (cid:27) (x)n (x)isthestress center points of the boundary elements, by specify-
i ij j
force, and n(x) is the normal vector to the boundary @D ing 2n boundary conditions. This is a linear equa-
at position x, pointing towards the interior of the solution tionsystemwith4nunknowns(f1,f2,u1,u2 forev-
domain D. (cid:17)(x0) accounts for the dependencyof the ve- ery boundaryelement) and 2nequations. Fast algo-
locityintegralonthelocationofx0: (cid:17)(x0)=1ifx0 2D, rithms for solving such boundary integral equations
(cid:17)(x0)= 12 ifx0 2@D. exist,e.g.[8].
b) Solveforthevelocityatarbitrarypositionsinsidethe
1Notethatwearelookingatthetwo-dimensionalStokesflowproblem.
This is the solution of the two-dimensional equation, written with Ein- domainDbyusingEquation(4).
stein’ssummationconvention(summingovermultiplyoccurringindices).
Solutionsinthreedimensionsalsoexist,astreatedin[5].Wedenotevector 2We omit the explicit analytical expressions for the integrals due to
componentsbysubscripts,i.e.u1andu2forthecomponentsofu. spaceconstraints.
Theboundarydecomposesintotheinlet,outletandwall
boundaryparts:@D=@D [@D [@D .Weprescribethe
i o w
velocityprofileu1,u23attheinlet(@Di)andassumestress
freeboundaryconditionsfortheoutlet(s)(i.e.f =f =0
1 2
for x 2 @D ). @D exhibits perfect slip without perme-
o w
ation: u2 = 0, f1 = 0. Note, that these boundary condi-
tionsdirectlytranslatetothethree-dimensionalcase. While
the Stokes equation itself is linear and reversible, we note
thatthegeneratedflowfieldswillvaryslightlyunderinter-
changing the identification of the inlet and outlet, because
of the prescribed boundary conditions. These differences, (a)Discretization,withboundaryelements
andassociatednormalvectors.
however,willlargelyberestrictedtotheimmediatevicini-
ties of the inlet and outlet for reasonable prescribed inlet
flows.
4. MEASUREMENTOFCROSSSECTIONS
Tomeasurecrosssectionsofastructureweproposetointe-
grateperpendiculartothecomputedflowfield.Wecompute
l
x(l)= vdl; (7)
(b)Normalizedvelocityfield.Lightgrayar-
Z0
rowsindicatetheprescribedinflowveloci-
ties.
wherev 2R2isaunitvectornormaltothefluidspeeduat
allpoints,x(L) 2 @DwithL 6= 0,andx(0)isthebound-
ary point for whichthe cross sectionis to be measured. L
isthenthepathlengthfromtheinitialtothefinalboundary
point. We setthe initial conditionofthe vectorto the nor-
mal vector for the boundary element of the starting point,
anddisallowdirectionalflipsthroughouttheintegration.
For the evaluation of (7) we use a modified version
of the variable order Runge-Kutta method by Cash and
Karp[9].Itismodifiedinthesensethataminimalstepsize
h isintroduced. Oncetheintegrationalgorithmtriesto
min
use a step size h (cid:20) h we stop the integration and de-
min
clare the current integrator state as the end point. In this (c)Crosssections.
way we stop the integration when either a boundary point
or a stagnant point (where the step size will tend to zero) Fig.1. Thebuddingdomain.
is reached. The combination of an integrator with adap-
tivestepsizeandtheboundaryintegralmethodthusyields
a method that will naturally adapt to the problem at hand.
If a lot is known about the shape of the boundary (i.e. we its boundary elements and the associated normal vectors.
haveahighresolutionimage)thisinformationwillbeused, Figures 1(b) and 2(a) show the normalized flow fields of
without requiring too many steps to integrate from an ini- theexamples,andthecomputertomography(CT)imageof
tialpointontheboundarytothetargetpoint. Forverythin the spinous process4. Figures 1(c) and 2(b) show their re-
structures(thicknessintheorderofafewpixels)themethod spectivemeasuresforthecrosssections. Ofspecialinterest
willautomaticallyresultinsubpixelaccuracy. isthebehaviorclosetothebranchingpointforthespinous
process. While the Stokes flow based approach produced
sensibleresults,amethodbasedforexampleonfindingthe
5. EXAMPLES
nearestpointontheoppositeboundarywouldleadtounrea-
This section presents examples for cross section measure- sonableresults.Thebuddingdomainexhibitsavortexinits
ments: a budding domain, and a spinous process of a ver- flowfield(seeFigure1(b)). Thepathlinesperpendicularto
tebra. We set (cid:22) = 1 and h = 0:005, and the relative theflowfieldofthebuddingdomaincapturetherectangular
min
errortolerance oftheintegratorto (cid:15) = 10(cid:0)4. Figures1(a) part of the domain well. The cavity shows the end points
showsthediscretizeddomainofthebuddingexample,with of several path lines at the stagnant center of the vortex.
3Overlinedvariablesdenotevariablesgiveninalocallydefinedcoordi- 4Theauthorswouldlike tothankProf. Yamamotoforprovidingthe
natesystemofaboundaryelement,wherethex2directioncoincideswith image. The boundaries were extracted by means of an active contours
thenormalvector. basedapproach.
6. CONCLUSIONANDFUTUREWORK
In this paper a Stokes flow based method for measuring
crosssectionsoftwo-dimensionalstructureswaspresented.
This method is more versatile than the Laplace equation
based approach. A greater variety of domains can be han-
dled; e.g.branchingandbuddingdomains. There isnore-
striction to domains that are enclosed by two simply con-
nected surfaces. A boundary integral approach was used
to deal easily with arbitrary shaped boundaries, yielding a
cleanmethod,withouttheneedfordiscretizationofthein-
terior of the domain. The method is extendable to three
dimensions, where one of the natural geometric measure-
ments would be cross-sectional surface area, and we will
expect to find lines of stagnant points (instead of isolated
stagnantpointsatthecenterofvortices).Futureworkcould
includetheextensiontothreedimensionalproblems,inves-
tigations on the handling of stagnant points/lines (the vor-
tices),andcenterlineconstruction.
(a)NormalizedvelocityfieldandCTimage.Light
grayarrowsindicatetheprescribedinflowveloci- 7. REFERENCES
ties.
[1] S. E. Jones, B. R. Buchbinder, and I. Aharon, “Three-
dimensional mapping of cortical thickness using Laplace’s
equation,” HumanBrainMapping,vol.11,pp.12–32,2000.
[2] A. Yezziand J.L. Prince, “APDE approachfor measuring
tissue thickness,” in Proceedings of the International Con-
ference on ComputerVision andPattern Recognition. IEEE,
2001,vol.1,pp.87–92.
[3] M.Bertalmio,A.L.Bertozzi,andG.Sapiro, “Navier-Stokes,
Fluid Dynamics, and Image and Video Inpainting,” in Pro-
ceedingsoftheInternationalConferenceonComputerVision
andPatternRecognition.IEEE,2001,vol.1,pp.355–362.
[4] G.G. Gu and M. A. Gennert, “Boundary element methods
forsolvingPoissonequationsincomputervisionproblems,”
inProceedingsoftheInternationalConferenceonComputer
VisionandPatternRecognition.IEEE,1991,pp.546–551.
[5] C.Pozrikidis, Boundaryintegralandsingularitymethodsfor
linearizedviscousflow, CambridgeTextsinAppliedMathe-
matics.CambridgeUniversityPress,1992.
[6] J.J.L.Higdon,“Stokesflowinarbitrarytwo-dimensionaldo-
mains: shearflowoverridgesandcavities,” JournalofFluid
(b)Crosssections.
Mechanics,vol.159,pp.195–226,1985.
Fig.2. Thespinousprocess. [7] A. Zeb, L. Elliott, D. B. Ingham, and D. Lesnic, “The
boundary element method for the solution of Stokes equa-
tions in two-dimensional domains,” Engineering Analysis
withBoundaryElements,vol.22,pp.317–326,1998.
[8] W. Ye, J. Kanapka, X. Wang, and J. White, “Efficient and
AccuracyImprovementsforFastStokes,aPrecorrected-FFT
Accelerated3-DStokesSolver,” inProceedingsofModeling
One could for example connect these path lines by requir- andSimulationofMicrosystems,1999,pp.502–505.
ing minimal change of direction at the joining point. We
[9] J. R. Cash and A. H. Karp, “A variable order Runge-Kutta
see that through the presence of the vortex, a measure for
methodforinitialvalueproblemswithrapidlyvaryingright-
thecrosssectionofthebuddingpartofthestructureaswell
handsides,”ACMtransactionsonmathematicalsoftware,vol.
asoftherectangularpartcanbegiven.Thiswouldnothave
16,no.3,pp.201–222,1990.
beenpossiblewithanapproachbasedonLaplace’sequation
(heatflow). Note,thattheobservationsmadeinthissection
willhavecorrespondencesinthethree-dimensionalcase.
