text:
function [eddy] = thresholdBU(cyc, block_bottom_index, block_top_index, block_left_index, block_right_index, ...
        ssh, extrema, extrema_lat_index, extrema_lon_index, thresh, threshold_step, last_step, previous, ...
        lat, lon, R, areamap, min_pixel_size, is_padding, cyc_ssh)
% THRESHOLDBU Get an eddy by bottom up method
%   cyc: 1 for anticyclonic and -1 for cyclonic
%   block_bottom(top/left/right)_index: index of the bottom/top/left/right of the block that will be used for
%   thresholding
%   ssh: ssh data(extended if is_padding is true)
%   extrema: logical index of ssh extrema
%   extrema_lat/lon_index: lat/lon index of the extremum that is being used to find an eddy
%   thresh: current threshold that is being used to find an eddy
%   threshold_step: the step to increase/decrease threshold value, based on eddy type
%   last_step: the last step was used for thresholding
%   previous: 2d logical array of the connected component that contains the extremum in the last thresholdBU call
%   lat: 1d array of latitudes of the ssh grid
%   lon: 1d array of longitudes of the ssh grid
%   R: the georasterref object for SSH grid
%   areamap: 1D or 2D array of reference to area of each pixel in SSH grid
%   min_pixel_size: minimum number of pixels for an eddy
%   is_padding: whether or not the SSH data is padded

    switch cyc
        case 1
            intensity = 'MaxIntensity';
        case -1
            intensity = 'MinIntensity';
        otherwise
            error('Invalid cyc');
    end

    if block_bottom_index < 1 || block_top_index > size(ssh, 1) || ...
            block_left_index < 1 || block_right_index > size(ssh, 2)
        edgeOfWorld = true;

        while block_bottom_index < 1 || block_top_index > size(ssh, 1) || ...
            block_left_index < 1 || block_right_index > size(ssh, 2)
            % Make sure that the block is inside the grid

            block_bottom_index = block_bottom_index + 1;
            block_top_index = block_top_index - 1;
            block_left_index = block_left_index + 1;
            block_right_index = block_right_index - 1;
            if block_top_index <= block_bottom_index + 2 || block_right_index <= block_left_index + 2
                % If the block is too small, just return an empty eddy
                eddy = new_eddy();
                return;
            end
        end

        block = ssh(block_bottom_index:block_top_index, block_left_index:block_right_index);

        extremaBlock = extrema(block_bottom_index:block_top_index, block_left_index:block_right_index);

    else
        edgeOfWorld = false;
        block = ssh(block_bottom_index:block_top_index, block_left_index:block_right_index);
        extremaBlock = extrema(block_bottom_index:block_top_index, block_left_index:block_right_index);
    end

    if isnan(last_step)
        step = threshold_step;
    else
        step = last_step;
    end

    iter = 1;
    while true
        iter = iter+1;
        if iter > 5000

            perim = imdilate(logical(current), ones(3)) & ~logical(current);
            if all(isnan(block(perim)))
                eddy = new_eddy();
                return;
            end
            disp('potential infinite loop')
        end

        bw = cyc .* block >= cyc .* thresh;
        labels = bwlabel(bw);

        extrema_label = labels(extrema_lat_index - block_bottom_index + 1, extrema_lon_index - block_left_index + 1);
        current = labels == extrema_label;
        currentExtrema = extremaBlock(current);

        existing_pixel_at_box_edge = outterRing(labels, extrema_label);

        if sum(currentExtrema) > 1 || ( edgeOfWorld && existing_pixel_at_box_edge)

            if step ~= threshold_step
                % Go back to last threshold
                thresh = thresh + cyc*step;
                step = threshold_step;
                thresh = thresh - cyc*step;
                continue;
            end

            if size(block, 1) ~= size(previous(block_bottom_index:block_top_index, block_left_index:block_right_index), 1)
                prevBlock = block(2:end-1, 2:end-1);
            else
                prevBlock = block;
            end

            perim = imdilate(logical(previous(block_bottom_index:block_top_index, block_left_index:block_right_index)), ones(3)) ...
                & ~logical(previous(block_bottom_index:block_top_index, block_left_index:block_right_index));
            nan = isnan(prevBlock(perim));
            if sum(nan) / length(nan) > .3
                %if more than half of your perimeter is land, then throw it out.
                eddy = new_eddy();
                return;
            end

            if sum(previous(:)) < min_pixel_size
                    eddy = new_eddy();
                    return;
            end

            perim = bwperim(previous);
            meanPerim = mean(ssh(logical(perim)));
            amp = cyc * (ssh(extrema_lat_index, extrema_lon_index)-meanPerim);

            stats = regionprops(previous, ssh, 'Area', 'Extrema',...
                'PixelIdxList', intensity, 'ConvexImage', 'PixelList', ...
                'Solidity', 'Extent', 'Orientation', 'MajorAxisLength', ...
                'MinorAxisLength');

            stats.Intensity = stats.(intensity);
            stats = rmfield(stats, intensity);

            if is_padding
                [idx, r, c] = extidx2original(stats.PixelIdxList, [length(lat) length(lon)], size(ssh));
                stats.PixelIdxList = idx;
            else
                [r, c] = ind2sub(size(ssh), stats.PixelIdxList);
            end

            stats.PixelList = [r, c];

            % Getting geodesic speed
            if is_padding
                geoSpeed = mean_geo_speed(ssh(:, 201:end-200), stats.PixelIdxList, lat, lon);
                if ~isempty(R)
                    [elat, elon] = weighted_centroid(cyc_ssh(:, 201:end-200), stats.PixelList, stats.PixelIdxList, R);
                else
                    [elat, elon] = weighted_centroid_irregular_grid(cyc_ssh(:, 201:end-200), stats.PixelList, stats.PixelIdxList, lat, lon);
                end
            else
                geoSpeed = mean_geo_speed(ssh, stats.PixelIdxList, lat, lon);
                if ~isempty(R)
                    [elat, elon] = weighted_centroid(cyc_ssh, stats.PixelList, stats.PixelIdxList, R);
                else
                    [elat, elon] = weighted_centroid_irregular_grid(cyc_ssh, stats.PixelList, stats.PixelIdxList, lat, lon);
                end
            end

            % weighted_centroid returns lon from 0-360, fix this
            % TODO: should we also fix lat lon -270 to 80?
            elon = (elon > 180).*(elon - 360) + (elon <= 180).*elon;

            % Getting surface area of the eddy
            if all(size(areamap) == [length(lat) length(lon)])
                % area is 2D array for areas of pixels at [lat, lon]
                sarea = sum(areamap(stats.PixelIdxList));
            elseif any(size(areamap) == [1 1]) && length(areamap) == length(lat)
                % Area is 1D array for areas of pixels at a specific latitude
                sarea = sum(areamap(stats.PixelList(:, 1)));
            else
                % Invalid areamap
                sarea = NaN;
            end

            eddy = new_eddy(rmfield(stats, 'PixelList'), amp, elat, elon, thresh, sarea, cyc, geoSpeed, 'ESv2');

            return
        end

        if existing_pixel_at_box_edge
            %disp('expanding size');
            eddy = thresholdBU(cyc, block_bottom_index-1, block_top_index+1, block_left_index-1,...
                block_right_index+1, ssh, extrema, extrema_lat_index, extrema_lon_index, thresh, threshold_step, ...
                step, previous, lat, lon, R, areamap, min_pixel_size, ...
                is_padding, cyc_ssh);
            return
        end

        previous(block_bottom_index:block_top_index, block_left_index:block_right_index) = current;

        step = step * 2; % double the step for less number of iterations
        thresh = thresh - cyc*step;
    end

end

function [res] = outterRing(box, val)
    ring = [box(1, 1:end)'; box(end, 1:end)'; box(1:end, 1); box(1:end, end)];
    res = any(ring == val);
end

author: jfaghm

repo: OceanEddies

save_path: github-repos/MATLAB/jfaghm-OceanEddies


text:
function select_segments_on_base_and_followup_trees(baseCerrFile,baseTreeFile,...
    followupTreeFile,vfFile,segmentsFile)
% function  select_segments_on_base_and_followup_trees(baseCerrFile,baseTreeFile,...
%     followupTreeFile,vfFile)
%
% Example
%
% APA, 5/4/2021

global airwayStateS stateS
% airwayStateS.baseAddNodes = 0;
% airwayStateS.baseShowInCerrViewer = 1;
% airwayStateS.baseRemoveNodes = 0;
% airwayStateS.followupAddNodes = 0;
% airwayStateS.followupShowInCerrViewer = 1;
% airwayStateS.followupRemoveNodes = 0;
% airwayStateS.baseSegmentStartNodes = [];
% airwayStateS.baseSegmentEndNodes = [];
% airwayStateS.basehStartNodes = [];
% airwayStateS.basehEndNodes = [];
% airwayStateS.followupSegmentStartNodes = [];
% airwayStateS.followupSegmentEndNodes = [];
% airwayStateS.followuphStartNodes = [];
% airwayStateS.followuphEndNodes = [];
% airwayStateS.hBaseSegment = [];
% airwayStateS.hFollowupSegment = [];
airwayStateS.selectBaseStartNode = 1;
airwayStateS.selectFollowupStartNode = 1;
% airwayStateS.hBaseStart
% airwayStateS.hBaseEnd
% airwayStateS.hFollowupStart
% airwayStateS.hFollowupEnd

% Create crosshairs
numAxis = length(stateS.handle.CERRAxis);
for i = 1:numAxis
    hCrosshairV(i) = plot(0,0, ...
        'parent', stateS.handle.CERRAxis(i), 'Marker','+', 'MarkerSize',10,...
        'color', [1 1 0], 'hittest', 'off','linewidth',3,'visible','off');
end

airwayStateS.hCrosshairV = hCrosshairV;

% Load base scan with dose
planC = loadPlanC(baseCerrFile,tempdir);
planC = updatePlanFields(planC);
planC = quality_assure_planC(baseCerrFile,planC);

% Get grid for vf
scanNum = 1;
indexS = planC{end};
[xUnifV,yUnifV,zUnifV] = getUniformScanXYZVals(planC{indexS.scan}(scanNum));
xFieldV = [xUnifV(1),xUnifV(2)-xUnifV(1),xUnifV(end)];
yFieldV = [yUnifV(end),yUnifV(1)-yUnifV(2),yUnifV(1)];

% % Load base tree
% load(baseTreeFile)
% elemDistBaseV = elemDistV;
% elemAreaBaseV = elemAreaV;
% elemBaseM = elemM;
% nodeXyzBaseM = nodeXyzM;

% Load base tree
load(baseTreeFile)
minDistBaseV = minDistV;
elemBaseM = elemM;
nodeXyzBaseM = nodeXyzM;
centerTreeXyzBaseM = centerTreeXyzM;
% nodeDistFromCarinaBaseV = nodeDistFromCarinaV;
treeS.nodeXyzM = nodeXyzM;
treeS.elemM = elemM;
xyzStartM = [nodeXyzM(elemM(:,1),1),nodeXyzM(elemM(:,1),2),nodeXyzM(elemM(:,1),3)];
xyzEndM = [nodeXyzM(elemM(:,2),1),nodeXyzM(elemM(:,2),2),nodeXyzM(elemM(:,2),3)];
elemDistV = sum((xyzStartM - xyzEndM).^2,2).^0.5;
airwayGraph = graph(elemM(:,1),elemM(:,2),elemDistV);
treeS.airwayGraph = airwayGraph;
treeS.minDistV = minDistV;
airwayStateS.baseTreeS = treeS;

% Load followup tree
load(followupTreeFile)
minDistFollowupV = minDistV;
elemFollowupM = elemM;
nodeXyzFollowupM = nodeXyzM;
nodeXyzFollowupSmoothM = centerTreeXyzM;
treeS.nodeXyzM = nodeXyzM;
treeS.elemM = elemM;
xyzStartM = [nodeXyzM(elemM(:,1),1),nodeXyzM(elemM(:,1),2),nodeXyzM(elemM(:,1),3)];
xyzEndM = [nodeXyzM(elemM(:,2),1),nodeXyzM(elemM(:,2),2),nodeXyzM(elemM(:,2),3)];
elemDistV = sum((xyzStartM - xyzEndM).^2,2).^0.5;
airwayGraph = graph(elemM(:,1),elemM(:,2),elemDistV);
treeS.airwayGraph = airwayGraph;
% Smooth radius calculation
minDistBaseV = smoothRadius(nodeXyzBaseM,centerTreeXyzBaseM,minDistBaseV);
minDistFollowupV = smoothRadius(nodeXyzFollowupM,nodeXyzFollowupM,minDistFollowupV);
treeS.minDistV = minDistFollowupV;
airwayStateS.followupTreeS = treeS;

% Get deformation at nodes
load(vfFile)

% Get location of base points on followup scan
xDeformV = finterp3(nodeXyzBaseM(:,1),nodeXyzBaseM(:,2),nodeXyzBaseM(:,3),...
    flip(vf(:,:,:,1),1),xFieldV,yFieldV,zUnifV);
yDeformV = finterp3(nodeXyzBaseM(:,1),nodeXyzBaseM(:,2),nodeXyzBaseM(:,3),...
    flip(vf(:,:,:,2),1),xFieldV,yFieldV,zUnifV);
