text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Q = INVERSEKINEMATIC_KUKA_KR5_ARC(robot, T)
%   Solves the inverse kinematic problem for the KUKA KR5 ARC robot
%   where:
%   robot stores the robot parameters.
%   T is an homogeneous transform that specifies the position/orientation
%   of the end effector.
%
%   A call to Q=INVERSEKINEMATIC__KUKA_KR5_ARC returns 8 possible solutions, thus,
%   Q is a 6x8 matrix where each column stores 6 feasible joint values.
%
%
%   Example code:
%
%   robot=load_robot('KUKA', 'KR160_R1570_nanoC');
%   q = [0 0 0 0 0 0];
%   T = directkinematic(robot, q);
%   %Call the inversekinematic for this robot
%   qinv = inversekinematic(robot, T);
%   check that all of them are feasible solutions!
%   and every Ti equals T
%   for i=1:8,
%        Ti = directkinematic(robot, qinv(:,i))
%   end
%	See also DIRECTKINEMATIC.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.
function q = inversekinematic_kuka_kr160_r1570_nanoC(robot, T)

%initialize q,
%eight possible solutions are generally feasible
q=zeros(6,8);

% %Evaluate the parameters
% theta = eval(robot.DH.theta);
d = eval(robot.DH.d);
L6=abs(d(6));

%T= [ nx ox ax Px;
%     ny oy ay Py;
%     nz oz az Pz];
Px=T(1,4);
Py=T(2,4);
Pz=T(3,4);

%Compute the position of the wrist, being W the Z component of the end effector's system
W = T(1:3,3);

% Pm: wrist position
Pm = [Px Py Pz]' + L6*W;

%first joint, two possible solutions admited:
% if q(1) is a solution, then q(1) + pi is also a solution
q1=atan2(Pm(2), Pm(1));

%solve for q2
q2_1=solve_for_theta2(robot, [q1 0 0 0 0 0 0], Pm);
%the other possible solution is q1 + pi
q2_2=solve_for_theta2(robot, [q1+pi 0 0 0 0 0 0], Pm);

%solve for q3
q3_1=solve_for_theta3(robot, [q1 0 0 0 0 0 0], Pm);
%solver for q3 for both cases
q3_2=solve_for_theta3(robot, [q1+pi 0 0 0 0 0 0], Pm);

%the next matrix doubles each column. For each two columns, two different
%configurations for theta4, theta5 and theta6 will be computed. These
%configurations are generally referred as wrist up and wrist down solution
q = [q1         q1         q1        q1       q1+pi   q1+pi   q1+pi   q1+pi;
     q2_1(1)    q2_1(1)    q2_1(2)   q2_1(2)  q2_2(1) q2_2(1) q2_2(2) q2_2(2);
     q3_1(1)    q3_1(1)    q3_1(2)   q3_1(2)  q3_2(1) q3_2(1) q3_2(2) q3_2(2);
     0          0          0         0         0      0       0       0;
     0          0          0         0         0      0       0       0;
     0          0          0         0         0      0       0       0];
%leave only the real part of the solutions
q=real(q);

%Note that in this robot, the joint q3 has a non-simmetrical range. In this
%case, the joint ranges from 60 deg to -219 deg, thus, the typical normalizing
%step is avoided in this angle (the next line is commented). When solving
%for the orientation, the solutions are normalized to the [-pi, pi] range
%only for the theta4, theta5 and theta6 joints.

%normalize q to [-pi, pi]
q(1,:) = normalize(q(1,:));
q(2,:) = normalize(q(2,:));

% solve for the last three joints
% for any of the possible combinations (theta1, theta2, theta3)
for i=1:2:size(q,2),
    qtemp = solucion_muneca_esferica(robot, q(:,i), T, 1,'geometric'); %wrist up
    qtemp(4:6)=normalize(qtemp(4:6));
    q(:,i)=qtemp;

    qtemp = solucion_muneca_esferica(robot, q(:,i), T, -1, 'geometric'); %wrist up
    qtemp(4:6)=normalize(qtemp(4:6));
    q(:,i+1)=qtemp;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve for second joint theta2, two different
% solutions are returned, corresponding
% to elbow up and down solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function q2 = solve_for_theta2(robot, q, Pm)

%Evaluate the parameters
theta = eval(robot.DH.theta);
d = eval(robot.DH.d);
a = eval(robot.DH.a);
alpha = eval(robot.DH.alpha);

%See geometry
L2=abs(a(2));
L3=abs(d(4));
A2 = abs(a(3));

%See geometry of the robot
%compute L4
L4 = sqrt(A2^2 + L3^2);

%The inverse kinematic problem can be solved as in the IRB 140 (for example)

%given q1 is known, compute first DH transformation
T01=dh(robot, q, 1);

%Express Pm in the reference system 1, for convenience
p1 = inv(T01)*[Pm; 1];

r = sqrt(p1(1)^2 + p1(2)^2);

beta = atan2(-p1(2), p1(1));
gamma = real(acos((L2^2+r^2-L4^2)/(2*r*L2)));

%return two possible solutions
%elbow up and elbow down
%the order here is important and is coordinated with the function
%solve_for_theta3
q2(1) = pi/2 - beta - gamma; %elbow up
q2(2) = pi/2 - beta + gamma; %elbow down

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve for third joint theta3, two different
% solutions are returned, corresponding
% to elbow up and down solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function q3 = solve_for_theta3(robot, q, Pm)

%Evaluate the parameters
theta = eval(robot.DH.theta);
d = eval(robot.DH.d);
a = eval(robot.DH.a);
alpha = eval(robot.DH.alpha);

%See geometry
L2=abs(a(2));
L3=abs(d(4));

A2 = abs(a(3));

%See geometry of the robot
%compute L4
L4 = sqrt(A2^2 + L3^2);

%the angle phi is fixed
phi=acos((A2^2+L4^2-L3^2)/(2*A2*L4));

%given q1 is known, compute first DH transformation
T01=dh(robot, q, 1);

%Express Pm in the reference system 1, for convenience
p1 = inv(T01)*[Pm; 1];

r = sqrt(p1(1)^2 + p1(2)^2);

eta = real(acos((L2^2 + L4^2 - r^2)/(2*L2*L4)));

%return two possible solutions
%elbow up and elbow down solutions
%the order here is important
q3(1) = phi - eta;
q3(2) = -2*pi +( phi + eta); %%%?????????????????????????????????????????

function q = solucion_muneca_esferica(robot, q, T, wrist, method)

