text:
function f = allenAtlasBrowser(f, templateVolume, annotationVolume, structureTree, slice_figure, save_location, save_suffix, plane)
% ------------------------------------------------
% Browser for the allen atlas ccf data in matlab.
% ------------------------------------------------
%
% Inputs templateVolume, annotationVolume, and structureTree are the data describing the atlas.
% The annotation volume should be the "by_index" version
%

% print instructions
display_controls(plane)

% create figure and adjust to user's screen size
% f = figure('Name','Atlas Viewer');
figure(f);
try; screen_size = get(0,'ScreenSize'); screen_size = [max(screen_size(3:4)) min(screen_size(3:4))]./[2560 1440];
catch; screen_size = [1900 1080]./[2560 1440];
end
set(f,'Position', [1050*screen_size(1) 660*screen_size(2) 880*screen_size(1) 650*screen_size(2)])
movegui(f,'onscreen')

% initialize user data variables held by the figure
ud.plane = plane;
ud.bregma = allenCCFbregma;
ud.currentSlice = ud.bregma(1);
ud.currentAngle = zeros(2,1);
ud.scrollMode = 0;

ud.transform_type = 'projective'; %can change to 'affine' or 'pwl'

ud.oldContour = [];
ud.showContour = false;
ud.showOverlay = false; ud.overlayAx = [];
ud.pointList = cell(1,3); ud.pointList{1} = zeros(0,3);
ud.pointHands = cell(1,3);
ud.probe_view_mode = false;
ud.currentProbe = 0; ud.ProbeColors = [1 1 1; 1 .75 0;  .3 1 1; .4 .6 .2; 1 .35 .65; .7 .7 1; .65 .4 .25; .7 .95 .3; .7 0 0; .5 0 .6; 1 .6 0];
ud.ProbeColor =  {'white','gold','turquoise','fern','bubble gum','overcast sky','rawhide', 'green apple','red','purple','orange'};
ud.getPoint_for_transform =false; ud.pointList_for_transform = zeros(0,2); ud.pointHands_for_transform = [];
ud.current_pointList_for_transform = zeros(0,2); ud.curr_slice_num = 1;
ud.clicked = false;
ud.showAtlas = false;
ud.viewColorAtlas = false;
ud.histology_overlay = 0;
ud.atlasAx = axes('Position', [0.05 0.05 0.9 0.9]);
ud.transform = [];
ud.transformed_slice_figure = [];
ud.slice_shift = 0;
ud.loaded_slice = 0;
ud.slice_at_shift_start = 1;
ud.text = [];

reference_image = squeeze(templateVolume(ud.currentSlice,:,:));
ud.im = plotTVslice(reference_image);
ud.ref_size = size(reference_image);
ud.ref = uint8(squeeze(templateVolume(ud.currentSlice,:,:)));
ud.curr_im = uint8(squeeze(templateVolume(ud.currentSlice,:,:)));
ud.curr_slice_trans = uint8(squeeze(templateVolume(ud.currentSlice,:,:)));
ud.im_annotation = squeeze(annotationVolume(ud.currentSlice,:,:));
ud.atlas_boundaries = zeros(ud.ref_size,'uint16');;
ud.offset_map = zeros(ud.ref_size);
ud.loaded = 0;

% create functions needed to interact with the figure
set(ud.im, 'ButtonDownFcn', @(f,k)atlasClickCallback(f, k, slice_figure, save_location));
ud.bregmaText = annotation('textbox', [0 0.95 0.4 0.05], ...
    'String', '[coords]', 'EdgeColor', 'none', 'Color', 'k');

ud.angleText = annotation('textbox', [.7 0.95 0.4 0.05], ...
    'EdgeColor', 'none', 'Color', 'k');

allData.tv = templateVolume;
allData.av = annotationVolume;
allData.st = structureTree;
hold(ud.atlasAx, 'on');
set(ud.atlasAx, 'HitTest', 'off');
set(f, 'UserData', ud);
set(f, 'KeyPressFcn', @(f,k)hotkeyFcn(f, slice_figure, k, allData, save_location, save_suffix));
set(f, 'WindowScrollWheelFcn', @(src,evt)updateSlice(f, evt, allData, slice_figure, save_location))
set(f, 'WindowButtonMotionFcn',@(f,k)fh_wbmfcn(f, allData, slice_figure, save_location)); % Set the motion detector.

% display user controls in the console
function display_controls(plane)
fprintf(1, '\n Controls: \n');
fprintf(1, '--------- \n');
fprintf(1, 'Navigation: \n');
switch plane
    case 'coronal'
        fprintf(1, 'up: scroll through D/V angles (for coronal sections)\n');
        fprintf(1, 'right: scroll through M/L angles  (for coronal sections)\n');
        fprintf(1, 'down: scroll through A/P atlas slices \n');
        fprintf(1, 'left: scroll for Slice Viewer \n');

        fprintf(1, 'scroll: move between slices or angles \n');
    case 'sagittal'
        fprintf(1, 'up: scroll through D/V angles (for sagittal sections)\n');
        fprintf(1, 'right: scroll through A/P angles  (for sagittal sections)\n');
        fprintf(1, 'down: scroll through M/L atlas slices \n');
        fprintf(1, 'left: scroll for Slice Viewer \n');

        fprintf(1, 'scroll: move between slices or angles \n');

    case 'transverse'
        fprintf(1, 'up: scroll through M/L angles (for transverse sections)\n');
        fprintf(1, 'right: scroll through A/P angles  (for transverse sections)\n');
        fprintf(1, 'down: scroll through D/V atlas slices \n');
        fprintf(1, 'left: scroll for Slice Viewer \n');

        fprintf(1, 'scroll: move between slices or angles \n');

end

fprintf(1, '\n Registration: \n');
fprintf(1, 't: toggle mode where clicks are logged for transform \n');
fprintf(1, 'h: toggle overlay of current histology slice \n');
fprintf(1, 'p: toggle mode where clicks are logged for probe or switch probes \n');
fprintf(1, 'n: add a new probe \n');
fprintf(1, 'x: save transform and current atlas location \n');
fprintf(1, 'l: load transform for current slice; press again to load probe points \n');
fprintf(1, 's: save current probe \n');
fprintf(1, 'd: delete most recent probe point or transform point \n');
fprintf(1, 'w: enable/disable probe viewer mode for current probe  \n');

fprintf(1, '\n Viewing modes: \n');
fprintf(1, 'o: toggle overlay of current region extent \n');
fprintf(1, 'a: toggle to viewing boundaries \n');
fprintf(1, 'v: toggle to color atlas mode \n');
fprintf(1, 'g: toggle gridlines \n');

fprintf(1, '\n space: display controls \n \n');

% ------------------------
% react to keyboard press
% ------------------------
function hotkeyFcn(f, slice_figure, keydata, allData, save_location, save_suffix)

% retrieve user data from figure
ud = get(f, 'UserData');
ud_slice = get(slice_figure, 'UserData');

key_letter = lower(keydata.Key);
switch key_letter
% space -- display controls
    case 'space'
        display_controls(ud.plane)
% o -- toggle showing brain region overlay
    case 'o'
        ud.showOverlay = ~ud.showOverlay;
        if ~ud.showOverlay
            delete(ud.overlayAx); ud.overlayAx = [];
            disp('Overlay OFF');
        elseif ~ud.viewColorAtlas; disp('Overlay on!');
        end
% g -- toggle showing Gridlines
    case 'g'
        if ~isfield(ud, 'gridlines') || isempty(ud.gridlines)
            axes(ud.atlasAx); hold on;
            gridY = 100:100:ud.ref_size(1); % assuming the size of the atlas for this for now
            gridX = 70:100:ud.ref_size(2);
            xl = xlim(); yl = ylim();
            gx = arrayfun(@(x)plot(x*[1 1], yl, 'w'), gridX, 'uni', false);
            gy = arrayfun(@(y)plot(xl, y*[1 1], 'w'), gridY, 'uni', false);
            ud.gridlines = [gx gy];
        elseif strcmp(get(ud.gridlines{1}, 'Visible'), 'on');
            cellfun(@(x)set(x, 'Visible', 'off'), ud.gridlines);
        elseif strcmp(get(ud.gridlines{1}, 'Visible'), 'off');
            cellfun(@(x)set(x, 'Visible', 'on'), ud.gridlines);
        end
% p -- toggle mode to register clicks as probe points
    case 'p'
        ud.probe_view_mode = 0;
        if ud.currentProbe < size(ud.pointList,1)
            ud.currentProbe = ud.currentProbe + 1;

            disp(['probe point mode -- selecting probe ' num2str(ud.currentProbe) ' (' ud.ProbeColor{ud.currentProbe} ')']);
            ud.getPoint_for_transform = false;

            % show Transformed Slice & Probage Viewer, if not already showing
            if ~ud.slice_at_shift_start; add = 1; else; add = 0; end
                slice_name = ud_slice.processed_image_names{ud.slice_at_shift_start+ud.slice_shift+add}(1:end-4);
                folder_transformations = fullfile(save_location, ['transformations' filesep]);
                try; load([folder_transformations slice_name '_transform_data.mat']);
                    try; figure(ud.transformed_slice_figure);
                    catch

                        ud.transformed_slice_figure = figure('Name','Transformed Slice & Probe Point Viewer');
                        try; screen_size = get(0,'ScreenSize'); screen_size = [max(screen_size(3:4)) min(screen_size(3:4))]./[2560 1440];
                        catch; screen_size = [1900 1080]./[2560 1440];
                        end

                        set(ud.transformed_slice_figure,'Position', [256*screen_size(1) 12*screen_size(2) 560*screen_size(1) 420*screen_size(2)]);
                        movegui(ud.transformed_slice_figure,'onscreen')

                        highlight_point = false;
                        transformed_sliceBrowser(ud.transformed_slice_figure, save_location, f, highlight_point, [], [], [], [], [], add);
                    end; figure(f);
                end
        else
            ud.currentProbe = 0;
            disp(['probe point mode OFF']);
        end
