classdef Gt6DTwinTestCreator < Gt6DTestCreator
    methods (Access = public)
        function obj = Gt6DTwinTestCreator(volsNum)
            obj = obj@Gt6DTestCreator(volsNum);
        end

        function grain = generateTwinsGeometry(obj, parameters, R_vector, center)
        % To be cleaned!
            grain = [];
            grain.phaseid = 1;
            grain.R_vector = R_vector;
            pixelsize = mean([parameters.labgeo.pixelsizeu, parameters.labgeo.pixelsizev]);
            grain.center = center .* pixelsize;

            lp = parameters.cryst(grain.phaseid).latticepar;
            twin_axis = [1 1 1];
            twin_angle = 60;
            spacegroup = 225;

            twin_vars = gtTwinOrientations(grain.R_vector, twin_angle, twin_axis, spacegroup, lp, false);

            for jj = length(twin_vars):-1:1
                disp('CHANGE THE SIGN OF THE R-VECTOR')
                my_twin_vars(jj).R_vector = -twin_vars(jj).R_vector;
                my_twin_vars(jj).g = twin_vars(jj).g;
                my_twin_vars(jj).center = grain.center;
            end
            twin_vars = my_twin_vars;

            for jj = 1:length(twin_vars)
              grain.R_vector = -twin_vars(jj).R_vector;
              grain = gtCalculateGrain(grain, parameters);
              twin_vars(jj).allblobs = grain.allblobs;
            end

            grain.twin_vars = twin_vars; % this will be the parent, then all the twins

            % parent grain - reset the normal stuff in grain
            grain.R_vector = R_vector;
            grain = gtCalculateGrain(grain, parameters);

            grain.proj.geom = obj.buildGrainGeom(parameters, grain, pixelsize);
            for ii = 1:numel(grain.twin_vars)
                grain.twin_vars(ii).proj.geom = obj.buildGrainGeom(parameters, grain.twin_vars(ii), pixelsize);
            end

              % If we want to show the result of forward simulation
%             fig = figure();
%             ax = axes('Parent', fig);
%             hold(ax, 'on');
%             for ii = 1:size(grain.proj.geom, 1)
%                 vect = grain.proj.geom(ii, 1:3);
%                 quiver3(ax, 0, 0, 0, vect(1), vect(2), vect(3));
%             end
%             hold(ax, 'off');
%             grid(ax, 'on');
%             zlim(ax, [-1, 1]);
        end

        function geoms = matchParentTwinGeometries(obj, grain)
            % Now we need to find out which geometries are in common with
            % the parent
            numVols = numel(obj.singleVolumes);

            detectorDistance = grain.proj.geom(1, 4:6);
            detectorDistance = norm(detectorDistance);

            parentGeom = grain.proj.geom(:, 1:3);
            numCols = size(parentGeom, 1);
            geoms = zeros(numCols, numVols-1);

            for ii = 1:numVols-1
                twinGeom = grain.twin_vars(ii).proj.geom(:, 1:3);
                indexesTwin = ones(size(twinGeom, 1), 1);

                for jj = 1:numCols
                    diffGeoms = parentGeom(jj .* indexesTwin, :) - twinGeom;
                    normDiffs = sqrt(sum(diffGeoms .* diffGeoms, 2));
                    [val, indx] = min(normDiffs);
                    angle = 2 * asind( val / 2);
                    if (angle < 5e-1)
                        geoms(jj, ii) = indx;
                        fprintf('Found: norm %f, angle %f (deg), offset %f\n', ...
                            val, angle, val * detectorDistance)
                    end
                end
                indexes = find(geoms(:, ii) ~= 0);
                [parentGeom(indexes, :), twinGeom(geoms(indexes, ii), :)]
            end

        end

        function [indexes, finalGeoms] = createTwins(obj, metaInfoGrains, grains)
            if (length(grains) < 2)
                error('CS:wrong_parameter', ...
                      'At least two twin geometries are needed');
            end
            if (length(grains) > 2)
                error('CS:wrong_parameter', ...
                      'More than two is not supported, yet.');
            end
            indexes = [];
            [indexes.common_id, indexes.common_indx1, indexes.common_indx2] ...
                = intersect(grains{1}.difspotID(metaInfoGrains{1}.selectedDiffspots), ...
                            grains{2}.difspotID(metaInfoGrains{2}.selectedDiffspots));
            [indexes.only_1_id, indexes.only_1_indx] ...
                = setdiff(grains{1}.difspotID, grains{2}.difspotID);
            [indexes.only_2_id, indexes.only_2_indx] ...
                = setdiff(grains{2}.difspotID, grains{1}.difspotID);

            finalGeoms = cell(2, 1);
            finalGeoms{1} = grains{1}.proj.geom(indexes.common_indx1, :);
            finalGeoms{2} = grains{1}.proj.geom(indexes.common_indx1, :);

            modifiedOnly1 = grains{1}.proj.geom(indexes.only_1_indx, :);
            finalGeoms{1} = [finalGeoms{1}; modifiedOnly1];
            modifiedOnly1(:, 4:6) = obj.findPerpendicularDirection(modifiedOnly1(:, 4:6));
            finalGeoms{2} = [finalGeoms{2}; modifiedOnly1];