switch method

     %algebraic solution
    case 'algebraic'
        T01=dh(robot, q, 1);
        T12=dh(robot, q, 2);
        T23=dh(robot, q, 3);

        Q=inv(T23)*inv(T12)*inv(T01)*T;

        %detect the degenerate case when q(5)=0, this leads to zeros
        % in Q13, Q23, Q31 and Q32 and Q33=1
        thresh=1e-12;
        %detect if q(5)==0
        % this happens when cos(q5) in the matrix Q is close to 1
        if abs(Q(3,3)-1)>thresh
            %normal solution
            if wrist==1 %wrist up
                q(4)=atan2(-Q(2,3),-Q(1,3));
                q(6)=atan2(-Q(3,2),Q(3,1));
                %q(5)=atan2(-Q(3,2)/sin(q(6)),Q(3,3));
            else %wrist down
                q(4)=atan2(-Q(2,3),-Q(1,3))-pi;
                q(6)=atan2(-Q(3,2),Q(3,1))+pi;
                %q(5)=atan2(-Q(3,2)/sin(q(6)),Q(3,3));
            end
            if abs(cos(q(6)+q(4)))>thresh
                cq5=(Q(1,1)+Q(2,2))/cos(q(4)+q(6))-1;
            end
            if abs(sin(q(6)+q(4)))>thresh
                cq5=(-Q(1,2)+Q(2,1))/sin(q(4)+q(6))-1;
            end
            if abs(sin(q(6)))>thresh
                sq5=-Q(3,2)/sin(q(6));
            end
            if abs(cos(q(6)))>thresh
                sq5=Q(3,1)/cos(q(6));
            end
            q(5)=atan2(sq5,cq5);

        else %degenerate solution, in this case, q4 cannot be determined,
             % so q(4)=0 is assigned
            if wrist==1 %wrist up
                q(4)=0;
                q(5)=0;
                q(6)=atan2(-Q(1,2)+Q(2,1),Q(1,1)+Q(2,2));
            else %wrist down
                q(4)=-pi;
                q(5)=0;
                q(6)=atan2(-Q(1,2)+Q(2,1),Q(1,1)+Q(2,2))+pi;
            end

        end

       %geometric solution
    case 'geometric'
        % T is the noa matrix defining the position/orientation of the end
        % effector's reference system
        vx6=T(1:3,1);
        vz5=T(1:3,3); % The vector a z6=T(1:3,3) is coincident with z5

        % Obtain the position and orientation of the system 3
        % using the already computed joints q1, q2 and q3
        T01=dh(robot, q, 1);
        T12=dh(robot, q, 2);
        T23=dh(robot, q, 3);
        T03=T01*T12*T23;

        vx3=T03(1:3,1);
        vy3=T03(1:3,2);
        vz3=T03(1:3,3);

        % find z4 normal to the plane formed by z3 and a
        vz4=cross(vz3, vz5);	% end effector's vector a: T(1:3,3)

        % in case of degenerate solution,
        % when vz3 and vz6 are parallel--> then z4=0 0 0, choose q(4)=0 as solution
        if norm(vz4) <= 0.000001
            if wrist == 1 %wrist up
                q(4)=0;
            else
                q(4)=-pi; %wrist down
            end
        else
            %this is the normal and most frequent solution
            cosq4=wrist*dot(vy3,vz4);
            sinq4=wrist*dot(-vx3,vz4);
            q(4)=atan2(sinq4, cosq4);
        end
        %propagate the value of q(4) to compute the system 4
        T34=dh(robot, q, 4);
        T04=T03*T34;
        vx4=T04(1:3,1);
        vy4=T04(1:3,2);

        % solve for q5
        cosq5=dot(-vy4,vz5);
        sinq5=dot(vx4,vz5);
        q(5)=atan2(sinq5, cosq5);

        %propagate now q(5) to compute T05
        T45=dh(robot, q, 5);
        T05=T04*T45;
        vx5=T05(1:3,1);
        vy5=T05(1:3,2);

        % solve for q6
        cosq6=dot(vx6,vx5);
        sinq6=dot(vx6,vy5);
        q(6)=atan2(sinq6, cosq6);

    otherwise
        disp('no method specified in solve_spherical_wrist');
end

repo: ARTE


text:
function G = grad3(V,T)
  % GRAD3 Compute the numerical gradient operator for tet meshes in 3d
  %
  % G = grad(V,F)
  %
  % Inputs:
  %   V  #vertices by dim list of mesh vertex positions
  %   T  #elements by 3 list of mesh tet indices
  % Outputs:
  %   G  #elements*dim by #V Gradient operator
  %
  % Example:
  %   L = cotmatrix(V,T);
  %   G  = grad3(V,T);
  %   vol = volume(V,T);
  %   GMG = -G'*repdiag(diag(sparse(vol)),3)*G;
  %

  warning('Deprecated. Call grad directly.');
  G = grad(V,T);
end

author: alecjacobson

repo: gptoolbox


text:
function [bioQualifiers,standardQualifiers] = getBioQualifiers()
% Get a list of possible bioqualifiers
%
% USAGE:
%    annotationQualifiers = getBioQualifiers()
%
% OUTPUT:
%    bioQualifiers:    A cell array of bioqualifier strings in descending order.
%    standardQualifiers:        A n x 2 cell array with the first column representing the standard identifiers and the second column representing the alternative name.

bioQualifiers = {'version','taxon','relation',...
                 'propertyBearer','property',...
                 'parthood','part','origin','occursIn',...
                 'isVersionOf','isRelatedTo','isPropertyOf',...
                 'isPartOf','isInstanceOf','isHomologTo','isEncodedBy',...
                 'isDescribedBy','isDerivedFrom','is','instance','identity',...
                 'hypernym','homolog','hasVersion',...
                 'hasTaxon','hasProperty','hasPart','hasInstance',...
                 'encodes','encoder','encodement',...
                 'description','container','class'};

standardQualifiers = {'is', 'identity';...
                      'isDerivedFrom', 'origin';...
                      'isDescribedBy', 'description';...
                      'isInstanceOf', 'class';...
                      'hasInstance', 'instance';...
                      'encodes', 'encodement';...
                      'hasPart', 'part';...
                      'hasProperty', 'property';...
                      'hasVersion', 'version';...
                      'isEncodedBy', 'encoder';...
                      'isHomologTo', 'homolog';...
                      'isPartOf', 'parthood';...
                      'isPropertyOf', 'propertyBearer';...
                      'isVersionOf', 'hypernym';...
                      'occursIn', 'container';...
                      'hasTaxon', 'taxon'};

author: opencobra

repo: cobratoolbox


text:
function det = load_det(varargin)
% Load detection results of different local/pairwise models

%% parse parameters
opts = struct;
% model type
opts.fix_ann.x_off = 0;
opts.fix_ann.y_off = 0;
opts.fix_ann.w = inf;
opts.fix_ann.h = inf;

opts.regression = struct;               % detection regression parameters
opts.regression.param = [0 0 1 1];
opts.regression.fix_ann = struct;       % fix annotation if needed
opts.regression.fix_ann.x_off = 0;
opts.regression.fix_ann.y_off = 0;
opts.regression.fix_ann.w = inf;
opts.regression.fix_ann.h = inf;
opts.regression.warp = 'none';

opts.nms = struct;
opts.nms.nmsIntersectionOverAreaThreshold = 0.3;
opts.nms.numBoundingBoxMax = inf;

opts.regression_before_nms = false;     % do regression before nms

opts.det = struct;
opts.det.modeltype = 'local';
opts.det.thres = -inf;
opts.det.path_format = '';
opts.det.scoretype = 'raw';
opts.det.as_filter = false;             % remove *weird* aspect ratio
opts.det.as_ratio = 1.5;

opts.im_path_format = '';
opts.im_set = [];

opts.viz = struct;                      % visulization
opts.viz.doviz = false;
opts.viz.max_det = 5;

opts.verbose = true;
opts.progress_part = 5;

opts = vl_argparse(opts, varargin);

det = struct;
numimage = length(opts.im_set);

fprintf('Loading detections');
progress_part_num = ceil(numimage/opts.progress_part);

for i = 1:numimage
    if (opts.verbose)
        if (i==numimage)
            fprintf('...100%%\n');
        else
            if ~mod(i,progress_part_num)
                fprintf('...%d%%', i*100/(progress_part_num*opts.progress_part));
            end
        end
    end

    if (~iscell(opts.im_set))
        idname = opts.im_set(i);
        det_path = sprintf(opts.det.path_format, idname);
        im_path = sprintf(opts.im_path_format , idname);

    else
        idname = opts.im_set{i};
        det_path = sprintf(opts.det.path_format, opts.im_set{i});
        im_path = sprintf(opts.im_path_format , opts.im_set{i});
    end
    det(i).path = det_path;
    det(i).impath = im_path;
    det(i).id = idname;

    im = imread(im_path);
    [img_h, img_w, ~] = size(im);

    BB = load_BB(opts.det.modeltype, det_path);
    if (isempty(BB))
        det(i).bb = BB;
        continue;
    end

