Mercurial > hg > human-echolocation

annotate private/EDB1readac.m @ 18:2d5f50205527 jabuilder_int tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Escape the trailing backslash as well
author Chris Cannam
date Tue, 30 Sep 2014 16:23:00 +0100
parents 90220f7249fc
children
rev   line source
tp@0 1 function [outputfile] = EDB1readac(CADfile,outputfile,planecornerstype,checkgeom)
tp@0 2 % 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 3 %
tp@0 4 % Input parameters:
tp@0 5 % CADfile (optional) The input file, with or without the .AC extension.
tp@0 6 % If this file is not specified, a file opening window will appear.
tp@0 7 % outputfile (optional) The output file. If not specified, it will be given the
tp@0 8 % same name as the CAD file with an '_cadgeo.mat' ending instead of the '.AC'.
tp@0 9 % planecornerstype (optional) Could have the value 'zero' or 'circ'.Default: 'circ'.
tp@0 10 % Affects the matrix planecorners, see below.
tp@0 11 % checkgeom (optional) If this parameter is given the value 'check', then a few checks
tp@0 12 % of the geometry consistency will be done: a check for duplicate corners
tp@0 13 % for redundant corners and for corners that are connected to only one plane.
tp@0 14 % Only warnings are given. As default, no check is done.
tp@0 15 % SHOWTEXT (global) If this global parameter has the value 3 or higher, informative text will
tp@0 16 % be printed on the screen.
tp@0 17 %
tp@0 18 % Output parameters:
tp@0 19 % outputfile The name of the outputfile, generated either automatically or specified as
tp@0 20 % an input parameter.
tp@0 21 %
tp@0 22 % Output parameters in the outputfile:
tp@0 23 % corners Matrix [ncorners,3] with the corner coordinates
tp@0 24 % cornernumbers Vector [ncorners,1] with the corner numbers used in the CAD-file. That is
tp@0 25 % in the outputfile, the corners will effectively be renumbered from 1 to ncorners
tp@0 26 % in the order they appeared in the CAD-file.
tp@0 27 % planecorners Matrix [planes,nmaxcornersperplane] with the corner numbers that make up each plane.
tp@0 28 % Since planes can have different numbers of corners, the number of columns is the
tp@0 29 % maximum number of corners that any plane has. NB! The corner numbers are in the
tp@0 30 % renumbered system, not in the CAD-file system.
tp@0 31 % planenames Matrix [nplanes,nmaxcharacters1] (sparse), with the planenames in the CAD file.
tp@0 32 % planeabstypes Matrix [nplanes,nmaxcharacters2] (sparse), with the absorber names in the CAD file.
tp@0 33 % planenumbers Vector [nplanes,1] with the plane numbers used in the CAD file.
tp@0 34 % cornercrossref Vector [nmaxnumberinCADfile,1] (sparse), which gives the matlab-file corner
tp@0 35 % numbers for each CAD-file corner number. For instance, if a CAD-file corner
tp@0 36 % number 1200 is translated to corner number 72 in the matlab-file, then
tp@0 37 % cornercrossref(1000) = 72.
tp@0 38 % Snumbers Vector [nsources,1] of the source numbers in the CAD-file. NB! Only CAD-files
tp@0 39 % v6 use source numbers; later versions use source names instead of numbers.
tp@0 40 % For CAD-file v7 and v8, Snumbers is empty.
tp@0 41 % Snames Matrix [nsources,2] of the source names in the CAD-file. NB! Only CAD-files v7
tp@0 42 % or later use source names (such as 'A0' etc). Older versions use source numbers,
tp@0 43 % in which case Snames is empty.
tp@0 44 % Sdirectivitynames Matrix [nsources,nmaxcharacters3] (sparse) with the source directivity names
tp@0 45 % in the CAD file.
tp@0 46 % Sdelays Vector [nsources,1] of the source delays in the CAD-file, in milliseconds.
tp@0 47 % Scoords Matrix [nsources,3] of the source coordinates.
tp@0 48 % Sdirections Matrix [nsources,3] of the aim point coordinates for the sources.
tp@0 49 % Srotations Vector [nsources,1] of the source rotations, in degrees.
tp@0 50 % Slevels Matrix [nsources,6/8] of the source levels for six (CAD-files v6) or eight
tp@0 51 % octave bands (CAD-files v7 or later).
tp@0 52 % Rnumbers Vector [nreceivers,1] of the receiver numbers in the CAD-file.
tp@0 53 % Rcoords Matrix [nreceivers,3] of the receiver coordinates.
tp@0 54 % CATTversionnumber 6,7 or 8
tp@0 55 % planeeqs Matrix [nplanes,4] of the plane equations as derived from the plane definitions.
tp@0 56 % Each row has the values [A B C D] of the plane equation on the form Ax + By + Cz = D
tp@0 57 % planenvecs Matrix [nplanes,3] of the normalized normal vectors of the planes.
tp@0 58 % ncornersperplanevec Vector [nplanes,1] which gives the number of corners for each plane.
tp@0 59 % planesatcorners Matrix [ncorners,4] which gives the plane numbers that each corner is connected to.
tp@0 60 % NB! Here it is assumed that each corner is connected to maximum four planes.
tp@0 61 % nplanespercorners Vector [ncorners,1] which gives the number of planes that each corner is connected to.
tp@0 62 % minvals Matrix [nplanes,3] which for each plane gives the smallest coordinate in the x-, y- and z-direction.
tp@0 63 % maxvals Matrix [nplanes,3] which for each plane gives the largest coordinate in the x-, y- and z-direction.
tp@0 64 % planehasindents Vector [nplanes,1] which for each plane gives the
tp@0 65 % number of indeting corners
tp@0 66 % indentingcorners Matrix [nplanes,max(ncornersperplanevec)] which for each plane gives the number to the first corner
tp@0 67 % in a corner triplet which identifies an indenting corner. The number of the corner is the
tp@0 68 % order given in planecorners for that plane.
tp@0 69 %
tp@0 70 % Uses the functions EDB1fistr, EDB1extrnums, EDB1cross
tp@0 71 %
tp@0 72 % ----------------------------------------------------------------------------------------------
tp@0 73 % This file is part of the Edge Diffraction Toolbox by Peter Svensson.
tp@0 74 %
tp@0 75 % The Edge Diffraction Toolbox is free software: you can redistribute it and/or modify
tp@0 76 % it under the terms of the GNU General Public License as published by the Free Software
tp@0 77 % Foundation, either version 3 of the License, or (at your option) any later version.
tp@0 78 %
tp@0 79 % The Edge Diffraction Toolbox is distributed in the hope that it will be useful,
tp@0 80 % but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
tp@0 81 % FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
tp@0 82 %
tp@0 83 % You should have received a copy of the GNU General Public License along with the
tp@0 84 % Edge Diffraction Toolbox. If not, see <http://www.gnu.org/licenses/>.
tp@0 85 % ----------------------------------------------------------------------------------------------
tp@0 86 % Peter Svensson (svensson@iet.ntnu.no) 18Jul09
tp@0 87 %
tp@0 88 % ----------------------------------------------------------------------------------------------
tp@0 89 % This File is created based on EDB1eadcad.m to support the Edge Diffraction Toolbox
tp@0 90 % ----------------------------------------------------------------------------------------------
tp@0 91 % Alexander Pohl (alexander.pohl@hcu-hamburg.de) 25May2011
tp@0 92 %
tp@0 93 % [outputfile] = EDB1readac(CADfile,outputfile,planecornerstype,checkgeom);
tp@0 94
tp@0 95 % Modified by Peter Svensson on 3Dec2011
tp@0 96
tp@0 97 global SHOWTEXT
tp@0 98
tp@0 99 if nargin == 0
tp@0 100 CADfile = '';
tp@0 101 outputfile = '';
tp@0 102 planecornerstype = '';
tp@0 103 checkgeom = '';
tp@0 104 elseif nargin == 1
tp@0 105 outputfile = '';
tp@0 106 planecornerstype = '';
tp@0 107 checkgeom = '';
tp@0 108 elseif nargin == 2
tp@0 109 planecornerstype = '';
tp@0 110 checkgeom = '';
tp@0 111 elseif nargin == 3
tp@0 112 checkgeom = '';
tp@0 113 end
tp@0 114
tp@0 115 geomacc = 1e-10;
tp@0 116
tp@0 117 %---------------------------------------------------------------
tp@0 118 % If no AC-file was specified, present a file opening window
tp@0 119
tp@0 120 if isempty(CADfile)
tp@0 121 [CADfile,filepath] = uigetfile('*.ac','Please select the AC3Dfile');
tp@0 122
tp@0 123 [filepath,temp1] = fileparts(filepath);
tp@0 124
tp@0 125 if ~isstr(CADfile)
tp@0 126 return
tp@0 127 end
tp@0 128
tp@0 129 [temp1,filestem] = fileparts(CADfile);
tp@0 130
tp@0 131
tp@0 132 CADfile = [[filepath,filesep],filestem,'.ac'];
tp@0 133 else
tp@0 134
tp@0 135 [filepath,filestem,fileext] = fileparts(CADfile);
tp@0 136 CADfile = [[filepath,filesep],filestem,fileext];
tp@0 137 end
tp@0 138
tp@0 139
tp@0 140 if exist(CADfile) ~= 2
tp@0 141 error(['ERROR: AC3D-file: ',CADfile,' can not be found'])
tp@0 142 end
tp@0 143
tp@0 144 %---------------------------------------------------------------
tp@0 145 % If no output file was specified, construct an automatic file name
tp@0 146
tp@0 147 if isempty(outputfile)
tp@0 148 outputfile = [[filepath,filesep],filestem,'_cadgeo.mat'];
tp@0 149 end
tp@0 150
tp@0 151 %---------------------------------------------------------------
tp@0 152 % Read in the entire file into the string B
tp@0 153 % Find the line starts and ends
tp@0 154
tp@0 155 fid = fopen([CADfile],'r');
tp@0 156 if fid == 0
tp@0 157 error(['ERROR: The AC3Dfile ',CADfile,' could not be opened.'])
tp@0 158 end
tp@0 159
tp@0 160 %---------------------------------------------------------------
tp@0 161 % read header "AC3Db"
tp@0 162 tline = fgetl(fid);
tp@0 163 if (~strcmp(tline,'AC3Db'))
tp@0 164 disp('AC3Db header missing');
tp@0 165 end
tp@0 166
tp@0 167 % read materials
tp@0 168 p_NumberOfMaterials=0;
tp@0 169 while(true)
tp@0 170 tline = fgetl(fid);
tp@0 171 tokens = readLineToToken(tline);
tp@0 172
tp@0 173 if(strcmp(tokens{1},'MATERIAL'))
tp@0 174 p_MaterialList{p_NumberOfMaterials+1}=tokens{2};
tp@0 175 p_NumberOfMaterials=p_NumberOfMaterials+1;
tp@0 176
tp@0 177 else
tp@0 178 break;
tp@0 179 end
tp@0 180 end
tp@0 181
tp@0 182 % read world object (all objects are kids of this)
tp@0 183 m_LocalMatrix = [1.0 0.0 0.0 0.0;
tp@0 184 0.0 1.0 0.0 0.0;
tp@0 185 0.0 0.0 1.0 0.0;
tp@0 186 0.0 0.0 0.0 1.0];
tp@0 187
tp@0 188 % Read OBJECT world
tp@0 189 tokens = readLineToToken(tline);
tp@0 190 if(~strcmp(tokens{1},'OBJECT'))
tp@0 191 disp('Section does not start with "object"');
tp@0 192 end
tp@0 193 if(~strcmp(tokens{2},'world'))
tp@0 194 disp('Section does not start with "world"');
tp@0 195 end
tp@0 196
tp@0 197 % Read kids
tp@0 198 tline = fgetl(fid);
tp@0 199 tokens = readLineToToken(tline);
tp@0 200 if(~strcmp(tokens{1},'kids'))
tp@0 201 disp('root object section does not start with kids');
tp@0 202 end
tp@0 203 p_Kids = str2num(tokens{2});
tp@0 204
tp@0 205 for i = 1:p_Kids
tp@0 206 p_Objects{i}=read_ac3d_object(fid, m_LocalMatrix);
tp@0 207 end
tp@0 208
tp@0 209 fclose(fid);
tp@0 210
tp@0 211 %% Export
tp@0 212 eval(['save a.mat p_Objects']);
tp@0 213 % Combine Vertices of EVERY Object to on corner list (incl. Counting PlanesPerPoly)
tp@0 214 for i = 1:p_Kids
tp@0 215 if(i==1)
tp@0 216 corners=p_Objects{i}{1}';
tp@0 217 p_CornersPerPoly = size(p_Objects{i}{1},2);
tp@0 218 p_PlanesPerPoly = size(p_Objects{i}{2},2);
tp@0 219 else
tp@0 220 corners = [corners; p_Objects{i}{1}'];
tp@0 221 p_PlanesPerPoly = [p_PlanesPerPoly size(p_Objects{i}{2},2)];
tp@0 222 p_CornersPerPoly = [p_CornersPerPoly size(p_Objects{i}{1},2)];
tp@0 223 end
tp@0 224 end
tp@0 225
tp@0 226 % Enumerate each corner from 1:end
tp@0 227 ncorners = size(corners,1);
tp@0 228 cornernumbers=(1:ncorners)';
tp@0 229
tp@0 230 %Number of Planes is sum over all PlanesPerPoly
tp@0 231 nplanes=sum(p_PlanesPerPoly);
tp@0 232
tp@0 233 ncornersperplanevec=zeros(nplanes,1);
tp@0 234 currentPlaneIndex=1;
tp@0 235 for i = 1:p_Kids
tp@0 236 for j=1:p_PlanesPerPoly(i);
tp@0 237 ncornersperplanevec(currentPlaneIndex)=length(p_Objects{i}{2}{j});
tp@0 238 currentPlaneIndex=currentPlaneIndex+1;
tp@0 239 end
tp@0 240 end
tp@0 241
tp@0 242 % We need to compute the maximum number of vertices for any plane
tp@0 243 % (others are filled up with 0)
tp@0 244 p_MaximumNumberOfVerticesOfPlane = max(ncornersperplanevec);
tp@0 245
tp@0 246 % Compute indices of corners for each plane
tp@0 247 % Simple for first Object, as here numbers are identical
tp@0 248 % For every further object, the number of previous used corners has be be
tp@0 249 % taken as offset
tp@0 250 planecorners=zeros(nplanes,p_MaximumNumberOfVerticesOfPlane);
tp@0 251 currentPlaneIndex=1;
tp@0 252 p_Offset=0;
tp@0 253 for i = 1:p_Kids
tp@0 254 for j=1:p_PlanesPerPoly(i);
tp@0 255 for k=1:length(p_Objects{i}{2}{j})
tp@0 256 planecorners(currentPlaneIndex,k)=p_Objects{i}{2}{j}(k)+p_Offset;
tp@0 257 end
tp@0 258 currentPlaneIndex=currentPlaneIndex+1;
tp@0 259 end
tp@0 260 p_Offset=p_Offset+p_CornersPerPoly(i);
tp@0 261 end
tp@0 262
tp@0 263 planenames=zeros(nplanes,30);
tp@0 264 currentPlaneIndex=1;
tp@0 265 for i = 1:p_Kids
tp@0 266 for j=1:p_PlanesPerPoly(i);
tp@0 267 p_Text = sparse(double(p_MaterialList{p_Objects{i}{4}{j}+1}));
tp@0 268
tp@0 269 if length(p_Text) > size(planenames,2)
tp@0 270 planenames = [planenames zeros(planenumbers(currentPlaneIndex),length(p_Text)-size(planenames,2))];
tp@0 271 end
tp@0 272 planenames(currentPlaneIndex,1:length(p_Text)) = p_Text;
tp@0 273 currentPlaneIndex=currentPlaneIndex+1;
tp@0 274 end
tp@0 275 end
tp@0 276
tp@0 277 planeabstypes = sparse(planenames(:,1:6));
tp@0 278
tp@0 279 % Addition on 2.12.2011 by Peter Svensson
tp@0 280 % The absorber names had an initial " which must be removed.
tp@0 281 % We remove columns until A-Z,a-z
tp@0 282
tp@0 283 if size(planeabstypes,1) > 1
tp@0 284 columnhasalphabet = 1 - sign( prod(sign(64-full(planeabstypes))+1));
tp@0 285 else
tp@0 286 columnhasalphabet = 1 - sign( (sign(64-full(planeabstypes))+1));
tp@0 287 end
tp@0 288 firstOKcolumn = find(columnhasalphabet);
tp@0 289 if isempty(firstOKcolumn)
tp@0 290 error('ERROR: All plane absorber types have names without any letters. They must start with a letter')
tp@0 291 end
tp@0 292 firstOKcolumn = firstOKcolumn(1);
tp@0 293 planeabstypes = planeabstypes(:,firstOKcolumn:end);
tp@0 294
tp@0 295
tp@0 296 planenumbers=1:nplanes;
tp@0 297
tp@0 298 eval(['save ',outputfile,' corners planecorners ncornersperplanevec'])
tp@0 299
tp@0 300 %% Begining from here, it is code based of EDB1readcad.m without changes - except SOURCE - RECIEVER Values, as not in AC3Dfile
tp@0 301 %---------------------------------------------------------------
tp@0 302 % Change the numbering in the CAD file to a contiguous one
tp@0 303
tp@0 304 % Corner numbers
tp@0 305 % First we make an inelegant crossreference vector giving the
tp@0 306 % resulting corner number in the position given by the CAD-file
tp@0 307 % corner number.
tp@0 308
tp@0 309 % PS 20120414: CHECK Corners with number zero are not handled correctly.
tp@0 310 % Insert code from EDB1readcad here.
tp@0 311
tp@0 312 maxnumb = max(cornernumbers);
tp@0 313 cornercrossref = sparse(zeros(1,maxnumb));
tp@0 314 cornercrossref(cornernumbers) = [1:ncorners];
tp@0 315
tp@0 316 [nplanes,ncornersperplane] = size(planecorners);
tp@0 317 iv = find(planecorners~=0);
tp@0 318 planecorners(iv) = full(cornercrossref(planecorners(iv)));
tp@0 319
tp@0 320 %---------------------------------------------------------------
tp@0 321 % Go through all planes. If there is a plane definition including
tp@0 322 % zeros, and planecornerstype == 'circ', expand it repeating the
tp@0 323 % same corner order again.
tp@0 324
tp@0 325 if isempty(planecornerstype)
tp@0 326 planecornerstype = 'circ';
tp@0 327 else
tp@0 328 planecornerstype = setstr(lower(planecornerstype(1)));
tp@0 329 if planecornerstype(1) == 'z'
tp@0 330 planecornerstype = 'zero';
tp@0 331 else
tp@0 332 planecornerstype = 'circ';
tp@0 333 end
tp@0 334 end
tp@0 335
tp@0 336 if planecornerstype == 'circ'
tp@0 337 for ii = 1:nplanes
tp@0 338 iv = find( planecorners(ii,:) ~= 0);
tp@0 339 ncornersatplane = length(iv);
tp@0 340 if ncornersatplane ~= ncornersperplane
tp@0 341 pattern = planecorners(ii,iv);
tp@0 342 nrepeatings = ceil(ncornersperplane/ncornersatplane);
tp@0 343 for jj = 1:nrepeatings-1
tp@0 344 pattern = [pattern planecorners(ii,iv)];
tp@0 345 end
tp@0 346 planecorners(ii,:) = pattern(1:ncornersperplane);
tp@0 347 end
tp@0 348 end
tp@0 349 end
tp@0 350
tp@0 351
tp@0 352 %---------------------------------------------------------------
tp@0 353 % Find the normal vectors to the planes using the cross products
tp@0 354 %
tp@0 355 % The normal vector is basically the cross product between two vectors
tp@0 356 % connecting three consecutive points. If there are indents though
tp@0 357 % they will cause reversed normal vectors, so one must go through all
tp@0 358 % combinations and choose the majority normal vector.
tp@0 359 %
tp@0 360 % 26mar09 Use the fact described above for detecting indention corners
tp@0 361
tp@0 362 planenvecs = zeros(nplanes,3);
tp@0 363 planehasindents = zeros(nplanes,1);
tp@0 364 indentingcorners = sparse(zeros(nplanes,max(ncornersperplanevec)));
tp@0 365
tp@0 366 for ii = 1:nplanes
tp@0 367
tp@0 368 iv = find(planecorners(ii,:)~=0);
tp@0 369 cornerlist = planecorners(ii,iv);
tp@0 370 iv = find(cornerlist == cornerlist(1));
tp@0 371 if length(iv) > 1
tp@0 372 cornerlist = cornerlist(1:iv(2)-1);
tp@0 373 end
tp@0 374 ncorners = length( cornerlist );
tp@0 375 cornerlist = [cornerlist cornerlist(1) cornerlist(2)];
tp@0 376
tp@0 377 nvectorlist = zeros(ncorners,3);
tp@0 378 nveclen = zeros(ncorners,1);
tp@0 379
tp@0 380 for jj = 1:ncorners
tp@0 381 co1numb = cornerlist(jj);
tp@0 382 co2numb = cornerlist(jj+1);
tp@0 383 co3numb = cornerlist(jj+2);
tp@0 384 vec1 = corners(co1numb,:) - corners(co2numb,:);
tp@0 385 vec2 = corners(co3numb,:) - corners(co2numb,:);
tp@0 386
tp@0 387 nvec = EDB1cross(vec1.',vec2.').';
tp@0 388 nveclen(jj) = norm(nvec);
tp@0 389 if nveclen(jj) > 0
tp@0 390 nvectorlist(jj,:) = nvec./nveclen(jj);
tp@0 391 end
tp@0 392 end
tp@0 393
tp@0 394 iv = find(nveclen < max(nveclen)*0.001);
tp@0 395 nvectorlist(iv,:) = [];
tp@0 396 nvecref = nvectorlist(1,:);
tp@0 397
tp@0 398 [n1,slask] = size(nvectorlist);
tp@0 399 nvecsigns = round(sum( (nvectorlist.').*(nvecref(ones(n1,1),:).') ));
tp@0 400
tp@0 401 if sum(nvecsigns) == 0
tp@0 402 disp(' ')
tp@0 403 error(['ERROR: Plane ',int2str(planenumbers(ii)),' (plane numbering as in the CAD file) seems to be twisted.'])
tp@0 404 end
tp@0 405
tp@0 406
tp@0 407 if abs(sum(nvecsigns)) ~= n1
tp@0 408 disp(['Plane ',int2str(ii),' has indents!'])
tp@0 409 nindents = (n1 - abs(sum(nvecsigns)))/2;
tp@0 410 planehasindents(ii) = nindents;
tp@0 411 ivindent = find(nvecsigns == -1);
tp@0 412 if length(ivindent) == nindents
tp@0 413 indentingcorners(ii,ivindent) = 1;
tp@0 414 else
tp@0 415 if length(ivindent) == ncorners - nindents
tp@0 416 ivindent = find(nvecsigns == 1);
tp@0 417 indentingcorners(ii,ivindent) = 1;
tp@0 418 else
tp@0 419 error(['ERROR: An unexpected problem. Please report to the developer'])
tp@0 420 end
tp@0 421 end
tp@0 422
tp@0 423 end
tp@0 424
tp@0 425 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 426
tp@0 427 if n1 > 2
tp@0 428 nvecdiff = sum(nvecdiff.'.^2).';
tp@0 429 else
tp@0 430 nvecdiff = norm(nvecdiff);
tp@0 431 end
tp@0 432
tp@0 433 %APO Change
tp@0 434 if any(nvecdiff>1e-3),
tp@0 435 nvecdiff
tp@0 436 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 437 elseif any(nvecdiff>1e-8)
tp@0 438 nvecdiff
tp@0 439 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 440 end
tp@0 441
tp@0 442 if ncorners > 5 & abs(sum(nvecsigns)) <= 1
tp@0 443 disp(['WARNING for plane number ',int2str(planenumbers(ii)),' in the CAD-file'])
tp@0 444 disp([' with the name ',strtrim(setstr(full(planenames(ii,:))))])
tp@0 445 disp(' The normal vector can not be determined for this plane because there are')
tp@0 446 disp(' the same number of inwards and outwards corners')
tp@0 447 disp(' This is a list of the possible normal vectors:')
tp@0 448 [nv1,slask] = size(nvectorlist);
tp@0 449 for kk = 1:nv1
tp@0 450 vecstr = [' ',int2str(kk),'. ',num2str(-nvectorlist(kk,1)),' ',num2str(-nvectorlist(kk,2)),' ',num2str(-nvectorlist(kk,3))];
tp@0 451 disp(vecstr)
tp@0 452 end
tp@0 453 disp(' ')
tp@0 454
tp@0 455 preferredsign = input(' Give the number of a correct normal vector for this plane please ');
tp@0 456 switchsign = nvecsigns.'./nvecsigns(preferredsign);
tp@0 457 nvectorlist = nvectorlist.*switchsign(:,ones(1,3));
tp@0 458 else
tp@0 459 mostcommonsign = sign(sum(nvecsigns));
tp@0 460
tp@0 461 switchsign = nvecsigns.'./mostcommonsign;
tp@0 462 nvectorlist = nvectorlist.*switchsign(:,ones(1,3));
tp@0 463 end
tp@0 464
tp@0 465 planenvecs(ii,:) = mean(nvectorlist);
tp@0 466 end
tp@0 467
tp@0 468 planenvecs = -planenvecs;
tp@0 469
tp@0 470 %---------------------------------------------------------------
tp@0 471 % Plane equations, as Ax + By + Cz = D for each plane
tp@0 472
tp@0 473 planeeqs = zeros(nplanes,4);
tp@0 474 planeeqs(:,1:3) = planenvecs;
tp@0 475 planeeqs(:,4) = sum( (planenvecs.').*(corners(planecorners(:,1),:).') ).';
tp@0 476
tp@0 477 %---------------------------------------------------------------
tp@0 478 % Useful data: planesatcorners, minvals and maxvals
tp@0 479
tp@0 480 [ncorners,slask] = size(corners);
tp@0 481 planesatcornerhits = zeros(ncorners,nplanes);
tp@0 482
tp@0 483 for ii = 1:nplanes
tp@0 484 cornerlist = planecorners(ii,1:ncornersperplanevec(ii));
tp@0 485 planesatcornerhits(cornerlist,ii) = planesatcornerhits(cornerlist,ii) + 1;
tp@0 486 end
tp@0 487
tp@0 488 maxplanespercorner = 0;
tp@0 489 for ii = 1:ncorners
tp@0 490 nplanes = length(find(planesatcornerhits(ii,:) ~= 0));
tp@0 491 if nplanes > maxplanespercorner
tp@0 492 maxplanespercorner = nplanes;
tp@0 493 end
tp@0 494 end
tp@0 495
tp@0 496 planesatcorners = zeros(ncorners,maxplanespercorner);
tp@0 497 nplanespercorners = zeros(ncorners,1);
tp@0 498 for ii = 1:ncorners
tp@0 499 iv = find(planesatcornerhits(ii,:)~=0);
tp@0 500 planesatcorners(ii,1:length(iv)) = iv;
tp@0 501 nplanespercorners(ii) = length(iv);
tp@0 502 end
tp@0 503
tp@0 504 % find cubic boxes inside which the planes are placed
tp@0 505
tp@0 506 [nplanes,slask] = size(planeeqs);
tp@0 507
tp@0 508 minvals = zeros(nplanes,3);
tp@0 509 maxvals = zeros(nplanes,3);
tp@0 510
tp@0 511 for ii = 1:nplanes
tp@0 512 cornerlist = planecorners(ii,:);
tp@0 513 cornerlist = cornerlist( find(cornerlist~= 0) );
tp@0 514 cornercoord = corners(cornerlist,:);
tp@0 515 minvals(ii,:) = min(cornercoord);
tp@0 516 maxvals(ii,:) = max(cornercoord);
tp@0 517 end
tp@0 518
tp@0 519 minvals = minvals - geomacc;
tp@0 520 maxvals = maxvals + geomacc;
tp@0 521
tp@0 522 %---------------------------------------------------------------
tp@0 523 % Check the geometry a bit
tp@0 524
tp@0 525 if isempty(checkgeom)
tp@0 526 checkgeom = 0;
tp@0 527 else
tp@0 528 if setstr(checkgeom(1)) == 'c'
tp@0 529 checkgeom = 1;
tp@0 530 else
tp@0 531 checkgeom = 0;
tp@0 532 end
tp@0 533 end
tp@0 534
tp@0 535 if checkgeom
tp@0 536
tp@0 537 badproblem = 0;
tp@0 538 disp(' ')
tp@0 539 disp(' Checking if there are corners with identical coordinates')
tp@0 540 for ii = 1:ncorners-1
tp@0 541 othercorners = [ii+1:ncorners];
tp@0 542 iv = find((corners(othercorners,1)==corners(ii,1)) & (corners(othercorners,2)==corners(ii,2)) & (corners(othercorners,3)==corners(ii,3)));
tp@0 543 if ~isempty(iv)
tp@0 544 disp([' WARNING: corner ',int2str(ii),' and ',int2str(othercorners(iv(1))),' have the same coordinates'])
tp@0 545 disp([' This should be fixed in the CAD-file or in the CATT-model'])
tp@0 546 badproblem = 1;
tp@0 547 end
tp@0 548 end
tp@0 549
tp@0 550 disp(' Checking if there are unused corners')
tp@0 551 iv = find(planesatcorners(:,1)==0);
tp@0 552 if ~isempty(iv)
tp@0 553 disp(' WARNING: The following corners in the CAD-file are not used:')
tp@0 554 printvec = int2str(cornernumbers(iv(1)));
tp@0 555 for iiprint = 2:length(iv)
tp@0 556 printvec = [printvec,' ',int2str(cornernumbers(iv(iiprint)))];
tp@0 557 end
tp@0 558 disp([' ',printvec])
tp@0 559 disp([' This is not a big problem'])
tp@0 560 planesatcorners(iv,:) = [];
tp@0 561 end
tp@0 562
tp@0 563 disp(' Checking if any corners seem redundant')
tp@0 564 iv = find(planesatcorners(:,2)==0);
tp@0 565 if ~isempty(iv)
tp@0 566 disp(' WARNING: The following corners in the CAD-file belong to only one plane:')
tp@0 567 printvec = int2str(cornernumbers(iv(1)));
tp@0 568 for iiprint = 2:length(iv)
tp@0 569 printvec = [printvec,' ',int2str(cornernumbers(iv(iiprint)))];
tp@0 570 end
tp@0 571 disp([' ',printvec])
tp@0 572 disp([' This should be fixed in the CAD-file or in the CATT-model'])
tp@0 573 badproblem = 1;
tp@0 574 end
tp@0 575
tp@0 576 if badproblem == 1
tp@0 577 disp(' ');
tp@0 578 error(['ERROR: Problems with the geometry defined in the CAD-file: ',CADfile])
tp@0 579 end
tp@0 580
tp@0 581 end
tp@0 582
tp@0 583 %---------------------------------------------------------------
tp@0 584 % Save the relevant variables in the output file
tp@0 585
tp@0 586 % NB! Sparse matrices can not get a non-double format
tp@0 587
tp@0 588 if ncorners < 256
tp@0 589 planecorners = uint8(planecorners);
tp@0 590 elseif ncorners < 65536
tp@0 591 planecorners = uint16(planecorners);
tp@0 592 end
tp@0 593
tp@0 594 if max(ncornersperplanevec) <= 255
tp@0 595 ncornersperplanevec = uint8(ncornersperplanevec);
tp@0 596 planehasindents = uint8(planehasindents);
tp@0 597 else
tp@0 598 ncornersperplanevec = uint16(ncornersperplanevec);
tp@0 599 planehasindents = uint16(planehasindents);
tp@0 600 end
tp@0 601
tp@0 602 Varlist = [' corners cornernumbers planecorners planenames planeabstypes '];
tp@0 603 Varlist = [Varlist,' planenumbers cornercrossref '];
tp@0 604 Varlist = [Varlist,' planeeqs planenvecs ncornersperplanevec planesatcorners nplanespercorners'];
tp@0 605 Varlist = [Varlist,' minvals maxvals planehasindents indentingcorners'];
tp@0 606 %Varlist = [Varlist,' Scoords Sdirections Srotations Slevels Rnumbers Rcoords CATTversionnumber'];
tp@0 607
tp@0 608 planenames = sparse(planenames+1-1);
tp@0 609 planeabstypes = sparse(planeabstypes+1-1);
tp@0 610 %Sdirectivitynames = sparse(Sdirectivitynames+1-1);
tp@0 611
tp@0 612 eval(['save ',outputfile,Varlist])
tp@0 613
tp@0 614 end
tp@0 615
tp@0 616 function [p_Plane, p_Norm, p_Mat] = read_ac3d_surface(fid, i_Vertices)
tp@0 617 tline = fgetl(fid);
tp@0 618 tokens = readLineToToken(tline);
tp@0 619 if(~strcmp(tokens{1},'SURF'))
tp@0 620 disp(['surface section does not start with SURF',tokens{1}]);
tp@0 621 end
tp@0 622 if(~strcmp(tokens{2},'0x10'))
tp@0 623 disp(['only polygons are accepted TODO?', tokens{2}]);
tp@0 624 end
tp@0 625
tp@0 626 tline = fgetl(fid);
tp@0 627 tokens = readLineToToken(tline);
tp@0 628 if(~strcmp(tokens{1},'mat'))
tp@0 629 disp('missing expected token "mat"');
tp@0 630 end
tp@0 631 p_Mat=str2num(tokens{2});
tp@0 632
tp@0 633 tline = fgetl(fid);
tp@0 634 tokens = readLineToToken(tline);
tp@0 635
tp@0 636 if(~strcmp(tokens{1},'refs'))
tp@0 637 disp('surface section missing refs');
tp@0 638 end
tp@0 639
tp@0 640 p_NumberOfRefs = str2num(tokens{2});
tp@0 641
tp@0 642 if(tokens{2}<3)
tp@0 643 disp(['too less points in surface:',tokens{2}]);
tp@0 644 return;
tp@0 645 end
tp@0 646
tp@0 647 p_Plane=zeros(p_NumberOfRefs,1);
tp@0 648
tp@0 649 for r=1:p_NumberOfRefs
tp@0 650 tline = fgetl(fid);
tp@0 651 tokens = readLineToToken(tline);
tp@0 652 p_Plane(r)=str2num(tokens{1})+1;
tp@0 653 end
tp@0 654
tp@0 655 v1 = i_Vertices(:,p_Plane(2))-i_Vertices(:,p_Plane(1));
tp@0 656 v2 = i_Vertices(:,p_Plane(3))-i_Vertices(:,p_Plane(1));
tp@0 657 p_Norm = cross(v1,v2);
tp@0 658 p_Norm = p_Norm./norm(p_Norm);
tp@0 659 end
tp@0 660
tp@0 661 function [token] = readLineToToken(tline)
tp@0 662 [token{1}, remain] = strtok(tline);
tp@0 663 i=1;
tp@0 664 while(remain~=0)
tp@0 665 i=i+1;
tp@0 666 tline = remain;
tp@0 667 [token{i}, remain] = strtok(tline);
tp@0 668 end
tp@0 669 end
tp@0 670
tp@0 671 function [p_Object] = read_ac3d_object(fid, i_LocalMatrix)
tp@0 672
tp@0 673 tline = fgetl(fid);
tp@0 674 tokens = readLineToToken(tline);
tp@0 675 if(~strcmp(tokens{1},'OBJECT'))
tp@0 676 disp('Section does not start with "object"');
tp@0 677 end
tp@0 678 if(strcmp(tokens{2},'poly'))
tp@0 679 p_Object = read_ac3d_poly(fid, i_LocalMatrix);
tp@0 680 elseif(strcmp(tokens{2},'group'))
tp@0 681 disp('GROUP not implemented yet');
tp@0 682 elseif(strcmp(tokens{2},'light'))
tp@0 683 disp('LIGHT not implemented yet');
tp@0 684 else
tp@0 685 disp('Section does not start with valid key');
tp@0 686 end
tp@0 687 end
tp@0 688 function [p_Polygon] = read_ac3d_poly(fid, i_LocalMatrix)
tp@0 689
tp@0 690 p_LocalMatrix = [1.0 0.0 0.0 0.0;
tp@0 691 0.0 1.0 0.0 0.0;
tp@0 692 0.0 0.0 1.0 0.0;
tp@0 693 0.0 0.0 0.0 1.0];
tp@0 694
tp@0 695 tline = fgetl(fid);
tp@0 696
tp@0 697 p_Vertices=0;
tp@0 698 while(ischar(tline))
tp@0 699 tokens = readLineToToken(tline);
tp@0 700
tp@0 701 if(strcmp(tokens{1},'numvert'))
tp@0 702
tp@0 703 p_ResultMatrix = i_LocalMatrix * p_LocalMatrix;
tp@0 704
tp@0 705 p_NumVertices=str2num(tokens{2});
tp@0 706 p_Vertices=zeros(3,p_NumVertices);
tp@0 707 for l=1:p_NumVertices
tp@0 708 tline = fgetl(fid);
tp@0 709 tokens = readLineToToken(tline);
tp@0 710 if(~ischar(tline))
tp@0 711 disp(['unexpected EOF while reading vertex (', num2str(l),'/',num2str(p_NumVertices)]);
tp@0 712 end
tp@0 713 if(size(tokens,2)~=3)
tp@0 714 disp(['error while reading vertex coords (', num2str(l),'/',num2str(p_NumVertices)]);
tp@0 715 end
tp@0 716
tp@0 717 p_Vec = p_ResultMatrix * [str2num(tokens{1}); str2num(tokens{2}); str2num(tokens{3}); 1];
tp@0 718 p_Vertices(:,l) = p_Vec(1:3);
tp@0 719 end
tp@0 720
tp@0 721 elseif(strcmp(tokens{1},'numsurf'))
tp@0 722 p_NumSurfaces=str2num(tokens{2});
tp@0 723 for s=1:p_NumSurfaces
tp@0 724 [p_Planes{s}, p_Norms{s}, p_Mats{s}] = read_ac3d_surface(fid, p_Vertices);
tp@0 725 end
tp@0 726 elseif(strcmp(tokens{1},'loc'))
tp@0 727 p_LocalMatrix(1,4)=str2num(tokens{2});
tp@0 728 p_LocalMatrix(2,4)=str2num(tokens{3});
tp@0 729 p_LocalMatrix(3,4)=str2num(tokens{4});
tp@0 730 elseif(strcmp(tokens{1},'rot'))
tp@0 731 p_LocalMatrix(1,1)=str2num(tokens{2});
tp@0 732 p_LocalMatrix(1,2)=str2num(tokens{3});
tp@0 733 p_LocalMatrix(1,3)=str2num(tokens{4});
tp@0 734 p_LocalMatrix(2,1)=str2num(tokens{5});
tp@0 735 p_LocalMatrix(2,2)=str2num(tokens{6});
tp@0 736 p_LocalMatrix(2,3)=str2num(tokens{7});
tp@0 737 p_LocalMatrix(3,1)=str2num(tokens{8});
tp@0 738 p_LocalMatrix(3,2)=str2num(tokens{9});
tp@0 739 p_LocalMatrix(3,3)=str2num(tokens{10});
tp@0 740 elseif(strcmp(tokens{1},'kids'))
tp@0 741 if(str2num(tokens{2})>2)
tp@0 742 disp(['Warning: Kids of Polygon not handled yet:',tokens{2}])
tp@0 743 p_Polygon=0;
tp@0 744 return;
tp@0 745 end
tp@0 746 p_Polygon{1}=p_Vertices;
tp@0 747 p_Polygon{2}=p_Planes;
tp@0 748 p_Polygon{3}=p_Norms;
tp@0 749 p_Polygon{4}=p_Mats;
tp@0 750 return;
tp@0 751 else
tp@0 752 if(strcmp(tokens{1},'name'))
tp@0 753 % Only name of poly, irrelevant -->SKIP
tp@0 754 elseif(strcmp(tokens{1},'crease'))
tp@0 755 else
tp@0 756 disp(['Warning: Ignored Entry in AC3D File: ', tokens{1}]);
tp@0 757 end
tp@0 758 end
tp@0 759 tline = fgetl(fid);
tp@0 760 end
tp@0 761 end