tp@0: function outputfile = EDB1edgeox(cadgeofile,outputfile,open_or_closed_model,int_or_ext_model,specorder,difforder,nedgesubs,firstskipcorner,planeseesplanestrategy)
tp@0: % EDB1edgeox - Calculates some plane- and edge-related geometrical parameters.
tp@0: % EDB1edgeo calculates some plane- and edge-related geometrical parameters
tp@0: % based on corners and planes in a .mat-file generated by EDB1readcad
tp@0: % but only data that is independent of the source and receiver.
tp@0: % The output is saved in a .mat-file.
tp@0: %
tp@0: % Input parameters:
tp@0: % cadgeofile (optional) The .mat file with the corners and plane data.
tp@0: % If it is not specified, a file opening window is presented.
tp@0: % outputfile (optional) The .mat file where the output data will be stored.
tp@0: % If it is not specified, an automatic file name is constructed
tp@0: % using the same filename stem as the cadgeofile, but ending
tp@0: % with '_edgeo'.
tp@0: % open_or_closed_model 'o' or 'c' defining whether the model is open or closed.
tp@0: % int_or_ext_model 'i' or 'e' defining whether the model is interior or exterior.
tp@0: % specorder The highest order of any reflection kind (specular and/or diffraction).
tp@0: % difforder (optional) The highest order of diffraction. If it is 0 or 1 then the parameter
tp@0: % edgeseespartialedge is not calculated. Default: 1
tp@0: % nedgesubs (optional) The number of segments that each edge will be
tp@0: % subdivided into for visibility/obstruction tests. Default: 2
tp@0: % NB! When nedgesubs = 2, only the two end points will be checked.
tp@0: % firstskipcorner (optional) All edges including at least one corner with this number or
tp@0: % higher will be excluded from the calculations. Default: 1e6
tp@0: % planeseesplanestrategy (optional) If this parameter is given the value 1, a plane-to-plane
tp@0: % visibility check is done by checking the plane-midpoint to plane-midpoint
tp@0: % for obstructions, which might be enough for some geometries.
tp@0: % SHOWTEXT (global) An integer value; if it is 4 or higher, then messages will be printed on the
tp@0: % screen during calculations.
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: % (Taken directly from the cadgeo-file:)
tp@0: % corners planecorners planeeqs planenvecs ncornersperplanevec
tp@0: % minvals maxvals planehasindents indentingcorners
tp@0: % (Taken directly from the setup-file:)
tp@0: % int_or_ext_model
tp@0: % (Calculated by EDB1edgeo:)
tp@0: % edgecorners Matrix [nedges,2] with the corner numbers that
tp@0: % define each edge, in the reference order.
tp@0: % edgestartcoords Matrix [nedges,3] with the coordinates of the
tp@0: % start corner for each edge.
tp@0: % edgeendcoords Matrix [nedges,3] with the coordinates of the
tp@0: % end corner for each edge.
tp@0: % edgelengthvec List [nedges,1] with the lengths of each edge.
tp@0: % planesatedge Matrix [nedges,2] with the plane numbers of the
tp@0: % two planes that connect to each edge. The first
tp@0: % plane is considered as the reference plane for
tp@0: % each edge. Rule: if the RH thumb is placed
tp@0: % along the edge, from start to end, and the hand
tp@0: % is rotated, the fingers should "come up
tp@0: % through" the reference plane.
tp@0: % closwedangvec List [nedges,1] with the angles in radians of
tp@0: % each edge - for the closed part of each edge.
tp@0: % edgenvecs Matrix [nedges,3] with the normal vectors of
tp@0: % the reference plane for each edge.
tp@0: % offedges List [nlength,1] where nlength = 0-nedges
tp@0: % containing the edge numbers that should not
tp@0: % used.
tp@0: % reflfactors List [nplanes,1] with the reflection factors
tp@0: % of each plane. The only supported values are
tp@0: % +1 for rigid planes, 0 for perfectlt absorbing
tp@0: % planes and -1 for pressure-release planes. The
tp@0: % pressure-release cases might not be implemented
tp@0: % fully.
tp@0: % planeisthin List [nplanes,1] with the value 1 or 0 stating
tp@0: % whether the plane is thin or not.
tp@0: % rearsideplane List [nplanes,1] with the plane number that is at
tp@0: % the rear side. If a plane is non-thin, the
tp@0: % value in this list is zero.
tp@0: % planehasindents List [nplanes,1] with the value 1 or 0 stating
tp@0: % whether a plane has indents or not.
tp@0: % indentingedgepairs Matrix [nindentingedgepairs,2] where each row
tp@0: % gives two connected edges that form an indent
tp@0: % in the shared plane. Consequently, those two
tp@0: % edges can never see each other.
tp@0: % canplaneobstruct List [nplanes,1] with the value 1 or 0 stating
tp@0: % whether a plane has the potential to obstruct
tp@0: % or not.
tp@0: % planeseesplane A matrix, [nplanes,nplanes], with the value 1
tp@0: % 0,-1,-2 stating:
tp@0: % 1 A plane is front of another plane and could
tp@0: % then potentially see it.
tp@0: % 0 A plane is completely behind another plane
tp@0: % but they are not aligned
tp@0: % -1 A plane is aligned with another plane
tp@0: % -2 A plane is back-to-back with another plane.
tp@0: % edgesatplane A matrix, [nplanes,nmaxedgesperplane] which for
tp@0: % row N contains the edge numbers that connect to
tp@0: % to plane N.
tp@0: % edgeseesplane A matrix, [nplanes,nedges], which contains the
tp@0: % values 0,1,-2 or -2, meaning:
tp@0: % 0 The edge is completely behind a plane
tp@0: % -1 The edge can never see the plane and the
tp@0: % edge is aligned with the plane
tp@0: % -2 The edge can never see the plane and the
tp@0: % edge belongs to the plane
tp@0: % 1 The edge has at least one point in front of the plane.
tp@0: %
tp@0: % Uses the functions EDB1coordtrans1
tp@0: % EDB1infrontofplane, EDB1strpend, 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 .
tp@0: % ----------------------------------------------------------------------------------------------
tp@0: % Peter Svensson (svensson@iet.ntnu.no) 20090709
tp@0: %
tp@0: % outputfile = EDB1edgeox(cadgeofile,outputfile,open_or_closed_model,int_or_ext_model,specorder,difforder,nedgesubs,firstskipcorner,planeseesplanestrategy);
tp@0:
tp@0: global SHOWTEXT
tp@0:
tp@0: geomacc = 1e-10;
tp@0:
tp@0: if nargin < 9
tp@0: planeseesplanestrategy = 0;
tp@0: if nargin < 8
tp@0: firstskipcorner = 1e6;
tp@0: if nargin < 7
tp@0: nedgesubs = 2;
tp@0: if nargin < 6
tp@0: difforder = 1;
tp@0: else
tp@0: if isempty(difforder), difforder = 1; end
tp@0: end
tp@0: else
tp@0: if isempty(nedgesubs), nedgesubs = 2; end
tp@0: end
tp@0: else
tp@0: if isempty(firstskipcorner), firstskipcorner = 1e6; end
tp@0: end
tp@0: end
tp@0:
tp@0: %---------------------------------------------------------------
tp@0: % If no cadgeofile was specified, present a file opening window
tp@0:
tp@0: if isempty(cadgeofile)
tp@0: [cadgeofile,filepath] = uigetfile('*.mat','Please select the cadgeofile');
tp@0: [filepath,temp1,temp2] = fileparts(filepath);
tp@0: if ~isstr(cadgeofile)
tp@0: return
tp@0: end
tp@0: [temp1,filestem,fileext,temp2] = fileparts(cadgeofile);
tp@0: cadgeofile = [filepath,filestem,'.mat'];
tp@0: else
tp@0: [filepath,filestem,fileext] = fileparts(cadgeofile);
tp@0: if ~isempty(filepath)
tp@0: cadgeofile = [[filepath,filesep],filestem,'.mat'];
tp@0: else
tp@0: cadgeofile = [filestem,'.mat'];
tp@0: end
tp@0: end
tp@0:
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: filestem = EDB1strpend(filestem,'_cadgeo');
tp@0: outputfile = [[filepath,filesep],filestem,'_eddata.mat'];
tp@0: end
tp@0:
tp@0: %---------------------------------------------------------------
tp@0:
tp@0: eval(['load ',cadgeofile])
tp@0:
tp@0: ncorners = length(cornernumbers);
tp@0: [nplanes,maxncornersperplane] = size(planecorners);
tp@0: minncornersperplane = min(ncornersperplanevec);
tp@0: nedgesperplanevec = double(ncornersperplanevec);
tp@0: ncornersperplanevec = double(ncornersperplanevec);
tp@0:
tp@0: switchoffcorners = find( cornernumbers >= firstskipcorner );
tp@0:
tp@0: onesvec2 = ones(nplanes,1);
tp@0: zerosvec1 = zeros(nplanes,1);
tp@0: zerosvec2 = int8(zeros(nplanes,nplanes));
tp@0:
tp@0: %----------------------------------------------------------------------------
tp@0: %
tp@0: % GEOMETRICAL ACOUSTICS PARAMETERS
tp@0: %
tp@0: %----------------------------------------------------------------------------
tp@0: %
tp@0: % Go through all planeabstypes. If it is:
tp@0: % 'SOFT' or 'soft' then reflfactors should be -1.
tp@0: % 'TOTA' or 'tota' then reflfactors should be 0 (short for TOTABS)
tp@0: % Otherwise it is assumed to be rigid, so reflfactors = 1.
tp@0:
tp@0: if SHOWTEXT >= 4
tp@0: disp(' Check abstypes')
tp@0: end
tp@0:
tp@0: [slask,nchars] = size(planeabstypes);
tp@0: reflfactors = ones(nplanes,1);
tp@0: if nchars >= 4
tp@0: planeabstypes = lower(full(planeabstypes(:,1:min([4,nchars]))));
tp@0:
tp@0: comptxt = 'soft';
tp@0: ivpotential = find(planeabstypes(:,1)==comptxt(1));
tp@0: if ~isempty(ivpotential)
tp@0: comptxt = 'oft';
tp@0: compmat = comptxt(ones(length(ivpotential),1),:);
tp@0: ivsoft = ivpotential(find(prod( (planeabstypes(ivpotential,2:4).*compmat).' ).'));
tp@0: reflfactors(ivsoft) = -1*ones(size(ivsoft));
tp@0: end
tp@0:
tp@0: comptxt = 'tota';
tp@0: ivpotential = find(planeabstypes(:,1)==comptxt(1));
tp@0: if ~isempty(ivpotential)
tp@0: comptxt = 'ota';
tp@0: compmat = comptxt(ones(length(ivpotential),1),:);
tp@0: ivtotabs = ivpotential(find(prod( double(planeabstypes(ivpotential,2:4)==compmat).' ).'));
tp@0: reflfactors(ivtotabs) = zeros(size(ivtotabs));
tp@0: end
tp@0: end
tp@0:
tp@0: %----------------------------------------------------------------------------
tp@0: %
tp@0: % EDGE PARAMETERS
tp@0: %
tp@0: %----------------------------------------------------------------------------
tp@0: %
tp@0: % Check that the planecorners matrix contains no zeros. In such case
tp@0: % add the circular repetition of coordinates.
tp@0:
tp@0: if SHOWTEXT >= 4
tp@0: disp(' Circular')
tp@0: end
tp@0: if sum(sum( planecorners == 0 )) ~= 0
tp@0: for ii = 1:nplanes
tp@0: iv = find( planecorners(ii,:) ~= 0);
tp@0: ncornersatplane = length(iv);
tp@0: if ncornersatplane ~= maxncornersperplane
tp@0: pattern = planecorners(ii,iv);
tp@0: nrepeatings = ceil(maxncornersperplane/ncornersatplane);
tp@0: for jj = 1:nrepeatings-1
tp@0: pattern = [pattern planecorners(ii,iv)];
tp@0: end
tp@0: planecorners(ii,:) = pattern(1:maxncornersperplane);
tp@0: end
tp@0: end
tp@0: end
tp@0:
tp@0: %--------------------------------------------------------------------------------
tp@0: %
tp@0: % Find all edges
tp@0: %
tp@0: % Go through all planes. All consecutive pairs of corners form an edge.
tp@0: % The list planesatedge gives the one or two planes that connect
tp@0: % to each edge. If the second plane number, for a certain edge
tp@0: % is zero, then that specific edge is a plate edge.
tp@0:
tp@0: edgesatplane = zeros( size(planecorners) );
tp@0:
tp@0: if SHOWTEXT >= 4
tp@0: disp(' Defining edges...')
tp@0: end
tp@0:
tp@0: nedgesguess = sum(ncornersperplanevec);
tp@0: edgecorners = zeros(nedgesguess,2);
tp@0: tempplanesatedge = zeros(nedgesguess,1);
tp@0: nedges = 0;
tp@0: thirdpoint = [];
tp@0: edgecounter = 0;
tp@0: multfac = round(10^(floor(log10(nedgesguess))+2));
tp@0:
tp@0: % First go through all planes, and construct edges from the list of corners
tp@0: % that define every plane. This way, many edges will occur twice, unless
tp@0: % a thin plane has no rear side (then the edge occurs a single time)
tp@0: % or if two thin planes are connected (then the edge occurs four times)
tp@0:
tp@0: for ii = 1:nplanes
tp@0: for jj = 1:nedgesperplanevec(ii)
tp@0: edgecounter = edgecounter + 1;
tp@0: corner1 = planecorners(ii,jj);
tp@0: if jj + 1 <= ncornersperplanevec(ii)
tp@0: corner2 = planecorners(ii,jj+1);
tp@0: else
tp@0: corner2 = planecorners(ii,1);
tp@0: end
tp@0: edgecorners(edgecounter,:) = [corner1 corner2];
tp@0: tempplanesatedge(edgecounter) = ii;
tp@0: end
tp@0:
tp@0: end
tp@0:
tp@0: if edgecounter < nedgesguess
tp@0: edgecorners(edgecounter+1:nedgesguess,:) = [];
tp@0: tempplanesatedge(edgecounter+1:nedgesguess,:) = [];
tp@0: end
tp@0:
tp@0: % To find all duplicates, the edge list is sorted in numerical order
tp@0: % If an edge is defined twice or four times, then this is an OK edge.
tp@0: % If an edge is defined a single time, it is connected to a plane which is not connected, and then the
tp@0: % edge should be switched off.
tp@0: % If an edge is defined three times, there must be something like a reflector array where one of
tp@0: % the rear planes is missing.
tp@0:
tp@0: [test,flipvec] = sort(edgecorners.');
tp@0: test = test.';
tp@0: flipvec = flipvec(1,:).';
tp@0:
tp@0: test = test(:,1)*multfac + test(:,2);
tp@0: [test,sortvec] = sort(test);
tp@0: tempplanesatedge = tempplanesatedge(sortvec);
tp@0: flipvec = flipvec(sortvec);
tp@0: planesatedge = [tempplanesatedge.*(flipvec==1) tempplanesatedge.*(flipvec==2)];
tp@0:
tp@0: % Check if there are loners; remove and give a warning.
tp@0:
tp@0: dtest = diff([0;test;test(length(test))+100]);
tp@0: ntest = length(dtest);
tp@0: loners = find(dtest(1:ntest-1)~=0 & dtest(2:ntest)~=0);
tp@0:
tp@0: if ~isempty(loners)
tp@0: disp('WARNING! Some edges had only one plane connected and they will be excluded.')
tp@0: disp(' Check if they should be corrected. ')
tp@0: disp(' Remember that reflectors must have a plane at the rear side too')
tp@0: disp(' Set SHOWTEXT to 4 or higher to get a list of these edges')
tp@0: if SHOWTEXT >= 4
tp@0: disp(' The edges are (EDB1 corner numbers for the two edge point are given):')
tp@0: nloners = length(loners);
tp@0: for jjdisp = 1:nloners
tp@0: indlon = loners(jjdisp);
tp@0: disp([' ',int2str( floor(test(indlon)/multfac) ),' ',int2str( test(indlon)-floor(test(indlon)/multfac)*multfac )])
tp@0: disp([' CATT: ',int2str( cornernumbers(floor(test(indlon)/multfac) )),' ',int2str( cornernumbers(test(indlon)-floor(test(indlon)/multfac)*multfac ))])
tp@0: end
tp@0: end
tp@0:
tp@0: test(loners) = [];
tp@0: planesatedge(loners,:) = [];
tp@0:
tp@0: end
tp@0:
tp@0: ntest = length(test);
tp@0: if ntest >= 2
tp@0: if test(ntest) ~= test(ntest-1)
tp@0: test(ntest) = [];
tp@0: planesatedge(ntest,:) = [];
tp@0: end
tp@0: end
tp@0:
tp@0: % Check if there are triplets
tp@0:
tp@0: triptest = test(1:2:length(test))-test(2:2:length(test));
tp@0: triplets = find(triptest~=0);
tp@0:
tp@0: if ~isempty(triplets)
tp@0: disp('ERROR: Triplets: some edges are defined three times which means that two thin planes have a correct')
tp@0: disp(' definition but some error on the rear side. You must correct these. ')
tp@0: disp(' Only the first one can be detected; it is:')
tp@0: disp(' (EDB1 corner numbers for the two edge points are given)')
tp@0: [floor(test(triplets(1)*2-1)/multfac) test(triplets(1)*2-1)-floor(test(triplets(1)*2-1)/multfac)*multfac]
tp@0: save ~/Documents/Temp/test.mat
tp@0: error
tp@0: end
tp@0:
tp@0: edgecounter = length(test);
tp@0:
tp@0: iv1 = find(planesatedge(:,1)~=0);
tp@0: iv2 = 1:edgecounter;
tp@0: iv2(iv1) = [];
tp@0:
tp@0: edgecorners = [floor(test(iv1)/multfac) test(iv1)-floor(test(iv1)/multfac)*multfac];
tp@0: planesatedge = planesatedge(iv1,:) + planesatedge(iv2,:);
tp@0:
tp@0: [nedges,slask] = size(edgecorners);
tp@0:
tp@0: zerosvec3 = zeros(nedges,1);
tp@0: zerosvec5 = zeros(nedges,3);
tp@0:
tp@0: if SHOWTEXT >= 4
tp@0: disp([' ',int2str(nedges),' edges found'])
tp@0: end
tp@0: thirdpoint = zerosvec5;
tp@0: for ii = 1:nedges
tp@0: refplane = planesatedge(ii,1);
tp@0: secondplane = planesatedge(ii,2);
tp@0: if secondplane ~= 0
tp@0: edgeco1 = corners(edgecorners(ii,1),:);
tp@0: edgeco2 = corners(edgecorners(ii,2),:);
tp@0: edgemidpoint = edgeco1 + (edgeco2-edgeco1)/2;
tp@0: intoplanevec = EDB1cross( (edgeco2-edgeco1).',(planenvecs(secondplane,:)).').';
tp@0: inplanepoint = edgemidpoint + intoplanevec*0.1;
tp@0: if sum(abs(inplanepoint)) == 0
tp@0: inplanepoint = edgemidpoint + intoplanevec*0.2;
tp@0: end
tp@0: thirdpoint(ii,:) = inplanepoint;
tp@0: end
tp@0: end
tp@0:
tp@0:
tp@0: %---------------------------------------------------------------
tp@0: % Calculate the closwedang values for all edges
tp@0:
tp@0: if SHOWTEXT >= 4
tp@0: disp(' wedge angles...')
tp@0: end
tp@0: closwedangvec = zerosvec3;
tp@0:
tp@0: ivec = find( planesatedge(:,1).*planesatedge(:,2) ~= 0);
tp@0: xwedgevec = [corners(edgecorners(ivec,1),:) corners(edgecorners(ivec,2),:)];
tp@0: nvec1vec = planenvecs( planesatedge(ivec,1),: );
tp@0: xsouvec = thirdpoint(ivec,:);
tp@0:
tp@0: for jj = 1:length(ivec)
tp@0:
tp@0: ii = ivec(jj);
tp@0: [slask,thetas,slask] = EDB1coordtrans1(xsouvec(jj,:),[xwedgevec(jj,1:3);xwedgevec(jj,4:6)],nvec1vec(jj,:));
tp@0: if thetas == 0
tp@0: closwedangvec(ii) = 0;
tp@0: else
tp@0: closwedangvec(ii) = 2*pi - thetas;
tp@0: end
tp@0:
tp@0: end
tp@0:
tp@0: %-------------------------------------------------------------------
tp@0: % Now we check for duplicates of the edge definitions
tp@0: % Some edge definitions might need to be swapped: if there are duplicate edges,
tp@0: % and one of them has closwedangvec = 0
tp@0: % then there is a mistake in the edge definition, so a swap will be made.
tp@0:
tp@0: edgecosinglelist = edgecorners(:,1)*multfac + edgecorners(:,2);
tp@0: nsteps = diff([0;edgecosinglelist]);
tp@0: ivduplicates = find(nsteps == 0);
tp@0: nduplicates = length(ivduplicates);
tp@0: if nduplicates > 0
tp@0: if SHOWTEXT >= 4
tp@0: disp([' ',int2str(nduplicates),' edges formed by connected thin plates'])
tp@0: end
tp@0: for ii = 1:nduplicates
tp@0: comb1 = ivduplicates(ii)-1;
tp@0: comb2 = comb1+1;
tp@0: if closwedangvec(comb1) == 0 | closwedangvec(comb2) == 0, % edge definitions should be swapped
tp@0: plane1 = planesatedge(comb1,2);
tp@0: plane2 = planesatedge(comb2,2);
tp@0: planesatedge(comb1,2) = plane2;
tp@0: planesatedge(comb2,2) = plane1;
tp@0: temp1 = thirdpoint(comb1,:);
tp@0: temp2 = thirdpoint(comb2,:);
tp@0: thirdpoint(comb1,:) = temp2;
tp@0: thirdpoint(comb2,:) = temp1;
tp@0: if SHOWTEXT >= 4
tp@0: disp([' Swapping an edge definition'])
tp@0: disp([' Planes ',int2str(planesatedge(comb1,1)),' and ',int2str(planesatedge(comb1,2))])
tp@0: disp([' CATT: ',int2str(planenumbers(planesatedge(comb1,1))),' and ',int2str(planenumbers(planesatedge(comb1,2)))])
tp@0: disp([' Planes ',int2str(planesatedge(comb2,1)),' and ',int2str(planesatedge(comb2,2))])
tp@0: disp([' CATT: ',int2str(planenumbers(planesatedge(comb2,1))),' and ',int2str(planenumbers(planesatedge(comb2,2)))])
tp@0: end
tp@0:
tp@0: xwedge = [corners(edgecorners(comb1,1),:);corners(edgecorners(comb1,2),:)];
tp@0: nvec1 = planenvecs( planesatedge(comb1,1),: );
tp@0: xsou = thirdpoint(comb1,:);
tp@0: [slask,thetas,slask] = EDB1coordtrans1(xsou,xwedge,nvec1);
tp@0: if thetas == 0
tp@0: closwedangvec(comb1) = 0;
tp@0: else
tp@0: closwedangvec(comb1) = 2*pi - thetas;
tp@0: end
tp@0:
tp@0: xwedge = [corners(edgecorners(comb2,1),:);corners(edgecorners(comb2,2),:)];
tp@0: nvec1 = planenvecs( planesatedge(comb2,1),: );
tp@0: xsou = thirdpoint(comb2,:);
tp@0: [slask,thetas,slask] = EDB1coordtrans1(xsou,xwedge,nvec1);
tp@0: if thetas == 0
tp@0: closwedangvec(comb2) = 0;
tp@0: else
tp@0: closwedangvec(comb2) = 2*pi - thetas;
tp@0: end
tp@0:
tp@0: end % if closwedangvec(comb1) == 0 |
tp@0:
tp@0: end % for ii = 1:nduplicates
tp@0: end
tp@0:
tp@0: if nplanes < 256
tp@0: planesatedge = uint8(planesatedge);
tp@0: elseif nplanes < 65535
tp@0: planesatedge = uint16(planesatedge);
tp@0: else
tp@0: planesatedge = uint32(planesatedge);
tp@0: end
tp@0:
tp@0: %-------------------------------------------------------------------
tp@0: % Find which edges each plane is connected to
tp@0:
tp@0: if SHOWTEXT >= 4
tp@0: disp(' find edgesatplane')
tp@0: end
tp@0:
tp@0: edgesatplane = zeros(nplanes,double(max(ncornersperplanevec)));
tp@0:
tp@0: for ii= 1:nplanes
tp@0: iv = find(planesatedge==ii);
tp@0: iv = iv - floor((iv-1)/nedges)*nedges;
tp@0: if ~isempty(iv)
tp@0: edgesatplane(ii,1:length(iv)) = iv.';
tp@0: end
tp@0: end
tp@0:
tp@0: % Check how many edges are defined for each plane. There must be at
tp@0: % least three edges per plane.
tp@0:
tp@0: tempnedgesperplanevec = sum(edgesatplane.'~=0).';
tp@0:
tp@0: iv = find(tempnedgesperplanevec<3);
tp@0:
tp@0: if ~isempty(iv)
tp@0:
tp@0: disp(' ')
tp@0: if SHOWTEXT >= 4
tp@0: for ii = 1:length(iv)
tp@0: disp([' Plane ',int2str(iv(ii)),' has only ',int2str(tempnedgesperplanevec(iv(ii))),' edges connected to it'])
tp@0: end
tp@0: error(['ERROR: The planes above have less than three edges.'])
tp@0: else
tp@0: error(['ERROR: Some planes have less than three edges. Set SHOWTEXT >= 4 to see a list of those planes.'])
tp@0: end
tp@0:
tp@0: end
tp@0:
tp@0:
tp@0:
tp@0: %-------------------------------------------------------------------
tp@0: % Make a special list of thin planes
tp@0:
tp@0: planeisthin = zerosvec1;
tp@0: rearsideplane = zerosvec1;
tp@0:
tp@0: iv = find(closwedangvec==0);
tp@0: if ~isempty(iv)
tp@0: planenumberlist = planesatedge(iv,1);
tp@0: planeisthin(planenumberlist) = planeisthin(planenumberlist) + 1;
tp@0: rearsideplane(planenumberlist) = planesatedge(iv,2);
tp@0: end
tp@0:
tp@0: iv = find(closwedangvec == 0 & planesatedge(:,2)~=0);
tp@0: if ~isempty(iv)
tp@0: planenumberlist = planesatedge(iv,2);
tp@0: planeisthin(planenumberlist) = planeisthin(planenumberlist) + 1;
tp@0: rearsideplane(planenumberlist) = planesatedge(iv,1);
tp@0: end
tp@0:
tp@0: planeisthin = sign(planeisthin);
tp@0: listofthinplanes = find(planeisthin);
tp@0:
tp@0: %---------------------------------------------------------------
tp@0: % Closwedangvec should be calculated into nyvec
tp@0:
tp@0: nyvec = pi./(2*pi - closwedangvec);
tp@0: integerny = ( abs(nyvec - round(nyvec)) < 1e-10 );
tp@0:
tp@0: %-----------------------------------------------------------
tp@0: % Construct some other useful edge variables
tp@0:
tp@0: if difforder >= 1
tp@0: edgestartcoords = zerosvec5;
tp@0: edgestartcoordsnudge = zerosvec5;
tp@0: edgeendcoords = zerosvec5;
tp@0: edgeendcoordsnudge = zerosvec5;
tp@0: edgemidcoords = zerosvec5;
tp@0: edgenvecs = zerosvec5;
tp@0:
tp@0: edgestartcoords = [corners(edgecorners(:,1),:)];
tp@0: edgeendcoords = [corners(edgecorners(:,2),:)];
tp@0: edgemidcoords = (edgestartcoords+edgeendcoords)/2;
tp@0: edgenvecs = planenvecs(planesatedge(:,1),:);
tp@0:
tp@0: edgenormvecs = edgeendcoords - edgestartcoords;
tp@0:
tp@0: edgestartcoordsnudge = edgestartcoords + 1e-10*edgenormvecs;
tp@0: edgeendcoordsnudge = edgeendcoords - 1e-10*edgenormvecs;
tp@0:
tp@0: edgelengthvec = sqrt(sum( ((edgenormvecs).^2).' )).';
tp@0: edgenormvecs = edgenormvecs./edgelengthvec(:,ones(1,3));
tp@0:
tp@0: else
tp@0: edgestartcoords = [];
tp@0: edgeendcoords = [];
tp@0: edgenvecs = [];
tp@0: edgelengthvec = [];
tp@0: end
tp@0:
tp@0:
tp@0: %---------------------------------------------------------------
tp@0: %
tp@0: % BACK TO SOME PURE PLANE RELATED PARAMETERS
tp@0: %
tp@0: %---------------------------------------------------------------
tp@0: %
tp@0: % First, make a list of which planes that are absorbing/non-reflective
tp@0:
tp@0: planeisabsorbing = zerosvec1;
tp@0:
tp@0: listofactiveplanes = find(reflfactors~=0);
tp@0: listofabsorbingplanes = find(reflfactors == 0);
tp@0: planeisabsorbing(listofabsorbingplanes) = ones(size(listofabsorbingplanes));
tp@0:
tp@0: %---------------------------------------------------------------
tp@0: % Help variables: planealignedwplane and planeconnectstoplane
tp@0: %
tp@0: % Preparatory for the check of which planes see which plane: check
tp@0: % if any planes are aligned with each other. If two planes are aligned
tp@0: % then planealignedwplane = 1 for that combination.
tp@0: % However, if two planes are back-to back (thin planes), then
tp@0: % planealignedwplane = 2 for that combination.
tp@0: %
tp@0: % Also make a matrix of which planes connect to which planes and via
tp@0: % which edges.
tp@0:
tp@0: planealignedwplane = zerosvec2;
tp@0:
tp@0: % First, check which planes are aligned with each other
tp@0:
tp@0: for ii = 1:nplanes-1
tp@0: if planeisabsorbing(ii) ~= 1
tp@0: oneplaneeq = planeeqs(ii,:);
tp@0: diffvec1 = oneplaneeq(ones(nplanes-ii,1),:) - planeeqs(ii+1:nplanes,:);
tp@0: diffvec2 = oneplaneeq(ones(nplanes-ii,1),:) + planeeqs(ii+1:nplanes,:);
tp@0: diffvec = min( [sum( diffvec1.'.^2 ).' sum( diffvec2.'.^2 ).'].' ).';
tp@0: alignedplanes = find(diffvec < geomacc) + ii;
tp@0: if ~isempty(alignedplanes)
tp@0: planealignedwplane(alignedplanes,ii) = int8(double(planealignedwplane(alignedplanes,ii)) + double((planeisabsorbing(alignedplanes)~=1)));
tp@0: end
tp@0: end
tp@0: end
tp@0:
tp@0: planealignedwplane = (sparse(sign(double(planealignedwplane) + double(planealignedwplane).')));
tp@0:
tp@0: % Second, check which planes are back to back
tp@0:
tp@0: if ~isempty(listofthinplanes)
tp@0: for ii = 1:length(listofthinplanes)
tp@0: plane = listofthinplanes(ii);
tp@0: rearplane = rearsideplane(plane);
tp@0: if planeisabsorbing(plane) ~= 1 & planeisabsorbing(rearplane) ~= 1
tp@0: planealignedwplane(plane,rearplane) = 2;
tp@0: planealignedwplane(rearplane,plane) = 2;
tp@0: end
tp@0: end
tp@0: end
tp@0:
tp@0: planeconnectstoplane = zerosvec2;
tp@0: clear zerosvec2
tp@0:
tp@0: for ii = 1:nplanes
tp@0: if planeisabsorbing(ii) ~= 1
tp@0: edgelist = edgesatplane(ii,:);
tp@0: edgelist = edgesatplane(ii,1:double(ncornersperplanevec(ii)));
tp@0: ivec = planesatedge(edgelist,:);
tp@0: ivec = reshape(ivec.',length(edgelist)*2,1);
tp@0: ivec = ivec(find(ivec~=ii));
tp@0: okplanes = find(planeisabsorbing(ivec)~=1);
tp@0: ivec = ivec(okplanes);
tp@0: if ~isempty(ivec)
tp@0: planeconnectstoplane(ii,ivec) = edgelist(okplanes);
tp@0: end
tp@0: end
tp@0: end
tp@0:
tp@0: %---------------------------------------------------------------
tp@0: %
tp@0: % planeseesplane
tp@0: %
tp@0: %---------------------------------------------------------------
tp@0: % Check which planes see which planes:
tp@0: % 1. If one of the planes has reflfactors = 0, then the two can't see each other.
tp@0: % 2. Reflective planes that are aligned with each other can not reach each other
tp@0: % 3. Reflective planes that are not aligned but have the same normal vectors
tp@0: % can not reach each other
tp@0: % 4. Planes that have all its corners behind another plane can not be seen by that plane.
tp@0: %
tp@0: % planeseesplane = 0 means that either, one of the planes is non-reflective or, that
tp@0: % two planes never can see each other, but they are not aligned with each other.
tp@0: % planeseesplane = 1 means that two reflective planes might see each other, but there could
tp@0: % be obstructions
tp@0: % planeseesplane = -1 means that two reflective planes are aligned (but they are not back-to-back)
tp@0: % and thus they can never see each other.
tp@0: % planeseesplane = -2 means that two reflective planes are back-to-back with each other and thus
tp@0: % they can never see each other.
tp@0:
tp@0: if SHOWTEXT >= 4
tp@0: disp(' Check which planes see which planes')
tp@0: end
tp@0:
tp@0: % First an auxilliary matrix which can be used for edgeseesplane later:
tp@0: % cornerinfrontofplane, size [nplanes,ncorners]
tp@0: % Values are the same as given by EDB1infrontofplane:
tp@0: % 1 means that a point is in front of the plane
tp@0: % 0 means that a point is aligned with a plane
tp@0: % -1 means that a point is behind a plane
tp@0: % All corners that belong to a plane will have the value 0.
tp@0:
tp@0: % Corner number is given by the col no.
tp@0: if ncorners < 256
tp@0: cornernumb = uint8([1:ncorners]);
tp@0: elseif ncorners < 65536
tp@0: cornernumb = uint16([1:ncorners]);
tp@0: else
tp@0: cornernumb = uint32([1:ncorners]);
tp@0: end
tp@0:
tp@0: % Plane numbers is given by the row no.
tp@0: if nplanes < 256
tp@0: planenumb = uint8([1:nplanes].');
tp@0: elseif nplanes < 65536
tp@0: planenumb = uint16([1:nplanes].');
tp@0: else
tp@0: planenumb = uint32([1:nplanes].');
tp@0: end
tp@0:
tp@0: cornerinfrontofplane = int8(EDB1infrontofplane(corners,planenvecs,planecorners,[],cornernumb,planenumb));
tp@0: clear planenumb cornernumb
tp@0: cornerinfrontofplane = reshape(cornerinfrontofplane,nplanes,ncorners);
tp@0:
tp@0: planeseesplane = int8(ones(nplanes,nplanes));
tp@0:
tp@0: % 1. If one of the planes has reflfactors = 0, then the two can't see each other.
tp@0:
tp@0: if ~isempty(listofabsorbingplanes)
tp@0: zerosvec = zeros(length(listofabsorbingplanes),nplanes);
tp@0: planeseesplane(listofabsorbingplanes,:) = zerosvec;
tp@0: planeseesplane(:,listofabsorbingplanes) = zerosvec.';
tp@0: end
tp@0:
tp@0: % 2. Reflective planes that are aligned with each other can not reach each other
tp@0:
tp@0: ivec = find( planealignedwplane == 1);
tp@0: if ~isempty(ivec)
tp@0: planeseesplane(ivec) = - 1*ones(size(ivec));
tp@0: end
tp@0:
tp@0: ivec = find( planealignedwplane == 2);
tp@0: if ~isempty(ivec)
tp@0: planeseesplane(ivec) = - 2*ones(size(ivec));
tp@0: end
tp@0:
tp@0: % 3. Reflective planes that have the same normal vectors can not reach each other
tp@0:
tp@0: numvec = [1:nplanes];
tp@0: for ii = 1:nplanes
tp@0: if planeisabsorbing(ii) ~= 1
tp@0: nvec1 = planenvecs(ii,:);
tp@0: ivec = find(planealignedwplane(ii,:)==0 & planeisabsorbing.'==0 & numvec>ii);
tp@0: if ~isempty(ivec)
tp@0: similarvec = abs( nvec1(ones(length(ivec),1),:) - planenvecs(ivec,:)) < geomacc;
tp@0: similarvec = prod( double(similarvec.') ).';
tp@0: ivec2 = ivec(find(similarvec==1));
tp@0: if ~isempty(ivec2)
tp@0: zerosvec = zeros(size(ivec2));
tp@0: planeseesplane(ii,ivec2) = zerosvec;
tp@0: planeseesplane(ivec2,ii) = zerosvec.';
tp@0: end
tp@0: end
tp@0: end
tp@0: end
tp@0:
tp@0: % 3.5 A plane does not see itself
tp@0:
tp@0: iv = [0:nplanes-1]*nplanes + [1:nplanes];
tp@0: planeseesplane(iv) = zeros(size(iv));
tp@0:
tp@0: % 4. Planes that have all its corners behind another plane can not be seen by that plane.
tp@0: % First, we construct a list referring to the complete matrix of size [nplanes,nplanes]
tp@0: % The index vector is called iv. Then we find which combinations that
tp@0: % could be cleared. After having cleared a number of combinations (in
tp@0: % iv, fromplanenumb and toplanenumb) we pick out out only the non-zero
tp@0: % indices in iv.
tp@0:
tp@0: if nplanes*nplanes < 128
tp@0: iv = int8([1:nplanes*nplanes].');
tp@0: elseif nplanes*nplanes < 32768
tp@0: iv = int16([1:nplanes*nplanes].');
tp@0: else
tp@0: iv = int32([1:nplanes*nplanes].');
tp@0: end
tp@0: iv = reshape(iv,nplanes,nplanes);
tp@0:
tp@0: % If there are any absorbing planes, we zero all plane-to-plane
tp@0: % combinations that involve such a plane
tp@0:
tp@0: if ~isempty(listofabsorbingplanes)
tp@0: nabsorbingplanes = length(listofabsorbingplanes);
tp@0: iv(listofabsorbingplanes,:) = zeros(nabsorbingplanes,nplanes);
tp@0: iv(:,listofabsorbingplanes) = zeros(nplanes,nabsorbingplanes);
tp@0: end
tp@0:
tp@0: % Check connecting planes first: if two planes are connected and their
tp@0: % shared edge has a closwedangvec > pi, then the two planes can see each
tp@0: % other. So, we can zero the plane-pair combinations where closwedangvec
tp@0: % <= pi but larger than zero.
tp@0:
tp@0: edgecombstozero = find(closwedangvec<=pi & closwedangvec > 0);
tp@0: if ~isempty(edgecombstozero)
tp@0: ivzerocombs = (double(planesatedge(edgecombstozero,1))-1)*nplanes + double(planesatedge(edgecombstozero,2));
tp@0: iv(ivzerocombs) = zeros(size(ivzerocombs));
tp@0: ivzerocombs = (double(planesatedge(edgecombstozero,2))-1)*nplanes + double(planesatedge(edgecombstozero,1));
tp@0: iv(ivzerocombs) = zeros(size(ivzerocombs));
tp@0: end
tp@0:
tp@0: % Now we should keep only the index numbers in iv that are still non-zero
tp@0: % We also take the chance to zero planeseesplane for the combinations that
tp@0: % could be cleared.
tp@0: ivclear = uint32(find(iv==0));
tp@0: planeseesplane(ivclear) = zeros(size(ivclear));
tp@0: iv(ivclear) = [];
tp@0:
tp@0: % Of the remaining plane-to-plane combinations, we pick out the ones
tp@0: % for which we need to check if the "to-plane" is in front of the "from-plane"
tp@0: iv = iv(find(planeseesplane(iv)==1 & planeconnectstoplane(iv)==0));
tp@0:
tp@0: % We create full matrices, fromplanenumb and toplanenumb, and then keep
tp@0: % only the combinations that remain in iv
tp@0: % to keep fromplanenumb and toplanenumb aligned with iv.
tp@0:
tp@0: % To-plane numbers is given by the row no.
tp@0: if nplanes < 256
tp@0: toplanenumb = uint8([1:nplanes].');
tp@0: elseif nplanes < 65536
tp@0: toplanenumb = uint16([1:nplanes].');
tp@0: else
tp@0: toplanenumb = uint32([1:nplanes].');
tp@0: end
tp@0: toplanenumb = reshape(toplanenumb(:,ones(1,nplanes)),nplanes*nplanes,1);
tp@0: toplanenumb = toplanenumb(iv);
tp@0:
tp@0: % From-plane numbers is given by the col no.
tp@0: if nplanes < 256
tp@0: fromplanenumb = uint8([1:nplanes]);
tp@0: elseif nplanes < 65536
tp@0: fromplanenumb = uint16([1:nplanes]);
tp@0: else
tp@0: fromplanenumb = uint32([1:nplanes]);
tp@0: end
tp@0: fromplanenumb = reshape(fromplanenumb(ones(nplanes,1),:),nplanes*nplanes,1);
tp@0: fromplanenumb = fromplanenumb(iv);
tp@0:
tp@0: % Now, if a plane has *all* its corners behind another plane, those two
tp@0: % planes can not see each other.
tp@0: % We check the corners of the "to-plane" (columns) and check if they are in front of
tp@0: % the "from-plane" (rows).
tp@0: % The ivreftocoplamatrix is the index number to the
tp@0: % cornerinfrontofplanematrix, which has the size [nplanes,ncorners].
tp@0: % NB! In order to save some space, we don't construct the ivreftocoplamatrix specifically.
tp@0:
tp@0: % First we check the first 3 corners since all planes have at least 3
tp@0: % corners.
tp@0:
tp@0: % Plane corner 1
tp@0: cornersbehind = cornerinfrontofplane(double(fromplanenumb) + ( double(planecorners( toplanenumb,1 ))-1 )*nplanes)<1;
tp@0:
tp@0: % Plane corner 2
tp@0: cornersbehind = cornersbehind &(cornerinfrontofplane(double(fromplanenumb) + ( double(planecorners( toplanenumb,2 ))-1 )*nplanes) < 1);
tp@0:
tp@0: % Plane corner 3
tp@0: cornersbehind = cornersbehind &(cornerinfrontofplane(double(fromplanenumb) + ( double(planecorners( toplanenumb,3 ))-1 )*nplanes) < 1);
tp@0:
tp@0: if minncornersperplane == 4 & maxncornersperplane == 4
tp@0: % Plane corner 4
tp@0: cornersbehind = cornersbehind &(cornerinfrontofplane(double(fromplanenumb) + ( double(planecorners( toplanenumb,4 ))-1 )*nplanes) < 1);
tp@0: elseif not(minncornersperplane == 3 & maxncornersperplane == 3)
tp@0: ivall = uint32([1:length(fromplanenumb)].');
tp@0: if minncornersperplane == 3
tp@0: iv3cornerplanes = find(ncornersperplanevec(toplanenumb)==3);
tp@0: ivall(iv3cornerplanes) = [];
tp@0: end
tp@0: if maxncornersperplane == 4
tp@0: % Plane corner 4
tp@0: cornersbehind(ivall) = cornersbehind(ivall) &(cornerinfrontofplane(double(fromplanenumb(ivall)) + ( double(planecorners( toplanenumb(ivall),4 ))-1 )*nplanes) < 1);
tp@0: else
tp@0: ivmorethan4 = find(ncornersperplanevec(toplanenumb(ivall))>4);
tp@0: temp = ivall(ivmorethan4);
tp@0: ivall(ivmorethan4) = [];
tp@0: ivmorethan4 = temp;
tp@0: cornersbehind(ivall) = cornersbehind(ivall) &(cornerinfrontofplane(double(fromplanenumb(ivall)) + ( double(planecorners( toplanenumb(ivall),4 ))-1 )*nplanes) < 1);
tp@0: for ii = 5:maxncornersperplane
tp@0: ivselection = find(ncornersperplanevec(toplanenumb(ivmorethan4))>=ii);
tp@0: cornersbehind(ivmorethan4(ivselection)) = cornersbehind(ivmorethan4(ivselection)) &(cornerinfrontofplane(double(fromplanenumb(ivmorethan4(ivselection))) + ( double(planecorners( toplanenumb(ivmorethan4(ivselection)),ii ))-1 )*nplanes) < 1);
tp@0: end
tp@0: clear ivselection
tp@0: end
tp@0: end
tp@0: clear toplanenumb fromplanenumb
tp@0:
tp@0: ivclear = iv(find(cornersbehind==1));
tp@0: clear iv
tp@0:
tp@0: planeseesplane(ivclear) = zeros(size(ivclear));
tp@0: clear ivclear
tp@0:
tp@0: temp1 = (planeseesplane~=0);
tp@0: temp2 = (planeseesplane.'~=0);
tp@0: planeseesplane = int8(double(planeseesplane).*(temp1.*temp2));
tp@0:
tp@0: % New addition:
tp@0: % If the user has asked for it, check plane-mid-point to plane-mid-point
tp@0: % for obstruction. If there are obstructions, set planeseesplane to 0 for
tp@0: % that combination.
tp@0:
tp@0: if planeseesplanestrategy == 1
tp@0: ivorig = uint32(find(planeseesplane==1));
tp@0: iv = ivorig;
tp@0: fromplane = ceil(double(iv)/nplanes);
tp@0: toplane = uint16(double(iv) - (fromplane-1)*nplanes);
tp@0: fromplane = uint16(fromplane);
tp@0: listofplanesthatcanobscure = find( sum(planeseesplane==1) );
tp@0:
tp@0: iv4planes = find(ncornersperplanevec == 4);
tp@0: planemidpoints = zeros(nplanes,3);
tp@0:
tp@0: if any(ncornersperplanevec~=4)
tp@0: disp('WARNING: planeseesplanestrategy 1 not implemented for models with non-4-corner planes')
tp@0: planeseesplanestrategy = 0;
tp@0: else
tp@0: xvalues = corners(planecorners.',1);
tp@0: xvalues = reshape(xvalues,4,length(xvalues)/4);
tp@0: yvalues = corners(planecorners.',2);
tp@0: yvalues = reshape(yvalues,4,length(yvalues)/4);
tp@0: zvalues = corners(planecorners.',3);
tp@0: zvalues = reshape(zvalues,4,length(zvalues)/4);
tp@0: planemidpoints = [mean(xvalues).' mean(yvalues).' mean(zvalues).'];
tp@0: end
tp@0: end
tp@0:
tp@0: if planeseesplanestrategy == 1
tp@0: for ii = listofplanesthatcanobscure
tp@0: planeobsc = ii;
tp@0: ivcheckvis1 = uint32( planeobsc + (double(fromplane)-1)*nplanes);
tp@0: ivcheckvis2 = uint32( planeobsc + (double(toplane)-1)*nplanes);
tp@0: iv3 = find(toplane~=planeobsc & fromplane~=planeobsc & ((planeseesplane(ivcheckvis1) == 1) | (planeseesplane(ivcheckvis2) == 1)) & (planeseesplane(ivcheckvis1) >= 0) & (planeseesplane(ivcheckvis2) >= 0) );
tp@0: tempcanplaneobstruct = zerosvec1.';
tp@0: tempcanplaneobstruct(planeobsc) = 1;
tp@0: [nonobstructedpaths,nobstructions,edgehits,cornerhits] = EDB1checkobstrpaths(planemidpoints(fromplane(iv3),:),planemidpoints(toplane(iv3),:),[],[],tempcanplaneobstruct,planeseesplane,...
tp@0: planeeqs,planenvecs,minvals,maxvals,planecorners,corners,ncornersperplanevec,rearsideplane);
tp@0: if length(iv3) > length(nonobstructedpaths)
tp@0: iv3(nonobstructedpaths) = [];
tp@0: iv(iv3) = [];
tp@0: toplane(iv3) = [];
tp@0: fromplane(iv3) = [];
tp@0: end
tp@0: end
tp@0: planeseesplane(ivorig) = 0;
tp@0: planeseesplane(iv) = 1;
tp@0: clear iv ivorig toplane fromplane ivcheckvis ivcheckvis2
tp@0: end
tp@0:
tp@0: %---------------------------------------------------------------
tp@0: %
tp@0: % EDGE RELATED PARAMETERS AGAIN
tp@0: %
tp@0: %---------------------------------------------------------------
tp@0: %
tp@0: % Which edges should be switched off?
tp@0: % (1) Go through the list edgecorners. If any edge contains one of the
tp@0: % corners, whose numbers are in switchoffcorners, then that edge should be
tp@0: % stored in the list offedges.
tp@0: % (2) If an edge has an integer ny-value, then it should be switched off.
tp@0: % (3) If one of the two connecting planes has reflfactors = 0, then it should
tp@0: % be switched off.
tp@0:
tp@0: % disp(' switching off possibly unwanted edges...')
tp@0:
tp@0: offedges = zerosvec3;
tp@0:
tp@0: % (1) Go through the list edgecorners. If any edge contains one of the
tp@0: % corners, whose numbers are in switchoffcorners, then that edge should be
tp@0: % stored in the list offedges.
tp@0:
tp@0: if ~isempty(switchoffcorners)
tp@0: for ii = 1:nedges
tp@0: corner1 = edgecorners(ii,1);
tp@0: corner2 = edgecorners(ii,2);
tp@0: remove = sum( corner1 == switchoffcorners ) + sum( corner2 == switchoffcorners );
tp@0: offedges(ii) = sign(remove);
tp@0: end
tp@0: end
tp@0:
tp@0: % (2) If an edge has an integer ny-value, then it should be switched off.
tp@0:
tp@0: offedges = offedges + integerny;
tp@0:
tp@0: % (3) If one of the two connecting planes has reflfactors = 0, then it should
tp@0: % be switched off.
tp@0:
tp@0: reflectingedges = reshape(reflfactors(planesatedge),nedges,2);
tp@0: reflectingedges = reflectingedges(:,1).*reflectingedges(:,2);
tp@0: offedges = offedges + (1-reflectingedges);
tp@0:
tp@0: offedges = find(offedges ~= 0);
tp@0:
tp@0: %--------------------------------------------------------------------------------
tp@0: %
tp@0: % edgeseesplane
tp@0: %
tp@0: %--------------------------------------------------------------------------------
tp@0: %
tp@0: % Now check which planes each edge can see.
tp@0: % 1. Set all combinations with reflfactors = 0 to -3
tp@0: % 2. Switch off all combinations with offedges
tp@0: % 3. For all edges that are aligned with a plane, set edgeseesplane to:
tp@0: % -1 if the edge doesn't belong to the plane
tp@0: % -2 if the edge does belong to the plane
tp@0: % 4. If at least one edge end point is in front of a plane, then the
tp@0: % plane is potentially visible (edgeseesplane = 1)
tp@0: % 5. If at least one corner point is in front of:
tp@0: % *one* of the two edge planes, when closwedang < pi
tp@0: % *both* edge planes, when closwedangvec > pi
tp@0:
tp@0: if SHOWTEXT >= 4
tp@0: disp(' Check which edges see which planes')
tp@0: end
tp@0:
tp@0: edgeseesplane = int8(ones(nplanes,nedges));
tp@0:
tp@0: % % % 20120414 PS: Why has this first stage been shut off? CHECK
tp@0: % 1. Set all edges belonging to planes with reflfactors = 0 to
tp@0: % edgeseesplane = -3.
tp@0:
tp@0: % % if ~isempty(listofabsorbingplanes)
tp@0: % % edgeseesplane(listofabsorbingplanes,:) = -3*ones(length(listofabsorbingplanes),nedges);
tp@0: % % end
tp@0:
tp@0: % 2. Switch off all combinations with offedges
tp@0:
tp@0: if ~isempty(offedges)
tp@0: edgeseesplane(:,offedges) = zeros(nplanes,length(offedges));
tp@0: end
tp@0:
tp@0: % 3. For all edges that belong to a plane, set edgeseesplane to -2
tp@0: % Also, for all other edges that are aligned with a plane, set
tp@0: % edgeseesplane to -1
tp@0:
tp@0: edgelist = [1:nedges].';
tp@0: edgelist(offedges) = [];
tp@0:
tp@0: plane1list = planesatedge(edgelist,1);
tp@0: plane2list = planesatedge(edgelist,2);
tp@0: indexvec = uint32( double(plane1list) + (edgelist-1)*nplanes);
tp@0: edgeseesplane(indexvec) = -2*(reflfactors(plane1list)~=0);
tp@0: indexvec = uint32( double(plane2list) + (edgelist-1)*nplanes);
tp@0: edgeseesplane(indexvec) = -2*(reflfactors(plane2list)~=0);
tp@0:
tp@0: for ii = 1:length(edgelist)
tp@0:
tp@0: aligners1 = find(planealignedwplane(:,plane1list(ii))==1);
tp@0: if ~isempty(aligners1)
tp@0: edgeseesplane(aligners1,edgelist(ii)) = -1*ones(size(aligners1));
tp@0: end
tp@0:
tp@0: aligners2 = find(planealignedwplane(:,plane2list(ii))==1);
tp@0: if ~isempty(aligners2)
tp@0: edgeseesplane(aligners2,edgelist(ii)) = -1*ones(size(aligners2));
tp@0: end
tp@0:
tp@0: end
tp@0:
tp@0: % 4. If at least one edge end point is in front of a plane, then the
tp@0: % plane is potentially visible.
tp@0: % Do this for all plane-edge combinations for which edgeseesplane = 1
tp@0: % up til now.
tp@0: % 5. If at least one corner point is in front of:
tp@0: % *one* of the two edge planes, when closwedang < pi
tp@0: % *both* edge planes, when closwedangvec > pi
tp@0:
tp@0: ntot = nplanes*nedges;
tp@0: ivclear = find(edgeseesplane~=1);
tp@0: % Edge number is given by the col no.
tp@0: if nedges < 256
tp@0: edgenumb = uint8([1:nedges]);
tp@0: elseif nedges < 65536
tp@0: edgenumb = uint16([1:nedges]);
tp@0: else
tp@0: edgenumb = uint32([1:nedges]);
tp@0: end
tp@0: edgenumb = reshape(edgenumb(ones(nplanes,1),:),nplanes*nedges,1);
tp@0: edgenumb(ivclear) = [];
tp@0: % Plane numbers is given by the row no.
tp@0: if nplanes < 256
tp@0: planenumb = uint8([1:nplanes].');
tp@0: elseif nplanes < 65536
tp@0: planenumb = uint16([1:nplanes].');
tp@0: else
tp@0: planenumb = uint32([1:nplanes].');
tp@0: end
tp@0: planenumb = reshape(planenumb(:,ones(1,nedges)),nplanes*nedges,1);
tp@0: planenumb(ivclear) = [];
tp@0:
tp@0: if nplanes*nedges < 128
tp@0: iv = int8([1:nplanes*nedges].');
tp@0: elseif nplanes*nedges < 32768
tp@0: iv = int16([1:nplanes*nedges].');
tp@0: else
tp@0: iv = int32([1:nplanes*nedges].');
tp@0: end
tp@0: iv(ivclear) = [];
tp@0:
tp@0: if ~isempty(edgenumb)
tp@0: % For the "edgecorners in front of plane"-test
tp@0: % we can use the cornerinfrontofplane matrix, but we need to construct
tp@0: % index vectors, based on the planenumb and edgenumb values.
tp@0:
tp@0: edgeinfrontofplane = cornerinfrontofplane( double(planenumb) + ( double( edgecorners(edgenumb,1) )-1 )*nplanes );
tp@0: edgeinfrontofplane = (edgeinfrontofplane==1) | (cornerinfrontofplane( double(planenumb) + ( double( edgecorners(edgenumb,2) )-1 )*nplanes )==1);
tp@0: edgeseesplane(iv) = int8(double(edgeseesplane(iv)).*double(edgeinfrontofplane));
tp@0: clear edgeinfrontofplane
tp@0: end
tp@0:
tp@0: if ~isempty(edgenumb)
tp@0: % For the "planecorners in front of edge-plane"-test
tp@0: % we can use the cornerinfrontofplane matrix, but we need to construct
tp@0: % index vectors, based on the planenumb and edgenumb values.
tp@0: %
tp@0: % First we split up the iv into two halves: the one with edges < pi and the
tp@0: % one with edges > pi;
tp@0:
tp@0: ivsmaller = uint32(find(closwedangvec(edgenumb)4);
tp@0: temp = ivall(ivmorethan4);
tp@0: ivall(ivmorethan4) = [];
tp@0: ivmorethan4 = temp;
tp@0: % Plane corner 4
tp@0: ivreftocoplamatrix = double( planesatedge(edgenumb(ivsmaller(ivall)),1) ) + ( double(planecorners( planenumb(ivsmaller(ivall)),4 ))-1 )*nplanes;
tp@0: planecornerinfrontofedge(ivall) = planecornerinfrontofedge(ivall) | (cornerinfrontofplane(ivreftocoplamatrix)==1);
tp@0: ivreftocoplamatrix = double( planesatedge(edgenumb(ivsmaller(ivall)),2) ) + ( double(planecorners( planenumb(ivsmaller(ivall)),4 ))-1 )*nplanes;
tp@0: planecornerinfrontofedge(ivall) = planecornerinfrontofedge(ivall) | (cornerinfrontofplane(ivreftocoplamatrix)==1);
tp@0: for ii = 5:maxncornersperplane
tp@0: % Plane corner 5,6,7,etc
tp@0: ivselection = find(ncornersperplanevec(planenumb(ivsmaller(ivmorethan4)))>=ii);
tp@0: ivreftocoplamatrix = double( planesatedge(edgenumb(ivsmaller(ivmorethan4)),1) ) + ( double(planecorners( planenumb(ivsmaller(ivmorethan4)),4 ))-1 )*nplanes;
tp@0: planecornerinfrontofedge(ivmorethan4) = planecornerinfrontofedge(ivmorethan4) | (cornerinfrontofplane(ivreftocoplamatrix)==1);
tp@0: ivreftocoplamatrix = double( planesatedge(edgenumb(ivsmaller(ivmorethan4)),2) ) + ( double(planecorners( planenumb(ivsmaller(ivmorethan4)),4 ))-1 )*nplanes;
tp@0: planecornerinfrontofedge(ivmorethan4) = planecornerinfrontofedge(ivmorethan4) | (cornerinfrontofplane(ivreftocoplamatrix)==1);
tp@0: end
tp@0: clear ivselection
tp@0: end
tp@0: end
tp@0: clear ivreftocoplamatrix
tp@0:
tp@0: ivclear = ivsmaller(find(planecornerinfrontofedge==0));
tp@0: edgeseesplane(iv(ivclear)) = 0;
tp@0: end
tp@0:
tp@0: if ~isempty(ivlarger)
tp@0: % Second the edges larger than pi (closwedang > pi): at least one plane corner needs
tp@0: % to be in front of *both* edge planes.
tp@0: % Plane corner 1
tp@0: ivreftocoplamatrix = double( planesatedge(edgenumb(ivlarger),1) ) + ( double(planecorners( planenumb(ivlarger),1 ))-1 )*nplanes;
tp@0: planecornerinfrontofedge = cornerinfrontofplane(ivreftocoplamatrix)==1;
tp@0: ivreftocoplamatrix = double( planesatedge(edgenumb(ivlarger),2) ) + ( double(planecorners( planenumb(ivlarger),1 ))-1 )*nplanes;
tp@0: planecornerinfrontofedge = planecornerinfrontofedge & (cornerinfrontofplane(ivreftocoplamatrix)==1);
tp@0: % Plane corner 2
tp@0: ivreftocoplamatrix = double( planesatedge(edgenumb(ivlarger),1) ) + ( double(planecorners( planenumb(ivlarger),2 ))-1 )*nplanes;
tp@0: temp = cornerinfrontofplane(ivreftocoplamatrix)==1;
tp@0: ivreftocoplamatrix = double( planesatedge(edgenumb(ivlarger),2) ) + ( double(planecorners( planenumb(ivlarger),2 ))-1 )*nplanes;
tp@0: planecornerinfrontofedge = planecornerinfrontofedge | (temp & (cornerinfrontofplane(ivreftocoplamatrix)==1));
tp@0: % Plane corner 3
tp@0: ivreftocoplamatrix = double( planesatedge(edgenumb(ivlarger),1) ) + ( double(planecorners( planenumb(ivlarger),3 ))-1 )*nplanes;
tp@0: temp = cornerinfrontofplane(ivreftocoplamatrix)==1;
tp@0: ivreftocoplamatrix = double( planesatedge(edgenumb(ivlarger),2) ) + ( double(planecorners( planenumb(ivlarger),3 ))-1 )*nplanes;
tp@0: planecornerinfrontofedge = planecornerinfrontofedge | (temp & (cornerinfrontofplane(ivreftocoplamatrix)==1));
tp@0:
tp@0: if minncornersperplane == 4 & maxncornersperplane == 4
tp@0: % Plane corner 4
tp@0: ivreftocoplamatrix = double( planesatedge(edgenumb(ivlarger),1) ) + ( double(planecorners( planenumb(ivlarger),4 ))-1 )*nplanes;
tp@0: temp = cornerinfrontofplane(ivreftocoplamatrix)==1;
tp@0: ivreftocoplamatrix = double( planesatedge(edgenumb(ivlarger),2) ) + ( double(planecorners( planenumb(ivlarger),4 ))-1 )*nplanes;
tp@0: planecornerinfrontofedge = planecornerinfrontofedge | (temp & (cornerinfrontofplane(ivreftocoplamatrix)==1));
tp@0: elseif not(minncornersperplane == 3 & maxncornersperplane == 3)
tp@0: ivall = uint32(1:length(ivlarger));
tp@0: if minncornersperplane == 3
tp@0: iv3cornerplanes = find(ncornersperplanevec(planenumb(ivlarger))==3);
tp@0: ivall(iv3cornerplanes) = [];
tp@0: end
tp@0: if maxncornersperplane == 4
tp@0: % Plane corner 4
tp@0: ivreftocoplamatrix = double( planesatedge(edgenumb(ivlarger(ivall)),1) ) + ( double(planecorners( planenumb(ivlarger(ivall)),4 ))-1 )*nplanes;
tp@0: temp = cornerinfrontofplane(ivreftocoplamatrix)==1;
tp@0: ivreftocoplamatrix = double( planesatedge(edgenumb(ivlarger(ivall)),2) ) + ( double(planecorners( planenumb(ivlarger(ivall)),4 ))-1 )*nplanes;
tp@0: planecornerinfrontofedge(ivall) = planecornerinfrontofedge(ivall) | (temp & (cornerinfrontofplane(ivreftocoplamatrix)==1));
tp@0: else
tp@0: ivmorethan4 = find(ncornersperplanevec(planenumb(ivlarger(ivall)))>4);
tp@0: temp = ivall(ivmorethan4);
tp@0: ivall(ivmorethan4) = [];
tp@0: ivmorethan4 = temp;
tp@0: % Plane corner 4
tp@0: ivreftocoplamatrix = double( planesatedge(edgenumb(ivlarger(ivall)),1) ) + ( double(planecorners( planenumb(ivlarger(ivall)),4 ))-1 )*nplanes;
tp@0: temp = cornerinfrontofplane(ivreftocoplamatrix)==1;
tp@0: ivreftocoplamatrix = double( planesatedge(edgenumb(ivlarger(ivall)),2) ) + ( double(planecorners( planenumb(ivlarger(ivall)),4 ))-1 )*nplanes;
tp@0: planecornerinfrontofedge(ivall) = planecornerinfrontofedge(ivall) | (temp & (cornerinfrontofplane(ivreftocoplamatrix)==1));
tp@0: for ii = 5:maxncornersperplane
tp@0: % Plane corner 5,6,7,etc
tp@0: ivselection = find(ncornersperplanevec(planenumb(ivlarger(ivmorethan4)))>=ii);
tp@0: ivreftocoplamatrix = double( planesatedge(edgenumb(ivlarger(ivmorethan4)),1) ) + ( double(planecorners( planenumb(ivlarger(ivmorethan4)),4 ))-1 )*nplanes;
tp@0: temp = cornerinfrontofplane(ivreftocoplamatrix)==1;
tp@0: ivreftocoplamatrix = double( planesatedge(edgenumb(ivlarger(ivmorethan4)),2) ) + ( double(planecorners( planenumb(ivlarger(ivmorethan4)),4 ))-1 )*nplanes;
tp@0: planecornerinfrontofedge(ivmorethan4) = planecornerinfrontofedge(ivmorethan4) | (temp & (cornerinfrontofplane(ivreftocoplamatrix)==1));
tp@0: end
tp@0: clear ivselection
tp@0: end
tp@0: end
tp@0:
tp@0: ivclear = ivlarger(find(planecornerinfrontofedge==0));
tp@0: edgeseesplane(iv(ivclear)) = 0;
tp@0: end
tp@0: end
tp@0:
tp@0: clear ivlarger ivsmaller ivreftocoplamatrix edgenumb planenumb iv
tp@0:
tp@0: %----------------------------------------------------------------------------
tp@0: %
tp@0: % indentingedgepairs
tp@0: %
tp@0: %--------------------------------------------------------------------------------
tp@0: %
tp@0: % If there are any planes with indents, make a list of all consecutive edge pairs that
tp@0: % form an indent, because such edge pairs could never see each other.
tp@0:
tp@0: indentingedgepairs = [];
tp@0:
tp@0: if any(planehasindents)
tp@0: iv = find(indentingcorners);
tp@0: indentingedgepairs = zeros(length(iv),2);
tp@0: indentingplanenumbers = mod(iv,nplanes);
tp@0: ivzeros = find(indentingplanenumbers==0);
tp@0: if ~isempty(ivzeros)
tp@0: indentingplanenumbers(ivzeros) = nplanes;
tp@0: end
tp@0: conumbers = (iv - indentingplanenumbers)/nplanes+1;
tp@0:
tp@0: % We expand planecorners temporarily, cyclically
tp@0: planecorners = [planecorners planecorners(:,1:2)];
tp@0:
tp@0: for ii = 1:length(iv)
tp@0: edge1cornerpair = sort(planecorners(indentingplanenumbers(ii),conumbers(ii):conumbers(ii)+1));
tp@0: edge2cornerpair = sort(planecorners(indentingplanenumbers(ii),conumbers(ii)+1:conumbers(ii)+2));
tp@0:
tp@0: [slask,edge1] = ismember(edge1cornerpair,edgecorners,'rows');
tp@0: [slask,edge2] = ismember(edge2cornerpair,edgecorners,'rows');
tp@0:
tp@0: edgepair = sort([edge1 edge2]);
tp@0: indentingedgepairs(ii,:) = edgepair;
tp@0:
tp@0: end
tp@0:
tp@0: % Remove the temporary expansion
tp@0:
tp@0: planecorners = planecorners(:,1:size(planecorners,2)-2);
tp@0: end
tp@0:
tp@0: %----------------------------------------------------------------------------
tp@0: %
tp@0: % canplaneobstruct
tp@0: %
tp@0: %--------------------------------------------------------------------------------
tp@0: %
tp@0: % Check which planes that are potentially obstructing:
tp@0: %
tp@0: % If **all** closwedangvec are < pi and we have an exterior problem =>
tp@0: % the scatterer is without indents.
tp@0: %
tp@0: % If **all** closwedangvec are > pi and we have an interior problem =>
tp@0: % the room is without any indents so there can be no obstructions.
tp@0: %
tp@0: % For exterior problems, all active planes can potentially obstruct.
tp@0: %
tp@0: % If these simple checks give a results that conflict with the input
tp@0: % parameters open_or_closed_model or int_or_ext_model, give an
tp@0: % error message.
tp@0:
tp@0: canplaneobstruct = [];
tp@0: if sum(closwedangvec= 2
tp@0: disp(' ')
tp@0: disp(' The model is an interior problem, and the room has no indents')
tp@0: disp(' so no obstruction tests will be necessary.')
tp@0: end
tp@0: if int_or_ext_model == 'e'
tp@0: disp( 'ERROR: In the setup file this model was defined as an exterior problem, but it is clearly an interior problem.')
tp@0: error(' Please check the orientation of the planes.')
tp@0: end
tp@0: canplaneobstruct = zeros(1,nplanes);
tp@0: elseif sum(closwedangvec>0) == 0
tp@0: if SHOWTEXT >= 2
tp@0: disp(' ')
tp@0: disp(' The model is an exterior problem (scattering)')
tp@0: disp(' with only thin planes.')
tp@0: end
tp@0: if int_or_ext_model == 'i'
tp@0: disp( 'ERROR: In the setup file this model was defined as an interior problem, but it is clearly an exterior problem.')
tp@0: error(' Please check the orientation of the planes.')
tp@0: end
tp@0: canplaneobstruct = ones(1,nplanes);
tp@0: elseif sum(closwedangvec>pi) == 0
tp@0: if SHOWTEXT >= 2
tp@0: disp(' ')
tp@0: disp(' The model is an exterior problem (scattering)')
tp@0: disp(' and the scatterer has no indents.')
tp@0: end
tp@0: if int_or_ext_model == 'i'
tp@0: disp( 'ERROR: In the setup file this model was defined as an interior problem, but it is clearly an exterior problem.')
tp@0: error(' Please check the orientation of the planes.')
tp@0: end
tp@0: canplaneobstruct = ones(1,nplanes);
tp@0: end
tp@0:
tp@0: % For an interior problem we know for sure that if a plane
tp@0: % has all other points in front of itself, it can not obstruct.
tp@0: % NB! We need to check only points that do not belong to planes
tp@0: % that are aligned with the plane, or planes that are non-active.
tp@0:
tp@0: if int_or_ext_model == 'e'
tp@0: canplaneobstruct = (reflfactors.'~=0);
tp@0: elseif isempty(canplaneobstruct)
tp@0:
tp@0: canplaneobstruct = (reflfactors.'~=0);
tp@0: maskvec1 = canplaneobstruct.';
tp@0: listofactiveplanes = find(canplaneobstruct);
tp@0: zerosvec = zeros(nplanes,1);
tp@0:
tp@0: for ii = 1:length(listofactiveplanes)
tp@0: iplane = listofactiveplanes(ii);
tp@0: maskvec2 = zerosvec;
tp@0: maskvec2(iplane) = 1;
tp@0: otherplanestocheck = find((double(planeseesplane(:,iplane))+maskvec2)~=1 & planeseesplane(:,iplane)~=-1 & maskvec1~=0);
tp@0: if ~isempty(otherplanestocheck),
tp@0: cornerstocheck = unique(planecorners(otherplanestocheck,:));
tp@0: cornerstocheck = setdiff(unique(cornerstocheck),planecorners(iplane,1:ncornersperplanevec(iplane)));
tp@0: pointinfront = EDB1infrontofplane(corners(cornerstocheck,:),planenvecs(iplane,:),corners(planecorners(iplane,1),:),corners(planecorners(iplane,2),:));
tp@0: if isempty(find(pointinfront==-1))
tp@0: canplaneobstruct(iplane) = 0;
tp@0: end
tp@0: else
tp@0: canplaneobstruct(iplane) = 0;
tp@0: end
tp@0:
tp@0: end
tp@0: end
tp@0:
tp@0: %---------------------------------------------------------------
tp@0: % For each thin plane, make sure that all edges have the same normal
tp@0: % vector. If not, switch the direction of that edge, and all the other
tp@0: % relevant parameters.
tp@0:
tp@0: if difforder >= 1
tp@0: ivthin = find(planeisthin);
tp@0:
tp@0: for ii = 1:length(ivthin)
tp@0: plane = ivthin(ii);
tp@0: if rearsideplane(plane) > plane
tp@0: edgelist = unique(edgesatplane(plane,1:nedgesperplanevec(plane)));
tp@0: iv = find(closwedangvec(edgelist,:)==0);
tp@0: if ~isempty(iv)
tp@0: edgelist = edgelist(iv);
tp@0: nedgestemp = length(edgelist);
tp@0: edgenveclist = edgenvecs(edgelist,:);
tp@0: refnvec = edgenveclist(1,:);
tp@0: nvecdiff = edgenveclist(2:nedgestemp,:) - refnvec(ones(nedgestemp-1,1),:);
tp@0: nvecdiff = sum(abs(nvecdiff.')).';
tp@0: ivswitch = find(nvecdiff)+1;
tp@0: if ~isempty(ivswitch)
tp@0: edgenumber = edgelist(ivswitch);
tp@0: edgecorners(edgenumber,:) = [edgecorners(edgenumber,2) edgecorners(edgenumber,1)];
tp@0: planesatedge(edgenumber,:) = [planesatedge(edgenumber,2) planesatedge(edgenumber,1)];
tp@0: edgenvecs(edgenumber,:) = -edgenvecs(edgenumber,:);
tp@0: edgestartcoords(edgenumber,:) = corners(edgecorners(edgenumber,1),:);
tp@0: edgeendcoords(edgenumber,:) = corners(edgecorners(edgenumber,2),:);
tp@0: tempvec = edgeendcoordsnudge(edgenumber,:);
tp@0: edgeendcoordsnudge(edgenumber,:) = edgestartcoordsnudge(edgenumber,:);
tp@0: edgestartcoordsnudge(edgenumber,:) = tempvec;
tp@0: end
tp@0: end
tp@0: end
tp@0:
tp@0: end
tp@0: end
tp@0:
tp@0: %----------------------------------------------------------------------------
tp@0: %
tp@0: % SAVE THE VARIABLES
tp@0: %
tp@0: %----------------------------------------------------------------------------
tp@0: %
tp@0: % Also the variables from cadgeo??
tp@0: % Yes: planeeqs planenvecs ncornersperplanevec planecorners corners
tp@0: % minvals maxvals
tp@0:
tp@0: edgeseesplane = int8(edgeseesplane);
tp@0: planeseesplane = int8(planeseesplane);
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: planeisthin = uint8(planeisthin);
tp@0:
tp@0: if max(ncornersperplanevec) <= 255
tp@0: ncornersperplanevec = uint8(ncornersperplanevec);
tp@0: else
tp@0: ncornersperplanevec = uint16(ncornersperplanevec);
tp@0: end
tp@0:
tp@0: Varlist = [ ' edgecorners planesatedge closwedangvec planehasindents indentingedgepairs'];
tp@0: Varlist = [Varlist,' planeisthin planeseesplane rearsideplane edgeseesplane canplaneobstruct'];
tp@0: Varlist = [Varlist,' corners planecorners planeeqs planenvecs ncornersperplanevec cornerinfrontofplane'];
tp@0: Varlist = [Varlist,' edgestartcoords edgeendcoords edgenvecs reflfactors edgesatplane edgelengthvec'];
tp@0: Varlist = [Varlist,' minvals maxvals offedges '];
tp@0: Varlist = [Varlist,' int_or_ext_model'];
tp@0: if difforder >= 1
tp@0: edgedata.edgestartcoordsnudge = edgestartcoordsnudge;
tp@0: edgedata.edgeendcoordsnudge = edgeendcoordsnudge;
tp@0: edgedata.edgenormvecs = edgenormvecs;
tp@0: Varlist = [Varlist, ' edgedata '];
tp@0: end
tp@0: eval(['save ',outputfile,Varlist])