text:
function [imX,imY,imZ] = dtiSaveImageSlicesOverlays(handles, overlayFGs, overlayROIs, queryParameters, fname, upSamp, curPosAcpc, bg)
%
%   dtiSaveImageSlicesOverlays(handles,[overlayFGs=[]],[overlayROIs=[]],[queryParameters=0],[fname],[upSamp])
%
%Author: Dougherty, Wandell
%Purpose:
%   Save out image slices.  Potentially we overlay the FGs and ROIs.
%   Potentially we query for the image processing parameters
%
% Examples:
%   dtiSaveImageSlicesOverlays(handles)       %Just slices, no overlays
%   dtiSaveImageSlicesOverlays(handles,1,0)   % Overlay fiber groups
%   dtiSaveImageSlicesOverlays(handles,1,0,1) % Overlay fiber groups and
%              query for imaging parameters

if ieNotDefined('overlayFGs'),  overlayFGs = []; end
if ieNotDefined('overlayROIs'), overlayROIs = []; end
if ieNotDefined('queryParameters'), queryParameters = 1; end
if(~exist('upSamp','var') | isempty(upSamp)) upSamp = 4; end
if(~exist('curPosAcpc','var') | isempty(curPosAcpc)) curPosAcpc = dtiGet(handles, 'acpcpos'); end
skipAxes = isnan(curPosAcpc);
curPosAcpc(skipAxes) = 0;
if(~exist('bg','var') | isempty(bg))
    bg = handles.bg(dtiGet(handles,'curbgnum'));
    bg.acpcToImgXform = dtiGet(handles, 'acpc2imgxform');
    bg.mmPerVox = bg.mmPerVoxel;
end

% May want to interpolate here...
curPosImgInt = round(mrAnatXformCoords(bg.acpcToImgXform, curPosAcpc));

if(size(bg.img,4)==3)
	% Special case for vector images
	z = squeeze(bg.img(:,:,curPosImgInt(3),:));
	y = squeeze(bg.img(:,curPosImgInt(2),:,:));
	x = squeeze(bg.img(curPosImgInt(1),:,:,:));
else
    z = repmat(squeeze(bg.img(:,:,curPosImgInt(3))),[1,1,3]);
    y = repmat(squeeze(bg.img(:,curPosImgInt(2),:)),[1,1,3]);
    x = repmat(squeeze(bg.img(curPosImgInt(1),:,:)),[1,1,3]);
end

curPosImg = mrAnatXformCoords(bg.acpcToImgXform, curPosAcpc);
if(~exist('fname','var')) fname = ''; end

% Maybe the user wanted to set the parameters
roiBlurSize = 0;
fgBlurSize = 0;
if queryParameters
    prompt={'Upsample factor (1,2,4,8,16):','ROI blur size:','Fiber blur size:'};
    def={num2str(upSamp),num2str(roiBlurSize),num2str(fgBlurSize)};
    dlgTitle='Image processing parameters';
    lineNo=1;
    answer=inputdlg(prompt,dlgTitle,lineNo,def);
    if isempty(answer), disp('dtiSaveImageSlice cancelled.'); return;
    else
        upSamp = str2num(answer{1});
        roiBlurSize = str2num(answer{2});
        fgBlurSize = str2num(answer{3});
    end
end

% Figure out the file
persistent imPath;
if(nargout==0 & isempty(fname))
    if(isempty(imPath)) imPath = fullfile(handles.defaultPath,'slice'); end
    [p,f,e] = fileparts(imPath);
    imPath = fullfile(p,f);
    [f, p] = uiputfile({'*.png'}, ['Save current view...'], imPath);
    if(isnumeric(f)), disp('dtiSaveImageSlices cancelled.'); fname = 'show';
    else [junk,f,e] = fileparts(f); fname = fullfile(p, f); end
end
imPath = fname;

% Start image processing code
upSamp = 2^round(log2(upSamp));

% Resizing the image
imX = imageResize(x,upSamp);
imY = imageResize(y,upSamp);
imZ = imageResize(z,upSamp);

% Which order?  Put down the ROI, then overlay the fibers on top of that.
if(~isempty(overlayROIs))
    [imX,imY,imZ] = dtiImageOverlayROIs(overlayROIs,upSamp,roiBlurSize,imX,imY,imZ,curPosImg,bg.acpcToImgXform);
