level_model_slice

PURPOSE ^

LEVEL_MODEL_SLICE - level 3D points for slicing at z=0

SYNOPSIS ^

function [out, out2, out3, out4] = level_model_slice(varargin)

DESCRIPTION ^

LEVEL_MODEL_SLICE - level 3D points for slicing at z=0
 
 XYZ = LEVEL_MODEL_SLICE(NODES, LEVEL)
 XYZ = LEVEL_MODEL_SLICE(NODES, LEVEL, SN)
 Transforms a list of 3D points such that level is a slice at z=0.
    NODES: n_nodes x 3 list of point coordinates 
    LEVEL: definition of the slice(s). Can be:
     - scalar: number of slices along the z axis ordered max to min
     - vector (N x 3): intercepts with the coordinate axis
     - struct:
       .centre (N x 3): centre point for transformation
        .normal_angle (1 x 4) [nx, ny, nz, theta]: normal vector and angle 
          of rotation around the axis defined by it (right hand rule).
         The model will be rotated such that the normal vector becomes 
         [0,0,1].         
        .rotation_matrix (3 x 3)
        An error is thrown if both normal_angle and rotation_matrix are present.
   SN     : slice number to use if level defines multiple slices (default: 1).
    XYZ     : NODES transformed such that level is a slice at z=0.

 MDL = LEVEL_MODEL_SLICE(MDL,__) 
 Modifies the nodes of an EIDORS fwd_model structure MDL. Can be used in
 place of NODES with every other syntax.

 {XYZ} = LEVEL_MODEL_SLICE(MDL, LEVEL, 'all') 
 {XYZ} = LEVEL_MODEL_SLICE(NODES, LEVEL, 'all') 
 Returns a cell array (one per slice) of transformed NODES (or MDL.nodes).

 N = LEVEL_MODEL_SLICE(LEVEL) returns the number of slices defined in LEVEL.

 [__, C, M] = LEVEL_MODEL_SLICE(NODES/MDL, __) returns the center C and 
 tranformation matrix M.

 [__, C, M, ifun] = LEVEL_MODEL_SLICE(NODES/MDL,__) also returns a function 
 for an inverse transformation.

 [C,M] = LEVEL_MODEL_SLICE(LEVEL, SN) (discouraged) returns the center C and
 transformation matrix M of slice SN (optional, default 1). An error is 
 thrown if level is scalar. Note that two outputs must always be requested for 
 this syntax.  

 [C,M] = LEVEL_MODEL_SLICE(NODES, LEVEL, SN) (discouraged) allows to obtain C
 and M when LEVEL is scalar. Note that two outputs must always be requested for
 this syntax.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [out, out2, out3, out4] = level_model_slice(varargin)
