function [ebsd_aus_pairs,ebsd_aus_triple,V_rec,F_rec,variants] = rec_aus_triple(ebsd,fcc2bcc_true,symmetries)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
symmetry1_string = symmetries{1}.mineral;

% Determine neighboring point pairs:
[V,F,I_FD] = spatialDecomposition([ebsd(symmetry1_string).prop.x(:), ebsd(symmetry1_string).prop.y(:)],ebsd(symmetry1_string).unitCell);
A_D = I_FD'*I_FD==1;
[Dl,Dr] = find(triu(A_D,1)); % lists of neighboring point pairs

% Reduce all orientations to quaternions:
O = quaternion(ebsd(symmetry1_string).orientations);
% Determine experimental misorientation set:
g1 = inv(O(Dl)).*O(Dr);

%% Set up the necessary symmetries and the simulated intermartensitic misorientations:

%The assumption is that we're dealing with misorientations between cubic symmetry:
cs = crystalSymmetry('432');%432 is used instead of m3m because it only gives 24 symmetry operators, which is enough

% Set up symmetry operators in the order that is assumed to produce a
% misorientation set similar to that used by Morito et al. and others:
cs_symmetries = rotation(cs.properGroup);

symrots_true(1,1) = cs_symmetries.rot(1);
symrots_true(2,1) = cs_symmetries.rot(21);
symrots_true(3,1) = cs_symmetries.rot(11);
symrots_true(4,1) = cs_symmetries.rot(16);
symrots_true(5,1) = cs_symmetries.rot(24);
symrots_true(6,1) = cs_symmetries.rot(8);
symrots_true(7,1) = cs_symmetries.rot(18);
symrots_true(8,1) = cs_symmetries.rot(10);
symrots_true(9,1) = cs_symmetries.rot(20);
symrots_true(10,1) = cs_symmetries.rot(3);
symrots_true(11,1) = cs_symmetries.rot(7);
symrots_true(12,1) = cs_symmetries.rot(23);
symrots_true(13,1) = cs_symmetries.rot(19);
symrots_true(14,1) = cs_symmetries.rot(2);
symrots_true(15,1) = cs_symmetries.rot(9);
symrots_true(16,1) = cs_symmetries.rot(22);
symrots_true(17,1) = cs_symmetries.rot(17);
symrots_true(18,1) = cs_symmetries.rot(12);
symrots_true(19,1) = cs_symmetries.rot(5);
symrots_true(20,1) = cs_symmetries.rot(15);
symrots_true(21,1) = cs_symmetries.rot(4);
symrots_true(22,1) = cs_symmetries.rot(14);
symrots_true(23,1) = cs_symmetries.rot(6);
symrots_true(24,1) = cs_symmetries.rot(13);

% Determine intermartensitic misorientations based on a predetermined orientation relationship (determined earlier by an iterative procedure):
bcc_trans =  inv(orientation(symrots_true*rotation(fcc2bcc_true),cs,1)) * (orientation(rotation(fcc2bcc_true(1)),cs,1));

% Reduce the simulated intermartensitic misorientations to quaternions:
g2 = quaternion(bcc_trans);

%% Determine intermartensitic misorientations

%Using a modified dot_outer_quat_cs, get the dot products, along with the index of the maximum dot product:
[ind,d] = dot_outer_quat_cs_modified(g2,g1,cs,cs);

%Calculate the angles corresponding to each dot product:
omega_1 = [real(2*acos(abs(d))) / degree];

% Get the minimum angle and index corresponding to the best-fitting
% simulated intermartensitic misorientation:
[omega_min_1, omega_min_1_I] = min(omega_1,[],2);

condition = omega_min_1<3 & omega_min_1_I~=1 & omega_min_1_I~=4;

% There are many misorientations corresponding to intragranular
% misorientations (omega_min_1_I == 1). There are also some misorientations
% that do not correspond to intermartensitic misorientations
% (omega_min_1>arbitrary_cutoff_value). Let us remove these from consideration:
cutoff_g1 = g1(condition);
cutoff_ind = ind(find(condition),:);
cutoff_omega_min_1 = omega_min_1(condition);
cutoff_omega_min_1_I = omega_min_1_I(condition);

%% Determine all necessary symmetry operators

