text:
function f = allenAtlasBrowser(f, templateVolume, annotationVolume, structureTree, 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
%

if nargin<4
    save_location = '';
    save_suffix = '';
end

fprintf(1, 'Controls: \n');
fprintf(1, '--------- \n');
fprintf(1, 'scroll: move between slices \n');
fprintf(1, 'g: add/remove gridlines \n');
fprintf(1, 'o: add/remove overlay of current region extent \n');
fprintf(1, 'a: toggle viewing boundaries \n');
fprintf(1, 'v: toggle color atlas mode \n');
fprintf(1, 'p: enable/disable mode where clicks are logged for probe or switch probes \n');
fprintf(1, 'n: trace a new probe \n');
fprintf(1, 'l: load saved probe points \n');
fprintf(1, 's: save current probe \n');
fprintf(1, 'w: enable/disable probe viewer mode for current probe  \n');
fprintf(1, 'd: delete most recent probe point \n');

fprintf(1, 'up: scroll through A/P angles \n');
fprintf(1, 'right: scroll through M/L angles \n');
fprintf(1, 'down: scroll through slices \n');

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,'Name','Atlas Viewer','Position', [1000*screen_size(1) 250*screen_size(2) 1100*screen_size(1) 850*screen_size(2)])
movegui(f,'onscreen')

ud.bregma = allenCCFbregma;

ud.plane = plane;
ud.currentSlice = ud.bregma(1);
ud.currentAngle = zeros(2,1);
ud.scrollMode = 0;
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.showAtlas = false;
ud.viewColorAtlas = false;
ud.histology_overlay = 0;
ud.atlasAx = axes('Position', [0.05 0.05 0.9 0.9]);
ud.angleText = annotation('textbox', [.7 0.95 0.4 0.05], ...
    'EdgeColor', 'none', 'Color', 'k');

reference_image = squeeze(templateVolume(ud.currentSlice,:,:));
ud.ref_size = size(reference_image);
ud.im = plotTVslice(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 = 10;

set(ud.im, 'ButtonDownFcn', @(f,k)atlasClickCallback(f, k));

ud.bregmaText = annotation('textbox', [0 0.95 0.4 0.05], ...
    'String', '[coords]', 'EdgeColor', 'none', 'Color', 'k');

allData.tv = templateVolume;
allData.av = annotationVolume;
allData.st = structureTree;

hold(ud.atlasAx, 'on');
set(ud.atlasAx, 'HitTest', 'on');

set(f, 'UserData', ud);

set(f, 'KeyPressFcn', @(f,k)hotkeyFcn(f, k, allData, save_location, save_suffix));
set(f, 'WindowScrollWheelFcn', @(src,evt)updateSlice(f, evt, allData))
set(f, 'WindowButtonMotionFcn',@(f,k)fh_wbmfcn(f, allData)); % Set the motion detector.

function hotkeyFcn(f, keydata, allData, save_location, save_suffix)

ud = get(f, 'UserData');

switch lower(keydata.Key)
    case 'g' % toggle showing Gridlines
        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
    case 'o' % toggle Overlay
        ud.showOverlay = ~ud.showOverlay;
        if ~ud.showOverlay
            delete(ud.overlayAx); ud.overlayAx = [];
            disp('Overlay OFF');
        elseif ~ud.viewColorAtlas; disp('Overlay on!');
        end

    case 'p' % toggle mode to register clicks as Probe points
        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} ')']);
        else
            ud.currentProbe = 0;
            disp(['probe point mode OFF']);
        end
    case 'l' % load saved probe points
         % 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;
                disp([num2str(size(ud.pointList,1)) ' probes loaded'])
        else
            disp('probe points already present')
        end
    case 'w'
        % remove overlay
        ud.showOverlay = 0;
        delete(ud.overlayAx); ud.overlayAx = [];

        % plot probe points
        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;
                end; disp(['activate probe view mode for probe ' num2str(ud.currentProbe)])

                for probe_point = 1:size(ud.pointHands{ud.currentProbe, 1}(:),1)
                    ud.pointHands{ud.currentProbe, 1}(probe_point) = scatter(ud.atlasAx, ...
                            ud.pointList{ud.currentProbe,1}(probe_point,1), ud.pointList{ud.currentProbe,1}(probe_point,2), 20, 'ro',...
                            'MarkerFaceColor', [ .1 .1 .1], 'MarkerEdgeColor', ud.ProbeColors(ud.currentProbe, :), 'MarkerFaceAlpha',.4,'LineWidth',2);
                end

                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);

                % 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); ud = get(f, 'UserData');

            else; disp('probe view mode OFF'); end

    case 'a' % toggle View boundaries

        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

    case 'v' % toggle View Color Atlas
        ud.viewColorAtlas = ~ud.viewColorAtlas;

        if ud.viewColorAtlas
            % turn off overlay first
            ud.showOverlay = 0;
            ref_mode = false;
            delete(ud.overlayAx); ud.overlayAx = [];
            set(ud.im, 'CData', ud.im_annotation)
            cmap = allen_ccf_colormap('2017');
            colormap(ud.atlasAx, cmap); caxis(ud.atlasAx, [1 size(cmap,1)]);
        else
            set(ud.im, 'CData', ud.ref)
            colormap(ud.atlasAx, 'gray'); caxis(ud.atlasAx, [0 400]);
        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 'n' % new probe
        new_num_probes = size(ud.pointList,1) + 1;
        if new_num_probes <= size(ud.ProbeColors,1)
            disp(['probe ' num2str(new_num_probes) ' added (' ud.ProbeColor{new_num_probes} ')']);
            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
    case 's' % save probe trajectory and points of each probe per histology image (and associated histology name/number)
        pointList.pointList = ud.pointList;
        pointList.pointHands = ud.pointHands;
        save(fullfile(save_location ,['probe_points' save_suffix]), 'pointList'); disp('probe points saved');

    case 'd' % delete current transform or most recent probe point
        if 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
set(f, 'UserData', ud);

function updateSlice(f, evt, allData)

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 A/P 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;
end

% update title/overlay with brain region
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 do normal thing
% ---------------------------------
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

    [name, acr, ann] = getPixelAnnotation(allData, pixel, ud.currentSlice);
    updateTitle(ud.atlasAx, name, acr)
    if ud.showOverlay
        updateOverlay(f, allData, ann);
    end
    ud.ref = uint8(reference_slice);
    set(f, 'UserData', ud);
% ---------------------------
% if angle, do angle thing
% ---------------------------
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);

end

if ud.viewColorAtlas
    ud.curr_im = ud.im_annotation;
else
    ud.curr_im = ud.ref;
end

% then update boundary overlay unless in color view mode
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), 'MarkerFaceColor',color,'MarkerEdgeColor', color, 'SizeData', 20)
    end
% disp(['mean distance from slice to 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 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/5,'blend','Scaling','none') )* 2;
        set(ud.im, 'CData', image_blend);
    end

    set(f, 'UserData', ud);

function atlasClickCallback(im, keydata)
f = get(get(im, 'Parent'), 'Parent');
ud = get(f, 'UserData');

if 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 Yn
    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, :) = 0;
    slice_name = 'no histology image';
    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, :), 'MarkerFaceAlpha',.4,'LineWidth',2);
    ud.pointHands{ud.currentProbe, 2}(end+1) = 0;

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

function fh_wbmfcn(f, allData)
% 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)
    end
end

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));

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]);

function updateOverlay(f, allData, ann)
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));