% w -- toggle mode to visualize probe trajectory
    case 'w'
        % remove color atlas
        ud.viewColorAtlas = false;
        set(ud.im, 'CData', ud.ref)
        colormap(ud.atlasAx, 'gray'); caxis(ud.atlasAx, [0 400]);

        % remove overlay
        ud.showOverlay = 0;
        delete(ud.overlayAx); ud.overlayAx = [];

        for probe_plotted = 1:size(ud.pointHands,1)
            set(ud.pointHands{probe_plotted,1}(:), 'Visible', 'off');
            set(ud.pointHands{probe_plotted,3}(:), 'Visible', 'off');
        end
        if ~ud.currentProbe
            ud.currentProbe = 1;
        end
            ud.probe_view_mode = ~ud.probe_view_mode;

            if ud.probe_view_mode && ~isempty(ud.pointList{ud.currentProbe,1})
                 % load probe points if none are already up
                if ~size(ud.pointList{1,1},1)
                    probe_points = load(fullfile(save_location, ['probe_points' save_suffix]));  disp('loading probe points')
                    ud.pointList = probe_points.pointList.pointList;
                    ud.pointHands = probe_points.pointList.pointHands;
                for probe = 1:size(ud.pointList,1)
                    for probe_point = 1:size(ud.pointList{probe, 1},1)
                        ud.pointHands{probe, 1}(probe_point) = scatter(ud.atlasAx, ...
                            ud.pointList{probe,1}(probe_point,1), ud.pointList{probe,1}(probe_point,2), 20, 'ro', ...
                        'MarkerFaceColor', ud.ProbeColors(probe, :),'MarkerEdgeColor', ud.ProbeColors(probe, :), ...
                        'LineWidth',3);
                    end
                    if probe ~= ud.currentProbe
                        set(ud.pointHands{probe, 1}(:), 'Visible', 'off');
                    end
                end

                end; disp(['activate probe view mode for probe ' num2str(ud.currentProbe)])

                % show probe points
                set(ud.pointHands{ud.currentProbe, 1}(:), 'Visible', 'on');

                curr_probe_handle = ud.pointHands{ud.currentProbe,1};
                set(curr_probe_handle, 'ButtonDownFcn', @(f,k)atlasClickCallback(f, k, slice_figure, save_location));

                % switch the order to z y x
                curr_probePoints = ud.pointList{ud.currentProbe,1}(:, [3 2 1]);

                % get line of best fit through points
                % m is the mean value of each dimension; p is the eigenvector for largest eigenvalue
                [m,p,s] = best_fit_line(curr_probePoints(:,1), curr_probePoints(:,2), curr_probePoints(:,3));
                atlas_res = 0.010;
                disp(['root mean square error to line of ' num2str(round(sqrt(s / size(curr_probePoints,1))*(atlas_res*1000),1)) ' micron'])

                min_y = min(ud.pointList{ud.currentProbe,1}(:,2));
                max_y = max(ud.pointList{ud.currentProbe,1}(:,2));
                min_x = m(3) + (min_y - m(2))  * p(3) / p(2);
                max_x = m(3) + (max_y - m(2))  * p(3) / p(2);
                max_z = m(1) + (max_y - m(2))  * p(1) / p(2);
                ud.pointHands{ud.currentProbe, 3} = plot([min_x max_x],[min_y max_y],'color',ud.ProbeColors(ud.currentProbe,:),'linestyle',':','linewidth',1.5);
                set(ud.pointHands{ud.currentProbe,3}, 'ButtonDownFcn', @(f,k)atlasClickCallback(f, k, slice_figure, save_location));

                % ensure proper orientation: want 0 at the top of the brain
                % and positive distance goes down into the brain
                if p(2)<0
                    p = -p;
                end

                % calculate slice angle along probe track -- do so by
                % constraining either ML angle or DV angle to zero
                position_at_x_center_point = m + (ud.ref_size(2)/2 - m(3)) * p / p(3);
                position_at_y_center_point = m + (ud.ref_size(1)/2 - m(2)) * p / p(2);

                angle_DV_if_constraining_ML = round(ud.ref_size(1)/2 / (ud.ref_size(1)/2 - m(2)) * (position_at_y_center_point(1)-m(1)) );
                angle_ML_if_constraining_DV = round(ud.ref_size(2)/2 / (ud.ref_size(2)/2 - m(3)) * (position_at_x_center_point(1)-m(1)) );

                if abs(angle_DV_if_constraining_ML) < abs(angle_ML_if_constraining_DV) || abs(angle_DV_if_constraining_ML) < 150
                    ud.currentSlice = round(position_at_y_center_point(1));
                    ud.currentAngle(1) = angle_DV_if_constraining_ML;
                    ud.currentAngle(2) = 0;
                else
                    ud.currentSlice = round(position_at_x_center_point(1));
                    ud.currentAngle(2) = angle_ML_if_constraining_DV;
                    ud.currentAngle(1) = 0;
                end

                % if the angles are too extreme, just show mean slice without angles
                if abs(ud.currentAngle(1))+abs(ud.currentAngle(2)) >= ud.currentSlice || ...
                        abs(ud.currentAngle(1))+abs(ud.currentAngle(2)) >= ud.ref_size(2)-ud.currentSlice
                    ud.currentAngle(1) = 0;
                    ud.currentAngle(2) = 0;
                    ud.currentSlice = round(m(1));
                    disp('probe angle not ideal for viewing in coronal slice -- angles set to 0')
                end

                % report estimated probe angle
                AP_angle = round(atand(angle_DV_if_constraining_ML/(ud.ref_size(1)/2)),1);
                ML_angle = round(atand((max_x - min_x)/(max_y - min_y)),1);
                position_at_bregma_depth = [ (m(3) + (0 - m(2)) * p(3) / p(2)) ((m(1) + (0 - m(2)) * p(1) / p(2)))];
                ML_position = round((position_at_bregma_depth(1)-ud.bregma(3))*atlas_res,2);
                AP_position = round((ud.bregma(1) - position_at_bregma_depth(2))*atlas_res,2);
                disp(' ');
                disp('---estimated probe insertion---')
                disp(['entry position at DV = 0: AP = ' num2str(AP_position) ' mm, ' ...
                                                'ML = ' num2str(ML_position) ' mm'])
                insertion_dist = sqrt((position_at_bregma_depth(1) - max_x)^2+(position_at_bregma_depth(2)-max_z)^2+(0-max_y)^2);
                disp(['insertion distance from the above position: ' num2str(round(insertion_dist*atlas_res,3)) ' mm'])
                direction_AP = {'posterior','anterior'};
                disp([num2str(abs(AP_angle)) ' degrees in the ' direction_AP{(AP_angle<0)+1} ' direction'])
                direction_ML = {'medial','lateral'};
                disp([num2str(abs(ML_angle)) ' degrees in the ' direction_ML{(ML_angle<0&ML_position<0)+1} ' direction']);
                disp(' ');

                % update slice
                update.VerticalScrollCount = 0; ud.scrollMode = 0; ud.histology_overlay = 0; set(f, 'UserData', ud);
                updateSlice(f, update, allData, slice_figure, save_location); ud = get(f, 'UserData');
                fill([5 5 250 250],[5 50 50 5],[0 0 0]);

                % show Transformed Slice & Probage Viewer, if not already showing
                try; figure(ud.transformed_slice_figure);
                catch
                    ud.transformed_slice_figure = figure('Name','Transformed Slice & Probe Point Viewer');
                    try; screen_size = get(0,'ScreenSize'); screen_size = screen_size(3:4)./[2560 1440];
                    catch; screen_size = [1900 1080]./[2560 1440];
                    end

                    set(ud.transformed_slice_figure,'Position', [256*screen_size(1) 37*screen_size(2) 560*screen_size(1) 420*screen_size(2)]);
                    movegui(ud.transformed_slice_figure,'onscreen');
                    highlight_point = false;
                    transformed_sliceBrowser(ud.transformed_slice_figure, save_location, f, highlight_point, [], [], [], [], [],0);
                end; figure(f);

                set(ud.im, 'CData', ud.ref);
                ud.curr_im = ud.ref; set(f, 'UserData', ud);
            else; disp('probe view mode OFF'); end

