MATLAB Spline Fitting for Volume Approximation

The following MATLAB script generates the images and volume approximation for the light bulb example.

The MATLAB script assumes that the photo bulb.jpg is located in your current working directory.

Software required:
MATLAB version 5 or 6
Symbolic Math Toolbox

Download the M-file and picture by clicking HERE.


%volumebulb.m LAST updated 4/17/02 LFR
%MATLAB script for light bulb volume approximation--Symbolic Toolbox Required
%
%---This script produces images and cubic spline approximation for light bulb volume.
%---The script places the image in a coordinate system with horizontal
%---and vertical axes scaled by pixel information of the image.
%---Points are selected along the outline of the bulb using the mouse.
%---CHOOSE points from LEFT to RIGHT
%---Press Q or q when finished selecting points.
%
%--Press Enter when ready to begin
clear
%see help
help volumebulb
pause
%++++++++++++++++++++++++++ 
%Set up XData and YData for scaling
dpi=72;
%--Read in the image information
bulb=imread('bulb.jpg');
h=image(bulb);
xdata=get(h,'XData');
xdata=[1:xdata(2)]/dpi;
n=length(xdata);
ydata=get(h,'YData');
ydata=[1:ydata(2)]/dpi;
m=length(ydata);
close
%++++++++++++++++++++++++++
%Declare symbolic variable s--needed for integration.
syms s
% --Set axes for appropriate aspect ratio.
axis([-2 6 -3.2 3.2])
hold on
%--Place image in the graph.
h=image(bulb,'XData',[xdata(1)-1.6,xdata(n)-1.6],'YData',[ydata(m)-2.0,ydata(1)-2.0]);
%--Plot x-and y-axes.
plot([-200/dpi 600/dpi],[0 0],'k','LineWidth',2)
plot([0 0],[-450/dpi 500/dpi],'k','LineWidth',2)
%addaxes;
hold on
title('When crosshairs appear, use mouse and click to outline the bulb.')
pause(3)
title('Press Q or q when done selecting points.')
over='N';jj=1;
while over=='N'
[xx(jj),yy(jj),button]=ginput(1);
if abs(button)=='q' | abs(button)=='Q';
over='Y';
else;
plot(xx(jj),yy(jj),'ro','erasemode','none');
jj=jj+1;
end
end

%--Begin spline calculation
k=jj-1; %set back one since last time through loop we quit!
xx=xx(1:k);yy=yy(1:k);
x=[xx(1):.1:xx(k)];%NOW curve will start at first point instead of zero
y=spline(xx,yy,x);
pp=spline(xx,yy);
axis(axis);
[B,c,L,K,D]=unmkpp(pp);
numbreaks=length(B);
for i=1:numbreaks-1
xxx=[B(i):.01:B(i+1)];
yyy=polyval(c(i,1:K),xxx-xx(i));
plot(xxx,yyy,'b','LineWidth',2,'erasemode','none')
%Generate strings for polynomial title
term=['p ' num2str(i) '(x)= '];
term1=[num2str(c(i,1),3) ' (x -' num2str(xx(i),3) ')^3 + '];
term2=[num2str(c(i,2),3) ' (x -' num2str(xx(i),3) ')^2 + '];
term3=[num2str(c(i,3),3) ' (x -' num2str(xx(i),3) ') + '];
term4=[num2str(c(i,4),3)];
title([term term1 term2 term3 term4])
pause(3)
end
%Calculate volume integrals
for i=1:numbreaks-1
volumeinc(i)=pi*int((poly2sym(c(i,1:K),(s-xx(i))))^2,s,B(i),B(i+1));
end
volume=sum(volumeinc);
volume=vpa(volume,8);
pause(1)
volstring=char(volume);
titlestring=['Revolve and compute volume.'];
title(titlestring)
pause(3)
%--Generate solid of revolution.
figure(2)
titlestring=['Volume = ' volstring];
title(titlestring)
axis([-2 6 -3.2 3.2 -3.2 3.2])
n=20; 
tt=(0:n)'*2*pi/n;
xsurf=x;ysurf=y;
surfl(ones(n+1,1)*xsurf,cos(tt)*ysurf,sin(tt)*ysurf);
title(titlestring)
%
%
%MATLAB script by
% Lila F. Roberts
% Mathematics and Computer Science Department
% Georgia Southern University
% lroberts@gasou.edu
%
%AND
%
% David R. Hill
% Mathematics Department
% Temple University
% dhill001@temple.edu
%
%Produced for Demos with Positive Impact: NSF-DUE 9952306
%
%Digital photo by Judy O'Neal, North Georgia College and State University