%     if (strcmp(detector,'globalmasked')~=0)
%         scale_factor = 224/size(im,2);
%         pad_size = floor((size(im,2)-size(im,1))/2);
%         BB_unpad = BB/scale_factor;
%         BB_unpad(:, 2) = BB_unpad(:, 2) - pad_size;
%         BB = BB_unpad;
%     end

    %thresholding lowscore detections
    BB = BB(BB(:,5)>opts.det.thres, :);
    if (isempty(BB)) det(i).bb = BB; continue; end

    %remove *weird* aspect ratio if needed
    if (opts.det.as_filter)
        w = BB(:,3);
        h = BB(:,4);
        as = w./h;
        BB = BB(as < opts.det.as_ratio && as > 1/opts.det.as_ratio, :);
    end
    if (isempty(BB)) det(i).bb = BB; continue; end

    if (~opts.regression_before_nms)
        top = selectBoundingBoxesNonMaxSup(BB(:,1:4), BB(:,5), opts.nms);
        BB = BB(top, :);
        BB = do_regression(BB, img_w, img_h, opts.regression);
    else
        BB = do_regression(BB, img_w, img_h, opts.regression);
        top = selectBoundingBoxesNonMaxSup(BB(:,1:4), BB(:,5), opts.nms);
        BB = BB(top, :);
    end

    if (opts.viz.doviz)
        fig = figure;
        imshow(im); hold on;
        for z=1:min(opts.viz.max_det, size(BB,1))
            rectangle('position', BB(z, 1:4), 'edgecolor', 'r');
            text(double(BB(z, 1))+5, double(BB(z, 2))+5, num2str(BB(z, 5)), 'color', 'yellow');
        end
        disp('Press any key to continue...');
        pause;
        close(fig);
    end

    det(i).bb = BB;
end
end

author: aosokin


text:
function aminz = i4vec_aminz ( n, a )

