Mercurial > hg > human-echolocation

annotate private/EDB1readcad.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] = EDB1readcad(CADfile,outputfile,planecornerstype,checkgeom)
tp@0 2 % EDB1readcad - Reads a file of type .CAD (made by e.g. CATT-Acoustic) and saves all the geometry data in a mat-file.
tp@0 3 % CAD-files v6,v7,v8 are read.
tp@0 4 %
tp@0 5 % Input parameters:
tp@0 6 % CADfile (optional) The input file, with or without the .CAD extension.
tp@0 7 % If this file is not specified, a file opening window will appear.
tp@0 8 % outputfile (optional) The output file. If not specified, it will be given the
tp@0 9 % same name as the CAD file with an '_cadgeo.mat' ending instead of the '.CAD'.
tp@0 10 % planecornerstype (optional) Could have the value 'zero' or 'circ'.Default: 'circ'.
tp@0 11 % Affects the matrix planecorners, see below.
tp@0 12 % checkgeom (optional) If this parameter is given the value 'check', then a few checks
tp@0 13 % of the geometry consistency will be done: a check for duplicate corners
tp@0 14 % for redundant corners and for corners that are connected to only one plane.
tp@0 15 % Only warnings are given. As default, no check is done.
tp@0 16 % SHOWTEXT (global) If this global parameter has the value 3 or higher, informative text will
tp@0 17 % be printed on the screen.
tp@0 18 %
tp@0 19 % Output parameters:
tp@0 20 % outputfile The name of the outputfile, generated either automatically or specified as
tp@0 21 % an input parameter.
tp@0 22 %
tp@0 23 % Output parameters in the outputfile:
tp@0 24 % corners Matrix [ncorners,3] with the corner coordinates
tp@0 25 % cornernumbers Vector [ncorners,1] with the corner numbers used in the CAD-file. That is
tp@0 26 % in the outputfile, the corners will effectively be renumbered from 1 to ncorners
tp@0 27 % in the order they appeared in the CAD-file.
tp@0 28 % planecorners Matrix [planes,nmaxcornersperplane] with the corner numbers that make up each plane.
tp@0 29 % Since planes can have different numbers of corners, the number of columns is the
tp@0 30 % maximum number of corners that any plane has. NB! The corner numbers are in the
tp@0 31 % renumbered system, not in the CAD-file system.
tp@0 32 % planenames Matrix [nplanes,nmaxcharacters1] (sparse), with the planenames in the CAD file.
tp@0 33 % planeabstypes Matrix [nplanes,nmaxcharacters2] (sparse), with the absorber names in the CAD file.
tp@0 34 % planenumbers Vector [nplanes,1] with the plane numbers used in the CAD file.
tp@0 35 % cornercrossref Vector [nmaxnumberinCADfile,1] (sparse), which gives the matlab-file corner
tp@0 36 % numbers for each CAD-file corner number. For instance, if a CAD-file corner
tp@0 37 % number 1200 is translated to corner number 72 in the matlab-file, then
tp@0 38 % cornercrossref(1000) = 72.
tp@0 39 % Snumbers Vector [nsources,1] of the source numbers in the CAD-file. NB! Only CAD-files
tp@0 40 % v6 use source numbers; later versions use source names instead of numbers.
tp@0 41 % For CAD-file v7 and v8, Snumbers is empty.
tp@0 42 % Snames Matrix [nsources,2] of the source names in the CAD-file. NB! Only CAD-files v7
tp@0 43 % or later use source names (such as 'A0' etc). Older versions use source numbers,
tp@0 44 % in which case Snames is empty.
tp@0 45 % Sdirectivitynames Matrix [nsources,nmaxcharacters3] (sparse) with the source directivity names
tp@0 46 % in the CAD file.
tp@0 47 % Sdelays Vector [nsources,1] of the source delays in the CAD-file, in milliseconds.
tp@0 48 % Scoords Matrix [nsources,3] of the source coordinates.
tp@0 49 % Sdirections Matrix [nsources,3] of the aim point coordinates for the sources.
tp@0 50 % Srotations Vector [nsources,1] of the source rotations, in degrees.
tp@0 51 % Slevels Matrix [nsources,6/8] of the source levels for six (CAD-files v6) or eight
tp@0 52 % octave bands (CAD-files v7 or later).
tp@0 53 % Rnumbers Vector [nreceivers,1] of the receiver numbers in the CAD-file.
tp@0 54 % Rcoords Matrix [nreceivers,3] of the receiver coordinates.
tp@0 55 % CATTversionnumber 6,7 or 8
tp@0 56 % planeeqs Matrix [nplanes,4] of the plane equations as derived from the plane definitions.
tp@0 57 % Each row has the values [A B C D] of the plane equation on the form Ax + By + Cz = D
tp@0 58 % planenvecs Matrix [nplanes,3] of the normalized normal vectors of the planes.
tp@0 59 % ncornersperplanevec Vector [nplanes,1] which gives the number of corners for each plane.
tp@0 60 % planesatcorners Matrix [ncorners,4] which gives the plane numbers that each corner is connected to.
tp@0 61 % NB! Here it is assumed that each corner is connected to maximum four planes.
tp@0 62 % nplanespercorners Vector [ncorners,1] which gives the number of planes that each corner is connected to.
tp@0 63 % minvals Matrix [nplanes,3] which for each plane gives the smallest coordinate in the x-, y- and z-direction.
tp@0 64 % maxvals Matrix [nplanes,3] which for each plane gives the largest coordinate in the x-, y- and z-direction.
tp@0 65 % planehasindents Vector [nplanes,1] which for each plane gives the
tp@0 66 % number of indeting corners
tp@0 67 % indentingcorners Matrix [nplanes,max(ncornersperplanevec)] which for each plane gives the number to the first corner
tp@0 68 % in a corner triplet which identifies an indenting corner. The number of the corner is the
tp@0 69 % order given in planecorners for that plane.
tp@0 70 %
tp@0 71 % Uses the functions EDB1extrnums, EDB1cross
tp@0 72 %
tp@0 73 % ----------------------------------------------------------------------------------------------
tp@0 74 % This file is part of the Edge Diffraction Toolbox by Peter Svensson.
tp@0 75 %
tp@0 76 % The Edge Diffraction Toolbox is free software: you can redistribute it and/or modify
tp@0 77 % it under the terms of the GNU General Public License as published by the Free Software
tp@0 78 % Foundation, either version 3 of the License, or (at your option) any later version.
tp@0 79 %
tp@0 80 % The Edge Diffraction Toolbox is distributed in the hope that it will be useful,
tp@0 81 % but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
tp@0 82 % FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
tp@0 83 %
tp@0 84 % You should have received a copy of the GNU General Public License along with the
tp@0 85 % Edge Diffraction Toolbox. If not, see <http://www.gnu.org/licenses/>.
tp@0 86 % ----------------------------------------------------------------------------------------------
tp@0 87 % Peter Svensson (svensson@iet.ntnu.no) 20090718
tp@0 88 %
tp@0 89 % [outputfile] = EDB1readcad(CADfile,outputfile,planecornerstype,checkgeom);
tp@0 90
tp@0 91 global SHOWTEXT
tp@0 92
tp@0 93 if nargin == 0
tp@0 94 CADfile = '';
tp@0 95 outputfile = '';
tp@0 96 planecornerstype = '';
tp@0 97 checkgeom = '';
tp@0 98 elseif nargin == 1
tp@0 99 outputfile = '';
tp@0 100 planecornerstype = '';
tp@0 101 checkgeom = '';
tp@0 102 elseif nargin == 2
tp@0 103 planecornerstype = '';
tp@0 104 checkgeom = '';
tp@0 105 elseif nargin == 3
tp@0 106 checkgeom = '';
tp@0 107 end
tp@0 108
tp@0 109 geomacc = 1e-10;
tp@0 110
tp@0 111 %---------------------------------------------------------------
tp@0 112 % If no CAD-file was specified, present a file opening window
tp@0 113
tp@0 114 if isempty(CADfile)
tp@0 115 [CADfile,filepath] = uigetfile('*.cad','Please select the CADfile');
tp@0 116 [filepath,temp1,temp2] = fileparts(filepath);
tp@0 117 if ~isstr(CADfile)
tp@0 118 return
tp@0 119 end
tp@0 120 [temp1,filestem,temp2] = fileparts(CADfile);
tp@0 121
tp@0 122 CADfile = [[filepath,filesep],filestem,'.cad'];
tp@0 123 else
tp@0 124 [filepath,filestem,fileext] = fileparts(CADfile);
tp@0 125 CADfile = [[filepath,filesep],filestem,fileext];
tp@0 126 end
tp@0 127
tp@0 128 if exist(CADfile) ~= 2
tp@0 129 error(['ERROR: CAD-file: ',CADfile,' can not be found'])
tp@0 130 end
tp@0 131
tp@0 132 %---------------------------------------------------------------
tp@0 133 % If no output file was specified, construct an automatic file name
tp@0 134
tp@0 135 if isempty(outputfile)
tp@0 136 outputfile = [[filepath,filesep],filestem,'_cadgeo.mat'];
tp@0 137 end
tp@0 138
tp@0 139 %---------------------------------------------------------------
tp@0 140 % Read in the entire file into the string B
tp@0 141 % Find the line starts and ends
tp@0 142
tp@0 143 fid = fopen([CADfile],'r');
tp@0 144 if fid == 0
tp@0 145 error(['ERROR: The CADfile ',CADfile,' could not be opened.'])
tp@0 146 end
tp@0 147 B = fread(fid,inf,'char').';
tp@0 148 fclose(fid);
tp@0 149
tp@0 150 % Cut out sequential spaces to save some space.
tp@0 151 iv = find(B(1:length(B)-1)==32 & B(2:length(B))==32);
tp@0 152 B(iv) = [];
tp@0 153
tp@0 154 % If the text file was generated on a Mac, each line has ASCII 13 at the end.
tp@0 155 % On a PC, each line ends with ASCII 13, and the next line starts with ASCII 10.
tp@0 156 % In Unix, each line ends with ASCII 10.
tp@0 157 %
tp@0 158 % Change this to a single ASCII 13
tp@0 159
tp@0 160 % New 17 Jul 09: use the builtin function regexprep instead:
tp@0 161
tp@0 162 % 1. Replace all ASCII 10 by ASCII 13
tp@0 163
tp@0 164 B = regexprep(setstr(B),'\n','\r');
tp@0 165
tp@0 166 % 2. Cut out sequential CRs (ASCII13) to save some space.
tp@0 167 iv = find(B(1:length(B)-1)==13 & B(2:length(B))==13);
tp@0 168 B(iv) = [];
tp@0 169
tp@0 170 % 3. Cut out sequences of CR, space, CR because they "survived" the
tp@0 171 % search for consecutive blanks.
tp@0 172 iv = find(B(1:length(B)-1)==32 & B(2:length(B))==13);
tp@0 173 B(iv) = [];
tp@0 174 iv = find(B(1:length(B)-1)==13 & B(2:length(B))==13);
tp@0 175 B(iv) = [];
tp@0 176
tp@0 177 % 4. Convert all text to lower case
tp@0 178 B = lower(setstr(B));
tp@0 179
tp@0 180 % Some special characters might give negative ASCII values
tp@0 181 % (Comment by PS 17Jul09: Will this ever happen??)
tp@0 182
tp@0 183 iv = find(B<0);
tp@0 184 if ~isempty(iv)
tp@0 185 B(iv) = B(iv)+256;
tp@0 186 end
tp@0 187
tp@0 188 %---------------------------------------------------------------
tp@0 189 % Look for the text 'CATT-Acoustic v' anywhere, and read the numerical
tp@0 190 % value right after this text.
tp@0 191
tp@0 192 iv = regexp(B,'catt-acoustic v');
tp@0 193
tp@0 194 if isempty(iv)
tp@0 195 CATTversionnumber = [];
tp@0 196 else
tp@0 197 CATTversionnumber = str2num(setstr(B(iv(1)+15)));
tp@0 198 end
tp@0 199
tp@0 200 %---------------------------------------------------------------
tp@0 201 % Find the four sections of the file:
tp@0 202 % %CORNERS %PLANES %SOURCES %RECEIVERS
tp@0 203 % (or small letters)
tp@0 204 % We assume that they come in this order
tp@0 205
tp@0 206 stringposCORNERS = regexp(B,'%corners');
tp@0 207 if isempty(stringposCORNERS)
tp@0 208 error('ERROR: The .cad file must contain the word %CORNERS');
tp@0 209 end
tp@0 210
tp@0 211 stringposPLANES = regexp(B,'%planes');
tp@0 212 if isempty(stringposPLANES)
tp@0 213 error('ERROR: The .cad file must contain the word %PLANES');
tp@0 214 end
tp@0 215
tp@0 216 stringposSOURCES = regexp(B,'%sources');
tp@0 217
tp@0 218 stringposRECEIVERS = regexp(B,'%receivers');
tp@0 219
tp@0 220 %---------------------------------------------------------------
tp@0 221 % Look for the corners definitions between the lines containing
tp@0 222 % %SOURCES and another keyword (%PLANES or %SOURCES or %RECEIVERS)
tp@0 223 % The lines should look like:
tp@0 224 % 1 -0.2000 0.34 2.3
tp@0 225
tp@0 226 if stringposPLANES < stringposCORNERS
tp@0 227 error('ERROR: Sorry, but in the .cad file, the %CORNERS section must come before the %PLANES section')
tp@0 228 end
tp@0 229
tp@0 230 C = textscan(B(stringposCORNERS:stringposPLANES-1),'%d%f%f%f','Headerlines',1);
tp@0 231 cornernumbers = (C{1});
tp@0 232 corners = [C{2} C{3} C{4}];
tp@0 233 ncorners = length(cornernumbers);
tp@0 234
tp@0 235 %---------------------------------------------------------------
tp@0 236 % Look for the planes definitions between the lines containing
tp@0 237 % %PLANES and %SOURCES, or the end
tp@0 238 % The lines should look like:
tp@0 239 % 1 / planename /RIGID
tp@0 240 % 1 4 3 2
tp@0 241 % 2 / /RIGID
tp@0 242
tp@0 243 if ~isempty(stringposSOURCES)
tp@0 244 if stringposSOURCES < stringposPLANES
tp@0 245 error('ERROR: Sorry, but in the .cad file, the %PLANES section must come before the %SOURCES section')
tp@0 246 end
tp@0 247 Str = B(stringposPLANES:stringposSOURCES-1);
tp@0 248 else
tp@0 249 Str = B(stringposPLANES:end);
tp@0 250 end
tp@0 251
tp@0 252 % To make it easier with textscan, we put the two types of line onto
tp@0 253 % a single line, with a new '/' inserted. Then all lines have the same
tp@0 254 % format.
tp@0 255
tp@0 256 % Find all CR and replace every second with a '/'
tp@0 257 ivendlines = find(Str==13);
tp@0 258 Str(ivendlines(2:2:end)) = '/';
tp@0 259
tp@0 260 C = textscan(Str,'%d%s%s%s','Headerlines',1,'Delimiter','/','BufSize',32768);
tp@0 261
tp@0 262 planenumbers = C{1};
tp@0 263 Str2 = C{2};
tp@0 264 Str3 = C{3};
tp@0 265 Str = C{4};
tp@0 266
tp@0 267 nplanes = size(Str,1);
tp@0 268 ncornersperplanevec = zeros(nplanes,1);
tp@0 269 planecorners = zeros(nplanes,4);
tp@0 270 planenames = zeros(nplanes,30);
tp@0 271 planeabstypes = zeros(nplanes,6);
tp@0 272
tp@0 273 for ii = 1:nplanes
tp@0 274 listofplanes = EDB1extrnums(Str{ii});
tp@0 275 ncornersperplanevec(ii) = length(listofplanes);
tp@0 276 if ncornersperplanevec(ii) > size(planecorners,2)
tp@0 277 planecorners = [planecorners zeros(nplanes, ncornersperplanevec(ii)-size(planecorners,2))];
tp@0 278 end
tp@0 279 planecorners(ii,1:ncornersperplanevec(ii)) = listofplanes;
tp@0 280
tp@0 281 txtstr2 = Str2{ii};
tp@0 282 if length(txtstr2) > size(planenames,2)
tp@0 283 planenames = [planenames zeros(nplanes,length(txtstr2)-size(planenames,2))];
tp@0 284 end
tp@0 285 planenames(ii,1:length(txtstr2)) = txtstr2;
tp@0 286
tp@0 287 txtstr3 = Str3{ii};
tp@0 288 if length(txtstr3) > size(planeabstypes,2)
tp@0 289 planeabstypes = [planeabstypes zeros(nplanes,length(txtstr3)-size(planeabstypes,2))];
tp@0 290 end
tp@0 291 planeabstypes(ii,1:length(txtstr3)) = txtstr3;
tp@0 292 end
tp@0 293
tp@0 294 if max(max(planecorners)) > max(cornernumbers)
tp@0 295 error('ERROR: One plane definition in the CAD-file used a higher corner number than was defined in the CORNERS section')
tp@0 296 end
tp@0 297
tp@0 298 %---------------------------------------------------------------
tp@0 299 % Look for the sources definitions between the lines containing
tp@0 300 % %SOURCES and %RECEIVERS, or the end
tp@0 301 % The lines should look like: (v6)
tp@0 302 % 0 OMNI:SD0
tp@0 303 % 0.0000000 0.0000000 1.0000000
tp@0 304 % 0.0000000 0.0000000 0.0000000
tp@0 305 % 85.0 88.0 91.0 94.0 97.0 100.0
tp@0 306 %
tp@0 307 % or: (v7)
tp@0 308 % A0 OMNI:SD0
tp@0 309 % 0.0000000 0.0000000 1.0000000
tp@0 310 % 0.0000000 0.0000000 0.0000000 15.00000
tp@0 311 % 85.0 88.0 91.0 94.0 97.0 100.0 : 103.0 106.0
tp@0 312
tp@0 313 if ~isempty(stringposSOURCES)
tp@0 314
tp@0 315 if ~isempty(stringposRECEIVERS)
tp@0 316 if stringposRECEIVERS < stringposSOURCES
tp@0 317 error('ERROR: Sorry, but in the .cad file, the %SOURCES section must come before the %RECEIVERS section')
tp@0 318 end
tp@0 319 Str = B(stringposSOURCES:stringposRECEIVERS-1);
tp@0 320 else
tp@0 321 Str = B(stringposSOURCES:end);
tp@0 322 end
tp@0 323
tp@0 324 % To make it easier with textscan, we put the four types of line onto
tp@0 325 % a single line, with new '/' inserted. Then all lines have the same
tp@0 326 % format.
tp@0 327 % Find all CR and replace number 2,3,4, 6,7,8, etc with a '/'
tp@0 328
tp@0 329 ivendlines = find(Str==13);
tp@0 330 ivkeep = [1:4:length(ivendlines)];
tp@0 331 ivremove = [1:length(ivendlines)];
tp@0 332 ivremove(ivkeep) = [];
tp@0 333 Str(ivendlines(ivremove)) = '/';
tp@0 334
tp@0 335 C = textscan(Str,'%s%s%s%s','Headerlines',1,'Delimiter','/');
tp@0 336
tp@0 337 nsources = size( C{1},1 );
tp@0 338 Scoords = zeros(nsources,3);
tp@0 339 Sdirections = zeros(nsources,3);
tp@0 340 Srotations = zeros(nsources,1);
tp@0 341 Slevels = zeros(nsources,8);
tp@0 342 Snumbers = zeros(nsources,1);
tp@0 343 Snames = zeros(nsources,4);
tp@0 344 Sdirectivitynames = ( zeros(nsources,30) );
tp@0 345 Sdelays = zeros(nsources,1);
tp@0 346
tp@0 347 Str1 = C{1};
tp@0 348 Str2 = C{2};
tp@0 349 Str3 = C{3};
tp@0 350 Str4 = C{4};
tp@0 351 for ii = 1:nsources
tp@0 352 Scoords(ii,:) = EDB1extrnums(Str2{ii});
tp@0 353 listofvalues3 = EDB1extrnums(Str3{ii});
tp@0 354 listofvalues4 = EDB1extrnums(Str4{ii});
tp@0 355 if ii==1
tp@0 356 if length(listofvalues3)==3
tp@0 357 CATTversionnumber = 6;
tp@0 358 elseif length(listofvalues3)==4
tp@0 359 CATTversionnumber = 7;
tp@0 360 else
tp@0 361 error('ERROR: The source data should have three or four values for the direction')
tp@0 362 end
tp@0 363 end
tp@0 364
tp@0 365 Sdirections(ii,1:3) = listofvalues3(1:3);
tp@0 366
tp@0 367 % The first field contains three items which are a mix of text and numerical values.
tp@0 368 % First we remove possible consecutive blanks.
tp@0 369 textstr = Str1{ii};
tp@0 370 iv = find(textstr(1:length(textstr)-1)==32 & textstr(2:length(textstr))==32);
tp@0 371 textstr(iv) = [];
tp@0 372 iv = find(textstr==32);
tp@0 373
tp@0 374
tp@0 375 if CATTversionnumber == 7
tp@0 376 Srotations(ii) = listofvalues3(4);
tp@0 377 Snames(ii,1:length(textstr(1:iv(1)-1))) = textstr(1:iv(1)-1);
tp@0 378 Sdirectivitynames(ii,1:iv(2)-iv(1)-1) = textstr(iv(1)+1:iv(2)-1);
tp@0 379 Sdelays(ii) = str2num( textstr(iv(2)+1:end) );
tp@0 380 else
tp@0 381 Snumbers(ii) = str2num( textstr(1:iv(1)-1) );
tp@0 382 Sdirectivitynames(ii,1:length(textstr(iv(1)+1:end))) = textstr(iv(1)+1:end);
tp@0 383 end
tp@0 384 Slevels(ii,1:length(listofvalues4)) = listofvalues4;
tp@0 385 end
tp@0 386 Snames = setstr(Snames);
tp@0 387 Sdirectivitynames = setstr(Sdirectivitynames);
tp@0 388
tp@0 389 else
tp@0 390
tp@0 391 Scoords = [];
tp@0 392 Sdirections = [];
tp@0 393 Srotations = [];
tp@0 394 Slevels = [];
tp@0 395 Snumbers = [];
tp@0 396 Snames = [];
tp@0 397 Sdirectivitynames = [];
tp@0 398 Sdelays = [];
tp@0 399
tp@0 400 end
tp@0 401
tp@0 402 %---------------------------------------------------------------
tp@0 403 % Look for the receivers definitions after the line containing
tp@0 404 % %RECEIVERS
tp@0 405 % The lines should look like:
tp@0 406 % 1 -0.2000 0.34 2.3
tp@0 407
tp@0 408 if ~isempty(stringposRECEIVERS)
tp@0 409 C = textscan(B(stringposRECEIVERS:end),'%d%f%f%f','Headerlines',1);
tp@0 410 Rnumbers = (C{1});
tp@0 411 Rcoords= [C{2} C{3} C{4}];
tp@0 412 else
tp@0 413 Rnumbers = [];
tp@0 414 Rcoords = [];
tp@0 415 end
tp@0 416
tp@0 417 %---------------------------------------------------------------
tp@0 418 % Change the numbering in the CAD file to a contiguous one
tp@0 419
tp@0 420 % Corner numbers
tp@0 421 % First we make an inelegant crossreference vector giving the
tp@0 422 % resulting corner number in the position given by the CAD-file
tp@0 423 % corner number.
tp@0 424
tp@0 425 maxnumb = max(cornernumbers);
tp@0 426 cornercrossref = sparse(zeros(1,maxnumb));
tp@0 427
tp@0 428 % Temporary fix: if one corner has the number zero, then we make a
tp@0 429 % separate crossref variable for that one.
tp@0 430 if cornernumbers(1) == 0
tp@0 431 cornercrossrefzero = 1;
tp@0 432 cornercrossref(cornernumbers(2:end)) = [2:ncorners];
tp@0 433 else
tp@0 434 cornercrossrefzero = [];
tp@0 435 cornercrossref(cornernumbers) = [1:ncorners];
tp@0 436 end
tp@0 437
tp@0 438 [nplanes,ncornersperplane] = size(planecorners);
tp@0 439
tp@0 440 planecorners_CATTnumbers = planecorners;
tp@0 441
tp@0 442 iv = find(planecorners_CATTnumbers~=0);
tp@0 443 planecorners(iv) = full(cornercrossref(planecorners(iv)));
tp@0 444 iv = find(planecorners_CATTnumbers==0);
tp@0 445 planecorners(iv) = cornercrossrefzero;
tp@0 446
tp@0 447 %---------------------------------------------------------------
tp@0 448 % Go through all planes. If there is a plane definition including
tp@0 449 % zeros, and planecornerstype == 'circ', expand it repeating the
tp@0 450 % same corner order again.
tp@0 451
tp@0 452 if isempty(planecornerstype)
tp@0 453 planecornerstype = 'circ';
tp@0 454 else
tp@0 455 planecornerstype = setstr(lower(planecornerstype(1)));
tp@0 456 if planecornerstype(1) == 'z'
tp@0 457 planecornerstype = 'zero';
tp@0 458 else
tp@0 459 planecornerstype = 'circ';
tp@0 460 end
tp@0 461 end
tp@0 462
tp@0 463 if planecornerstype == 'circ'
tp@0 464 for ii = 1:nplanes
tp@0 465 iv = find( planecorners(ii,:) ~= 0);
tp@0 466 ncornersatplane = length(iv);
tp@0 467 if ncornersatplane ~= ncornersperplane
tp@0 468 pattern = planecorners(ii,iv);
tp@0 469 nrepeatings = ceil(ncornersperplane/ncornersatplane);
tp@0 470 for jj = 1:nrepeatings-1
tp@0 471 pattern = [pattern planecorners(ii,iv)];
tp@0 472 end
tp@0 473 planecorners(ii,:) = pattern(1:ncornersperplane);
tp@0 474 end
tp@0 475 end
tp@0 476 end
tp@0 477
tp@0 478 %---------------------------------------------------------------
tp@0 479 % Find the normal vectors to the planes using the cross products
tp@0 480 %
tp@0 481 % The normal vector is basically the cross product between two vectors
tp@0 482 % connecting three consecutive points. If there are indents though
tp@0 483 % they will cause reversed normal vectors, so one must go through all
tp@0 484 % combinations and choose the majority normal vector.
tp@0 485 %
tp@0 486 % 26mar09 Use the fact described above for detecting indention corners
tp@0 487
tp@0 488 planenvecs = zeros(nplanes,3);
tp@0 489 planehasindents = zeros(nplanes,1);
tp@0 490 indentingcorners = sparse(zeros(nplanes,max(ncornersperplanevec)));
tp@0 491
tp@0 492 for ii = 1:nplanes
tp@0 493
tp@0 494 iv = find(planecorners(ii,:)~=0);
tp@0 495 cornerlist = planecorners(ii,iv);
tp@0 496 iv = find(cornerlist == cornerlist(1));
tp@0 497 if length(iv) > 1
tp@0 498 cornerlist = cornerlist(1:iv(2)-1);
tp@0 499 end
tp@0 500 ncorners = length( cornerlist );
tp@0 501 cornerlist = [cornerlist cornerlist(1) cornerlist(2)];
tp@0 502
tp@0 503 nvectorlist = zeros(ncorners,3);
tp@0 504 nveclen = zeros(ncorners,1);
tp@0 505
tp@0 506 for jj = 1:ncorners
tp@0 507 co1numb = cornerlist(jj);
tp@0 508 co2numb = cornerlist(jj+1);
tp@0 509 co3numb = cornerlist(jj+2);
tp@0 510 vec1 = corners(co1numb,:) - corners(co2numb,:);
tp@0 511 vec2 = corners(co3numb,:) - corners(co2numb,:);
tp@0 512 nvec = EDB1cross(vec1.',vec2.').';
tp@0 513 nveclen(jj) = norm(nvec);
tp@0 514 if nveclen(jj) > 0
tp@0 515 nvectorlist(jj,:) = nvec./nveclen(jj);
tp@0 516 end
tp@0 517 end
tp@0 518
tp@0 519 iv = find(nveclen < max(nveclen)*0.001);
tp@0 520 nvectorlist(iv,:) = [];
tp@0 521 nvecref = nvectorlist(1,:);
tp@0 522
tp@0 523 [n1,slask] = size(nvectorlist);
tp@0 524 nvecsigns = round(sum( (nvectorlist.').*(nvecref(ones(n1,1),:).') ));
tp@0 525
tp@0 526 if sum(nvecsigns) == 0
tp@0 527 disp(' ')
tp@0 528 error(['ERROR: Plane ',int2str(planenumbers(ii)),' (plane numbering as in the CAD file) seems to be twisted.'])
tp@0 529 end
tp@0 530
tp@0 531
tp@0 532 if abs(sum(nvecsigns)) ~= n1
tp@0 533 disp(['Plane ',int2str(ii),' has indents!'])
tp@0 534 nindents = (n1 - abs(sum(nvecsigns)))/2;
tp@0 535 planehasindents(ii) = nindents;
tp@0 536 ivindent = find(nvecsigns == -1);
tp@0 537 if length(ivindent) == nindents
tp@0 538 indentingcorners(ii,ivindent) = 1;
tp@0 539 else
tp@0 540 if length(ivindent) == ncorners - nindents
tp@0 541 ivindent = find(nvecsigns == 1);
tp@0 542 indentingcorners(ii,ivindent) = 1;
tp@0 543 else
tp@0 544 error(['ERROR: An unexpected problem. Please report to the developer'])
tp@0 545 end
tp@0 546 end
tp@0 547 end
tp@0 548
tp@0 549 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 550
tp@0 551 if n1 > 2
tp@0 552 nvecdiff = sum(nvecdiff.'.^2).';
tp@0 553 else
tp@0 554 nvecdiff = norm(nvecdiff);
tp@0 555 end
tp@0 556 if any(nvecdiff>1e-4),
tp@0 557 nvecdiff
tp@0 558 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 559 elseif any(nvecdiff>1e-8)
tp@0 560 nvecdiff
tp@0 561 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 562 end
tp@0 563
tp@0 564 if ncorners > 5 & abs(sum(nvecsigns)) <= 1
tp@0 565 disp(['WARNING for plane number ',int2str(planenumbers(ii)),' in the CAD-file'])
tp@0 566 disp([' with the name ',strtrim(setstr(full(planenames(ii,:))))])
tp@0 567 disp(' The normal vector can not be determined for this plane because there are')
tp@0 568 disp(' the same number of inwards and outwards corners')
tp@0 569 disp(' This is a list of the possible normal vectors:')
tp@0 570 [nv1,slask] = size(nvectorlist);
tp@0 571 for kk = 1:nv1
tp@0 572 vecstr = [' ',int2str(kk),'. ',num2str(-nvectorlist(kk,1)),' ',num2str(-nvectorlist(kk,2)),' ',num2str(-nvectorlist(kk,3))];
tp@0 573 disp(vecstr)
tp@0 574 end
tp@0 575 disp(' ')
tp@0 576
tp@0 577 preferredsign = input(' Give the number of a correct normal vector for this plane please ');
tp@0 578 switchsign = nvecsigns.'./nvecsigns(preferredsign);
tp@0 579 nvectorlist = nvectorlist.*switchsign(:,ones(1,3));
tp@0 580 else
tp@0 581 mostcommonsign = sign(sum(nvecsigns));
tp@0 582
tp@0 583 switchsign = nvecsigns.'./mostcommonsign;
tp@0 584 nvectorlist = nvectorlist.*switchsign(:,ones(1,3));
tp@0 585 end
tp@0 586
tp@0 587 planenvecs(ii,:) = mean(nvectorlist);
tp@0 588 end
tp@0 589
tp@0 590 planenvecs = -planenvecs;
tp@0 591
tp@0 592 %---------------------------------------------------------------
tp@0 593 % Plane equations, as Ax + By + Cz = D for each plane
tp@0 594
tp@0 595 planeeqs = zeros(nplanes,4);
tp@0 596 planeeqs(:,1:3) = planenvecs;
tp@0 597 planeeqs(:,4) = sum( (planenvecs.').*(corners(planecorners(:,1),:).') ).';
tp@0 598
tp@0 599 %---------------------------------------------------------------
tp@0 600 % Useful data: planesatcorners, minvals and maxvals
tp@0 601
tp@0 602 [ncorners,slask] = size(corners);
tp@0 603 planesatcornerhits = zeros(ncorners,nplanes);
tp@0 604
tp@0 605 for ii = 1:nplanes
tp@0 606 cornerlist = planecorners(ii,1:ncornersperplanevec(ii));
tp@0 607 planesatcornerhits(cornerlist,ii) = planesatcornerhits(cornerlist,ii) + 1;
tp@0 608 end
tp@0 609
tp@0 610 maxplanespercorner = 0;
tp@0 611 for ii = 1:ncorners
tp@0 612 nplanes = length(find(planesatcornerhits(ii,:) ~= 0));
tp@0 613 if nplanes > maxplanespercorner
tp@0 614 maxplanespercorner = nplanes;
tp@0 615 end
tp@0 616 end
tp@0 617
tp@0 618 planesatcorners = zeros(ncorners,maxplanespercorner);
tp@0 619 nplanespercorners = zeros(ncorners,1);
tp@0 620 for ii = 1:ncorners
tp@0 621 iv = find(planesatcornerhits(ii,:)~=0);
tp@0 622 planesatcorners(ii,1:length(iv)) = iv;
tp@0 623 nplanespercorners(ii) = length(iv);
tp@0 624 end
tp@0 625
tp@0 626 % find cubic boxes inside which the planes are placed
tp@0 627
tp@0 628 [nplanes,slask] = size(planeeqs);
tp@0 629
tp@0 630 minvals = zeros(nplanes,3);
tp@0 631 maxvals = zeros(nplanes,3);
tp@0 632
tp@0 633 for ii = 1:nplanes
tp@0 634 cornerlist = planecorners(ii,:);
tp@0 635 cornerlist = cornerlist( find(cornerlist~= 0) );
tp@0 636 cornercoord = corners(cornerlist,:);
tp@0 637 minvals(ii,:) = min(cornercoord);
tp@0 638 maxvals(ii,:) = max(cornercoord);
tp@0 639 end
tp@0 640
tp@0 641 minvals = minvals - geomacc;
tp@0 642 maxvals = maxvals + geomacc;
tp@0 643
tp@0 644 %---------------------------------------------------------------
tp@0 645 % Check the geometry a bit
tp@0 646
tp@0 647 if isempty(checkgeom)
tp@0 648 checkgeom = 0;
tp@0 649 else
tp@0 650 if setstr(checkgeom(1)) == 'c'
tp@0 651 checkgeom = 1;
tp@0 652 else
tp@0 653 checkgeom = 0;
tp@0 654 end
tp@0 655 end
tp@0 656
tp@0 657 if checkgeom
tp@0 658
tp@0 659 badproblem = 0;
tp@0 660 disp(' ')
tp@0 661 disp(' Checking if there are corners with identical coordinates')
tp@0 662 for ii = 1:ncorners-1
tp@0 663 othercorners = [ii+1:ncorners];
tp@0 664 iv = find((corners(othercorners,1)==corners(ii,1)) & (corners(othercorners,2)==corners(ii,2)) & (corners(othercorners,3)==corners(ii,3)));
tp@0 665 if ~isempty(iv)
tp@0 666 disp([' WARNING: corner ',int2str(ii),' and ',int2str(othercorners(iv(1))),' have the same coordinates'])
tp@0 667 disp([' This should be fixed in the CAD-file or in the CATT-model'])
tp@0 668 badproblem = 1;
tp@0 669 end
tp@0 670 end
tp@0 671
tp@0 672 disp(' Checking if there are unused corners')
tp@0 673 iv = find(planesatcorners(:,1)==0);
tp@0 674 if ~isempty(iv)
tp@0 675 disp(' WARNING: The following corners in the CAD-file are not used:')
tp@0 676 printvec = int2str(cornernumbers(iv(1)));
tp@0 677 for iiprint = 2:length(iv)
tp@0 678 printvec = [printvec,' ',int2str(cornernumbers(iv(iiprint)))];
tp@0 679 end
tp@0 680 disp([' ',printvec])
tp@0 681 disp([' This is not a big problem'])
tp@0 682 planesatcorners(iv,:) = [];
tp@0 683 end
tp@0 684
tp@0 685 disp(' Checking if any corners seem redundant')
tp@0 686 iv = find(planesatcorners(:,2)==0);
tp@0 687 if ~isempty(iv)
tp@0 688 disp(' WARNING: The following corners in the CAD-file belong to only one plane:')
tp@0 689 printvec = int2str(cornernumbers(iv(1)));
tp@0 690 for iiprint = 2:length(iv)
tp@0 691 printvec = [printvec,' ',int2str(cornernumbers(iv(iiprint)))];
tp@0 692 end
tp@0 693 disp([' ',printvec])
tp@0 694 disp([' This should be fixed in the CAD-file or in the CATT-model'])
tp@0 695 badproblem = 1;
tp@0 696 end
tp@0 697
tp@0 698 if badproblem == 1
tp@0 699 disp(' ');
tp@0 700 error(['ERROR: Problems with the geometry defined in the CAD-file: ',CADfile])
tp@0 701 end
tp@0 702
tp@0 703 end
tp@0 704
tp@0 705 %---------------------------------------------------------------
tp@0 706 % Save the relevant variables in the output file
tp@0 707
tp@0 708 % NB! Sparse matrices can not get a non-double format
tp@0 709
tp@0 710 if ncorners < 256
tp@0 711 planecorners = uint8(planecorners);
tp@0 712 elseif ncorners < 65536
tp@0 713 planecorners = uint16(planecorners);
tp@0 714 end
tp@0 715
tp@0 716 if max(ncornersperplanevec) <= 255
tp@0 717 ncornersperplanevec = uint8(ncornersperplanevec);
tp@0 718 planehasindents = uint8(planehasindents);
tp@0 719 else
tp@0 720 ncornersperplanevec = uint16(ncornersperplanevec);
tp@0 721 planehasindents = uint16(planehasindents);
tp@0 722 end
tp@0 723
tp@0 724 Varlist = [' corners cornernumbers planecorners planenames planeabstypes '];
tp@0 725 Varlist = [Varlist,' planenumbers cornercrossref cornercrossrefzero Snumbers Snames Sdirectivitynames Sdelays'];
tp@0 726 Varlist = [Varlist,' Scoords Sdirections Srotations Slevels Rnumbers Rcoords CATTversionnumber'];
tp@0 727 Varlist = [Varlist,' planeeqs planenvecs ncornersperplanevec planesatcorners nplanespercorners'];
tp@0 728 Varlist = [Varlist,' minvals maxvals planehasindents indentingcorners'];
tp@0 729
tp@0 730 planenames = sparse(planenames+1-1);
tp@0 731 planeabstypes = sparse(planeabstypes+1-1);
tp@0 732 Sdirectivitynames = sparse(Sdirectivitynames+1-1);
tp@0 733
tp@0 734 eval(['save ',outputfile,Varlist])