end
if(~isempty(overlayFGs))
    [imX,imY,imZ] = dtiImageOverlayFGs(overlayFGs,upSamp,fgBlurSize,imX,imY,imZ,curPosImg,bg.acpcToImgXform);
end

imX = permute(imX,[2,1,3]);
imY = permute(imY,[2,1,3]);
imZ = permute(imZ,[2,1,3]);

for(ii=1:3)
    imX(:,:,ii) = flipud(imX(:,:,ii));
    imY(:,:,ii) = flipud(imY(:,:,ii));
    imZ(:,:,ii) = flipud(imZ(:,:,ii));
end

% Add a scale bar
cmBarLen = round(1/bg.mmPerVox(1) * 10 * upSamp);
cbarVal = max(imX(:));
imX(end-5, end-cmBarLen-5:end-5, :) = cbarVal;
imX(end-4, end-cmBarLen-5:end-5, :) = cbarVal;
cmBarLen = round(1/bg.mmPerVox(2) * 10 * upSamp);
cbarVal = max(imY(:));
imY(end-5, end-cmBarLen-5:end-5, :) = cbarVal;
imY(end-4, end-cmBarLen-5:end-5, :) = cbarVal;
cmBarLen = round(1/bg.mmPerVox(3) * 10 * upSamp);
cbarVal = max(imZ(:));
imZ(end-5, end-cmBarLen-5:end-5, :) = cbarVal;
imZ(end-4, end-cmBarLen-5:end-5, :) = cbarVal;

if(~isempty(fname))
    if(strcmp(fname,'show'))
        figure; image(imX./max(imX(:))); axis image; truesize;
        figure; image(imY./max(imY(:))); axis image; truesize;
        figure; image(imZ./max(imZ(:))); axis image; truesize;
    else
        imwrite(imX, [fname '_X.png']);
        disp(['Wrote X slice to ' fname '_X.png']);
        imwrite(imY, [fname '_Y.png']);
        disp(['Wrote Y slice to ' fname '_Y.png']);
        imwrite(imZ, [fname '_Z.png']);
        disp(['Wrote Z slice to ' fname '_Z.png']);
    end
end
return;

%------------------------------------------------
function imOut = imageResize(im,m)
%
%   imOut = imageResize(im,m)
%
%Author: Wandell, Dougherty
%Purpose:
%   Resize an RGB image.  Uses the local upSample/upConv code.
%

sz = size(im);
m = round(log2(m));
for ii=1:3, imOut(:,:,ii) = upSample(im(:,:,ii),m); end

% Old, slow code
% sz = size(im);
% newSize = [sz(1), sz(2)]*m;
% imOut = zeros(newSize);
% for ii=1:3
%     imOut(:,:,ii) = imresize(im(:,:,ii),newSize,'bilinear');
% end
% return;

return;

%-----------------------------------------------------
function [imX,imY,imZ] = dtiImageOverlayFGs(overlayFGs,upSamp,blurSize,imX,imY,imZ,curPosImg,acpcToImgXform)
%
%        [imX,imY,imZ] = dtiImageOverlayFGs(handles,upSamp,blurSize,imX,imY,imZ,curPosImg,acpcToImgXform)
%
%Author: Dougherty, Wandell
%Purpose:
%   Add the fiber group positions to the image
%   It might be possible to upsample differently in the three dimensions.
%   But for now, we only upsample the whole data set with a single upSamp
%   value.
%

if ieNotDefined('upSamp'), upSamp = 1; end
if ieNotDefined('blurSize'), blurSize = 0; end

