Interpolating Functions to Approximate Volumes

This notebook utilizes digital photos and the method of cross sections to approximate the volume of a sea worm.  A digital photograph is read and the "edges" of the worm are approximated using an interpolating function (by default a third order interpolation function).  The cross sections are constructed and displayed in an animation.  The volume is approximated using a sum of the volumes of cross-sectional disks and also by using an integral involving the interpolating functions.

Read in graphics packages and turn off spelling checks (turning off spelling checks is merely a convenience).

In[1]:=

<<Graphics`Graphics3D` <<Graphics`Shapes` <<Graphics`FilledPlot` <<Graphics`Colors` <<Graphics`ParametricPlot3D` Off[General :: spell1] Off[General :: spell]

Import the digital photo.  This demo does not account for forshortening so it is useful to include in a photo some article such as a ruler that would give a sense of scale.  In this example, the photo is in JPEG format with a resolution of 72 pixels/inch.  This will be used in the calculations.  NOTE:  If you use your own photograph or if you download the image, you will need to change the path before executing.

In[8]:=

worm = Import["C:\Documents and Settings\user\Desktop\Volumeapps\seaworm_lr.jpg"]

Out[8]=

⁃Graphics⁃

In[9]:=

Show[worm]

[Graphics:HTMLFiles/seaworm_mma_5.gif]

Out[9]=

⁃Graphics⁃

Select points along the top edge of the worm in the photograph.  Click on the picture and while holding the Control Key, click the edge while moving from left to right (if you look carefully you will see white pixels corresponding to the selected points).  After selecting the points, choose Edit->Copy to copy the coordinates on the clipboard.  Then paste into the notebook input cell for the variable "toppoints".

In[10]:=

RowBox[{RowBox[{toppoints, =, RowBox[{{, RowBox[{RowBox[{{, RowBox[{126.997, ,, 192}], }}], ,, ... Box[{{, RowBox[{563.001, ,, 145}], }}], ,, RowBox[{{, RowBox[{575.001, ,, 136}], }}]}], }}]}], ;}]

Compute interpolating function and find maximum and minimum domain and range values.  These values will be used later when constructing graphs.

In[11]:=

top = Interpolation[toppoints] ;

In[12]:=

firstcoordstop = Table[toppoints[[i, 1]], {i, 1, Length[toppoints]}] ; secondcoordstop = Table ... topminx = Min[firstcoordstop] ; topmaxy = Max[secondcoordstop] ; topminy = Min[secondcoordstop] ;

Select points along the bottom edge of the photograph using the same approach as for the top points.  Copy and paste into the input for variable "bottompoints."

In[18]:=

RowBox[{RowBox[{bottompoints, =, RowBox[{{, RowBox[{RowBox[{{, RowBox[{117.997, ,, 180}], }}], ... Box[{{, RowBox[{545.001, ,, 140}], }}], ,, RowBox[{{, RowBox[{559.001, ,, 130}], }}]}], }}]}], ;}]

Compute interpolating function and find maximum and minimum domain and range values.  These values will be used later when constructing graphs.

In[19]:=

bottom = Interpolation[bottompoints] ;

In[20]:=

firstcoordsbottom = Table[bottompoints[[i, 1]], {i, 1, Length[bottompoints]}] ; secondcoordsbo ... irstcoordsbottom] ; bottommaxy = Max[secondcoordsbottom] ; bottomminy = Min[secondcoordsbottom] ;

Compute upper and lower bounds for the graphs.

In[26]:=

upperx = Min[bottommaxx, topmaxx] ; lowerx = Max[bottomminx, topminx] ;

Set PlotRange for the 3D plots.  The parameter "addabit" gives a little space around the actual graph.  It can be adjusted to suit the user.

In[28]:=

addabit = 10 ;

In[29]:=

plotrange3d = {{lowerx - addabit, upperx + addabit}, {-Min[bottommaxy, topmaxy] - addabit, Max ...  topmaxy] + addabit}, {-Min[bottommaxy, topmaxy] - addabit, Max[bottommaxy, topmaxy] + addabit}} ;

Show the interpolating function that fits the points.

In[30]:=

wormgraph = Show[Plot[{top[x], bottom[x]}, {x, lowerx, upperx}, AxesAutomatic, PlotRan ...  1, 0], Thickness[.005]}, {CMYColor[0, 1, 0], Thickness[.005]}}, DisplayFunctionIdentity]]

Out[30]=

⁃Graphics⁃

In[31]:=

all = Show[{worm, wormgraph}, AspectRatioAutomatic, ImageSize500, DisplayFunction$DisplayFunction]

[Graphics:HTMLFiles/seaworm_mma_19.gif]

Out[31]=

⁃Graphics⁃

Set up the partition.  

In[32]:=

n = 21 ; Δx = (upperx - lowerx)/(n - 1) ; xgrid[i_] = lowerx + (i - 1) Δx ; partitio ... id[i]]}, {xgrid[i], 0, top[xgrid[i]]}}]}], ImageSize300, DisplayFunctionIdentity]

In[36]:=

partitionall := Table[partition[i], {i, 1, n - 1}]

Locate the centers of the circular cross sections and the radius function which is computed as 1/2 (top(x) - bottom(x)).

In[37]:=

centers[x_] = {x, (top[x] + bottom[x])/2} ; radii[x_] = (top[x] - bottom[x])/2 ;

Set up the ends of the cross sectional pieces and store as a graphics array.

In[39]:=

