gtMosaicityCmaps.m

function maps = gtMosaicityCmaps(dct_vol, maxgrain, spotInfo, pairInfo, etaCorr)
% maps=gtMosaicityCmaps(dct_vol, maxgrain, spotInfo, pairInfo, etaCorr)
% ---------------------------------------------------------------------
% function to read average omega extension per grain data from the
% spotpairs table and make colormaps.  
% make volume and radius cmaps,based on the volume passed in
%
% data(grainid,:)=[Nspots,mean(extent) med(extent) mean(ext-extent)
% med(ext-extent) volume radius]

maxImage = max(spotInfo.EndImage);

imspread = spotInfo.EndImage - spotInfo.StartImage;
extimspread = spotInfo.ExtEndImage - spotInfo.ExtStartImage;
if any(imspread < 0)
    imspread(imspread<0) = imspread(imspread<0) + maxImage;
end
if any(extimspread < 0)
    extimspread(extimspread<0) = extimspread(extimspread<0) + maxImage;
end

if ~exist('maxgrain','var') || isempty(maxgrain)
    maxgrain = max(dct_vol(:));
end

data = [];

for gID = 1:maxgrain
    listA = pairInfo.difAID(pairInfo.grainID==gID);
    listB = pairInfo.difBID(pairInfo.grainID==gID);
    eta   = pairInfo.eta(pairInfo.grainID==gID);
    theta = pairInfo.theta(pairInfo.grainID==gID);
    
    eta   = [eta; eta+180];
    theta = [theta; theta];
   
    list = [listA; listB];
    tmp    = [];
    exttmp = [];
    if (etaCorr)
        tan2eta=tand(eta).^2;
        sin2th=sind(theta).^2;
        ny=sqrt(tan2eta.*(1-sin2th)./(1+tan2eta));
        ny(isnan(ny))=1; % if eta=90,cor->NaN
    end
    for ii = 1:length(list)
        tmp(ii) = imspread(spotInfo.difspotID==list(ii));
        exttmp(ii) = extimspread(spotInfo.difspotID==list(ii));

        if (etaCorr)
            %sin(eta) correction for the eta angle...  is this reasonable?
            % yes,but better is
            tmp(ii)=tmp(ii)*ny(ii);
            exttmp(ii)=exttmp(ii)*ny(ii);
        end
    end

    data(gID,1)  = mean(tmp);
    data(gID,2)  = median(tmp);
    data(gID,3)  = mean(exttmp);
    data(gID,4)  = median(exttmp);
    n{gID}       = length(list);
    ids{gID}     = list';
    imsp{gID}    = tmp';
    extimsp{gID} = exttmp';
end

% define color map
maps.defineMap = @(x) [0 0 0; repmat(data(:,x),1,3)/max(data(:,x))];
% add a colour for grains with no data; x is the map
maps.addNaNColor = @(x,color) repmat(color,length(find(isnan(x(:,1)))),1);

% add data to the maps structure for future reference
maps.data     = data;
maps.data_key = 'nspots mean med meanext medext';% volume radius';
maps.n        = n;
maps.ids      = ids;
maps.imsp     = imsp;
maps.extimsp  = extimsp;

end