tp@0: function [outputfile] = EDB1readac(CADfile,outputfile,planecornerstype,checkgeom)
tp@0: % EDB1readac - Reads a file of type .AC (made by e.g. Invis AC3D) and saves all the geometry data in a mat-file. 
tp@0: %
tp@0: % Input parameters:
tp@0: %  CADfile     	(optional) The input file, with or without the .AC extension.
tp@0: %					If this file is not specified, a file opening window will appear.
tp@0: %  outputfile      (optional) The output file. If not specified, it will be given the
tp@0: %                   same name as the CAD file with an '_cadgeo.mat' ending instead of the '.AC'.
tp@0: %  planecornerstype (optional) Could have the value 'zero' or 'circ'.Default: 'circ'.
tp@0: %               	Affects the matrix planecorners, see below.
tp@0: %  checkgeom		(optional) If this parameter is given the value 'check', then a few checks
tp@0: %                   of the geometry consistency will be done: a check for duplicate corners
tp@0: %                   for redundant corners and for corners that are connected to only one plane.
tp@0: %					Only warnings are given. As default, no check is done.
tp@0: %  SHOWTEXT (global)	If this global parameter has the value 3 or higher, informative text will
tp@0: %                   be printed on the screen.
tp@0: %
tp@0: % Output parameters:
tp@0: %	outputfile		The name of the outputfile, generated either automatically or specified as 
tp@0: %					an input parameter.
tp@0: %
tp@0: % Output parameters in the outputfile:
tp@0: % 	corners         Matrix [ncorners,3] with the corner coordinates
tp@0: %   cornernumbers   Vector [ncorners,1] with the corner numbers used in the CAD-file. That is
tp@0: %                   in the outputfile, the corners will effectively be renumbered from 1 to ncorners
tp@0: %                   in the order they appeared in the CAD-file. 
tp@0: %   planecorners    Matrix [planes,nmaxcornersperplane] with the corner numbers that make up each plane.
tp@0: %                   Since planes can have different numbers of corners, the number of columns is the 
tp@0: %                   maximum number of corners that any plane has. NB! The corner numbers are in the
tp@0: %                   renumbered system, not in the CAD-file system.
tp@0: %   planenames      Matrix [nplanes,nmaxcharacters1] (sparse), with the planenames in the CAD file.
tp@0: %   planeabstypes   Matrix [nplanes,nmaxcharacters2] (sparse), with the absorber names in the CAD file. 
tp@0: % 	planenumbers    Vector [nplanes,1] with the plane numbers used in the CAD file.
tp@0: %   cornercrossref  Vector [nmaxnumberinCADfile,1] (sparse), which gives the matlab-file corner
tp@0: %                   numbers for each CAD-file corner number. For instance, if a CAD-file corner
tp@0: %                   number 1200 is translated to corner number 72 in the matlab-file, then
tp@0: %                   cornercrossref(1000) = 72.
tp@0: %   Snumbers        Vector [nsources,1] of the source numbers in the CAD-file. NB! Only CAD-files
tp@0: %                   v6 use source numbers; later versions use source names instead of numbers. 
tp@0: %                   For CAD-file v7 and v8, Snumbers is empty.
tp@0: %   Snames          Matrix [nsources,2] of the source names in the CAD-file. NB! Only CAD-files v7
tp@0: %                   or later use source names (such as 'A0' etc). Older versions use source numbers, 
tp@0: %                   in which case Snames is empty.
tp@0: %   Sdirectivitynames   Matrix [nsources,nmaxcharacters3] (sparse) with the source directivity names 
tp@0: %                   in the CAD file. 
tp@0: %   Sdelays         Vector [nsources,1] of the source delays in the CAD-file, in milliseconds.
tp@0: % 	Scoords         Matrix [nsources,3] of the source coordinates.
tp@0: %   Sdirections     Matrix [nsources,3] of the aim point coordinates for the sources.
tp@0: %   Srotations      Vector [nsources,1] of the source rotations, in degrees.
tp@0: %   Slevels         Matrix [nsources,6/8] of the source levels for six (CAD-files v6) or eight
tp@0: %                   octave bands (CAD-files v7 or later).
tp@0: %   Rnumbers        Vector [nreceivers,1] of the receiver numbers in the CAD-file.
tp@0: %   Rcoords         Matrix [nreceivers,3] of the receiver coordinates.
tp@0: %   CATTversionnumber   6,7 or 8
tp@0: %	planeeqs        Matrix [nplanes,4] of the plane equations as derived  from the plane definitions. 
tp@0: %                   Each row has the values [A B C D] of the plane equation on the form Ax + By + Cz = D
tp@0: %   planenvecs      Matrix [nplanes,3] of the normalized normal vectors of the planes.
tp@0: %   ncornersperplanevec     Vector [nplanes,1] which gives the number of corners for each plane.
tp@0: %   planesatcorners     Matrix [ncorners,4] which gives the plane numbers that each corner is connected to.
tp@0: %                   NB! Here it is assumed that each corner is connected to maximum four planes.
tp@0: %   nplanespercorners   Vector [ncorners,1] which gives the number of planes that each corner is connected to.
tp@0: % 	minvals         Matrix [nplanes,3] which for each plane gives the smallest coordinate in the x-, y- and z-direction.
tp@0: % 	maxvals         Matrix [nplanes,3] which for each plane gives the largest coordinate in the x-, y- and z-direction.
tp@0: %   planehasindents Vector [nplanes,1] which for each plane gives the
tp@0: %                   number of indeting corners
tp@0: %   indentingcorners  Matrix [nplanes,max(ncornersperplanevec)] which for each plane gives the number to the first corner
tp@0: %                   in a corner triplet which identifies an indenting corner. The number of the corner is the
tp@0: %                   order given in planecorners for that plane.
tp@0: %
tp@0: % Uses the functions EDB1fistr, EDB1extrnums, EDB1cross
tp@0: % 
tp@0: % ----------------------------------------------------------------------------------------------
tp@0: %   This file is part of the Edge Diffraction Toolbox by Peter Svensson.                       
tp@0: %                                                                                              
tp@0: %   The Edge Diffraction Toolbox is free software: you can redistribute it and/or modify       
tp@0: %   it under the terms of the GNU General Public License as published by the Free Software     
tp@0: %   Foundation, either version 3 of the License, or (at your option) any later version.        
tp@0: %                                                                                              
tp@0: %   The Edge Diffraction Toolbox is distributed in the hope that it will be useful,       
tp@0: %   but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS  
tp@0: %   FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.             
tp@0: %                                                                                              
tp@0: %   You should have received a copy of the GNU General Public License along with the           
tp@0: %   Edge Diffraction Toolbox. If not, see <http://www.gnu.org/licenses/>.                 
tp@0: % ----------------------------------------------------------------------------------------------
tp@0: % Peter Svensson (svensson@iet.ntnu.no) 18Jul09
tp@0: %
tp@0: % ----------------------------------------------------------------------------------------------
tp@0: %   This File is created based on EDB1eadcad.m to support the Edge Diffraction Toolbox
tp@0: % ----------------------------------------------------------------------------------------------
tp@0: % Alexander Pohl (alexander.pohl@hcu-hamburg.de) 25May2011
tp@0: %
tp@0: % [outputfile] = EDB1readac(CADfile,outputfile,planecornerstype,checkgeom);
tp@0: 
tp@0: % Modified by Peter Svensson on 3Dec2011
tp@0: 
tp@0: global SHOWTEXT
tp@0: 
tp@0: if nargin == 0
tp@0: 	CADfile = '';
tp@0: 	outputfile = '';
tp@0: 	planecornerstype = '';
tp@0: 	checkgeom = '';
tp@0: elseif nargin == 1
tp@0: 	outputfile = '';
tp@0: 	planecornerstype = '';
tp@0: 	checkgeom = '';
tp@0: elseif nargin == 2
tp@0: 	planecornerstype = '';
tp@0: 	checkgeom = '';
tp@0: elseif nargin == 3
tp@0: 	checkgeom = '';
tp@0: end
tp@0: 
tp@0: geomacc = 1e-10;
tp@0: 
tp@0: %---------------------------------------------------------------
tp@0: % If no AC-file was specified, present a file opening window
tp@0: 
tp@0: if isempty(CADfile)
tp@0: 	[CADfile,filepath] = uigetfile('*.ac','Please select the AC3Dfile');
tp@0:     
tp@0:     [filepath,temp1] = fileparts(filepath);
tp@0: 	
tp@0: 	if ~isstr(CADfile)
tp@0: 		return
tp@0: 	end
tp@0: 
tp@0:     [temp1,filestem] = fileparts(CADfile);
tp@0: 
tp@0: 
tp@0:     CADfile = [[filepath,filesep],filestem,'.ac'];
tp@0: else
tp@0: 
tp@0:     [filepath,filestem,fileext] = fileparts(CADfile);
tp@0:     CADfile = [[filepath,filesep],filestem,fileext];
tp@0: end
tp@0: 
tp@0: 
tp@0: if exist(CADfile) ~= 2
tp@0: 	error(['ERROR: AC3D-file: ',CADfile,' can not be found'])
tp@0: end
tp@0: 
tp@0: %---------------------------------------------------------------
tp@0: % If no output file was specified, construct an automatic file name
tp@0: 
tp@0: if isempty(outputfile)
tp@0: 	outputfile = [[filepath,filesep],filestem,'_cadgeo.mat'];
tp@0: end
tp@0: 
tp@0: %---------------------------------------------------------------
tp@0: % Read in the entire file into the string B
tp@0: % Find the line starts and ends
tp@0: 
tp@0: fid = fopen([CADfile],'r');
tp@0: if fid == 0
tp@0: 	error(['ERROR: The AC3Dfile ',CADfile,' could not be opened.'])
tp@0: end
tp@0: 
tp@0: %---------------------------------------------------------------
tp@0: % read header "AC3Db"
tp@0: tline = fgetl(fid);
tp@0: if (~strcmp(tline,'AC3Db'))
tp@0:     disp('AC3Db header missing');
tp@0: end
tp@0: 
tp@0: % read materials
tp@0: p_NumberOfMaterials=0;
tp@0: while(true)
tp@0:     tline = fgetl(fid);
tp@0:     tokens = readLineToToken(tline);
tp@0:     
tp@0:     if(strcmp(tokens{1},'MATERIAL'))
tp@0:         p_MaterialList{p_NumberOfMaterials+1}=tokens{2};
tp@0:         p_NumberOfMaterials=p_NumberOfMaterials+1;
tp@0:         
tp@0:     else
tp@0:         break;
tp@0:     end  
tp@0: end
tp@0: 
tp@0: % read world object (all objects are kids of this)
tp@0: m_LocalMatrix = [1.0 0.0 0.0 0.0; 
tp@0:                  0.0 1.0 0.0 0.0; 
tp@0:                  0.0 0.0 1.0 0.0; 
tp@0:                  0.0 0.0 0.0 1.0];
tp@0: 
tp@0: % Read OBJECT world
tp@0: tokens = readLineToToken(tline);
tp@0: if(~strcmp(tokens{1},'OBJECT'))    
tp@0:     disp('Section does not start with "object"');
tp@0: end
tp@0: if(~strcmp(tokens{2},'world'))    
tp@0:     disp('Section does not start with "world"');
tp@0: end
tp@0: 
tp@0: % Read kids
tp@0: tline = fgetl(fid);
tp@0: tokens = readLineToToken(tline);
tp@0: if(~strcmp(tokens{1},'kids'))    
tp@0:     disp('root object section does not start with kids');
tp@0: end
tp@0: p_Kids = str2num(tokens{2});
tp@0: 
tp@0: for i = 1:p_Kids
tp@0:     p_Objects{i}=read_ac3d_object(fid, m_LocalMatrix);
tp@0: end
tp@0: 
tp@0: fclose(fid);
tp@0: 
tp@0:  %% Export
tp@0: eval(['save a.mat p_Objects']);
tp@0: % Combine Vertices of EVERY Object to on corner list (incl. Counting PlanesPerPoly)
tp@0: for i = 1:p_Kids
tp@0:     if(i==1)
tp@0:         corners=p_Objects{i}{1}';
tp@0:         p_CornersPerPoly = size(p_Objects{i}{1},2);
tp@0:         p_PlanesPerPoly = size(p_Objects{i}{2},2);
tp@0:     else
tp@0:         corners = [corners; p_Objects{i}{1}'];
tp@0:         p_PlanesPerPoly = [p_PlanesPerPoly size(p_Objects{i}{2},2)];
tp@0:         p_CornersPerPoly = [p_CornersPerPoly size(p_Objects{i}{1},2)];
tp@0:     end
tp@0: end
tp@0: 
tp@0: % Enumerate each corner from 1:end
tp@0: ncorners = size(corners,1);
tp@0: cornernumbers=(1:ncorners)';
tp@0: 
tp@0: %Number of Planes is sum over all PlanesPerPoly
tp@0: nplanes=sum(p_PlanesPerPoly);
tp@0: 
tp@0: ncornersperplanevec=zeros(nplanes,1);
tp@0: currentPlaneIndex=1;
tp@0: for i = 1:p_Kids
tp@0:     for j=1:p_PlanesPerPoly(i);
tp@0:         ncornersperplanevec(currentPlaneIndex)=length(p_Objects{i}{2}{j});
tp@0:         currentPlaneIndex=currentPlaneIndex+1;
tp@0:     end
tp@0: end
tp@0: 
tp@0: % We need to compute the maximum number of vertices for any plane 
tp@0: % (others are filled up with 0)
tp@0: p_MaximumNumberOfVerticesOfPlane = max(ncornersperplanevec);
tp@0: 
tp@0: % Compute indices of corners for each plane
tp@0: % Simple for first Object, as here numbers are identical
tp@0: % For every further object, the number of previous used corners has be be
tp@0: % taken as offset
tp@0: planecorners=zeros(nplanes,p_MaximumNumberOfVerticesOfPlane);
tp@0: currentPlaneIndex=1;
tp@0: p_Offset=0;
tp@0: for i = 1:p_Kids
tp@0:     for j=1:p_PlanesPerPoly(i);
tp@0:         for k=1:length(p_Objects{i}{2}{j})
tp@0:             planecorners(currentPlaneIndex,k)=p_Objects{i}{2}{j}(k)+p_Offset;
tp@0:         end
tp@0:         currentPlaneIndex=currentPlaneIndex+1;
tp@0:     end
tp@0:     p_Offset=p_Offset+p_CornersPerPoly(i);
tp@0: end
tp@0: 
tp@0: planenames=zeros(nplanes,30);
tp@0: currentPlaneIndex=1;
tp@0: for i = 1:p_Kids
tp@0:     for j=1:p_PlanesPerPoly(i);
tp@0:         p_Text = sparse(double(p_MaterialList{p_Objects{i}{4}{j}+1}));
tp@0:         
tp@0:         if length(p_Text) > size(planenames,2)
tp@0:            planenames = [planenames zeros(planenumbers(currentPlaneIndex),length(p_Text)-size(planenames,2))]; 
tp@0:         end
tp@0:         planenames(currentPlaneIndex,1:length(p_Text)) = p_Text;
tp@0:         currentPlaneIndex=currentPlaneIndex+1;
tp@0:     end
tp@0: end
tp@0: 
tp@0: planeabstypes = sparse(planenames(:,1:6));
tp@0: 
tp@0: % Addition on 2.12.2011 by Peter Svensson
tp@0: % The absorber names had an initial " which must be removed.
tp@0: % We remove columns until A-Z,a-z
tp@0: 
tp@0: if size(planeabstypes,1) > 1
tp@0:     columnhasalphabet = 1 - sign( prod(sign(64-full(planeabstypes))+1));
tp@0: else
tp@0:     columnhasalphabet = 1 - sign( (sign(64-full(planeabstypes))+1));    
tp@0: end
tp@0: firstOKcolumn = find(columnhasalphabet);
tp@0: if isempty(firstOKcolumn)
tp@0:     error('ERROR: All plane absorber types have names without any letters. They must start with a letter')
tp@0: end
tp@0: firstOKcolumn = firstOKcolumn(1);
tp@0: planeabstypes = planeabstypes(:,firstOKcolumn:end);
tp@0: 
tp@0: 
tp@0: planenumbers=1:nplanes;
tp@0: 
tp@0: eval(['save ',outputfile,' corners planecorners ncornersperplanevec'])
tp@0: 
tp@0: %% Begining from here, it is code based of EDB1readcad.m without changes - except SOURCE - RECIEVER Values, as not in AC3Dfile
tp@0: %---------------------------------------------------------------
tp@0: % Change the numbering in the CAD file to a contiguous one
tp@0: 
tp@0: % Corner numbers
tp@0: % First we make an inelegant crossreference vector giving the
tp@0: % resulting corner number in the position given by the CAD-file
tp@0: % corner number.
tp@0: 
tp@0: % PS 20120414: CHECK Corners with number zero are not handled correctly.
tp@0: % Insert code from EDB1readcad here.
tp@0: 
tp@0: maxnumb = max(cornernumbers);
tp@0: cornercrossref = sparse(zeros(1,maxnumb));
tp@0: cornercrossref(cornernumbers) = [1:ncorners];
tp@0: 
tp@0: [nplanes,ncornersperplane] = size(planecorners);
tp@0: iv = find(planecorners~=0);
tp@0: planecorners(iv) = full(cornercrossref(planecorners(iv)));
tp@0: 
tp@0: %---------------------------------------------------------------
tp@0: % Go through all planes. If there is a plane definition including
tp@0: % zeros, and planecornerstype == 'circ', expand it repeating the
tp@0: % same corner order again.
tp@0: 
tp@0: if isempty(planecornerstype)
tp@0: 	planecornerstype = 'circ';
tp@0: else
tp@0: 	planecornerstype = setstr(lower(planecornerstype(1)));
tp@0: 	if planecornerstype(1) == 'z'
tp@0: 		planecornerstype = 'zero';
tp@0: 	else
tp@0: 		planecornerstype = 'circ';
tp@0: 	end
tp@0: end
tp@0: 
tp@0: if planecornerstype == 'circ'
tp@0: 	for ii = 1:nplanes
tp@0: 		iv = find( planecorners(ii,:) ~= 0);
tp@0: 		ncornersatplane = length(iv);
tp@0: 		if ncornersatplane ~= ncornersperplane
tp@0: 			pattern = planecorners(ii,iv);
tp@0: 			nrepeatings = ceil(ncornersperplane/ncornersatplane);
tp@0: 			for jj = 1:nrepeatings-1
tp@0: 				pattern = [pattern planecorners(ii,iv)];
tp@0: 			end
tp@0: 			planecorners(ii,:) = pattern(1:ncornersperplane);
tp@0: 		end
tp@0: 	end
tp@0: end
tp@0: 
tp@0: 
tp@0: %---------------------------------------------------------------
tp@0: % Find the normal vectors to the planes using the cross products
tp@0: %
tp@0: % The normal vector is basically the cross product between two vectors
tp@0: % connecting three consecutive points. If there are indents though
tp@0: % they will cause reversed normal vectors, so one must go through all
tp@0: % combinations and choose the majority normal vector.
tp@0: %
tp@0: % 26mar09  Use the fact described above for detecting indention corners
tp@0: 
tp@0: planenvecs = zeros(nplanes,3);
tp@0: planehasindents = zeros(nplanes,1);
tp@0: indentingcorners = sparse(zeros(nplanes,max(ncornersperplanevec)));
tp@0: 
tp@0: for ii = 1:nplanes
tp@0:     
tp@0: 	iv = find(planecorners(ii,:)~=0);
tp@0: 	cornerlist = planecorners(ii,iv);
tp@0: 	iv = find(cornerlist == cornerlist(1));
tp@0: 	if length(iv) > 1
tp@0: 		cornerlist = cornerlist(1:iv(2)-1);
tp@0: 	end
tp@0: 	ncorners = length( cornerlist );
tp@0: 	cornerlist = [cornerlist cornerlist(1) cornerlist(2)];
tp@0: 
tp@0: 	nvectorlist = zeros(ncorners,3);
tp@0: 	nveclen = zeros(ncorners,1);	
tp@0: 	
tp@0: 	for jj = 1:ncorners
tp@0: 		co1numb = cornerlist(jj);
tp@0: 		co2numb = cornerlist(jj+1);
tp@0: 		co3numb = cornerlist(jj+2);
tp@0: 		vec1 = corners(co1numb,:) - corners(co2numb,:);
tp@0: 		vec2 = corners(co3numb,:) - corners(co2numb,:);
tp@0: 
tp@0: 		nvec = EDB1cross(vec1.',vec2.').';
tp@0: 		nveclen(jj) = norm(nvec);
tp@0: 		if nveclen(jj) > 0
tp@0: 			nvectorlist(jj,:) = nvec./nveclen(jj);
tp@0: 		end
tp@0: 	end
tp@0: 	
tp@0: 	iv = find(nveclen < max(nveclen)*0.001);
tp@0: 	nvectorlist(iv,:) = [];
tp@0: 	nvecref = nvectorlist(1,:);
tp@0: 
tp@0:     [n1,slask] = size(nvectorlist);
tp@0: 	nvecsigns = round(sum(     (nvectorlist.').*(nvecref(ones(n1,1),:).')      ));
tp@0:     
tp@0:     if sum(nvecsigns) == 0
tp@0:         disp(' ')        
tp@0:        error(['ERROR: Plane ',int2str(planenumbers(ii)),' (plane numbering as in the CAD file) seems to be twisted.'])        
tp@0:     end
tp@0:     
tp@0:     
tp@0:     if abs(sum(nvecsigns)) ~= n1
tp@0:        disp(['Plane ',int2str(ii),' has indents!']) 
tp@0:        nindents = (n1 - abs(sum(nvecsigns)))/2;
tp@0:        planehasindents(ii) = nindents;
tp@0:        ivindent = find(nvecsigns == -1);
tp@0:        if length(ivindent) == nindents
tp@0:             indentingcorners(ii,ivindent) = 1;
tp@0:        else
tp@0:            if length(ivindent) == ncorners - nindents
tp@0:                ivindent = find(nvecsigns == 1);
tp@0:                 indentingcorners(ii,ivindent) = 1;           
tp@0:            else
tp@0:               error(['ERROR: An unexpected problem. Please report to the developer']) 
tp@0:            end
tp@0:        end
tp@0: 
tp@0:     end
tp@0:     
tp@0:     nvecdiff = [nvectorlist(2:n1,1).*nvecsigns(2:n1).' nvectorlist(2:n1,2).*nvecsigns(2:n1).' nvectorlist(2:n1,3).*nvecsigns(2:n1).'] - nvecref(ones(n1-1,1),:);
tp@0:         
tp@0:     if n1 > 2
tp@0:         nvecdiff = sum(nvecdiff.'.^2).';    
tp@0:     else
tp@0:         nvecdiff = norm(nvecdiff);
tp@0:     end
tp@0: 
tp@0:     %APO Change
tp@0:     if any(nvecdiff>1e-3),        
tp@0:         nvecdiff
tp@0:         error(['ERROR: Normal vectors for plane ',int2str(planenumbers(ii)),' (in the CAD file, = ',int2str(ii),' in the EDB1 file), get different normal vectors for consecutive corner triplets. Check the geometry in the CAD-file'])
tp@0:     elseif any(nvecdiff>1e-8)
tp@0:         nvecdiff
tp@0:         disp(['WARNING: Normal vectors for plane ',int2str(planenumbers(ii)),' (in the CAD file, = ',int2str(ii),' in the EDB1 file), get somewhat different normal vectors for consecutive corner triplets. Check the geometry in the CAD-file'])
tp@0:     end
tp@0:     
tp@0: 	if ncorners > 5 & abs(sum(nvecsigns)) <= 1
tp@0: 		disp(['WARNING for plane number ',int2str(planenumbers(ii)),' in the CAD-file'])
tp@0: 		disp(['   with the name ',strtrim(setstr(full(planenames(ii,:))))])
tp@0: 		disp('   The normal vector can not be determined for this plane because there are')
tp@0: 		disp('   the same number of inwards and outwards corners')
tp@0: 		disp('   This is a list of the possible normal vectors:')
tp@0: 		[nv1,slask] = size(nvectorlist);
tp@0: 		for kk = 1:nv1
tp@0: 			vecstr = ['   ',int2str(kk),'. ',num2str(-nvectorlist(kk,1)),' ',num2str(-nvectorlist(kk,2)),' ',num2str(-nvectorlist(kk,3))];
tp@0: 			disp(vecstr)
tp@0: 		end
tp@0: 		disp(' ')
tp@0: 	
tp@0:       preferredsign = input('   Give the number of a correct normal vector for this plane please ');
tp@0:       switchsign = nvecsigns.'./nvecsigns(preferredsign);
tp@0: 		nvectorlist = nvectorlist.*switchsign(:,ones(1,3));	
tp@0: 	else
tp@0: 		mostcommonsign = sign(sum(nvecsigns));
tp@0: 
tp@0: 		switchsign = nvecsigns.'./mostcommonsign;
tp@0: 		nvectorlist = nvectorlist.*switchsign(:,ones(1,3));
tp@0:     end
tp@0: 
tp@0: 	planenvecs(ii,:) = mean(nvectorlist);
tp@0: end
tp@0: 
tp@0: planenvecs = -planenvecs;
tp@0: 
tp@0: %---------------------------------------------------------------
tp@0: % Plane equations, as Ax + By + Cz = D for each plane
tp@0: 
tp@0: planeeqs = zeros(nplanes,4);
tp@0: planeeqs(:,1:3) = planenvecs;
tp@0: planeeqs(:,4) = sum( (planenvecs.').*(corners(planecorners(:,1),:).')  ).';
tp@0: 
tp@0: %---------------------------------------------------------------
tp@0: % Useful data: planesatcorners, minvals and maxvals
tp@0: 
tp@0: [ncorners,slask] = size(corners);
tp@0: planesatcornerhits = zeros(ncorners,nplanes);
tp@0: 
tp@0: for ii = 1:nplanes
tp@0: 	cornerlist = planecorners(ii,1:ncornersperplanevec(ii));
tp@0: 	planesatcornerhits(cornerlist,ii) = planesatcornerhits(cornerlist,ii) + 1;
tp@0: end
tp@0: 
tp@0: maxplanespercorner = 0;
tp@0: for ii = 1:ncorners
tp@0: 	nplanes = length(find(planesatcornerhits(ii,:) ~= 0));
tp@0: 	if nplanes > maxplanespercorner
tp@0: 		maxplanespercorner = nplanes;
tp@0: 	end	
tp@0: end
tp@0: 
tp@0: planesatcorners = zeros(ncorners,maxplanespercorner);
tp@0: nplanespercorners = zeros(ncorners,1);
tp@0: for ii = 1:ncorners
tp@0: 	iv = find(planesatcornerhits(ii,:)~=0);
tp@0: 	planesatcorners(ii,1:length(iv)) = iv;
tp@0: 	nplanespercorners(ii) = length(iv);
tp@0: end
tp@0: 
tp@0: % find cubic boxes inside which the planes are placed
tp@0: 
tp@0: [nplanes,slask] = size(planeeqs);
tp@0: 
tp@0: minvals = zeros(nplanes,3);
tp@0: maxvals = zeros(nplanes,3);
tp@0: 
tp@0: for ii = 1:nplanes
tp@0: 	cornerlist = planecorners(ii,:);
tp@0: 	cornerlist = cornerlist( find(cornerlist~= 0) );
tp@0: 	cornercoord = corners(cornerlist,:);
tp@0: 	minvals(ii,:) = min(cornercoord);
tp@0: 	maxvals(ii,:) = max(cornercoord);
tp@0: end
tp@0: 
tp@0: minvals = minvals - geomacc;
tp@0: maxvals = maxvals + geomacc;
tp@0: 
tp@0: %---------------------------------------------------------------
tp@0: % Check the geometry a bit
tp@0: 
tp@0: if isempty(checkgeom)
tp@0: 	checkgeom = 0;
tp@0: else
tp@0: 	if setstr(checkgeom(1)) == 'c'
tp@0: 		checkgeom = 1;
tp@0: 	else
tp@0: 		checkgeom = 0;
tp@0: 	end
tp@0: end
tp@0: 
tp@0: if checkgeom
tp@0: 
tp@0: 	badproblem = 0;
tp@0: 	disp(' ')
tp@0: 	disp('   Checking if there are corners with identical coordinates')
tp@0: 	for ii = 1:ncorners-1
tp@0: 		othercorners = [ii+1:ncorners];
tp@0: 		iv = find((corners(othercorners,1)==corners(ii,1)) & (corners(othercorners,2)==corners(ii,2)) & (corners(othercorners,3)==corners(ii,3)));
tp@0: 		if ~isempty(iv)
tp@0: 			disp(['      WARNING: corner ',int2str(ii),' and ',int2str(othercorners(iv(1))),' have the same coordinates'])
tp@0: 			disp(['      This should be fixed in the CAD-file or in the CATT-model'])
tp@0: 			badproblem = 1;
tp@0: 		end
tp@0: 	end
tp@0: 
tp@0: 	disp('   Checking if there are unused corners')
tp@0: 	iv = find(planesatcorners(:,1)==0);
tp@0: 	if ~isempty(iv)
tp@0: 		disp('      WARNING: The following corners in the CAD-file are not used:')
tp@0: 		printvec = int2str(cornernumbers(iv(1)));
tp@0: 		for iiprint = 2:length(iv)
tp@0: 			printvec = [printvec,' ',int2str(cornernumbers(iv(iiprint)))];
tp@0: 		end
tp@0: 		disp(['            ',printvec])
tp@0: 		disp(['      This is not a big problem'])
tp@0: 		planesatcorners(iv,:) = [];
tp@0: 	end
tp@0: 
tp@0: 	disp('   Checking if any corners seem redundant')
tp@0: 	iv = find(planesatcorners(:,2)==0);
tp@0: 	if ~isempty(iv)
tp@0: 		disp('      WARNING: The following corners in the CAD-file belong to only one plane:')
tp@0: 		printvec = int2str(cornernumbers(iv(1)));
tp@0: 		for iiprint = 2:length(iv)
tp@0: 			printvec = [printvec,' ',int2str(cornernumbers(iv(iiprint)))];
tp@0: 		end
tp@0: 		disp(['            ',printvec])
tp@0: 		disp(['      This should be fixed in the CAD-file or in the CATT-model'])
tp@0: 		badproblem = 1;
tp@0: 	end
tp@0: 
tp@0: 	if badproblem == 1
tp@0: 		disp(' ');
tp@0: 		error(['ERROR: Problems with the geometry defined in the CAD-file: ',CADfile])
tp@0: 	end
tp@0: 
tp@0: end
tp@0: 
tp@0: %---------------------------------------------------------------
tp@0: % Save the relevant variables in the output file
tp@0: 
tp@0: % NB! Sparse matrices can not get a non-double format
tp@0: 
tp@0: if ncorners < 256
tp@0:     planecorners = uint8(planecorners);
tp@0: elseif ncorners < 65536
tp@0:     planecorners = uint16(planecorners);    
tp@0: end   
tp@0: 
tp@0: if max(ncornersperplanevec) <= 255
tp@0:     ncornersperplanevec = uint8(ncornersperplanevec);
tp@0:     planehasindents = uint8(planehasindents);
tp@0: else
tp@0:     ncornersperplanevec = uint16(ncornersperplanevec);
tp@0:     planehasindents = uint16(planehasindents);
tp@0: end
tp@0: 
tp@0: Varlist = [' corners cornernumbers planecorners planenames planeabstypes '];
tp@0: Varlist = [Varlist,' planenumbers cornercrossref '];
tp@0: Varlist = [Varlist,' planeeqs planenvecs ncornersperplanevec planesatcorners nplanespercorners'];
tp@0: Varlist = [Varlist,' minvals maxvals planehasindents indentingcorners'];
tp@0: %Varlist = [Varlist,' Scoords Sdirections Srotations Slevels Rnumbers Rcoords CATTversionnumber'];
tp@0: 
tp@0: planenames = sparse(planenames+1-1);
tp@0: planeabstypes = sparse(planeabstypes+1-1);
tp@0: %Sdirectivitynames = sparse(Sdirectivitynames+1-1);
tp@0: 
tp@0: eval(['save ',outputfile,Varlist])
tp@0: 
tp@0: end
tp@0: 
tp@0: function [p_Plane, p_Norm, p_Mat] = read_ac3d_surface(fid, i_Vertices)
tp@0:     tline = fgetl(fid);
tp@0:     tokens = readLineToToken(tline);
tp@0:     if(~strcmp(tokens{1},'SURF'))    
tp@0:         disp(['surface section does not start with SURF',tokens{1}]);
tp@0:     end
tp@0:     if(~strcmp(tokens{2},'0x10'))    
tp@0:         disp(['only polygons are accepted TODO?', tokens{2}]);
tp@0:     end
tp@0:     
tp@0:     tline = fgetl(fid);
tp@0:     tokens = readLineToToken(tline);
tp@0:     if(~strcmp(tokens{1},'mat'))    
tp@0:         disp('missing expected token "mat"');
tp@0:     end
tp@0:     p_Mat=str2num(tokens{2});
tp@0:     
tp@0:     tline = fgetl(fid);
tp@0:     tokens = readLineToToken(tline);
tp@0:     
tp@0:     if(~strcmp(tokens{1},'refs'))    
tp@0:         disp('surface section missing refs');
tp@0:     end
tp@0:     
tp@0:     p_NumberOfRefs = str2num(tokens{2});
tp@0:     
tp@0:     if(tokens{2}<3)
tp@0:         disp(['too less points in surface:',tokens{2}]);
tp@0:         return;
tp@0:     end
tp@0:     
tp@0:     p_Plane=zeros(p_NumberOfRefs,1);
tp@0:     
tp@0:     for r=1:p_NumberOfRefs
tp@0:         tline = fgetl(fid);
tp@0:         tokens = readLineToToken(tline);
tp@0:         p_Plane(r)=str2num(tokens{1})+1; 
tp@0:     end
tp@0:     
tp@0:     v1 = i_Vertices(:,p_Plane(2))-i_Vertices(:,p_Plane(1));
tp@0:     v2 = i_Vertices(:,p_Plane(3))-i_Vertices(:,p_Plane(1));
tp@0:     p_Norm = cross(v1,v2);
tp@0:     p_Norm = p_Norm./norm(p_Norm);
tp@0: end
tp@0: 
tp@0: function [token] = readLineToToken(tline)
tp@0: 	[token{1}, remain] = strtok(tline);
tp@0:     i=1;
tp@0:     while(remain~=0)
tp@0:         i=i+1;
tp@0:         tline = remain;
tp@0:         [token{i}, remain] = strtok(tline);
tp@0:     end
tp@0: end
tp@0: 
tp@0: function [p_Object] = read_ac3d_object(fid, i_LocalMatrix)
tp@0: 
tp@0:     tline = fgetl(fid);
tp@0:     tokens = readLineToToken(tline);
tp@0:     if(~strcmp(tokens{1},'OBJECT'))    
tp@0:         disp('Section does not start with "object"');
tp@0:     end
tp@0:     if(strcmp(tokens{2},'poly'))    
tp@0:         p_Object = read_ac3d_poly(fid, i_LocalMatrix);
tp@0:     elseif(strcmp(tokens{2},'group'))    
tp@0:         disp('GROUP not implemented yet');
tp@0:     elseif(strcmp(tokens{2},'light'))    
tp@0:         disp('LIGHT not implemented yet');
tp@0:     else        
tp@0:         disp('Section does not start with valid key');
tp@0:     end
tp@0: end
tp@0: function [p_Polygon] = read_ac3d_poly(fid, i_LocalMatrix)
tp@0:     
tp@0:     p_LocalMatrix = [1.0 0.0 0.0 0.0; 
tp@0:                      0.0 1.0 0.0 0.0; 
tp@0:                      0.0 0.0 1.0 0.0; 
tp@0:                      0.0 0.0 0.0 1.0];
tp@0:     
tp@0:     tline = fgetl(fid);
tp@0:     
tp@0:     p_Vertices=0;
tp@0:     while(ischar(tline))
tp@0:         tokens = readLineToToken(tline);
tp@0: 
tp@0:         if(strcmp(tokens{1},'numvert'))
tp@0:             
tp@0:             p_ResultMatrix = i_LocalMatrix * p_LocalMatrix;
tp@0:             
tp@0:             p_NumVertices=str2num(tokens{2});
tp@0:             p_Vertices=zeros(3,p_NumVertices);
tp@0:             for l=1:p_NumVertices
tp@0:                tline = fgetl(fid);
tp@0:                tokens = readLineToToken(tline);
tp@0:                if(~ischar(tline))
tp@0:                    disp(['unexpected EOF while reading vertex (', num2str(l),'/',num2str(p_NumVertices)]);
tp@0:                end
tp@0:                if(size(tokens,2)~=3)
tp@0:                    disp(['error while reading vertex coords (', num2str(l),'/',num2str(p_NumVertices)]);
tp@0:                end
tp@0:                
tp@0:                p_Vec = p_ResultMatrix * [str2num(tokens{1}); str2num(tokens{2}); str2num(tokens{3}); 1];
tp@0:                p_Vertices(:,l) = p_Vec(1:3);
tp@0:             end
tp@0:             
tp@0:         elseif(strcmp(tokens{1},'numsurf'))
tp@0:             p_NumSurfaces=str2num(tokens{2});
tp@0:             for s=1:p_NumSurfaces
tp@0:                 [p_Planes{s}, p_Norms{s}, p_Mats{s}] = read_ac3d_surface(fid, p_Vertices);
tp@0:             end
tp@0:         elseif(strcmp(tokens{1},'loc'))
tp@0:             p_LocalMatrix(1,4)=str2num(tokens{2});
tp@0:             p_LocalMatrix(2,4)=str2num(tokens{3});
tp@0:             p_LocalMatrix(3,4)=str2num(tokens{4});
tp@0:         elseif(strcmp(tokens{1},'rot'))
tp@0:             p_LocalMatrix(1,1)=str2num(tokens{2});
tp@0:             p_LocalMatrix(1,2)=str2num(tokens{3});
tp@0:             p_LocalMatrix(1,3)=str2num(tokens{4});
tp@0:             p_LocalMatrix(2,1)=str2num(tokens{5});
tp@0:             p_LocalMatrix(2,2)=str2num(tokens{6});
tp@0:             p_LocalMatrix(2,3)=str2num(tokens{7});
tp@0:             p_LocalMatrix(3,1)=str2num(tokens{8});
tp@0:             p_LocalMatrix(3,2)=str2num(tokens{9});
tp@0:             p_LocalMatrix(3,3)=str2num(tokens{10});
tp@0:         elseif(strcmp(tokens{1},'kids'))
tp@0:             if(str2num(tokens{2})>2)
tp@0:                 disp(['Warning: Kids of Polygon not handled yet:',tokens{2}])
tp@0:                 p_Polygon=0;
tp@0:                 return;
tp@0:             end
tp@0:             p_Polygon{1}=p_Vertices;
tp@0:             p_Polygon{2}=p_Planes;
tp@0:             p_Polygon{3}=p_Norms;
tp@0:             p_Polygon{4}=p_Mats;
tp@0:             return;
tp@0:         else
tp@0:             if(strcmp(tokens{1},'name'))
tp@0:                 % Only name of poly, irrelevant -->SKIP
tp@0:             elseif(strcmp(tokens{1},'crease'))
tp@0:             else
tp@0:                 disp(['Warning: Ignored Entry in AC3D File: ', tokens{1}]);
tp@0:             end
tp@0:         end
tp@0:         tline = fgetl(fid);
tp@0:     end
tp@0: end