% t -- toggle mode to register clicks as Points
    case 't'
        ud.getPoint_for_transform = ~ud.getPoint_for_transform;
        ud.loaded = false;

        if ud.getPoint_for_transform
            disp('transform point mode on');

            ud.currentProbe = 0;

            % launch transform point mode
            if ~size(ud.current_pointList_for_transform,1)% ) (ud.curr_slice_num ~= (ud.slice_at_shift_start+ud.slice_shift) ||  && ~ud.loaded
                ud.curr_slice_num = ud.slice_at_shift_start+ud.slice_shift; %ud_slice.slice_num;
                ud.current_pointList_for_transform = zeros(0,2);
                set(ud.pointHands_for_transform(:), 'Visible', 'off');
                num_hist_points = size(ud_slice.pointList,1);
                template_point = 1; template_points_shown = 0;
                updateBoundaries(f,ud, allData); ud = get(f, 'UserData');
            end
        else; disp('transform point mode OFF');
        end
% a -- toggle viewing of annotation boundaries
    case 'a'

        ud.showAtlas = ~ud.showAtlas;

        if ud.showAtlas % superimpose boundaries
            updateBoundaries(f,ud, allData);
        else % return to image
            set(ud.im, 'CData', ud.curr_im)
        end

% v -- toggle viewing of Color Atlas
    case 'v'
        ud.viewColorAtlas = ~ud.viewColorAtlas;
        ud.histology_overlay = 0;

        if ud.viewColorAtlas
            % remove overlay
            ud.showOverlay = 0;
            ref_mode = false;
            delete(ud.overlayAx); ud.overlayAx = [];
            set(ud.im, 'CData', ud.im_annotation)
            ud.curr_im = ud.im_annotation;
            cmap = allen_ccf_colormap('2017');
            colormap(ud.atlasAx, cmap); caxis(ud.atlasAx, [1 size(cmap,1)]);
            fill([5 5 250 250],[5 50 50 5],[1 1 1],'edgecolor',[1 1 1]);
        else
            set(ud.im, 'CData', ud.ref)
            ud.curr_im = ud.ref;
            colormap(ud.atlasAx, 'gray'); caxis(ud.atlasAx, [0 400]);
            fill([5 5 250 250],[5 50 50 5],[0 0 0]);
        end

        if ud.showAtlas % superimpose boundaries
            updateBoundaries(f,ud, allData);
        end

    case 'uparrow' % scroll angles along M/L axis
        ud.scrollMode = 1;
        if strcmp(ud.plane,'coronal')
            disp('switch scroll mode -- tilt D/V')
        elseif strcmp(ud.plane,'sagittal')
            disp('switch scroll mode -- tilt D/V')
        elseif strcmp(ud.plane,'transverse')
            disp('switch scroll mode -- tilt M/L')
        end
    case 'rightarrow' % scroll angles along A/P axis
        ud.scrollMode = 2;
        if strcmp(ud.plane,'coronal')
            disp('switch scroll mode -- tilt M/L')
        elseif strcmp(ud.plane,'sagittal')
            disp('switch scroll mode -- tilt A/P')
        elseif strcmp(ud.plane,'transverse')
            disp('switch scroll mode -- tilt A/P')
        end

    case 'downarrow' % scroll along A/P axis
        ud.scrollMode = 0;
        if strcmp(ud.plane,'coronal')
            disp('switch scroll mode -- scroll along A/P axis')
        elseif strcmp(ud.plane,'sagittal')
            disp('switch scroll mode -- scroll along M/L axis')
        elseif strcmp(ud.plane,'transverse')
            disp('switch scroll mode -- scroll along D/V axis')
        end
    case 'leftarrow' % scroll along A/P axis
        ud.scrollMode = 3;
        if ~ud.slice_at_shift_start
            ud.slice_at_shift_start = ud_slice.slice_num;
            ud.slice_shift = 0;
        end
        disp('switch scroll mode -- scroll along slice images')
% h -- toggle viewing of histology / histology overlay
    case 'h'
        disp('  ');
        % remove color atlas
        ud.viewColorAtlas = false;
        set(ud.im, 'CData', ud.ref)
        colormap(ud.atlasAx, 'gray'); caxis(ud.atlasAx, [0 400]);
        % remove overlay
        ud.showOverlay = 0;
        ref_mode = false;
        delete(ud.overlayAx); ud.overlayAx = [];
        % toggle which mode is active
        if ~ud.clicked || ~ud.histology_overlay;
            ud.histology_overlay = ud.histology_overlay + 1 - 3*(ud.histology_overlay==2);
        end
        ud.clicked = false;
        % get clicked points and slice info from the slice figure
        slice_points = ud_slice.pointList;
        slice_name = ud_slice.processed_image_names{ud.slice_at_shift_start+ud.slice_shift}(1:end-4);
        folder_transformations = [save_location filesep 'transformations' filesep];
        if size(ud.current_pointList_for_transform,1)  && size(slice_points,1) && ud.slice_at_shift_start+ud.slice_shift == ud_slice.slice_num && ~ud.probe_view_mode
            key_letter = 'x'; % save transform automatically
        end
        % if in one of the transformation modes, perform the transform
        if (ud.histology_overlay == 1 || ud.histology_overlay == 2) && ...
                ( (size(ud.current_pointList_for_transform,1) && size(slice_points,1)) || ud.loaded)
            try
                if ud.slice_shift > 0
                    ud.curr_slice_trans = imread([folder_transformations slice_name '_transformed.tif']);
                else
                    set(ud.text,'Visible','off');
                    fill([5 5 250 250],[5 50 50 5],[0 0 0]); ud.text(end+1) = text(5,15,['Slice ' num2str(ud.slice_at_shift_start+ud.slice_shift)],'color','white');

                    reference_points = ud.current_pointList_for_transform;
                    slice_points = ud_slice.pointList;

                    current_slice_image = flip(get(ud_slice.im, 'CData'),1);

                    % ** this is where the transform happens, using the
                    % clicked points from the reference and histology images
                    if ~ud.loaded  % use loaded version if 'l' was just pressed
                        ud.transform = fitgeotrans(slice_points,reference_points,ud.transform_type); %can use 'affine', 'projective', 'polynomial', or 'pwl'
                    end
                    R = imref2d(size(ud.ref));
                    ud.curr_slice_trans = imwarp(current_slice_image, ud.transform, 'OutputView',R);
                end

                image_blend =  imfuse(uint8(ud.ref*.6), ud.curr_slice_trans(:,:,:),'blend','Scaling','none');
                if ud.histology_overlay == 2 % 2 ~ blend
                    disp('Slice + Reference mode!');
                    set(ud.im, 'CData', image_blend);
                    ud.curr_im = image_blend;
                else % 1 ~ just see slice
                    disp('Slice mode!');
                    set(ud.im, 'CData', ud.curr_slice_trans);
                    ud.curr_im = ud.curr_slice_trans;
                end
            % if wrong number of points clicked
            catch
                ref_mode = true;
                disp(['Unable to transform -- ' num2str(size(ud_slice.pointList,1)) ...
                     ' slice points and ' num2str(size(ud.current_pointList_for_transform,1)) ' reference points']);
                key_letter = 'h';
            end
        end
        % if not doing transform, just show reference atlas
        if ud.histology_overlay == 0 || ref_mode
            ud.histology_overlay = 0;
            disp('Reference mode!');
            set(ud.im, 'CData', ud.ref);
            ud.curr_im = ud.ref; set(f, 'UserData', ud);
        end
        if ud.showAtlas
            updateBoundaries(f,ud, allData);
        end
% n -- start marking a new probe
    case 'n'
        new_num_probes = size(ud.pointList,1) + 1;
        ud.getPoint_for_transform = false;
        if new_num_probes <= size(ud.ProbeColors,1)
            disp(['probe ' num2str(new_num_probes) ' added! (' ud.ProbeColor{new_num_probes} ')']);
%         ud.probe_view_mode = 0;
        probe_point_list = cell(new_num_probes,1); probe_hands_list = cell(new_num_probes,3);
        for prev_probe = 1:new_num_probes-1
            probe_point_list{prev_probe,1} = ud.pointList{prev_probe,1};
            probe_point_list{prev_probe,2} = ud.pointList{prev_probe,2};
            probe_point_list{prev_probe,3} = ud.pointList{prev_probe,3};
            probe_hands_list{prev_probe, 1} = ud.pointHands{prev_probe, 1};
            probe_hands_list{prev_probe, 2} = ud.pointHands{prev_probe, 2};
            probe_hands_list{prev_probe, 3} = ud.pointHands{prev_probe, 3};
        end; probe_point_list{new_num_probes,1} = zeros(0,3);
        ud.pointList = probe_point_list; ud.pointHands = probe_hands_list;
        ud.currentProbe = new_num_probes;
        end
 % s -- save probe trajectory and points of each probe per histology image (and associated histology name/number)
    case 's'
        pointList.pointList = ud.pointList;
        for probe = 1:size(ud.pointList,1)
            set(ud.pointHands{probe,3}, 'ButtonDownFcn', '');
            set(ud.pointHands{probe,1}, 'ButtonDownFcn', '');
        end
        pointList.pointHands = ud.pointHands;
        warning('off', 'MATLAB:Figure:FigureSavedToMATFile');
        save(fullfile(save_location, ['probe_points' save_suffix]), 'pointList'); disp('probe points saved');
% l -- load transform and current slice position and angle
    case 'l'
        slice_name = ud_slice.processed_image_names{ud_slice.slice_num}(1:end-4);
        folder_transformations = fullfile(save_location, ['transformations' filesep]);
        ud.clicked = false;
        try
        if ud.loaded_slice+ud.slice_shift ~= ud_slice.slice_num

            ud.curr_slice_num = ud_slice.slice_num;

            % remove overlay
            ud.showOverlay = 0;
            delete(ud.overlayAx); ud.overlayAx = [];

            % load transform data
            transform_data = load(fullfile(folder_transformations, [slice_name '_transform_data.mat']));
            transform_data = transform_data.save_transform;

            % load new transform
            ud.transform = transform_data.transform;

            if ~isempty(transform_data.transform_points{1}) && ~isempty(transform_data.transform_points{2})
                ud.current_pointList_for_transform = transform_data.transform_points{1};
                ud_slice.pointList = transform_data.transform_points{2};
                set(slice_figure, 'UserData', ud_slice);
            end

            % load allen ref location
            ud.currentSlice = transform_data.allen_location{1}; ud.currentAngle = transform_data.allen_location{2};

            % create transformed histology image
            current_slice_image = flip(get(ud_slice.im, 'CData'),1); R = imref2d(size(ud.ref));
            ud.curr_slice_trans = imwarp(current_slice_image, ud.transform, 'OutputView',R);

            % update figure
            update.VerticalScrollCount = 0; temp_scroll_mode = ud.scrollMode; ud.scrollMode = 4; set(f, 'UserData', ud);
            updateSlice(f, update, allData, slice_figure, save_location); ud = get(f, 'UserData');
            ud.scrollMode = temp_scroll_mode;
            ud.loaded = true;

            ud.slice_at_shift_start = ud_slice.slice_num; ud.slice_shift = 0;
            ud.loaded_slice = ud_slice.slice_num;

            if ~isempty(ud.text)
                set(ud.text,'Visible','off');
                fill([5 5 250 250],[5 50 50 5],[0 0 0]); ud.text(end+1) = text(5,15,['Slice ' num2str(ud.slice_at_shift_start+ud.slice_shift)],'color','white');
            end

            disp('transform loaded -- press ''l'' again now to load probe points');
        else % load probe points
            if ~size(ud.pointList{1,1},1)
                probe_points = load(fullfile(save_location, ['probe_points' save_suffix]));  disp('probe points loaded')
                ud.pointList = probe_points.pointList.pointList;
                ud.pointHands = probe_points.pointList.pointHands;
            else
                disp('probe points not loaded -- there are already some current probe points')
            end
        end

        for probe = 1:size(ud.pointList,1)
            % create point plot handles anew
            try; set(ud.pointHands{probe, 1}(:),'Visible','off'); end
            for probe_point = 1:size(ud.pointList{probe, 1},1)
                ud.pointHands{probe, 1}(probe_point) = scatter(ud.atlasAx, ...
                    ud.pointList{probe,1}(probe_point,1), ud.pointList{probe,1}(probe_point,2), 20, 'ro', ...
                'MarkerFaceColor', [.1 .1 .1],'MarkerEdgeColor', ud.ProbeColors(probe, :), ...
                'LineWidth',2);

            % set the point plot from the current slice visible
            slice_point_belongs_to = ud.pointHands{probe, 2}(probe_point);
            if slice_point_belongs_to == ud_slice.slice_num
                set(ud.pointHands{probe, 1}(probe_point), 'Visible', 'on');
            else
                set(ud.pointHands{probe, 1}(probe_point), 'Visible', 'off');
            end

            end
            % turn off best fit line if applicable
            set(ud.pointHands{probe, 3}(:),'Visible','off'); ud.pointHands{probe, 3} = [];
        end

        ud.slice_shift = 0;
        catch;
            disp(['loading failed']); end
% d -- delete current transform or most recent probe point
    case 'd'
        if ud.getPoint_for_transform
%             ud.current_pointList_for_transform = zeros(0,2); set(ud.pointHands_for_transform(:), 'Visible', 'off');
%             ud.pointHands_for_transform = []; ud_slice.pointList = []; set(slice_figure, 'UserData', ud_slice);
%             disp('current transform erased');

            % Try to delete only the most recent point
            ud.current_pointList_for_transform = ud.current_pointList_for_transform(1:end-1,:);
            ud_slice.pointList = ud_slice.pointList(1:end-1,:);
            set(slice_figure, 'UserData', ud_slice);
            if ud.pointHands_for_transform
                % remove circle for most revent point
                set(ud.pointHands_for_transform(end), 'Visible', 'off');
                ud.pointHands_for_transform = ud.pointHands_for_transform(1:end-1);
                if ud.pointHands_for_transform
                    % recolor points
                    set(ud.pointHands_for_transform(end), 'color', [0 .9 0]);
                end
            end

            disp('transform point deleted');

        elseif ud.currentProbe
            ud.pointList{ud.currentProbe,1} = ud.pointList{ud.currentProbe,1}(1:end-1,:);
            ud.pointList{ud.currentProbe,2} = ud.pointList{ud.currentProbe,2}(1:end-1,:);
            ud.pointList{ud.currentProbe,3} = ud.pointList{ud.currentProbe,3}(:,1:end-1);
            set(ud.pointHands{ud.currentProbe, 1}(end), 'Visible', 'off');
            ud.pointHands{ud.currentProbe, 1} = ud.pointHands{ud.currentProbe, 1}(1:end-1);
            ud.pointHands{ud.currentProbe, 2} = ud.pointHands{ud.currentProbe, 2}(1:end-1);
            disp('probe point deleted')
        end
end
% x -- save transform and current slice position and angle
if strcmp(key_letter,'x')

        % find or create folder location for transformations
        try
        folder_transformations = fullfile(save_location, ['transformations' filesep]);
        if ~exist(folder_transformations)
            mkdir(folder_transformations)
        end

        % find name of slice
        if ud.slice_shift || ud.slice_at_shift_start ~= ud_slice.slice_num
            slice_name = ud_slice.processed_image_names{(ud.slice_at_shift_start+ud.slice_shift)}(1:end-4);
        else
            slice_name = ud_slice.processed_image_names{ud_slice.slice_num}(1:end-4);
        end

        if isempty(ud.current_pointList_for_transform)
            ud.transform = [];
        end
        % store transform, if applicable
        save_transform.transform = ud.transform;

        % store transform points
        transform_points = cell(2,1); transform_points{1} = ud.current_pointList_for_transform;
        if ~isempty(ud_slice.pointList)
            transform_points{2} = ud_slice.pointList;
        end

        save_transform.transform_points = transform_points;
        % store reference location
        allen_location = cell(2,1); allen_location{1} = ud.currentSlice; allen_location{2} = ud.currentAngle;
        save_transform.allen_location = allen_location;
        % save all this
        save(fullfile(folder_transformations, [slice_name '_transform_data.mat']), 'save_transform');
        disp('atlas location saved')

        % save transformed histology image
        current_slice_image = imread(fullfile(save_location, [slice_name '.tif']));
%         current_slice_image = flip(get(ud_slice.im, 'CData'),1);
        R = imref2d(size(ud.ref));
        curr_slice_trans = imwarp(current_slice_image, ud.transform, 'OutputView',R);
        imwrite(curr_slice_trans, fullfile(folder_transformations, [slice_name '_transformed.tif']))

        disp('transform saved')
        catch
            disp('transform not saved')
        end
end

set(f, 'UserData', ud);

% -----------------------------------------
% Update slice (from scrolling or loading)
% -----------------------------------------
function updateSlice(f, evt, allData, slice_figure, save_location)

ud = get(f, 'UserData');

% scroll through slices
if ud.scrollMode==0
    ud.currentSlice = ud.currentSlice+evt.VerticalScrollCount*3;

    if ud.currentSlice>size(allData.tv,1); ud.currentSlice = 1; end %wrap around
    if ud.currentSlice<1; ud.currentSlice = size(allData.tv,1); end %wrap around

% scroll through D/V angles
elseif ud.scrollMode==1 %&&  abs(ud.currentAngle(1)) < 130
  ud.currentAngle(1) = ud.currentAngle(1)+evt.VerticalScrollCount*3;

% scroll through M/L angles
elseif ud.scrollMode==2 %&&  abs(ud.currentAngle(2)) < 130
  ud.currentAngle(2) = ud.currentAngle(2)+evt.VerticalScrollCount*3;

% scroll through slices (left arrow pressed)
elseif ud.scrollMode == 3
  set(ud.pointHands_for_transform(:), 'Visible', 'off');
  ud.showOverlay = 0;
  delete(ud.overlayAx); ud.overlayAx = [];
  ud_slice = get(slice_figure, 'UserData');

  try
    ud.slice_shift = ud.slice_shift-evt.VerticalScrollCount;
    slice_name = ud_slice.processed_image_names{ud.slice_at_shift_start+ud.slice_shift}(1:end-4);
  catch
    ud.slice_shift = ud.slice_shift+evt.VerticalScrollCount;
    slice_name = ud_slice.processed_image_names{ud.slice_at_shift_start+ud.slice_shift}(1:end-4);
  end
  folder_transformations = fullfile(save_location, ['transformations' filesep]);

    % set probe points from other slices invisible and from this slice visible
    for probe = 1:size(ud.pointList,1)
        set(ud.pointHands{probe, 3}(:),'Visible','off'); ud.pointHands{probe, 3} = [];
        for probe_point = 1:size(ud.pointList{probe,1},1)
            slice_point_belongs_to = ud.pointList{probe, 2}(probe_point);
            if slice_point_belongs_to == ud.slice_at_shift_start+ud.slice_shift
                set(ud.pointHands{probe, 1}(probe_point), 'Visible', 'on');
            else
                set(ud.pointHands{probe, 1}(probe_point), 'Visible', 'off');
            end

        end
    end

    ud.current_pointList_for_transform = zeros(0,2);
    try; load([folder_transformations slice_name '_transform_data.mat']);

        ud.clicked = false;

        % load transform data
        transform_data = load(fullfile(folder_transformations, [slice_name '_transform_data.mat']));
        transform_data = transform_data.save_transform;

        % load new transform
        ud.transform = transform_data.transform;

        if ~isempty(transform_data.transform_points{1}) && ~isempty(transform_data.transform_points{2})
            ud.current_pointList_for_transform = transform_data.transform_points{1};
            ud_slice.pointList = transform_data.transform_points{2};
        else
            ud_slice.pointList = [];
        end
        set(slice_figure, 'UserData', ud_slice);

        % load allen ref location
        ud.currentSlice = transform_data.allen_location{1}; ud.currentAngle = transform_data.allen_location{2};

        % create transformed histology image
        ud.curr_slice_trans = imread([folder_transformations slice_name '_transformed.tif']);

        % update figure
        update.VerticalScrollCount = 0; set(f, 'UserData', ud);
        ud.loaded = true;

        ud.curr_slice_num = ud.slice_at_shift_start+ud.slice_shift;

        ud.histology_overlay = 1;

        set(ud.text,'Visible','off');
        fill([5 5 250 250],[5 50 50 5],[0 0 0]); ud.text(end+1) = text(5,15,['Slice ' num2str(ud.slice_at_shift_start+ud.slice_shift)],'color','white');
    catch;
        % if no transform, just show reference
        ud.histology_overlay = 0;
        ud.current_pointList_for_transform = zeros(0,2);
        set(ud.im, 'CData', ud.ref);
        ud.curr_im = ud.ref; set(f, 'UserData', ud);
        set(ud.text,'Visible','off');
        fill([5 5 250 250],[5 50 50 5],[0 0 0]); ud.text(end+1) = text(5,15,['Slice ' num2str(ud.slice_at_shift_start+ud.slice_shift) ' - no transform'],'color','white');
    end

end

% update coordinates at the top
pixel = getPixel(ud.atlasAx);
updateStereotaxCoords(ud.currentSlice, pixel, ud.bregma, ud.bregmaText, ud.angleText, ud.currentSlice, ud.currentAngle(1), ud.currentAngle(2), ud.ref_size, ud.plane);

% ----------------------------------------
% if no angle, just change reference slice
% ----------------------------------------
if ud.currentAngle(1) == 0 && ud.currentAngle(2) == 0

    reference_slice = squeeze(allData.tv(ud.currentSlice,:,:));
    ud.im_annotation = squeeze(allData.av(ud.currentSlice,:,:));

    if ud.viewColorAtlas
        set(ud.im, 'CData', ud.im_annotation);
    else
        set(ud.im, 'CData', reference_slice);
    end

   % update title/overlay with brain region
    [name, acr, ann] = getPixelAnnotation(allData, pixel, ud.currentSlice);
    updateTitle(ud.atlasAx, name, acr)
    if ud.showOverlay
        updateOverlay(f, allData, ann, slice_figure, save_location);
    end
    ud.ref = uint8(reference_slice);
    set(ud.pointHands_for_transform(:), 'Visible', 'off');
    ud.offset_map = zeros(ud.ref_size);
    set(f, 'UserData', ud);

% ---------------------------
% if angle, angle the atlas
% ---------------------------
else

  image_size = size(squeeze(allData.av(ud.currentSlice,:,:)));
  angle_slice = zeros(image_size);

  if ud.currentAngle(1)==0; offset_DV = 0;
  else; offset_DV = -ud.currentAngle(1):sign(ud.currentAngle(1)):ud.currentAngle(1);
  end; start_index_DV = 1;

  % loop through AP offsets
  num_DV_iters_add_ind = floor( (image_size(1) - floor( image_size(1) / length(offset_DV))*length(offset_DV)) / 2);
  for curr_DV_iter = 1:length(offset_DV)
      cur_offset_DV = offset_DV(curr_DV_iter);
      if cur_offset_DV == ud.currentAngle(1)
          end_index_DV = image_size(1);
      elseif curr_DV_iter <= num_DV_iters_add_ind  || length(offset_DV - curr_DV_iter) < num_DV_iters_add_ind
          end_index_DV = start_index_DV + floor( image_size(1) / length(offset_DV));
      else
          end_index_DV = start_index_DV + floor( image_size(1) / length(offset_DV)) - 1;
      end

       if ud.currentAngle(2)==0;  offset_ML = 0;
       else; offset_ML = -ud.currentAngle(2):sign(ud.currentAngle(2)):ud.currentAngle(2);
       end; start_index_ML = 1;
    % nested: loop through ML offsets
  num_ML_iters_add_ind = floor( (image_size(2) - floor( image_size(2) / length(offset_ML))*length(offset_ML)) / 2);
  for curr_ML_iter = 1:length(offset_ML)
      cur_offset_ML = offset_ML(curr_ML_iter);
      if cur_offset_ML == ud.currentAngle(2)
          end_index_ML = image_size(2);
      elseif curr_ML_iter <= num_ML_iters_add_ind  || length(offset_ML - curr_ML_iter) < num_ML_iters_add_ind
          end_index_ML = start_index_ML + floor( image_size(2) / length(offset_ML));
      else
          end_index_ML = start_index_ML + floor( image_size(2) / length(offset_ML)) - 1;
      end

      % update current slice
      try
     angle_slice(start_index_DV:end_index_DV, start_index_ML:end_index_ML) = ...
         squeeze(allData.tv(ud.currentSlice + cur_offset_DV + cur_offset_ML,start_index_DV:end_index_DV,start_index_ML:end_index_ML));
      catch
          disp('')
      end

      ud.im_annotation(start_index_DV:end_index_DV,start_index_ML:end_index_ML) = squeeze(allData.av(ud.currentSlice + cur_offset_DV + cur_offset_ML,...
                                                            start_index_DV:end_index_DV,start_index_ML:end_index_ML));
      ud.offset_map(start_index_DV:end_index_DV, start_index_ML:end_index_ML) = cur_offset_DV + cur_offset_ML;

      start_index_ML = end_index_ML + 1;
    end
      start_index_DV = end_index_DV + 1;
  end
  if ud.viewColorAtlas
      set(ud.im, 'CData', ud.im_annotation);
  elseif ~ud.showAtlas
      set(ud.im, 'CData', angle_slice);
  end

  ud.ref = uint8(angle_slice);
  set(ud.pointHands_for_transform(:), 'Visible', 'off');
end

% in all cases. update histology overlay
if ud.histology_overlay == 1 || ud.histology_overlay == 2
    updateHistology(f,ud); ud = get(f, 'UserData');
else
    if ud.viewColorAtlas
        ud.curr_im = ud.im_annotation;
    else
        ud.curr_im = ud.ref;
    end
end

% then update boundary overlay
if ud.showAtlas
    updateBoundaries(f,ud, allData)
end

% in normal mode, show only the probe points for that slice
% in probe view mode, color depends on distance from point
if ud.probe_view_mode && ud.currentProbe
    probe = ud.currentProbe; mean_distance = 0;
    for probe_point = 1:size(ud.pointList{probe,1},1)
        distance = abs(ud.pointList{probe,1}(probe_point,3) - ...
            (ud.currentSlice + ud.offset_map(ud.pointList{probe,1}(probe_point,2),ud.pointList{probe,1}(probe_point,1))));
        mean_distance = mean_distance + distance / size(ud.pointList{probe,1},1);
        if distance < 50
            color = abs((ud.ProbeColors(probe,:) * (50 - distance) + [0 0 0] * distance) / 50);
        else
            color = [.33 .33 .33];
        end
            set(ud.pointHands{probe, 1}(probe_point),'MarkerEdgeColor', color, 'SizeData', 20)
    end
% disp(['mean distance from this slice to probe points is ' num2str(round(mean_distance*10)) ' microns'])
    if mean_distance < 50
        color = abs((ud.ProbeColors(probe,:) * (50 - mean_distance) + [0 0 0] * mean_distance) / 50);
    else
        color = [.33 .33 .33];
    end
    set(ud.pointHands{probe, 3}, 'Color',color)

end
  set(f, 'UserData', ud);

% ---------------------------------------------------------------
% update the image shown if histology is currently being overlaid
% ---------------------------------------------------------------
function updateHistology(f, ud)
    if ud.histology_overlay == 2
        image_blend =  imfuse(uint8(ud.ref*.6), ud.curr_slice_trans(:,:,:),'blend','Scaling','none');
        set(ud.im, 'CData', image_blend);
        ud.curr_im = image_blend;
    elseif ud.histology_overlay == 1
        set(ud.im, 'CData', ud.curr_slice_trans);
        ud.curr_im = ud.curr_slice_trans;
    end
    set(f, 'UserData', ud);

% -------------------------------------------------
% update the position of the region boundary image
% -------------------------------------------------
function updateBoundaries(f, ud, allData)
    if ud.currentAngle(1) == 0 && ud.currentAngle(2) == 0
        curr_annotation = squeeze(allData.av(ud.currentSlice,:,:));
    else
        curr_annotation = ud.im_annotation;
    end

    atlas_vert_1 = double(curr_annotation(1:end-2,:));
    atlas_vert_2 = double(curr_annotation(3:end,:));
    atlas_vert_offset = abs( atlas_vert_1 - atlas_vert_2 ) > 0;
    shifted_atlas_vert1 = zeros(size(curr_annotation(:,:)));
    shifted_atlas_vert1(3:end,:) = atlas_vert_offset;
    shifted_atlas_vert2 = zeros(size(curr_annotation(:,:)));
    shifted_atlas_vert2(1:end-2,:) = atlas_vert_offset;

    atlas_horz_1 = double(curr_annotation(:,1:end-2));
    atlas_horz_2 = double(curr_annotation(:,3:end));
    atlas_horz_offset = abs( atlas_horz_1 - atlas_horz_2 )>0;
    shifted_atlas_horz1 = zeros(size(curr_annotation(:,:)));
    shifted_atlas_horz1(:,3:end) = atlas_horz_offset;
    shifted_atlas_horz2 = zeros(size(curr_annotation(:,:)));
    shifted_atlas_horz2(:,1:end-2) = atlas_horz_offset;

    shifted_atlas = shifted_atlas_horz1 + shifted_atlas_horz2 + shifted_atlas_vert1 + shifted_atlas_vert2;

    atlas_boundaries = (shifted_atlas>0); ud.atlas_boundaries = atlas_boundaries;

    if ud.showAtlas
        image_blend =  uint8( imfuse(ud.curr_im, atlas_boundaries/3.5*(1+.35*isa(ud.curr_im,'uint16')),'blend','Scaling','none') )* 2;
        set(ud.im, 'CData', image_blend);
    end

    set(f, 'UserData', ud);

% ----------------
% react to clicks
% ----------------
function atlasClickCallback(im, keydata, slice_figure, save_location)
f = get(get(im, 'Parent'), 'Parent');
ud = get(f, 'UserData');
ud_slice = get(slice_figure, 'UserData');

% probe view mode
if ud.probe_view_mode && ud.currentProbe
    clickX = round(keydata.IntersectionPoint(1));
    clickY = round(keydata.IntersectionPoint(2));
    if ud.showOverlay
        clickY = size(ud.ref,1) - clickY;
    end
    clickZ = ud.currentSlice + ud.offset_map(clickY,clickX);

    % find the probe point closest to this clicked point
    [min_dist, point_ind] = min( sqrt(sum(([clickX clickY clickZ] - ud.pointList{ud.currentProbe}).^2,2)));

    % find the slice corresponding to that point
    relevant_slice = ud.pointList{ud.currentProbe,2}(point_ind);
    highlight_point = true;
    ud.transformed_slice_figure = transformed_sliceBrowser(ud.transformed_slice_figure, save_location, f, ...
        highlight_point, relevant_slice, min_dist, clickX, clickY, point_ind,0);
    figure(f);

% selecting probe points mode
elseif ud.currentProbe > 0
    clickX = round(keydata.IntersectionPoint(1));
    clickY = round(keydata.IntersectionPoint(2));
    if ud.showOverlay; clickY = size(ud.ref,1) - clickY; end % overlay inverts Y
    clickZ = ud.currentSlice + ud.offset_map(clickY,clickX);

    ud.pointList{ud.currentProbe,1}(end+1, :) = [clickX, clickY, clickZ];
    ud.pointList{ud.currentProbe,2}(end+1, :) = ud.slice_at_shift_start + ud.slice_shift;
    slice_name = ud_slice.processed_image_names{ud.slice_at_shift_start + ud.slice_shift}(1:end-4);
    ud.pointList{ud.currentProbe,3}{end+1} = slice_name;

    ud.pointHands{ud.currentProbe, 1}(end+1) = scatter(ud.atlasAx, clickX, clickY, 20, 'ro', ...
                'MarkerFaceColor', [.1 .1 .1],'MarkerEdgeColor', ud.ProbeColors(ud.currentProbe, :), ...
                'LineWidth',2);

    ud.pointHands{ud.currentProbe, 2}(end+1) = ud_slice.slice_num + ud.slice_shift;

% transform mode
elseif ud.getPoint_for_transform
    clickX = round(keydata.IntersectionPoint(1));
    clickY = round(keydata.IntersectionPoint(2));
    if ud.showOverlay; clickY = size(ud.ref,1) - clickY; end

    if ud.curr_slice_num ~= ud.slice_at_shift_start+ud.slice_shift
        if ~ud.loaded
            ud.current_pointList_for_transform = zeros(0,2);
            disp('transforming new slice');
        end
    end
    ud.pointList_for_transform(end+1, :) = [clickX, clickY];
    ud.current_pointList_for_transform(end+1, :) = [clickX, clickY];
    set(ud.pointHands_for_transform(:), 'color', [.7 .3 .3]);
    ud.pointHands_for_transform(end+1) = plot(ud.atlasAx, clickX, clickY, 'ro', 'color', [0 .9 0],'LineWidth',2,'markers',4);

    ud.slice_at_shift_start = ud_slice.slice_num;
    ud.slice_shift = 0;
    ud.curr_slice_num = ud.slice_at_shift_start+ud.slice_shift;
    ud.loaded = 0;
    ud.clicked = true;
end
set(f, 'UserData', ud);

% ------------------------
% react to mouse hovering
% ------------------------
function fh_wbmfcn(f, allData, slice_figure, save_location)
% WindowButtonMotionFcn for the figure.

ud = get(f, 'UserData');
ax = ud.atlasAx;
pixel = getPixel(ax);

%get offset due to angling
if 0<pixel(1) && pixel(1)<=ud.ref_size(1) && 0<pixel(2) && pixel(2)<=ud.ref_size(2)
    offset = ud.offset_map(pixel(1),pixel(2));
else; offset = 0;
end

% show bregma coords
updateStereotaxCoords(ud.currentSlice + offset, pixel, ud.bregma, ud.bregmaText, ud.angleText, ud.currentSlice, ud.currentAngle(1), ud.currentAngle(2), ud.ref_size, ud.plane);

% get annotation for this pixel
[name, acr, ann] = getPixelAnnotation(allData, pixel, ud.currentSlice+offset);

updateTitle(ax, name, acr);

if ~isempty(name)
    if ud.showContour
        if ~isempty(ud.oldContour)
            delete(ud.oldContour);
        end
        [~,ch] = contour(squeeze(allData.av(ud.currentSlice,:,:)==ann), 1, 'r');
        ud.oldContour = ch;
        set(f, 'UserData', ud);
    end

    if ud.showOverlay
        updateOverlay(f, allData, ann, slice_figure, save_location)
    end
end

% ---------------------------------------------
% update the coordinates shown in the top left
% ---------------------------------------------
function updateStereotaxCoords(currentSlice, pixel, bregma, bregmaText, angleText, slice_num, ap_angle, ml_angle, ref_size, plane)
atlasRes = 0.010; % mm
if strcmp(plane,'coronal')
    ap = -(currentSlice-bregma(1))*atlasRes;
    dv = (pixel(1)-bregma(2))*atlasRes;
    ml = (pixel(2)-bregma(3))*atlasRes;
    set(angleText, 'String', ['Slice ' num2str(bregma(1) - slice_num) ', DV angle ' num2str(round(atand(ap_angle/(ref_size(1)/2)),1)) '^{\circ}, ML angle ' num2str(round(atand(ml_angle/(ref_size(2)/2)),1)) '^{\circ}']);
elseif strcmp(plane,'sagittal')
    ap = -(pixel(2)-bregma(1))*atlasRes;
    dv = (pixel(1)-bregma(2))*atlasRes;
    ml = -(currentSlice-bregma(3))*atlasRes;
    set(angleText, 'String', ['Slice ' num2str(bregma(1) - slice_num) ', DV angle ' num2str(round(atand(ap_angle/(ref_size(1)/2)),1)) '^{\circ}, AP angle ' num2str(round(atand(ml_angle/(ref_size(2)/2)),1)) '^{\circ}']);
elseif strcmp(plane,'transverse')
    ap = -(pixel(2)-bregma(1))*atlasRes;
    dv = (currentSlice-bregma(2))*atlasRes;
    ml = -(pixel(1)-bregma(3))*atlasRes;
    set(angleText, 'String', ['Slice ' num2str(bregma(1) - slice_num) ', ML angle ' num2str(round(atand(ap_angle/(ref_size(1)/2)),1)) '^{\circ}, AP angle ' num2str(round(atand(ml_angle/(ref_size(2)/2)),1)) '^{\circ}']);
end
set(bregmaText, 'String', sprintf('%.2f AP, %.2f DV, %.2f ML', ap, dv, ml));

% ---------------------------------
% update the current mouse location
% ---------------------------------
function pixel = getPixel(ax)

currPoint = get(ax,'currentpoint');  % The current point w.r.t the axis.

Cx = currPoint(1,1); Cy = currPoint(1,2);
pixel = round([Cy Cx]);

% ---------------------------------
% update the overlaid brain region
% ---------------------------------
function updateOverlay(f, allData, ann, slice_figure, save_location)
ud = get(f, 'UserData');
if isempty(ud.overlayAx) % first time
    if ud.currentAngle(1) == 0 && ud.currentAngle(2) == 0
        avo = plotAVoverlay(fliplr(squeeze(allData.av(ud.currentSlice,:,:))'), ann, ud.atlasAx);
    else
        avo = plotAVoverlay(fliplr(ud.im_annotation)', ann, ud.atlasAx);
    end
    ud.overlayAx = avo;
    set(ud.overlayAx, 'HitTest', 'off');
    set(f, 'UserData', ud);
else
    ovIm = get(ud.overlayAx, 'Children');
%     set(ovIm, 'HitTest', 'off');

    if ud.currentAngle(1) == 0 && ud.currentAngle(2) == 0
        thisSlice = squeeze(allData.av(ud.currentSlice,:,:));
    else
        thisSlice = ud.im_annotation;
    end

    set(ovIm, 'CData', flipud(thisSlice));
    plotAVoverlay(fliplr(thisSlice'), ann, ud.atlasAx, ud.overlayAx);
    set(ovIm, 'ButtonDownFcn', @(f,k)atlasClickCallback(f, k, slice_figure, save_location));

end

% ---------------------------------
% find the region being hovered on
% ---------------------------------
function [name, acr, ann] = getPixelAnnotation(allData, pixel, currentSlice)
if pixel(1)>0&&pixel(1)<size(allData.av,2) && pixel(2)>0&&pixel(2)<=size(allData.av,3)
    ann = allData.av(currentSlice,pixel(1),pixel(2));
    name = allData.st.safe_name(ann);
    acr = allData.st.acronym(ann);
else
    ann = []; name = []; acr = [];
end

% ---------------------------------
% update the title, showing region
% ---------------------------------
function updateTitle(ax, name, acr)
if ~isempty(name)
    title(ax, [name{1} ' (' acr{1} ')']);
else
    title(ax, 'not found');
end

author: cortex-lab

repo: allenCCF

save_path: github-repos/MATLAB/cortex-lab-allenCCF


text:
% IN ORDER TO SIMULATE THE PROGRAM:
%   A) FIRST, LOAD A ROBOT
%       robot = load_robot('abb','irb140');
%   B) NEXT, LOAD SOME EQUIPMENT.
%       robot.equipment{1} = load_robot('equipment','tables/table_small');
%       OR
%       robot.equipment{1} = load_robot('equipment','bumper_cutting');
%   C) NOW, LOAD AN END TOOL
%       robot.tool= load_robot('equipment','end_tools/parallel_gripper_0');
%   D) FINALLY, LOAD A PIECE TO GRAB BY THE ROBOT
%       robot.piece=load_robot('equipment','cylinders/cylinder_tiny');
%
%   E) IF NECESSARY, CHANGE THE POSITION AND ORIENTATION OF THE ROBOT'S
%   BASE
%       robot.piece.T0= [1 0 0 -0.35;
%                        0 1 0 -0.55;
%                        0 0 1 0.2;
%                        0 0 0 1];
%
% during the simulation, call simulation_open_tool; to open the tool and
% simulation_close_tool; to close it.
% To grip the piece, call simulation_grip_piece; and
% simulation_release_piece to release it.
% The call to each function must be correct, thus, typically the correct
% sequence is:

% simulation_open_tool;
% approach the piece to grab.
% simulation_close_tool;
% simulation_grip_piece; --> the piece will be drawn with the robot
% move to a different place
% simulation_open_tool;
% simulation_release_piece

function basic_robot_program

global RT_tp1 RT_tp2 RT_tp3 RT_tp4 TD_gripper

TD_gripper=[1,[[0,0,0.125],[1,0,0,0]],[0.1,[0,0,0.100],[1,0,0,0],0,0,0]];
RT_tp1=[[0.541,0.1171,0.713],[0.727812,-0.115363,0.676004,-0.000468],[0,0,-1,1],[9E+09,9E+09,9E+09,9E+09,9E+09,9E+09]];
RT_tp2=[[0.515,-0.200,0.712],[0.7071,0,0.7071,0],[-1,-2,1,0],[9E+09,9E+09,9E+09,9E+09,9E+09,9E+09]];
RT_tp3=[[0.515,0,0.912],[0.7071,0,0.7071,0],[0,-1,0,0],[9E+09,9E+09,9E+09,9E+09,9E+09,9E+09]];
RT_tp4=[[0.515,0,0.512],[0.7071,0,0.7071,0],[0,0,-1,0],[9E+09,9E+09,9E+09,9E+09,9E+09,9E+09]];

main
end

function main

global RT_tp1 RT_tp2 RT_tp3 RT_tp4 TD_gripper

MoveJ(RT_tp1, 'vmax' , 'fine' , TD_gripper, 'wobj0');
MoveJ(RT_tp2, 'vmax' , 'fine' , TD_gripper, 'wobj0');
MoveJ(RT_tp3,'vmax','fine', TD_gripper, 'wobj0');
MoveL(RT_tp4,'vmax','fine', TD_gripper, 'wobj0');
end

repo: ARTE


text:
%kvocoASF 'Alternating Sequential Filter of type Opening(Closing(Opening()))'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros vocoASF.pane file
%
% Parameters:
% InputFile: i 'Input Image ', required: 'input image'
% InputFile: str 'Struct Element', required: 'Structuring Element'
% Integer: n 'Number of open/close/opens', default: 5: 'number of open/close/opens'
% OutputFile: o 'Output Image ', required: 'resulting output image'
%
% Example: o = kvocoASF({i, str}, {'i','';'str','';'n',5;'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% vocoASF - Alternating Sequential Filter of type Opening(Closing(Opening()))
%
%  DESCRIPTION
% .I vocoASF
% will do N successive open/close/opens by the same structuring element.
% Non-sequential filtering can be selected, so that no iteration is done. Otherwise, there will be iteration from 1 to n.
%
% The output image data type will be the same as that of the highest data type
% of the input image.
%
%
%
%  EXAMPLES
% vocoASF -i feath.xv -str b.str -n 5 -o outimage.xv
%
% Will do 5 successive open/close/opens (vocoASF) of the input image "feath.xv" by the structuring element "b.str" and the resulting image will be written to  "outimage.xv".
%
%  "SEE ALSO"
% vcast1(1)
%
%  RESTRICTIONS
% .I vocoASF
% can be defined for all data types supported by Khoros, but at the moment it has been implemented just for the bit and unsigned char types.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993-1997 Junior Barrera, Roberto Lotufo.  All rights reserved.
%

function varargout = kvocoASF(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kvocoASF(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'str', '__input';'n', 5;'o', '__output'};
maxval={0,0,2,0};
minval={0,0,2,0};
istoggle=[0,0,0,0];
was_set=istoggle * 0;
paramtype={'InputFile','InputFile','Integer','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'vocoASF"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function NewFemFile = fem_mergelayers(FemFile, MergeLayers, MergedLabel)
% FEM_MERGELAYERS: Merge multiple layers from a FEM model
%
% USAGE: NewFemFile = fem_mergelayers(FemFile, MergeLayers=[ask], MergedLabel=[ask])

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2020

% ===== PARSE INPUTS =====
if (nargin < 3) || isempty(MergedLabel)
    MergedLabel = [];
end
if (nargin < 2) || isempty(MergeLayers)
    MergeLayers = [];
end
NewFemFile = [];

% ===== INPUTS =====
% Load input file
isProgress = bst_progress('isVisible');
if ~isProgress
    bst_progress('start', 'Merge layers', 'Loading FEM mesh...');
end
FemMat = load(FemFile);
% Get subject
[sSubject, iSubject] = bst_get('SurfaceFile', FemFile);

% Check input layers
if (length(FemMat.TissueLabels) == 1)
    error('There is only one layer in the input file.');
elseif ~isempty(MergeLayers) && any(~ismember(MergeLayers, FemMat.TissueLabels))
    error('The requested layers are not available in the file.');
end
% Ask layers to merge
if isempty(MergeLayers)
    isSelect = ismember(FemMat.TissueLabels, {'white', 'gray', 'csf'});
    isSelect = java_dialog('checkbox', 'Select layers to merge:', 'Merge layers', [], FemMat.TissueLabels, isSelect);
    if isempty(isSelect)
        if ~isProgress
            bst_progress('stop');
        end
        return;
    end
    MergeLayers = FemMat.TissueLabels(isSelect == 1);
end
% Check number of layers selected
if (length(MergeLayers) < 2)
    error('You must select at least two layers.');
end
% Get merged tissue indices
iMerged = cellfun(@(c)find(strcmpi(c, FemMat.TissueLabels)), MergeLayers);

% Ask new layer name
if isempty(MergedLabel)
    MergedLabel = java_dialog('input', 'Please enter a name for the new merged layer:', 'Merge layers', [], 'brain');
    if isempty(MergedLabel)
        if ~isProgress
            bst_progress('stop');
        end
        return;
    end
end
% Check that new layer name doesn't exist
nTissueOld = length(FemMat.TissueLabels);
nTissueNew = nTissueOld - length(iMerged) + 1;
if ismember(MergedLabel, FemMat.TissueLabels(setdiff(nTissueOld, iMerged)))
    error(['Label "' MergedLabel '" already exists in file.']);
end

% ===== RELABEL TISSUES =====
% Update tissue labels
FemMat.TissueLabels{iMerged(1)} = MergedLabel;
FemMat.TissueLabels(iMerged(2:end)) = [];
% Relabel all elements
iRelabel = 1:nTissueNew;
iRelabel(iMerged) = iMerged(1);
iKept = setdiff(1:nTissueOld, iMerged(2:end));
iRelabel(iKept) = 1:length(iKept);
FemMat.Tissue = reshape(iRelabel(FemMat.Tissue), [], 1);

% ===== SAVE NEW FILE =====
if ~isProgress
    bst_progress('start', 'Merge layers', 'Saving new FEM mesh...');
end
% Update comment
FemMat.Comment = strrep(FemMat.Comment, sprintf('%d lay', nTissueOld), sprintf('%d lay', nTissueNew));
% Add history
FemMat = bst_history('add', FemMat, 'merge', ['Merged layers:', sprintf(' %s', MergeLayers{:})]);
% Add file tag
NewFemFile = file_unique(strrep(FemFile, '.mat', '_merge.mat'));
% Save new file
bst_save(NewFemFile, FemMat, 'v7');
% Reference file in database
db_add_surface(iSubject, NewFemFile, FemMat.Comment);

if ~isProgress
    bst_progress('stop');
end

author: brainstorm-tools


text:
function y = cvx_isconstant( x )
y = testall( @cvx_isconstant, x );

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function Dual_TVL1_optic_flow( I0,I1,     u1, u2, nx, ny,     tau,    lambda, theta,   warps,  epsilon,verbose,iflagMedian)
size = round(nx * ny);
l_t = lambda * theta;
if isa(u1,'cell') && isa(u2,'cell')
   u1 = cell2mat(u1);
   u2 = cell2mat(u2);
end
MAX_ITERATIONS= 300;
PRESMOOTHING_SIGMA =0.8;
GRAD_IS_ZERO = 1E-10;
[I1x, I1y] = centered_gradient(I1, 0, 0, nx, ny);
  for  i = 1: size
        p11(i) =0.0;
        p12(i) = 0.0;
        p21(i) =0.0;
        p22(i) = 0.0;
  end
  rho_c = zeros(1,size);
  grad = zeros(1,size);
  for  warpings = 1: warps
      I1w = bicubic_interpolation_warp(I1,  u1, u2, 0,  nx, ny, 1);
      I1wx =  bicubic_interpolation_warp(I1x, u1, u2, 0, nx, ny, 1);
     I1wy = bicubic_interpolation_warp(I1y, u1, u2, 0, nx, ny, 1);
      for  i = 1:size
           Ix2 = I1wx(i) * I1wx(i);
             Iy2 = I1wy(i) * I1wy(i);

            grad(i) = (Ix2 + Iy2);
             sizes = [length(u2),length(I1w),length(I1wx),length(u1),length(I1wy),length(I0)];
if i <= min(sizes)
            rho_c(i) = (I1w(i) - I1wx(i) * u1(i)- I1wy(i) * u2(i) - I0(i));
end
      end
     n = 0;
       error = Inf;
        while (error > epsilon * epsilon && n < MAX_ITERATIONS)
            n = n +1;

            for  i = 1 : size
                 sizes = [length(u2),length(I1w),length(I1wx),length(u1),length(I1wy),length(I0)];
if i <= min(sizes)
               rho = rho_c(i)+ (I1wx(i) * u1(i) + I1wy(i) * u2(i));

               d1 = 0.0;  d2 = 0.0;

                if (rho < - l_t * grad(i))
                    d1 = l_t * I1wx(i);
                    d2 = l_t * I1wy(i);
                 else
                    if (rho > l_t * grad(i))
                        d1 = -l_t * I1wx(i);
                        d2 = -l_t * I1wy(i);
                     else
                        if (grad(i) < GRAD_IS_ZERO)
                            d1 = 0;
                            d2 = 0;
                        else
                             fi = -rho/grad(i);
                            d1 = fi * I1wx(i);
                            d2 = fi * I1wy(i);
                        end

                    end
                end
end
if(i <= length(u1))
                v1(i) = u1(i) + d1;
end
if(i <= length(u2))
                v2(i) = u2(i) + d2;
end
            end
           div_p1 =  divergence(p11, p12, 0, nx ,ny);
           div_p2 = divergence(p21, p22, 0, nx ,ny);
            error = 0.0;

            for  i = 1:  size-1
                if(i <= length(u1))
                u1k = u1(i);
                u2k = u2(i);

                u1(i) = v1(i) + theta * div_p1(i);
                u2(i) = v2(i) + theta * div_p2(i);

                error = error + (u1(i) - u1k) * (u1(i) - u1k) +(u2(i) - u2k) * (u2(i) - u2k);
                end
            end
            error =  error / size;
           [u1x, u1y ] = forward_gradient(u1, 0, 0, nx ,ny);
           [u2x, u2y] =  forward_gradient(u2, 0, 0, nx ,ny);
          for (i = 1:  size-1)
                 taut = tau / theta;
                 if i < length(u1x)
                g1   = hypot(u1x(i), u1y(i));
                g2   = hypot(u2x(i), u2y(i));
                 end
               ng1  = 1.0 + taut * g1;
                 ng2  = 1.0 + taut * g2;

                p11(i) = (p11(i) + taut * u1x(i)) / ng1;
                p12(i) = (p12(i) + taut * u1y(i)) / ng1;
                p21(i) = (p21(i) + taut * u2x(i)) / ng2;
                p22(i) = (p22(i) + taut * u2y(i)) / ng2;
          end

        end

        %applay the filter and continue
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
classdef PrecomputedVariableTest < handle

    properties (Access = protected)
        testName
        variablesToStore
        storedVar
        computedVar
        computation
        computerType
        testResultsName
    end

    methods (Access = public)

        function obj = PrecomputedVariableTest(cParams)
            obj.init(cParams);
            obj.computeVariable(cParams);
            obj.selectComputedVar();
            obj.loadStoredVariable();
        end

        function error = computeError(obj)
            d = numel(obj.variablesToStore);
            err = ones(d,1);
            for ivar = 1:d
                sV = obj.storedVar{ivar};
                cV = obj.computedVar{ivar};
                err(ivar) = norm(sV - cV)/norm(sV);
            end
            error = norm(err);
        end

        function overwriteResults(obj)
            % NOTE: change the variable name "x" to the one in the .mat
            % file so that it works
            variableNames = who('-file', obj.testResultsName);
            d = numel(variableNames);
            for ivar = 1:d
                varName = variableNames{ivar};
                eval([varName ' = obj.computedVar{ivar}' ])
                save(obj.testResultsName, varName, '-append')
            end
        end

        function comp = getComputation(obj)
            comp = obj.computation;
        end
    end

    methods (Access = private)

        function init(obj, cParams)
            obj.testName         = cParams.testName;
            obj.variablesToStore = cParams.variablesToStore;
            obj.computerType     = cParams.computerType;
            obj.testResultsName  = cParams.testName;
            if isfield(cParams, 'testResultsName')
                obj.testResultsName = cParams.testResultsName;
            end
        end

        function computeVariable(obj, s)
            s.testName = obj.testName;
            testComputer = TestComputer.create(obj.computerType, s);
            testComputer.compute();
            obj.computation = testComputer.computation;
        end

        function selectComputedVar(obj)
            toStore = obj.variablesToStore;
            comp = obj.computation;
            if isprop(comp, 'variables')
                vars = obj.computation.variables;
            elseif isprop(comp, 'designVariable')
                vars = struct(obj.computation.designVariable);
                valX = toStore{1};
                [vars.(valX)] = vars.value;
                vars = rmfield(vars, 'value');
            end
            fnms = fieldnames(vars);
            totalComputationVariables = numel(fnms);
            totalStoredVariables = size(toStore,2);
            for i = 1:totalComputationVariables
                calcVar = fnms(i);
                for j = 1:totalStoredVariables
                    storVar = toStore(j);
                    if strcmp(calcVar, storVar)
                        obj.computedVar{j} = vars.(fnms{i});
                        break
                    end
                end
            end
        end

        function loadStoredVariable(obj)
            file2load = obj.testResultsName;
            load(file2load);
            for icell = 1:numel(obj.variablesToStore)
              obj.storedVar{icell} = eval(obj.variablesToStore{icell});
            end
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function om = set_params(om, st, name, idx, params, vals)
%SET_PARAMS  Modifies parameters for variable, cost or constraint in model
%
%   This method can be used to modify parameters for an existing variable,
%   constraint or cost in the model.
%
%   OM.SET_PARAMS(SET_TYPE, NAME, PARAMS, VALS)
%   OM.SET_PARAMS(SET_TYPE, NAME, IDX, PARAMS, VALS)
%
%   Inputs:
%       SET_TYPE : one of 'var', 'lin', 'nle', 'nli', 'nlc', 'qdc' for
%           variables, linear constraints, nonlinear equality constraints,
%           nonlinear inequality constraints, general nonlinear costs,
%           and quadratic costs, respectively
%       NAME : name of set
%       IDX : index of named set (for an indexed set)
%       PARAMS : can be one of three options:
%           1 - 'all', indicating that VALS is a cell array whose elements
%               correspond to the input parameters of the respective
%               add_*() method
%           2 - the name of a PARAM, VAL is the value of that parameter
%           3 - a cell array of PARAM names, VALS is a cell array of
%               corresponding values
%           Note: Changing the dimension of a 'var' is not allowed
%               and changing the #1 ('all') is the only option for 'nle', 'nli', and 'nlc'
%       VALS : new value or cell array of new values for PARAMS
%
%   Valid PARAM names:
%       var - N, v0, vl, vu, vt
%       lin - A, l, u, vs
%       nle - N, fcn, hess, include, vs
%       nli - N, fcn, hess, include, vs
%       nlc - N, fcn, vs
%       qdc - Q, c, k, vs
%
%   Examples:
%       om.set_params('var', 'Pg', 'v0', Pg0);
%       om.set_params('lin', 'y', {2,3}, {'l', 'u'}, {l, u});
%       om.set_params('nle', 'Pmis', 'all', {N, @fcn, @hess, vs});
%
%   See also OPT_MODEL, ADD_VAR, ADD_LIN_CONSTRAINT, ADD_NLN_CONSTRAINT
%            ADD_QUAD_COST and ADD_NLN_COST.

%   MP-Opt-Model
%   Copyright (c) 2008-2020, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MP-Opt-Model.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://github.com/MATPOWER/mp-opt-model for more info.

if nargin < 6
    vals = params;
    params = idx;
    idx = {};
end

%% create default list of parameters to update based on set type & inputs
switch st
    case 'var'
        default_params = {'N', 'v0', 'vl', 'vu', 'vt'};
    case 'lin'
        default_params = {'A', 'l', 'u', 'vs'};
    case {'nle', 'nli'}
        default_params = {'N', 'fcn', 'hess', 'vs'};
    case 'nlc'
        default_params = {'N', 'fcn', 'vs'};
    case 'qdc'
        default_params = {'Q', 'c', 'k', 'vs'};
    otherwise
        error('@opt_model/set_params: ''%s'' is not a valid SET_TYPE', st);
end

%% standardize provided arguments in cell arrays params, vals
is_all = 0;     %% flag to indicate all params for set are being replaced
if ischar(params)
    if strcmp(params, 'all')
        is_all = 1;
        np = length(vals);      %% number of parameter values provided
        params = default_params(1:np);
    else
        np = 1;                 %% number of parameter values provided
        params = {params};
        vals = {vals};
    end
else
    np = length(vals);          %% number of parameter values provided
end

if ~isempty(idx)
    %% calls to substruct() are relatively expensive, so we pre-build the
    %% struct for addressing cell array fields
    %% sc = substruct('.', name, '{}', idx);
    sc = struct('type', {'.', '{}'}, 'subs', {name, idx});  %% cell array field
end

switch st
    case 'var'
        %% get current parameters
        [v0, vl, vu, vt] = om.params_var(name, idx);
        N0 = om.getN('var', name, idx);
        p = struct('N', N0, 'v0', v0, 'vl', vl, 'vu', vu, 'vt', vt);    %% current parameters
        u = struct('N',  0, 'v0',  0, 'vl',  0, 'vu',  0, 'vt',  0);    %% which ones to update

        %% replace with new parameters
        for k = 1:np
            p.(params{k}) = vals{k};
            u.(params{k}) = 1;
        end
        N = p.N;

        %% set missing default params for 'all'
        if is_all
            if np < 5
                p.vt = 'C';
                u.vt = 1;               %% update vt
                if np < 4
                    p.vu = Inf(N, 1);
                    u.vu = 1;           %% update vu
                    if np < 3
                        p.vl = -Inf(N, 1);
                        u.vl = 1;       %% update vl
                        if np < 2
                            p.v0 = zeros(N, 1);
                            u.v0 = 1;   %% update v0
                        end
                    end
                end
            end
        end

        %% check consistency of parameters
        %% no dimension change
        if N ~= N0
            error('@opt_model/set_params: dimension change for ''%s'' ''%s'' not allowed', st, nameidxstr(name, idx));
        end

        %% check sizes of new values of v0, vl, vu, vt
        for pn = {'v0', 'vl', 'vu', 'vt'}
            if u.(pn{1})
                nn = length(p.(pn{1}));
                if nn ~= N
                    if nn == 0
                        switch pn{1}
                            case 'v0'
                                p.(pn{1}) = zeros(N, 0);
                            case 'vl'
                                p.(pn{1}) = -Inf(N, 0);
                            case 'vu'
                                p.(pn{1}) =  Inf(N, 0);
                            case 'vt'
                                p.(pn{1}) = 'C';
                        end
                    elseif nn == 1
                        if pn{1} ~= 'vt'
                            p.(pn{1}) = p.(pn{1}) * ones(N, 1);   %% expand from scalar
                        end
                    else
                        error('@opt_model/set_params: parameter ''%s'' ''%s'' ''%s'' should have length %d (or 1)', st, nameidxstr(name, idx), pn{1}, N);
                    end
                end
            end
        end

        %% assign new parameters
        if isempty(idx)     %% simple named set
            for k = 2:length(default_params)
                pn = default_params{k};     %% param name
                if u.(pn)   %% assign new val for this parameter
                    om.var.data.(pn).(name) = p.(pn);
                end
            end
        else                %% indexed named set
            for k = 2:length(default_params)
                pn = default_params{k};     %% param name
                if u.(pn)   %% assign new val for this parameter
                    om.var.data.(pn) = subsasgn(om.var.data.(pn), sc, p.(pn));
                end
            end
        end
    case 'lin'
        %% get current parameters
        [A, l, u, vs] = om.params_lin_constraint(name, idx);
        [N0, M0] = size(A);
        if isempty(vs), vs = {vs}; end
        p = struct('A', A, 'l', l, 'u', u, 'vs', vs);   %% current parameters
        u = struct('A', 0, 'l', 0, 'u', 0, 'vs',  0);   %% which ones to update

        %% replace with new parameters
        for k = 1:np
            p.(params{k}) = vals{k};
            u.(params{k}) = 1;
        end

        %% set missing default params for 'all'
        [N, M] = size(p.A);
        if is_all
            u.A = 1;            %% always update A
            u.l = 1;            %% alwaus update l
            if np < 4
                p.vs = {};
                u.vs = 1;       %% update vs
                if np < 3
                    p.u = Inf(N, 1);
                    u.u = 1;    %% update u
                end
            end
        end

        %% check consistency of parameters
        %% no dimension change unless 'all'
        if N ~= N0 && ~is_all
            error('@opt_model/set_params: dimension change for ''%s'' ''%s'' not allowed except for ''all''', st, nameidxstr(name, idx));
        end

        %% check sizes of new values of l, u
        for pn = {'l', 'u'}
            if u.(pn{1})
                nn = length(p.(pn{1}));
                if nn ~= N
                    if nn == 0
                        switch pn{1}
                            case 'l'
                                p.(pn{1}) = -Inf(N, 0);
                            case 'u'
                                p.(pn{1}) =  Inf(N, 0);
                        end
                    elseif nn == 1
                        p.(pn{1}) = p.(pn{1}) * ones(N, 1);   %% expand from scalar
                    else
                        error('@opt_model/set_params: parameter ''%s'' ''%s'' ''%s'' should have length %d (or 1)', st, nameidxstr(name, idx), pn{1}, N);
                    end
                end
            end
        end

        %% check consistency of A and vs
        if u.A || u.vs
            p.vs = om.varsets_cell2struct(p.vs);
            nv = om.varsets_len(p.vs);      %% number of variables
            if M ~= nv
                error('@opt_model/set_params: for ''%s'' ''%s'' number of columns of ''A'' (%d) must be consistent with ''vs'' (%d)', st, nameidxstr(name, idx), M, nv);
            end
        end

        %% assign new parameters
        if isempty(idx)     %% simple named set
            for k = 1:length(default_params)
                pn = default_params{k};     %% param name
                if u.(pn)   %% assign new val for this parameter
                    om.lin.data.(pn).(name) = p.(pn);
                end
            end
        else                %% indexed named set
            for k = 1:length(default_params)
                pn = default_params{k};     %% param name
                if u.(pn)   %% assign new val for this parameter
                    om.lin.data.(pn) = subsasgn(om.lin.data.(pn), sc, p.(pn));
                end
            end
        end
    case 'qdc'
        %% get current parameters
        [Q, c, kk, vs] = om.params_quad_cost(name, idx);
        [MQ0, NQ0] = size(Q);
        [Mc0, Nc0] = size(c);
        Nk0 = length(kk);
        nx0 = max([MQ0 Mc0 Nk0]);
        N0 = om.getN('qdc', name, idx);
        if isempty(vs), vs = {vs}; end
        p = struct('Q', Q, 'c', c, 'k', kk, 'vs', vs);  %% current parameters
        u = struct('Q', 0, 'c', 0, 'k',  0, 'vs',  0);  %% which ones to update

        %% replace with new parameters
        for k = 1:np
            p.(params{k}) = vals{k};
            u.(params{k}) = 1;
        end

        %% set missing default params for 'all'
        [MQ, NQ] = size(p.Q);
        [Mc, Nc] = size(p.c);
        Nk = length(p.k);
        nx = max([MQ Mc Nk]);
        if NQ == 1 || (isempty(p.Q) && (Nk > 1 || k == 0))
            %% Q is a column vector (cost is element-wise, i.e. a vector)
            %% OR Q is empty and k is either a vector or zero
            N = nx;
        else            %% Q is a square matrix (cost is a scalar)
            N = 1;
        end
        if is_all
            u.Q = 1;                %% always update Q
            if np < 4
                p.vs = {};
                u.vs = 1;           %% update vs
                if np < 3
                    p.k = 0;
                    u.k = 1;        %% update k
                    if np < 2
                        p.c = [];
                        u.c = 1;    %% update c
                    end
                end
            end
        end

        %% check consistency of parameters
        %% no dimension change unless 'all'
        if (N ~= N0 || nx ~= nx0) && ~is_all
            error('@opt_model/set_params: dimension change for ''%s'' ''%s'' not allowed except for ''all''', st, nameidxstr(name, idx));
        end

        %% Q and c can't both be empty
        if ~MQ && ~Mc
            error('@opt_model/set_params: ''%s'' ''%s'' : ''Q'' and ''c'' cannot both be empty', st, nameidxstr(name, idx));
        end

        %% check sizes of new values of Q, c, k
        if ~isempty(p.Q) && (MQ ~= nx || (MQ ~= NQ && NQ ~= 1) )
            error('@opt_model/set_params: ''%s'' ''%s'' : ''%s'' is expected to be (%d x %d)', st, nameidxstr(name, idx), 'Q', MQ, NQ);
        end
        if ~isempty(p.c) && Mc ~= nx
            error('@opt_model/set_params: ''%s'' ''%s'' : ''%s'' is expected to be (%d x %d)', st, nameidxstr(name, idx), 'c', Mc, 1);
        end
        if ~isempty(p.k) && any(p.k) && Nk ~= N && Nk ~= 1
            error('@opt_model/set_params: ''%s'' ''%s'' : ''%s'' is expected to be (%d x %d)', st, nameidxstr(name, idx), 'k', N, 1);
        end

        %% check consistency of Q, c, k and vs
        if u.Q || u.c || u.vs
            p.vs = om.varsets_cell2struct(p.vs);
            nv = om.varsets_len(p.vs);      %% number of variables
            if nx ~= nv
                error('@opt_model/set_params: for ''%s'' ''%s'' dimensions of ''Q'', ''c'', ''k'' (%d) must be consistent with ''vs'' (%d)', st, nameidxstr(name, idx), nx, nv);
            end
        end

        %% assign new parameters
        if isempty(idx)     %% simple named set
            for k = 1:length(default_params)
                pn = default_params{k};     %% param name
                if u.(pn)   %% assign new val for this parameter
                    om.qdc.data.(pn).(name) = p.(pn);
                end
            end
        else                %% indexed named set
            for k = 1:length(default_params)
                pn = default_params{k};     %% param name
                if u.(pn)   %% assign new val for this parameter
                    om.qdc.data.(pn) = subsasgn(om.qdc.data.(pn), sc, p.(pn));
                end
            end
        end
    case {'nle', 'nli'}
        %% get current parameters
        if isempty(idx)
            [N0, fcn, hess, vs, include] = om.params_nln_constraint(st(3) == 'e', name, idx);
        else
            [N0, fcn, hess, vs] = om.params_nln_constraint(st(3) == 'e', name, idx);
            include = '';
        end
        if isempty(vs), vs = {vs}; end
        p = struct('N', N0, 'fcn', fcn, 'hess', hess, 'vs', vs);    %% current parameters
        u = struct('N',  0, 'fcn',   0, 'hess',    0, 'vs',  0);    %% which ones to update

        %% replace with new parameters
        for k = 1:np
            p.(params{k}) = vals{k};
            u.(params{k}) = 1;
        end
        N = p.N;

        %% set missing default params for 'all'
        if is_all
            u.N    = 1;         %% always update N
            u.fcn  = 1;         %% alwaus update fcn
            u.hess = 1;         %% alwaus update hess
            if np < 4
                p.vs = {};
                u.vs = 1;       %% update vs
            end
        end

        %% check consistency of parameters
        %% no dimension change unless 'all'
        if N ~= N0 && ~is_all
            error('@opt_model/set_params: dimension change for ''%s'' ''%s'' not allowed except for ''all''', st, nameidxstr(name, idx));
        end

        %% included constraints not yet implemented
        if ~isempty(include)
            error('@opt_model/set_params: modifications for ''%s'' ''%s'' not (yet) supported since it includes evaluation of other constraints', st, nameidxstr(name, idx));
        end

        %% convert vs to struct
        if u.vs
            p.vs = om.varsets_cell2struct(p.vs);
        end

        %% assign new parameters
        if isempty(idx)     %% simple named set
            for k = 2:length(default_params)
                pn = default_params{k};     %% param name
                if u.(pn)   %% assign new val for this parameter
                    om.(st).data.(pn).(name) = p.(pn);
                end
            end
        else                %% indexed named set
            for k = 2:length(default_params)
                pn = default_params{k};     %% param name
                if u.(pn)   %% assign new val for this parameter
                    om.(st).data.(pn) = subsasgn(om.(st).data.(pn), sc, p.(pn));
                end
            end
        end
    case 'nlc'
        %% get current parameters
        [N0, fcn, vs] = om.params_nln_cost(name, idx);
        if isempty(vs), vs = {vs}; end
        p = struct('N', N0, 'fcn', fcn, 'vs', vs);  %% current parameters
        u = struct('N',  0, 'fcn',   0, 'vs',  0);  %% which ones to update

        %% replace with new parameters
        for k = 1:np
            p.(params{k}) = vals{k};
            u.(params{k}) = 1;
        end
        N = p.N;

        %% set missing default params for 'all'
        if is_all
            u.N   = 1;          %% always update N
            u.fcn = 1;          %% alwaus update fcn
            if np < 3
                p.vs = {};
                u.vs = 1;       %% update vs
            end
        end

        %% check consistency of parameters
        %% no dimension change unless 'all'
        if N ~= N0 && ~is_all
            error('@opt_model/set_params: dimension change for ''%s'' ''%s'' not allowed except for ''all''', st, nameidxstr(name, idx));
        end

        %% vector valued costs not yet implemented
        if N ~= 1
            error('@opt_model/set_params: vector value for ''%s'' ''%s'' not yet implemented', st, nameidxstr(name, idx));
        end

        %% assign new parameters
        if isempty(idx)     %% simple named set
            for k = 2:length(default_params)
                pn = default_params{k};     %% param name
                if u.(pn)   %% assign new val for this parameter
                    om.nlc.data.(pn).(name) = p.(pn);
                end
            end
        else                %% indexed named set
            for k = 2:length(default_params)
                pn = default_params{k};     %% param name
                if u.(pn)   %% assign new val for this parameter
                    om.nlc.data.(pn) = subsasgn(om.nlc.data.(pn), sc, p.(pn));
                end
            end
        end
    otherwise
        error('@opt_model/set_params: ''%s'' is not a valid SET_TYPE', st);
end

%% update dimensions and indexing, if necessary
dN = N - N0;
if is_all && dN
    %% calls to substruct() are relatively expensive, so we pre-build the
    %% struct for addressing num array fields
    %% sn = substruct('.', name, '()', idx);
    sn = struct('type', {'.', '()'}, 'subs', {'', 1});  %% num array field
    om_ff = om.(st);
    update = 0;             %% not yet reached set being updated
    update_i1 = 0;          %% flag to indicate whether to update i1
    for k = 1:om_ff.NS
        o = om_ff.order(k);
        if ~update && strcmp(o.name, name) && isequal(o.idx, idx)
            update = 1;     %% arrived at set being updated
        end
        if update
            if isempty(o.idx)   %% simple named set
                if update_i1
                    om_ff.idx.i1.(o.name) = om_ff.idx.i1.(o.name) + dN;
                else
                    om_ff.idx.N.(o.name) = om_ff.idx.N.(o.name) + dN;
                end
                om_ff.idx.iN.(o.name) = om_ff.idx.iN.(o.name) + dN;
            else                %% indexed named set
                sn(1).subs = o.name;
                sn(2).subs = o.idx;
                if update_i1
                    v = subsref(om_ff.idx.i1, sn);
                    om_ff.idx.i1 = subsasgn(om_ff.idx.i1, sn, v + dN);
                else
                    v = subsref(om_ff.idx.N, sn);
                    om_ff.idx.N = subsasgn(om_ff.idx.N, sn, v + dN);
                end
                v = subsref(om_ff.idx.iN, sn);
                om_ff.idx.iN = subsasgn(om_ff.idx.iN, sn, v + dN);
            end
            update_i1 = 1;  %% update i1 from here on out
        end
    end
    om_ff.N = om_ff.N + dN;
    om.(st) = om_ff;
end

function str = nameidxstr(name, idx)
str = sprintf('%s%s', name, idxstr(idx));

function str = idxstr(idx)
if isempty(idx)
    str = '';
elseif length(idx) == 1
    str = sprintf('(%d)', idx{1});
else
    str = ['(' sprintf('%d', idx{1}) sprintf(',%d', idx{2:end}) ')'];
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function [V3D] = Vandermonde3D(N, r, s, t);

% function [V3D] = Vandermonde3D(N, r, s, t);
% Purpose : Initialize the 3D Vandermonde Matrix, V_{ij} = phi_j(r_i, s_i, t_i);

V3D = zeros(length(r),(N+1)*(N+2)*(N+3)/6);

% Transfer to (a,b) coordinates
[a, b, c] = rsttoabc(r, s, t);

% build the Vandermonde matrix
sk = 1;

for i=0:N % old ordering
  for j=0:N - i
    for k=0:N - i - j
      V3D(:,sk) = Simplex3DP(a,b,c,i,j,k);
      sk = sk+1;
    end
  end
end
return;

repo: Aero-matlab


text:
function [ SimilarityMatrix ] = similarityNeighbor( x, n, ~)
%similarityNeighbor Create a link for the n'th neighbour

    sz = size(x,1);
    SimilarityMatrix = eye(sz);

    i = 1:sz-n;
    SimilarityMatrix(sub2ind([sz, sz], i+n,i)) = 1;
    SimilarityMatrix(sub2ind([sz, sz], i,i+n)) = 1;

    DiagMask = ones(size(x, 1)) - eye(size(x,1));
    SimilarityMatrix = SimilarityMatrix .* DiagMask;
    SimilarityMatrix = SimilarityMatrix + eye(size(x, 1));

end

author: TadasBaltrusaitis

repo: OpenFace


text:
mex -lopencv_core -lopencv_imgproc -L/usr/local/lib -I/usr/local/include/ resize.cpp MxArray.cpp
mex -lopencv_core -lopencv_imgproc -L/usr/local/lib -I/usr/local/include/ calcIIF.cpp MxArray.cpp
mex im2colstep.c

author: flyers

repo: drone-tracking


text:
node = [0,0,0; 1,0,0; 0,1,0; 0,0,1];
elem = [1 2 3 4];


text:
function [H,G] = make_bank_3(h,nbands)

% make_bank_3      Creates Analysis/Synthesis Filters for Pseudo-QMF Filter Banks
%
%                  Standard type of cosine modulation where the phase reference is (L-1)/2
% Arguments:
% h                Impulse response of prototype filter
% ndands           Number of subbands
%
% by Lee, Gan, and Kuo, 2008
% Subband Adaptive Filtering: Theory and Implementation
% Publisher: John Wiley and Sons, Ltd

flen = max(size(h));
H = zeros(nbands,flen);

% Note: k starts from 1, thus (k-0.5) and (k-1) in the cosine modulations

n = 0:flen-1;
for k=1:nbands
    H(k,:) = h.*cos((pi/nbands)*(k-0.5)*(n-(flen-1)/2)+ (pi/2)*(k+0.5)*(flen-1-nbands)/nbands)*2;
end

% Synthesis filters are time-reversed versions of the analysis filters

% Gk(z) = z^(-L+1)H(z^-1)
G = fliplr(H);

author: CharlesThaCat

repo: acoustic-interference-cancellation


text:
function r = plus(p,q)
% this function adds the x, y, and z values of the structures that I use in
% the FVtool.
%
% SYNOPSIS:
%
%
% PARAMETERS:
%
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

% Copyright (c) 2012-2016 Ali Akbar Eftekhari
% See the license file

if (isa(p, 'FaceVariable')&&isa(q, 'FaceVariable'))
    r=p;
    r.xvalue = p.xvalue+q.xvalue;
    r.yvalue = p.yvalue+q.yvalue;
    r.zvalue = p.zvalue+q.zvalue;
elseif isa(p, 'FaceVariable')
    r=p;
    r.xvalue = p.xvalue+q;
    r.yvalue = p.yvalue+q;
    r.zvalue = p.zvalue+q;
else
    r=q;
    r.xvalue = p+q.xvalue;
    r.yvalue = p+q.yvalue;
    r.zvalue = p+q.zvalue;
end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function i4mat_transpose_print_some ( m, n, a, ilo, jlo, ihi, jhi, title )

%*****************************************************************************80
%
%% I4MAT_TRANSPOSE_PRINT_SOME prints some of an I4MAT, transposed.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 June 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns.
%
%    Input, integer A(M,N), an M by N matrix to be printed.
%
%    Input, integer ILO, JLO, the first row and column to print.
%
%    Input, integer IHI, JHI, the last row and column to print.
%
%    Input, string TITLE, an optional title.
%
  incx = 10;

  if ( 0 < s_len_trim ( title ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, '%s\n', title );
  end

  for i2lo = max ( ilo, 1 ) : incx : min ( ihi, m )

    i2hi = i2lo + incx - 1;
    i2hi = min ( i2hi, m );
    i2hi = min ( i2hi, ihi );

    inc = i2hi + 1 - i2lo;

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Row: ' );
    for i = i2lo : i2hi
      fprintf ( 1, '%7d  ', i );
    end
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Col\n' );
    fprintf ( 1, '\n' );

    j2lo = max ( jlo, 1 );
    j2hi = min ( jhi, n );

    for j = j2lo : j2hi

      fprintf ( 1, '%5d  ', j );
      for i2 = 1 : inc
        i = i2lo - 1 + i2;
        fprintf ( 1, '%7d  ', a(i,j) );
      end
      fprintf ( 1, '\n' );

    end

  end

  fprintf ( 1, '\n' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function result___ = get_function_symbol(expression___)
% internal: some function_handle expressions have a function symbol (an @name expression), and this function obtains it

% note: we are using funny names here to bypass potential name conflicts within the eval() clause further below
if ~isa(expression___,'function_handle')
    error('the expression has no associated function symbol.'); end

string___ = char(expression___);
if string___(1) == '@'
    % we are dealing with a lambda function
    if is_symbolic_lambda(expression___)
        result___ = eval(string___(27:end-21));
    else
        error('cannot derive a function symbol from a non-symbolic lambda function.');
    end
else
    % we are dealing with a regular function handle
    result___ = expression___;
end

repo: REST


text:
function u = pend_control(x)

load('pend_data', 'K', 'U', 'domain');

p = [x(2) x(4)];
u = tpcontroller(p, x, K, U, domain);

author: Sable

repo: mcbench-benchmarks