%% Set plotting preferences
setMTEXpref('xAxisDirection','east');
setMTEXpref('zAxisDirection','intoPlane');
setMTEXpref('figSize','normal');

%% Calculate grains
[grains,ebsd.grainId,ebsd.mis2mean] = calcGrains(ebsd('Iron'),'angle',10*degree,'boundary','tight');

%% Plot IPF for raw data
oM = ipdfHSVOrientationMapping(ebsd('indexed'));
oM.inversePoleFigureDirection = xvector;
color_1 = oM.orientation2color(ebsd('indexed').orientations);

figure;  
plot(ebsd('indexed'),color_1);
hold on;
plot(grains.boundary);

%% Plot KAM for RAW data
nearest_neighbours=1; %Number of nearest neighbours for KAM
misorientation_threshold=2; %Maximum misorientation allowed inside the kernel in degrees

%Raw data
kam_raw=KAM(ebsd,'threshold',misorientation_threshold*degree,'order',nearest_neighbours)./degree; %First nearest neighbours with analysis range 0-2 degrees

figure;
plot(ebsd,kam_raw);
mtexColorMap(parula);
title('KAM, raw data');
mtexColorbar
CLim(gcm,[0,1]);
hold on;
plot(grains.boundary);

%% Smooth ebsd data with half-quadratic
%{
Filters available in MTEX:
F = meanFilter;
F.numNeighbours = 1; % 3x3 window

F = medianFilter;
F.numNeighbours = 3; % this corresponds to a 7x7 window

F = KuwaharaFilter;
F.numNeighbours = 1; %3x3 window, avoid large windows

F = splineFilter;

F = halfQuadraticFilter; %Rather slow
F.alpha = 0.15;
%}

%With these params in beta the denoising is very minor in newest beta
F = halfQuadraticFilter;
F.alpha = 0.15; %Decrease alpha to reduce filtering.
fprintf('Smoothing ebsd data. \n');
ebsd_smoothed = smooth(ebsd('indexed'),F,'fill'); %Filter and Fill all non-indexed pixels
ebsd_grid=gridify(ebsd_smoothed('Iron'));

%{
%Other params for HQ from 5.0.1
 properties
    weight = @(t,eps,th) (t<=th)./(t.^2+eps^2).^(1/2); % weight function handle b^*(t,epsilon)
    eps = 1e-1;                             % parameter for the weight function
    tol = 1e-5                              % tolerance for gradient descent
    threshold = 5*degree                    % threshold for subgrain boundaries
%}

%Remove orientations from non-indexed top area
%{
  %Syntax to read black-white mask
  mask=imread(['d:\temp\crop_mask.png']);
  mask=im2bw(mask,0.5);
%}
index_mask=find(mask==0);
ebsd_grid(index_mask).orientations=orientation(nanquaternion,ebsd_grid.CS); %Assign nan-orientation to masked area 

%% Plot IPF for filtered data
oM = ipdfHSVOrientationMapping(ebsd_grid('indexed'));
oM.inversePoleFigureDirection = xvector;
color_2 = oM.orientation2color(ebsd_grid('indexed').orientations);

figure; 
plot(ebsd_grid('indexed'),color_2);
hold on;
plot(grains.boundary);

%% Plot KAM for filtered data
kam_filtered=KAM(ebsd_grid,'threshold',misorientation_threshold*degree,'order',nearest_neighbours)./degree; %First nearest neighbours with analysis range 0-2 degrees

figure;
plot(ebsd_grid,kam_filtered);
mtexColorMap(parula);
title('KAM, filtered data');
mtexColorbar
CLim(gcm,[0,1]);
hold on;
plot(grains.boundary);

%% Iterative filtering
%2nd pass HQ
F = halfQuadraticFilter;
F.alpha = 0.15;
ebsd_grid_smoothed=smooth(ebsd_grid('indexed'),F,'fill');

