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