zDeformV = finterp3(nodeXyzBaseM(:,1),nodeXyzBaseM(:,2),nodeXyzBaseM(:,3),...
    flip(vf(:,:,:,3),1),xFieldV,yFieldV,zUnifV);

nodeXyzInterpM = nodeXyzBaseM + [xDeformV(:), yDeformV(:), zDeformV(:)];

airwayStateS.nodeXyzInterpM = nodeXyzInterpM;

% Calculate radius change
distTolerance = 1; %cm
radiusDiffV = calculateRadiusChange(nodeXyzBaseM,nodeXyzFollowupM,...
    minDistBaseV,minDistFollowupV,distTolerance,...
    vf,xFieldV,yFieldV,zUnifV);

% Calculate dose
doseNum = 1;
% nodeDoseV = getDoseAt(doseNum, centerTreeXyzBaseM(:,1), centerTreeXyzBaseM(:,2),...
%     centerTreeXyzBaseM(:,3), planC);
nodeOrigDoseV = getDoseAt(doseNum, nodeXyzBaseM(:,1), nodeXyzBaseM(:,2),...
    nodeXyzBaseM(:,3), planC);

airwayStateS.baseTreeS.nodeOrigDoseV = nodeOrigDoseV;

% Initialize figure to show airway trees
hFig = figure('CloseRequestFcn','airwaySegmentCallback(''closeRequest'')',...
    'numbertitle','off');
% baseAxis = subplot(1,2,1);
% followupAxis = subplot(1,2,2);
% set(baseAxis,'NextPlot','add')
% set(followupAxis,'NextPlot','add')

baseAxis = axes('parent', hFig, 'units', 'normalized','NextPlot','add', ...
    'position', [0.2 0.1 0.38 0.8], 'color', [1,1,1]);
followupAxis = axes('parent', hFig, 'units', 'normalized','NextPlot','add', ...
    'position', [0.6 0.1 0.38 0.8], 'color', [1,1,1]);
hSegList = uicontrol(hFig,'units','normalized','Position',[0.02 0.55, 0.15, 0.4],...
    'Style','list','callback','airwaySegmentCallback(''SEGMENT_CLICKED'')');

hSegAdd = uicontrol(hFig,'units','normalized','Position',[0.02 0.5, 0.03, 0.05],...
    'Style','push','String','+','callback','airwaySegmentCallback(''ADD_SEGMENT'')');
hSegRemove = uicontrol(hFig,'units','normalized','Position',[0.06 0.5, 0.03, 0.05],...
    'Style','push','String','-','callback','airwaySegmentCallback(''REMOVE_SEGMENT'')');
hSegRename = uicontrol(hFig,'units','normalized','Position',[0.1 0.5, 0.07, 0.05],...
    'Style','push','String','Name','callback','airwaySegmentCallback(''RENAME_SEGMENT'')');

hLoadSegments = uicontrol(hFig,'units','normalized','Position',[0.02 0.4, 0.15, 0.05],...
    'Style','push','String','Load','callback','airwaySegmentCallback(''LOAD_SEGMENTS'')');
hSaveSegments = uicontrol(hFig,'units','normalized','Position',[0.02 0.33, 0.15, 0.05],...
    'Style','push','String','Save','callback','airwaySegmentCallback(''SAVE_SEGMENTS'')');
hSegViewer = uicontrol(hFig,'units','normalized','Position',[0.02 0.26, 0.15, 0.05],...
    'Style','push','String','Viewer','callback','airwaySegmentCallback(''SHOW_IN_CERR_VIEWER'')');
hSelectSegNodes = uicontrol(hFig,'units','normalized','Position',[0.02 0.19, 0.15, 0.05],...
    'Style','push','String','Select seg nodes','callback','airwaySegmentCallback(''SELECT_SEGMENT_NODES'')');
hLockBaseTreeToggle = uicontrol(hFig,'units','normalized','Position',[0.02 0.12, 0.15, 0.05],...
    'Style','toggle','String','Lock base tree','callback','airwaySegmentCallback(''LOCK_BASE_TREE'')');
airwayStateS.hAxisBase = baseAxis;
airwayStateS.hAxisFollowup = followupAxis;
airwayStateS.hSegList = hSegList;
airwayStateS.hSegAdd = hSegAdd;
airwayStateS.hSegRemove = hSegRemove;
airwayStateS.hSegRename = hSegRename;
airwayStateS.hSegViewer = hSegViewer;
airwayStateS.hLoadSegments = hLoadSegments;
airwayStateS.hSaveSegments = hSaveSegments;
airwayStateS.hSelectSegNodes = hSelectSegNodes;
airwayStateS.hLockBaseTreeToggle = hLockBaseTreeToggle;
airwayStateS.viewerMode = 1;
airwayStateS.segSelectMode = 0;
segS = struct('segName','segment','nodeXyzM',[],...
    'startNode',[],'endNode',[],'allNodes',[],...
    'hPlot',[],'assocScanUID','');
airwayStateS.followupSegmentS = segS;
airwayStateS.followupSegmentS(:) = [];
airwayStateS.baseSegmentS = segS;
airwayStateS.baseSegmentS(:) = [];
airwayStateS.totalSegments = 0;
airwayStateS.currentSegment = 0;

basePt = plot3(baseAxis,0,0,0,'marker','o','markerSize',12,...
    'MarkerEdgeColor',[0,0,0],'visible','off','hittest','off');
airwayStateS.hBaseStart = plot3(baseAxis,0,0,0,'marker','o','markerSize',12,...
    'MarkerEdgeColor','g','visible','off','color','g','hittest','off');
airwayStateS.hBaseEnd = plot3(baseAxis,0,0,0,'marker','o','markerSize',12,...
    'MarkerEdgeColor','r','visible','off','color','r','hittest','off');

set(baseAxis,'zDir','reverse')
hPlotBase = scatter3(baseAxis,nodeXyzBaseM(:,1),nodeXyzBaseM(:,2),nodeXyzBaseM(:,3),...
'sizeData',12,'MarkerEdgeColor','None',...
'MarkerFaceColor','flat','MarkerFaceAlpha',0.5,'CData',[0.5,0.5,0.5]);
airwayStateS.hPlotBase = hPlotBase;
airwayStateS.basePt = basePt;
%set(hp,'ButtonDownFcn',{buttonDownFcn,varargin{:},nodeValV});
axis(baseAxis,'equal')
view(baseAxis,3)
xlabel(baseAxis,'R - L','fontsize',10)
ylabel(baseAxis,'A - P','fontsize',10)
zlabel(baseAxis,'I - S','fontsize',10)
grid(baseAxis,'on')
title(baseAxis,'Baseline')

% Load followup tree
% load(followupTreeFile)
% baseFig = figure;
% hold on,
% baseAxis = gca;
followupPt = plot3(followupAxis,0,0,0,'marker','o','markerSize',12,...
    'MarkerEdgeColor',[0,0,0],'visible','off','hittest','off');
airwayStateS.hFollowupStart = plot3(followupAxis,0,0,0,'marker','o','markerSize',12,...
    'MarkerEdgeColor','g','visible','off','color','g','hittest','off');
airwayStateS.hFollowupEnd = plot3(followupAxis,0,0,0,'marker','o','markerSize',12,...
    'MarkerEdgeColor','r','visible','off','color','r','hittest','off');
airwayStateS.followupPt = followupPt;
set(followupAxis,'zDir','reverse')
hPlotFollowup = scatter3(followupAxis,nodeXyzFollowupM(:,1),...
    nodeXyzFollowupM(:,2),nodeXyzFollowupM(:,3),...
'sizeData',12,'MarkerEdgeColor','None',...
'MarkerFaceColor','flat','MarkerFaceAlpha',0.5,'CData',[0.5,0.5,0.5]);
airwayStateS.hPlotFollowup = hPlotFollowup;

axis(followupAxis,'equal')
view(followupAxis,3)
xlabel(followupAxis,'R - L','fontsize',10)
ylabel(followupAxis,'A - P','fontsize',10)
zlabel(followupAxis,'I - S','fontsize',10)
grid(followupAxis,'on')
title(followupAxis,'Followup')

airwaySegmentCallback('SHOW_IN_CERR_VIEWER')

if exist('segmentsFile','var') && ~isempty(segmentsFile)
    airwaySegmentCallback('LOAD_SEGMENTS',segmentsFile);
end

% % Callback to add nodes
% buttonDownFcn = @pickStatrStopNodes;
% set(hPlotFollowup,'ButtonDownFcn',{buttonDownFcn,followupAxis});

% % Callback to show nodes in Viewer
% buttonDownFcn = @showNodeInCerrViewer;
% set(hPlotBase,'ButtonDownFcn',{buttonDownFcn,vf,xFieldV,yFieldV,zUnifV,...
%     basePt,followupPt})

% % Callback to remove nodes
% buttonDownFcn = @removeStatrStopNodes;
% set(hPlotFollowup,'ButtonDownFcn',{buttonDownFcn,followupAxis});

% Set buttondown function

%% Plots

% Create contextmenu
% baseMenu = uicontextmenu('Callback', 'airwayAxisMenu(''init_start_stop_view_nodes'')',...
%     'userdata', {'base',baseAxis,hPlotBase,basePt,vf,xFieldV,yFieldV,zUnifV,followupPt,...
%     minDistBaseV, radiusDiffV, nodeOrigDoseV}, 'parent', hFig);
baseMenu = uicontextmenu('Callback', 'airwayAxisMenu(''init_start_stop_view_nodes'')',...
    'userdata', {'base',baseAxis,hPlotBase,basePt,nodeXyzInterpM,followupPt,...
    minDistBaseV, radiusDiffV, nodeOrigDoseV}, 'parent', hFig);
% followupMenu = uicontextmenu('Callback', 'airwayAxisMenu(''init_start_stop_view_nodes'')',...
%     'userdata', {'followup',followupAxis,hPlotFollowup,followupPt}, 'parent', hFig);
set(baseAxis, 'UIContextMenu', baseMenu);
%set(followupAxis, 'UIContextMenu', followupMenu);

% buttonDownFcn = @showNodeInCerrViewer;
% set(hPlotBase,'ButtonDownFcn',{buttonDownFcn,'base',basePt,nodeXyzInterpM,...
%    followupPt})
% set(hPlotFollowup,'ButtonDownFcn',{buttonDownFcn,'followup',followupPt})

% followupMenu = uicontextmenu('Callback', 'airwayAxisMenu(''choose_start_stop'')',...
%     'userdata', {}, 'parent', followupFig);
% set(baseAxis, 'UIContextMenu', followupMenu);

% m1 = uimenu(baseCtxMenu,'Text','Airway Radius');
% m2 = uimenu(baseCtxMenu,'Text','RT Dose');
% m3 = uimenu(baseCtxMenu,'Text','Change in Airway Radius');

% % Plot elements color-coded by opening radius
% radiusDiffV = (elemDistV(minIndV)-elemDistBaseV)./elemDistBaseV*100;
% radiusDiffV(minDistV > 0.5) = NaN;
% minDist = min(radiusDiffV);
% maxDist = max(radiusDiffV);
% maxMinusMinDist = maxDist - minDist;
% cmapM = CERRColorMap('jetmod');
% cmapSiz = size(cmapM,1)-1;
% cmapIndV = round((radiusDiffV - minDist) / maxMinusMinDist * cmapSiz) + 1;
%
% xyzStartM = [nodeXyzBaseM(elemBaseM(:,1),1),nodeXyzBaseM(elemBaseM(:,1),2),nodeXyzBaseM(elemBaseM(:,1),3)];
% xyzEndM = [nodeXyzBaseM(elemBaseM(:,2),1),nodeXyzBaseM(elemBaseM(:,2),2),nodeXyzBaseM(elemBaseM(:,2),3)];
%
% figure ,hold on,
% set(gca,'zDir','reverse')
% for i = 1:size(xyzStartM,1)
%     if ~isnan(elemDistV(minIndV(i))) && ~isnan(elemDistBaseV(i)) && ~isnan(cmapIndV(i))
%         xV = [xyzStartM(i,1), xyzEndM(i,1)];
%         yV = [xyzStartM(i,2), xyzEndM(i,2)];
%         zV = [xyzStartM(i,3), xyzEndM(i,3)];
%         plot3(xV,yV,zV,'color',[cmapM(cmapIndV(i),:),1],...
%             'linewidth',3,'ButtonDownFcn',@showNodeOnFollowupTree); %@gotoNodeInCerrViewer
%         %plot3(xV,yV,zV,'b.')
%     end
% end
% set(gca,'cLim',[minDist,maxDist])
% colormap(cmapM)
% ticksV = linspace(minDist,maxDist,10);
% tickC = {};
% for i = 1:length(ticksV)
%     tickC{i} = num2str(ticksV(i));
% end
% hCbar = colorbar('Ticks',ticksV,'TickLabels',tickC);
% axis('equal')
% view(3)
% xlabel('x','fontsize',18)
% ylabel('y','fontsize',18)
% zlabel('z','fontsize',18)
% grid('on')
% title(titleStr)

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
filename='Bridge_hexahedra_coarse';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance';'perimeter'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
optimizer = 'MMA'; incrementFactor = 1;
filterType = 'P1';

nsteps = 5;
Vfrac_final = 0.05;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-2;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function api = WizardProgressApi( api, firstPage )
%IMPORTGUIPROGRESS Summary of this function goes here
%   Detailed explanation goes here