%Kuwahara
F = KuwaharaFilter;
F.numNeighbours = 1; %3x3 window
ebsd_grid_smoothed_2=smooth(ebsd_grid_smoothed('indexed'),F,'fill');
ebsd_grid_smoothed_2(index_mask).orientations=orientation(nanquaternion,ebsd_grid.CS); %Assign nan-orientation to masked area 

%% Plot IPF for iterative filtering
oM = ipdfHSVOrientationMapping(ebsd_grid_smoothed_2('indexed'));
oM.inversePoleFigureDirection = xvector;
color_3 = oM.orientation2color(ebsd_grid_smoothed_2('indexed').orientations);

figure; 
plot(ebsd_grid_smoothed_2('indexed'),color_3);
hold on;
plot(grains.boundary);

%% Plot KAM for iterative filtering
nearest_neighbours=1; %Number of nearest neighbours for KAM
misorientation_threshold=2; %Maximum misorientation allowed inside the kernel in degrees

kam_filtered_iter=KAM(ebsd_grid_smoothed_2,'threshold',misorientation_threshold*degree,'order',nearest_neighbours)./degree; %First nearest neighbours with analysis range 0-2 degrees

figure;
plot(ebsd_grid_smoothed_2,kam_filtered_iter);
mtexColorMap(parula);
title('KAM, iterative filtered data');
mtexColorbar
CLim(gcm,[0,1]);
hold on;
plot(grains.boundary);

%% Define misorientation profile between two dislocation cells
id_point_1=69495;
id_point_2=73611;

%Take EBSD line sectoin and calculate misorientation
%xy1 linesec start, xy2 linesec end

%Reference point
xy1=[ebsd_grid(id_point_1).x ebsd_grid(id_point_1).y];
xy2=[ebsd_grid(id_point_2).x ebsd_grid(id_point_2).y];
lineSec=[xy1;xy2];

[ebsd_line,ebsd_line_dist] = spatialProfile(ebsd_grid_smoothed_2('Iron'),lineSec); %Change ebsd variable here to plot other filtering options
ref_ori=ebsd_line(1).orientations;
ebsd_line_ids=ebsd_line.id;
misori_to_origin=angle(ref_ori,ebsd_line.orientations)./degree;

%Plot linesegment overlay on KAM map
figure;
title(['KAM, iterative filtering data']);
plot(ebsd_grid_smoothed_2,kam_filtered_iter);
mtexColorMap(parula);
CLim(gcm,[0,1]);
mtexColorbar;
hold on;
plot(grains.boundary);

%Plot line-segment
plot(lineSec(:,1),lineSec(:,2),'-w','linewidth',2);
plot(lineSec(:,1),lineSec(:,2),'-k','linewidth',1);
scatter(xy1(1,1),xy1(1,2),50,'w','filled','o');
scatter(xy1(1,1),xy1(1,2),30,'k','filled','o');

%% Plot misorientation profiles
figure('Color',[1 1 1]);
plot(ebsd_line_dist,misori_to_origin,'k','linestyle','-','linewidth',1);
hold on;
xlabel('Distance along line (µm)'); ylabel('Rotation (deg)');
xlim([nanmin(ebsd_line_dist) nanmax(ebsd_line_dist)]);
ylabel('Misorientation (deg)');

yyaxis right;
plot(ebsd_line_dist,kam_filtered_iter(ebsd_line_ids),'b','linestyle','-','linewidth',1);
ylabel('Misorientation (deg)');
legend('Misorientation to reference','KAM');

%Misorientation angle
grod_angle=angle(ebsd_line(1).orientations,ebsd_line(end).orientations)./degree

%Misorientation axis in specimen coordinates
grod_axis=axis(ebsd_line(1).orientations,ebsd_line(end).orientations)

%Misorientation axis in crystal coordinates
misor=inv(ebsd_line(1).orientations.*ebsd_line(end).orientations);
grod_axis_cs=axis(misor,'noSymmetry')