0002 %LEVEL_MODEL_SLICE - level 3D points for slicing at z=0
0003 %
0004 % XYZ = LEVEL_MODEL_SLICE(NODES, LEVEL)
0005 % XYZ = LEVEL_MODEL_SLICE(NODES, LEVEL, SN)
0006 % Transforms a list of 3D points such that level is a slice at z=0.
0007 %    NODES: n_nodes x 3 list of point coordinates
0008 %    LEVEL: definition of the slice(s). Can be:
0009 %     - scalar: number of slices along the z axis ordered max to min
0010 %     - vector (N x 3): intercepts with the coordinate axis
0011 %     - struct:
0012 %       .centre (N x 3): centre point for transformation
0013 %        .normal_angle (1 x 4) [nx, ny, nz, theta]: normal vector and angle
0014 %          of rotation around the axis defined by it (right hand rule).
0015 %         The model will be rotated such that the normal vector becomes
0016 %         [0,0,1].
0017 %        .rotation_matrix (3 x 3)
0018 %        An error is thrown if both normal_angle and rotation_matrix are present.
0019 %   SN     : slice number to use if level defines multiple slices (default: 1).
0020 %    XYZ     : NODES transformed such that level is a slice at z=0.
0021 %
0022 % MDL = LEVEL_MODEL_SLICE(MDL,__)
0023 % Modifies the nodes of an EIDORS fwd_model structure MDL. Can be used in
0024 % place of NODES with every other syntax.
0025 %
0026 % {XYZ} = LEVEL_MODEL_SLICE(MDL, LEVEL, 'all')
0027 % {XYZ} = LEVEL_MODEL_SLICE(NODES, LEVEL, 'all')
0028 % Returns a cell array (one per slice) of transformed NODES (or MDL.nodes).
0029 %
0030 % N = LEVEL_MODEL_SLICE(LEVEL) returns the number of slices defined in LEVEL.
0031 %
0032 % [__, C, M] = LEVEL_MODEL_SLICE(NODES/MDL, __) returns the center C and
0033 % tranformation matrix M.
0034 %
0035 % [__, C, M, ifun] = LEVEL_MODEL_SLICE(NODES/MDL,__) also returns a function
0036 % for an inverse transformation.
0037 %
0038 % [C,M] = LEVEL_MODEL_SLICE(LEVEL, SN) (discouraged) returns the center C and
0039 % transformation matrix M of slice SN (optional, default 1). An error is
0040 % thrown if level is scalar. Note that two outputs must always be requested for
0041 % this syntax.
0042 %
0043 % [C,M] = LEVEL_MODEL_SLICE(NODES, LEVEL, SN) (discouraged) allows to obtain C
0044 % and M when LEVEL is scalar. Note that two outputs must always be requested for
0045 % this syntax.
0046 %
0047 
0048 
0049 % (C) 2006-2022 Andy Adler and Bartek Grychtol
0050 % License: GPL version 2 or version 3
0051 % $Id: level_model_slice.m 7045 2024-11-30 21:10:50Z bgrychtol $
0052 
0053 if nargin==1 && ischar(varargin{1}) && strcmp(varargin{1}, 'UNIT_TEST')
0054   do_unit_test; return
0055 end
0056 
0057 if nargin == 0 || isempty(varargin{1})
0058    error('EIDORS:WrongInput', 'level_model_slice requires non-empty input');
0059 end
0060 
0061 if nargout <= 1 || nargout >= 3
0062     if nargin == 1
0063         % N = LEVEL_MODEL_SLICE(LEVEL) returns the number of slices defined in LEVEL
0064         out = get_number_of_slices(varargin{1});
0065     else
0066         % XYZ = LEVEL_MODEL_SLICE(NODES, LEVEL)
0067         % XYZ = LEVEL_MODEL_SLICE(NODES, LEVEL, SN)
0068         if nargin == 3 && ischar(varargin{3}) && strcmp(varargin{3},'all')
0069            n_slices = get_number_of_slices(varargin{2});
0070            if isstruct(varargin{1}), varargin{1}=varargin{1}.nodes; end
0071            for i = n_slices:-1:1
0072               out(i,1) = cell(1);
0073               [out{i,1}, C(i,:), M] = level_model_slice(varargin{1},varargin{2}, i);
0074            end
0075         else % nargin==3
0076             [C, M] = matrix_from_level(varargin{:});
0077             SHIFT = [C(1), C(2), C(3)] * M';
0078             SHIFT(3) = 0;
0079             % accept fwd_model as input
0080             if isstruct(varargin{1})
0081                out = varargin{1};
0082                % (nodes - C) * M' + SHIFT
0083                out.nodes = bsxfun(@plus, bsxfun(@minus, out.nodes, C ) * M', SHIFT);
0084             else
0085                % (nodes - C) * M' + SHIFT
0086                out = bsxfun(@plus, bsxfun(@minus, varargin{1}, C ) * M', SHIFT);
0087             end
0088             if nargout >= 3
0089                 out2 = C;
0090                 out3 = M;
0091             end
0092             if nargout == 4
0093                out4 = @(nodes) (nodes - SHIFT) / M' + C;
0094             end
0095         end
0096     end
0097 elseif nargout == 2
0098     if nargin < 3
0099         % [C,M] = LEVEL_MODEL_SLICE(LEVEL)
0100         % [C,M] = LEVEL_MODEL_SLICE(LEVEL, SN)
0101         [out, out2] = matrix_from_level(varargin{:});
0102     else
0103         % [C,M] = LEVEL_MODEL_SLICE(NODES, LEVEL, SN)
0104         [out, out2] = matrix_from_number(varargin{:});
0105     end
0106 end
0107 
0108 function N = get_number_of_slices(level)  
0109     if isnumeric(level) && isscalar(level)
0110         check_scalar_level(level);
0111       N = level;
0112     elseif isnumeric(level)
0113         N = size(level,1);
0114     elseif isstruct(level) && isscalar(level)
0115         N = size(level.centre, 1);
0116     else
0117         error('level must be numeric or a scalar struct');
0118     end
0119 
0120 function [C,M] = matrix_from_level(varargin)
0121     nodes = []; sn = 1;
0122     switch nargin
0123         case 3 % nodes, level, sn
0124             nodes = varargin{1}; level = varargin{2}; sn = varargin{3};
0125         case 2 % level & sn or nodes & level
0126             if isnumeric(varargin{2}) && isscalar(varargin{2}) 
0127                 level = varargin{1}; sn = varargin{2};
0128             else
0129                 nodes = varargin{1}; level = varargin{2};
0130             end
0131         case 1 % level
0132             level = varargin{1};
0133     end
0134 
0135     if isnumeric(level) && isscalar(level)
0136         if isempty(nodes)
0137             error('Need nodes to calculate with scalar numeric level');
0138         else
0139             [C, M] = matrix_from_number(nodes, level, sn);
0140         end
0141     elseif isnumeric(level)
0142         [C, M] = matrix_from_intercept(level, sn);
0143     elseif isstruct(level) && isscalar(level)
0144         [C, M] = matrix_from_struct(level, sn);
0145     else
0146         error('level must be numeric or a scalar struct');
0147     end
0148 
0149     
0150 function [C, M] = matrix_from_number(nodes, num_levels, sn)
0151     check_scalar_level(num_levels);
0152     if isstruct(nodes)
0153         nodes = nodes.nodes;
0154     end
0155     M = diag([1,1,1]);
0156     
0157     zmax= max(nodes(:,3));
0158     zmin= min(nodes(:,3));
0159     levels = linspace(zmax,zmin, num_levels+2);
0160     levels = levels(2:end-1);
0161     
0162     C = (max(nodes) + min(nodes))/2;
0163     C(3) = levels(sn);
0164     
0165     
0166     
0167   
0168 function [C, M] = matrix_from_intercept(level, sn)   
0169    if nargin < 2
0170      sn = 1;
0171    end
0172    level = level(sn, :);
0173    R = eye(3);
0174    if length(level)==4
0175      theta = -level(4); % we rotate the model, user thinks of the image
0176      R = [cos(theta), -sin(theta), 0; sin(theta), cos(theta), 0; 0,0,1;];
0177    end
0178    
0179    % Infinities tend to cause issues -> replace with realmax
0180    % Don't need to worry about the sign of the inf
0181    level( isinf(level) | isnan(level) ) = realmax;
0182    level( level==0 ) =     1e-10; %eps;
0183 
0184    % Step 1: Choose a centre point in the plane
0185    %  Weight the point by it's inv axis coords
0186    invlev= 1./level;
0187    ctr= invlev / sum( invlev.^2 );
0188 
0189    % Step 2: Choose basis vectors in the plane
0190    %  First is the axis furthest from ctr
0191    [jnk, s_ax]= sort( - abs(level - ctr) );
0192    v1= [0,0,0]; v1(s_ax(1))= level(s_ax(1));
0193    v1= v1 - ctr;
0194    v1= v1 / norm(v1);
0195 
0196    % Step 3: Get off-plane vector, by cross product
0197    v2= [0,0,0]; v2(s_ax(2))= level(s_ax(2));
0198    v2= v2 - ctr;
0199    v2= v2 / norm(v2);
0200    v3= cross(v1,v2);
0201 
0202    % Step 4: Get orthonormal basis. Replace v2
0203    v2= cross(v1,v3);
0204 
0205    % Step 5: Get bases to point in 'positive directions'
0206    v1= v1 * (1-2*(sum(v1)<0));
0207    v2= v2 * (1-2*(sum(v2)<0));
0208    v3= v3 * (1-2*(sum(v3)<0));
0209    
0210    C = ctr;
0211    M = R * [v1;v2;v3];
0212 %   NODE= [v1;v2;v3] * (vtx' - ctr'*ones(1,nn) );
0213   
0214     
0215 function [C, M] = matrix_from_struct(level, sn)
0216   if nargin < 2
0217     sn = 1;
0218   end
0219   if ~isfield(level,'centre')
0220       error('Struct input ''level'' not understood.');
0221   end
0222   ctr = level.centre;
0223   if numel(ctr) == 3
0224       C = ctr; % make row
0225   elseif size(ctr,2) == 3 && size(ctr, 1) ~= 3
0226       C = ctr(:,sn);
0227   else
0228       C = ctr(sn,:);
0229   end
0230   C = C(:)'; % make row
0231   
0232   if isfield(level, 'rotation_matrix')
0233       M = level.rotation_matrix;
0234       return % no questions asked
0235   end
0236   
0237   if isfield(level, 'normal') && ~isfield(level,'normal_angle')
0238       level.normal_angle = level.normal; %quietly accept
0239   end
0240  
0241   if ~isfield(level,'normal_angle')
0242       error('Struct input ''level'' not understood.');
0243   end
0244   
0245   A = level.normal_angle(1:3);
0246   normA = norm(A);
0247   if normA == 0
0248     error('Normal vector has zero norm')
0249   end
0250   A = A(:) / normA; % normalize and make column
0251   B = [0,0,1]';
0252   R = eye(3);
0253   if length(level.normal_angle) == 4
0254     theta = level.normal_angle(4);
0255     R = [cos(theta), -sin(theta), 0; sin(theta), cos(theta), 0; 0,0,1;];
0256   end
0257    
0258   cross_AB = cross(A,B);
0259   if norm(cross_AB) == 0 % parallel vectors
0260     if A(3) > 0
0261       M = R*eye(3);
0262     elseif A(3) < 0
0263       M = R*diag([1,-1,-1]);
0264     end
0265     return
0266   end
0267   
0268   if 0
0269       % from https://math.stackexchange.com/questions/180418/
0270       % calculate-rotation-matrix-to-align-vector-a-to-vector-b-in-3d
0271       % works well, but leaves a rotation that's hard to control
0272       dot_AB = A(3);
0273       ssc = @(v) [0 -v(3) v(2); v(3) 0 -v(1); -v(2) v(1) 0];
0274       M = eye(3) + ssc(cross_AB) + ...
0275           ssc(cross_AB)^2*(1-dot_AB)/(norm(cross_AB)^2);
0276   else
0277       % here we define a nice coordinate system on the plane, and then
0278       % transform it to become the identity (ie. invert).
0279       
0280       % project each coordinate axis on the plane
0281       I = eye(3);
0282       for i = 1:3
0283           proj(:,i) = I(:,i) - (dot(I(:,i),A) / normA^2) * A;
0284       end
0285       norm_proj = vecnorm(proj);
0286       max_norm = max(norm_proj);
0287       
0288       % choose what happens based on which projection is longest
0289       M = zeros(3);
0290       M(:,3) = A;
0291       if norm_proj(3) == max_norm
0292           % projection of z becomes y
0293           M(:,2) = proj(:,3);
0294           M(:,1) = cross(M(:,2),M(:,3));
0295       elseif norm_proj(2) == max_norm
0296           % projection of y becomes y
0297           M(:,2) = proj(:,2);
0298           M(:,1) = cross(M(:,2),M(:,3));
0299       else
0300           % projection of x becomes x
0301           M(:,1) = proj(:,1);
0302           M(:,2) = cross(M(:,3),M(:,1));
0303       end
0304       M = inv(M);  
0305   end
0306   M = R*M ;  
0307   
0308 
0309 function check_scalar_level(num_levels)
0310    if num_levels <= 0
0311       error('EIDORS:WrongInput', 'Number of levels must be greater than 0.')
0312    end
0313 
0314   
0315 function do_unit_test
0316   unit_test_cmp('Scalar level', level_model_slice(5),5);
0317   unit_test_cmp('Vector level', level_model_slice([inf, inf, 3]),1);
0318   unit_test_cmp('Matrix level', level_model_slice([inf, inf, 3; inf, inf, 5]),2);
0319   LVL = [inf, inf, 3; inf, inf, 5];
0320   LVL(:,4) = rand(1,2);
0321   unit_test_cmp('Matrix level angle', level_model_slice(LVL),2);
0322   lvl.centre = [ 0,0,0];
0323   unit_test_cmp('Struct level', level_model_slice(lvl),1);
0324   lvl.centre(2,:) = [0,0,3];
0325   unit_test_cmp('Struct level', level_model_slice(lvl),2);
0326   lvl.normal_angle = rand(3,1)-.5;
0327   lvl.normal_angle = lvl.normal_angle / norm(lvl.normal_angle);
0328   lvl.normal_angle(4) = 2*pi * rand(1);
0329   [C, M] = level_model_slice(lvl);
0330   unit_test_cmp('Normal -> Matrix', M*lvl.normal_angle(1:3),[0,0,1]', 5*eps);
0331   [C, M] = level_model_slice([inf 5 inf]);
0332   lvl.normal_angle = [0,0,-1];
0333   [C, M] = level_model_slice(lvl);
0334   unit_test_cmp('Parallel normal', M,diag([1,-1,-1]));
0335   
0336   lvl.normal_angle = [0,0,0];
0337   try
0338     [C, M] = level_model_slice(lvl);
0339     fprintf('TEST: %20s = %4s \n','Expect error', 'Fail');  
0340   catch
0341     fprintf('TEST: %20s = %4s \n','Expect error', 'OK'); 
0342   end
0343   imdl = mk_common_model('n3r2',[16,2]); fmdl = imdl.fwd_model;
0344   nodes = level_model_slice(fmdl,[inf, inf, 1; inf, inf, 2],1);
Generated on Sun 29-Dec-2024 11:41:59 by m2html © 2005