%             diffCenter = grains{2}.center - grains{1}.center;

            modifiedOnly2 = grains{2}.proj.geom(indexes.only_2_indx, :);
%             modifiedOnly2(:, 4:6) = modifiedOnly2(:, 4:6) ...
%                                     + diffCenter(ones(size(modifiedOnly2, 1), 1), :);
            finalGeoms{2} = [finalGeoms{2}; modifiedOnly2];
            modifiedOnly2(:, 4:6) = obj.findPerpendicularDirection(modifiedOnly2(:, 4:6));
            finalGeoms{1} = [finalGeoms{1}; modifiedOnly2];

            fprintf('Common spots (num %04d):\n', length(indexes.common_id))
            fprintf('  %d\n', indexes.common_id);
            fprintf('First volume''s spots (num %04d):\n', length(indexes.only_1_id))
            fprintf('  %d\n', indexes.only_1_id);
            fprintf('Second volume''s spots (num %04d):\n', length(indexes.only_2_id))
            fprintf('  %d\n', indexes.only_2_id);
        end

        function [detectorImages, trueDetectorImages] ...
                            = getSyntheticTwinDetectorImages(obj, grain, geoms, maxIndivSpots)
        % GETSYNTHETICTWINDETECTORIMAGES

%             [~, finalGeoms] = obj.createTwins(metaInfoGrains, grains);

            numVols = numel(obj.singleVolumes);
            finalGeoms = cell(1, numVols);

            if (numVols > 2)
                error('Gt6DTwinTestCreator:wrong_parameter', ...
                    'More than 2 volumes, not supported yet')
            end

            jj = 1;
            parentGeom = grain.proj.geom;
            twinGeom = grain.twin_vars(jj).proj.geom;

            commonIndx = find(geoms(:, jj) ~= 0);
            onlyParentIndx = setdiff(1:size(parentGeom, jj), commonIndx);
            onlyTwinIndx = setdiff(1:size(twinGeom, jj), geoms(commonIndx, jj));

            if (exist('maxIndivSpots', 'var'))
                if (numel(onlyParentIndx) > maxIndivSpots)
                    onlyParentIndx = onlyParentIndx(randperm(numel(onlyParentIndx), maxIndivSpots));
                end
                if (numel(onlyTwinIndx) > maxIndivSpots)
                    onlyTwinIndx = onlyTwinIndx(randperm(numel(onlyTwinIndx), maxIndivSpots));
                end
            end

            perpendParent = [ ...
                obj.findPerpendicularDirection(twinGeom(onlyTwinIndx, 1:3)), ...
                twinGeom(onlyTwinIndx, 4:end) ...
                ];
            finalGeoms{1} = [ ...
                parentGeom(onlyParentIndx, :); ...
                parentGeom(commonIndx, :); ...
                perpendParent ...
                ];

            perpendTwin = [ ...
                obj.findPerpendicularDirection(parentGeom(onlyParentIndx, 1:3)), ...
                parentGeom(onlyParentIndx, 4:end) ...
                ];
            finalGeoms{jj+1} = [ ...
                perpendTwin; ...
                parentGeom(commonIndx, :); ...
                twinGeom(onlyTwinIndx, :) ...
                ];

            for ii = 1:numel(finalGeoms)
                obj.proj.addProjectionGeometry( finalGeoms{ii} );
            end

            % Visualization of diffraction directions
            f = figure();
            ax = axes('parent', f);

            hold(ax, 'on');
            zeroCols = zeros(numel(onlyParentIndx), 1);
            quiver3(ax, zeroCols, zeroCols, zeroCols, ...
                parentGeom(onlyParentIndx, 1), parentGeom(onlyParentIndx, 2), parentGeom(onlyParentIndx, 3))

            zeroCols = zeros(numel(commonIndx), 1);
            quiver3(ax, zeroCols, zeroCols, zeroCols, ...
                parentGeom(commonIndx, 1), parentGeom(commonIndx, 2), parentGeom(commonIndx, 3))

            zeroCols = zeros(numel(onlyTwinIndx), 1);
            quiver3(ax, zeroCols, zeroCols, zeroCols, ...
                twinGeom(onlyTwinIndx, 1), twinGeom(onlyTwinIndx, 2), twinGeom(onlyTwinIndx, 3))
            hold(ax, 'off');

            grid(ax, 'on');
            zlim(ax, [-1 1]);

            [detectorImages, trueDetectorImages] = obj.getDetectorImages();
        end

        function [detectorImages, trueDetectorImages] ...
                            = getTrueTwinDetectorImages(obj, metaInfoGrains, grains)
        %
            numGrains = length(grains);
            volumeOversize = 1;
            spotsOversize = 1.2;

            if (numGrains < 2)
                error('CS:wrong_parameter', ...
                      'At least two twin geometries are needed');
            end
            if (numGrains > 2)
                error('CS:wrong_parameter', ...
                      'More than two is not supported, yet.');
            end

            % Let's create initial geometry, and decide which spots go in
            [indexes, finalGeoms] = obj.createTwins(metaInfoGrains, grains);

            % Shift geometries (relative to centers), create main voloume
            % and give it a decent size
            obj.fixCenterToGeometriesAndAdd(finalGeoms, grains, volumeOversize);

            % First crop/pad the bigger images
            for ii = 1:numGrains
                sizes(ii, :) = size(grains{ii}.proj.stack);
            end
            maxSpotBBox = [ max([sizes(:, 1); round(obj.proj.proj_size(1) * spotsOversize)]), ...
                            max([sizes(:, 3); round(obj.proj.proj_size(2) * spotsOversize)])];
            offsets(:, 1) = floor((maxSpotBBox(1) - sizes(:, 1)) / 2);
            offsets(:, 2) = zeros(numGrains, 1);
            offsets(:, 3) = floor((maxSpotBBox(2) - sizes(:, 3)) / 2);

            modifiedStacks = cell(numGrains, 1);
            for ii = 1:numGrains
                modifiedStacks{ii} = zeros( maxSpotBBox(1), ...
                                            size(grains{ii}.proj.stack, 2), ...
                                            maxSpotBBox(2));

                modifiedStacks{ii} = gtPlaceSubVolume(modifiedStacks{ii}, ...
                                                  grains{ii}.proj.stack, ...
                                                  offsets(ii, :), ...
                                                  0, 'assign');
            end

            detectorImages = [modifiedStacks{1}(:, indexes.common_indx1, :), ...
                              modifiedStacks{1}(:, indexes.only_1_indx, :), ...
                              modifiedStacks{2}(:, indexes.only_2_indx, :)];

            trueDetectorImages = permute(detectorImages, [1 3 2]);
        end
    end

    methods (Access = protected)
        function proj_geom = buildGrainGeom(obj, parameters, grain, pixelsize)
            u  = grain.allblobs.uv(:, 1);
            v  = grain.allblobs.uv(:, 2);

            spotsCommProps = getDiffractionSpotsCommonProperties(obj);

            radiuses = sqrt((u - parameters.labgeo.detrefu).^2 + (v - parameters.labgeo.detrefv).^2);

            %%% Get the indices of all spots which fully intersect the active
            % region of the detector
            onDet = Gt6DTwinTestCreator.findSpotsOnDetector(parameters, spotsCommProps, u, v, radiuses);

            rotComp = gtMathsRotationMatrixComp(parameters.labgeo.rotdir, 'row');
            getRotationTensor = @(x)gtMathsRotationTensor(-x, rotComp);

            detposcorner = gtGeoDetOrig(parameters.labgeo) ./ pixelsize;

            proj_geom = [];
            for ii = numel(onDet):-1:1
                currentReflection = onDet(ii);
                % use theoretically predited directions (default)
                omega = grain.allblobs.omega(currentReflection);
                omegaRotTens = getRotationTensor(omega);
                pl = grain.allblobs.pl(currentReflection, :);
                theta = grain.allblobs.theta(currentReflection);
                % theoretically predicted beam direction
                diffBeamDir = parameters.labgeo.beamdir * omegaRotTens + 2 * sind(theta) * pl;

                % position of Detector