state = {'off','on'};

gui = localCreateGUIControls();

% local progress api
progressApi.currentPage = firstPage;
progressApi.prevPage = {};
progressApi.leavePageCallback  = [];
progressApi.gotoPageCallback   = [];

% global available functions
api.Progress.enableNext   = @enableNext;
api.Progress.enablePrev   = @enablePrev;
api.Progress.enableFinish = @enableFinish;
api.Progress.enablePlot   = @enablePlot;

% some helper
api.initProgressApi       = @showPage;
% api.updateProgressApi     = @updateProgressApi;
api.setLeavePageCallback  = @setLeavePageCallback;
api.setGotoPageCallback   = @setGotoPageCallback;

api.nextCallback  = @localNextCallback;

% set callbacks
set(gui.hPlot   ,'CallBack',@localPlotDataCallback );
set(gui.hPrev   ,'Callback',@localPrevCallback   );
set(gui.hNext   ,'Callback',@localNextCallback   );
set(gui.hFinish ,'CallBack',@localFinishCallback   );

% initalize buttons
api.Progress.enableNext(false);
api.Progress.enablePrev(false);
api.Progress.enablePlot(false);
api.Progress.enableFinish(false);

  function showPage()

    api.clearPage();

    set(api.hPanel,'visible','off');

    progressApi.currentPage(api);
    try
      progressApi.gotoPageCallback();
    catch
    end
    set(api.hPanel,'visible','on');

  end

  function setLeavePageCallback(callback)

    progressApi.leavePageCallback = callback;

  end

  function setGotoPageCallback(callback)

    progressApi.gotoPageCallback = callback;

  end

  function localPrevCallback(varargin)

    if numel(progressApi.prevPage) > 0

      progressApi.currentPage = progressApi.prevPage{end};
      progressApi.prevPage(end) = [];

      enablePrev(numel(progressApi.prevPage)>0);

      showPage();

    end

  end

  function localNextCallback(varargin)

    nextPage = progressApi.leavePageCallback();

    if ~isempty(nextPage)

      progressApi.prevPage{end+1} = progressApi.currentPage;
      progressApi.currentPage = nextPage;

      enablePrev(numel(progressApi.prevPage)>0);

      showPage();

    end

  end

  function localPlotDataCallback(varargin)

    data = api.getDataTransformed();
    figure

    if isa(data,'EBSD') && strcmp('pageSS',char(progressApi.currentPage)) % && numel(data.indexedPhasesId) == 1

      counts = accumarray(data.phaseId,1);
      [~,maxPhaseId] = max(counts(2:end));
      ind = data.phaseId == (1+maxPhaseId);
      ori = data(ind).orientations;
      plot(data(ind),ori);
      set(gcf,'name',['IPF Z map of phase ' data.CSList{maxPhaseId+1}.mineral])

      figure
      h = ori.CS.basicHKL;
      plotPDF(ori,h(1:3),'contourf','antipodal');
      set(gcf,'name',['pole figures of phase ' data.CSList{maxPhaseId+1}.mineral])

    else

      plot(data,'silent');
    end

  end

  function localFinishCallback(varargin)

    e = WizardFinish(api);

    if ~e, close(api.hFigure); end

  end

  function gui = localCreateGUIControls()

    parent = api.hFigure;
    w    = api.Spacings.Width;
    m    = api.Spacings.Margin;

    bW   = api.Spacings.ButtonWidth;
    bH   = api.Spacings.ButtonHeight;
    pH   = api.Spacings.BottomHeight;
    fs   = api.Spacings.FontSize;

    pos = @(offset) [offset pH-bH-m bW bH];

    % navigation
    next =  uicontrol(...
      'Parent',parent,...
      'String','Next >>',...
      'FontSize',fs,...
      'Position',pos(w-2*bW-2*m));

    prev = uicontrol(...
      'Parent',parent,...
      'String','<< Previous',...
      'FontSize',fs,...
      'Position',pos(w-3*bW-2*m));

    finish = uicontrol(...
      'Parent',parent,...
      'String','Finish',...
      'Enable','off',...
      'FontSize',fs,...
      'Position',pos(w-bW-m));

    plot = uicontrol(...
      'Parent',parent,...
      'String','Plot',...
      'FontSize',fs,...
      'Position',pos(m),...
      'Visible','on');

    uipanel(...
      'BorderType','line',...
      'BorderWidth',1,...
      'units','pixel',...
      'parent',parent,...
      'FontSize',fs,...
      'HighlightColor',[0 0 0],...
      'Position',[m pH w-2*m 1]);

    gui.hNext    = next;
    gui.hPrev    = prev;
    gui.hFinish  = finish;
    gui.hPlot    = plot;

  end

  function enableNext(b)

    set(gui.hNext,'enable',state{1+b});

  end

  function enablePrev(b)

    set(gui.hPrev,'enable',state{1+b});

  end

  function enableFinish(b)

    set(gui.hFinish,'enable',state{1+b});

  end

  function enablePlot(b)

    set(gui.hPlot,'enable',state{1+b});

  end

end

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
%  myConstantParamFilter - filtering procedure for state-space models with
%  constant parameters
%
%  ::
%
%
%    [LogLik,Incr,retcode,Filters]=myConstantParamFilter(syst,y,U,z,options)
%
%  Args:
%
%     - **syst** [struct]: structure provided by dsge.filter
%
%     - **y** [matrix]: matrix of data provided by dsge.filter
%
%     - **U** [matrix]: matrix of trends provided by dsge.filter
%
%     - **z** [matrix]: matrix of deterministic terms provided by dsge.filter
%
%     - **options** [struct]: options provided by dsge.filter
%
%  Returns:
%     :
%
%     - **LogLik** [numeric]: value of the log likelihood
%
%     - **Incr** [vector]: contributions to the likelihood in each period
%
%     - **retcode** [numeric]: flag equal to 0 if there is no problem
%
%     - **Filters** [struct]: structure containing all the filtering
%       information
%
%  Note:
%
%     - If the function is passed through a rise/dsge object, then it should
%       be called as ff=filter(m,'kf_user_algo',@myConstantParamFilter).
%
%     See also: myKnownRegimesFilter
%

author: jmaih


text:
function [b, mpe] = back_mpe(engine, bfuture, f, ev2, t)

if f.t ~= t
  error('mixed up time stamps')
end
if t==1
  error('should call back1_mpe')
end

maximize = 1;
bnet = bnet_from_engine(engine);
ss = bnet.nnodes_per_slice;

int = engine.interface;
D = engine.in_clq;
C = engine.out_clq;
phiD = marginalize_pot(bfuture.clpot{D}, int, maximize);
phiD = set_domain_pot(phiD, int+ss); % shift to slice 2
phiC = marginalize_pot(f.clpot{C}, int+ss, maximize);
ratio = divide_by_pot(phiD, phiC);
f.clpot{C} = multiply_by_pot(f.clpot{C}, ratio);

[mpe, b.clpot] = find_max_config(engine.jtree_engine, f.clpot, f.seppot, ev2);
mpe = mpe((1:ss)+ss); % extract values for slice 2
for c=1:length(b.clpot)
  [b.clpot{c}, ll(c)] = normalize_pot(b.clpot{c});
end
b.t = t;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [mellin,beta]=fmt(X,fmin,fmax,N);
%FMT    Fast Fourier Mellin Transform.
%       [MELLIN,BETA]=FMT(X,FMIN,FMAX,N) computes the Fast Mellin
%       Transform of signal X.
%
%       X : signal in time (Nx=length(X)).
%       FMIN,FMAX : respectively lower and upper frequency bounds of
%        the analyzed signal. These parameters fix the equivalent
%        frequency bandwidth (expressed in Hz). When unspecified, you
%        have to enter them at the command line from the plot of the
%        spectrum. FMIN and FMAX must be >0 and <=0.5.
%       N : number of analyzed voices. N must be even
%				 (default : automatically determined).
%       MELLIN : the N-points Mellin transform of signal S.
%       BETA : the N-points Mellin variable.
%
%       Examples :
%        sig=altes(128,0.05,0.45);
%	 [MELLIN,BETA]=fmt(sig,0.05,0.5,128);
%        plot(BETA,real(MELLIN));
%
%       See also IFMT, FFT, IFFT.