%*****************************************************************************80
%
%% I4VEC_AMINZ returns the smallest nonzero magnitude in an I4VEC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries to be checked.
%
%    Input, integer A(N), the vector to be checked.
%
%    Output, integer AMINZ, the value of the smallest nonzero magnitude.
%    If all entries are zero, AMINZ is 0.
%
  range = find ( a(1:n) ~= 0 );

  if ( isempty ( range ) )
    aminz = 0;
  else
    aminz = min ( abs ( a(range) ) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef FeaturesMatcher < handle
    %FEATURESMATCHER  Feature matchers class
    %
    % See also: cv.Stitcher, cv.FeaturesFinder, cv.DescriptorMatcher
    %

    properties (SetAccess = private)
        % Object ID
        id
    end

    methods
        function this = FeaturesMatcher(matcherType, varargin)
            %FEATURESMATCHER  Constructor
            %
            %     obj = cv.FeaturesMatcher(matcherType)
            %     obj = cv.FeaturesMatcher(matcherType, 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __matcherType__ One of:
            %   * __BestOf2NearestMatcher__ A "best of 2 nearest" matcher.
            %     Features matcher which finds two best matches for each
            %     feature and leaves the best one only if the ratio between
            %     descriptor distances is greater than the threshold
            %     `MatchConf`.
            %   * __BestOf2NearestRangeMatcher__
            %   * __AffineBestOf2NearestMatcher__ A "best of 2 nearest"
            %     matcher that expects affine trasformation between images.
            %     Features matcher similar to `BestOf2NearestMatcher` which
            %     finds two best matches for each feature and leaves the best
            %     one only if the ratio between descriptor distances is
            %     greater than the threshold `MatchConf`. Unlike
            %     `BestOf2NearestMatcher` this matcher uses affine
            %     transformation (affine trasformation estimate will be placed
            %     in `matches_info`).
            %
            % ## Options
            % The following are options accepted by all matchers:
            %
            % * __TryUseGPU__ Should try to use GPU or not. default false
            % * __MatchConf__ Match distances ration threshold. Confidence for
            %   feature matching step. default 0.3
            % * __NumMatchesThresh1__ Minimum number of matches required for
            %   the 2D projective transform estimation used in the inliers
            %   classification step. default 6
            % * __NumMatchesThresh2__ Minimum number of matches required for
            %   the 2D projective transform re-estimation on inliers.
            %   default 6
            %
            % The following are options for the various algorithms:
            %
            % ### `BestOf2NearestRangeMatcher`
            % * __RangeWidth__ Range width. default 5
            %
            % ### `AffineBestOf2NearestMatcher`
            % * __FullAffine__ whether to use full affine transformation with
            %   6 degress of freedom (cv.estimateAffine2D) or reduced
            %   transformation with 4 degrees of freedom
            %   (cv.estimateAffinePartial2D) using only rotation, translation
            %   and uniform scaling. default false
            %
            % See also: cv.FeaturesMatcher.match
            %
            this.id = FeaturesMatcher_(0, 'new', matcherType, varargin{:});
        end

        function delete(this)
            %DELETE  Destructor
            %
            %     obj.delete()
            %
            % See also: cv.FeaturesMatcher
            %
            if isempty(this.id), return; end
            FeaturesMatcher_(this.id, 'delete');
        end

        function typename = typeid(this)
            %TYPEID  Name of the C++ type (RTTI)
            %
            %     typename = obj.typeid()
            %
            % ## Output
            % * __typename__ Name of C++ type
            %
            typename = FeaturesMatcher_(this.id, 'typeid');
        end
    end

    %% FeaturesMatcher
    methods
        function collectGarbage(this)
            %COLLECTGARBAGE  Frees unused memory allocated before if there is any
            %
            %     obj.collectGarbage()
            %
            % See also: cv.FeaturesMatcher.FeaturesMatcher
            %
            FeaturesMatcher_(this.id, 'collectGarbage');
        end

        function tf = isThreadSafe(this)
            %ISTHREADSAFE  Check if matcher is thread safe
            %
            %     tf = obj.isThreadSafe()
            %
            % ## Output
            % * __tf__ True, if it's possible to use the same matcher instance
            %   in parallel, false otherwise.
            %
            % See also: cv.FeaturesMatcher.FeaturesMatcher
            %
            tf = FeaturesMatcher_(this.id, 'isThreadSafe');
        end

        function matches_info = match(this, features1, features2)
            %MATCH  Performs images matching
            %
            %     matches_info = obj.match(features1, features2)
            %
            % ## Input
            % * __features1__ First image features. See cv.FeaturesFinder.
            % * __features2__ Second image features. See cv.FeaturesFinder.
            %
            % ## Output
            % * **matches_info** Found matches. Structure containing
            %   information about matches between two images. It's assumed
            %   that there is a transformation between those images.
            %   Transformation may be homography or affine transformation
            %   based on selected matcher. Struct with the following fields:
            %   * **src_img_idx** Images indices (optional).
            %   * **dst_img_idx** Images indices (optional).
            %   * __matches__ Matches. A 1-by-N structure array with the
            %     following fields:
            %     `{'queryIdx', 'trainIdx', 'imgIdx', 'distance'}`
            %   * **inliers_mask** Geometrically consistent matches mask.
            %   * **num_inliers** Number of geometrically consistent matches.
            %   * __H__ Estimated transformation.
            %   * __confidence__ Confidence two images are from the same
            %     panorama.
            %
            % See also: cv.FeaturesMatcher.match_pairwise,
            %  cv.FeaturesFinder.find
            %
            matches_info = FeaturesMatcher_(this.id, 'match', features1, features2);
        end

        function pairwise_matches = match_pairwise(this, features, varargin)
            %MATCH_PAIRWISE  Performs images matching
            %
            %     pairwise_matches = obj.match_pairwise(features)
            %     pairwise_matches = obj.match_pairwise(features, 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __features__ Features of the source images.
            %   See cv.FeaturesFinder.
            %
            % ## Output
            % * **pairwise_matches** Found pairwise matches.
            %
            % ## Options
            % * __Mask__ Mask indicating which image pairs must be matched.
            %   default empty
            %
            % The function is parallelized with the TBB library.
            %
            % See also: cv.FeaturesMatcher.match, cv.FeaturesFinder.find
            %
            pairwise_matches = FeaturesMatcher_(this.id, 'match_pairwise', features, varargin{:});
        end
    end

    %% Auxiliary methods
    methods (Static)
        function str = matchesGraphAsString(pairwise_matches, conf_threshold)
            %MATCHESGRAPHASSTRING  Covert matches to graph
            %
            %     str = cv.FeaturesMatcher.matchesGraphAsString(pairwise_matches, conf_threshold)
            %
            % ## Input
            % * **pairwise_matches** Pairwise matches.
            % * **conf_threshold** Threshold for two images are from the same
            %   panorama confidence.
            %
            % ## Output
            % * __str__ matches graph represented in DOT language. Labels
            %   description: `Nm` is number of matches, `Ni` is number of
            %   inliers, `C` is confidence.
            %
            % Returns matches graph representation in DOT language.
            %
            str = FeaturesMatcher_(0, 'matchesGraphAsString', pairwise_matches, conf_threshold);
        end

        function indices = leaveBiggestComponent(features, pairwise_matches, conf_threshold)
            %LEAVEBIGGESTCOMPONENT  Leave biggest component
            %
            %     indices = cv.FeaturesMatcher.leaveBiggestComponent(features, pairwise_matches, conf_threshold)
            %
            % ## Input
            % * __features__ Features of the source images.
            % * **pairwise_matches** Pairwise matches.
            % * **conf_threshold** Threshold for two images are from the same
            %   panorama confidence.
            %
            % ## Output
            % * __indices__ array of image indices (0-based).
            %
            indices = FeaturesMatcher_(0, 'leaveBiggestComponent', features, pairwise_matches, conf_threshold);
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function M = euclideanfactory(m, n)
% Returns a manifold struct to optimize over m-by-n matrices.
%
% function M = euclideanfactory(m, n)
%
% Returns M, a structure describing the Euclidean space of m-by-n matrices
% equipped with the standard Frobenius distance and associated trace inner
% product as a manifold for Manopt.

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Dec. 30, 2012.
% Contributors:
% Change log:
%  July 5, 2013 (NB): added egred2rgrad, ehess2rhess, mat, vec, tangent.

    if ~exist('n', 'var') || isempty(n)
        n = 1;
    end

    M.name = @() sprintf('Euclidean space R^(%dx%d)', m, n);

    M.dim = @() m*n;

    M.inner = @(x, d1, d2) d1(:).'*d2(:);

    M.norm = @(x, d) norm(d, 'fro');

    M.dist = @(x, y) norm(x-y, 'fro');

    M.typicaldist = @() sqrt(m*n);

    M.proj = @(x, d) d;

    M.egrad2rgrad = @(x, g) g;

    M.ehess2rhess = @(x, eg, eh, d) eh;

    M.tangent = M.proj;

    M.exp = @exp;
    function y = exp(x, d, t)
        if nargin == 3
            y = x + t*d;
        else
            y = x + d;
        end
    end

    M.retr = M.exp;

	M.log = @(x, y) y-x;

    M.hash = @(x) ['z' hashmd5(x(:))];

    M.rand = @() randn(m, n);

    M.randvec = @randvec;
    function u = randvec(x) %#ok<INUSD>
        u = randn(m, n);
        u = u / norm(u, 'fro');
    end

    M.lincomb = @lincomb;
    function v = lincomb(x, a1, d1, a2, d2) %#ok<INUSL>
        if nargin == 3
            v = a1*d1;
        elseif nargin == 5
            v = a1*d1 + a2*d2;
        else
            error('Bad usage of euclidean.lincomb');
        end
    end

    M.zerovec = @(x) zeros(m, n);

    M.transp = @(x1, x2, d) d;

    M.pairmean = @(x1, x2) .5*(x1+x2);

    M.vec = @(x, u_mat) u_mat(:);
    M.mat = @(x, u_vec) reshape(u_vec, [m, n]);
    M.vecmatareisometries = @() true;

end

author: andrewssobral

repo: lrslibrary


text:
function bcdof = BoundaryCondition(typeBC,coordinates)

%--------------------------------------------------------------------------
%   Purpose:
%           To determine the boundary conditions degree of freedom
%   Synopsis:
%           bcdof = BoundaryCondition(typeBC,coordinates)
%
%   Variable Description:
%           bcdof - boundary condition degree of freedom
%                   dof's are (UZ,RX,RY)
%           typeBC - string which gives the state of boundary condition
%           coordinates - geometric coordinates of nodes
%
%--------------------------------------------------------------------------

L1 = find(coordinates(:,2)==min(coordinates(:,2))) ; % at y = 0 (along X-axes)
L2 = find(coordinates(:,1)==max(coordinates(:,1))) ; % at x = a (along Y-axes)
L3 = find(coordinates(:,2)==max(coordinates(:,2))) ; % at y = b (along X-axes)
L4 = find(coordinates(:,1)==min(coordinates(:,1))) ; % at x = 0 (along Y-axes)
n = length(L1) ;

switch typeBC
    case 'ss-ss-ss-ss'
        disp('plate is simply supported at all the edges')
        dofL1 = zeros(1,2*n) ;
        dofL2 = zeros(1,2*n) ;
        dofL3 = zeros(1,2*n) ;
        dofL4 = zeros(1,2*n) ;
        for i = 1:n
            i1 = 2*(i-1)+2 ;
            i2 = i1-1 ;
            dofL1(i1) = 3*L1(i);
            dofL1(i2) = 3*L1(i)-2 ;
            dofL3(i1) = 3*L3(i) ;
            dofL3(i2) = 3*L3(i)-2 ;
        end
        for i = 1:n
            i1 = 2*(i-1)+2 ;
            i2 = i1-1 ;
            dofL2(i1) = 3*L2(i)-1 ;
            dofL2(i2) = 3*L2(i)-2 ;
            dofL4(i1) = 3*L4(i)-1 ;
            dofL4(i2) = 3*L4(i)-2 ;
        end
        L1UL3 = union(dofL1,dofL3) ;
        L2UL4 = union(dofL2,dofL4) ;
        bcdof = union(L1UL3,L2UL4) ;

    case 'c-c-c-c'
        disp('plate is clamped at all the edges')
        dofL1 = zeros(1,2*n) ;
        dofL2 = zeros(1,2*n) ;
        dofL3 = zeros(1,2*n) ;
        dofL4 = zeros(1,2*n) ;
        for i = 1:n
            i1 = 2*(i-1)+2 ;
            i2 = i1-1 ;
            dofL1(i1) = 3*L1(i)-1;
            dofL1(i2) = 3*L1(i)-2 ;
            dofL3(i1) = 3*L3(i)-1 ;
            dofL3(i2) = 3*L3(i)-2 ;
        end
        for i = 1:n
            i1 = 2*(i-1)+2 ;
            i2 = i1-1 ;
            dofL2(i1) = 3*L2(i) ;
            dofL2(i2) = 3*L2(i)-2 ;
            dofL4(i1) = 3*L4(i) ;
            dofL4(i2) = 3*L4(i)-2 ;
        end
        L1UL3 = union(dofL1,dofL3) ;
        L2UL4 = union(dofL2,dofL4) ;
        bcdof = union(L1UL3,L2UL4) ;
end

author: Sable

repo: mcbench-benchmarks


text:
function [fg, fg_unclassified]=dtiFindMoriTracts(dt6File, outFile, fgFile, Atlas, showFig, saveQuench, saveMrDiffusion, useJhuFa, useInterhemisphericSplit, useRoiBasedApproach)
%Categorizes fibers based on Mori white matter atlas.
%
%  dtiFindMoriTracts(dt6File, [outFile=fullfile(fileparts(dt6File), 'fibers', 'MoriGroups.mat')], fgFile,...
%      [Atlas='MNI_JHU_tracts_prob.nii.gz'], [showFig=false],[saveQuench=false],...
%      [saveMrDiffusion=false], [useJhuFa=false],...
%      [useInterhemisphericSplit=true], [useRoiBasedApproach=true]);
%
% By default whole brain tractography is performed and then
% Mori-classified. Otherwise, provide a file name fgFile for a fiber group
% if willing to Mori-classify an existing fiber group.
%
% Input parameters:
% fgFile                   - A file with previously tracked elsewhere fibers
%                          to be categorized.
% Atlas                    - probabilistic atlas defining probabilities for
%                          each voxel to be passed by a fiber within each of
%                          atlas fiber groups. We usually use Mori atlas
%                          supplied with fsl: MNI_JHU_tracts_prob.nii.gz.
%                          This atlas is not symmetric. For a "symmetrified"
%                          atlas use 'MNI_JHU_tracts_prob_Symmetric.nii.gz'
%                          but we strongly recommend using the original
%                          atlas.
% useInterhemisphericSplit - cut fibers crossing between hemispheres with a
%                          midsaggital plane below z=-10. This is to get
%                          rid of
% useRoiBasedApproach      - use the approach describing in Zhang (2008)
%                           Neuroimage 42. For each of the 20 Mori Groups 2
%                           critical ROIs are computed by spatially
%                           transforming ROIs provided in
%                           templates/MNI_JHU_tracts_ROIs. A fiber becomes
%                           a candidate to being labeled as a part of a
%                           given Mori group if this fiber "passes through"
%                           both critical ROIs for that Mori group. Our
%                           modification of Zhang (2008) approach: In case
%                           a single fiber is a candidate for >1 Mori
%                           group, respective cumulative probabilities are
%                           computed with probabilistic Mori atlas, then
%                           compared. useRoiBasedApproach can take the
%                           following values: (1) 'false' (to not use the
%                           approach); (2) 'true'(use the approach; the
%                           minimal distance from a fiber to ROI to count
%                           as "a fiber  is crossing the ROI" minDist=2mm);
%                           (3) a scalar value for minDist in mm; (4) a
%                           vector with the first value being minDist and
%                           the second value being a flag 1/0 for whether
%                           ROIs should be recomputed (and overwritten) for
%                           this subject (1, by default). This is useful
%                           because sometimes if you  are rerunning
%                           dtiFindMoriTracts with different parameters,
%                           you do not need to recompute ROIs. E.g., to
%                           avoid recomputing  ROIs and use minDist of 4mm
%                           one would pass [useRoiBasedApproach=[4 0]];
%
% Output parameters:
% fg              - fibers structure containing all fibers assigned to
%                   one of Mori Groups. Respective group labeles are stored
%                   in fg.subgroups field.
% fg_unclassified - fiber structure containing the rest of the (not Mori) fibers.
%
% Example:
%    dt6File = '/biac3/wandell4/data/reading_longitude/dti_y1234/at040918/dti06trilinrt/dt6.mat';
%    fg = dtiFindMoriTracts(dt6File);
%
% See also: dtiSplitInterhemisphericFibers.m
%
% Usage Notes:
%    Text explaining the usage of this code for publication can be found at
%    the end of this code along with references.
%
% (c) Stanford University, Vistalab

% HISTORY:
% 2008.10.08 RFD wrote it. 2009.01.08 EIR added an option of Mori-categorizing
%            a provided FG (do not supply a fgFile parameter, or supply []
%            if want whole brain tractography performed)
% 2009.01.09 RFD removed the feature that discarded fibers that scored too
%            high on more than one atlas group. This was causing too many
%            good fibers to be removed, especially for groups 19 & 20
%            (arcuate).
% 2009.01.30 EIR modified to write out Mori labels as fg.subgroup
% 2009.02.01 EIR added an option to use symmetrified Mori atlas
% 2009.05.31 EIR modified to output fghandle in addition to fg
% 2009.07.21 EIR modified to output unclassified fibers
% 2009.08.28 EIR added interhemispheric fiber split below ACPC and removed
%                fghandle output (because the parent FG is modified to add
%                split fibers.
% 2009.09.27 ER  added useRoiBasedApproach option which implements Zhang
%                (2008) Neuroimage 42 method.
% 2009.10.2 ER:  no longer need to fit Quench state file labels to 8 grp
%                limit.

%% Check INPUTS

if ~exist('showFig', 'var')|| isempty(showFig)
    showFig= false;
end
if ~exist('saveQuench', 'var')|| isempty(saveQuench)
    saveQuench = false;
end
if ~exist('saveMrDiffusion', 'var')|| isempty(saveMrDiffusion)
    saveMrDiffusion= false;
end
if ~exist('useJhuFa', 'var')|| isempty(useJhuFa)
    useJhuFa = false; %JhuFa should never be used
end
if ~exist('useInterhemisphericSplit', 'var')|| isempty(useInterhemisphericSplit)
    useInterhemisphericSplit=true;
end
if ~exist('useRoiBasedApproach', 'var')|| isempty(useRoiBasedApproach)
    useRoiBasedApproach=true;
end

if useRoiBasedApproach==false
    recomputeROIs = false;
elseif length(useRoiBasedApproach)<2
    recomputeROIs=1;
else
    recomputeROIs=useRoiBasedApproach(2);
end

if recomputeROIs
    display('dtiFindMoriTracts: You chose to recompute ROIs');
end

% E.g., to avoid recomputing  ROIs and use minDist of 4mm one would pass [useRoiBasedApproach=[4 0]];
if isnumeric(useRoiBasedApproach)
    minDist = useRoiBasedApproach(1);
    useRoiBasedApproach = 'true';
else
    minDist=2; %.89;
end
display(['Fibers that get as close to the ROIs as ' num2str(minDist) 'mm will become candidates for the Mori Groups']);

if(~exist('Atlas','var') || isempty(Atlas))
    % Default scenario: use original Mori Atlas
    Atlas='MNI_JHU_tracts_prob.nii.gz';
end

if(~exist('fgFile','var') || isempty(fgFile))
    % Default scenario: perform whole brain tractogrpahy
    wholeBrainFlag=1;
else
    wholeBrainFlag=0;
end

if(~exist('outFile','var') || isempty(outFile))
    bd = fileparts(dt6File);
    if(isempty(bd)), bd = pwd; end
    outFile = fullfile(bd,'fibers', 'MoriGroups.mat');
end

%%

% Load a dt6 file
dt = dtiLoadDt6(dt6File);

tdir = fullfile(fileparts(which('mrDiffusion.m')), 'templates');
spm_defaults; global defaults; params = defaults.normalise.estimate;
if(useJhuFa)
    template = fullfile(tdir,'MNI_JHU_FA.nii.gz');
    alignIm = dtiComputeFA(dt.dt6);
    params.cutoff = 19;
    params.reg = 0.09;
else
    % Spatially normalize it with the MNI (ICBM) template
    template = fullfile(tdir,'MNI_JHU_T2.nii.gz');
    alignIm = mrAnatHistogramClip(double(dt.b0),0.3,0.99);
end

[sn, Vtemplate, invDef] = mrAnatComputeSpmSpatialNorm(alignIm, dt.xformToAcpc, template, params);

% check the normalization
mm = diag(chol(Vtemplate.mat(1:3,1:3)'*Vtemplate.mat(1:3,1:3)))';
bb = mrAnatXformCoords(Vtemplate.mat,[1 1 1; Vtemplate.dim]);
alignIm_sn = mrAnatResliceSpm(alignIm, sn, bb, [2 2 2], [1 1 1 0 0 0], 0);
tIm = mrAnatResliceSpm(double(Vtemplate.dat), inv(Vtemplate.mat), bb, [2 2 2], [1 1 1 0 0 0], 0);
im(:,:,:,1) = uint8(tIm);
im(:,:,:,2) = uint8(round(clip(alignIm_sn)*255));
im(:,:,:,3) = im(:,:,:,2);

if(showFig)

    showMontage(im);
else
    if(~exist(fileparts(outFile),'dir')), mkdir(fileparts(outFile)); end
    imwrite(makeMontage(im),[outFile(1:end-4) '_snCheck.png']);
end

% Load the Mori atlas maps and the corresponding label files
% ldir = fileparts(which('dtiGetBrainlabel.m'));
moriTracts = niftiRead(fullfile(tdir, Atlas));
% 15 is a subregion of 19 and 16 a subregion of 20. To better separate them,
% we subtract 19 from 15 and 20 from 16.
moriTracts.data(:,:,:,15) = moriTracts.data(:,:,:,15)-moriTracts.data(:,:,:,19);
moriTracts.data(:,:,:,16) = moriTracts.data(:,:,:,16)-moriTracts.data(:,:,:,20);

labels = readTab(fullfile(tdir,'MNI_JHU_tracts_prob.txt'),',',false);
labels = labels(1:20,2);
% If you wanted to inverse-normalize the maps to this subject's brain:
% invDef.outMat = moriTracts.qto_ijk;
% bb = mrAnatXformCoords(dt.xformToAcpc,[1 1 1; size(dt.b0)]);
% tprob = mrAnatResliceSpm(tprob, invDef, bb, dt.mmPerVoxel, [1 1 1 0 0 0]);

if wholeBrainFlag
    % Track all white matter fibers in the native subject space. We do this by
    % seeding all voxels with high FA (>0.3).
    faThresh = 0.30;
    opts.stepSizeMm = 1;
    opts.faThresh = .2; %0.15
    opts.lengthThreshMm = [50 250];
    opts.angleThresh = 50;
    opts.wPuncture = 0.2;
    opts.whichAlgorithm = 1;
    opts.whichInterp = 1;
    opts.seedVoxelOffsets = [0.25 0.75];
    opts.offsetJitter = 0.1;
    fa = dtiComputeFA(dt.dt6);
    fa(fa>1) = 1; fa(fa<0) = 0;
    roiAll = dtiNewRoi('all');
    mask = dtiCleanImageMask(fa>=faThresh);
    [x,y,z] = ind2sub(size(mask), find(mask));
    clear mask fa;
    roiAll.coords = mrAnatXformCoords(dt.xformToAcpc, [x,y,z]);
    clear x y z;
    fg = dtiFiberTrack(dt.dt6, roiAll.coords, dt.mmPerVoxel, dt.xformToAcpc, 'wholeBrain', opts);
    clear roiAll
    % wholeBrainFGFile=fullfile(fileparts(dt6File), 'fibers', 'all.mat');
    % dtiWriteFiberGroup(fg, wholeBrainFGFile);
    % fghandle.parent=wholeBrainFGFile;
else
    % get fg from a file
    fg = fgRead(fgFile);
    %fg = dtiLoadFiberGroup(fgFile);
    %fghandle.parent=fgFile;
end

%Cut the fibers below acpc, to disentangle CST& ATL crossing at the pons
%level.
if useInterhemisphericSplit
    fgname=fg.name;
    [fg]=dtiSplitInterhemisphericFibers(fg, dt, -10);
    fg.name=fgname; %To avoid "bilaterally split" extra long comment in the fg name.
end

%Because we know that we are looking for major Mori tracts, we should
%keep only 5 points or longer fibers (otherwise spline in contrack
%crashes).  This is important for Blue matter procedures planned, but wont hurt otherwise (for Mori detection) either.
if sum(cellfun(@length, fg.fibers)<5)~=0
    if isfield(fg, 'subgroup')&&~isempty(fg.subgroup)
        fg.subgroup(cellfun(@length, fg.fibers)<5)=[];
    end
    if isfield(fg, 'seeds')&&~isempty(fg.seeds)
        fg.seeds(cellfun(@length, fg.fibers)<5, :)=[];
    end
    fg.fibers(cellfun(@length, fg.fibers)<5)=[];
    fprintf('dtiFindMoriTracts: Removing %s fibers with 5 points or less \n', num2str(sum(cellfun(@length, fg.fibers)<5)));
end

% Warp the fibers in 'fg' to the MNI standard space:
fg_sn = dtiXformFiberCoords(fg, invDef);

% moriTracts.data is a an XxYxZx20 array contianing the 20 Mori probability
% atlases (range is 0-100 where 100 represents p(1)).
sz = size(moriTracts.data);

% fg_sn fiber coords are in MNI space- now convert them to atlas space by
% applying the affine xform from the atlas NIFTI header. SInce the atlas is
% already in MNI space, this transform will just account for any
% translation and scale differences between the atlas maps and the MNI
% template used to compute our sn.
fgCoords = mrAnatXformCoords(moriTracts.qto_ijk, horzcat(fg_sn.fibers{:}));
clear fg_sn;   % what we need from fg_sn is now stored in fgCoords
fgLen = cellfun('size',fg.fibers,2);

% Now loop over the 20 atlases and get the atlas probability score for each
% fiber point. We collapse the scores across all pints in a fiber by taking
% the mean. Below, we will use these 20 mean scores to categorize the fibers.
% TO DO: consider doing something more sophisticated than taking the mean.
fp = zeros(sz(4),numel(fg.fibers));
for(ii=1:sz(4))
    % Get the Mori atlas score for each point in the fibers using
    % trilinear interpolation.
    p = myCinterp3(double(moriTracts.data(:,:,:,ii))/100, sz([1,2]), sz(3), fgCoords(:,[2,1,3]));
    % The previous line interpolated one giant array with all fiber points
    % concatenated. The next loop will separate the coordinates back into
    % fibers and take the mean score for the points within each fiber.
    fiberCoord = 1;
    for(jj=1:numel(fg.fibers))
        fp(ii,jj) = nanmean(p([fiberCoord:fiberCoord+fgLen(jj)-1]));
        fiberCoord = fiberCoord+fgLen(jj);
    end
end
clear p fgCoords;

if useRoiBasedApproach
    %Warp Mori ROIs to individual space; collect candidates for each fiber
    %group based on protocol of 2 or > ROIs a fiber should travel thru.
    %The following ROIs are saved within
    %trunk/mrDiffusion/templates/MNI_JHU_tracts_ROIs folder and are created using MNI template as
    %described in Wakana et al.(2007) Neuroimage 36 with a single modification:
    %For SLFt Roi2, they recommend drawing the ROI at the AC level, whereas we
    %use a lisce just inferior of CC splenium. The reason for this modification
    %is that Wakana et al. ACPC aligned images appear different from MNI images
    %(the latter we use for defininng ROIs). If defining SLFt-Roi2 on a slice
    %actually at the AC level (althought highly consistently across human
    %raters), many SLFt fibers were not correctly labeled as they extend
    %laterally into temporal lobe just above the aforementioned ROI plane.

    moriRois={'ATR_roi1_L.nii.gz',  'ATR_roi2_L.nii.gz'; 'ATR_roi1_R.nii.gz', 'ATR_roi2_R.nii.gz'; ...
        'CST_roi1_L.nii.gz', 'CST_roi2_L.nii.gz'; 'CST_roi1_R.nii.gz',  'CST_roi2_R.nii.gz'; ...
        'CGC_roi1_L.nii.gz', 'CGC_roi2_L.nii.gz'; 'CGC_roi1_R.nii.gz', 'CGC_roi2_R.nii.gz'; ...
        'HCC_roi1_L.nii.gz', 'HCC_roi2_L.nii.gz'; 'HCC_roi1_R.nii.gz', 'HCC_roi2_R.nii.gz';...
        'FP_R.nii.gz', 'FP_L.nii.gz'; ...
        'FA_L.nii.gz', 'FA_R.nii.gz'; ...
        'IFO_roi1_L.nii.gz', 'IFO_roi2_L.nii.gz'; 'IFO_roi2_R.nii.gz', 'IFO_roi1_R.nii.gz'; ...
        'ILF_roi1_L.nii.gz', 'ILF_roi2_L.nii.gz'; 'ILF_roi1_R.nii.gz', 'ILF_roi2_R.nii.gz'; ...
        'SLF_roi1_L.nii.gz', 'SLF_roi2_L.nii.gz'; 'SLF_roi1_R.nii.gz', 'SLF_roi2_R.nii.gz'; ...
        'UNC_roi1_L.nii.gz', 'UNC_roi2_L.nii.gz'; 'UNC_roi1_R.nii.gz', 'UNC_roi2_R.nii.gz'; ...
        'SLF_roi1_L.nii.gz', 'SLFt_roi2_L.nii.gz'; 'SLF_roi1_R.nii.gz', 'SLFt_roi2_R.nii.gz'};

    midSaggitalRoi = dtiRoiMakePlane([0, dt.bb(1, 2), dt.bb(1, 3); 0 , dt.bb(2, 2) , dt.bb(2, 3)], 'midsaggital', 'g');

    keep1 = zeros(length(fg.fibers), size(moriRois, 1)); keep2=zeros(length(fg.fibers), size(moriRois, 1));

    [fgOut,contentiousFibers,InterHemisphericFibers] = dtiIntersectFibersWithRoi([], 'not', [], midSaggitalRoi, fg);
    %NOTICE: ~keep3 (not "keep3") will mark fibers that DO NOT cross midSaggitalRoi.
    keep3 = repmat(InterHemisphericFibers, [1 size(moriRois, 1)]);

    fgCopy=fg; fgCopy.subgroup=[];
    for roiID=1:size(moriRois, 1)

        ROI_img_file=fullfile(tdir, 'MNI_JHU_tracts_ROIs',  [moriRois{roiID, 1}]);
        if recomputeROIs
            [RoiFileName, invDef, roi]=dtiCreateRoiFromMniNifti(dt6File, ROI_img_file, invDef, true);
        else
            RoiFileName=fullfile(fileparts(dt6File), 'ROIs',  [prefix(prefix(ROI_img_file, 'short'), 'short') '.mat']);
            load(RoiFileName);

        end
        [fgOut,contentiousFibers, keep1(:, roiID)] = dtiIntersectFibersWithRoi([], 'and', minDist, roi, fg);
        keepID1=find(keep1(:, roiID));

        ROI_img_file=fullfile(tdir, 'MNI_JHU_tracts_ROIs',  [moriRois{roiID, 2}]);
        if recomputeROIs
            [RoiFileName, invDef, roi]=dtiCreateRoiFromMniNifti(dt6File, ROI_img_file, invDef, true);
        else
            RoiFileName=fullfile(fileparts(dt6File), 'ROIs',  [prefix(prefix(ROI_img_file, 'short'), 'short') '.mat']);
            load(RoiFileName);
        end

        %To speed up the function, we intersect with the second ROI not all the
        %fibers, but only those that passed first ROI.
        fgCopy.fibers=fg.fibers(keepID1(keepID1>0));
        [a,b, keep2given1] = dtiIntersectFibersWithRoi([], 'and', minDist, roi, fgCopy);
        keep2(keepID1(keep2given1), roiID)=true;

    end
    clear fgOut contentiousFibers keepID
    %Note: forceps major and minor should NOT have interhemipsheric fibers
    %excluded
    keep3(:, 9:10)=keep3(:, 9:10).*0;
    fp(~(keep1'&keep2'&~keep3'))=0;
    %Also note: Tracts that cross through slf_t rois should be automatically
    %classified as slf_t, without considering their probs.
    fp(19, (keep1(:, 19)'&keep2(:, 19)'&~keep3(:, 19)'))=max(fp(:));
    fp(20, (keep1(:, 20)'&keep2(:, 20)'&~keep3(:, 20)'))=max(fp(:));

end

% We have a set of atlas scores for each each fiber. To categorize the
% fibers, we will find the atlas with the highest score (using 'sort').
[atlasScore,atlasInd] = sort(fp,1,'descend');
% Eliminate fibers that don't match any of the atlases very well:
unclassified=atlasScore(1,:)==0; %ER 09.15.09 dropped this threshold to zero
goodEnough = atlasScore(1,:)~=0; % RFD 09.01.09: removed "& atlasScore(1,:)>sum(atlasScore(2:end,:),1);"

for ii=1:sz(4)
    curAtlasFibers{ii} = find(atlasInd(1,:)==ii & goodEnough);
    %    if(showFig)
    %        fc = round(mrAnatXformCoords(inv(dt.xformToAcpc), horzcat(fg.fibers{curAtlasFibers{ii}})));
    %        im = dt.b0;
    %        inds = sub2ind(size(dt.b0),fc(:,1),fc(:,2),fc(:,3));
    %        im(inds) = 1;
    %        makeMontage3(im,dt.b0,dt.b0);
    %        set(gcf,'name',labels{ii},'NumberTitle','off');
    %    end
end

% We now have a cell array (curAtlasFibers) that contains 20 arrays, each
% listing the fiber indices for the corresponding atlas group. E.g.,
% curAtlasFibers{3} is a list of indices into fg.fibers that specify the
% fibers belonging to group 3.

%Create a FG for unclassified fibers
fg_unclassified=fg;
fg_unclassified.name=[fg.name ' not Mori Groups'];
fg_unclassified.fibers = fg.fibers(unclassified); %prepare fg output
if ~isempty(fg.seeds)
    fg_unclassified.seeds = fg.seeds(unclassified);
end
fg_unclassified.subgroup=zeros(size(fg.fibers(unclassified)))'+(1+sz(4));
fg_unclassified.subgroupNames(1)=struct('subgroupIndex', 1+sz(4), 'subgroupName', 'NotMori');

% Modify fg.fibers to discard the fibers that didn't make it into any of
% the atlas groups:
fg.name=[fg.name ' Mori Groups'];
%fghandle.name=fg.name;
fg.fibers = fg.fibers([curAtlasFibers{:}]); %prepare fg output
%fghandle.ids=horzcat(curAtlasFibers{:}); %prepare fghandle output

if ~isempty(fg.seeds)
    fg.seeds = fg.seeds([curAtlasFibers{:}],:);
end

% We changed the size of fg.fibers by discarding the uncategorized fibers,
% so we need to create a new array to categorize the fibers. This time we
% make an array with one entry corresponding to each fiber, with integer
% values indicating to which atlas group the corresponding fiber belongs.
fg.subgroup = zeros(1,numel(fg.fibers));
curInd = 1;
for(ii=1:numel(curAtlasFibers))
    fg.subgroup(curInd:curInd+numel(curAtlasFibers{ii})-1) = ii;
    %fghandle.subgroup(curInd:curInd+numel(curAtlasFibers{ii})-1) = ii;
    curInd = curInd+numel(curAtlasFibers{ii});
    %Save labels for the fiber subgroups within the file
    fg.subgroupNames(ii)=struct('subgroupIndex', ii, 'subgroupName', labels(ii));
    %fghandle.subgroupNames(ii)=struct('subgroupIndex', ii, 'subgroupName', labels(ii));
end

% Save in mrDiffusion format:
if (saveMrDiffusion)
    if(~exist(fileparts(outFile),'dir')), mkdir(fileparts(outFile)); end
    dtiWriteFiberGroup( fg, outFile);
end

if(saveQuench)
    if(~exist(fileparts(outFile),'dir')), mkdir(fileparts(outFile)); end
    % Also save a Quench pdb file/state file to show the fiber groups.
    % Should we use dtiWriteFibersPdb?
    dtiWriteFibersPdb(fg,dt.xformToAcpc, [outFile(1:end-4) '.pdb']);
    % merge some of the groups to fit into the Quench 8-group limit:
    fgInds = fg.subgroup;
    dtiQuenchSaveFibersState(fgInds, [outFile(1:end-4) '.qst']);
end

return;

%%

% To run this on a bunch of subjects:
bd = '/biac3/wandell4/data/reading_longitude/dti_y1234';
[dt6Files,subCodes,subDirs,subIn] = findSubjects(fullfile(bd,'*'), 'dti06rt');
subIn = unique(subIn);
doThese = [1:numel(dt6Files)]; % strmatch('pt0',subCodes)';
for(ii=doThese)
    fiberDir = fullfile(fileparts(dt6Files{ii}),'Mori');
    if(~exist(fiberDir,'dir')), mkdir(fiberDir); end
    outBase = fullfile(fiberDir,'MoriGroups');
    if(~exist([outBase '.mat'],'file'))
        fprintf('Processing %s...\n',dt6Files{ii});
        dtiFindMoriTracts(dt6Files{ii},outBase);
    end
end

subs = []; doThese = [];
for(ii=1:numel(subIn))
    cur = strmatch([subIn{ii} '0'],subCodes);
    if(numel(cur)==4)
        subs = [subs cur'];
    end
end
for(ii=subs)
    fiberDir = fullfile(fileparts(dt6Files{ii}),'Mori');
    if(exist(fiberDir,'dir')&&exist(fullfile(fiberDir,'MoriGroups.mat')))
        %disp([subCodes{ii} ' is finished.']);
    else
        disp([subCodes{ii} ' is MISSING.']);
        doThese = [doThese ii];
    end
end

% To analyze a bunch of subjects:
outDir = '/biac3/wandell4/data/reading_longitude/moriGroupAnalysis';
bd = '/biac3/wandell4/data/reading_longitude/dti_y1234';
[dt6Files,subCodes,subDirs,subLetters] = findSubjects(fullfile(bd,'*'), 'dti06rt');
nGroups = 20;
nSubs = numel(dt6Files);
sgVol = zeros(nSubs,nGroups);
for(ii=1:nSubs)
    fname = fullfile(fileparts(dt6Files{ii}),'Mori','MoriGroups.mat');
    fprintf('Processing %s (%d of %d)...\n',fname,ii,numel(dt6Files));
    fg = dtiReadFibers(fname);
    for(jj=1:nGroups)
        sgCoords = fg.fibers(fg.subgroup==jj);
        % crude volume measure:
        sgVol(ii,jj) = size(unique(round(horzcat(sgCoords{:}))','rows'),1);
    end
end

% Get the labels for the 20 groups:
labels = readTab(which('MNI_JHU_tracts_prob.txt'),',',false);
labels = labels(:,2);
save(fullfile(outDir,'moriGroupSum.mat'),'bd','labels','sgVol','dt6Files','subCodes','subLetters');

outDir = '/biac3/wandell4/data/reading_longitude/moriGroupAnalysis';
load(fullfile(outDir,'moriGroupSum.mat'));
nGroups = size(sgVol,2);
nSubs = size(sgVol,1);
mn = mean(sgVol)/1000;
sd = std(sgVol)/1000;
rng = [max(sgVol); min(sgVol)]/1000;
for(ii=1:nGroups)
    fprintf('%50s:\t %0.2fcc (%0.3f, %0.2f-%0.2f)\n',labels{ii},mn(ii),sd(ii),rng(1,ii),rng(2,ii));
end

[bd, colNames, sc, subYr] = dtiGetBehavioralData(subCodes);

yr = 4;
bdVars = {'Passage Comprehension.', 'Rapid Naming.', 'Word attack ss.', 'Phonological Awareness.', 'Calculation.','WISC Full-Scale IQ.','DTI Age.','Sex (1=male)'};
for(ii=1:numel(bdVars))
    for(jj=1:2:nGroups-1)
        if(bdVars{ii}(end)=='.'), bVar = sprintf('%s%d',bdVars{ii},yr);
        else bVar = bdVars{ii}; end
        s = subYr==yr;
        x = bd(s,strmatch(bVar,colNames,'exact'));
        if(~isempty(x))
            y = sgVol(s,jj+1)-sgVol(s,jj);  % sgVol(s,jj);
            gv = ~isnan(x)&~isnan(y);
            x = x(gv); y = y(gv);
            %figure(34); plot(x,y,'.');
            [p,r,df] = myStatTest(x,y,'r');
            fprintf('%25s vs. %50s:\t r=%+0.3f (p=%0.2g, df=%d)\n',bVar,labels{jj},r,p,df);
        end
    end
end

%% TEXT FOR PUBLICATION

% The fiber tracts from the whole brain tractography were automatically
% classified into twenty fiber structures as defined in the JHU
% white-matter tractography atlas (Wakana et al., 2007) using a modified
% form of the reference ROI approach ( [Wakana et al., 2007] , [Zhang et
% al., 2008] , [Hua et al., 2008] and [Zhang et al., 2010] ). Specifically,
% we manually defined reference ROIs (rROIs) describing two waypoints for
% each of the 20 major white matter tracts described in (Wakana et al.,
% 2007). These ROIs were drawn in MNI space on the ICBM-DTI-81 atlas by two
% experts. The rROIs were warped from MNI space into each individual's
% diffusion space and fibers were retained if they passed through any pair
% of rROIs. Note that some fibers passed through the rROIs for more than
% one major tract. Thus, we applied an additional inclusion criterion by
% warping each fiber to MNI space and measuring the approximate overlap
% between the fiber points and each of the major tracts of the
% probabilistic JHU tractography atlas (Wakana et al., 2007). The fiber was
% classified as representing the JHU tract with which it had the highest
% degree of overlap.
%

%% REFERENCES:
%
% Wakana, Setsu et al. 2007. “Reproducibility of quantitative tractography
% methods applied to cerebral white matter.” NeuroImage 36(3): 630-44.
%
% Zhang, Weihong, Alessandro Olivi, Samuel J Hertig, Peter van Zijl, et
% al. 2008. “Automated fiber tracking of human brain white matter using
% diffusion tensor imaging.” NeuroImage 42(2): 771-7.
%
% Hua, K, J Zhang, S Wakana, H Jiang, et al. 2008. “Tract probability maps
% in stereotaxic spaces: analyses of white matter anatomy and
% tract-specific quantification.” Neuroimage 39(1): 336-347.
%
% Zhang et al., 2010 Y. Zhang, J. Zhang, K. Oishi, A. Faria, H. Jiang, X.
% Li, K. Akhter, P. Rosa-Neto, G. Pike, A. Evans, A. Toga, R. Woods, J.
% Mazziotta, M. Miller, P. van Zijl and S. Mori, Atlas-guided tract
% reconstruction for automated and comprehensive examination of the white
% matter anatomy, . Neuroimage, 52 4 (2010), pp. 1289–1301 (May).

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft