00001 function [err,PatchHandles] = ShowPlane (EqPlane,Opt)
00002 %
00003 % [err,PatchHandles] = ShowPlane (EqPlane,Opt)
00004 %
00005 %Purpose:
00006 % This function displays a plane with equation ax + by + cz +d = 0
00007 % The plane is defined by one sqaure patch
00008 % The planes span the current axis limits of the figure.
00009 %
00010 %Input Parameters:
00011 % EqPlane is a Nx4 matrix containing the equation of the plane containing each
00012 % triplet in Triplets. The plane passing by triplet i is speicifed in
00013 % EqPlane(i,:) the plane would be
00014 % EqPlane(i,1)x + EqPlane(i,2)y + EqPlane(i,3)z + EqPlane(i,4) = 0
00015 % Opt is the options structure
00016 % .Fig is a handle to the figure you want the planes displayed in,
00017 % default is the current figure.
00018 % .WriteIV if this string is not empty, the planes that are displayed
00019 % on the graph are written to an inventor format file
00020 % .units 'mm' or 'tesscon'. default is tesscon
00021 % if you specify mm, then the xyz coordinates are transformed to tesscon
00022 % before writing them out. (*319.7). The idea is to write all inventor
00023 % files in tesscon units. This option will only be used if WriteIV is not empty.
00024 % .OvrWrite (0/1) default is 0, flag for overwriting existing .iv file
00025 %
00026 %Output Parameters:
00027 % err : 0 No Problem
00028 % : 1 Mucho Problems
00029 %
00030 % PatchHandles : the handle to the patches displayed on the figure
00031 %
00032 %More Info :
00033 % see Plane_Equation
00034 %
00035 %
00036 %
00037 % Author : Ziad Saad
00038 % Date : Thu Oct 22 20:19:36 CDT 1998
00039
00040
00041 %Define the function name for easy referencing
00042 FuncName = 'ShowPlane';
00043
00044 %initailize return variables
00045 err = 1;
00046
00047
00048 %check on the size of input data
00049 if (nargin == 1),
00050 Opt.Fig = [];
00051 Opt.OvrWrite = 0;
00052 Opt.WriteIV = '';
00053 end
00054
00055 if (~isfield(Opt,'OvrWrite') | isempty(Opt.OvrWrite)), Opt.OvrWrite = 0; end
00056
00057 if (size(EqPlane,2) ~= 4), err = ErrEval(FuncName,'Err_Bad size for EqPlane'); return; end
00058
00059 Nplanes = size(EqPlane,1);
00060 Nnodes = 4.* Nplanes;
00061
00062 %pop up a figure
00063 if (~isfield(Opt,'Fig') | isempty(Opt.Fig)),
00064 Opt.Fig = gcf;
00065 end
00066
00067 figure(Opt.Fig);
00068
00069 %get the axis roperties of the figure
00070 Xlim = get(gca,'Xlim');
00071 Ylim = get(gca,'ylim');
00072 Zlim = get(gca,'Zlim');
00073
00074 inode = 0;
00075 Node = zeros(Nnodes,3);
00076 Pat = zeros(Nplanes,4);
00077 ztmp = zeros(1,4);
00078
00079 for (i=1:1:Nplanes),
00080 %using the XY limits, find the corrspondign z values
00081 if (EqPlane(i,3) ~= 0),
00082 ztmp(1) = (-EqPlane(i,4) - EqPlane(i,1).*Xlim(1) - EqPlane(i,2).*Ylim(1)) ./ EqPlane(i,3);
00083 ztmp(2) = (-EqPlane(i,4) - EqPlane(i,1).*Xlim(2) - EqPlane(i,2).*Ylim(1)) ./ EqPlane(i,3);
00084 ztmp(3) = (-EqPlane(i,4) - EqPlane(i,1).*Xlim(2) - EqPlane(i,2).*Ylim(2)) ./ EqPlane(i,3);
00085 ztmp(4) = (-EqPlane(i,4) - EqPlane(i,1).*Xlim(1) - EqPlane(i,2).*Ylim(2)) ./ EqPlane(i,3);
00086 %form the four points on the plane
00087 Node(inode+1,:) = [Xlim(1) Ylim(1) ztmp(1)];
00088 Node(inode+2,:) = [Xlim(2) Ylim(1) ztmp(2)];
00089 Node(inode+3,:) = [Xlim(2) Ylim(2) ztmp(3)];
00090 Node(inode+4,:) = [Xlim(1) Ylim(2) ztmp(4)];
00091
00092 elseif (EqPlane(i,2) ~= 0),
00093 ytmp(1) = (-EqPlane(i,4) - EqPlane(i,1).*Xlim(1) - EqPlane(i,3).*Zlim(1)) ./ EqPlane(i,2);
00094 ytmp(2) = (-EqPlane(i,4) - EqPlane(i,1).*Xlim(2) - EqPlane(i,3).*Zlim(1)) ./ EqPlane(i,2);
00095 ytmp(3) = (-EqPlane(i,4) - EqPlane(i,1).*Xlim(2) - EqPlane(i,3).*Zlim(2)) ./ EqPlane(i,2);
00096 ytmp(4) = (-EqPlane(i,4) - EqPlane(i,1).*Xlim(1) - EqPlane(i,3).*Zlim(2)) ./ EqPlane(i,2);
00097 %form the four points on the plane
00098 Node(inode+1,:) = [Xlim(1) ytmp(1) Zlim(1)];
00099 Node(inode+2,:) = [Xlim(2) ytmp(2) Zlim(1)];
00100 Node(inode+3,:) = [Xlim(2) ytmp(3) Zlim(2)];
00101 Node(inode+4,:) = [Xlim(1) ytmp(4) Zlim(2)];
00102 elseif (EqPlane(i,1) ~= 0),
00103 xtmp(1) = (-EqPlane(i,4) - EqPlane(i,2).*Ylim(1) - EqPlane(i,3).*Zlim(1)) ./ EqPlane(i,1);
00104 xtmp(2) = (-EqPlane(i,4) - EqPlane(i,2).*Ylim(2) - EqPlane(i,3).*Zlim(1)) ./ EqPlane(i,1);
00105 xtmp(3) = (-EqPlane(i,4) - EqPlane(i,2).*Ylim(2) - EqPlane(i,3).*Zlim(2)) ./ EqPlane(i,1);
00106 xtmp(4) = (-EqPlane(i,4) - EqPlane(i,2).*Ylim(1) - EqPlane(i,3).*Zlim(2)) ./ EqPlane(i,1);
00107 %form the four points on the plane
00108 Node(inode+1,:) = [xtmp(1) Ylim(1) Zlim(1)];
00109 Node(inode+2,:) = [xtmp(2) Ylim(2) Zlim(1)];
00110 Node(inode+3,:) = [xtmp(3) Ylim(2) Zlim(2)];
00111 Node(inode+4,:) = [xtmp(4) Ylim(1) Zlim(2)];
00112 end
00113
00114 %verify that all points are on plane for debugging only
00115 %sum(EqPlane.*[Node(inode+1,:) 1])
00116 %sum(EqPlane.*[Node(inode+2,:) 1])
00117 %sum(EqPlane.*[Node(inode+3,:) 1])
00118 %sum(EqPlane.*[Node(inode+4,:) 1])
00119
00120 %form the faceset connection, for the patch
00121 Pat(i,:) = [inode+1 inode+2 inode+3 inode+4];
00122
00123 inode = inode + 4;
00124
00125 end
00126
00127 %Now display those patches
00128 vc = 0:1./Nnodes:(1-1./Nnodes);
00129 tcolor = [vc' (1-vc)' vc'];
00130
00131 PatchHandles = patch('vertices',Node,'faces',Pat,...
00132 'FaceVertexCData',tcolor,'FaceColor','flat');
00133
00134 %write to file ?
00135 if (isfield(Opt,'WriteIV') & ~isempty(Opt.WriteIV)),
00136 fprintf (1,'Saving patches to iv file %s ...\n',Opt.WriteIV);
00137 Opt.OptIV.BaseCol = tcolor;
00138 if (~isfield(Opt,'units') | isempty(Opt.units)), Opt.units = 'tesscon'; end
00139 Opt.OptIV.units = Opt.units;
00140 Opt.OptIV.OvrWrite = Opt.OvrWrite;
00141 Opt.OptIV.verbose = 0;
00142 [err] = WriteInv21Surf(Opt.WriteIV,Node,Pat,Opt.OptIV); %exclude redundant last node
00143 end
00144
00145 err = 0;
00146 return;
00147