%                 detPos = detposcorner ...
%                     + spotsCommProps.stackUSize/2 * parameters.labgeo.detdiru * omegaRotTens ...
%                     + spotsCommProps.stackVSize/2 * parameters.labgeo.detdirv * omegaRotTens;
%                 detPos = detPos * omegaRotTens - grain.center ./ pixelsize;
                detPos = diffBeamDir * norm(detposcorner);

                detdiru = parameters.labgeo.detdiru * omegaRotTens;
                detdirv = parameters.labgeo.detdirv * omegaRotTens;

                proj_geom(ii, :) = [diffBeamDir, detPos, detdiru, detdirv];
            end
        end

        function spotsProps = getDiffractionSpotsCommonProperties(obj)
            %%% Get the average difspot X and Y sizes for the current spots in the grain,
            spotsProps = [];
            spotsProps.Xsize = mean(size(obj.fullVol)); % average X Size of difspots
            spotsProps.Ysize = mean(size(obj.fullVol)); % average Y Size of difspots

            % define theshold values which will be used during the database search later
            % - limit to 5 pixels minimum width
            spotsProps.XsizeMin = max(5, spotsProps.Xsize - 1);
            spotsProps.XsizeMax = spotsProps.Xsize + 1;
            spotsProps.YsizeMin = max(5, spotsProps.Ysize - 1);
            spotsProps.YsizeMax = spotsProps.Ysize + 1;

            % allow for a lateral deviation of spot postions based on grain orientation statics
            % dif_props.Offset = fsim.Rdist_factor * gr.stat.Rdiststd * norm(detrefpos_pix)
            % !!! detrefpos is not necessarily the distance; depends on definition by user
            spotsProps.Offset = 1;

            diagonal = round(norm(size(obj.fullVol)));
            % take the biggest spot and add a relative extra margin (multiplicative factor oversize)
            spotsProps.stackUSize = diagonal;
            % ...same for vertical direction
            spotsProps.stackVSize = diagonal;
        end
    end

    methods (Access = protected, Static)
        function ondet = findSpotsOnDetector(parameters, dif_props, u, v, radius)
            %%% Get the indices of all spots which fully intersect the active region of the detector
            props = [ u > dif_props.Xsize/2, u < (parameters.acq.xdet - dif_props.Xsize/2), ...
                      v > dif_props.Ysize/2, v < (parameters.acq.ydet - dif_props.Ysize/2), ...
                      radius < parameters.acq.maxradius];
            ondet = find(all(props, 2));

            % Exclude spots which overlap with segmentation boundingbox
            seg = find(all([( (u + dif_props.Xsize/2 > parameters.seg.bbox(1)) ...
                            & (u - dif_props.Xsize/2 < parameters.seg.bbox(1) + parameters.seg.bbox(3))), ...
                            ( (v + dif_props.Ysize/2 > parameters.seg.bbox(2)) ...
                            & (v - dif_props.Ysize/2 < parameters.seg.bbox(2) + parameters.seg.bbox(4)))], 2));

            ondet = setdiff(ondet, seg); % indices of spots intersection active area of detector
        end
    end
end