%% test identification of continuously rounded particles
% Base: octahedron
% Roundness: RelAddSph varies from 0 to 1
% Orientation: Set of random orientations
clear all
clc
% close all

%% orientations to track
rng(1)
nSample = 20;
x = randn(3, nSample);
dirMatrix = bsxfun(@times, 1./sqrt(sum(x.*x, 1)), x);
% plot(dirMatrix(1,:), dirMatrix(2,:), 'b.')
% plot3(dirMatrix(1,:), dirMatrix(2,:), dirMatrix(3,:), 'b.')

%% relative roughness to track
rrVector = [0 0.05 0.1 0.2 0.4 0.6 0.8 0.9 0.95 1];

%% array of random background
% not used
rbxMatrix = rand(length(rrVector), size(dirMatrix, 2))*0;
rbyMatrix = rand(length(rrVector), size(dirMatrix, 2))*0;

%% sizes to track
sizeVector = [100 125 150 200 250 300 400 500];

%% options
caseName = 'localSqRough';
  cdef = obtainCrystal('generic_potash_alum', 3);
  hrep.A = cdef.A;
%   hrep.h = cdef.groupMappingMatrix * [1 1.45 1.15]';
  hrep.h = cdef.groupMappingMatrix * [1 1.5 1.2]';
%   hrep.h = cdef.groupMappingMatrix * [1 3 3]';
optionStruct.hrepBase = hrep;
optionStruct.writeBaseImage = false;
optionStruct.writeBlurrImage = false;
optionStruct.writeEdgeImage = false;
optionStruct.blurrSize = 50;
optionStruct.blurrStd = 12;
optionStruct.noise = 0.0002;
optionStruct.rrDiff = 0.1*0;

%% run
clear data
for iSize = 1:length(sizeVector)
    disp(['Processing size: ' num2str(sizeVector(iSize))]);
    optionStruct.particleSize = sizeVector(iSize);
    data(iSize) = runSimulation(dirMatrix, rrVector, {rbxMatrix, rbyMatrix}, optionStruct);
end

%% prepara data
[mwGrid rrGrid] = meshgrid(sizeVector, rrVector);
mwMeanMatrix = mwGrid*0;
mwStdMatrix = mwGrid*0;
rrMeanMatrix = mwGrid*0;
rrStdMatrix = mwGrid*0;
for iSize = 1:length(sizeVector)
    mwMeanMatrix(:, iSize) = data(iSize).mwMeanVector;
    mwStdMatrix(:, iSize) = data(iSize).mwStdVector;
    rrMeanMatrix(:, iSize) = data(iSize).rrMeanVector;
    rrStdMatrix(:, iSize) = data(iSize).rrStdVector;
end

figure,surf(mwGrid, rrGrid, rrStdMatrix)
figure,surf(mwGrid, rrGrid, rrMeanMatrix)
figure,surf(mwGrid, rrGrid, mwStdMatrix./mwGrid)


save(['data' filesep caseName '.mat'])
%% save
% save(['data' filesep caseName '.mat'], 'data');

%% eval
% figure,hold on
% plot([0 1], [0 1], 'k-', 'color', [1 1 1]*0.5, 'lineWidth', 5);
% plot(data.rrVector, data.rrMeanVector, 'r-', 'lineWidth', 3);
% plot(data.rrVector, data.rrMeanVector + data.rrStdVector, 'g-', 'lineWidth', 2);
% plot(data.rrVector, data.rrMeanVector - data.rrStdVector, 'g-', 'lineWidth', 2);
% plot(data.rrVector, data.rrMaxVector, 'm-', 'lineWidth', 1);
% plot(data.rrVector, data.rrMinVector, 'm-', 'lineWidth', 1);
% print(['data' filesep caseName '.eps'], '-depsc')
% 
% data.mwMatrix