%Using omega_min_1_I and its row number, we can pick the corresponding index from ind. First,
%we need to alter omega_min_1_I to a corresponding linear index for ind:
linearind_right_stuff = sub2ind(size(cutoff_ind),... %size of ind matrix
    [1:length(cutoff_omega_min_1)]',... % the row of each omega_min_1_I value
    [cutoff_omega_min_1_I]); % the index of the best-fitting intermartensitic misorientation on each row

% dot_outer_quat_cs_modified works by picking the best option out of a
% 24x24 matrix, from which we need the column value. All we get from ind is a linear index, however. %The rem command gives the division remainder for each linear index when divided by 24, giving the column value:
ind_necessary_2 = rem(cutoff_ind(linearind_right_stuff),24);
%The zero value given by the rem command actually corresponds to the 24th
%symmetry operator:
ind_necessary_2(ind_necessary_2 == 0) = 24;

%To obtain the final necessary symmetry operator, we need to go back to the list of indexes given by dot_outer_quat_modified. This time we need the row,
%which we can get by the command:
ind_necessary_3 = ceil(cutoff_ind(linearind_right_stuff)/24);

%The operation necessary to rotate a given orientation from the fcc to the
%bcc orientation is given by using symmetry operator 3:
candidate = rotation(fcc2bcc_true(1))*symrots_true(ind_necessary_3);
candidate = reshape(orientation(candidate,cs,1),length(candidate),1);

%The inverse of the candidate rotations can now be used to rotate each
%corresponding bcc orientation to the fcc orientation:
fcc_parents = orientation(O(Dr(condition)).*inv(candidate),symmetries{2},1);

candidate = rotation(fcc2bcc_true(1))*symrots_true(ind_necessary_2);
candidate = reshape(orientation(candidate,cs,1),length(candidate),1);
fcc_parents = [fcc_parents,orientation(O(Dl(condition)).*inv(candidate),symmetries{2},1)];

% We get the correct spatial map indexes for each rotated orientation:
Dr_restricted = Dr(condition);
Dl_restricted = Dl(condition);

miso = quaternion(inv(fcc_parents(:,1)).*fcc_parents(:,2));
[~,d_2] = dot_outer_quat_cs_modified(miso,quaternion.id,cs,cs);
angle_comparison_fccparents = [real(2*acos(abs(d_2))) / degree];
fcc_parents = fcc_parents(angle_comparison_fccparents<3,:);

Dr_restricted = Dr_restricted(angle_comparison_fccparents<3);
Dl_restricted = Dl_restricted(angle_comparison_fccparents<3);