%       P. Goncalves 9-95 - O. Lemoine, June 1996.
%       Copyright (c) 1995 Rice University - CNRS (France) 1996.
%
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2 of the License, or
%  (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

if (nargin == 0),
 error('At least one parameter required');
end;

[xrow,xcol] = size(X);

if (nargin==2),
 disp('FMIN will not be taken into account. Determine it with FMAX');
 disp('     from the following plot of the spectrum.');
elseif nargin==3,
 N=[];
elseif (nargin==4 & rem(N,2)~=0),
 error('N must be even');
end;
if (xcol==0)|(xcol>2),
 error('X must have one or two columns');
end

Mt = length(X);
Z  = hilbert(real(X));
M  = (Mt+rem(Mt,2))/2;

if nargin<=2,                                   % fmin,fmax,N unspecified
 STF = fft(fftshift(X)); Nstf=length(STF);
 sp = (abs(STF(1:Nstf/2))).^2; Maxsp=max(sp);
 f = linspace(0,0.5,Nstf/2+1) ; f = f(1:Nstf/2);
 plot(f,sp) ; grid;
 xlabel('Normalized frequency');
 title('Analyzed signal energy spectrum');
 indmin=min(find(sp>Maxsp/1000));
 indmax=max(find(sp>Maxsp/1000));
 fmindflt=max([0.001 0.05*fix(f(indmin)/0.05)]);
 fmaxdflt=0.05*ceil(f(indmax)/0.05);
 txtmin=['Lower frequency bound [',num2str(fmindflt),'] : '];
 txtmax=['Upper frequency bound [',num2str(fmaxdflt),'] : '];
 fmin = input(txtmin); fmax = input(txtmax);
 if fmin==[], fmin=fmindflt; end
 if fmax==[], fmax=fmaxdflt; end
end

if fmin >= fmax
 error('FMAX must be greater or equal to FMIN');
elseif fmin<=0.0 | fmin>0.5,
 error('FMIN must be > 0 and <= 0.5');
elseif fmax<=0.0 | fmax>0.5,
 error('FMAX must be > 0 and <= 0.5');
end

B = fmax-fmin;       		% Bandwidth of the signal X
R = B/((fmin+fmax)/2);		% Relative bandwidth of X

Nq= ceil((B*Mt*(1+2/R)*log((1+R/2)/(1-R/2)))/2);
Nmin = Nq-rem(Nq,2);
Ndflt = 2^nextpow2(Nmin);
if nargin<=2,
 Ntxt=['Number of frequency samples (>=',num2str(Nmin),') [',num2str(Ndflt),'] : '];
 N = input(Ntxt);
end

if N~=[],
 if (N<Nmin),
  dispstr=['Warning : the number of analyzed voices (N) should be >= ',num2str(Nmin)];
  disp(dispstr);
 end
else
 N=Ndflt;
end

% Geometric sampling of the analyzed spectrum
No2 = N/2;
k = (1:No2);
q = (fmax/fmin)^(1/(No2-1));
t = (1:Mt)-M-1;
geo_f  = fmin*(exp((k-1).*log(q)));
tfmatx = zeros(Mt,N);
tfmatx = exp(-2*j*pi*t'*geo_f);
ZS = Z.'*tfmatx;
ZS(No2+1:N) = zeros(1,N-No2);

% Mellin transform computation of the analyzed signal
p = 0:(N-1);
L = log(fmin)/log(q);
mellin = N*log(q)*fftshift(ifft(ZS)).*exp(j*2*pi*L*(p/N-1/2));
beta   = (p/N-1/2)./log(q);

% Normalization
SP = fft(hilbert(real(X)));
indmin = 1+round(fmin*(xrow-2));
indmax = 1+round(fmax*(xrow-2));
SPana = SP(indmin:indmax);
nu = (indmin:indmax)'/N;
SPp = SPana./nu;
Normsig = sqrt(SPp'*SPana);

mellin = mellin*Normsig/norm(mellin);

author: HeLiangHIT


text:
function [cm2] = ft22cm2(ft2)
% Convert area from square feet to square centimeters.
% Chad A. Greene 2012
cm2 = ft2*929.0304;

author: Sable

repo: mcbench-benchmarks


text:
function linplus_test571 ( )

%*****************************************************************************80
%
%% TEST571 tests R8SP_INDICATOR.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2009
%
%  Author:
%
%    John Burkardt
%
  m = 7;
  n = 5;
  nz_num = 10;
  col = [ 2, 5, 1, 5, 1, 2, 3, 4, 4, 1 ];
  row = [ 1, 1, 2, 2, 4, 4, 4, 5, 6, 7 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST571\n' );
  fprintf ( 1, '  R8SP_INDICATOR sets up a R8SP indicator matrix;\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix rows M =    %d\n', m );
  fprintf ( 1, '  Matrix columns N = %d\n', n );
  fprintf ( 1, '  Matrix nonzeros =  %d\n', nz_num );

  a = r8sp_indicator ( m, n, nz_num, row, col );

  r8sp_print ( m, n, nz_num, row, col, a, '  The R8SP indicator matrix:' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [] = closematlabpool
   poolobj = gcp('nocreate');
   delete(poolobj);
end

author: ravenxrz

repo: Mathematical-Modeling


text:
function [ wval, dwdx ] = r8poly_lagrange_factor ( npol, xpol, xval )

%*****************************************************************************80
%
%% R8POLY_LAGRANGE_FACTOR evaluates the polynomial Lagrange factor at a point.
%
%  Formula:
%
%    W(X) = Product ( 1 <= I <= NPOL ) ( X - XPOL(I) )
%
%  Discussion:
%
%    Suppose F(X) is at least N times continuously differentiable in the
%    interval [A,B].  Pick NPOL distinct points XPOL(I) in [A,B] and compute
%    the interpolating polynomial P(X) of order NPOL ( and degree NPOL-1)
%    which passes through all the points ( XPOL(I), F(XPOL(I)) ).
%    Then in the interval [A,B], the maximum error
%
%      abs ( F(X) - P(X) )
%
%    is bounded by:
%
%      C * FNMAX * W(X)
%
%    where
%
%      C is a constant,
%      FNMAX is the maximum value of the NPOL-th derivative of F in [A,B],
%      W(X) is the Lagrange factor.
%
%    Thus, the value of W(X) is useful as part of an estimated bound
%    for the interpolation error.
%
%    Note that the Chebyshev abscissas have the property that they minimize
%    the value of W(X) over the interval [A,B].  Hence, if the abscissas may
%    be chosen arbitrarily, the Chebyshev abscissas have this advantage over
%    other choices.
%
%    For a set of points XPOL(I), 1 <= I <= NPOL, the IPOL-th Lagrange basis
%    polynomial L(IPOL)(X), has the property:
%
%      L(IPOL)( XPOL(J) ) = delta ( IPOL, J )
%
%    and may be expressed as:
%
%      L(IPOL)(X) = W(X) / ( ( X - XPOL(IPOL) ) * W'(XPOL(IPOL)) )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NPOL, the number of abscissas.
%    NPOL must be at least 1.
%
%    Input, real XPOL(NPOL), the abscissas, which should
%    be distinct.
%
%    Input, real XVAL, the point at which the Lagrange
%    factor is to be evaluated.
%
%    Output, real WVAL, the value of the Lagrange factor at XVAL.
%
%    Output, real DWDX, the derivative of W with respect to XVAL.
%
  wval = prod ( xval - xpol(1:npol) );

  dwdx = 0.0;

  for i = 1 : npol

    term = 1.0;

    for j = 1 : npol
      if ( i ~= j )
        term = term * ( xval - xpol(j) );
      end
    end

    dwdx = dwdx + term;

  end

  return
end


text:
function planC = getplanCDownSample(planC, optS, index)

indexS = planC{end};

for i = 1: size(planC{indexS.scan}.scanArray,3)
    planC{indexS.scan}.scanInfo(1,i).grid1Units = index*planC{indexS.scan}.scanInfo(1,i).grid1Units;
    planC{indexS.scan}.scanInfo(1,i).grid2Units = index*planC{indexS.scan}.scanInfo(1,i).grid2Units;
    planC{indexS.scan}.scanInfo(1,i).sizeOfDimension2 = size(planC{indexS.scan}.scanArray,2)/index;
    planC{indexS.scan}.scanInfo(1,i).sizeOfDimension1 = size(planC{indexS.scan}.scanArray,1)/index;
end

planC{indexS.scan}.scanArray = getDownsample3(planC{indexS.scan}.scanArray,index,1);

planC{indexS.scan}.scanArraySuperior = [];
planC{indexS.scan}.scanArrayInferior = [];
planC{indexS.scan}.uniformScanInfo = [];

for i=1:length(planC{indexS.structures})
    planC{indexS.structures}(i).rasterSegments =[];
end

planC = getRasterSegs(planC, optS);

planC = setUniformizedData(planC, optS);

try
	for i=1:length(planC{indexS.dose})
        planC = clearCache(planC, i);
	end
end


text:
function out = simbin(mat1,mat2,type,mask)

% Computes 106 measures of similarity and dissimilarity (distance) between
% two binary matrices. Matrices can be any dimensions, but must have the
% same dimensions. A mask can be used to indicate the relevant elements of
% the matrices (needed when the measure takes into account elements present
% in neither binary object). See further information below (including next
% to each measure).
%
% IMPORTANT: naming conventions for certain measures (e.g., the 5 Sokal and
% Sneath measures) are inconsistent across different researchers, different
% formulas for the same 'measure' exist in the literature, and different
% measures have similar names (e.g., Anderberg's vs. Anderberg's D). Thus,
% you should check that the formula used below corresponds to the
% particular measure that you wish to compute.
%
% Usage: out = overlap_similarity(mat1,mat2,type,mask)
%           mat1: 1st binary matrix
%           mat2: 2nd binary matrix
%           type: measure of (dis)similarity to be used
%           mask: binary mask indicating relevant matrix elements
%
%
% Further Information:
%
% All measures are symmetric (i.e., it does not matter which matrix is mat1
% vs. mat 2) except:
%                    AD, AMPLE, BUB1, BUB2, cole, diceasym1, diceasym2,
%                    digby, fleiss, forbes, GKlambda, GKtau, GK2, Ko2, LH,
%                    MP, Pcs, Pmsc, Pei1, Pei2, PHI, PH, Q0, RDeprob, SS4,
%                    SS5, tarantula, YQ, YQD, YY
%
% Computation of these measures takes into account (in some way) the
% elements present in neither binary object (i.e., d) (which may not be
% optimal if you want a straightforward measure of overlap):
%                    AD, AMPLE, baulieu, bc, Bshape, BUB1, BUB2, CK, cole,
%                    dennis, digby, disper, eyraud, faith, FD, fleiss,
%                    GKlambda, GKtau, GK1, GK2, GK3, GK4, GLS, goodall,
%                    gower, GW, hamann, HD, inprod, Kcam, Kpoai, LH,
%                    meanman, michael, MK, MP, pattern, Pcs, Pmsc, Pei1,
%                    Pei2, PHI, PH, Q0, RDeprob, Rpattern, RR, RT, scott,
%                    SM, Spattern, SS1, SS3, SS4, SS5, Stau, stiles,
%                    tarantula, tarwid, var, YQ, YQD, YY
%
% Positive monotomic relationship: bc, pattern
%                                  BSeuc, meanman, mirkin, var
%                                  BUB1, BUB2
%                                  DK, johnson
%                                  disper, MK, Stau
%                                  int, RR
%                                  hamann, inprod, SM
%                                  chord, hellinger
%                                  faith, tversky
%
% Negative monotomic relationship: BB with savage
%                                  BLWMN with dice
%                                  BSeuc, meanman, mirkin, and var with hamann, inprod, and SM
%                                  Bpattern with bc and pattern
%                                  Beuc with RT
%                                  steffensen with int and RR
%                                  YQ with YQD
%
% TAG:       MEASURE NAME:
% AMPLE      Analyzing Method Patterns to Locate Errors (Dallmeier et al, 2005) (similarity)
% AD         Anderberg's D (similarity)
% anderberg  Anderberg 1973 (similarity)
% baulieu    Baulieu 1989 (dissimilarity) (also called size difference)
% BB         Braun and Blanquet 1932 (similarity)
% bc         bc (dissimilarity)
% benini     Benini 1901 (similarity)
% Beuc       Binary euclidean distance (dissimilarity)
% BLWMN      Binary Lance and Williams nonmetric (dissimilarity) (also called Sokal and Sneath nonmetric or Bray and Curtis)
% Bpattern   Browsing pattern (similarity)
% BSeuc      Binary squared euclidean distance (dissimilarity) (also called Hamming, city block, Manhattan, or Filkov)
% Bshape     Binary shape (dissimilarity)
% BUB1       Baroni-Urbani and Buser 1976 1 (similarity)
% BUB2       Baroni-Urbani and Buser 1976 2 (similarity)
% chiY       Chi Square with correction of Yates
% chord      Chord (dissimilarity)
% CK         Cohen's kappa (similarity)
% cole       Cole 1949 (similarity)
% dennis     Dennis (similarity)
% dice       Dice (similarity) (also called Nei and Lei's genetic coefficient, Czekanowski, or Sorenson, or Bray)
% diceasym1  Dice asymmetric 1 (similarity)
% diceasym2  Dice asymmetric 2 (similarity)
% digby      Digby 1983 (similarity)
% disper     Dispersion (similarity)
% DK         Driver and Kroeber 1932 (similarity) (also called Kulczynski 1927 2)
% eyraud     Eyraud 1936 (similarity)
% fager      Fager (similarity)
% faith      Faith (similarity)
% FaMc       Fager and McGowan 1963 (similarity)
% FD         Forbes 1907's D (similarity)
% fleiss     Fleiss 1975 (similarity)
% forbes     Forbes 1925 (similarity)
% fossum     Fossum (similarity) (also called Jones and Curtis)
% gilbert    Gilbert's coefficient (similarity)
% gini       Gini 1912 index (similarity)
% GK1        Goodman and Kruskal 1954 1 (similarity)
% GK2        Goodman and Kruskal 1954 2 (similarity)
% GK3        Goodman and Kruskal 1954 3 (similarity)
% GK4        Goodman and Kruskal 1954 4 (similarity)
% GKlambda   Goodman and Kruskal's lambda (similarity)
% GKtau      Goodman and Kruskal's tau (similarity)
% goodall    Goodall 1967's angular transformation of simple matching coefficient (similarity) (also called Austin and Colwell)
% gower      Gower (similarity)
% GW         Gilbert and Wells 1966 (similarity)
% hamann     Hamann 1961 (similarity)
% HD         Hawkins and Dotson 1975 (similarity)
% hellinger  Hellinger (dissimilarity)
% inprod     Inner product (similarity)
% int        Intersection of two masks (similarity)
% jaccard    Jaccard 1908 (similarity) (also called similarity ratio or Tanimioto)
% johnson    Johnson 1967 (similarity) (also called McConnaughey 1964)
% Kcam       Kuhns 1965's coefficient of arithmetic means (similarity)
% Ko1        Koppen 1884 (similarity)
% Ko2        Koppen 1870 (similarity)
% Kpoai      Kuhns 1965's proportion of overlap above independence (similarity)
% kulczynski Kulczynski 1927 1 (similarity)
% meanman    Mean Manhattan (dissimilarity) (also called Canberra or Sneath total difference)
% michael    Michael 1920 (similarity)
% mirkin     Mirkin (dissimilarity)
% MK         Maron and Kuhns 1960 (similarity)
% modgini    Modified gini index (similarity)
% mountford  Mountford 1962 (similarity)
% MP         Maxwell and Pilliner 1968 (similarity)
% ochiai     Ochiai 1957 (similarity) (also called Driver and Kroeber or Otsuka)
% pattern    Pattern (dissimilarity)
% Pcs        Pearson 1905 chi-square (similarity)
% Pei1       Peirce 1884 1 (similarity)
% Pei2       Peirce 1884 2 (similarity)
% PH         Pearson and Heron 1913 (similarity)
% PHI        Fourfold point correlation (similarity) (also called Pearson and Heron 1913)
% Pmsc       Pearson 1905's coefficient of mean square contingency (similarity)
% Q0         Q0 (Batagelj and Bren 1995) (dissimilarity)
% Rcost      R cost
% RDeprob    Relative decrease of error probability (similarity)
% Rpattern   Retrieval pattern (similarity)
% RR         Russell and Rao 1940 (similarity)
% RT         Rogers and Tanimoto 1960 (similarity)
% savage     Savage (dissimilarity)
% Scost      S cost
% scott      Scott 1955 (similarity)
% simpson    Simpson 1943's ecological coexistence coefficient (similarity) (also called overlap)
% SM         Sokal and Michener 1958's simple matching (similarity) (also called Rand or Kendall)
% soergel    Soergel distance
% sorensen   Sorensen 1948 (similarity)
% sorgenfrei Sorgenfrei 1958 (similarity) (also called Fowlkes-Mallows)
% Spattern   Sneath pattern difference (dissimilarity)
% SS1        Sokal and Sneath 1963 1 (similarity)
% SS2        Sokal and Sneath 1963 2 (similarity)
% SS3        Sokal and Sneath 1963 3 (similarity)
% SS4        Sokal and Sneath 1963 4 (similarity) (also called Anderberg)
% SS5        Sokal and Sneath 1963 5 (similarity) (also called Ochiai 2)
% Stau       Stuart's tau (similarity)
% steffensen Steffensen 1934 (dissimilarity)
% stiles     Stiles (similarity)
% Tcost      T combined cost
% tarantula  Tarantula (Jones et al., 2002 & 2005) (similarity)
% tarwid     Tarwid 1960 (similarity)
% tversky    Tversky 1977's feature contrast model (similarity)
% Ucost      U cost
% unigram    Unigram subtuples
% US         Upholt 1977's S
% US         Upholt 1977's F
% var        Variance (dissimilarity)
% YQ         Yule 1911's Q coefficient of association (similarity)
% YQD        Yule's Q (Yule and Kendall, 1957) distance (dissimilarity)
% YY         Yule 1912's Y (or omega) coefficient of colligation (similarity)
%
%
% Author: Jeffrey M. Spielberg (jspielb2@gmail.com)
% Version: 05.06.16
%
% WARNING: This is a beta version. There no known bugs, but only limited
% testing has been perfomed. This software comes with no warranty (even the
% implied warranty of merchantability or fitness for a particular purpose).
% Therefore, USE AT YOUR OWN RISK!!!
%
% Copyleft 2014-2016. Software can be modified and redistributed, but
% modifed, redistributed versions must have the same rights

mat1 = double(mat1>0);
mat2 = double(mat2>0);

if nargin==2||isempty(type)
    type = 'dice';
    mask = ones(size(mat1));
elseif nargin<4
    mask = ones(size(mat1));
elseif nargin==4
    mask = double(mask>0);
end

mat1  = mat1.*mask;
mat2  = mat2.*mask;
nmat1 = not(mat1).*mask;
nmat2 = not(mat2).*mask;

a = sumn(mat1.*mat2);
b = sumn(nmat1.*mat2);
c = sumn(mat1.*nmat2);
d = sumn(nmat1.*nmat2);
tot = a+b+c+d;
q = (a+b)*(a+c)*(b+d)*(c+d);

switch type
    case 'AD'           % Anderberg's D (similarity) (measures the reduction in error probability when one item is used to predict the other) (range: 0:1)
        t1  = max(a,b)+max(c,d)+max(a,c)+max(b,d);
        t2  = max((a+c),(b+d))+max((a+d),(c+d));
        out = (t1-t2)/(2*tot);
    case 'AMPLE'        % Analyzing Method Patterns to Locate Errors (Dallmeier et al, 2005) (similarity) (undefined when a+b and/or c+d are 0)
        out = abs((a/(a+b))-(c/(c+d)));
    case 'anderberg'    % Anderberg 1973 (similarity)
        out = (8*a)/((8*a)+b+c);
    case 'baulieu'      % Baulieu 1989 (dissimilarity) (also called size difference) (range: 0:Inf)
        out = ((b-c)^2)/(tot^2);
    case 'BB'           % Braun and Blanquet 1932 (similarity) (range: 0:1)
        out = a/max((a+b),(a+c));
    case 'bc'           % bc (dissimilarity)
        out = (4*b*c)/(tot^2);
    case 'benini'       % Benini 1901 (similarity)
        out = (a-((a+b)*(a+c)))/(a+min(b,c)-((a+b)*(a+c)));
    case 'Beuc'         % Binary euclidean distance (dissimilarity) (range: 0:Inf)
        out = sqrt(b+c);
    case 'BLWMN'        % Binary Lance and Williams nonmetric (dissimilarity) (also called Sokal and Sneath nonmetric or Bray and Curtis) (range: 0:1)
        out = (b+c)/((2*a)+b+c);
    case 'Bpattern'     % Browsing pattern (similarity)
        out = a-(b*c);
    case 'BSeuc'        % Binary squared euclidean distance (dissimilarity) (also called Hamming, city block, Manhattan, or Filkov) (range: 0:Inf)
        out = b+c;
    case 'Bshape'       % Binary shape (dissimilarity) (range: no limits)
        out = ((tot*(b+c))-((b-c)^2))/(tot^2);
    case 'BUB1'         % Baroni-Urbani and Buser 1976 1 (similarity) (range: 0:1)
        out = (sqrt(a*d)+a)/(sqrt(a*d)+a+b+c);
    case 'BUB2'         % Baroni-Urbani and Buser 1976 2 (similarity)
        out = (sqrt(a*d)+a-b-c)/(sqrt(a*d)+a+b+c);
    case 'chiY'         % Chi Square with correction of Yates
        out = (tot*(abs((a*d)-(b*c))-(tot/2))^2)/q;
    case 'chord'        % Chord (dissimilarity)
        out = sqrt(2*(1-(a/sqrt((a+b)*(a+c)))));
    case 'CK'           % Cohen's kappa (similarity)
        out = (2*((a*d)-(b*c)))/(((a+b)*(b+d))+((a+c)*(c+d)));
    case 'cole'         % Cole 1949 (similarity) (undefined when a+b, b+d, a+c, and/or c+d are 0)
        out = ((a*d)-(b*c))/min(((a+b)*(a+c)),((b+d)*(c+d)));
    case 'dennis'       % Dennis (similarity)
        out = ((a*d)-(b*c))/sqrt(tot*(a+b)*(a+c));
    case 'dice'         % Dice (similarity) (also called Nei and Lei's genetic coefficient, Czekanowski, Sorenson, Bray, Upholt's F, Gower and Legendre's T, or percent positive agreement) (range: 0:1)
        out = (2*a)/((2*a)+b+c);
    case 'diceasym1'    % Dice asymmetric 1 (similarity)
        out = a/(a+c);
    case 'diceasym2'    % Dice asymmetric 2 (similarity)
        out = a/(a+b);
    case 'digby'        % Digby 1983 (similarity) (undefined when a+b, b+d, a+c, and/or c+d are 0)
        out = (((a*d)^(3/4))-((b*c)^(3/4)))/(((a*d)^(3/4))+((b*c)^(3/4)));
    case 'disper'       % Dispersion (similarity) (range: -1:1)
        out = ((a*d)-(b*c))/(tot^2);
    case 'DK'           % Driver and Kroeber 1932 (similarity) (also called Kulczynski 1927 2) (measures the conditional probability that both items are positive) (range: 0:1)
        out = (a/2)*((1/(a+b))+(1/(a+c)));
    case 'eyraud'       % Eyraud 1936 (similarity) (undefined when a+b, b+d, a+c, and/or c+d are 0)
        out = (a-((a+b)*(a+c)))/q;
        %out = ((tot^2)*((tot*a)-((a+b)*(a+c))))/q; % Alternative formula
    case 'fager'        % Fager (similarity)
        out = (a/sqrt((a+b)*(a+c)))-(1/(2*sqrt(min((a+b),(a+c)))));
        %out = (a/(((a+b)*(a+c))^2))-(max(b,c)/2); % Alternative formula
    case 'faith'        % Faith (similarity)
        out = (a+(0.5*d))/tot;
    case 'FaMc'         % Fager and McGowan 1963 (similarity)
        out = (a/sqrt((a+b)*(a+c)))-((a+max(b,c))/2);
    case 'FD'           % Forbes 1907's D (similarity) (also called Kocher and Wong)
        out = (a*tot)/((a+b)*(a+c));
    case 'fleiss'       % Fleiss 1975 (similarity) (undefined when a+b, b+d, a+c, and/or c+d are 0)
        out = (((a*d)-(b*c))*(((a+b)*(b+d))+((a+c)*(c+d))))/(2*q);
    case 'forbes'       % Forbes 1925 (similarity) (also called Loevinger's H) (undefined when a+b, b+d, a+c, and/or c+d are 0)
        out = ((tot*a)-((a+b)*(a+c)))/((tot*min((a+b),(a+c)))-((a+b)*(a+c)));
        %out = ((tot*a)-((a+b)*(a+c)))/((tot*min(b,c))-((a+b)*(a+c))); % Alternative formula
    case 'fossum'       % Fossum (similarity) (also called Jones and Curtis)
        out = (tot*((a-0.5)^2))/((a+b)*(a+c));
    case 'gilbert'      % Gilbert's coefficient (similarity) (undefined under certain circumstances)
        out = (a-((a+b)*(a+c)))/(a+b+c-((a+b)*(a+c)));
    case 'gini'         % Gini 1912 index (similarity) (undefined under certain circumstances)
        out = (a-((a+b)*(a+c)))/sqrt((1-((a+b)^2))*(1-((a+c)^2)));
    case 'GKlambda'     % Goodman and Kruskal's lambda (similarity) (measures the proportional reduction in error using one item to predict another) (range: 0:1)
        t1  = max(a,b)+max(c,d)+max(a,c)+max(b,d);
        t2  = max((a+c),(b+d))+max((a+d),(c+d));
        out = (t1-t2)/((2*tot)-t2);
    case 'GKtau'        % Goodman and Kruskal's tau (similarity) (undefined when a+b and/or c+d are 0)
        out = (((((a-((a+b)*(a+c)))^2)+((b-((a+b)*(b+d)))^2))/(a+b))/(1-((a+c)^2)-((b+d)^2)))+(((((c-((a+c)*(c+d)))^2)+((d-((b+d)*(c+d)))^2))/(c+d))/(1-((a+c)^2)-((b+d)^2)));
    case 'GK1'          % Goodman and Kruskal 1 1954 (similarity)
        out = ((2*min(a,d))-b-c)/((2*min(a,d))+b+c);
    case 'GK2'          % Goodman and Kruskal 2 1954 (similarity)
        out = (max(a,c)+max(b,d)-max((a+b),(c+d)))/(1-max((a+b),(c+d)));
    case 'GK3'          % Goodman and Kruskal 3 1954 (similarity)
        out = (a+d-max(a,d)-((b+c)/2))/(1-max(a,d)-((b+c)/2));
    case 'GK4'          % Goodman and Kruskal 4 1954 (similarity)
        out = ((max(a,b)+max(c,d)+max(a,c)+max(b,d))/(2-max((a+c),(b+d))-max((a+b),(c+d))))-((max((a+c),(b+d))+max((a+b),(c+d)))/(2-max((a+c),(b+d))-max((a+b),(c+d))));
    case 'goodall'      % Goodall 1967's angular transformation of simple matching coefficient (similarity) (also called Austin and Colwell)
        out = asin(sqrt((a+d)/tot))/(50*pi);
    case 'gower'        % Gower (similarity) (undefined when a+b, b+d, a+c, and/or c+d are 0)
        out = (a+d)/sqrt(q);
    case 'GW'           % Gilbert and Wells 1966 (similarity) (undefined when a+b and/or a+c are 0)
        out = log((a*tot)/((a+b)*(a+c)));
    case 'hamann'       % Hamann 1961 (similarity) (range: -1:1)
        out = ((a+d)-(b+c))/tot;
    case 'HD'           % Hawkins and Dotson 1975 (similarity)
        out = 0.5*((a/(a+b+c))+(d/(b+c+d)));
    case 'hellinger'    % Hellinger (dissimilarity)
        out = 2*sqrt(1-(a/sqrt((a+b)*(a+c))));
    case 'inprod'       % Inner product (similarity)
        out = a+d;
    case 'int'          % Intersection of two masks (similarity)
        out = a;
    case 'jaccard'      % Jaccard 1908 (similarity) (also called similarity ratio or Tanimioto) (range: 0:1)
        out = a/(a+b+c);
    case 'johnson'      % Johnson 1967 (similarity) (also called McConnaughey 1964)
        out = (a/(a+b))+(a/(a+c));
    case 'Kcam'         % Kuhns 1965's coefficient of arithmetic means (similarity)
        out = (2*((a*d)-(b*c)))/(tot*((2*a)+b+c));
    case 'Ko1'          % Koppen 1884 (similarity)
        out = a+((b+c)/2);
    case 'Ko2'          % Koppen 1870 (similarity) (undefined when a+b and/or 1-a-b are 0)
        out = (((a+b)*(1-a-b))-c)/((a+b)*(1-a-b));
    case 'Kpoai'        % Kuhns 1965's proportion of overlap above independence (similarity) (undefined when a+c and/or a+c are 0)
        out = ((a*d)-(b*c))/(tot*(a/((a+b)*(a+c)))*((2*a)+b+c-(((a+b)*(a+c))/tot)));
    case 'kulczynski'   % Kulczynski 1927 1 (similarity) (range: 0:Inf)
        out = a/(b+c);
    case 'meanman'      % Mean Manhattan (dissimilarity) (also called Canberra or Sneath total difference)
        out = (b+c)/tot;
    case 'michael'      % Michael 1920 (similarity) (range: 0:1)
        out = (4*((a*d)-(b*c)))/(((a+d)^2)+((b+c)^2));
    case 'mirkin'       % Mirkin (dissimilarity)
        out = 2*(b+c);
    case 'MK'           % Maron and Kuhns 1960 (similarity)
        out = ((a*d)-(b*c))/tot;
    case 'modgini'      % Modified gini index (similarity)
        out = (a-((a+b)*(a+c)))/(1-(abs(b-c)/2)-((a+b)*(a+c)));
    case 'mountford'    % Mountford 1962 (similarity)
        out = a/((0.5*((a*b)+(a*c)))+(b*c));
    case 'MP'           % Maxwell and Pilliner 1968 (similarity) (undefined when a+b, b+d, a+c, and/or c+d are 0)
        out = (2*((a*d)-(b*c)))/q;
    case 'ochiai'       % Ochiai 1957 (similarity) (also called Driver and Kroeber or Otsuka) (binary form of cosine similarity) (range: 0:1)
        out = a/sqrt((a+b)*(a+c));
    case 'pattern'      % Pattern (dissimilarity) (range: 0:1)
        out = (b*c)/(tot^2);
    case 'Pcs'          % Pearson 1905 chi-square (similarity) (undefined when a+b, b+d, a+c, and/or c+d are 0)
        out = (tot*(((a*d)-(b*c))^2))/q;
    case 'Pmsc'         % Pearson 1905's coefficient of mean square contingency (similarity) (undefined when a+b, b+d, a+c, and/or c+d are 0)
        ch2 = (tot*(((a*d)-(b*c))^2))/q;
        out = sqrt(ch2/(tot+ch2));
    case 'Pei1'         % Peirce 1884 1 (similarity) (undefined when a+c and/or b+d are 0)
        out = ((a*b)-(b*c))/((a+c)*(b+d));
    case 'Pei2'         % Peirce 1884 2 (similarity) (undefined when a+c and/or b+d are 0)
        out = ((a*b)+(b*c))/((a*b)+(2*b*c)+(c*d));
    case 'PHI'          % Fourfold point correlation (similarity) (also called Pearson and Heron 1913) (binary form of the Pearson product-moment correlation) range: 0:1) (undefined when a+b, b+d, a+c, and/or c+d are 0)
        out = ((a*d)-(b*c))/sqrt(q);
    case 'PH'           % Pearson and Heron 1913 (similarity) (undefined when a+b, b+d, a+c, and/or c+d are 0)
        out = cos((180*sqrt(b*c))/(sqrt(a*d)+sqrt(b*c)));
    case 'Q0'           % Q0 (Batagelj and Bren 1995) (dissimilarity) (range: 0:Inf) (undefined when a+b, b+d, a+c, and/or c+d are 0)
        out = (b*c)/(a*d);
    case 'Rcost'        % R cost
        out = log(1+(a/(a+b)))*log(1+(a/(a+c)));
    case 'RDeprob'      % Relative decrease of error probability (similarity)
        out = (max(a,b)+max(c,d)-max((a+c),(b+d)))/(1-max((a+c),(b+d)));
    case 'Rpattern'     % Retrieval pattern (similarity)
        out = a*d;
    case 'RR'           % Russell and Rao 1940 (similarity) (range: 0:1)
        out = a/tot;
    case 'RT'           % Rogers and Tanimoto 1960 (similarity) (binary dot product) (range: 0:1)
        out = (a+d)/(a+d+(2*(b+c)));
    case 'savage'       % Savage (dissimilarity)
        out = 1-(a/(a+max(b,c)));
    case 'Scost'        % S cost
        out = log(1+min(b,c)/(a+1))^-.5;
    case 'scott'        % Scott 1955 (similarity)
        out = ((4*((a*d)-(b*c)))-((b-c)^2))/(((2*a)+b+c)+(b+c+(2*d)));
    case 'simpson'      % Simpson 1943's ecological coexistence coefficient (similarity) (also called overlap) (range: 0:1)
        out = a/min((a+b),(a+c));
    case 'SM'           % Sokal and Michener 1958's simple matching (similarity) (also called Rand or Kendall) (range: 0:1)
        out = (a+d)/tot;
    case 'soergel'      % Soergel distance
        out = (b+c)/(b+c+d);
    case 'sorensen'     % Sorensen 1948 (similarity)
        out = (4*a)/((4*a)+b+c);
    case 'sorgenfrei'   % Sorgenfrei 1958 (similarity) (also called Fowlkes-Mallows)
        out = (a^2)/((a+b)*(a+c));
    case 'Spattern'     % Sneath pattern difference (dissimilarity)
        out = (2*sqrt(b*c))/tot;
    case 'SS1'          % Sokal and Sneath 1963 1 (similarity) (also called Gower and Legendre's S)
        out = (2*(a+d))/((2*(a+d))+b+c);
    case 'SS2'          % Sokal and Sneath 1963 2 (similarity)
        out = a/(a+(2*(b+c)));
    case 'SS3'          % Sokal and Sneath 1963 3 (similarity) (range: 0:Inf)
        out = (a+d)/(b+c);
    case 'SS4'          % Sokal and Sneath 1963 4 (similarity) (also called Anderberg) (measures the conditional probability that both items are in the same state [present or absent]) (range: 0:1) (undefined when a+b, b+d, a+c, and/or c+d are 0)
        out = ((a/(a+b))+(a/(a+c))+(d/(b+d))+(d/(c+d)))/4;
    case 'SS5'          % Sokal and Sneath 1963 5 (similarity) (also called Ochiai 2) (range: 0:1) (undefined when a+b, b+d, a+c, and/or c+d are 0)
        out = (a*d)/sqrt(q);
    case 'Stau'         % Stuart's tau (similarity)
        out = 2*((a*d)-(b*c));
    case 'steffensen'   % Steffensen 1934 (dissimilarity)
        out = (a-((a+b)*(a+c)))/2;
    case 'stiles'       % Stiles (similarity) (undefined when a+b, b+d, a+c, and/or c+d are 0 and other circumstances)
        out = log10((tot*((abs((a*d)-(b*c))-(tot/2))^2))/q);
    case 'tarantula'    % Tarantula (Jones et al., 2002 & 2005) (similarity) (undefined when a+b and/or c+d are 0)
        out = (a/(a+b))/((a/(a+b))+(c/(c+d)));
    case 'tarwid'       % Tarwid 1960 (similarity)
        out = ((tot*a)-((a+b)*(a+c)))/((tot*a)+((a+b)*(a+c)));
    case 'Tcost'        % T combined cost
        R = log(1+(a/(a+b)))*log(1+(a/(a+c)));
        S = log(1+min(b,c)/(a+1))^-.5;
        U = log(1+((min(b,c)+a)/(max(b,c)+a)));
        out = sqrt(U*S*R);
    case 'tversky'      % Tversky 1977's feature contrast model (similarity)
        out = a-b-c;
    case 'Ucost'        % U cost
        out = log(1+((min(b,c)+a)/(max(b,c)+a)));
    case 'unigram'      % Unigram subtuples
        out = log((a*d)/(b*c))-(3.29*sqrt((1/a)+(1/b)+(1/c)+(1/d)));
    case 'UF'           % Upholt 1977's F (also called Gower and Legendre's T)
        out = (2*a)/((2*a)+b+c);
    case 'US'           % Upholt's S 1977
        F = (2*a)/((2*a)+b+c);
        out = (0.5*(-F+sqrt((F^2)+(8*F))))^(1/tot);
    case 'var'          % Variance (dissimilarity) (range: 0:Inf)
        out = (b+c)/(4*tot);
    case 'YQ'           % Yule 1911's Q coefficient of association (similarity) (undefined when a+b, b+d, a+c, and/or c+d are 0)
        out = ((a*d)-(b*c))/((a*d)+(b*c));
    case 'YQD'          % Yule's Q distance (dissimilarity) (range: -1:1) (undefined when a+b, b+d, a+c, and/or c+d are 0)
        out = (b*c)/((a*d)+(b*c));
    case 'YY'           % Yule 1912's Y (or omega) coefficient of colligation (similarity) (range: -1:1) (undefined when a+b, b+d, a+c, and/or c+d are 0)
        out = (sqrt(a*d)-sqrt(b*c))/(sqrt(a*d)+sqrt(b*c));
end

function out = sumn(in)
out = sum(in(:));

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function phase=resolveLogScaleEccToPhase(Ecc,phaseZero,EccZero,noRings)

%   phase=resolveLogScaleEccToPhase(Ecc,phaseZero,EccZero,noRings)
%   PURPOSE: find the phase an eccentricity reflects when the ringmapping was logscaled.
%   (inverse of resolveLogScalePhaseToEcc)
%
%--------------------------------------------------------------------------
%   USAGE:  phase       = the phase that needs to be translated,
%           phaseZero   = the lowest phase (or reference phase)
%           EccZero     = the Ecc the phaseZero reflects
%           noRings     = number of stimuli in the ringmapping, usually 8,
%                         importent because noRing makes 2pi
%--------------------------------------------------------------------------
%   HISTORY
%
%   2005.03.14 by Mark Schira mark@ski.org

tau=noRings/pi/2;
phase=(log2(Ecc/EccZero)/tau)+phaseZero;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function test_bug1397

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_preprocessing ft_appenddata

% the following code was obtained from http://www.fieldtriptoolbox.org/tutorial/coherence
% on Wed Mar 28 15:36:40 CEST 2012

% find the interesting epochs of data
cfg = [];
% MODIFICATION, use trialfun handle and other path to the data
cfg.trialfun                  = @trialfun_left;
cfg.dataset                   = dccnpath('/home/common/matlab/fieldtrip/data/ftp/test/ctf/SubjectCMC.ds');
cfg = ft_definetrial(cfg);

% MODIFICATION, use only 10 trials
cfg.trl = cfg.trl(1:10,:);

% MODIFICATION, the following should not affect the problem
%
% % detect EOG artifacts in the MEG data
% cfg.continuous                = 'yes';
% cfg.artfctdef.eog.padding     = 0;
% cfg.artfctdef.eog.bpfilter    = 'no';
% cfg.artfctdef.eog.detrend     = 'yes';
% cfg.artfctdef.eog.hilbert     = 'no';
% cfg.artfctdef.eog.rectify     = 'yes';
% cfg.artfctdef.eog.cutoff      = 2.5;
% cfg.artfctdef.eog.interactive = 'no';
% cfg = ft_artifact_eog(cfg);
%
% % detect jump artifacts in the MEG data
% cfg.artfctdef.jump.interactive = 'no';
% cfg.padding                    = 5;
% cfg = ft_artifact_jump(cfg);
%
% % detect muscle artifacts in the MEG data
% cfg.artfctdef.muscle.cutoff      = 8;
% cfg.artfctdef.muscle.interactive = 'no';
% cfg = ft_artifact_muscle(cfg);
%
% % reject the epochs that contain artifacts
% cfg.artfctdef.reject          = 'complete';
% cfg = ft_rejectartifact(cfg);

% preprocess the MEG data
cfg.demean                    = 'yes';
cfg.dftfilter                 = 'yes';
cfg.channel                   = {'MEG'};
cfg.continuous                = 'yes';
meg = ft_preprocessing(cfg);

cfg              = [];
cfg.dataset      = meg.cfg.dataset;
cfg.trl          = meg.cfg.trl;
cfg.continuous   = 'yes';
cfg.demean       = 'yes';
cfg.dftfilter    = 'yes';
cfg.channel      = {'EMGlft' 'EMGrgt'};
cfg.hpfilter     = 'yes';
cfg.hpfreq       = 10;
cfg.rectify      = 'yes';
emg = ft_preprocessing(cfg);

% see http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=1397
% the reported problem in fieldtrip-20120302 was
%
% ??? Error using ==> ft_appenddata at 266
% there is a difference in the time axes of the input data

data = ft_appenddata([], meg, emg);
assert(numel(data.label)==153);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION to avoid external dependencies of this test
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function trl = trialfun_left(cfg)

% read in the triggers and create a trial-matrix
% consisting of 1-second data segments, in which
% left ECR-muscle is active.

event = ft_read_event(cfg.dataset);
trig  = [event(find(strcmp('backpanel trigger', {event.type}))).value];
indx  = [event(find(strcmp('backpanel trigger', {event.type}))).sample];

%left-condition
sel = [find(trig==1028):find(trig==1029)];

trig = trig(sel);
indx = indx(sel);

trl = [];
for j = 1:length(trig)-1
  trg1 = trig(j);
  trg2 = trig(j+1);
  if trg1<=100 && trg2==2080
    trlok      = [[indx(j)+1:1200:indx(j+1)-1200]' [indx(j)+1200:1200:indx(j+1)]'];
    trlok(:,3) = [0:-1200:-1200*(size(trlok,1)-1)]';
    trl        = [trl; trlok];
  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% Example script how to use PIVlab from the commandline
% Just run this script to see what it does.
% You can adjust the settings in "s" and "p", specify a mask and a region of interest
clc; clear all

nr_of_cores = 1; % integer, 1 means single core, greater than 1 means parallel
if nr_of_cores > 1
	try
		local_cluster=parcluster('local'); % single node
		corenum =  local_cluster.NumWorkers ; % fix : get the number of cores available
	catch
		warning('on');
		warning('parallel local cluster can not be created, assigning number of cores to 1');
		nr_of_cores = 1;
	end
end
%% Create list of images inside user specified directory
directory= fullfile(fileparts(mfilename('fullpath')) , 'Examples') ; %directory containing the images you want to analyze
% default directory: PIVlab/Examples
addpath(fileparts(mfilename('fullpath')));
disp(fileparts(mfilename('fullpath')))

suffix='*.jpg'; %*.bmp or *.tif or *.jpg or *.tiff or *.jpeg
disp(['Looking for ' suffix ' files in the selected directory.']);
direc = dir([directory,filesep,suffix]); filenames={};
[filenames{1:length(direc),1}] = deal(direc.name);
filenames = sortrows(filenames); %sort all image files
amount = length(filenames);

%% Standard PIV Settings
s = cell(15,2); % To make it more readable, let's create a "settings table"
%Parameter                          %Setting           %Options
s{1,1}= 'Int. area 1';              s{1,2}=64;         % window size of first pass
s{2,1}= 'Step size 1';              s{2,2}=32;         % step of first pass
s{3,1}= 'Subpix. finder';           s{3,2}=1;          % 1 = 3point Gauss, 2 = 2D Gauss
s{4,1}= 'Mask';                     s{4,2}=[];         % If needed, generate via: imagesc(image); [temp,Mask{1,1},Mask{1,2}]=roipoly;
s{5,1}= 'ROI';                      s{5,2}=[];         % Region of interest: [x,y,width,height] in pixels, may be left empty
s{6,1}= 'Nr. of passes';            s{6,2}=2;          % 1-4 nr. of passes
s{7,1}= 'Int. area 2';              s{7,2}=32;         % second pass window size
s{8,1}= 'Int. area 3';              s{8,2}=16;         % third pass window size
s{9,1}= 'Int. area 4';              s{9,2}=16;         % fourth pass window size
s{10,1}='Window deformation';       s{10,2}='*linear'; % '*spline' is more accurate, but slower
s{11,1}='Repeated Correlation';     s{11,2}=0;         % 0 or 1 : Repeat the correlation four times and multiply the correlation matrices.
s{12,1}='Disable Autocorrelation';  s{12,2}=0;         % 0 or 1 : Disable Autocorrelation in the first pass.
s{13,1}='Correlation style';        s{13,2}=0;         % 0 or 1 : Use circular correlation (0) or linear correlation (1).
s{14,1}='Repeat last pass';   s{14,2}=0; % 0 or 1 : Repeat the last pass of a multipass analyis
s{15,1}='Last pass quality slope';   s{15,2}=0.025; % Repetitions of last pass will stop when the average difference to the previous pass is less than this number.

%% Standard image preprocessing settings
p = cell(10,1);
%Parameter                       %Setting           %Options
p{1,1}= 'ROI';                   p{1,2}=s{5,2};     % same as in PIV settings
p{2,1}= 'CLAHE';                 p{2,2}=1;          % 1 = enable CLAHE (contrast enhancement), 0 = disable
p{3,1}= 'CLAHE size';            p{3,2}=50;         % CLAHE window size
p{4,1}= 'Highpass';              p{4,2}=0;          % 1 = enable highpass, 0 = disable
p{5,1}= 'Highpass size';         p{5,2}=15;         % highpass size
p{6,1}= 'Clipping';              p{6,2}=0;          % 1 = enable clipping, 0 = disable
p{7,1}= 'Wiener';                p{7,2}=0;          % 1 = enable Wiener2 adaptive denoise filter, 0 = disable
p{8,1}= 'Wiener size';           p{8,2}=3;          % Wiener2 window size
p{9,1}= 'Minimum intensity';     p{9,2}=0.0;        % Minimum intensity of input image (0 = no change)
p{10,1}='Maximum intensity';     p{10,2}=1.0;       % Maximum intensity on input image (1 = no change)

%% other settings
pairwise = 1; % 0 for [A+B], [B+C], [C+D]... sequencing style, and 1 for [A+B], [C+D], [E+F]... sequencing style

%% PIV analysis loop
if pairwise == 1
	if mod(amount,2) == 1 %Uneven number of images?
		disp('Image folder should contain an even number of images.')
		%remove last image from list
		amount=amount-1;
		filenames(size(filenames,1))=[];
	end

	disp(['Found ' num2str(amount) ' images (' num2str(amount/2) ' image pairs).'])
	x=cell(amount/2,1);
	y=x;
	u=x;
	v=x;
else
	disp(['Found ' num2str(amount) ' images'])
	x=cell(amount-1,1);
	y=x;
	u=x;
	v=x;
end

typevector=x; %typevector will be 1 for regular vectors, 0 for masked areas
correlation_map=x; % correlation coefficient

%% Pre-load the image names out side of the parallelizable loop
slicedfilename1=cell(0);
slicedfilename2=cell(0);
j = 1;
for i=1:1+pairwise:amount-1
	slicedfilename1{j}=filenames{i}; % begin
	slicedfilename2{j}=filenames{i+1}; % end
	j = j+1;
end

%% Main PIV analysis loop:
% parallel
if nr_of_cores > 1

	if pivparpool('size')<nr_of_cores
		pivparpool('open',nr_of_cores);
	end

	parfor i=1:size(slicedfilename1,2)  % index must increment by 1

		[x{i}, y{i}, u{i}, v{i}, typevector{i},correlation_map{i}] = ...
			piv_analysis(directory, slicedfilename1{i}, slicedfilename2{i},p,s,nr_of_cores,false);
	end
else % sequential loop

	for i=1:size(slicedfilename1,2)  % index must increment by 1

		[x{i}, y{i}, u{i}, v{i}, typevector{i},correlation_map{i}] = ...
			piv_analysis(directory, slicedfilename1{i}, slicedfilename2{i},p,s,nr_of_cores,true);

		disp([int2str((i+1)/amount*100) ' %']);

	end
end

%% PIV postprocessing loop
% Standard image post processing settings

r = cell(6,1);
%Parameter     %Setting                                     %Options
r{1,1}= 'Calibration factor, 1 for uncalibrated data';      r{1,2}=1;                   % Calibration factor for u
r{2,1}= 'Calibration factor, 1 for uncalibrated data';      r{2,2}=1;                   % Calibration factor for v
r{3,1}= 'Valid velocities [u_min; u_max; v_min; v_max]';    r{3,2}=[-50; 50; -50; 50];  % Maximum allowed velocities, for uncalibrated data: maximum displacement in pixels
r{4,1}= 'Stdev check?';                                     r{4,2}=1;                   % 1 = enable global standard deviation test
r{5,1}= 'Stdev threshold';                                  r{5,2}=7;                   % Threshold for the stdev test
r{6,1}= 'Local median check?';                              r{6,2}=1;                   % 1 = enable local median test
r{7,1}= 'Local median threshold';                           r{7,2}=3;                   % Threshold for the local median test

u_filt=cell(size(u));
v_filt=cell(size(v));
typevector_filt=typevector;

if nr_of_cores >1 % parallel

	if pivparpool('size')<nr_of_cores
		pivparpool('open',nr_of_cores);
	end

	parfor PIVresult=1:size(x,1)

		[u_filt{PIVresult,1}, v_filt{PIVresult,1},typevector_filt{PIVresult,1}]= ...
			post_proc_wrapper(u{PIVresult,1},v{PIVresult,1},typevector{PIVresult,1},r,true);

	end

else % sequential loop

	for PIVresult=1:size(x,1)

		[u_filt{PIVresult,1}, v_filt{PIVresult,1},typevector_filt{PIVresult,1}]= ...
			post_proc_wrapper(u{PIVresult,1},v{PIVresult,1},typevector{PIVresult,1},r,true);

	end

end

%% clean up parallel pool, and cluster
%{
if nr_of_cores >1 % parallel
	poolobj = gcp('nocreate'); % GET the current parallel pool
	if ~isempty(poolobj ); delete(poolobj );end
	clear local_cluster;
end
%}
%%
save(fullfile(directory, [filenames{1} '_' filenames{end} '_' num2str(amount) '_frames_result_.mat']));

%%
clearvars -except p s r x y u v typevector directory filenames u_filt v_filt typevector_filt correlation_map
disp('DONE.')

function [u_filt, v_filt,typevector_filt] = post_proc_wrapper(u,v,typevector,post_proc_setting,paint_nan)
% wrapper function for PIVlab_postproc

% INPUT
% u, v: u and v components of vector fields
% typevector: type vector
% post_proc_setting: post processing setting
% paint_nan: bool, whether to interpolate missing data

% OUTPUT
% u_filt, v_filt: post-processed u and v components of vector fields
% typevector_filt: post-processed type vector

[u_filt,v_filt] = PIVlab_postproc(u,v, ...
	post_proc_setting{1,2},...
	post_proc_setting{2,2},...
	post_proc_setting{3,2},...
	post_proc_setting{4,2},...
	post_proc_setting{5,2},...
	post_proc_setting{6,2},...
	post_proc_setting{7,2});

typevector_filt = typevector; % initiate
typevector_filt(isnan(u_filt))=2;
typevector_filt(isnan(v_filt))=2;
typevector_filt(typevector==0)=0; %restores typevector for mask

if paint_nan
	u_filt=inpaint_nans(u_filt,4);
	v_filt=inpaint_nans(v_filt,4);
end

end

author: Shrediquette

repo: PIVlab

save_path: github-repos/MATLAB/Shrediquette-PIVlab


text:
function sessions = selectSessions(studyDir,nSessions);
%
% sessions = selectSessions(studyDir,[nSessions]);
%
% Select multiple sessions. studyDir is the path
% to the "study", i.e., the path where most sessions
% are kept. Defaults to the parent of pwd. Also provides
% the option to manually select different sessions, if
% they're located all over the place, using uigetmanydirs.
%
% Returns a cell containing full path to each session.
%
% I wrote this on top of uigetmanydirs b/c I found that
% interface to be pretty slow and annoying (the default
% uigetdir for windows is always low and cut-off on my
% machine), plus I usually keep my session folders in one
% place anyway.
%
% nSessions: optional argument specifying # of
% sessions to select. If omitted or 0, takes
% all selected. Otherwise, takes first nSessions.
%
% ras 03/05
if notDefined('studyDir')
    studyDir = fileparts(pwd);
end

if notDefined('nSessions')
    nSessions = 0;
end

sessions = [];

% get names of directories
w = dir(studyDir);
names = {w.name};
dirCheck = [w.isdir];
names = names(dirCheck==1); % only dirs
names = setdiff(names,{'.' '..'});

% first option will be to manually locate sesions
% outside the guessed study dir
names = [{'(Find Other Session)'} names];

[whichSessions, ok] = listdlg('PromptString','Use coords from which sessions?',...
                        'ListSize',[400 600],...
                        'ListString',names,...
                        'InitialValue',1,...
                        'OKString','OK');
if ~ok  return;  end

% if subject requested to find other session, allow to
% add that to list:
if ismember(whichSessions,1)
    if nSessions==1
        while 1
            sess = uigetdir(pwd, 'Select a session to import');
            if exist(fullfile(sess, 'Gray', 'coords.mat'), 'file')
                sessions = {sess};
                break;
            else
                h = warndlg(['Selected session doesn''t have a ' ...
                           'Gray/coords.mat file. Please choose again.']);
                set(h, 'CloseRequestFcn', 'uiresume; closereq; ');
                uiwait
            end
        end
    else
        sessions = uigetmanydirs('Gray/coords.mat');
    end
end

% also get any selected from the listdlg, and add the
% parent dir to get a full path:
sessInds = whichSessions(whichSessions>1);
tmpSessions = names(sessInds);
for i = 1:length(tmpSessions)
    tmpSessions{i} = fullfile(studyDir,tmpSessions{i});
end
sessions = [sessions tmpSessions];

if nSessions > 0
    nSessions = min(length(sessions),nSessions);
    sessions = sessions(1:nSessions);
end

return


text:
function dt6 = dtiMakeDt6FromFsl(b0FileName, t1FileName)
%
% dt6 = dtiMakeDt6FromFsl([b0FileName], [t1FileName])
%
% Takes a set of images from FSL's DTI image files and saves them out as
% in a dt6 matlab file. You should specify the FSL S0 (usually b=0) file as
% input (a file dialog will pop-up if you don't).
%
% The t1 file is optional- if it is omitted, then the b=0 will be used as
% the anatomical reference.
%
% HISTORY:
% 2005.08.04 RFD (bob@white.stanford.edu) wrote it (based on dtiMakeDt6)
% 2006.04.27 RFD: adapted it for the Phillips data at MUSC. The old version
% (customized for data from U Washington) is now saved as
% dtiMakeDt6FromFsl_UW.

slice = [0,0,0]; % slice to display for checking coregistration
bb = [-80,80; -120,90; -60,90]';
mmDt = [1 1 1];
mmAnat = [1 1 1];
autoAlign = true;
showFigs = true;
tensorInterpSpline = false;

if ~exist('b0FileName','var') | isempty(b0FileName)
    [f, p] = uigetfile({'*.nii.gz','NIFTI gz';'*.*','All files'}, 'Select the FSL S0 file...');
    if(isnumeric(f)) error('Need an S0 file to continue...'); end
    b0FileName = fullfile(p, f);
    disp(b0FileName);
end
b0 = niftiRead(b0FileName);
% Fix ill-specified quaternion xforms so that they do something reasonable.
if(all(b0.qto_ijk(1:3,4) == [0 0 0]') | any(b0.qto_ijk(1:3,4)>1000))
    % *** HACK! [0.95 1.1 0.7] is near the ac for one data set
    if(det(b0.qto_ijk(1:3,1:3))<0)
        flipDti = 1;
        b0.data = flipdim(b0.data,1);
    end
    sz = size(b0.data);
    origin = round(sz(1:3)./2.*[0.95 1.1 0.7]);
    b0.qto_ijk = affineBuild(origin'+0.5, [0 0 0], 1./b0.pixdim, [0 0 0]);
    b0.qto_xyz = inv(b0.qto_ijk);
end
b0.data = double(b0.data);
% Clip vals <0 (these should be junk for a b=0 image)
b0.data(b0.data<0) = 0;

if ~exist('t1FileName','var')
    [f, p] = uigetfile({'*.nii.gz','NIFTI gz';'*.hdr','Analyze files (*.hdr)';'*.*','All files'},...
        'Select the T1 file...', b0FileName);
    if(isnumeric(f))
        t1FileName = [];
    else
        t1FileName = fullfile(p, f);
        disp(t1FileName);
    end
end
if(~isempty(t1FileName))
    t1 = niftiRead(t1FileName);
    % Fix ill-specified quaternion xforms so that they do something reasonable.
    if(all(t1.qto_ijk(1:3,4) == [0 0 0]'))
        sz = size(t1.data);
        origin = round(sz(1:3)./2.*[0.95 1.1 0.9]);
        t1.qto_ijk(1:3,4) = origin'+0.5;
        t1.qto_xyz = inv(t1.qto_ijk);
    end
    t1.data = double(t1.data);
    t1.data(t1.data<0) = 0;
else
    t1 = b0;
end

[datapath,basename] = fileparts(b0FileName);
us = strfind(basename,'_');
basename = basename(1:us(end));

if(~exist('outPathName','var') | isempty(outPathName))
    outPathName = fullfile(datapath, [basename 'dt6.mat']);
end
if(exist(outPathName,'file') | exist([outPathName '.mat'],'file'))
  disp('output file exists- please rename...');
  [f,p] = uiputfile('*.mat', 'Select output file...', outPathName);
  if(isnumeric(f)) error('User cancelled.'); end
  outPathName = fullfile(p,f);
end

disp(['Data will be saved in ' outPathName '...']);

% *** TO DO: fix this- here we assume that the data use a NIFTI-style
% header and this it's quaternion xform is specified correctly.
t1.acpcXform = t1.qto_xyz;
b0.acpcXform = b0.qto_xyz;

h = figure;
showFigure(h, t1, b0, bb, slice);
%set(h,'Position', [10, 50, 1000, 1000]);
set(h, 'PaperPositionMode', 'auto');
print(h, '-dpng', '-r90', [outPathName,'_headerAlign.png']);
if(~showFigs) close(h); end

if(autoAlign)
    spm_defaults;
    if(~exist('defaults','var')) global defaults; end
    disp('Coregistering (using spm tools)...');
    img = t1.data;
    img = mrAnatHistogramClip(img, 0.50, 0.995);
    VG.uint8 = uint8(round(img*255));
    VG.mat = t1.acpcXform;

    img = b0.data;
    img = mrAnatHistogramClip(img,0.40, 0.99);
    VF.uint8 = uint8(round(img*255));
    VF.mat = b0.acpcXform;
    p = defaults.coreg.estimate;
    % NOTE: there seems to be a consistent 1mm (1/2 of B0 voxel) translation
    % between the t1 and B0.
    transRot = spm_coreg(VG,VF,p);
    transRot(1:3) = transRot(1:3)+mmDt/2;
    b0.acpcXform = inv(VF.mat\spm_matrix(transRot(:)'));
    fig = figure;
    showFigure(fig, t1, b0, bb, [0,0,0], ['auto aligned']);
    set(fig, 'PaperPositionMode', 'auto');
    print(fig, '-dpng', '-r90', [outPathName,'_autoAlign.png']);
end

% Now build and save the dt6 data
eigVal = zeros([size(b0.data) 3]);
eigVec = zeros([size(b0.data) 3 3]);
for(ii=1:3)
    tmp = niftiRead(fullfile(datapath, [basename 'L' num2str(ii) '.nii.gz']));
    if(flipDti) tmp.data = flipdim(tmp.data,1); end
    eigVal(:,:,:,ii) = double(tmp.data);
end
for(ii=1:3)
    tmp = niftiRead(fullfile(datapath, [basename 'V' num2str(ii) '.nii.gz']));
    if(flipDti) tmp.data = flipdim(tmp.data,1); end
    eigVec(:,:,:,:,ii) = double(tmp.data);
    %eigVec(:,:,:,2,ii) = -eigVec(:,:,:,2,ii);
    if(flipDti) eigVec(:,:,:,1,ii) = -eigVec(:,:,:,1,ii); end
end

dt6 = dtiRebuildTensor(eigVec, eigVal);
clear eigVal eigVec;
%dt6_tmp = dtiLoadTensorElements(fullfile(fileparts(b0.baseFname), 'TensorElements.float.'));
%%dt6_tmp = permute(dt6_tmp,[2 1 3 4]);
%for(ii=1:6)
%    dt6(:,:,:,ii) = applyCannonicalXform(dt6_tmp(:,:,:,ii), b0.cannonical, b0.mmPerVox);
%end
%dt6 = permute(dt6,[2 1 3 4]);
%clear dt6_tmp;

% Reslice everything
disp('Interpolating tensors...');
%dt6A = dtiResliceTensorAffine(dt6, inv(b0.acpcXform), b0.cannonical_mmPerVox, bb, mmDt);
%dt6A = permute(dt6A,[2 1 3 4]);
if(tensorInterpSpline) tensorInterp = [7 7 7 0 0 0];
else tensorInterp = [1 1 1 0 0 0]; end
dt6 = mrAnatResliceSpm(dt6, inv(b0.acpcXform), bb, mmDt, tensorInterp, showFigs);
dt6(isnan(dt6)) = 0;
% NOTE: we want to apply b0.acpcXform to the tensors, even though we
% resliced them with inv(b0.acpcXform). inv(b0.acpcXform) maps from
% the new space to the old- the correct mapping for the interpolation,
% since we interpolate by creating a grid in the new space and fill it by
% pulling data from the old space. But the tensor reorientation should be
% done using the old-to-new space mapping (b0.acpcXform).
rigidXform = dtiFiniteStrainDecompose(b0.acpcXform);
[t,r] = affineDecompose([rigidXform,[0 0 0]';[0 0 0 1]]);
if(all(all(rigidXform==eye(3))))
    disp('No PPD correction needed- rotation is zero.');
else
    fprintf('Applying PPD rotation [%0.4f %0.4f %0.4f]...\n',r);
    dt6 = dtiXformTensors(dt6, rigidXform);
end
%figure; imagesc(makeMontage(dt6_new(:,:,:,1))); axis equal; colormap gray

disp('Interpolating B0...');
[b0_img,newXform] = mrAnatResliceSpm(b0.data, inv(b0.acpcXform), bb, mmDt, [7 7 7 0 0 0], showFigs);
b0_img(b0_img<0) = 0;
b0_img(isnan(b0_img)) = 0;

% create and apply a brain mask to the dt6 data
brainMask = mrAnatHistogramClip(b0_img, 0.4, 0.99);
brainMask = brainMask > 0.1;
% [c,v] = mrAnatHistogramSmooth(b0_img,256,0.05);
% % derivative of the thresholded derivative gives us the locations where a
% % rising or falling trend begins. In general, we can just take the first
% % rising trend (ie. just after the background noise) to find the brain.
% % We then back off from that by 25% (an empirically determined heuristic).
% peakStart = find(diff(diff(c)>0)>0);
% thresh = v(peakStart(1)+2);
% thresh = thresh.*0.75;
% brainMask = b0_img > thresh;
brainMask = dtiCleanImageMask(brainMask, 9);
figure; imagesc(makeMontage(brainMask)); axis equal tight; colormap(gray);
title('Brain Mask (based on b=0 image)');
dt6(repmat(~brainMask, [1,1,1,6])) = 0;

disp('Interpolating T1...');
anat.img = mrAnatResliceSpm(t1.data, inv(t1.acpcXform), bb, mmAnat, [], showFigs);
anat.img(anat.img<0) = 0;
anat.img(isnan(anat.img)) = 0;
anat.img = int16(round(anat.img));

mmPerVox = mmDt;
notes = 'Built from FSL data';
b0 = int16(round(b0_img));
anat.mmPerVox = mmAnat;

origin = bb(1,[1,2,3])-mmDt./2;
xformToAnat = [diag(mmDt./mmAnat) [0 0 0]'; [0 0 0 1]];
xformToAcPc = [diag(mmDt) origin'; [0 0 0 1]]; %swapXY*[diag(mm) origin'; [0 0 0 1]];
anat.xformToAcPc = [diag(mmAnat) origin'; [0 0 0 1]];
% rp = t1.talairach.refPoints;
% % Why swap coords like this? Well, the real reason has to do with the
% % history of ComputeTalairach. We could figure this out in a more general
% % way using t1.talairach.refPoints.mat, which tells us how the coords were
% % transformed in ComputeTalairach. In the end, the following reordering
% % works for our data, since the pre-processing is consistent.
% mm = t1.cannonical_mmPerVox([2 3 1]);
anat.talScale.notes = ['Scale factors to go from subject mm to Talairach mm.' ...
                       'Eg. talCoord = anat.talScale.sac * imgCoord.' ...
                       'The scales are (in order) superior, inferior, left, right ' ...
                       'and anterior of AC, ac-to-pc and posterior of pc.'];
% anat.talScale.sac = 72./sqrt(sum(((rp.acXYZ-rp.sacXYZ).*mm).^2));
% anat.talScale.iac = 42./sqrt(sum(((rp.acXYZ-rp.iacXYZ).*mm).^2));
% anat.talScale.lac = 62./sqrt(sum(((rp.acXYZ-rp.lacXYZ).*mm).^2));
% anat.talScale.rac = 62./sqrt(sum(((rp.acXYZ-rp.racXYZ).*mm).^2));
% anat.talScale.aac = 68./sqrt(sum(((rp.acXYZ-rp.aacXYZ).*mm).^2));
% anat.talScale.acpc = 24./sqrt(sum(((rp.acXYZ-rp.pcXYZ).*mm).^2));
% % The PPC is referenced to the PC. It is at Talairach (0,-102,0) and the PC
% % is at (0,-24,0), so it is 78 mm beyond the PC.
% anat.talScale.ppc = 78./sqrt(sum(((rp.pcXYZ-rp.ppcXYZ).*mm).^2));
anat.talScale.sac = 1.0;
anat.talScale.iac = 1.0;
anat.talScale.lac = 1.0;
anat.talScale.rac = 1.0;
anat.talScale.aac = 1.0;
anat.talScale.acpc = 1.0;
anat.talScale.ppc = 1.0;

disp(['Saving to ' outPathName '...']);
l = license('inuse'); h = hostid;
created = ['Created at ' datestr(now,31) ' by ' l(1).user ' on ' h{1} '.'];
save(outPathName, 'dt6', 'mmPerVox', 'notes', 'xformToAnat', 'xformToAcPc', 'b0', 'anat', 'created');
if(nargout==0)
    clear dt6;
end
return;

function showFigure(fig, t1, b0, bb, slice, figName)
if(~exist('figName','var')) figName = 'Interpolated slices'; end
% Get X,Y and Z (L-R, A-P, S-I) slices from T1 and b0 volumes
[t1Xsl] = dtiGetSlice(t1.acpcXform,t1.data,3,slice(3),bb);
[t1Ysl] = dtiGetSlice(t1.acpcXform,t1.data,2,slice(2),bb);
[t1Zsl] = dtiGetSlice(t1.acpcXform,t1.data,1,slice(1),bb);
[b0Xsl] = dtiGetSlice(b0.acpcXform,b0.data,3,slice(3),bb);
[b0Ysl] = dtiGetSlice(b0.acpcXform,b0.data,2,slice(2),bb);
[b0Zsl] = dtiGetSlice(b0.acpcXform,b0.data,1,slice(1),bb);
% Max values for image scaling
t1mv = max([t1Xsl(:); t1Ysl(:); t1Zsl(:)])+0.000001;
b0mv = max([b0Xsl(:); b0Ysl(:); b0Zsl(:)])+0.000001;
% Create XxYx3 RGB images for each of the axis slices. The green and
% blue channels are from the T1, the red channel is an average of T1
% and b=0.
Xsl(:,:,1) = t1Xsl./t1mv.*.5 + b0Xsl./b0mv.*.5;
Xsl(:,:,2) = t1Xsl./t1mv.*.5; Xsl(:,:,3) = t1Xsl./t1mv.*.5;
Ysl(:,:,1) = t1Ysl./t1mv.*.5 + b0Ysl./b0mv.*.5;
Ysl(:,:,2) = t1Ysl./t1mv.*.5; Ysl(:,:,3) = t1Ysl./t1mv.*.5;
Zsl(:,:,1) = t1Zsl./t1mv.*.5 + b0Zsl./b0mv.*.5;
Zsl(:,:,2) = t1Zsl./t1mv.*.5; Zsl(:,:,3) = t1Zsl./t1mv.*.5;

% Show T1 slices
figure(fig); set(fig, 'NumberTitle', 'off', 'Name', figName);
figure(fig); subplot(3,3,1); imagesc(bb(:,1), bb(:,2), t1Xsl);
colormap(gray); axis equal tight xy;
figure(fig); subplot(3,3,2); imagesc(bb(:,3), bb(:,1), t1Ysl);
colormap(gray); axis equal tight xy;
figure(fig); subplot(3,3,3); imagesc(bb(:,3), bb(:,2), t1Zsl);
colormap(gray); axis equal tight xy;
axis equal tight;

% Show b=0 slices
figure(fig); subplot(3,3,4); imagesc(bb(:,1), bb(:,2), b0Xsl);
colormap(gray); axis equal tight xy;
figure(fig); subplot(3,3,5); imagesc(bb(:,3), bb(:,1), b0Ysl);
colormap(gray); axis equal tight xy;
figure(fig); subplot(3,3,6); imagesc(bb(:,3), bb(:,2), b0Zsl);
colormap(gray); axis equal tight xy;
axis equal tight;

% Show combined slices
figure(fig); subplot(3,3,7); imagesc(bb(:,1), bb(:,2), Xsl);
axis equal tight xy;
figure(fig); subplot(3,3,8); imagesc(bb(:,3), bb(:,1), Ysl);
axis equal tight xy;
figure(fig); subplot(3,3,9); imagesc(bb(:,3), bb(:,2), Zsl);
axis equal tight xy;
axis equal tight;

return;