fpX = []; fpColorX = [];
fpY = []; fpColorY = [];
fpZ = []; fpColorZ = [];
for(grpNum=1:length(overlayFGs))
    if(overlayFGs(grpNum).visible)
        fp = horzcat(overlayFGs(grpNum).fibers{:});
        fp = mrAnatXformCoords(acpcToImgXform, fp');
        fiberColor = overlayFGs(grpNum).colorRgb/255;
        if(~isempty(fp))
            % Select those fiber points that are in this slice
            sfpX = fp(round(fp(:,1))==round(curPosImg(1)), [2,3])-1;
            sfpX = unique(round(sfpX.*upSamp),'rows');
            if(~isempty(sfpX))
                [mem,loc] = ismember(sfpX, fpX, 'rows');
                if(any(mem))
                    overlap = loc(loc>0);
                    for(ii=1:3)
                        fpColorX(overlap,ii) = fpColorX(overlap,ii)*0.5 + fiberColor(ii);
                    end
                end
                fpX = vertcat(fpX,sfpX(~mem,:));
                fpColorX = vertcat(fpColorX, repmat(fiberColor, sum(~mem), 1));
            end

            sfpY = fp(round(fp(:,2))==round(curPosImg(2)), [1,3])-1;
            sfpY = unique(round(sfpY.*upSamp),'rows');
            if(~isempty(sfpY))
                [mem,loc] = ismember(sfpY, fpY, 'rows');
                if(any(mem))
                    overlap = loc(loc>0);
                    for(ii=1:3)
                        fpColorY(overlap,ii) = fpColorY(overlap,ii)*0.5 + fiberColor(ii);
                    end
                end
                fpY = vertcat(fpY,sfpY(~mem,:));
                fpColorY = vertcat(fpColorY, repmat(fiberColor, sum(~mem), 1));
            end

            sfpZ = fp(round(fp(:,3))==round(curPosImg(3)), [1,2])-1;
            sfpZ = unique(round(sfpZ.*upSamp),'rows');
            if(~isempty(sfpZ))
                [mem,loc] = ismember(sfpZ, fpZ, 'rows');
                if(any(mem))
                    overlap = loc(loc>0);
                    for(ii=1:3)
                        fpColorZ(overlap,ii) = fpColorZ(overlap,ii)*0.5 + fiberColor(ii);
                    end
                end
                fpZ = vertcat(fpZ,sfpZ(~mem,:));
                fpColorZ = vertcat(fpColorZ, repmat(fiberColor, sum(~mem), 1));
            end
        end
    end
end

imX = dtiAddImageOverlay(imX, fpX, fpColorX, blurSize);
imY = dtiAddImageOverlay(imY, fpY, fpColorY, blurSize);
imZ = dtiAddImageOverlay(imZ, fpZ, fpColorZ, blurSize);

return;

%-----------------------------------------------------
function [imX,imY,imZ] = dtiImageOverlayROIs(overlayROIs,upSamp,blurSize,imX,imY,imZ,curPosImg,acpcToImgXform)
%
%   [imX,imY,imZ] = dtiImageOverlayROIs(overlayROIs,upSamp,blurSize,imX,imY,imZ,curPosImg,acpcToImgXform)
%
%Author: Dougherty, Wandell
%Purpose:
%   Add the ROI positions to the image
%

if ieNotDefined('upSamp'), error('upSamp required.'); end
if ieNotDefined('blurSize'), error('blurSize required.');; end

for(roiNum=1:length(overlayROIs))
    if(overlayROIs(roiNum).visible)
        roiPos = mrAnatXformCoords(acpcToImgXform, overlayROIs(roiNum).coords);
        roiColor =  dtiRoiGetColor(overlayROIs(roiNum));

        sroiPos = roiPos(round(roiPos(:,1))==round(curPosImg(1)),[2,3])';
        sroiPos = sroiPos.*upSamp;
        %sroiPosDilate = zeros(size(sroiPos,1)*upSamp,3);
        if(~isempty(sroiPos))
            % dilate up and to the left
            sroiPosDilate = [];
            for(ii=0:upSamp)
                for(jj=0:upSamp)
                    sroiPosDilate = [sroiPosDilate, [sroiPos(1,:)-ii; sroiPos(2,:)-jj]];
                end
            end
            imX = dtiAddImageOverlay(imX, round(sroiPosDilate)', roiColor, blurSize);
        end

        sroiPos = roiPos(round(roiPos(:,2))==round(curPosImg(2)),[1,3])';
        sroiPos = sroiPos.*upSamp;
        if(~isempty(sroiPos))
            sroiPosDilate = [];
            for(ii=0:upSamp)
                for(jj=0:upSamp)
                    sroiPosDilate = [sroiPosDilate, [sroiPos(1,:)-ii; sroiPos(2,:)-jj]];
                end
            end
            imY = dtiAddImageOverlay(imY, round(sroiPosDilate)', roiColor, blurSize);
        end

        sroiPos = roiPos(round(roiPos(:,3))==round(curPosImg(3)),[1,2])';
        sroiPos = sroiPos.*upSamp;
        if(~isempty(sroiPos))
            sroiPosDilate = [];
            for(ii=0:upSamp)
                for(jj=0:upSamp)
                    sroiPosDilate = [sroiPosDilate, [sroiPos(1,:)-ii; sroiPos(2,:)-jj]];
                end
            end
            imZ = dtiAddImageOverlay(imZ, round(sroiPosDilate)', roiColor, blurSize);
        end

    end
end

return;

%-----------------------------------------------------------
function newPts = dtiOverlayCoordXform(ax,upSamp,pts)
% Transform from acpc coordinates into image coordinates
newPts(2,:) = (-(ax(1,1) - pts(2,:)) + 1) *upSamp;
newPts(1,:) =  (ax(2,2) + (pts(1,:)  + 1)) *upSamp;
return;


text:
function BCType = CorrectBCTable_v2(EToV,VX,VY,BCType,fd,BCcode)

% function BCType = CorrectBCTable(EToV,BCType,fd,BCcode);
% Purpose: Setup BCType for boundary conditions in 2D
%
%    EToV   : Element-To-Vertice table
%    VX, VY : (x,y)-coordinates of mesh vertices
%    BCType : Table with types of faces for BC's
%    fd     : handle to distance function
%    BCcode : Integer for specific boundary type
%
% By Allan P. Engsig-Karup

Globals2D;

VNUM = [1 2;2 3;3 1]; % face orientations

pxc = 0.5*(VX(EToV)+VX(EToV(:,[2 3 1])));
pyc = 0.5*(VY(EToV)+VY(EToV(:,[2 3 1])));
dc  = abs(fd([pxc(:) pyc(:)])); % distances to boundaries from face centers
tol = 1e-4; % tolerance
idx = find(dc<tol);
BCType(idx) = BCcode;
return

repo: Aero-matlab


text:
function variance = gamma_variance ( a, b, c )

%*****************************************************************************80
%
%% GAMMA_VARIANCE returns the variance of the Gamma PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, C, the parameters of the PDF.
%    0.0 < B,
%    0.0 < C.
%
%    Output, real VARIANCE, the variance of the PDF.
%
  variance = b * b * c;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%% Copyright (C) 2016 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 3 of the License,
%% or (at your option) any later version.
%%
%% This software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod  @@sym hypot (@var{x}, @var{y})
%% @defmethodx @@sym hypot (@var{x}, @var{y}, @var{z}, @dots{})
%% Return hypoteneuse (distance) from symbolic expressions.
%%
%% Example of computing distance:
%% @example
%% @group
%% syms x y real
%% syms z
%% hypot (x, y, z)
%%   @result{} (sym)
%%          ________________
%%         ╱  2    2      2
%%       ╲╱  x  + y  + │z│
%% @end group
%% @end example
%%
%% Another example involving complex numbers:
%% @example
%% @group
%% hypot (sym([12 2]), [3+4i 1+2i])
%%   @result{} (sym) [13  3]  (1×2 matrix)
%% @end group
%% @end example
%%
%% @seealso{@@sym/atan2}
%% @end defmethod

function h = hypot(varargin)

  % two inputs:
  %h = sqrt(abs(sym(x)).^2 + abs(sym(y)).^2);

  two = sym(2);
  L = cellfun(@(x) power(abs(sym(x)), two), varargin, 'UniformOutput', false);
  s = L{1};
  for i=2:length(L);
    s = s + L{i};
  end
  h = sqrt(s);

end

%!assert (isequal (hypot (sym(3), 4), sym(5)))

%!test
%! % compare to @double (note Matlab hypot only takes 2 inputs)
%! A = hypot (hypot ([1 2 3], [4 5 6]), [7 8 9]);
%! B = double (hypot (sym([1 2 3]), [4 5 6], [7 8 9]));
%! assert (A, B, -eps)

%!test
%! % compare to @double, with complex
%! A = hypot ([1+2i 3+4i], [1 3+1i]);
%! B = double (hypot (sym([1+2i 3+4i]), [1 3+1i]));
%! assert (A, B, -eps)

%!test
%! % matrices
%! x = sym([1 -2; 0 3]);
%! y = sym([0 0; 8 4]);
%! A = hypot (x, y);
%! B = sym([1 2; 8 5]);
%! assert (isequal (A, B))

repo: octsympy