[C,ia,ic] = unique(Dr_restricted);
z = ebsd(symmetry1_string);

    zz = fcc_parents(ia,1);

    %We are ready to create the partially reconstructed orientation map:
    ebsd_aus_pairs = ebsd;
    ebsd_aus_pairs(z(C).id).phase = 2;
    ebsd_aus_pairs(z(C).id).orientations = zz;
    
    [C,ia,ic] = unique(Dl_restricted);
    zz = fcc_parents(ia,2);
    ebsd_aus_pairs(z(C).id).phase = 2;
    ebsd_aus_pairs(z(C).id).orientations = zz;
    
    ebsd_aus_triple = ebsd;
    
    % How many times does a point appear in Dr?
    n_length = accumarray(Dr_restricted,1);
    
    % What neighbors does a point have?
    n = accumarray(Dr_restricted,1:length(Dr_restricted),[],@(x) {x});
    indices = cell2mat(n(find(n_length==3))')';
    
    %So we have three contenders for parent orientation for each point. They
    %should all be equal (no angular deviation). Let's see if this is true:
    zzz = inv(fcc_parents(indices(:,1))).*fcc_parents(indices(:,2)).*inv(inv(fcc_parents(indices(:,1))).*fcc_parents(indices(:,3)));
    
    % Determine indices where previously said checks out:
    ind_to_keep = (zzz.angle/degree)<5;
    ind_to_keep = indices(ind_to_keep,:);
    ind_to_keep = reshape(ind_to_keep,length(ind_to_keep)*3,1);
    
    % Set the austenite orientations for triple points:
    ebsd_aus_triple(z(Dr_restricted(ind_to_keep)).id).phase = 2;
    ebsd_aus_triple(z(Dr_restricted(ind_to_keep)).id).orientations = fcc_parents(ind_to_keep);
    
    ebsd_aus_triple(z(Dl_restricted(ind_to_keep)).id).phase = 2;
    ebsd_aus_triple(z(Dl_restricted(ind_to_keep)).id).orientations = fcc_parents(ind_to_keep);
%% PAG boundaries

connect = omega_min_1<5;
A_Db = sparse(double(Dl(~connect)),double(Dr(~connect)),true,length(ebsd(symmetry1_string)),length(ebsd(symmetry1_string)));
A_Db = A_Db | A_Db';
[ix,iy] = find(A_Db);
D_Fsub  = diag(sum(abs(I_FD(:,ix)) & abs(I_FD(:,iy)),2)>0);
I_FDint = D_Fsub*I_FD;
isBoundary = full(any(I_FDint,2));
F_2 = F(isBoundary,:);
I_FDint = I_FDint.'; I_FDint = I_FDint(:,isBoundary).';
[inUse,~,F_2] = unique(F_2);
V_PAG = V(inUse,:);
F_PAG = reshape(F_2,[],2);

% Packet boundaries
connect = omega_min_1<5 & omega_min_1_I <6;
A_Db = sparse(double(Dl(~connect)),double(Dr(~connect)),true,length(ebsd(symmetry1_string)),length(ebsd(symmetry1_string)));
A_Db = A_Db | A_Db';
[ix,iy] = find(A_Db);
D_Fsub  = diag(sum(abs(I_FD(:,ix)) & abs(I_FD(:,iy)),2)>0);
I_FDint = D_Fsub*I_FD;
isBoundary = full(any(I_FDint,2));
F_2 = F(isBoundary,:);
I_FDint = I_FDint.'; I_FDint = I_FDint(:,isBoundary).';
[inUse,~,F_2] = unique(F_2);
V_packet = V(inUse,:);
F_packet = reshape(F_2,[],2);

% Block boundaries
connect = omega_min_1<5 & omega_min_1_I==1;
A_Db = sparse(double(Dl(~connect)),double(Dr(~connect)),true,length(ebsd(symmetry1_string)),length(ebsd(symmetry1_string)));
A_Db = A_Db | A_Db';
[ix,iy] = find(A_Db);
D_Fsub  = diag(sum(abs(I_FD(:,ix)) & abs(I_FD(:,iy)),2)>0);
I_FDint = D_Fsub*I_FD;
isBoundary = full(any(I_FDint,2));
F_2 = F(isBoundary,:);
I_FDint = I_FDint.'; I_FDint = I_FDint(:,isBoundary).';
[inUse,~,F_2] = unique(F_2);
V_block = V(inUse,:);
F_block = reshape(F_2,[],2);

V_rec{1} = V_block;
V_rec{2} = V_packet;
V_rec{3} = V_PAG;
F_rec{1} = F_block;
F_rec{2} = F_packet;
F_rec{3} = F_PAG;

cutoff_omega_min_2_I = omega_min_1_I(omega_min_1<3 & omega_min_1_I~=1);
[hist_omega_I, bins_omega_I] = hist(cutoff_omega_min_2_I,1:24);
variants  = [hist_omega_I', bins_omega_I'];

function [ind,d] = dot_outer_quat_cs_modified(g2,g1,cs,ss)
% quick version that ignores inversion

cs_symmetries = rotation(cs.properGroup);

symrots_true(1,1) = cs_symmetries.rot(1);
symrots_true(2,1) = cs_symmetries.rot(21);
symrots_true(3,1) = cs_symmetries.rot(11);
symrots_true(4,1) = cs_symmetries.rot(16);
symrots_true(5,1) = cs_symmetries.rot(24);
symrots_true(6,1) = cs_symmetries.rot(8);
symrots_true(7,1) = cs_symmetries.rot(18);
symrots_true(8,1) = cs_symmetries.rot(10);
symrots_true(9,1) = cs_symmetries.rot(20);
symrots_true(10,1) = cs_symmetries.rot(3);
symrots_true(11,1) = cs_symmetries.rot(7);
symrots_true(12,1) = cs_symmetries.rot(23);
symrots_true(13,1) = cs_symmetries.rot(19);
symrots_true(14,1) = cs_symmetries.rot(2);
symrots_true(15,1) = cs_symmetries.rot(9);
symrots_true(16,1) = cs_symmetries.rot(22);
symrots_true(17,1) = cs_symmetries.rot(17);
symrots_true(18,1) = cs_symmetries.rot(12);
symrots_true(19,1) = cs_symmetries.rot(5);
symrots_true(20,1) = cs_symmetries.rot(15);
symrots_true(21,1) = cs_symmetries.rot(4);
symrots_true(22,1) = cs_symmetries.rot(14);
symrots_true(23,1) = cs_symmetries.rot(6);
symrots_true(24,1) = cs_symmetries.rot(13);


% g2rot = symmetrise_modified(quaternion(g2),cs,ss).'; % g2 x CS x SS
g2rot = (inv(symrots_true)*((g2*symrots_true).'));
% g2rot = (inv(symrots_true)*((g2*symrots_true).'));
g2rot = reshape(g2rot,24 * 24,[]);
g2rot = g2rot.';

q1 = [g1.a(:) g1.b(:) g1.c(:) g1.d(:)];
a2 = g2rot.a; b2 = g2rot.b; c2 = g2rot.c; d2 = g2rot.d;

% this is implicite dot_outer
d = abs(q1 * [a2(:,1).';b2(:,1).';c2(:,1).';d2(:,1).']); 
% allocate matrix for indices:
zalathar = ones(length(d),24);

for k=2:length(cs)*length(ss)
d_compare = d;
    d = max(d,abs(q1 * [a2(:,k).';b2(:,k).';c2(:,k).';d2(:,k).'])); % g1 x g2 x CS * SS
zalathar(d_compare<d) = k;
end
ind = zalathar;
end

function varargout = symmetrise_modified(q,CS,SS,varargin)
% symmetrcially equivalent orientations
%
% Input
%  q  - @quaternion
%  CS - crystal @symmetry
%  SS - specimen @symmetry
%
% Output
%  q - symmetrically equivalent orientations CS x SS x q
%
% See also
% CrystalSymmetries

% q = (q * CS).'; % CS x M
% if nargin>2 && length(SS)>1
%   q = SS * q;     % SS x (CS X M)
%   lSS = length(SS);
% else
%   lSS = 1;
% end

% Symmetry order correction according to Morito (this is not strictly necessary unless variant pair selection is of interest):
cs_symmetries = rotation(CS.properGroup);

symrots_true(1,1) = cs_symmetries.rot(1);
symrots_true(2,1) = cs_symmetries.rot(21);
symrots_true(3,1) = cs_symmetries.rot(11);
symrots_true(4,1) = cs_symmetries.rot(16);
symrots_true(5,1) = cs_symmetries.rot(24);
symrots_true(6,1) = cs_symmetries.rot(8);
symrots_true(7,1) = cs_symmetries.rot(18);
symrots_true(8,1) = cs_symmetries.rot(10);
symrots_true(9,1) = cs_symmetries.rot(20);
symrots_true(10,1) = cs_symmetries.rot(3);
symrots_true(11,1) = cs_symmetries.rot(7);
symrots_true(12,1) = cs_symmetries.rot(23);
symrots_true(13,1) = cs_symmetries.rot(19);
symrots_true(14,1) = cs_symmetries.rot(2);
symrots_true(15,1) = cs_symmetries.rot(9);
symrots_true(16,1) = cs_symmetries.rot(22);
symrots_true(17,1) = cs_symmetries.rot(17);
symrots_true(18,1) = cs_symmetries.rot(12);
symrots_true(19,1) = cs_symmetries.rot(5);
symrots_true(20,1) = cs_symmetries.rot(15);
symrots_true(21,1) = cs_symmetries.rot(4);
symrots_true(22,1) = cs_symmetries.rot(14);
symrots_true(23,1) = cs_symmetries.rot(6);
symrots_true(24,1) = cs_symmetries.rot(13);

q = (q * symrots_true).'; % CS x M
q = symrots_true * q;     % SS x (CS X M)

varargout{1} = reshape(q,24 * 24,[]); % (CSxSS) x M
end

end