rightends = Table[ParametricPlot3D[{xgrid[i], u * radii[(xgrid[i - 1] + xgrid[i])/2] * Cos[t], ... ; {20, 20}, DisplayFunctionIdentity, BoxedFalse, ImageSize350], {i, 2, n}]

Out[39]=

{⁃Graphics3D⁃, ⁃Graphics3D⁃, ⁃Graphics3D⁃, ⁃Graphics ... 9;Graphics3D⁃, ⁃Graphics3D⁃, ⁃Graphics3D⁃, ⁃Graphics3D⁃}

In[40]:=

leftends = Table[ParametricPlot3D[{xgrid[i - 1], u * radii[(xgrid[i - 1] + xgrid[i])/2] * Cos[ ... ; {20, 20}, BoxedFalse, DisplayFunctionIdentity, ImageSize350], {i, 2, n}]

Out[40]=

{⁃Graphics3D⁃, ⁃Graphics3D⁃, ⁃Graphics3D⁃, ⁃Graphics ... 9;Graphics3D⁃, ⁃Graphics3D⁃, ⁃Graphics3D⁃, ⁃Graphics3D⁃}

Set up the circular outside edges of the cross sections.

In[41]:=

xsections := Table[ParametricPlot3D[{u, radii[(xgrid[i - 1] + xgrid[i])/2] * Cos[t], radii[(xg ... ; {20, 20}, BoxedFalse, DisplayFunctionIdentity, ImageSize350], {i, 2, n}]

Represent the filled region in the vertical plane using rectangles.

In[42]:=

numpts = 50 ; dx = (upperx - lowerx)/numpts ; xrect[i_] = lowerx + i * dx ; vertices[k_] := {{ ... ect[k], 0, bottom[xrect[k]]}, {xrect[k], 0, top[xrect[k]]}, {xrect[k - 1], 0, top[xrect[k - 1]]}}

In[46]:=

region = Show[Graphics3D[{SurfaceColor[Magenta], EdgeForm[], Table[Polygon[vertices[k]], {k, 1, numpts}]}]]

[Graphics:HTMLFiles/seaworm_mma_31.gif]

Out[46]=

⁃Graphics3D⁃

Show the region and the partition.

In[47]:=

Show[{region, partitionall}, ViewPoint {1, -3, .5}, BoxedFalse, ImageSize350]

[Graphics:HTMLFiles/seaworm_mma_34.gif]

Out[47]=

⁃Graphics3D⁃

Generate a 3-dimensional model of the worm.

In[48]:=

worm3d = Show[{leftends[[1]], rightends[[20]], ParametricPlot3D[{x, radii[x] * Cos[t], radii[x ... n$DisplayFunction, ViewPoint {2, -3, .5}, AxesFalse, BoxedFalse] ;

[Graphics:HTMLFiles/seaworm_mma_37.gif]

Generate the approximating cross sections as a sequence of frames that can be animated.

In[49]:=

For[m = 1, m<n, m ++, Show[{region, partitionall, Table[{xsections[[k]], rightends[[k]], le ... 3, 0}, PlotRangeplotrange3d, DisplayFunction$DisplayFunction, BoxedFalse]]

[Graphics:HTMLFiles/seaworm_mma_59.gif]

Approximate the volume using NIntegrate.  Because the interpolating functions are piecewise, NIntegrate may not always give convergence of the integral.  This can be avoided by limiting the number of points at which Mathematica computes function values.  In this example, the option MaxPoints is set to 100000000.  The integral is scaled by 1/72^3 because the resolution of the photograph is 72 pixels per inch.  

In[50]:=

volumeintegral = 1/72^3 * (NIntegrate[Pi * (radii[x]/72)^2, {x, lowerx, upperx}, MaxPoints100000000])

Out[50]=

0.0000545774

Convert to ml. 1 sqare inch is approximately 16.387 ml.

In[51]:=

RowBox[{volumeml, =, RowBox[{volumeintegral, *, 16.387}]}]

Out[51]=

0.00089436

Approximate the volume using a sum of the volumes of the approximating cross sections and convert to ml.

In[52]:=

volumesum = Sum[Pi * (radii[(xgrid[i - 1] + xgrid[i])/2]/72)^2 * Δx/72, {i, 2, n}]

Out[52]=

0.283187

In[53]:=

RowBox[{volumesumml, =, RowBox[{volumesum, *, 16.387}]}]

Out[53]=

4.64059

The codes below are not necessary for the demonstration but might be helpful if you wish to generate demos using a different photograph.  The two code sequences illustrate the selecton of the points along top and bottom and the generation of the interpolation functions.

In[54]:=

For[k = 1, k<Length[toppoints], k ++, If[k<Length[toppoints] - 1, Show[{worm, Table[Grap ... ; {Thickness[.005], CMYColor[0, 1, 0]}, DisplayFunctionIdentity]}, ImageSize350]]]

[Graphics:HTMLFiles/seaworm_mma_109.gif]

In[55]:=

For[k = 1, k<Length[bottompoints], k ++, If[k<Length[bottompoints] - 1, Show[{worm, Tabl ... ; {Thickness[.005], CMYColor[0, 1, 0]}, DisplayFunctionIdentity]}, ImageSize350]]]

[Graphics:HTMLFiles/seaworm_mma_142.gif]

This Mathematica notebook was developed for Demos with Positive Impact (NSF DUE-9952306) by
Lila F. Roberts  (lila.roberts@gcsu.edu)
Georgia College and State University

Special thanks to James P. Braselton (jbraselton@georgiasouthern.edu) for Mathematica advice and assistance.


Created by Mathematica  (December 31, 2003)