end

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

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:
function [sp, K, scores] = update_sp_indexing(sp, K, h, w, opts)
% After spagglom_sub.m has been run for a while (phase1, obtaining the
% refined superpixelation), The 'sp' variable contains lots of useless
% empty superpixels and the superpixels are made of multiple parts. If this
% information is not required (by default it's not), this function can be
% run to simplify the 'sp' variable and accordingly update neighbor
% variable 'K'.

cc = 1; % a counter
sp_old = sp;
sp = [];
key = [];

% Build each superpixel
for q = 1:length(sp_old)
    if size(sp_old{q}.pixels,1) > 0 % remove empty sp by ignoring them
        sp{cc} = sp_old{q}; % copy the old information
        sp{cc}.parts = [cc]; % each sp has one part, which is the same as its index
        key(end+1,1:2) = [q, cc]; % contruct mapping from old to new indexing
        cc = cc + 1;
    end
end

sp_amount = length(sp);

% Spply mapping to update sp labels
K = arrayfun(@(x)key(key(:,1) == x, 2), K); % changes each value of K by lookup using table matrix 'key'
for q = 1:length(sp)
    sp{q}.neighbors = arrayfun(@(x)key(key(:,1) == x, 2), sp{q}.neighbors); % same for neighbors
end

% List pixel coordinates in index format
for sus = 1:sp_amount
    sp{sus}.spind = sub2ind([h,w], double(sp{sus}.pixels(:,1)), double(sp{sus}.pixels(:,2)));
end

% Recalculate scores
scores = similarity_scores(sp, K, opts); % rows of this variable have been permutated. Instead of solving the permutation, we just recalculate the values.

author: Cloud-CV

repo: object-proposals


text:
function tc = tc_visualizeGlm(tc, parent);
%
% tc = tc_visualizeGlm(tc, parent);
%
% Provide graphics that visualize the
% setup and result of applying a general
% linear model to a time course. Requires
% that you first run tc_applyGlm on
% the tc time course struct.
%
%
% ras 04/05.
if nargin<1,    tc = get(gcf,'UserData');     end
if nargin<2,    parent = tc.ui.plot;          end
if parent==gcf | parent==get(gcf, 'CurrentAxes')
    % make a uipanel to fit on the target
    parent = uipanel('Parent', parent, ...
        'Units', 'normalized', ...
        'BackgroundColor', get(gcf, 'Color'), ...
        'Position', [0 0 1 1]);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% apply a GLM if needed
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ~isfield(tc, 'glm')
    tc = tc_applyGlm(tc);
end
X = tc.glm.designMatrix;
Y = tc.wholeTc(:);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% clean up existing objects in figure
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
otherAxes = findobj('Type', 'axes','Parent', parent);
delete(otherAxes);
otherUiControls = findobj('Type', 'uicontrol', 'Parent',parent);
delete(otherUiControls);
axes('Parent', parent); % shift focus to parent uipanel
delete(gca);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Visualize results: diff't options for deconvolved
% data and non-deconvolved (using HRF)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isequal(tc.glm.type, 'selective averaging')
    tc = tc_plotDeconvolvedTCs(tc, parent);
    return
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% if we got here, the GLM is non-deconvolved, using HRF   %
% show hemodynamic response function used for GLM         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
axes('Parent', parent, 'Position', [.05 .6 .15 .2]); % [.13 .58 .21 .34]
plot(tc.glm.hrf,'k','LineWidth',2);
if tc.params.grid==1
    grid on
end
xlabel('Time, frames')
ylabel('Arbitrary Units')
if isnumeric(tc.params.glmHRF)
    opts = {sprintf('Mean trial for conditions \n %s',num2str(tc.params.snrConds)), ...
        'Boynton gamma function', 'SPM difference-of-gammas' ...
        'Dale & Buckner ''97'};
    hrfName = opts{tc.params.glmHRF};
else
    hrfName = tc.params.glmHRF; hrfName(hrfName=='_') = ' ';
end
title({'HRF function used: ' hrfName}, 'FontWeight', 'bold')
axis tight, axis square, set(gca, 'Box', 'off');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% show design matrix                     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
axes('Parent', parent, 'Position', [.3 .55 .2 .34]); % [.13 .11 .21 .34]
hImg = imagesc(X);
colormap autumn
nConds = sum(tc.trials.condNums>0);
tickPts = [1:nConds];
tickLabels = {'Individual Conditions' 'DC Predictors for each run'};
set(gca,'XTick',tickPts);
xlabel('Predictors (Conditions + DC)')
ylabel('Time, Frames')
title('Design Matrix', 'FontWeight', 'bold')
set(hImg, 'ButtonDownFcn', 'zoom');
colorbar vert

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% bar beta values for selected conditions %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
axes('Parent', parent, 'Position', [.6 .55 .33 .34]);  %[.57 .58 .33 .34]
sel = setdiff( find(tc_selectedConds(tc)), 1);
tmp = tc.trials.condColors(sel);
for i = 1:length(sel), col(i,:) = tmp{i}; end
xstr = tc_condInitials(tc.trials.condNames(sel));
% for i = 1:length(sel), xstr{i} = num2str(tc.trials.condNums(sel(i))); end
mybar(tc.glm.betas(sel-1), tc.glm.sems(sel-1), xstr, [], col);
% xlabel Predictors
% ylabel('% Signal Change')
set(gca, 'Box', 'off');
if tc.params.grid, grid on; else grid off; end
title('Beta Values', 'FontWeight', 'bold');
ylabel('\beta', 'FontWeight', 'bold', 'Rotation', 0, 'FontSize', 14);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% show time course + selected predictors %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if tc.params.legend, ysz = .58; else ysz = .72; end
tc.glm.tcAxes = axes('Parent', parent, 'Position', [.1 .1 ysz .3]);

% construct matrix of predictor functions for each condition:
% (these predictors will include the dc components as well, so we
% won't have them as separate traces):
dc = nConds+1:size(tc.glm.betas, 2);  % indices of DC components
for c = 1:nConds
	predictors(:,c) = X(:,[c dc]) * [tc.glm.betas(:,[c dc])'];
end

% plot
t = [1:length(tc.wholeTc)] .* tc.TR;
hold on
hY = plot(t, tc.wholeTc(:), 'k-', 'LineWidth', 2);
hPred = plot(t, predictors(:,sel-1), 'LineStyle', '-', 'LineWidth', 1.5);
hTot = plot(t, X*tc.glm.betas', 'k-', 'LineWidth', 2);
hRes = plot(t, tc.glm.residual(:), 'k--', 'LineWidth', 1.5);
plot(t, zeros(size(t)), 'k:');
set(hPred, 'Visible', 'off');
set(hRes, 'Visible', 'off');
set(gca, 'Box', 'off');
if tc.params.grid, grid on; else grid off; end
setLineColors([{'k'} tc.trials.condColors(sel) {[.2 .4 .8] [.8 0 0]}]);
xlabel('Time, sec')
ylabel('% Signal')
title('Time course + Scaled Predictors', 'FontWeight', 'bold')
zoom

% for longer time courses, may want a UI control
% scrollbar(tc.glm.tcAxes, tc.wholeTc);
scale = 300;  % max seconds to nicely plot TC
scrollbar(gca, scale);

% create toggles for the different traces
if tc.params.legend, xx = .72; else xx = .9; end
cb = ['tmp = {''off'' ''on''}; val = get(gcbo, ''Value'')+1; ' ...
      'set(get(gcbo,''UserData''), ''Visible'', tmp{val}); ' ...
      'clear tmp val '];

h3 = uicontrol('Style', 'checkbox', 'Units', 'normalized', ...
               'Position', [xx .4 .08 .04], 'String', 'Time Course', ...
               'BackgroundColor', 'w', 'UserData', hY, ...
               'Value', 1, 'Callback', cb);
h4 = uicontrol('Style', 'checkbox', 'Units', 'normalized', ...
               'Position', [xx .32 .08 .04], 'String', 'Predictors', ...
               'BackgroundColor', 'w', 'UserData', hPred, ...
               'Value', 0, 'Callback', cb);
h5 = uicontrol('Style', 'checkbox', 'Units', 'normalized', ...
               'Position', [xx .24 .08 .04], 'String', 'Best Fit', ...
               'BackgroundColor', 'w', 'UserData', hTot, ...
               'Value', 1, 'Callback', cb);
h6 = uicontrol('Style', 'checkbox', 'Units', 'normalized', ...
               'Position', [xx .16 .08 .04], 'String', 'Residual', ...
               'BackgroundColor', 'w', 'UserData', hRes, ...
               'Value', 0, 'Callback', cb);

% lastly, show % variance explained
varEx = sprintf('%2.1f%% Variance Explained', tc.glm.varianceExplained * 100);
uicontrol('Style', 'text', 'Units', 'normalized', ...
          'Position', [xx-.1 .08 .18 .04], 'String', varEx, ...
          'BackgroundColor', 'w','FontSize', 15);

return
% /-------------------------------------------------------------------/ %

% /-------------------------------------------------------------------/ %
function tc = tc_plotDeconvolvedTCs(tc, parent)
% plot the deconvolved time courses in the same way as for the
% default plot, combining mean amplitudes on one side and mean time
% courses on the other.
sel = find(tc_selectedConds(tc));
% sel = sel(sel>1); % no baseline estimated
nConds = length(sel);
frameWindow = unique(round(tc.timeWindow./tc.TR));
prestim = -1 * frameWindow(1);
peakFrames = unique(round(tc.peakPeriod./tc.TR));
bslFrames = unique(round(tc.bslPeriod./tc.TR));
peakFrames = find(ismember(frameWindow,peakFrames));
bslFrames = find(ismember(frameWindow,bslFrames));

%%%%%%%%%%%%%%%%%%%%%%%%
% plot mean amplitudes %
%%%%%%%%%%%%%%%%%%%%%%%%
axes('Position', [.1 .2 .35 .6], 'Parent', parent);
% exp7_plotAmps(tc);
hold on

h1 = gca;

lineWidth = 2;
labels = {};

colors = tc.trials.condColors(sel);
X = tc.trials.condNums(sel);
Y = tc.glm.amps(sel-1);
E = tc.glm.amp_sems(sel-1);
mybar(Y, E, tc_condInitials(tc.trials.condNames(sel)), [], colors);
axis tight
set(gca, 'Box', 'off')
if tc.params.grid==1
    grid on
end

% set line width
htmp = findobj('Type','line','Parent',gca);
set(htmp,'LineWidth',lineWidth);

% add labels
set(gca,'XTick',[1:nConds]);
if tc.params.legend==0, set(gca, 'XTickLabel', labels); end
xlabel('Condition', 'FontWeight', 'bold', 'FontAngle', 'italic');
ylabel('Mean Amplitude, % Signal', 'FontWeight', 'bold', ...
    'FontAngle', 'italic');
title('Deconvolved Amplitudes', 'FontWeight', 'bold');

% set axes to frame bars nicely
% axis auto;
AX = axis;
AX(1:2) = [0 nConds+1];
if isfield(tc.params,'axisBounds') & ~isempty(tc.params.axisBounds)
    AX(3:4) = tc.params.axisBounds(3:4);
end
axis(AX);

%%%%%%%%%%%%%%%%%%%%%
% mean time courses %
%%%%%%%%%%%%%%%%%%%%%
h2 = axes('Position', [.5 .2 .4 .6], 'Parent', parent);
hold on

for i = sel-1
    htmp = errorbar(tc.timeWindow, tc.glm.betas(:,i), tc.glm.sems(:,i));
    set(htmp, 'Color', tc.trials.condColors{i+1}, 'LineWidth', 2);
end

% indicate the peak and baseline periods, if selected
if tc.params.showPkBsl==1
    AX = axis;
    plot(tc.bslPeriod, repmat(AX(3),size(tc.bslPeriod)), ...
        'k', 'LineWidth', 3.5);
    plot(tc.peakPeriod, repmat(AX(4),size(tc.peakPeriod)), ...
        'r', 'LineWidth', 3.5);
end

if tc.params.grid==1
    grid on
end

if isfield(tc.params,'axisBounds') & ~isempty(tc.params.axisBounds)
    axis(h2, tc.params.axisBounds);
end

xlabel('Trial time, secs', 'FontWeight', 'bold', 'FontAngle', 'italic');
ylabel('% Signal', 'FontWeight', 'bold', 'FontAngle', 'italic');
title('Deconvolved Time Courses', 'FontWeight', 'bold');

AX = axis;
txt = sprintf('Variance Explained: %3.1f%%', 100*tc.glm.varianceExplained);
text(AX(1) + .1*diff(AX(1:2)), AX(3) + .9*diff(AX(3:4)), txt, ...
           'FontSize', 12, 'FontWeight', 'bold');

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function curSliceOri = getCurSliceOri(view)
%
% curSliceOri = getCurSliceOri(view)
%
% Gets current slice orientation from the sliceOriButtons handles
try,
  curSliceOri = findSelectedButton(view.ui.sliceButtons);
catch,
  curSliceOri = findSelectedButton(view.ui.sliceOriButtons);
end


text:
function tri2tgf(input_name,output_name)
  % TRI2TGF Script to convert a triangle mesh to a tgf graph (skeleton file)
  %
  % tri2tgf(input_name,output_name)
  %
  % Inputs:
  %  input_name  name of file contaning input vertices and faces of surface,
  %    should end in .off or .obj
  %  output_name  name of output file to be written with 3D tet mesh, should
  %    end in .mesh

  % read input mesh
  [V,F] = load_mesh(input_name);
  % Find all edges in mesh, note internal edges are repeated
  E = sort( ...
    [F(:,1) (:,2); ...
     F(:,2) (:,3); ...
     F(:,3) (:,1)]')';
  % remove duplicates
  E = unique(E,'rows');
  writeTGF(output_name,V,E);

end

author: alecjacobson

repo: gptoolbox


text:
function [Im]=Walsh_DWT_Decoding(Header)
%Walsh-Harmard Transform with Tow level Discrete Wavelete Transform for image Decompression
% Designed by  Mohammed M. Siddeq
% Date 2012-2-22
% Email :- mamadmmx76@yahoo.com
%
% this program is used for Decompress grayscale images by using :
% INPUT\  Header : this parameter contains all infmration about compressed file
%
%OUTPUT\  Im : Decoded image from "Header"

% --- Very important for you ----------------------------------------------
%Befroe apply this function you need the following programs:
% 1- "Arith_code.m" and "Arith_Decode.m"
% 2- "Walsh2D_Transform.m"
%Note\ these programs free for download from the following website:
%     www.MathWorks.com  -> Auther = Siddeq

%--------------------------------------------------------------------------
% Example (1) --- for compression grayscale images-------------------------
%    I = imread('D:\images\image1.bmp'); % read an image
%    [Header]=Walsh_DWT_Coding(I,[0.025, 0.025],'db3',2); % Apply function
%
%
%    [Im]=Walsh_DWT_Decoding(Header); % for decoding;
%    imshow(uint8(Im)); % show final decoded image

% Example (2) -------------------------------------------------------------
%    I = imread('D:\images\imageG.bmp'); %-> read an image
%    [Header]=Walsh_DWT_Coding(I,[0.025, 0.025],'db5',3); %-> Apply function
%    save('C:\imageG.WWT', 'Header'); -> save compressed file....
%
%   X=load('C:\imageG.WWT', '-mat'); %-> load compressed file....
%   H=X.Header;
%   Im=Walsh_DWT_Decoding(H); %-> for decoding;
%   imshow(Im); %-> show final decoded image
%
%----------------------------------- Good Luck -------------------------

%--Read Header file for Decoding.....
S_2=double(Header(1).S_2);
S_=double(Header(2).S_);
Code_Arr=double(Header(3).Code_Arr);
Arr_H1=double(Header(4).Arr_H1);
Arr_H2=double(Header(5).Arr_H2);
Code_Data1=double(Header(6).Code_Data1);
Data1_H1=double(Header(7).Data1_H1);
Data1_H2=double(Header(8).Data1_H2);
Code_Data2=double(Header(9).Code_Data2);
Data2_H1=double(Header(10).Data2_H1);
Data2_H2=double(Header(11).Data2_H2);
Code_Data3=double(Header(12).Code_Data3);
Data3_H1=double(Header(13).Data3_H1);
Data3_H2=double(Header(14).Data3_H2);
DC_Values=double(Header(15).DC_Values);
M=double(Header(16).M);
Para=double(Header(17).Para);
F1=double(Header(18).F1);
F2=double(Header(19).F2);
wave_name=Header(20).wave_name;
%---------------------------------------------------------------

Arr_Size=8; %window size 2x2 or 4x4 or 8x8 used by Walsh-Transformed ................

%--------Compute the Quantization Matrix
Q(1:Arr_Size,1:Arr_Size)=1; % initilize quantization matrix........
L=M.*F1;
for i=1:Arr_Size
    for j=1:Arr_Size
        Q(i,j)=round((L)+(i+j).*Para);
    end;
end;

%----------Using Arihtmetic Decoding for extract original data.....
[Arr]=Arith_Decode(Code_Arr,Arr_H1,Arr_H2); % Low-Freqeuncy Decoding...

% High-Frequency Decoding......................................
if (Code_Data1==0)
  Data1(1:S_(1)*S_(2))=0;
else
  [Data1]=Arith_Decode(Code_Data1,Data1_H1,Data1_H2);
end;

if (Code_Data2==0)
 Data2(1:S_(1)*S_(2))=0;
else
 [Data2]=Arith_Decode(Code_Data2,Data2_H1,Data2_H2);
end;

if (Code_Data3==0)
  Data3(1:S_(1)*S_(2))=0;
else
  [Data3]=Arith_Decode(Code_Data3,Data3_H1,Data3_H2);
end;
%-------------------------------------------------------------------------

%Retrun "SAVE_" matrix from one-dimensional array "Arr"
POS=1;
SAVE_(1:S_2(1),1:S_2(2))=0;
sizeofArr=size(Arr);
for j=1:S_2(2)
  if (POS<=sizeofArr(2)) % to be ensure the postions is not exceeding
    for i=1:S_2(1)
      if (POS<=sizeofArr(2)) % to be ensure again not exceeding
        SAVE_(i,j)=Arr(POS);
      end;
      POS=POS+1;
    end;
  end;
end;

% Return "DC_Values" to matrix "SAVE_".....
SAVE_(1:S_2(1),1)=DC_Values(1:S_2(1));

% Deocde - HH2
% Return all original data from array
k=1;

for j=1:S_(2)
    for i=1:S_(1)
        HH2(i,j)=Data3(k);
        k=k+1;
    end;
end;
%-----------------------------------------------------------------------
% Decode HL2
% Return all original data from array
k=1;
for j=1:S_(2)
    for i=1:S_(1)
        HL2(i,j)=Data2(k);
        k=k+1;
    end;
end;
%-----------------------------------------------------------------------
% Decode LH2
% Return all original data from array
k=1;
for j=1:S_(2)
    for i=1:S_(1)
        LH2(i,j)=Data1(k);
        k=k+1;
    end;
end;

%------------------- Reconstruct "LL2" from matrix "SAVE_"
L=1;i=1;
while (i<=S_(1))
    j=1;
    while(j<=S_(2))
        Z=1;
        for k1=1:Arr_Size
            for k2=1:Arr_Size
                T(k1,k2)=SAVE_(L,Z);
                Z=Z+1;
            end;
        end;
        L=L+1;
        %------Apply inverse Quantization and inverse Walsh Transform------
        T=T.*Q;
        X=Walsh2D_Transform(T,Arr_Size,'i');
        D(i:i+Arr_Size-1,j:j+Arr_Size-1)=X(1:Arr_Size,1:Arr_Size);
        j=j+Arr_Size;
    end;
    i=i+Arr_Size;
end;
%-------------------------------------------------------------------------

%Adjustment for D matrix to be same size alike "LH2" or others
RELL(1:S_(1),1:S_(2))=D(1:S_(1),1:S_(2));
%RELL=RELL.*1
% Apply inverse Quantization for each sub-band (just high-Freqeuncy)
LH2=round(LH2.*(M*F2)); HL2=round(HL2.*(M.*F2)); HH2=round(HH2.*(M.*F2));

% Inverse -Wavelet Transform Stage-1
S_=size(LH2); LL2=0;
LL2(1:S_(1),1:S_(2))=RELL(1:S_(1),1:S_(2));
L1=idwt2(LL2,LH2,HL2,HH2,wave_name);
L1=L1*Para;
% Inverse -Wavelet Transform Stage-2
S_=size(L1);
HL(1:S_(1),1:S_(2))=0; % all main sub-bands are zeros

%L=0;
%L(1:S_(1),1:S_(2))=L1(1:S_(1),1:S_(2));

Im=idwt2(L1,HL,HL,HL,wave_name);% Decoded image
Im=uint8(Im);
'Decompression ..OK'
end

author: Sable

repo: mcbench-benchmarks


text:
function obj = G_func_pTop_new(in1,in2,in3,in4)
coef_f0_q_sym2 = in1(:,2);
coef_f0_q_sym3 = in1(:,3);
coef_f0_q_sym4 = in1(:,4);
coef_f0_q_sym11 = in1(:,11);
coef_f1_q_sym1 = in2(:,1);
coef_f1_q_sym5 = in2(:,5);
coef_f1_q_sym6 = in2(:,6);
coef_f1_q_sym7 = in2(:,7);
coef_f0_q_sym18 = in1(:,18);
coef_f1_q_sym8 = in2(:,8);
coef_f1_q_sym9 = in2(:,9);
coef_f2_q_sym1 = in3(:,1);
coef_f0_q_sym22 = in1(:,22);
coef_f0_q_sym24 = in1(:,24);
coef_f2_q_sym5 = in3(:,5);
coef_f2_q_sym6 = in3(:,6);
coef_f2_q_sym7 = in3(:,7);
coef_f2_q_sym8 = in3(:,8);
coef_f2_q_sym9 = in3(:,9);
coef_f3_q_sym1 = in4(:,1);
coef_f3_q_sym5 = in4(:,5);
coef_f3_q_sym6 = in4(:,6);
coef_f3_q_sym7 = in4(:,7);
coef_f3_q_sym8 = in4(:,8);
coef_f3_q_sym9 = in4(:,9);
coef_f1_q_sym10 = in2(:,10);
coef_f1_q_sym11 = in2(:,11);
coef_f1_q_sym18 = in2(:,18);
coef_f1_q_sym22 = in2(:,22);
coef_f1_q_sym24 = in2(:,24);
coef_f2_q_sym10 = in3(:,10);
coef_f2_q_sym11 = in3(:,11);
coef_f2_q_sym18 = in3(:,18);
coef_f2_q_sym22 = in3(:,22);
coef_f2_q_sym24 = in3(:,24);
coef_f3_q_sym10 = in4(:,10);
coef_f3_q_sym11 = in4(:,11);
coef_f3_q_sym18 = in4(:,18);
coef_f3_q_sym22 = in4(:,22);
coef_f3_q_sym24 = in4(:,24);
t2 = coef_f1_q_sym1.*2.0;
t3 = coef_f2_q_sym1.*2.0;
t4 = coef_f3_q_sym1.*2.0;
obj = reshape([coef_f0_q_sym11-coef_f1_q_sym18,-coef_f2_q_sym18,-coef_f3_q_sym18,-coef_f1_q_sym22,coef_f0_q_sym11-coef_f2_q_sym22,-coef_f3_q_sym22,-coef_f1_q_sym24,-coef_f2_q_sym24,coef_f0_q_sym11-coef_f3_q_sym24,coef_f0_q_sym2-coef_f1_q_sym5+coef_f0_q_sym18+coef_f1_q_sym11+t2,-coef_f2_q_sym5+coef_f2_q_sym11+t3,-coef_f3_q_sym5+coef_f3_q_sym11+t4,coef_f0_q_sym3-coef_f1_q_sym6+coef_f0_q_sym22,coef_f0_q_sym2+coef_f0_q_sym18-coef_f2_q_sym6,-coef_f3_q_sym6,coef_f0_q_sym4-coef_f1_q_sym7+coef_f0_q_sym24,-coef_f2_q_sym7,coef_f0_q_sym2+coef_f0_q_sym18-coef_f3_q_sym7,-coef_f1_q_sym8+coef_f1_q_sym11+t2,coef_f0_q_sym3+coef_f0_q_sym22-coef_f2_q_sym8+coef_f2_q_sym11+t3,-coef_f3_q_sym8+coef_f3_q_sym11+t4,-coef_f1_q_sym9,coef_f0_q_sym4+coef_f0_q_sym24-coef_f2_q_sym9,coef_f0_q_sym3+coef_f0_q_sym22-coef_f3_q_sym9,-coef_f1_q_sym10+coef_f1_q_sym11+t2,-coef_f2_q_sym10+coef_f2_q_sym11+t3,coef_f0_q_sym4+coef_f0_q_sym24-coef_f3_q_sym10+coef_f3_q_sym11+t4],[3,9]);

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function varargout = zettl(X)
% Zettl function
%
%   ZETTL([x1, x2]) returns the value of the Zettle function at the
%   specified points. [x1] and [x2] may be vectors. The search domain
%   is
%
%               -5 < x_i < 5
%
%   The global minimum is
%
%               f(x1, x2) = f(-0.0299, 0) = -0.003791

% Author: Rody P.S. Oldenhuis
% Delft University of Technology
% E-mail: oldenhuis@dds.nl
% Last edited 20/Jul/2009

    % if no input is given, return dimensions, bounds and minimum
    if (nargin == 0)
        varargout{1} = 2;        % # dims
        varargout{2} = [-5, -5]; % LB
        varargout{3} = [+5, +5]; % UB
        varargout{4} = [-2.989597760285287e-002, 0]; % solution
        varargout{5} = -3.791237220468656e-003; % function value at solution

    % otherwise, output function value
    else
        % keep all values in the search domain
        X(X < -5) = inf;        X(X > 5) = inf;

        % split input vector X into x1, x2
        if size(X, 1) == 2
            x1 = X(1, :);        x2 = X(2, :);
        else
            x1 = X(:, 1);        x2 = X(:, 2);
        end

        % output function value
        varargout{1} = (x1.^2 + x2.^2 - 2*x1).^2 + x1/4;
    end

end


text:
function s = logsumexp(x, dim)
% Returns log(sum(exp(x),dim)) while avoiding numerical underflow.
% Default is dim = 1 (columns).
% Written by Mo Chen (mochen@ie.cuhk.edu.hk). March 2009.

% Downloaded from:
% http://www.mathworks.com/matlabcentral/fileexchange/28899-logsumexp/content/logsumexp.m
%
% Copyright (c) 2010, Michael Chen
% All rights reserved.

% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:

%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution

% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.

if nargin == 1,
    % Determine which dimension sum will use
    dim = find(size(x)~=1,1);
    if isempty(dim), dim = 1; end
end

% subtract the largest in each column
y = max(x,[],dim);
x = bsxfun(@minus,x,y);
s = y + log(sum(exp(x),dim));
i = find(~isfinite(y));
if ~isempty(i)
    s(i) = y(i);
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [trg] = read_eep_trg(fn);

% READ_EEP_TRG reads triggers from an EEProbe *.trg file
%
% This function returns an Nx1 array with the N triggers
%
% [trg] = read_eep_trg(filename)
%
% trg(i).time   ... trigger latency in ms
% trg(i).offset ... byte offset
% trg(i).code   ... trigger code (string)
% trg(i).type   ... numeric value of trg.code
%
% where i as number between 1 and N (the number of triggers found)
%
% An EEProbe trigger file is formatted like
%   0.00195312 256
%   0.000     10350  __
%   17.033   2242926   1
%   20.535   2701934   5
%   21.096   2775406  13
%   21.098   2775662   8
%   ...
% where the first column is the trigger latency in seconds, the second
% column is the byte offset in the file and the third column is the triggercode.
% The triggercode can be numeric or a string. The first line of the file contains the
% sample duration.
%
% Author: Robert Oostenveld, Aalborg University, Denmark, 11 March 2003
%
% See also READ_EEP_CNT, READ_EEP_REJ, READ_EEP_AVR
%

% Copyright (C) 2002, Robert Oostenveld
%                     Aalborg University, Denmark
%                     http://www.smi.auc.dk/~roberto/
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

% $Log: not supported by cvs2svn $
% Revision 1.1  2004/11/26 13:17:02  jwiskerke
% Added m-files without binary code in maple distribution.
%
% Revision 1.2  2003/10/24 13:34:41  Maarten-Jan Hoeve
% Added GNU Licence and updated revision history
%
% Revision 1.1.1.2  2003/10/17 09:55:20  mvelde
% updated: consistent copyrights, arguments/data labels, fixed some typos
%
% Revision 1.1.1.1  2003/03/11 15:24:51  roberto
% updated help and copyrights
% ANT Software BV, The Netherlands, www.ant-neuro.com / info@ant-neuro.com
%

trg = [];

fid = fopen(fn, 'rb');
if fid<0
   return
end

header = fgetl(fid);
while ~feof(fid)
  tmp = fscanf(fid, '%f %d %s', 3);
  if ~isempty(tmp)
    new.time   = tmp(1)*1000;			% in ms
    new.offset = tmp(2)+1;			% offset 1
    new.code   = char(tmp(3:end));		% string
    new.type   = str2double(new.code);		% numeric
    trg = [trg; new];
  end
end

fclose(fid);

author: spm


text:
function [caNodeIndices, vResolution_] = ex_CreateIndexCatalog(mCatalog, mPolygon, bMap, nGriddingMode, nNumberEvents, fRadius, fSizeRectHorizontal, fSizeRectDepth)
    % Creates a cell-array with subcatalogs for every grid node defined by mPolygon
    %
    % [caNodeIndices] = ex_CreateIndexCatalog(mCatalog, mPolygon, bMap, nGriddingMode,
    %                                                  nNumberEvents, fRadius, fSizeRectHorizontal, fSizeRectDepth)
    %
    % Creates a cell-array with subcatalogs for every grid node defined by mPolygon. These subcatalogs
    %   contain only indices to the earthquake "rows" in mCatalog.
    %
    % Input parameters:
    %   mCatalog              Earthquake catalog
    %   mPolygon              Polygon (defined by ex_selectgrid)
    %   bMap                  Calculate cell-array for a map (true) or a cross-section (false)
    %   nGriddingMode         Mode of creating grid node subcatalogs
    %                         0: Constant number of events
    %                         1: Constant radius
    %                         2: Rectangular grid node samples
    %   nNumberEvents         Number of events per grid node (nGriddingMode == 0)
    %   fRadius               Radius of grid node sample (nGriddingMode == 1)
    %   fSizeRectHorizontal   Latitude/horizontal size of rectangle (nGriddingMode == 2)
    %   fSizeRectDepth        Longitude/depth size of rectangle (nGriddingMode == 2)
    %
    % Output parameters:
    %   caNodeIndices         Cell-array with index-catalogs per grid node of mPolygon
    %  vResolution_             Cell-array with distance to grid node of mPolygon
    %
    % Danijel Schorlemmer
    % June 17, 2002
    %
    % updates
    % Thomas van Stiphout August 14, 2006
    %  returns the vector vResolution that contains
    % either the maximal radius or the number of earthquake per gridnode
    %

    report_this_filefun();

    % Create the catalogs for each node with pointers to the overall catalog
    nNumberNodes_ = length(mPolygon(:,1));
    caNodeIndices = cell(nNumberNodes_, 1);
    vResolution_ = cell(nNumberNodes_, 1);
    % If cross-section calculate the length along cross-section
    if ~bMap
        nRow_ = mCatalog.Count;
        vXSecX_ = mCatalog(:,nColumn_);  % length along x-section
        vXSecY_ = (-1) * mCatalog.Depth;  % depth of hypocenters
    end
    % Loop over all points of the polygon
    for nNode_ = 1:nNumberNodes_
        % Get the grid node coordinates
        fX_ = mPolygon(nNode_, 1);
        fY_ = mPolygon(nNode_, 2);
        if (nGriddingMode == 0) | (nGriddingMode == 1)  % Fixed radius or fixed number
            % Calculate distance from center point
            if bMap
                %             vDistances_=deg2km(distance(mCatalog.Latitude,mCatalog.Longitude,ones(size(mCatalog,1),1)*fY_,ones(size(mCatalog,1),1)*fX_));
                vDistances_ = sqrt(   ((mCatalog.Longitude-fX_)*cosd(fY_)*111).^2   +     ((mCatalog.Latitude-fY_)*111).^2)  ;
            else
                vDistances_ = sqrt(((vXSecX_ - fX_)).^2 + ((vXSecY_ - fY_)).^2);
            end
            if nGriddingMode == 0 % Fixed number
                if mCatalog.Count == 0
                    caNodeIndices{nNode_} = [];
                elseif nNumberEvents > mCatalog.Count
                    caNodeIndices{nNode_} = vIndices(1:mCatalog.Count);
                    %                 Determine the Resolution i.e., radius
                    vResolution_{nNode_}=vDistances_(caNodeIndices{nNode_});
                else
                    % Use first nNumberEvents events
                    [vTmp, vIndices] = sort(vDistances_);
                    caNodeIndices{nNode_} = vIndices(1:nNumberEvents);
                    %                 Determine the Resolution i.e., radius
                    vResolution_{nNode_}=vDistances_(caNodeIndices{nNode_});
                end
            else % Fixed radius
                % Use all events within fRadius
                caNodeIndices{nNode_} = find(vDistances_ <= fRadius);
                %             Determine the Resolution i.e., radius
                vResolution_{nNode_}=size(caNodeIndices{nNode_},1);
            end
        else % Rectangular gridding (nGriddingMode == 2)
            if bMap
                vSel_ = ((mCatalog.Longitude >= (fX_ - fSizeRectHorizontal/2)) & (mCatalog.Longitude < (fX_ + fSizeRectHorizontal/2)) & ...
                    (mCatalog.Latitude >= (fY_ - fSizeRectDepth/2)) & (mCatalog.Latitude < (fY_ + fSizeRectDepth/2)));
            else
                vSel_ = ((vXSecX_ >= (fX_ - fSizeRectHorizontal/2)) & (vXSecX_ < (fX_ + fSizeRectHorizontal/2)) & ...
                    (vXSecY_ >= (fY_ - fSizeRectDepth/2)) & (vXSecY_ < (fY_ + fSizeRectDepth/2)));
            end
            caNodeIndices{nNode_} = find(vSel_ == 1);
            vResolution_{nNode_}=size(caNodeIndices{nNode_},1);
        end
    end % of for nNode_
end

author: CelsoReyes


text:
%% THE OPENMAS CYCLE PROCESSOR (OMAS_process.m) %%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function computes the progressive cycles of the simulation, for the
% environment, obstacles, agents and objects withing the simulation space.

% Author: James A. Douthwaite 06/10/2016

% //////////////////////// MAIN WRAPPER SCRIPT ////////////////////////////
function [DATA,SIM,EVENTS,objectIndex]  = OMAS_process(META,objectIndex)
% INPUTS:
% objectIndex - The cell array of object classes
% OUTPUTS:
% DATA        - The output DATA structure with timeseries data
% SIM         - The terminal META structure copy

%% PRE-SIMULATION DECLARATIONS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
META.phase = 'META'; % Declare the simuation active phase
EVENTS = [];         % Reset the EVENT log container

% PREPARE THE OUTPUT DATA CONTAINER
[DATA] = GetOutputStructure(META);

%%%%%%%%%%%%%%%%%%%%%% BEGIN TIME-STEP INTERATIONS %%%%%%%%%%%%%%%%%%%%%%%%
fprintf('[%s]\n[%s]\tLAUNCHING SIMULATION...\n[%s]\n',META.phase,META.phase,META.phase);
% EXECUTE OMAS MAIN CYCLE
[META,objectIndex,DATA,EVENTS] = ProcessGlobalTimeSeries(META,objectIndex,DATA,EVENTS);
% DATA.computationTime = toc; % Get total elapsed simulation time
% fprintf('\n[%s]\tSIMULATION COMPLETE (time elapsed: %ss)\n',META.phase,num2str(DATA.computationTime));
fprintf('[%s]\n[%s]\t...SIMULATION COMPLETE\n[%s]\n',META.phase,META.phase,META.phase);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% FINALISE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
META.phase = 'OUTPUT';                                   % Change to OUTPUT phase
% clearvars -except META EVENTS DATA objectIndex         % Clear loose simulation variables
fprintf('[%s]\tDumping initial simulation variables to output directory\n[%s]\n',META.phase,META.phase);
% sendOutputToFiles(META,EVENTS,DATA,objectIndex);
fprintf('[%s]\tClosing simulation.\n',META.phase);
SIM = META;
end

% //////////////////////// MAIN CYCLE FUNCTIONS ///////////////////////////
% PROCESS TIME VECTOR
function [META,objectIndex,DATA,EVENTS] = ProcessGlobalTimeSeries(META,objectIndex,DATA,EVENTS)
% This function computes the simulation loop cycle across the given
% descrete META.timevector:
% INPUTS:
% META - The initialised META data-structure
% objectIndex - The cell vector of objects (and/or agents)
% DATA        - The initialised output DATA structure
% EVENTS      - The initialised EVENTS structure list

% PROCESS TIME CYCLE
step = META.TIME.currentStep;
exitFlag = 0;
while step <= META.TIME.numSteps

    % 0. /////////////////// CHECK AGENT IDLE STATUS //////////////////////
    % If there are agents in the simulation, then break if all are idle,
    % else allow run to full duration
    idleCondition = sum([META.OBJECTS([META.OBJECTS.type] == OMAS_objectType.agent).idleStatus]) == META.totalAgents && META.totalAgents > 0;
    if idleCondition == 1 && ~exitFlag
        exitFlag = 1;
        referenceTime = META.TIME.currentTime;
        fprintf('[%s]\t... All agents are idle, aborting in t=%0.2fs.\n',META.phase,META.TIME.idleTimeOut);
    elseif exitFlag == 1
        timeOut = META.TIME.idleTimeOut - (META.TIME.currentTime - referenceTime);
        if timeOut <= 0
           fprintf('[%s]\t... Aborting.\n',META.phase);
           break
        end
        fprintf('[%s]\t... All agents are idle, aborting in t=%0.2fs.\n',META.phase,timeOut);
    end
    % /////////////////////////////////////////////////////////////////////

    % 1. ////////////// UPDATE TIMESTEP PARAMETERS (@t=k) /////////////////
    META.TIME.currentTime = META.TIME.timeVector(step);                    % The current sim-time
    META.TIME.currentStep = step;                                          % The current sim-step
    if META.verbosity >= 1
        fprintf('[%s]\tStep: %.0f\tTime: %.2fs\n',META.phase,META.TIME.currentStep,META.TIME.currentTime);
    end
    % /////////////////////////////////////////////////////////////////////

    % 2. ///////////////// RECORD GLOBAL STATE (@t=k) /////////////////////
    for ID1 = 1:META.totalObjects
        % Collect the META.OBJECTS.state data (as the simulations understanding of the global states)
        % and save to the output DATA.globalTrajectories set, must be done synchronously).
        DATA.globalTrajectories(DATA.stateIndex(1,ID1):DATA.stateIndex(2,ID1),META.TIME.currentStep) = META.OBJECTS(ID1).globalState;
    end
    % /////////////////////////////////////////////////////////////////////

    % 3. ////////////// UPDATE THE GLOBAL STATES (@t=k) ////////////////////
    for ID1 = 1:META.totalObjects                                          % For each META object element
        % Update META.OBJECT the structure
%         META.OBJECTS(ID1) = OMAS_updateGlobalStates_mex(...
%             META,...
%             objectIndex{ID1}.objectID,...
%             objectIndex{ID1}.GetGLOBAL('velocity'),...
%             objectIndex{ID1}.GetGLOBAL('quaternion'),...
%             objectIndex{ID1}.GetGLOBAL('idleStatus'));
        META.OBJECTS(ID1) = OMAS_updateGlobalStates(...
            META,...
            objectIndex{ID1}.objectID,...
            objectIndex{ID1}.GetGLOBAL('velocity'),...
            objectIndex{ID1}.GetGLOBAL('quaternion'),...
            objectIndex{ID1}.GetGLOBAL('idleStatus'));
    end
    % /////////////////////////////////////////////////////////////////////

%     % 4. ////////// UPDATE THE VISUAL REPRESENTATIONS (@t=k) //////////////
%     for ID1 = 1:META.totalObjects
%         if ~isstruct(META.OBJECTS(ID1).geometry)
%             continue
%         end
%         % MODIFY THE ASSOCIATED GEOMETRY
%         META.OBJECTS(ID1) = updateVisuals(META.OBJECTS(ID1));
%     end

    % 5. //////// UPDATE SIMULATION/ENVIRONMENTAL META DATA (@t=k) ////////
    [META,metaEVENTS] = UpdateSimulationMeta(META,objectIndex);            % Update META snapshot with equivilant objectIndex.state
    % LOG THE META EVENTS
    if ~isempty(metaEVENTS)                                                % If META events occurred this timestep
        EVENTS = vertcat(EVENTS,metaEVENTS);                               % Append to global EVENTS
    end
    % /////////////////////////////////////////////////////////////////////

    % 6. //// COMPUTE AGENT CYCLES (UPDATE GLOBAL REPRESENTATION (@t=k) ///
    if META.threadPool ~= 0
        objectSnapshot = objectIndex;                                      % Make a temporary record of the object set
        parfor (ID1 = 1:META.totalObjects)
            % MOVE THROUGH OBJECT INDEX AND UPDATE EACH AGENT
            [objectIndex{ID1},objectEVENTS] = UpdateObjects(META,objectSnapshot,objectIndex{ID1});
            % LOG THE OBJECT EVENTS
            if ~isempty(objectEVENTS)
                EVENTS = vertcat(EVENTS,objectEVENTS);
            end
        end
    else
        for ID1 = 1:META.totalObjects
            % MOVE THROUGH OBJECT INDEX AND UPDATE EACH AGENT
            [objectIndex{ID1},objectEVENTS] = UpdateObjects(META,objectIndex,objectIndex{ID1}); % Update objectIndex snapshot with new META data
            % LOG THE OBJECT EVENTS
            if ~isempty(objectEVENTS)                                      % If objectEVENTS occur in this timestep
                EVENTS = vertcat(EVENTS,objectEVENTS);                     % Append to global EVENTS
            end
        end
    end
    % /////////////////////////////////////////////////////////////////////

    % 7. /// THE 'OBJECT.VIRTUAL' PROPERTIES IS NOW UPDATED FOR (t=k+1) ///
    step = step + 1;
end
% CREATE TERMINAL VALUES, FOR CLARITY
META.TIME.endStep = META.TIME.currentStep;
META.TIME.endTime = META.TIME.currentTime;
end
% UPDATE THE OBJECT META PROPERTIES
function [SIM,metaEVENTS]	= UpdateSimulationMeta(SIM,objectIndex)
% In this function we wish to update the global separations and event
% conditions based on the new 'globalState' properties.

% INPUTS:
% SIM         - A local copy of the global META structure
% objectIndex - The cell array of object classes.
% OUTPUTS:
% SIM         - The updated META structure
% metaEVENTS  - A vector of new event objects.

% We want to move through the object set an update both the object being
% updated, but also the object it is being updated against.

% //// ASSESS THE COLLISIONS AND RELATIVE POSITIONS (SYMMETRICAL CHECK) ///
collisionLogicals = zeros(SIM.totalObjects);
warningLogicals   = zeros(SIM.totalObjects);
for entityA = 1:SIM.totalObjects                                           % Object A's position in the META.OBJECTS

    % NOTE:
    % - We only need to consider the objects whos separations have yet to be
    %   evaluated. If we update 1&2 and 2&1 simultaneously, then we must only
    %   examine the unique ID permutations to assess their conditions.

    for entityB = (entityA+1):(SIM.totalObjects)                           % Object B's position in the META.OBJECTS
        % THE OBJECTS ASSOCIATED objectIndex
        object_A = objectIndex{SIM.globalIDvector == SIM.OBJECTS(entityA).objectID};
        object_B = objectIndex{SIM.globalIDvector == SIM.OBJECTS(entityB).objectID};  % Determine the associated geometries

        % UPDATE RELATIVE POSITIONS  (centroid separation based)
        SIM.OBJECTS(entityA).relativePositions(entityB,:) = (SIM.OBJECTS(entityB).globalState(1:3,1) - SIM.OBJECTS(entityA).globalState(1:3,1))';   % AB vector
        SIM.OBJECTS(entityB).relativePositions(entityA,:) = (SIM.OBJECTS(entityA).globalState(1:3,1) - SIM.OBJECTS(entityB).globalState(1:3,1))';   % BA vector

        % GET THE WARNING CONDITIONS (centroid separation based)
        separationDistance = norm(SIM.OBJECTS(entityA).relativePositions(entityB,:)) - (SIM.OBJECTS(entityA).radius + SIM.OBJECTS(entityB).radius);
        warningCondition = SIM.warningDistance >= separationDistance;
        warningLogicals(entityA,entityB) = warningCondition;
        warningLogicals(entityB,entityA) = warningCondition;               % Assign the warning condition

        % EVALUATE COLLISIONS BETWEEN THE OBJECTS
        collisionCondition = OMAS_collisionDetection(...
            SIM.OBJECTS(entityA),object_A.GEOMETRY,...
            SIM.OBJECTS(entityB),object_B.GEOMETRY,...
            SIM.conditionTolerance);
        collisionLogicals(entityA,entityB) = collisionCondition;
        collisionLogicals(entityB,entityA) = collisionCondition;           % Assign the collision condition

        % ERROR CHECKING (Should never happen)
        if separationDistance > 0 && collisionCondition
            error('[ERROR] A collision occurred at distance %f without violating the minimum separation of %f',...
                    norm(SIM.OBJECTS(entityA).relativePositions(entityB,:)),(SIM.OBJECTS(entityA).radius + SIM.OBJECTS(entityB).radius));
        end
    end
end

% //////////// ASSESS DETECTION CONDITIONS (ASYMMETRICAL CHECK) ///////////
% DEFAULT DETECTION CONDITION
detectionLogicals = zeros(SIM.totalObjects);
for entityA = 1:SIM.totalObjects
    % Check objects have the capacity to observe other objects
    if SIM.OBJECTS(entityA).type ~= OMAS_objectType.agent
        continue
    end

    % Check all agent observations
    for entityB = 1:SIM.totalObjects                                       % Object B's position in the META.OBJECTS
        % SKIP SELF-REFERENCE
        if SIM.OBJECTS(entityA).objectID == SIM.OBJECTS(entityB).objectID
            continue                                                       % Only agents can generate notifications
        end
        % NOTE:
        % - We also need to determine which agents are able to observe other
        %   agents in the field. However, this is NOT A SYMMETRIC operation
        %   as A may be able to detect B without B detecting A.

        % Detection check #1 - Authorisation
        isWaypoint = SIM.OBJECTS(entityB).type == OMAS_objectType.waypoint;
        if isWaypoint
            isAuthorised = objectIndex{SIM.globalIDvector == SIM.OBJECTS(entityB).objectID}.IDAssociationCheck(SIM.OBJECTS(entityA).objectID);
            if ~isAuthorised
                continue
            end
        end

        % Detection check #2 - Spherical
        isDetected = OMAS_geometry.intersect_spheres(...
            SIM.OBJECTS(entityA).globalState(1:3),SIM.OBJECTS(entityA).detectionRadius + 0.5*SIM.conditionTolerance,...
            SIM.OBJECTS(entityB).globalState(1:3),SIM.OBJECTS(entityB).radius + 0.5*SIM.conditionTolerance);

        % Object is not detected
        if ~isDetected
            continue
        end

        % This check is sufficient for way-points
        if isWaypoint
            detectionLogicals(entityA,entityB) = 1;
            continue
        end

        % ///////////////// OBJECT MAY BE OBSERVED ////////////////////////
        % NOTE:
        % - If the distance(less the radius) is less than detection radius
        %   then its worth checking for more a more complex representation.
        % - We want to check if any vertices/edges can be observed. If any
        %   edges are visible, then they are to be sent to the first object.

        % Get the geometry of the second object (GEOMETRY OF B)
        geometryB = objectIndex{SIM.globalIDvector == SIM.OBJECTS(entityB).objectID}.GEOMETRY;

        % Check if a better check is necessary
        hasGeometry = isstruct(geometryB) && size(geometryB.vertices,1) > 0;
        if isDetected && ~hasGeometry
            detectionLogicals(entityA,entityB) = 1;                        % Detection occurred
            break
        end
        isDetected = 0; % Reset, potential false positive when checking against polygon

        % Assess the polygon
        relativeR = SIM.OBJECTS(entityA).R'*SIM.OBJECTS(entityB).R;        % Rotate second geometry orientated relative to
        % Design the vertex set
        geometryB.vertices = geometryB.vertices*relativeR + SIM.OBJECTS(entityB).globalState(1:3)';

        % Assess each face and vertex for detection
        for face = 1:size(geometryB.vertices,1)
            % GET THE MEMBER IDs OF THE FACE
            faceMembers  = geometryB.faces(face,:);
            faceVertices = geometryB.vertices(faceMembers,:);

            % CHECK THE GEOMETRY AGAINST THE SPHERICAL CONSTRAINT (with mex)
%             isDetected = CheckSphereTriangleIntersection_mex(...
%                 SIM.OBJECTS(entityA).globalState(1:3),...
%                 SIM.OBJECTS(entityA).detectionRadius,...
%                 faceVertices(1,:)',faceVertices(2,:)',faceVertices(3,:)');
            % CHECK THE GEOMETRY AGAINST THE SPHERICAL CONSTRAINT
            isDetected = CheckSphereTriangleIntersection(...
                SIM.OBJECTS(entityA).globalState(1:3),...
                SIM.OBJECTS(entityA).detectionRadius,...
                faceVertices(1,:)',faceVertices(2,:)',faceVertices(3,:)');

            % IF THE FACE IS DETECTED
            if isDetected
                detectionLogicals(entityA,entityB) = 1;
                break
            end
        end

        % NOTE:
        % - If any faces are within the constraint, then the detection
        %   logical is set true and then exited.
        % - Otherwise, if the none of the contraints are met, then it
        %   remains undetected.
    end
end

% NOTES:
% - The fields geometry based fields of SIM.OBJECTS are updated to this
%   time-step. The event conditions can now be re-evaluted.
% - CHECKS MUST BE ASYMMETRICAL - A may detect B without B detecting A

metaEVENTS = [];                                                           % Reset meta Events container
for entityA = 1:SIM.totalObjects
    % CHECK OBJECTS HAVE THE CAPACITY TO OBSERVE OTHER OBJECTS
    if SIM.OBJECTS(entityA).type ~= OMAS_objectType.agent
        continue
    end

    for entityB = 1:SIM.totalObjects
        % OMIT SELF-CHECK
        if SIM.OBJECTS(entityA).objectID == SIM.OBJECTS(entityB).objectID
            continue                                                       % Only agents can generate notifications
        end

        % //////////// ASSESS EVENT CONDITIONS FOR THE AGENTS /////////////
        % NOTES:
        % - The agent set is only capable of generating events associated
        %   with the detections, warnings and collisions.

        % ///////////////// ASSESS DETECTIONS CONDITIONS //////////////////
        % DETECTION EVENT LOGIC
        ABdetectionConstraint       =  detectionLogicals(entityA,entityB);
        novelDetectionConstraint    =  SIM.OBJECTS(entityA).objectStatus(entityB,eventType.detection) == 0;
        detectionEventCondition     =  ABdetectionConstraint &&  novelDetectionConstraint;   % Is within a range and a detection event has not been issued.
        detectionEventNullCondition = ~ABdetectionConstraint && ~novelDetectionConstraint;   % Is outside a range and a detection condition is still toggled.
        % EVALUATE OCCURANCE
        if detectionEventCondition
            % UPDATE DETECTION STATUS
            SIM.OBJECTS(entityA).objectStatus(entityB,eventType.detection) = 1;              % Ammend status of the META object
            % GENERATE DETECTION EVENT
            [EVENT] = OMAS_eventHandler(SIM.TIME.currentTime,SIM.OBJECTS(entityA),SIM.OBJECTS(entityB),eventType.detection);
            metaEVENTS = vertcat(metaEVENTS,EVENT);                                          % Add new META EVENTS to global structure
        elseif detectionEventNullCondition
            % AMEND THE META OBJECT
            SIM.OBJECTS(entityA).objectStatus(entityB,eventType.detection) = 0;              % Ammend status of the META object
            % GENERATE THE DETECTION NULLIFICATION EVENT
            [EVENT] = OMAS_eventHandler(SIM.TIME.currentTime,SIM.OBJECTS(entityA),SIM.OBJECTS(entityB),eventType.null_detection);
            metaEVENTS = vertcat(metaEVENTS,EVENT);                                          % Add new META EVENTS to global structure
        end

        % ///////////////// ASSESS WAYPOINT CONDITIONS ////////////////////
        if SIM.OBJECTS(entityB).type == OMAS_objectType.waypoint
            % Check A is allowed to get B
            waypointObj = objectIndex{SIM.globalIDvector == SIM.OBJECTS(entityB).objectID};
            hasAssociation = waypointObj.IDAssociationCheck(SIM.OBJECTS(entityA).objectID);
            % WAYPOINT EVENT LOGIC
            ABwaypointCondition        =  collisionLogicals(entityA,entityB) && hasAssociation;
            novelWaypointGetConstraint =  SIM.OBJECTS(entityA).objectStatus(entityB,eventType.waypoint) == 0;
            waypointGetEventCondition  =  ABwaypointCondition &&  novelWaypointGetConstraint;    % Is within a range and a collision event has not been issued.
            waypointEventNullCondition = ~ABwaypointCondition && ~novelWaypointGetConstraint;    % Is outside a range and a waypoint-get condition is still toggled.
            % EVALUATE OCCURANCE
            if waypointGetEventCondition
                % AMEND THE META OBJECT
                SIM.OBJECTS(entityA).objectStatus(entityB,eventType.waypoint) = 1;
                % GENERATE THE WAYPOINT ACHIEVED EVENT
                [EVENT] = OMAS_eventHandler(SIM.TIME.currentTime,SIM.OBJECTS(entityA),SIM.OBJECTS(entityB),eventType.waypoint);
                metaEVENTS = vertcat(metaEVENTS,EVENT);                                          % Add new META EVENTS to global structure
            elseif waypointEventNullCondition
                % AMEND THE META OBJECT
                SIM.OBJECTS(entityA).objectStatus(entityB,eventType.waypoint) = 0;
                % GENERATE THE WAYPOINT NULLIFICATION EVENT
                %[EVENT] = OMAS_eventHandler(SIM.TIME.currentTime,SIM.OBJECTS(entityA),SIM.OBJECTS(entityB),eventType.null_waypoint);
                %metaEVENTS = vertcat(metaEVENTS,EVENT);                                        % Add new META EVENTS to global structure
            end
            continue    % Skip collision check
        end

        % ////////// CONFIRM THE OBJECT CAN BE COLLIDED WITH //////////////
        if SIM.OBJECTS(entityB).hitBox == OMAS_hitBoxType.none
            % Check hit-box definition:
            % -If no hitbox is assigned, no need to evaluate warnings or
            % collisions.
            continue
        end

        % /////////////////// ASSESS WARNING CONDITIONS ///////////////////
        % WARNING EVENT LOGIC
        ABwarningConstraint       =  warningLogicals(entityA,entityB);
        novelWarningConstraint    =  SIM.OBJECTS(entityA).objectStatus(entityB,eventType.warning) == 0;
        warningEventCondition     =  ABwarningConstraint &&  novelWarningConstraint;         % Is within a range and a warning event has not been issued.
        warningEventNullCondition = ~ABwarningConstraint && ~novelWarningConstraint;         % Is outside a range and a warning condition is still toggled.
        % EVALUATE OCCURANCE
        if warningEventCondition
            % UPDATE WARNING STATUS
            SIM.OBJECTS(entityA).objectStatus(entityB,eventType.warning) = 1;
            % GENERATE THE WARNING EVENT
            [EVENT] = OMAS_eventHandler(SIM.TIME.currentTime,SIM.OBJECTS(entityA),SIM.OBJECTS(entityB),eventType.warning);
            metaEVENTS = vertcat(metaEVENTS,EVENT);                                          % Add new META EVENTS to global structure
        elseif warningEventNullCondition
            % UPDATE WARNING-NULL STATUS
            SIM.OBJECTS(entityA).objectStatus(entityB,eventType.warning) = 0;
            % GENERATE THE WARNING NULLIFICATION EVENT
            [EVENT] = OMAS_eventHandler(SIM.TIME.currentTime,SIM.OBJECTS(entityA),SIM.OBJECTS(entityB),eventType.null_warning);
            metaEVENTS = vertcat(metaEVENTS,EVENT);                                          % Add new META EVENTS to global structure
        end

        % ////////////////// ASSESS COLLISION CONDITIONS //////////////////
        % COLLISION EVENT LOGIC
        ABcollisionConstraint       =  collisionLogicals(entityA,entityB);
        novelCollisionConstraint    =  SIM.OBJECTS(entityA).objectStatus(entityB,eventType.collision) == 0;
        collisionEventCondition     =  ABcollisionConstraint &&  novelCollisionConstraint;   % Is within a range and a collision event has not been issued.
        collisionEventNullCondition = ~ABcollisionConstraint && ~novelCollisionConstraint;   % Is outside a range and a collision condition is still toggled.
        % EVALUATE OCCURANCE
        if collisionEventCondition
            % UPDATE COLLISIONL STATUS
            SIM.OBJECTS(entityA).objectStatus(entityB,eventType.collision) = 1;
            % GENERATE THE COLLISION EVENT
            [EVENT] = OMAS_eventHandler(SIM.TIME.currentTime,SIM.OBJECTS(entityA),SIM.OBJECTS(entityB),eventType.collision);
            metaEVENTS = vertcat(metaEVENTS,EVENT);                                          % Add new META EVENTS to global structure
        elseif collisionEventNullCondition
        	% AMEND THE META OBJECT
            SIM.OBJECTS(entityA).objectStatus(entityB,eventType.collision) = 0;
            % GENERATE THE COLLISION NULLIFICATION EVENT
            [EVENT] = OMAS_eventHandler(SIM.TIME.currentTime,SIM.OBJECTS(entityA),SIM.OBJECTS(entityB),eventType.null_collision);
            metaEVENTS = vertcat(metaEVENTS,EVENT);                                          % Add new META EVENTS to global structure
        end
    end
end

end
% UPDATE THE OBJECT PROPERTIES
function [referenceObject,objectEVENTS] = UpdateObjects(SIM,objectIndex,referenceObject)
% This function updates a referenceObject class against the rest of the
% objectIndex, independantly of the SIM.OBJECTS META data.
% INPUTS:
% SIM             - Current META structure
% objectIndex     - The current object list
% referenceObject - The agent object being updated

% OUTPUTS:
% referenceObject - The updated agent class

objectEVENTS = []; % Container for object based events

% GET THE OBJECTS EQUIVALENT SIM OBJECT
SIMfirstObject = SIM.OBJECTS(SIM.globalIDvector == referenceObject.objectID);

% If the observing object is "2D"
if referenceObject.Is3D()
    dimensionIndices = 1:3;
else
    dimensionIndices = 1:2;
end

if SIM.verbosity >= 2
    fprintf('[UPDATE]\tCycling ID:%d\t name: %s\n',SIMfirstObject.objectID,referenceObject.name);

end

% Default timing parameters
ENV = SIM.TIME;
ENV.outputPath = SIM.outputPath;

if SIM.verbosity >= 2
    fprintf('[UPDATE]\tCycling ID:%d\t name: %s\n',SIMfirstObject.objectID,referenceObject.name);
end

% SWITCH BEHAVIOUR BASED ON OBJECT TYPE
switch SIMfirstObject.type
    case OMAS_objectType.agent
        % AGENT - SIMULATION/ENVIROMENTAL FEEDBACK REQUIRED %%%%%%%%%%%%%%
        % Container for agent detection packet
        observationPacket = [];

        % //// PARSE THE ENVIRONMENTAL OBJECTS OBSERVABLE TO THE AGENT ////
        for ID2 = 1:SIM.totalObjects
            % GET THE SIM OBJECT OF THE EVALUATION OBJECT
            SIMsecondObject = SIM.OBJECTS(ID2);

            % Skip condition - Get the current status of this object
            isDetected = 1 == SIMfirstObject.objectStatus(ID2,eventType.detection);
            isSameObject = SIMfirstObject.objectID == SIMsecondObject.objectID;
            if ~isDetected || isSameObject
                continue
            end

            % GET THE EQUIVALENT OBJECT TO EVALUATE AGAINST THE REFERENCE
            secondObject = objectIndex{SIM.globalIDvector == SIMsecondObject.objectID};

            % ///////////////// ELSE AGENT IS DETECTED ////////////////////
            % ////// (GENERATE A LOCALISED PROXIMITY DESCIPTION) //////////
            % The second agent is within global seperation conditions. A
            % communication packet is therefore assembled containing the
            % data on the second object. rotated into the agent-local
            % coordinate frame.

            % GET THE RELATIVE POSITIONS OF THE AGENT AND OBJECT IN THE
            % GLOBAL COORDINATE SYSTEM (x_rel = x_B - x_A)
            relativePosition = SIMfirstObject.relativePositions(ID2,:)';
            relativeVelocity = SIMsecondObject.globalState(4:6) - SIMfirstObject.globalState(4:6);
            % ROTATE THE GLOBAL STATE OF THE OBJECT INTO THE AGENT FRAME
            observedPosition = SIMfirstObject.R*relativePosition;            % Rotate the from the global into the body frame of the simReference
            observedVelocity = SIMfirstObject.R*relativeVelocity;
            % SPHERICAL REPRESENTATION
            observedRange     = norm(observedPosition);
            observedElevation = asin(observedPosition(3)/observedRange);
            observedHeading   = atan2(observedPosition(2),observedPosition(1));
            % PRESENT SIZE PROPERTIES
            observedRadius    = SIMfirstObject.radius;
            observedAngularWidth = 2*asin(observedRadius/(observedRange + observedRadius));

            % OBJECT GEOMETRY PREPARATION
            % In order for the agent to observe the object correctly. The
            % geometry must be translated and rotated into the relative
            % frame of the detecting agent.
            % ASSUMPTION:
            % - The geometry is normalised to the body axes of the object.
            % TO DO:
            % - If the object is on the edge of vision.. only part of the
            %   object will be visable to the agent. a process must be in
            %   place to create a subset of the geometry.

            % THE GEOMETRIC PARAMETERS
            if size(secondObject.GEOMETRY.vertices,1) > 0
                % Rotate evaluation object geometry into the global space,
                % then rotate it back into the local space of the reference.
                % THE RELATIVE ROTATIONS OF THE SECOND BODY
                relativeR = SIMfirstObject.R'*SIMsecondObject.R;
                % THE CONTAINER FOR THE RELATIVE GEOMETRY
                observedGeometry = struct('vertices',secondObject.GEOMETRY.vertices*relativeR + observedPosition',...
                                          'normals', secondObject.GEOMETRY.normals*relativeR,...
                                          'faces',   secondObject.GEOMETRY.faces,...
                                          'centroid',secondObject.GEOMETRY.centroid + observedPosition');
                % [TO-DO] Reduce the structure sent to match the exposed geometry

%                 % PROCESS THE VISIBLE GEOMETRY (The sub-set of the geometry that is visible to the agent.)
%                 observedGeometry = OMAS_restrictedGeometry(zeros(3,1),SIMfirstObject.detectionRadius,relativeGeometry);
            else
                % THE OBSERVED GEOMETRY IS EMPTY
                observedGeometry = secondObject.GEOMETRY;                  % Pass the empty structure
            end

            % OBJECT PRIORITY (IF WAYPOINT)
            observedPriority = NaN;
            if SIMsecondObject.type == OMAS_objectType.waypoint
                association = objectIndex{ID2}.GetAgentAssociation(referenceObject);
                observedPriority = association.priority;
            end

            % A DEBUG STUCTURE
            DEBUG = struct('globalPosition',SIMsecondObject.globalState(dimensionIndices),...     % Added for simplicity
                           'globalVelocity',SIMsecondObject.globalState(dimensionIndices + 3),... % Added for simplicity
                           'priority',observedPriority);                                          % The objects global priority for that agent

            % ///////////////// PREPARE DETECTION PACKET //////////////////
            detectionObject = struct('objectID',SIMsecondObject.objectID,...            % The object ID (observed)
                                     'name',SIMsecondObject.name,...                    % The object name tag (observed)
                                     'type',SIMsecondObject.type,...                    % The objects sim-type enum
                                     'radius',observedRadius,...                        % The objects true size
                                     'position',observedPosition(dimensionIndices,1),...% The apparent position in the relative frame
                                     'velocity',observedVelocity(dimensionIndices,1),...% The apparent velocity in the relative frame
                                     'range',observedRange,...                          % The apparent range
                                     'elevation',observedElevation,...                  % The apparent inclination angle
                                     'heading',observedHeading,...                      % The apparent Azimuth angle
                                     'width',observedAngularWidth,...                   % The apparent angular width at that range
                                     'geometry',observedGeometry,...                    % The observable geometrical components
                                     'colour',SIMsecondObject.colour,...                % Finally the simulation's colourID
                                     'DEBUG',DEBUG);                                    % Pass additional parameters (not otherwise known)
            observationPacket = vertcat(observationPacket,detectionObject);             % Append object to packet to agent
        end

        % /////////////// COMPUTE THE LOCAL AGENT CYCLE ///////////////////
        % Given the detection object defined for this agent, compute the
        % agents cycle with this information.
        try
            % SEND OBJECT OBSERVERATION PACKET TO AGENT
            referenceObject = referenceObject.main(ENV,observationPacket);                      % Hand the current sim TIME object and complete observation packet
        catch objectProcessError
            warning('[%s]: Cycle failed, there was a problem with the object (%s) file, please check: %s',...
                     SIM.phase,referenceObject.name,objectProcessError.stack(1).name);
            rethrow(objectProcessError);
        end
    case OMAS_objectType.obstacle
        % PASSIVE - NO FEEDBACK REQUIRED %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        referenceObject = referenceObject.main(ENV);
    case OMAS_objectType.waypoint
        % WAYPOINTS ARE CONSIDERED PASSIVE - NO FEEDBACK REQUIRED %%%%%%%%%
        referenceObject = referenceObject.main(ENV);
    otherwise
        % CONSIDERED PASSIVE - NO FEEDBACK REQUIRED %%%%%%%%%%%%%%%%%%%%%%%
        % DO NOT UPDATE
end
end

% /////////////////// SIMULATION OUTPUT OPERATIONS ////////////////////////
% GET INITIAL OUTPUT DATA STRUCTURE
function [DATA] = GetOutputStructure(SIM)
% This function assembles the initial output DATA structure from the
% simulation input variables and it back as a simulation output container.
% INPUT:
% SIM  - A local reference to the META structure
% OUPUT:
% DATA - The output data structure

% CALCULATE THE INITIAL DATA PARAMETERS
globalStateNum = size(SIM.OBJECTS(1).globalState,1);
systemStates = globalStateNum*SIM.totalObjects;                            % The total number of system states
% DEFINE THE GLOBAL STATE INDICES
indexSet = zeros(2,length(SIM.OBJECTS));                                   % Indices are defined [start;end]*n
indexSet(1,:) = (0:globalStateNum:systemStates-globalStateNum) + 1;        % Declare the state indices
indexSet(2,:) = globalStateNum:globalStateNum:systemStates;

% GENERATE THE OUTPUT STRUCTURE FROM THE SIMULATION INPUT PARAMETERS
DATA = struct('outputPath',[SIM.outputPath,'DATA.mat'],...
               'outputDir',SIM.outputPath,...
            'totalObjects',SIM.totalObjects,...
             'totalAgents',SIM.totalAgents,...
          'totalObstacles',SIM.totalObstacles,...
          'totalWaypoints',SIM.totalWaypoints,...
              'timeVector',SIM.TIME.timeVector,...
                      'dt',SIM.TIME.dt,...
              'stateIndex',indexSet,...
      'globalTrajectories',NaN(systemStates,SIM.TIME.numSteps));         % Prepare the output container
%       'globalTrajectories',zeros(systemStates,SIM.TIME.numSteps));         % Prepare the output container
end

repo: OpenMAS


text:
function[lat,lon]=uv2latlon(varargin)
%UV2LATLON  Integrates horizontal velocity to give latitude and longitude.
%
%   [LAT,LON]=UV2LATLON(NUM,U,V,LATO,LONO) where NUM is the data in DATENUM
%   format, integates velocities U and V from intial location (LATO,LONO)
%   to form a trajectory on the earth described by LAT and LON.
%
%   LATO and LONO are an initial latitude and longitude in degrees,
%   and U and V are eastward and northward velocity components in cm/s.
%
%   NUM, U, and V may be column vectors, in which case LATO and LONO are
%   both scalars.  Alternatively U and V may be matrices with time oriented
%   in rows.  NUM is then either an array of length SIZE(U,1), or a matrix
%   of the same size as U and V, while LATO and LONO are either scalars or
%   arrays of length SIZE(U,2).
%
%   [LAT,LON]=UV2LATLON(NUM,CV,LATO,LONO), where CV is the complex-valued
%   velocity CV=U+SQRT(-1)*V, also works.
%
%   UV2LATLON is inverted by LATLON2UV.
%   ___________________________________________________________________
%
%   Cell array input / output
%
%   UV2LATLON returns cell array output given cell array input.
%
%   That is, if NUM, U, and V, are all cell arrays of length K, containing
%   K different numerical arrays, then the output will also be cell arrays
%   of length K.
%
%   In this case LATO and LONO are either scalars, or arrays of length K.
%   ___________________________________________________________________
%
%   Algorithm
%
%   UV2LATLON works by converting the velocity at the current point into
%   a vector in 3D Cartesian coordinates, using SPHERE2UVW, finding the new
%   location in 3D space, and then converting this back into latitude
%   and longitude using XYZ2LATLON, and then interating for the next point.
%
%   By default, UV2LATLON uses a forward integration from an initial point.
%   This inverts LATLON2UV(...,'forward') to a high degree of precision for
%   typical drifter and float velocity values and sampling intervals.
%
%   UV2LATLON(NUM,CV,LATF,LONF,'backward') instead uses a *backward*
%   integration, and LATF and LONF are now the *final* points rather than
%   the initial points.  This inverts LATLON2UV(...,'backward').
%   ___________________________________________________________________
%
%   See also XY2LATLON, LATLON2XY, LATLON2UV.
%
%   'uv2latlon --t' runs a test.
%
%   Usage: [lat,lon]=uv2latlon(num,u,v,lato,lono);
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2014 J.M. Lilly --- type 'help jlab_license' for details

if strcmpi(varargin{1}, '--t')
    uv2latlon_test,return
end

if ischar(varargin{end})
    str=varargin{end};
    varargin=varargin(1:end-1);
else
    str='forward';
end

na=length(varargin);
num=varargin{1};
if na==4
    if iscell(varargin{2})
        u=cellreal(varargin{2});
        v=cellimag(varargin{2});
    else
        u=real(varargin{2});
        v=imag(varargin{2});
    end
    varargin=varargin(3:end);
else
    u=varargin{2};
    v=varargin{3};
    varargin=varargin(4:end);
end
lato=varargin{1};
lono=varargin{2};

if ~iscell(num)
    if ~aresame(size(num),size(u))
        num=vrep(num(:),size(u,2),2);
    end
    if length(lato)==1
        lato=lato+zeros(size(u(1,:)));
    end
    if length(lono)==1
        lono=lono+zeros(size(u(1,:)));
    end
    %vsize(num,lato,lono,u,v)
    [lat,lon]=uv2latlon_one(num,lato,lono,u,v,str);
else
    if length(lato)==1
        lato=lato+zeros(size(u));
    end
    if length(lono)==1
        lono=lono+zeros(size(u));
    end
    for i=1:length(num)
        [lat{i,1},lon{i,1}]=uv2latlon_one(num{i},lato(i),lono(i),u{i},v{i},str);
    end
end

function[lat,lon]=uv2latlon_one(num,lato,lono,u,v,str)

switch str(1:3)
    case 'for'
        [lat,lon]=uv2latlon_integrator(num,lato,lono,u,v);
    case 'bac'
        [lat,lon]=uv2latlon_integrator(num,lato,lono,-flipud(u),-flipud(v));
        lat=flipud(lat);
        lon=flipud(lon);
%     case 'cen'
%         [lat1,lon1]=uv2latlon_integrator(num,lato(1),lono(1),u,v);
%         [lat2,lon2]=uv2latlon_integrator(num,lato(2),lono(2),-flipud(u),-flipud(v));
%         lat2=flipud(lat2);
%         lon2=flipud(lon2);
%         [x1,y1,z1]=latlon2xyz(lat1,lon1);
%         [x2,y2,z2]=latlon2xyz(lat2,lon2);
%         [lat,lon]=xyz2latlon((x1+x2)/2,(y1+y2)/2,(z1+z2)/2);
end

function[lat,lon]=uv2latlon_integrator(num,lato,lono,u,v)

dt=num(2)-num(1);
[xo,yo,zo]=latlon2xyz(lato,lono);
[lat,lon]=vzeros(size(u));
lat(1,:)=lato;
lon(1,:)=lono;

u=u(1:end-1,:);
v=v(1:end-1,:);

% if strcmpi(str(1:3),'for')
%     u=u(1:end-1,:);
%     u=u(1:end-1,:);
% elseif strcmpi(str(1:3),'for')
%     u=u(1:end-1,:);
%     u=u(1:end-1,:);
%
% u=1/2*(u+circshift(u,-1));
% u=u(1:end-1,:);
%
% v=1/2*(v+circshift(v,-1));
% v=v(1:end-1,:);

u=u*(dt*24*3600)/100/1000;%Now it's km
v=v*(dt*24*3600)/100/1000;%Now it's km

%vsize(lat,lon,u,v)
for i=1:size(u,1)
    [dx,dy,dz]=sphere2uvw(lat(i,:),lon(i,:),0,u(i,:),v(i,:));
    xo=xo+dx;
    yo=yo+dy;
    zo=zo+dz;
    [lat(i+1,:),lon(i+1,:)]=xyz2latlon(xo,yo,zo);
end

function[]=uv2latlon_test
load npg2006

use npg2006
vindex(num,lat,lon,1:max(find(~isnan(lat))),1);
cv=latlon2uv(num,lat,lon,'forward');
[lat2,lon2]=uv2latlon(num,cv,lat(1),lon(1),'forward');

bool1=allall(abs(frac(lon-lon2,lon))<5e-6);
bool2=allall(abs(frac(lat-lat2,lat))<5e-6);
reporttest('UV2LATLON inverts LATLON2UV to within 5e-6 for NPG2006 data, forward algorithm',bool1&&bool2)

[lat2,lon2]=uv2latlon([num num],[cv cv],[lat(1) lat(1)],[lon(1) lon(1)],'forward');
bool1=allall(abs(frac([lon lon]-lon2,[lon lon]))<5e-6);
bool2=allall(abs(frac([lat lat]-lat2,[lat lat]))<5e-6);
reporttest('UV2LATLON inverts LATLON2UV to within 5e-6 for NPG2006 data, forward algorithm, matrix input',bool1&&bool2)

use npg2006
vindex(num,lat,lon,1:max(find(~isnan(lat))),1);
cv=latlon2uv(num,lat,lon,'backward');
[lat3,lon3]=uv2latlon(num,cv,lat(end),lon(end),'backward');

bool1=allall(abs(frac(lon-lon3,lon))<5e-6);
bool2=allall(abs(frac(lat-lat3,lat))<5e-6);
reporttest('UV2LATLON inverts LATLON2UV to within 5e-6 for NPG2006 data, backward algorithm',bool1&&bool2)

% use npg2006
% vindex(num,lat,lon,1:max(find(~isnan(lat))),1);
% cv=latlon2uv(num,lat,lon);
% [lat4,lon4]=uv2latlon(num,cv,[lat(1) lat(end)],[lon(1) lon(end)],'central');
%
% bool1=allall(abs(frac(lon-lon4,lon))<5e-6);
% bool2=allall(abs(frac(lat-lat4,lat))<5e-6);
% reporttest('UV2LATLON inverts LATLON2UV to within 5e-6 for NPG2006 data, central algorithm',bool1&&bool2)

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab