text:
function test_example_nirs_fingertapping

% MEM 4gb
% WALLTIME 00:10:00

%
%% Analyzing NIRS data recorded during unilateral finger- and foot-tapping
%
% This example script demonstrates the analysis of data that is shared by Sujin Bak, Jinwoo Park, Jaeyoung Shin, and Jichai Jeong: _[Open Access fNIRS Dataset for Classification of Unilateral Finger- and Foot-Tapping](https://doi.org/10.6084/m9.figshare.9783755.v2)_.
%
% The following links point to the shared data, to the PDF manuscript that explains the shared data, and to a GitHub repository from one of the authors that contains some example analyses on the data.
%
%* <https://doi.org/10.6084/m9.figshare.9783755.v2>
%* <https://doi.org/10.3390/electronics8121486>
%* <https://github.com/JaeyoungShin/fNIRS-dataset>
%
% The data includes the concentration changes of oxygenated/reduced hemoglobin ∆HbO/HbR, trigger information, and fNIRS channel information. The recordings were done for 30 subjects that were performing a finger- and foot-tapping task.
%
% The data were recorded by a three-wavelength continuous-time multi-channel fNIRS system (LIGHTNIRS, Shimadzu, Kyoto, Japan) consisting of eight light sources (Tx) and eight detectors (Rx). Four each of Tx and Rx were placed around C3 on the left hemisphere, and the rest were placed around C4 on the right hemisphere.
%
% The temporal structure of each trial is
%
%* instruction for 2s
%* fixation cross, this indicates that the task has to be executed, for 10s
%* "stop" for 2s
%* fixation cross indicating rest period, for 15-17 seconds
%
% The total duration of a trial, including the instruction, ranges from 2+10+2+15=29 to 2+10+2+17=31 seconds.
%
%% # Building a MATLAB analysis script
%

% for now, getting the data from the below URL and unzipping it, works (JM added 20211201)
t = tempdir;
filenames = unzip('https://figshare.com/ndownloader/files/18069143', t);

% filenames = {
%   'v2/fNIRS 01.mat'
%   'v2/fNIRS 02.mat'
%   'v2/fNIRS 03.mat'
%   'v2/fNIRS 04.mat'
%   'v2/fNIRS 05.mat'
%   'v2/fNIRS 06.mat'
%   'v2/fNIRS 07.mat'
%   'v2/fNIRS 08.mat'
%   'v2/fNIRS 09.mat'
%   'v2/fNIRS 10.mat'
%   'v2/fNIRS 11.mat'
%   'v2/fNIRS 12.mat'
%   'v2/fNIRS 13.mat'
%   'v2/fNIRS 14.mat'
%   'v2/fNIRS 15.mat'
%   'v2/fNIRS 16.mat'
%   'v2/fNIRS 17.mat'
%   'v2/fNIRS 18.mat'
%   'v2/fNIRS 19.mat'
%   'v2/fNIRS 20.mat'
%   'v2/fNIRS 21.mat'
%   'v2/fNIRS 22.mat'
%   'v2/fNIRS 23.mat'
%   'v2/fNIRS 24.mat'
%   'v2/fNIRS 25.mat'
%   'v2/fNIRS 26.mat'
%   'v2/fNIRS 27.mat'
%   'v2/fNIRS 28.mat'
%   'v2/fNIRS 29.mat'
%   'v2/fNIRS 30.mat'
%   };

% this allows switching easily between subjects, and eventually looping over all 30 subjects
filename = filenames{1};

%% # Exploring the MATLAB files that hold the NIRS data
%
% The data is shared by the authors in the form of MATLAB files. Each file contains a bunch of separate variables, you can load it like this to represent it as a structure:
%
nirs = load(filename);

% >> nirs
% nirs =
%   struct with fields:
%      mrk: [1x1 struct]
%      mnt: [1x1 struct]
%      nfo: [1x1 struct]
%      ch1: [30003x1 double]
%      ch2: [30003x1 double]
%      ch3: [30003x1 double]
%      ch4: [30003x1 double]
%      ch5: [30003x1 double]
%      ch6: [30003x1 double]
%      ch7: [30003x1 double]
%      ch8: [30003x1 double]
%      ch9: [30003x1 double]
%     ch10: [30003x1 double]
%     ch11: [30003x1 double]
%     ch12: [30003x1 double]
%     ch13: [30003x1 double]
%     ch14: [30003x1 double]
%     ch15: [30003x1 double]
%     ch16: [30003x1 double]
%     ch17: [30003x1 double]
%     ch18: [30003x1 double]
%     ch19: [30003x1 double]
%     ch20: [30003x1 double]
%     ch21: [30003x1 double]
%     ch22: [30003x1 double]
%     ch23: [30003x1 double]
%     ch24: [30003x1 double]
%     ch25: [30003x1 double]
%     ch26: [30003x1 double]
%     ch27: [30003x1 double]
%     ch28: [30003x1 double]
%     ch29: [30003x1 double]
%     ch30: [30003x1 double]
%     ch31: [30003x1 double]
%     ch32: [30003x1 double]
%     ch33: [30003x1 double]
%     ch34: [30003x1 double]
%     ch35: [30003x1 double]
%     ch36: [30003x1 double]
%     ch37: [30003x1 double]
%     ch38: [30003x1 double]
%     ch39: [30003x1 double]
%     ch40: [30003x1 double]
%      dat: [1x1 struct]

% Contrary to the description in the accompanying publication (see table 1 in the PDF manuscript) and in the GitHub repository, the MATLAB files do _not_ contain the variables |cntHb|, |clab|, etc. However, it is not so hard to make sense of the data: each channel is represented as a vector, there is header information, there is information about the optode montage, and there is information about the events.
%
% >> nirs.nfo
% ans =
%   struct with fields:
%             fs: 13.3333
%           clab: {1x40 cell}
%              T: 30003
%        nEpochs: 1
%         length: 2.2502e+03
%         format: 'DOUBLE'
%     resolution: [1x40 double]
%           file: 'D:\Experimental data\NIRS motor imagery\matdata\fNIRS 01'
%        nEvents: 75
%       nClasses: 3
%      className: {'RIGHT'  'LEFT'  'FOOT'}

%
% >> nirs.mnt
% ans =
%  struct with fields:
%               x: [40x1 double]
%               y: [40x1 double]
%          pos_3d: [3x40 double]
%            clab: {1x40 cell}
%             box: [2x41 double]
%          box_sz: [2x41 double]
%       scale_box: [2x1 double]
%    scale_box_sz: [2x1 double]

%
% >> nirs.mrk
% ans =
%  struct with fields:
%        event: [1x1 struct]
%         time: [1x75 double]
%            y: [3x75 double]
%    className: {'RIGHT'  'LEFT'  'FOOT'}

%% Convering the MATLAB structure to a FieldTrip raw data structure
%
% Since the data is not stored on disk in a [dataformat](/faq/dataformat/) that FieldTrip can directly read, we will circumvent the FieldTrip reading functions as outlined in [this frequenly asked question](/faq/how_can_i_import_my_own_dataformat/#circumvent-the-fieldtrip-reading-functions).
%
% We start with constructing a MATLAB data structure according to **[ft_datatype_raw](https://github.com/fieldtrip/fieldtrip/blob/release/utilities/ft_datatype_raw.m)**, as if it were produced by **[ft_preprocessing](https://github.com/fieldtrip/fieldtrip/blob/release/ft_preprocessing.m)**.
%
data_raw = [];
data_raw.label = nirs.nfo.clab(:);
data_raw.trial{1} = [
  % transpose each channel, we want the data matrix to be Nchans * Nsamples
  nirs.ch1'
  nirs.ch2'
  nirs.ch3'
  nirs.ch4'
  nirs.ch5'
  nirs.ch6'
  nirs.ch7'
  nirs.ch8'
  nirs.ch9'
  nirs.ch10'
  nirs.ch11'
  nirs.ch12'
  nirs.ch13'
  nirs.ch14'
  nirs.ch15'
  nirs.ch16'
  nirs.ch17'
  nirs.ch18'
  nirs.ch19'
  nirs.ch20'
  nirs.ch21'
  nirs.ch22'
  nirs.ch23'
  nirs.ch24'
  nirs.ch25'
  nirs.ch26'
  nirs.ch27'
  nirs.ch28'
  nirs.ch29'
  nirs.ch30'
  nirs.ch31'
  nirs.ch32'
  nirs.ch33'
  nirs.ch34'
  nirs.ch35'
  nirs.ch36'
  nirs.ch37'
  nirs.ch38'
  nirs.ch39'
  nirs.ch40'
  ];
data_raw.fsample = nirs.nfo.fs;
data_raw.time{1} = ((1:nirs.nfo.T)-1)/nirs.nfo.fs;

%% # Convert the optode montage to a layout for plotting
%
% In FieldTrip we can have a [detailled description](/faq/how_are_electrodes_magnetometers_or_gradiometers_described/#the-definition-of-nirs-sensors) of the NIRS optode placement and how optodes are combined to form channels. However, in this dataset this information is not complete and we cannot make a complete |opto| structure.
%
% This is how far we can get
%
opto = [];
opto.label = nirs.mnt.clab(:);
opto.chanpos = nirs.mnt.pos_3d'; % these are all nan

% But this is not a complete description of the channel and sensor information according to **[ft_dataype_sens](https://github.com/fieldtrip/fieldtrip/blob/release/utilities/ft_datatype_sens.m)**. However, a full sensor definitioon is also not required: the opical densities have already been converted in HbO and HbR prior to sharing, so we only care about channel positions for plotting.
%
% For the plotting of channel level data (see also [this tutorial](/tutorial/plotting/)) we need a 2D layout. That is explained in detail in [this tutorial](/tutorial/layout/). If the |opto| definition had included 2D or 3D channel positions, then we could have used **[ft_prepare_layout](https://github.com/fieldtrip/fieldtrip/blob/release/ft_prepare_layout.m)** but now we will manually construct the layout structure.
%
layout = [];
layout.label = nirs.mnt.clab(:);
layout.pos = [nirs.mnt.x(:) nirs.mnt.y(:)]; % these are all nan
layout.pos = nirs.mnt.box';
layout.width = nirs.mnt.box_sz(1,:)';
layout.height = nirs.mnt.box_sz(2,:)';
layout.outline = {};
layout.mask = {};

% I don't think that the 41st position refers to an actual channel, so let's remove that.
%
layout.pos(41, :) = [];
layout.width(41)  = [];
layout.height(41) = [];

% We can plot the layout using **[ft_plot_layout](https://github.com/fieldtrip/fieldtrip/blob/release/plotting/ft_plot_layout.m)** and inspect it in detail. That is something I have done multiple times, going back and forth in the script to ensure that my interpretation of the |mnt| field matches the data.
%
figure
ft_plot_layout(layout, 'label', true)

%
% This shows four groups, with two groups of oxy channels on the left of the figure (for the left and right hemisphere), and the corresponding deoxy channels on the right of the figure. This matches with the schepatic display in figure 1 of the PDF manuscript, and with the displayed data in figure 3 of the PDF manuscript.
%
% To improve the plotting of topographies, we can also make an outline and mask. The outline comproses extra lines that are added to the figure, often we include a circle (representing the head) with a triangle at the top (representing the nose). The mask comprises a series of polygons that are used in topographic interpolation and masking the interpolated data that is outside the head, but also by masking the interpolated data that is in between the grid on the left and right hemisphere. This will become more clear at the end of this example script.
%
pad = 0.2;
section = {1:10 11:20 21:30 31:40};
for i=1:4
  sel = section{i};
  % upper left
  corner1x = min(layout.pos(sel,1)-layout.width (sel)/2) - pad;
  corner1y = max(layout.pos(sel,2)+layout.height(sel)/2) + pad;

  % upper right
  corner2x = max(layout.pos(sel,1)+layout.width (sel)/2) + pad;
  corner2y = max(layout.pos(sel,2)+layout.height(sel)/2) + pad;

  % lower right
  corner3x = max(layout.pos(sel,1)+layout.width (sel)/2) + pad;
  corner3y = min(layout.pos(sel,2)-layout.height(sel)/2) - pad;

  % lower left
  corner4x = min(layout.pos(sel,1)-layout.width (sel)/2) - pad;
  corner4y = min(layout.pos(sel,2)-layout.height(sel)/2) - pad;

  layout.outline{i} = [
    corner1x corner1y
    corner2x corner2y
    corner3x corner3y
    corner4x corner4y
    corner1x corner1y % this closes the outline
    ];

  layout.mask{i} = [
    corner1x corner1y
    corner2x corner2y
    corner3x corner3y
    corner4x corner4y
    corner1x corner1y % this closes the outline
    ];
end

figure
ft_plot_layout(layout, 'label', true, 'outline', true, 'mask', true)

%
% The layout now also shows the outline, as thick black lines. If you disable the outline, you can also see the mask as dotted lines.
%
%% # Parsing the events and segmenting the data in trials
%
% In general in FieldTrip events have to be aligned with samples. The first sample of a recording on disk is referred to as sample 1.
%
% The |mrk| field contains times
%
% >> nirs.mrk.time
% ans =
%   Columns 1 through 9
%         1950    30975    60975    91950    ...

% These numbers are larger than the total number of samples, so they do not directly map onto samples. I started off assuming that the marker time is expressed in milliseconds, but that turned out not to be correct. Therefore I experimented with the following piece of code until I had something that made sense.
%
% This timestamp to sample ratio results in the events being nicely spaced over the whole experiment
%
ratio = 75;

% Note that `1/data.fsample = 0.075`, i.e., with the 13.333 Hz sampling rate each sample happens to be 75 ms. When working with unknown data it is not uncommon to have to deal with a puzzle like this, it reminds me of [Kabbalah](https://en.wikipedia.org/wiki/Kabbalah).
%
event = [];
for i=1:nirs.nfo.nEvents
  event(i).type = 'marker';
  event(i).value = nirs.mrk.className{nirs.mrk.event.desc(i)};
  event(i).sample = round(nirs.mrk.time(i)/ratio) + 1;
  event(i).duration = nan;
  event(i).offset = nan;
end

cfg = [];
cfg.event = event; % store the events in the configuration to have them displayed
cfg.viewmode = 'vertical';
cfg.preproc.demean = 'yes';
cfg.blocksize = nirs.nfo.length; % show the data over the whole length
ft_databrowser(cfg, data_raw);

%
%% # Apply a band-pass filter on the continuous data
%
cfg = [];
cfg.bpfilter = 'yes';
cfg.bpfreq = [0.02 0.8];
cfg.bpfilford = 2;
data_filt = ft_preprocessing(cfg, data_raw);

%% # Segment the continuous data into trials
%
% To segment the data we have to create a so-called "trial definition". This would normally be done using **[ft_definetrial](https://github.com/fieldtrip/fieldtrip/blob/release/ft_definetrials.m)** and using the events that were read from the data on disk (had it been a supported file format). The trial definition specifies the begin and end sample of each trial, and the offset, i.e., how many samples the time defined as t=0 is shifted relative to the data segment. Furthermore, it can contain additional columns with trial specific information, such as condition codes.
%
% The temporal structure of each trial is
%
%* instruction for 2s
%* fixation cross, this indicates that the task has to be executed, for 10s
%* "stop" for 2s
%* fixation cross indicating rest period, for 15-17 seconds
%
% The total duration of a trial, including the instruction, ranges from 2+10+2+15=29 to 2+10+2+17=31 seconds.
%
% We will use the standard FieldTrip |event| structure to make the trials, this is also what **[ft_read_event](https://github.com/fieldtrip/fieldtrip/blob/release/fileio/ft_read_event.m)** returns for supported data formats.
%
taskonset = [event.sample];  % let's assume that the marker is at task onset
begsample = taskonset - 2*data_raw.fsample;
endsample = taskonset + (10+2+15)*data_raw.fsample;
offset = - 2*data_raw.fsample * ones(size(begsample)); % each trial has the same 2 seconds offset

% The begin sample, end sample and offset should be integers:
%
begsample = round(begsample);
endsample = round(endsample);
offset = round(offset);

% and we map the right, left, and foot condition onto condition codes 1, 2, and 3.
%
[dummy, condition] = ismember({event.value}, {'RIGHT'  'LEFT'  'FOOT'});

trl = [begsample(:) endsample(:) offset(:) condition(:)];

cfg = [];
cfg.trl = trl;
data_segmented = ft_redefinetrial(cfg, data_filt);

%% # Visualize the segmented data
%
cfg = [];
cfg.event = event;
cfg.viewmode = 'vertical';
cfg.preproc.demean = 'no'; % it is filtered by now
ft_databrowser(cfg, data_segmented);

%
%% # Average the segmented data, conditional on each condition
%
cfg = [];
cfg.trials = data_segmented.trialinfo==1;
timelockRight = ft_timelockanalysis(cfg, data_segmented);
cfg.trials = data_segmented.trialinfo==2;
timelockLeft = ft_timelockanalysis(cfg, data_segmented);
cfg.trials = data_segmented.trialinfo==3;
timelockFoot = ft_timelockanalysis(cfg, data_segmented);

%% # Plot the averaged responses
%
%
cfg = [];
cfg.layout = layout;
cfg.showoutline = 'yes';
cfg.baseline = [-2 0];
cfg.showlabels = 'yes';
cfg.interactive = 'no'; % use MATLAB for zooming in
ft_multiplotER(cfg, timelockRight, timelockLeft, timelockFoot);

%
% Although some of the more noisy channels dominate the figure due to the automatic vertical scaling, if you zoom in and pay attention to |Ch4oxy|, |Ch5oxy|, and |Ch6oxy| over the left hemisphere, and |Ch14oxy| and |Ch16oxy| over the right hemisphere, you can see that the left hemisphere shows more activity during the right fingertapping task (blue) and the right hemisphere shows more activity following the left fingertapping task (red).
%
% This [example script](/example/nirs_layout/) explains in more detail how for a NIRS dataset the oxy and deoxy channels can either be plotted on top of each other (which is spatially consistent with how the data is recorded) or side-by-side. Plotting them side-by-side makes the results easier to interpret, especially if you have multiple responses (like we have here for left, right and foot conditions). For plotting topographic distributions you also have to ensure that your 2D channel layout does not have overlapping channels that have bvery different numbers.
%
%% # Plot the topography of the averaged responses
%
% This step benefits from the layout having an outline and especially a mask, otherwise the interpolation will cause the whole image to be color coded, including the pieces in between the hemispheres and in between the oxy (on the left) and deoxy channels (on the right).
%
cfg = [];
cfg.layout = layout;
cfg.interactive = 'yes';
cfg.xlim = [7 15];
cfg.ylim = [-0.004 0.004];
figure; ft_topoplotER(cfg, timelockRight);
figure; ft_topoplotER(cfg, timelockLeft);
figure; ft_topoplotER(cfg, timelockFoot);

%
%
% From the topographic arrangements it is clear that the left hemisphere responds with an increase in HbO to the right fingertapping task and the right hemisphere to the left fingertapping task; the HbR (on the right side of the figure) shows the opposite pattern. Also interesting is that both left and right hemisphere show a decrease in HbO and an increase in HbR during the tapping of the foot.

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function test_suite = testBadSinTest
initTestSuite;

function testSinPi
% Example of a failing test case.  The test writer should have used
% assertAlmostEqual here.
assertEqual(sin(pi), 0);

author: ricedsp


text:
function savestack(img, fname, varargin)
% Save a 3D image into a file or a series of files.
%
%   savestack(IMG, FNAME)
%   Saves the 3D image given in IMG into the file(s) given by FNAME.
%   IMG should be either a 3D array of grayscale or intensity values, or a
%   4D array of color values (coded as uint8).
%   FNAME is either a single file name, or a file name pattern that is used
%   to save the stack into a file series.
%
%   If file name contains '??', '##' or '%0xd' (x being an integer), then
%   the image is saved into a series of files, with increasing indices.
%
%   savestack(IMG, MAP, FNAME)
%   Saves a grayscale image together with a colormap. IMG should be a uint8
%   or uint16 image, and MAP should be a valid Matlab colormap (N-by-3
%   array coded as double, between 0 and 1)
%
%   savestack(..., FNAME, INDICES)
%   Also specifies number associated with each slice. Can be useful when
%   exporting a portion of a 3D stack, after a crop or a downsampling in Z
%   direction.
%
%   savestack(..., OPTIONS)
%   Uses options to write each slice. See imwrite to details.
%
%   savestack(..., VERBOSITY)
%   Also specifies verbosity. VERBOSITY can be either 'verbose' or 'quiet'.
%
%
%   Examples:
%     savestack(img, 'imgBundle.tif');        % save as stack
%     savestack(img, 'imgBundle???.tif');     % save as image series
%     savestack(img, 'imgBundle###.tif');     % save as image series
%     savestack(img, 'imgBundle%03d.tif');    % save as image series
%
%
%   See also:
%     readstack, imwrite
%

%   ---------
%   author: David Legland, david.legland@inra.fr
%   INRA - Cepia Software Platform
%   created the 10/09/2003.
%   http://www.pfl-cepia.inra.fr/index.php?page=slicer

%   HISTORY
%   17/02/2004 change indices of slices, to start at -00 and not at -01
%   18/10/2004 adapt to save in other formats than tif, and keeping
%       possibility to specify options
%   02/01/2005 correct bug: if save a bundle TIF on existing file,
%       remove previous file.
%   01/06/2006 allow wildcard '###', ensure more portability with
%       windows, add verbosity options
%   28/11/2008 update doc, remove unused variable 'slice'
%   30/10/2019 uses TIFF lib to accelerate write time

%% Default values

% default verbosity is true
verbose = true;

%% Process input arguments

% check number of input arguments: at least two, at most five
if nargin < 1 || nargin > 5
    error('savestack requires at least two arguments, and at most 5');
end

% image size
dim = size(img);

% check if colormap is specified
map = [];
if isnumeric(fname)
    map = fname;
    if nargin < 3
        error('Should specify file name as third input');
    end
    fname = varargin{1};
    varargin(1) = [];
end

% check for verbosity options
if ~isempty(varargin)
    var = varargin{end};
    if ischar(var)
       if strcmp(var, 'quiet') || strcmp(var, 'silent')
            verbose = false;
            varargin = varargin{1:end-1};
       end
    end
end

% check image number of dimensions
nd = length(dim);
if nd == 3
    isGrayscale = true;
    sliceIndices = 0:dim(3)-1;

elseif nd == 4
    if ~isempty(map)
        error('Can not save a 4D image array with map');
    end
    isGrayscale = false;
    sliceIndices = 0:dim(4)-1;
end

% determine optional slice indices (default: 0 -> dim(3))
if ~isempty(varargin) && isnumeric(varargin{1})
    disp('use specified slice indices');
    sliceIndices = varargin{1};
    varargin(1) = [];
end

%% Compute file name pattern

% replace wildcard '??' or '???' by wildcard '%02d' or '%03d'
pos = strfind(fname, '?');
npos = length(pos);
if npos > 0
    fname = strrep(fname, repmat('?', [1 npos]), ['%0' num2str(npos) 'd']);
end

% replace wildcard '##' or '###' by wildcard '%02d' or '%03d'
pos = strfind(fname, '#');
npos = length(pos);
if npos > 0
    fname = strrep(fname, repmat('#', [1 npos]), ['%0' num2str(npos) 'd']);
end

% binary flag indicating if the stack shouldbe save as a single file, or as
% a collection of 2D slices.
saveAsStack = ~contains(fname, '%0');

%% Save image file(s)

if saveAsStack
    % save one file containing all slices of image
    if verbose
        disp('save a stack');
    end

    % setup TIFF tags shared by all images
    tagStruct.ImageWidth = size(img, 2);
    tagStruct.ImageLength = size(img, 1);
    if isa(img, 'uint8')
        tagStruct.BitsPerSample = 8;
    elseif isa(img, 'uint16')
        tagStruct.BitsPerSample = 16;
    else
        error(['Can not manage data with class: ' class(img)]);
    end
    if isGrayscale
        tagStruct.SamplesPerPixel = 1;
        tagStruct.Photometric = Tiff.Photometric.MinIsBlack;

        % eventually update Tag structure to understand color map
        if ~isempty(map)
            tagStruct.Photometric = Tiff.Photometric.Palette;
            tagStruct.ColorMap = map;
        end
    else
        tagStruct.SamplesPerPixel = 3;
        tagStruct.Photometric = Tiff.Photometric.RGB;
    end
    tagStruct.Compression = Tiff.Compression.PackBits;
    tagStruct.PlanarConfiguration = Tiff.PlanarConfiguration.Chunky;
    tagStruct.Software = 'MATLAB_MatImage';

    % create a TIFF file, and setup necessary tags
    t = Tiff(fname, 'w');
    setTag(t, tagStruct);

    if isGrayscale
        % save a grayscale stack

        % write first slice
        write(t, img(:,:,1));

        % append other slices
        for i = 2:dim(3)
            writeDirectory(t);
            setTag(t, tagStruct);
            write(t, img(:,:,i));
        end

    else
        % save a color stack

        % write first slice
        write(t, img(:,:,:,1));

        % append other slices
        for i = 2:dim(4)
            writeDirectory(t);
            setTag(t, tagStruct);
            write(t, img(:,:,:,i));
        end
    end

    close(t);

else
    % save a series of file, one file per slice of images
    if verbose
        disp('save slices');
    end

    % indentify the position of output file name that will contain indices
    pos = strfind(fname, '%0');
    pos = pos(end);

    % extract different parts of the file name
    basename = fname(1:pos-1);
    endname = fname(pos+4:end);

    % string format to compute image name
    % -> basename + indSlice + endName
    format = ['%s%0' fname(pos+2) 'd%s'];

    % save each slice of the image (gray-scale->3D or color->4D)
    if isGrayscale
        if isempty(map)
            % save grayscale slice images
            for i = 1:dim(3)
                index = sliceIndices(i);
                fileName = sprintf(format, basename, index, endname);
                imwrite(img(:,:,i), fileName, ...
                    'WriteMode', 'overwrite', varargin{:});
            end

        else
            % save grayscale slice images with colormap
            for i = 1:dim(3)
                index = sliceIndices(i);
                fileName = sprintf(format, basename, index, endname);
                imwrite(img(:,:,i), map, fileName, ...
                    'WriteMode', 'overwrite', varargin{:});
            end

        end
    else
        % save color slice images
        for i = 1:dim(4)
            index = sliceIndices(i);
            fileName = sprintf(format, basename, index, endname);
            imwrite(img(:,:,:,i), fileName, ...
                'WriteMode', 'overwrite', varargin{:});
        end
    end
end

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function [datafit,Amps,freqs,Fval,sig]=fitlinesc(data,params,p,plt,f0,fscanbw)
% fits significant sine waves to data (continuous data).
%
% Usage: [datafit,Amps,freqs,Fval,sig]=fitlinesc(data,params,p,plt,f0)
%
%  Inputs:
% Note that units of Fs, fpass have to be consistent.
%       data        (data in [N,C] i.e. time x channels/trials or a single
%       vector) - required.
%       params      structure containing parameters - params has the
%       following fields: tapers, Fs, fpass, pad
%           tapers : precalculated tapers from dpss or in the one of the following
%                    forms:
%                   (1) A numeric vector [TW K] where TW is the
%                       time-bandwidth product and K is the number of
%                       tapers to be used (less than or equal to
%                       2TW-1).
%                   (2) A numeric vector [W T p] where W is the
%                       bandwidth, T is the duration of the data and p
%                       is an integer such that 2TW-p tapers are used. In
%                       this form there is no default i.e. to specify
%                       the bandwidth, you have to specify T and p as
%                       well. Note that the units of W and T have to be
%                       consistent: if W is in Hz, T must be in seconds
%                       and vice versa. Note that these units must also
%                       be consistent with the units of params.Fs: W can
%                       be in Hz if and only if params.Fs is in Hz.
%                       The default is to use form 1 with TW=3 and K=5
%
%	        Fs 	        (sampling frequency) -- optional. Defaults to 1.
%               fpass       (frequency band to be used in the calculation in the form
%                                   [fmin fmax])- optional.
%                                   Default all frequencies between 0 and Fs/2
%	        pad		    (padding factor for the FFT) - optional (can take values -1,0,1,2...).
%                    -1 corresponds to no padding, 0 corresponds to padding
%                    to the next highest power of 2 etc.
%			      	 e.g. For N = 500, if PAD = -1, we do not pad; if PAD = 0, we pad the FFT
%			      	 to 512 points, if pad=1, we pad to 1024 points etc.
%			      	 Defaults to 0.
%	    p		    (P-value to calculate error bars for) - optional.
%                           Defaults to 0.05/N where N is data length.
%       plt         (y/n for plot and no plot respectively) - plots the
%       Fratio at all frequencies if y
%       f0          frequencies at which you want to remove the
%                   lines - if unspecified the program
%                   will compute the significant lines
%       fscanbw     bandwidth centered on ea. f0 to scan for significant
%                   lines (TM)
%
%  Outputs:
%       datafit        (linear superposition of fitted sine waves)
%       Amps           (amplitudes at significant frequencies)
%       freqs          (significant frequencies)
%       Fval           (Fstatistic at all frequencies)
%       sig            (significance level for F distribution p value of p)
data=change_row_to_column(data);
[N,C]=size(data);
if nargin < 2 || isempty(params); params=[]; end;
[tapers,pad,Fs,fpass,err,trialave,params]=getparams(params);
clear pad fpass err trialave;
if nargin < 3 || isempty(p);p=0.05/N;end;
if nargin < 4 || isempty(plt); plt='n'; end;
if nargin < 5; f0=[]; end;
if nargin < 6; fscanbw = []; end

params.tapers=dpsschk(tapers,N,Fs); % calculate the tapers
[Fval,A,f,sig] = ftestc(data,params,p,plt);
if isempty(f0) && isempty(fscanbw)
    fmax=chron_findpeaks(Fval,sig);
    freqs=cell(1,C);
    Amps=cell(1,C);
    datafit=data;
    for ch=1:C;
        fsig=f(fmax(ch).loc);
        freqs{ch}=fsig;
        Amps{ch}=A(fmax(ch).loc,ch);
        Nf=length(fsig);
        %       fprintf('The significant lines for channel %d and the amplitudes are \n',ch);
        %        for nf=1:Nf;
        %            fprintf('%12.8f\n',fsig(nf));
        %            fprintf('%12.8f\n',real(A(fmax(ch).loc(nf),ch)));
        %            fprintf('%12.8f\n',imag(A(fmax(ch).loc(nf),ch)));
        %            fprintf('\n');
        %        end;
        datafit(:,ch)=exp(i*2*pi*(0:N-1)'*fsig/Fs)*A(fmax(ch).loc,ch)+exp(-i*2*pi*(0:N-1)'*fsig/Fs)*conj(A(fmax(ch).loc,ch));
    end;
else

    fsig = [];

    if ~isempty(fscanbw)

        % for each line f0(n), scan f0+-BW/2 for largest significant peak of Fval

        k = 1;
        for n=1:length(f0)
            % extract scan range around f0 ( f0 +- fscanbw/2 )
            [dummy ridx(1)] = min(abs(f-(f0(n)-fscanbw/2)));
            [dummy ridx(2)] = min(abs(f-(f0(n)+fscanbw/2)));

            Fvalscan = Fval(ridx(1):ridx(2));
            Fvalscan(Fvalscan<sig)=0;
            if any(Fvalscan)
                % if there's a significant line, pull the max one
                [dummy rmaxidx] = max(Fvalscan);
                indx(k) = ridx(1)+rmaxidx-1;
                fsig(k) = f(indx(k));
                k = k+1;
            end

        end
    else

        % remove exact lines

        indx = zeros( size(f0) );
        for n=1:length(f0);
            [fsig,indx(n)]=min(abs(f-f0(n)));
        end;
        fsig=f(indx);
    end

    if isempty(fsig)
        % no significant lines
        datafit = zeros(N,C);
        Amps  = cell(1,C);
        freqs = cell(1,C);
    else
        for ch=1:C;
            freqs{ch}=fsig;
            Amps{ch}=A(indx,ch);
            Nf=length(fsig);
            %        fprintf('For channel %d the amplitudes and the Fstatistic at f=%f are \n',ch,f0);
            %        fprintf('Fstatistic = %12.8f Fthreshold = %12.8f\n',Fval(indx),sig);
            %        fprintf('Real part of amplitude = %12.8f\n',real(A(indx,ch)));
            %        fprintf('Imaginary part of amplitude = %12.8f\n',imag(A(indx,ch)));
            datafit(:,ch)=exp(i*2*pi*(0:N-1)'*fsig/Fs)*A(indx,ch)+exp(-i*2*pi*(0:N-1)'*fsig/Fs)*conj(A(indx,ch));
        end;
    end
end;

author: VisLab

repo: EEG-Clean-Tools


text:
function [V2D] = Vandermonde2D(N,r,s)

% function [V2D] = Vandermonde2D(N, r, s);
% Purpose : Initialize the 2D Vandermonde Matrix,  V_{ij} = phi_j(r_i, s_i);

V2D = zeros(length(r),(N+1)*(N+2)/2);

% Transfer to (a,b) coordinates
[a, b] = rstoab(r, s);

% build the Vandermonde matrix
sk = 1;
for i=0:N
  for j=0:N - i
    V2D(:,sk) = Simplex2DP(a,b,i,j);
    sk = sk+1;
  end
end
return;

repo: Aero-matlab


text:
classdef factorization_dense_lu < factorization_generic
%FACTORIZATION_DENSE_LU dense LU factorization: p*A = L*U

% Copyright 2009, Timothy A. Davis, University of Florida

    methods

        function F = factorization_dense_lu (A)
            %FACTORIZATION_DENSE_LU dense LU: p*A = L*U
            n = size (A,2) ;
            [F.L U p] = lu (A, 'vector') ;
            assert (nnz (diag (U)) == n, 'Matrix is rank deficient.') ;
            F.p = sparse (1:n, p, 1) ;
            F.U = U ;
            F.A = A ;
        end

        function disp (F)
            %DISP displays a dense LU factorization
            fprintf ('  dense LU factorization: p*A = L*U\n') ;
            fprintf ('  A:\n') ; disp (F.A) ;
            fprintf ('  L:\n') ; disp (F.L) ;
            fprintf ('  U:\n') ; disp (F.U) ;
            fprintf ('  p:\n') ; disp (F.p) ;
            fprintf ('  is_inverse: %d\n', F.is_inverse) ;
        end

        function x = mldivide_subclass (F,b)
            %MLDIVIDE_SUBCLASS x=A\b using dense LU
            % x = U \ (L \ (p*b)) ;
            if (issparse (b))
                b = full (b) ;
            end
            opL.LT = true ;
            opU.UT = true ;
            x = linsolve (F.U, linsolve (F.L, F.p*b, opL), opU) ;
        end

        function x = mrdivide_subclass (b,F)
            %MRDIVIDE_SUBCLASS x = b/A using dense LU
            % x = (P' * (L' \ (U' \ b')))'
            bT = b' ;
            if (issparse (bT))
                bT = full (bT) ;
            end
            opUT.UT = true ;
            opUT.TRANSA = true ;
            opLT.LT = true ;
            opLT.TRANSA = true ;
            x = (F.p' * linsolve (F.L, linsolve (F.U, bT, opUT), opLT))' ;
        end
    end
end

repo: REST


text:
function L = obslikepoisson(X,hmm)
%
% Evaluate likelihood of data given observation model, for one continuous trial
%
% INPUT
% X          N by ndim data matrix
% hmm        hmm data structure
%
% OUTPUT
% B          Likelihood of N data points
%
% Author: Cam Higgins, OHBA, University of Oxford

K = hmm.K;
[T,ndim] = size(X);

L = zeros(T,K);

for k=1:K
    constterm = -gammaln(X + 1);
    % note the expectation of log(lambda) is -log(lambda_b) + psigamma(lambda_a)
    num = (X.*(repmat(psi(hmm.state(k).W.W_shape),T,1) - log(hmm.state(k).W.W_rate))) - ...
        repmat(hmm.state(k).W.W_mean,T,1);
    L(:,k) = sum(num + constterm,2);
end
L = exp(L);
end

author: OHBA-analysis

repo: HMM-MAR


text:
%% Factorial
function [f] = find_factorial(n)
% calculate the factorial of a positive integer n
% factorial(n)can be used directly as a default function of Matlab

integerTest= ~mod(n,1); %it returns 0 if value is not an integer.
if integerTest== 0 || n < 0;   % checking n to be positive & integer
    disp('Error! your number muss be positive and integer');
else
    f = 1;
    for i = 1:n
        f = f*i;
    end
    disp(['factorial of ',num2str(n),' is: ',num2str(f)]);
end
end

author: TheAlgorithms

repo: MATLAB-Octave


text:
function marginal = marginal_nodes(engine, query)
% MARGINAL_NODES Compute the marginal on the specified query nodes (gaussian)
% marginal = marginal_nodes(engine, query)

% Compute sum_{Hsum} Pr(Hkeep, Hsum | o)
H = engine.hnodes;
bnet = bnet_from_engine(engine);
ns = bnet.node_sizes;
Hkeep = myintersect(H, query);
Hsum = mysetdiff(H, Hkeep);

[marginal.mu, marginal.Sigma] = marginalize_gaussian(engine.Hmu, engine.HSigma, Hkeep, Hsum, ns);
marginal.domain = query;
marginal.T = 1;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [tri, inds] = triangulateFaces(faces)
%TRIANGULATEFACES Convert face array to an array of triangular faces
%
%   TRI = triangulateFaces(FACES)
%   Returns a 3-columns array of indices, based on the data stored in the
%   argument FACES:
%   - if FACES is a N-by-3 array, returns the same array
%   - if FACES is a N-by-4 array, returns an array with 2*N rows and 3
%       columns, splitting each square into 2 triangles (uses first and
%       third vertex of each square as diagonal).
%   - if FACES is a cell array, split each face into a set of triangles,
%       and returns the union of all triangles. Faces are assumed to be
%       convex.
%
%   [TRI INDS] = triangulateFaces(FACES)
%   Also returns original face index of each new triangular face. INDS has
%   the same number of rows as TRI, and has values between 1 and the
%   number of rows of the original FACES array.
%
%
%   Example
%     % create a basic shape
%     [n e f] = createCubeOctahedron;
%     % draw with plain faces
%     figure;
%     drawMesh(n, f);
%     % draw as a triangulation
%     tri = triangulateFaces(f);
%     figure;
%     patch('vertices', n, 'faces', tri, 'facecolor', 'r');
%
%   See also
%   meshes3d, drawMesh, mergeCoplanarFaces
%

% ------
% Author: David Legland
% e-mail: david.legland@nantes.inra.fr
% Created: 2008-09-08,    using Matlab 7.4.0.287 (R2007a)
% Copyright 2008 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas.

%% Tri mesh case: return original set of faces

if isnumeric(faces) && size(faces, 2) == 3
    tri = faces;
    if nargout > 1
        inds = (1:size(faces, 1))';
    end
    return;
end

%% Square faces: split each square into 2 triangles

if isnumeric(faces) && size(faces, 2) == 4
    nf = size(faces, 1);
    tri = zeros(nf * 2, 3);
    tri(1:2:end, :) = faces(:, [1 2 3]);
    tri(2:2:end, :) = faces(:, [1 3 4]);

    if nargout > 1
        inds = kron(1:size(faces, 1), ones(1,2))';
    end

    return;
end

%% Pentagonal faces (for dodecahedron...): split into 3 triangles

if isnumeric(faces) && size(faces, 2) == 5
    nf = size(faces, 1);
    tri = zeros(nf * 3, 3);
    tri(1:3:end, :) = faces(:, [1 2 3]);
    tri(2:3:end, :) = faces(:, [1 3 4]);
    tri(3:3:end, :) = faces(:, [1 4 5]);

    if nargout > 1
        inds = kron(1:size(faces, 1), ones(1,2))';
    end

    return;
end

%% Faces as cell array

% number of faces
nf  = length(faces);

% compute total number of triangles
ni = zeros(nf, 1);
for i = 1:nf
    % as many triangles as the number of vertices minus 1
    ni(i) = length(faces{i}) - 2;
end
nt = sum(ni);

% allocate memory for triangle array
tri = zeros(nt, 3);
inds = zeros(nt, 1);

% convert faces to triangles
t = 1;
for i = 1:nf
    face = faces{i};
    nv = length(face);
    v0 = face(1);
    for j = 3:nv
        tri(t, :) = [v0 face(j-1) face(j)];
        inds(t) = i;
        t = t + 1;
    end
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function [ce_flux_p, ce_flux_ps, ce_flux_s, ce_flux_sn, ce_flux_n] = interpolateElectrolyteConcetrationFluxes(ce,param)
%	interpolateElectrolyteConcetrationFluxes interpolates the electrolyte concentration flux at the edges of control volumes using harmonic mean.

%   This file is part of the LIONSIMBA Toolbox
%
%	Official web-site: 	http://sisdin.unipv.it/labsisdin/lionsimba.php
% 	Official GitHUB: 	https://github.com/lionsimbatoolbox/LIONSIMBA
%
%   LIONSIMBA: A Matlab framework based on a finite volume model suitable for Li-ion battery design, simulation, and control
%   Copyright (C) 2016-2018 :Marcello Torchio, Lalo Magni, Davide Raimondo,
%                            University of Pavia, 27100, Pavia, Italy
%                            Bhushan Gopaluni, Univ. of British Columbia,
%                            Vancouver, BC V6T 1Z3, Canada
%                            Richard D. Braatz,
%                            Massachusetts Institute of Technology,
%                            Cambridge, Massachusetts 02142, USA
%
%   Main code contributors to LIONSIMBA 2.0:
%                           Ian Campbell, Krishnakumar Gopalakrishnan,
%                           Imperial college London, London, UK
%
%   LIONSIMBA is a free Matlab-based software distributed with an MIT
%   license.

% Fluxes within the positive electrode
ce_flux_p = (ce(2:param.Np)-ce(1:param.Np-1))/(param.deltax_p*param.len_p);

% Fluxes at the separator-positive interface
ce_flux_ps = (ce(param.Np+1)-ce(param.Np)) / ((param.deltax_p*param.len_p/2+param.deltax_s*param.len_s/2));

% Fluxes within the separator
ce_flux_s = (ce(param.Np+2:param.Np+param.Ns)-ce(param.Np+1:param.Np+param.Ns-1))/(param.deltax_s*param.len_s);

% Fluxes at the separator-negative interface
ce_flux_sn = (ce(param.Np+param.Ns+1)-ce(param.Np+param.Ns)) / ((param.deltax_n*param.len_n/2+param.deltax_s*param.len_s/2));

% Fluxes within the negative electrode
ce_flux_n = (ce(param.Np+param.Ns+2:end)-ce(param.Np+param.Ns+1:end-1))/(param.deltax_n*param.len_n);

end

author: lionsimbatoolbox

repo: LIONSIMBA


text:
%STYLIZATION  Stylization filter
%
%     dst = cv.stylization(src)
%     dst = cv.stylization(src, 'OptionName',optionValue, ...)
%
% ## Input
% * __src__ Input 8-bit 3-channel image.
%
% ## Output
% * __dst__ Output image with the same size and type as `src`.
%
% ## Options
% * __SigmaS__ Range between 0 to 200. default 60
% * __SigmaR__ Range between 0 to 1. default 0.45
% * __FlipChannels__ whether to flip the order of color channels in input
%   `src` and output `dst`, between MATLAB's RGB order and OpenCV's BGR
%   (input: RGB->BGR, output: BGR->RGB). default false
%
% Stylization aims to produce digital imagery with a wide variety of effects
% not focused on photorealism. Edge-aware filters are ideal for stylization,
% as they can abstract regions of low contrast while preserving, or enhancing,
% high-contrast features.
%
% See also: cv.pencilSketch
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function x = f1_abscissas_ab ( a, b, n )

%*****************************************************************************80
%
%% F1_ABSCISSAS_AB computes Fejer type 1 abscissas for the interval [A,B].
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 July 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Philip Davis, Philip Rabinowitz,
%    Methods of Numerical Integration,
%    Second Edition,
%    Dover, 2007,
%    ISBN: 0486453391,
%    LC: QA299.3.D28.
%
%    Walter Gautschi,
%    Numerical Quadrature in the Presence of a Singularity,
%    SIAM Journal on Numerical Analysis,
%    Volume 4, Number 3, 1967, pages 357-362.
%
%    Joerg Waldvogel,
%    Fast Construction of the Fejer and Clenshaw-Curtis Quadrature Rules,
%    BIT Numerical Mathematics,
%    Volume 43, Number 1, 2003, pages 1-18.
%
%  Parameters:
%
%    Input, real A, B, the endpoints of the interval.
%
%    Input, integer N, the order of the rule.
%
%    Output, real X(N), the abscissas.
%
  if ( n == 1 )
    x(1) = 0.5 * ( b + a );
    return
  end

  for i = 1 : n
    theta(i) = ( 2 * n - 2 * i + 1 ) * pi ...
             / ( 2 * n             );
  end

  x(1:n) = 0.5 * ( ( b + a ) + ( b - a ) * cos ( theta(1:n) ) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function makeGif( M, fileName, prm )
% Writes a matlab movie to an animated GIF.
%
% USAGE
%  makeGif( M, fileName, prm )
%
% INPUTS
%  M          - Matlab movie
%  fileName   - file name of the output gif file
%  prm
%   .fps        - number of frames per second
%   .nColor     - number of indexes colors
%   .scale      - scale to resize the frames at
%   .loop       - number of times to repeat the movie (Inf possible)
%
% OUTPUTS
%
% EXAMPLE
%  load( 'images.mat' );
%  M = playMovie( video, [], 1 );
%  makeGif( M, 'mouse.gif', struct('scale',0.5) );
%
% See also PLAYMOVIE, MONTAGE2
%
% Piotr's Computer Vision Matlab Toolbox      Version 2.0
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

dfs = {'fps',30, 'nColor',256, 'scale',1, 'loop', 0};
prm = getPrmDflt( prm, dfs );
fps=prm.fps;  nColor=prm.nColor; scale=prm.scale; loop=prm.loop;

if scale~=1
  for i=1:length(M); M(i).cdata=imresize(M(i).cdata,scale, 'bicubic'); end
end

if ndims(M(1).cdata)==3
  for i=1:length(M)
    [M(i).cdata,M(i).colormap] = rgb2ind(M(i).cdata, nColor);
  end
else
  for i=1:length(M)
    M(i).colormap = repmat(0:1/256:1,[3 1])';
  end
end

imwrite(M(1).cdata,M(1).colormap,fileName,'gif','LoopCount',loop);
for i=2:length(M)
  imwrite(M(i).cdata,M(i).colormap,fileName,'gif','DelayTime',1/fps,...
    'WriteMode','append');
end

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
% Coupled nonlinear PDE's
% Buckley Leverett equation
% dependent variables: pressure and water saturation
% Prepared for educational purposes by ** AAE **
% Works fine; slow!
% Written by Ali A. Eftekhari
% Last checked: June 2021
clc
%% define the geometry
well_radius = 0.1; % [m]
x_well_radius = 10*well_radius; % area close to the well
Nx_mid = 50; % number of cells in x direction
Ny_mid = 50; % number of cells in y direction
Lx = 50; % [m]
Ly = 50; % [m]
x = [0:well_radius:x_well_radius linspace(x_well_radius+well_radius, Lx-x_well_radius, Nx_mid) Lx-x_well_radius+well_radius:well_radius:Lx];
y = [0:well_radius:x_well_radius linspace(x_well_radius+well_radius, Ly-x_well_radius, Ny_mid) Ly-x_well_radius+well_radius:well_radius:Ly];
m = createMesh2D(x, y); % creates a 2D mesh
Nx = length(x)-1;
Ny = length(y)-1;
%% define the physical parametrs
krw0 = 1.0;
kro0 = 0.76;
nw = 2.4;
no = 2.0;
sor=0.12;
swc=0.09;
sws=@(sw)((sw>swc).*(sw<1-sor).*(sw-swc)/(1-sor-swc)+(sw>=1-sor).*ones(size(sw)));
kro=@(sw)((sw>=swc).*kro0.*(1-sws(sw)).^no+(sw<swc).*(1+(kro0-1)/swc*sw));
krw=@(sw)((sw<=1-sor).*krw0.*sws(sw).^nw+(sw>1-sor).*(-(1-krw0)/sor.*(1.0-sw)+1.0));
dkrwdsw=@(sw)((sw<=1-sor).*nw.*krw0.*(1/(1-sor-swc)).*sws(sw).^(nw-1)+(sw>1-sor)*((1-krw0)/sor));
dkrodsw=@(sw)((sw>=swc).*(-kro0*no*(1-sws(sw)).^(no-1))/(-swc-sor+1)+(sw<swc).*((kro0-1)/swc));
p0 = 100e5; % [bar] pressure
pin = 150e5; % [bar] injection pressure at the left boundary
u_in= 1.0/(24*3600); % [m/s] equal to 1 m/day
sw0 = swc+0.1; % initial water saturation
sw_in = 1;
mu_oil = 2e-3; % [Pa.s] oil viscosity
mu_water = 1e-3; % [Pa.s] water viscosity
% reservoir
k0 = 0.1e-12; % [m^2] average reservoir permeability
phi0 = 0.2; % average porosity
clx=1.2;
cly=0.2;
V_dp=0.1; % Dykstra-Parsons coef.
perm_val= field2d(Nx,Ny,k0,V_dp,clx,cly);
k=createCellVariable(m, perm_val);
phi=createCellVariable(m, phi0);
lw = geometricMean(k)/mu_water;
lo = geometricMean(k)/mu_oil;
Lw = @(sw)(krw(sw));
Lo = @(sw)(k/mu_oil*kro(sw));
dLwdsw = @(sw)(k/mu_water*dkrwdsw(sw));
dLodsw = @(sw)(k/mu_oil*dkrodsw(sw));
%% Define the boundaries
BCp = createBC(m); % Neumann BC for pressure
BCs = createBC(m); % Neumann BC for saturation
% bottom left boundary pressure gradient
BCp.left.a(1)=(krw(sw_in)*lw.xvalue(1,1)+kro(sw_in)*lo.xvalue(1,1)); BCp.left.b(1)=0; BCp.left.c(1)=-u_in;
BCp.bottom.a(1)=(krw(sw_in)*lw.yvalue(1,1)+kro(sw_in)*lo.yvalue(1,1)); BCp.bottom.b(1)=0; BCp.bottom.c(1)=-u_in;
% change the top right boandary to constant pressure (Dirichlet)
BCp.right.a(end)=0; BCp.right.b(end)=1; BCp.right.c(end)=p0;
BCp.top.a(end)=0; BCp.top.b(end)=1; BCp.top.c(end)=p0;
% change the bottom left boundary to constant saturation (Dirichlet)
BCs.left.a(1)=0; BCs.left.b(1)=1; BCs.left.c(1)=1;
BCs.bottom.a(1)=0; BCs.bottom.b(1)=1; BCs.bottom.c(1)=1;
%% define the time step and solver properties
% dt = 1000; % [s] time step
dt=(Lx/Nx)/u_in/10; % [s]

eps_p = 1e-5; % pressure accuracy
eps_sw = 1e-5; % saturation accuracy
%% define the variables
sw_old = createCellVariable(m, sw0, BCs);
p_old = createCellVariable(m, p0, BCp);
sw = sw_old;
p = p_old;
uw = -gradientTerm(p_old); % an estimation of the water velocity
%% start the main loop
% generate intial pressure profile (necessary to initialize the fully
% implicit solver)
dp_alwd= 100.0; % Pa
dsw_alwd= 0.1;
t = 0;
oil_init = domainInt(1-sw_old);
R(1) = 0;
i=0;
v_domain = domainInt(phi);
n_pv = 4; % injecting 4 pv of water
t_end = n_pv*v_domain/(u_in*well_radius*2); % final simulation time
while (t<t_end)
    error_p = 1e5;
    error_sw = 1e5;
    % Implicit loop
    loop_count=0;
    while ((error_p>eps_p) || (error_sw>eps_sw))
        loop_count=loop_count+1;
        if loop_count>10
            break
        end
        % calculate parameters
        pgrad = gradientTerm(p);
        sw_face = upwindMean(sw, -pgrad); % average value of water saturation
        labdao = lo.*funceval(kro, sw_face);
        labdaw = lw.*funceval(krw, sw_face);
        dlabdaodsw = lo.*funceval(dkrodsw, sw_face);
        dlabdawdsw = lw.*funceval(dkrwdsw, sw_face);
        labda = labdao+labdaw;
        dlabdadsw = dlabdaodsw+dlabdawdsw;
        % compute [Jacobian] matrices
        Mdiffp1 = diffusionTerm(-labda);
        Mdiffp2 = diffusionTerm(-labdaw);
        Mconvsw1 = convectionUpwindTerm(-dlabdadsw.*pgrad);
        Mconvsw2 = convectionUpwindTerm(-dlabdawdsw.*pgrad);
        [Mtranssw2, RHStrans2] = transientTerm(sw_old, dt, phi);
        % Compute RHS values
        RHS1 = divergenceTerm(-dlabdadsw.*sw_face.*pgrad);
        RHS2 = divergenceTerm(-dlabdawdsw.*sw_face.*pgrad);
        % include boundary conditions
        [Mbcp, RHSbcp] = boundaryCondition(BCp);
        [Mbcsw, RHSbcsw] = boundaryCondition(BCs);
        % Couple the equations; BC goes into the block on the main diagonal
        M = [Mdiffp1+Mbcp Mconvsw1; Mdiffp2 Mconvsw2+Mtranssw2+Mbcsw];
        RHS = [RHS1+RHSbcp; RHS2+RHStrans2+RHSbcsw];
        % solve the linear system of equations
        x = M\RHS;
        % x = agmg(M, RHS, [], 1e-10, 500, [], [p.value(:); sw.value(:)]);
        % separate the variables from the solution
        p_new = reshapeCell(m,full(x(1:(Nx+2)*(Ny+2))));
        sw_new = reshapeCell(m,full(x((Nx+2)*(Ny+2)+1:end)));
        % calculate error values
        error_p = max(abs((p_new(:)-p.value(:))./p_new(:)));
        error_sw = max(abs(sw_new(:)-sw.value(:)));
        % assign new values of p and sw
        p.value = p_new;
        sw.value = sw_new;
    end
    if loop_count>10
      p=p_old;
      sw=sw_old;
      dt=dt/5;
      continue
    end
    dsw=max(abs(sw_new(:)-sw_old.value(:))./sw_new(:));
    t=t+dt
    dt=min([dt*(dsw_alwd/dsw), 2*dt, t_end-t]);
    p_old = p;
    sw_old = sw;
    oil = domainInt(1-sw);
    i=i+1;
    R(i) = (oil_init-oil)/oil_init;
    p_inj(i) = p.value(2,2);
    t_series(i) = t;
    figure(1);visualizeCells(sw); drawnow;
end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
%function [pval, t_orig, crit_t, est_alpha, seed_state]=mult_comp_perm_t1(data,n_perm,tail,alpha_level,mu,reports,seed_state)
%
% mult_comp_perm_t1-One sample/paired sample permutation test based on a
% t-statistic.  This function can perform the test on one variable or
% simultaneously on multiple variables.  When applying the test to multiple
% variables, the "tmax" method is used for adjusting the p-values of each
% variable for multiple comparisons (Blair & Karnisky, 1993).  Like
% Bonferroni correction, this method adjusts p-values in a way that controls
% the family-wise error rate.  However, the permutation method will be more
% powerful than Bonferroni correction when different variables in the test
% are correlated.
%
% Required Input:
%  data   - 2D matrix of data (Observation x Variable)
%
% Optional Inputs:
%  n_perm      - Number of permutations used to estimate the distribution of
%                the null hypothesis. If the number of observations is less
%                than or equal to 12, all possible permutations are used and
%                this optional input has no effect.  If the number of
%                observations is greater than 12, n_perm specifies the
%                number of random permutations computed. Manly (1997)
%                suggests using at least 1000 permutations for an alpha level
%                of 0.05 and at least 5000 permutations for an alpha level of
%                0.01. {default=5000}
%  alpha_level - Desired family-wise alpha level. Note, because of the finite
%                number of possible permutations, the exact desired family-wise
%                alpha may not be possible. Thus, the closest approximation
%                is used and output as est_alpha. {default=.05}
%  tail        - [1, 0, or -1] If tail=1, the alternative hypothesis is that the
%                mean of the data is greater than 0 (upper tailed test).  If tail=0,
%                the alternative hypothesis is that the mean of the data is different
%                than 0 (two tailed test).  If tail=-1, the alternative hypothesis
%                is that the mean of the data is less than 0 (lower tailed test).
%                {default: 0}
%  mu          - The mean of the null hypothesis.  Must be a scalar or
%                1 x n vector (where n=the number of variables). {default: 0}
%  reports     - [0 or 1] If 0, function proceeds with no command line
%                reports. Otherwise, function reports what it is doing to
%                the command line. {default: 1}
%  seed_state  - The initial state of the random number generating stream
%                (see MATLAB documentation for "randstream"). If you pass
%                a value from a previous run of this function, it should
%                reproduce the exact same values. Note, this input only has
%                an effect if you have more than 12 observations.
%
% Outputs:
%  pval       - p-value (adjusted for multiple comparisons) of each
%               variable
%  t_orig     - t-score for each variable
%  crit_t     - Lower and upper critical t-scores for given alpha level.
%               t-scores that exceed these values significantly deviate from
%               the null hypothesis.  For upper tailed tests, the lower
%               critical t-score is NaN. The opposite is true of lower
%               tailed tests.
%  est_alpha  - The estimated family-wise alpha level of the test.  With
%               permutation tests, a finite number of p-values are possible.
%               This function tries to use an alpha level that is as close
%               as possible to the desired alpha level.  However, if the
%               sample size is small, a very limited number of p-values are
%               possible and the desired family-wise alpha level may be
%               impossible to approximately achieve.
%  seed_state - The initial state of the random number generating stream
%               (see MATLAB documentation for "randstream") used to
%               generate the permutations. You can use this to reproduce
%               the output of this function.  If the number of observations
%               is less than or equal to 12, seed_state='exact' since
%               all possible permutations are used in lieu of random
%               permutations.
%
% Note:
% -Unlike a parametric test (e.g., an ANOVA), a discrete set of p-values
% are possible (at most the number of possible permutations).  Since the
% number of possible permutations grows exponentially with the number of
% participants, this is only an issue for small sample sizes (e.g., 6
% participants).  When you have such a small sample size, the
% limited number of p-values may make the test overly conservative (e.g.,
% you might be forced to use an alpha level of .0286 since it is the biggest
% possible alpha level less than .05).
%
% -The null hypothesis of the permutation test is that the data come from a
% distribution that is symmetric around the mean of the null hypothesis.
% This is probably a generally reasonable assumption for
% paired-sample/repeated measures tests, but it might not be appropriate for
% one-sample tests.
%
%
% One Sample Example:
% >> data=randn(16,5);  %5 variables, 16 observations
% >> data(:,1:2)=data(:,1:2)+1; %mean of first two variables is 1
% >> [pval, t_orig, crit_t, est_alpha, seed_state]=mult_comp_perm_t1(data,50000);
% >> disp(pval); %adjusted p-values
%
% Paired-Sample/Repated Measures Example:
% >> dataA=randn(16,5);  %data from Condition A (5 variables, 16 observations)
% >> dataA(:,1:2)=dataA(:,1:2)+1; %mean of first two variables is 1
% >> dataB=randn(16,5); %data from Condition B (all variables have mean of 0)
% >> dif=dataA-dataB; %difference between conditions
% >> [pval, t_orig, crit_t, est_alpha, seed_state]=mult_comp_perm_t1(dif,50000);
% >> disp(pval); %adjusted p-values
%
% References:
% Blair, R.C. & Karniski, W. (1993) An alternative method for significance
% testing of waveform difference potentials. Psychophysiology.
%
% Manly, B.F.J. (1997) Randomization, Bootstrap, and Monte Carlo Methods in
% Biology. 2nd ed. Chapman and Hall, London.
%
%
%
% For a review on permutation tests and other contemporary techniques for
% correcting for multiple comparisons see:
%
%   Groppe, D.M., Urbach, T.P., & Kutas, M. (2011) Mass univariate analysis
% of event-related brain potentials/fields I: A critical tutorial review.
% Psychophysiology, 48(12) pp. 1711-1725, DOI: 10.1111/j.1469-8986.2011.01273.x
% http://www.cogsci.ucsd.edu/~dgroppe/PUBLICATIONS/mass_uni_preprint1.pdf
%
%
% Author:
% David Groppe
% Dec, 2010
% Kutaslab, San Diego
%

function [pval, t_orig, crit_t, est_alpha, seed_state]=mult_comp_perm_t1(data,n_perm,tail,alpha_level,mu,reports,seed_state)

if nargin<1,
    error('You need to provide data.');
end

if nargin<2,
    n_perm=5000;
end

if nargin<3,
    tail=0;
elseif (tail~=0) && (tail~=1) && (tail~=-1),
    error('Argument ''tail'' needs to be 0,1, or -1.');
end

if nargin<4,
    alpha_level=0.05;
elseif (alpha_level>=1) || (alpha_level<=0),
   error('Argument ''alpha_level'' needs to be a number between 0 and 1.');
end

if nargin<5,
    mu=0;
end

if nargin<6,
    reports=1;
end

defaultStream=RandStream.getDefaultStream; %random # generator state
if (nargin<7) || isempty(seed_state),
    %Store state of random number generator
    seed_state=defaultStream.State;
else
    defaultStream.State=seed_state; %reset random number generator to saved state
end

[n_obs n_var]=size(data);
if n_obs<2,
    error('You need data from at least two observations to perform a hypothesis test.')
end

if n_obs<7,
    n_psbl_prms=2^n_obs;
    if reports,
        watchit(sprintf(['Due to the very limited number of observations,' ...
            ' the total number of possible permutations is small.\nThus only a limited number of p-values (at most %d) are possible and the test might be overly conservative.'], ...
            n_psbl_prms));
    end
end

%% Remove null hypothesis mean from data
if isscalar(mu),
    data=data-mu;
elseif isvector(mu)
    s_mu=size(mu);
    if s_mu(1)>1,
        mu=mu';
        s_mu=size(mu);
    end
    if s_mu(2)~=n_var,
        error('mu needs to be a scalar or a 1 x %d vector (%d is the number of variables).',n_var,n_var);
    end
    data=data-repmat(mu,n_obs,1);
else
    error('mu needs to be a scalar or a 1 x %d vector (%d is the number of variables).',n_var,n_var);
end

if reports,
    fprintf('mult_comp_perm_t1: Number of variables: %d\n',n_var);
    fprintf('mult_comp_perm_t1: Number of observations: %d\n',n_obs);
    fprintf('t-score degrees of freedom: %d\n',n_obs-1);
end

%% Set up permutation test
if n_obs<=12,
    n_perm=2^n_obs; %total number of possible permutations
    exact=1;
    seed_state='exact';
    if reports,
        fprintf('Due to the limited number of observations, all possible permutations of the data will be computed instead of random permutations.\n');
    end
else
    exact=0;
end

if reports,
    fprintf('Executing permutation test with %d permutations...\n',n_perm);
    fprintf('Permutations completed: ');
end

%% Compute permutations

%Constant factor for computing t, speeds up computing t to precalculate
%now
sqrt_nXnM1=sqrt(n_obs*(n_obs-1));
if exact,
    %Use all possible permutations
    mxt=zeros(1,n_perm);
    for perm=1:n_perm
        if ~rem(perm,100)
            if reports,
                if ~rem(perm-100,1000)
                    fprintf('%d',perm);
                else
                    fprintf(', %d',perm);
                end
                if ~rem(perm,1000)
                    fprintf('\n');
                end
            end
        end
        %set sign of each participant's data
        if exist('de2bi.m','file') %% ?? check
            temp=de2bi(perm-1);
        else %% check ??
            temp=perm-1; %% check ??
        end %% check ??
        n_temp=length(temp);
        sn=-ones(n_obs,1);
        sn(1:n_temp,1)=2*temp-1;
        sn_mtrx=repmat(sn,1,n_var);
        d_perm=data.*sn_mtrx;

        %computes t-score of permuted data across all channels and time points
        sm=sum(d_perm,1);
        mn=sm/n_obs;
        sm_sqrs=sum(d_perm.^2,1)-(sm.^2)/n_obs;
        stder=sqrt(sm_sqrs)/sqrt_nXnM1;
        t=mn./stder;

        %get most extreme t-score
        [dummy mxt_id]=max(abs(t));
        mxt(perm)=t(mxt_id); %get the most extreme t-value with its sign (+ or -)
    end
else
    %Use random permutations
    mxt=zeros(1,n_perm*2);
    for perm=1:n_perm
        if ~rem(perm,100)
            if reports,
                if ~rem(perm-100,1000)
                    fprintf('%d',perm);
                else
                    fprintf(', %d',perm);
                end
                if ~rem(perm,1000)
                    fprintf('\n');
                end
            end
        end
        %randomly set sign of each participant's data
        sn=(rand(n_obs,1)>.5)*2-1;
        sn_mtrx=repmat(sn,1,n_var);

        d_perm=data.*sn_mtrx;

        %computes t-score of permuted data across all channels and time points
        sm=sum(d_perm,1);
        mn=sm/n_obs;
        sm_sqrs=sum(d_perm.^2,1)-(sm.^2)/n_obs;
        stder=sqrt(sm_sqrs)/sqrt_nXnM1;
        t=mn./stder;

        %get most extreme t-score (sign isn't immportant since we asumme
        %symmetric distribution of null hypothesis for one sample test)
        mxt(perm)=max(abs(t));
    end
    mxt(n_perm+1:2*n_perm)=-mxt(1:n_perm); %add the negative of all values since we assumme
    %null hypothesis distribution is symmetric
end

%End permutations completed line
if reports && rem(perm,1000)
    fprintf('\n');
end

%% Computes t-scores of observations at all variables and time points
sm=sum(data,1);
mn=sm/n_obs;
sm_sqrs=sum(data.^2,1)-(sm.^2)/n_obs;
stder=sqrt(sm_sqrs)/sqrt_nXnM1;
t_orig=mn./stder;

%% Compute p-values
pval=zeros(1,n_var);
for t=1:n_var,
    if tail==0,
        pval(t)=mean(mxt>=abs(t_orig(t)))*2;
    elseif tail==1,
        pval(t)=mean(mxt>=t_orig(t));
    elseif tail==-1,
        pval(t)=mean(mxt<=t_orig(t));
    end
end

%% Compute critical t-scores for specified alpha level,
if tail==0,
    %two-tailed
    crit_t(1)=prctile(mxt,100*alpha_level/2);
    crit_t(2)=-crit_t(1);
    est_alpha=mean(mxt>=crit_t(2))*2;
elseif tail==1,
    %upper tailed
    crit_t(1)=NaN;
    crit_t(2)=prctile(mxt,100-100*alpha_level);
    est_alpha=mean(mxt>=crit_t(2));
else
    %tail=-1, lower tailed
    crit_t(1)=prctile(mxt,alpha_level*100);
    est_alpha=mean(mxt<=crit_t(1));
    crit_t(2)=NaN;
end
if reports,
    fprintf('Desired family-wise alpha level: %f\n',alpha_level);
    fprintf('Estimated actual family-wise alpha level for returned values of crit_t: %f\n',est_alpha);
end

function watchit(msg)
%function watchit(msg)
%
% Displays a warning message on the Matlab command line.  Used by
% several Mass Univariate ERP Toolbox functions.
%

disp(' ');
disp('****************** Warning ******************');
disp(msg);
disp(' ');

author: Sable

repo: mcbench-benchmarks


text:
function msm_to_mm_coordinate_complex_hermitian ( output_filename, a )

%*****************************************************************************80
%
%% MSM_TO_MM_COORDINATE_COMPLEX_HERMITIAN writes a "matrix coordinate complex hermitian" Matrix Market file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 November 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the name of the file to which the information
%    should be written.
%
%    Input, sparse matrix A, the NROW by NCOL matrix, stored in MATLAB sparse
%    matrix format, which is to be written to the file.
%
  [ nrow, ncol ] = size ( a );
  nnzeros = nnz ( a );

  fid = fopen ( output_filename, 'wt+' );

  if ( fid < 0 );
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MSM_TO_MM_COORDINATE_COMPLEX_HERMITIAN - Fatal error!\n' );
    fprintf ( 1, '  Cannot open the output file.\n' );
    error ( 'MSM_TO_MM_COORDINATE_COMPLEX_HERMITIAN - Fatal error!');
  end;

  fprintf ( fid, '%%%%MatrixMarket matrix coordinate complex hermitian\n');
  fprintf ( fid, '%%%%  Created by MSM_TO_MM_COORDINATE_COMPLEX_HERMITIAN.M\n' );
  fprintf ( fid, '  %d  %d  %d\n', nrow, ncol, nnzeros );

  for j = 1 : ncol
    [ rows, temp, vals ] = find ( a(:,j) );
    sz = size ( rows );
    for k2 = 1 : sz
      if ( j <= rows(k2) )
        fprintf ( fid, '  %d  %d  %f  %f\n', rows(k2), j, real ( vals(k2) ), imag ( vals(k2) ) );
      end
    end

  end

  fclose ( fid );

  return
end


text:
filename='Cantileverbeam_Tetrahedra_Linear_Structured_SYM';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'holes';
cost = {'compliance'};
weights = [1];
constraint = {'volume'};
optimizer = 'HAMILTON-JACOBI';
incrementFactor = 1;
designVariable = 'LevelSet';
filterType = 'P1';
constraint_case = 'INEQUALITY';
line_search_initiator = 'INCREASING LAST STEP';
showBC = true;

HJiter0 = 1;
e2 = 100;
N_holes = [12 5 5];
R_holes = 0.7;
phase_holes = [0 0 0];
nsteps = 10;
Vfrac_final = 0.1;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

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

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function traveltime(dbName,varargin)

%TRAVELTIME make travel time plots.
% TRAVELTIME(dbName) creates travel time plots for database dbName. The top
% plot for P and S plots ditance vs. traveltime in seconds adjusted by a
% velocity reduction.
%
% TRAVELTIME(dbName,[Vp Vs]) use the specified velocities to reduce the
% times on the travel time plots. The default values are Vp = 7 and Vs = 4.
%
% see also ttimes.dbload

% Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
% $Date$
% $Revision$

% GET ARGUMENTS
if length(varargin)==1
    VpVs = varargin{1};
else
    VpVs = [7 4];
end

% LOAD DATABASE
[origin,site,arrival,ray] = ttimes.dbload(dbName);

% TRAVEL TIME PLOTS
figure('Position',[0 0 1100 850],'Color','w');
set(gcf,'DefaultAxesFontSize',14);
set(gcf,'DefaultAxesLineWidth',0.25);

% P WAVES
f = find(strcmp(arrival.iphase,'P'));
h1 = subplot(2,2,1);
scatter(ray.flatDist(f),arrival.travelTime(f)-ray.flatDist(f)/VpVs(1),30,ray.originDepth(f),'filled','MarkerEdgeColor','k');
hold on; box on; grid on;
xlabel('Distance (km)');
ylabel(['Traveltime - distance/' num2str(VpVs(1)) 'km/s (s)']);
title('P wave travel times');
xlim1 = get(gca,'xlim');
ylim1 = get(gca,'ylim');
%
h3 = subplot(2,2,3);
scatter(ray.flatDist(f),arrival.timeres(f),30,ray.originDepth(f),'filled','MarkerEdgeColor','k');
hold on; box on; grid on;
xlabel('Distance (km)');
ylabel('Time residual (s)');
title('P wave travel time residuals');
xlim3 = get(gca,'xlim');
ylim3 = get(gca,'ylim');

% S WAVES
f = find(strcmp(arrival.iphase,'S'));
h2 = subplot(2,2,2);
scatter(ray.flatDist(f),arrival.travelTime(f)-ray.flatDist(f)/VpVs(2),30,ray.originDepth(f),'filled','MarkerEdgeColor','k');
hold on; box on; grid on;
xlabel('Distance (km)');
ylabel(['Traveltime - distance/' num2str(VpVs(2)) 'km/s (s)']);
title('S wave travel times');
xlim2 = get(gca,'xlim');
ylim2 = get(gca,'ylim');
%
h4 = subplot(2,2,4);
scatter(ray.flatDist(f),arrival.timeres(f),30,ray.originDepth(f),'filled','MarkerEdgeColor','k');
hold on; box on; grid on;
xlabel('Distance (km)');
ylabel('Time residual (s)');
title('S wave travel time residuals');
xlim4 = get(gca,'xlim');
ylim4 = get(gca,'ylim');
%
h = colorbar('Location','east');
cmap = hot;
%colormap(flipud(cmap));
colormap(cmap);
position = get(h,'Position');
position(3) = position(3)/2;
position(4) = position(4)/2;
set(h,'Position',position);
set(h,'YDir','reverse');
set(h,'FontSize',9);
hh = get(h,'YLabel');
set(hh,'String','Depth (km)');
%
xLim = [0 max(ray.flatDist)];
yLim12 = [ min([ylim1 ylim2]) max([ylim1 ylim2])];
yLim34 = [ min([ylim3 ylim4]) max([ylim3 ylim4])];
set(h1,'xlim',xLim,'ylim',yLim12);
set(h2,'xlim',xLim,'ylim',yLim12);
set(h3,'xlim',xLim,'ylim',yLim34);
set(h4,'xlim',xLim,'ylim',yLim34);
%
set(gcf, 'paperorientation', 'landscape');
set(gcf, 'paperposition', [.5 .5 10 7.5] );
print(gcf, '-dpsc2', 'FIG_tt_curve.ps');

author: geoscience-community-codes

repo: GISMO


text:
function measPara = convertToMeasPara( drecksMDH, measPara )
%CONVERTTOMEASPARA convert from drecksMDH 2.0-struct to measPara-struct,
% but measPara-struct now contains just all needed variables
%
% (c) Thomas Kuestner
% ---------------------------------------------------------------------

if(drecksMDH.Seq.Is3D)
    measPara.dimensionSeq = '3D';
else
    measPara.dimensionSeq = '2D';
end
measPara.sliceThickness = drecksMDH.Geo.FOV(3)/drecksMDH.Geo.MatrixSize(3);
measPara.lines = drecksMDH.Geo.MatrixSize(2);
measPara.partitions = drecksMDH.Geo.MatrixSize(3);
measPara.baseResolution = drecksMDH.Geo.MatrixSize(1);
measPara.FOVread = drecksMDH.Geo.FOV(1);
measPara.FOVphase = drecksMDH.Geo.FOV(2);
measPara.sequenceName = drecksMDH.Seq.Sequence;
if(isfield(drecksMDH.Geo,'FFTLength'))
    measPara.ftlenPhase = drecksMDH.Geo.FFTLength(1);
    measPara.ftlenPar = drecksMDH.Geo.FFTLength(2);
else
    measPara.ftlenPhase = [];
    measPara.ftlenPar = [];
end
measPara.imageWidth = [];
measPara.imageHeight = [];

if(isfield(drecksMDH,'Accel') && isfield(drecksMDH.Accel,'EspressoDir'))
    dir_mapping = {'y', 'z', 'off'};
    measPara.ESPReSSoDirection = dir_mapping{drecksMDH.Accel.EspressoDir};
    pfn_mapping = [81, 90, 99, 108, 117; 1, 0.5, 0.625, 0.75, 0.875];
    measPara.ESPReSSo = pfn_mapping(2,pfn_mapping(1,:) == drecksMDH.Accel.EspressoFactor);
else
    measPara.ESPReSSoDirection = 'off';
    measPara.ESPReSSo = 1;
end

if(isfield(drecksMDH,'Wip'))
    if(isfield(drecksMDH.Wip,'SamplingFactor'))
        measPara.CSAcceleration = drecksMDH.Wip.SamplingFactor;
    end
    if(isfield(drecksMDH.Wip,'Phases'))
       	measPara.dim(4) = drecksMDH.Wip.Phases;
    end
    if(isfield(drecksMDH.Wip,'FullySampled'))
        measPara.CSFullySampled = drecksMDH.Wip.FullySampled;
    end
    if(isfield(drecksMDH.Wip,'Mask'))
        switch drecksMDH.Wip.Mask(1:2)
            case '06', measPara.CSSamplingType = 'Poisson';
            case '07', measPara.CSSamplingType = 'Random';
            case '08', measPara.CSSamplingType = 'Gaussian';
        end
        if(length(drecksMDH.Wip.Mask) > 3)
            switch drecksMDH.Wip.Mask(4:5)
                case '01', measPara.CSVDmap = 'None';
                case '02', measPara.CSVDmap = 'Point';
                case '03', measPara.CSVDmap = 'Block';
                case '04', measPara.CSVDmap = 'Ellipse';
                case '05', measPara.CSVDmap = 'Ring';
            end
    %         switch drecksMDH.Wip.Mask(7:8)
    %             case '14', measPara.CSBodyRegion = 'None';
    %             case '15', measPara.CSBodyRegion = 'Head';
    %             case '16', measPara.CSBodyRegion = 'Thorax';
    %             case '17', measPara.CSBodyRegion = 'Abdomen';
    %             case '18', measPara.CSBodyRegion = 'Pelvis';
    %         end
        end
    end
end

if(isfield(drecksMDH.Seq,'Bandwidth'))
    measPara.bandwidthPerPixel = drecksMDH.Seq.Bandwidth;
end

if(isfield(drecksMDH.Geo,'NormSag'))
    if(drecksMDH.Geo.NormSag)
        measPara.orientation    = 'sagittal';
    elseif(drecksMDH.Geo.NormCor)
        measPara.orientation    = 'coronal';
    elseif(drecksMDH.Geo.NormTra)
        measPara.orientation    = 'transversal';
    else
        measPara.orientation = 'rotated';
    end

    % pixel positions
    if(isfield(drecksMDH.Geo,'Shift'))
        measPara.position = drecksMDH.Geo.Shift;
    else
        measPara.position = zeros(1,3);
    end
    measPara.pixelSpacing = [drecksMDH.Geo.FOV(2)/drecksMDH.Geo.MatrixSize(2), drecksMDH.Geo.FOV(1)/drecksMDH.Geo.MatrixSize(1), drecksMDH.Geo.FOV(3)/drecksMDH.Geo.MatrixSize(3)]; % y-x-z/slice
    center = [drecksMDH.Geo.MatrixSize(2)/2, drecksMDH.Geo.MatrixSize(1)/2, drecksMDH.Geo.MatrixSize(3)/2];

    switch measPara.orientation
        case 'sagittal'
            measPara.xPos = measPara.position(1) + ([1:drecksMDH.Geo.MatrixSize(3)].' - center(3)) .* measPara.pixelSpacing(3);
            measPara.yPos = measPara.position(2) + ([1:drecksMDH.Geo.MatrixSize(2)].' - center(1)) .* measPara.pixelSpacing(1);
            measPara.zPos = measPara.position(3) + ([1:drecksMDH.Geo.MatrixSize(1)].' - center(2)) .* measPara.pixelSpacing(2);
        case 'coronal'
            measPara.xPos = measPara.position(1) + ([1:drecksMDH.Geo.MatrixSize(1)].' - center(2)) .* measPara.pixelSpacing(2);
            measPara.yPos = measPara.position(2) + ([1:drecksMDH.Geo.MatrixSize(3)].' - center(3)) .* measPara.pixelSpacing(3);
            measPara.zPos = measPara.position(3) + ([1:drecksMDH.Geo.MatrixSize(2)].' - center(1)) .* measPara.pixelSpacing(1);
        case 'transversal'
            measPara.xPos = measPara.position(1) + ([1:drecksMDH.Geo.MatrixSize(1)].' - center(2)) .* measPara.pixelSpacing(2);
            measPara.yPos = measPara.position(2) + ([1:drecksMDH.Geo.MatrixSize(2)].' - center(1)) .* measPara.pixelSpacing(1);
            measPara.zPos = measPara.position(3) + ([1:drecksMDH.Geo.MatrixSize(3)].' - center(3)) .* measPara.pixelSpacing(3);
        otherwise
            %
    end
end

% correct values
if(~isfield(measPara,'dim'))
    measPara.dim = zeros(1,5);
    measPara.LCall = zeros(1,4);
end
measPara.dim(1) = measPara.lines;
measPara.dim(2) = drecksMDH.Geo.OverSampling(1) * measPara.baseResolution; % due to CS_Trufi prescans
measPara.dim(3) = measPara.partitions;
% if(strcmp(measPara.sequenceName, 'CS_Trufi') && ~isempty(iLC)) % number of kernels (SetLC) = cardiac phases = time domain
%     measPara.dim(4) = double(max(iLC(:,10)) - min(iLC(:,10)) + 1); % SetLC
% end
measPara.LCall(2) = drecksMDH.LC.Averages;
if(measPara.dim(3) <= 1)
    measPara.dimension = '2D';
else
    if(measPara.dim(4) > 1)
        measPara.dimension = '4D';
    else
        measPara.dimension = '3D';
    end
end
measPara.dim(6) = 1;
if(isfield(drecksMDH,'Wip'))
    if(isfield(drecksMDH.Wip,'NSamples') && drecksMDH.Wip.NSamples > 0)
        measPara.dimension = '5D';
        measPara.dim(6) = drecksMDH.Wip.NSamples;
    elseif(isfield(drecksMDH.Wip,'NCardGates') && drecksMDH.Wip.NCardGates > 0)
        measPara.dimension = '5D';
        measPara.dim(6) = drecksMDH.Wip.NCardGates;
    end
end
if(isfield(measPara, 'sequenceName') && strcmp(measPara.sequenceName,'CS_Retro'))
    measPara.CSAcceleration = 2; % something unequal to 1
end

% local variables, just needed for further calculations
FOVphasePercentage = measPara.FOVphase/measPara.FOVread;
phaseOversampling = drecksMDH.Geo.OverSampling(2) - 1;
sliceOversampling = drecksMDH.Geo.OverSampling(3) - 1;

%% oversampling and anisotropy

% OVERSAMPLING PHASE INCORRECT!!!!!
% !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

% initialize anisotropic correction (y-x-z/slice-t)
measPara.aniso = [measPara.dim(1), measPara.dim(2), measPara.dim(3), measPara.dim(4), measPara.dim(6)];
% oversampling correction exactly possible => no case switch necessary later
if(~isfield(measPara,'oversampling'))
    measPara.oversampling = cell(2,3);
end
measPara.oversampling{2,2} = true;
measPara.oversampling{2,3} = true;

% readout oversampling
readoutDiff = (drecksMDH.Geo.OverSampling(1) - 1) * measPara.baseResolution;
measPara.oversampling{1,1} = readoutDiff/2+1:measPara.dim(2)-readoutDiff/2;
if(isempty(measPara.imageWidth) || measPara.imageWidth == 0)
    measPara.imageWidth = length(measPara.oversampling{1,1});
end

% correct anisotropy if oversampling correction was done during recon
if(measPara.oversampling{2,1})
    measPara.aniso(2) = length(measPara.oversampling{1,1});
end

% phase anistropy
if(isempty(measPara.ftlenPhase) || measPara.ftlenPhase == 0 || measPara.ftlenPhase ~= round(measPara.baseResolution * FOVphasePercentage * (1+phaseOversampling)))
    measPara.ftlenPhase = round(measPara.baseResolution * FOVphasePercentage * (1+phaseOversampling));
    if(measPara.ftlenPhase > measPara.baseResolution && FOVphasePercentage < 1) % can occur that y>x
        measPara.ftlenPhase = measPara.baseResolution;
    end
end
if(drecksMDH.Geo.PhaseRes(1) ~= 1)
    measPara.aniso(1) = measPara.ftlenPhase;
end
% phase oversampling
if(isempty(measPara.imageHeight) || measPara.imageHeight == 0)
    diffSize = abs(measPara.ftlenPhase - round(FOVphasePercentage * measPara.baseResolution));
    measPara.imageHeight = measPara.ftlenPhase;
else
    diffSize = abs(measPara.imageHeight - round(FOVphasePercentage * measPara.baseResolution));
end
if(mod(diffSize,2) == 0)
    oversamplimits = [round(diffSize/2)+1, round(diffSize/2)];
else
    oversamplimits = [round(diffSize/2)+1, 0];
    oversamplimits(2) = abs(diffSize-(oversamplimits(1)-1));
end
measPara.oversampling{1,2} = oversamplimits(1):(measPara.aniso(1)-oversamplimits(2));

if(strcmp(measPara.dimension,'3D') || strcmp(measPara.dimension,'4D') ||  strcmp(measPara.dimension,'5D'))
    % slice anisotropy
    if(isempty(measPara.ftlenPar) || measPara.ftlenPar == 0 || measPara.ftlenPar ~= round(drecksMDH.Geo.ImagesPerSlab * (1+sliceOversampling)))
        measPara.ftlenPar = round(drecksMDH.Geo.ImagesPerSlab * (1+sliceOversampling));
    end
    if(drecksMDH.Geo.PhaseRes(2) ~= 1)
        measPara.aniso(3) = measPara.ftlenPar;
    end
    % slice oversampling
    diffSize = abs(measPara.ftlenPar - drecksMDH.Geo.ImagesPerSlab);
    if(mod(diffSize,2) == 0)
        oversamplimits = [round(diffSize/2)+1, round(diffSize/2)];
    else
        oversamplimits = [round(diffSize/2)+1, 0];
        oversamplimits(2) = abs(diffSize-(oversamplimits(1)-1));
    end
    measPara.oversampling{1,3} = oversamplimits(1):(measPara.aniso(3)-oversamplimits(2));
%     measPara.imageDepth = measPara.imagesPerSlab;  %% ATTENTION !!!!!!!!!!! depth does not support time yet
else
    measPara.aniso(3) = measPara.LCall(1);
    measPara.oversampling{1,3} = 1:measPara.LCall(1);
%     measPara.imageDepth = measPara.LCall(1);     %% ATTENTION !!!!!!!!!!! depth does not support time yet
end

end

author: thomaskuestner


text:
 function [data M m] =scale_func(data,M,m)
%
% function data =rescale(data)
%
% This function rescale the input data between -1 and 1
%
% INPUT
%
% data: the data to rescale
% max: the maximum value of the ouput data
% min: the minimum value of the output data
%
% OUTPUT
%
% data: the rescaled data
[Nb_s Nb_b]=size(data);
if nargin==1
    M=max(data,[],1);
    m=min(data,[],1);
end

data = 2*(data-repmat(m,Nb_s,1))./(repmat(M-m,Nb_s,1))-1;

author: BehnoodRasti


text:
function f=crand(p1,p2);
%CRAND   Random complex numbers for testing.
%   Usage: f=tester_crand(p1,p2);

f=rand(p1,p2)-.5+i*(rand(p1,p2)-.5);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [t,a,k]=BSplineInterpDeriv(t,a,k,numDerivs)
%%BSPLINEINTERPDERIV Given a set of b-spline interpolation weights a and
%            knots t for scalar interpolation, compute the weights and
%            knots to interpolate the numDerivs derivative value using
%            b-splines interpolation.
%
%INPUTS: t The NX1 or 1XN set of knots for the interpolation function. The
%          first and last k-1 knots are outside of the ends of the region
%          with data or mark the ends of the region with data. These vslues
%          must be real.
%        a The npXnumSets collection of b-spline coefficients for
%          interpolating over a certain region for numSets sets of
%          interpolation problems. Such coefficients could be obtained from
%          the function BSplinePolyCoeffs along with t if one wishes to
%          interpolate over gridded data. These values can be complex.
%        k The order of the B-spline polynomials. The order of the
%          interpolation accuracy is k-numDerivs-1.
% numDerivs The number of derivatives to take. This must be less than k-1.
%
%OUTPUTS: a The (np-numDerivs)XnumSets modified coefficients.
%         t The (N-2*numDerivs)X1 modified knots.
%         k The modified B-spline order.
%
%This function implements Equation 12b in Chapter X of [1]. The
%modifications to t simply reflect the modification to the overall length
%of a.
%
%The first and last elements in t are not used in this function. However,
%the equations in [1] assume that those elements are present as they end up
%being used when computing integrals.
%
%EXAMPLE 1:
%Here, we approximate the derivative of a fourth-order function using
%piecewise third-order b-splines and plot the results.
% f=@(t)(t.^4-2*t.^2+t);
% tau=[-2;-1;0;0.5;1;1.5;2];
% y=f(tau);
% k=4;%k-1=3 is the order....
% [a,t]=BSplinePolyFit(tau,y,k);
%
% x=linspace(-2,2,500);
% numDerivs=1;
% [t,a,k]=BSplineInterpDeriv(t,a,k,numDerivs);
% ypDeriv=BSplineInterpVal(x,t,a,k);
% yDeriv=1-4*x+4*x.^3;%First derivative
%
% figure(1)
% clf
% hold on
% plot(x,yDeriv,'-k','linewidth',4)
% plot(x,ypDeriv,'--r','linewidth',2)
%
%EXAMPLE 2:
%This is similar to the previous example, except multiple sets of
%coefficients are computed and evaluated at once.
% f1=@(t)(t.^4-2*t.^2+t);
% f2=@(t)(t.^3-2*t.^2+t);
% tau=[-2;-1;0;0.5;1;1.5;2];
% y=[f1(tau),f2(tau)];
% k=4;%k-1=3 is the order....
% [a,t]=BSplinePolyFit(tau,y,k);
%
% x=linspace(-2,2,500).';
% numDerivs=1;
% [t,a,k]=BSplineInterpDeriv(t,a,k,numDerivs);
% ypDeriv=BSplineInterpVal(x,t,a,k);
% yDeriv=[1-4*x+4*x.^3,1-4*x+3*x.^2];%First derivative
%
% figure(1)
% clf
% hold on
% plot(x,yDeriv(:,1),'-k','linewidth',4)
% plot(x,ypDeriv(:,1),'--r','linewidth',2)
%
% figure(2)
% clf
% hold on
% plot(x,yDeriv(:,2),'-k','linewidth',4)
% plot(x,ypDeriv(:,2),'--r','linewidth',2)
%
%EXAMPLE 3:
%This is an example where a complex function parameterized by a single real
%value is differentiated and interpolated.
% numPoints=300;
% xMin=0;
% xMax=5;
% x=linspace(xMin,xMax,numPoints);
% f=@(x)exp(-1j*2*x).*(4*x.^3+x.*exp(1j*2*x.^2)+12*x+x.^3.*exp(-1j*2*x.^2));
%
% fDx=@(x)exp(-2*1j*x.*(1+x)).*(exp(4*1j*x.^2).*(1+2*1j*x.*(-1+2*x))+x.^2.*(3-2*1j*x.*(1+2*x))+4*exp(2*1j*x.^2).*(3+x.*(-6*1j+(3-2*1j*x).*x)));
%
% fx=f(x);
% k=5;
% [a,t]=BSplinePolyFit(x,fx(:),k);
% numDerivs=1;
% [t,a,k]=BSplineInterpDeriv(t,a,k,numDerivs);
%
% %Interpolate on a finer grid.
% numPoints=2000;
% x=linspace(xMin,xMax,numPoints);
% fx=fDx(x);
% fxInterp=BSplineInterpVal(x,t,a,k);
%
% %Plot the function values and the interpolated values for the real and
% %complex parts.
% figure(1)
% clf
% hold on
% plot(x,real(fx),'-k','linewidth',2)
% plot(x,imag(fx),'--k','linewidth',2)
% plot(x,real(fxInterp),'-r')
% plot(x,imag(fxInterp),'--r')
%
% %Plot the interpolation error.
% figure(2)
% clf
% hold on
% plot(x,real(fxInterp(:))-real(fx(:)))
% plot(x,imag(fxInterp(:))-imag(fx(:)))
%
%REFERENCES:
%[1] C. de Boor, A Practical Guide to Splines. New York: Springer-Verlag,
%    1978.
%
%April 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

numA=size(a,1);

for curDeriv=1:numDerivs
    aNext=a;
    for idxA=(curDeriv+1):numA
        denom=(t(idxA+k-curDeriv)-t(idxA))/(k-curDeriv);

        if(denom==0)
            aNext(idxA,:)=0;
        else
            aNext(idxA,:)=(a(idxA,:)-a(idxA-1,:))/denom;
        end
    end

    a=aNext;
end
a=a((1+numDerivs):end,:);
t=t((1+numDerivs):(end-numDerivs));
k=k-numDerivs;
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [Knots_n] = linear_3d_grid_from_displ(T, Spacing, volsz)
    ksz = ceil(volsz ./ Spacing) + 1;
    Knots_n =...
        T([1, 1+Spacing(1):Spacing(1):(ksz(1)-1)*Spacing(1), volsz(1)], ...
          [1, 1+Spacing(2):Spacing(2):(ksz(2)-1)*Spacing(2), volsz(2)], ...
          [1, 1+Spacing(3):Spacing(3):(ksz(3)-1)*Spacing(3), volsz(3)], :);
end

repo: pTVreg


text:
function f = p39_f ( n, x )

%*****************************************************************************80
%
%% P39_F evaluates the objective function for problem 39.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 January 2001
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Zbigniew Michalewicz,
%    Genetic Algorithms + Data Structures = Evolution Programs,
%    Third Edition,
%    Springer Verlag, 1996,
%    ISBN: 3-540-60676-9,
%    LC: QA76.618.M53.
%
%  Parameters:
%
%    Input, integer N, the number of variables.
%
%    Input, real X(N), the argument of the objective function.
%
%    Output, real F, the value of the objective function.
%
  f = x(1) * x(1) + 2.0 * x(2) * x(2) ...
    - 0.3 * cos ( 3.0 * pi * x(1) ) ...
    * cos ( 4.0 * pi * x(2) ) + 0.3;

  return
end


text:
function [u,sigma,eqn,info] = PoissonRT0(node,elem,bdFlag,pde,option)
%% POISSONRT0 Poisson equation: lowest order RT element.
%
%  [u,sigma] = PoissonRT0(node,elem,bdFlag,pde) produces an approximation of
%  the Poisson equation
%
%       -div(d*grad(u))=f  in \Omega, with
%       Dirichlet boundary condition u=g_D on \Gamma_D,
%       Neumann boundary condition   d*grad(u)*n=g_N on \Gamma_N
%
%  in the mixed formulation:
%
%  Find (\sigma , u) in H_{g_N,\Gamma_N}(div,\Omega)\times L^2(\Omega) s.t.
%
%  (d^-1\sigma,\tau) + (div \tau, u)  = <\tau*n,g_D>_{\Gamma_D}
%  \forall \tau in H_{0,\Gamma_N}(div,\Omega)
%  (div \sigma, v)                =  -(f,v)
%  \forall v in L^2(\Omega)
%
%  where
%  H_{g,\Gamma}(div,\Omega) = {\sigma \in H(div,\Omega); \sigma*n = g
%  on \Gamma \subset \partial\Omega }.
%
%  The unknown sigma = d*grad(u) is approximated using the lowest order
%  Raviart-Thomas element and u by piecewise constant element (with basis 1).
%
%  [u,sigma] = PoissonRT0(node,elem,bdFlag,pde,option) specifies options
%   - option.solver
%     'direct': the built in direct solver \ (mldivide)
%     'tri':    triangular preconditioner
%     'uzawapcg': PCG for the Schur complement equation
%     'none': only assemble the matrix equation but not solve
%
%   The default setting is to use the direct solver for small size problems
%   and transforming based multigrid solvers for large size problems.
%
%  Example
%
%    squarePoissonRT0
%
% Created by Ming Wang. Reorganized by Long Chen. Change basis for u.
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

%% Preprocess
if ~exist('bdFlag','var'), bdFlag = []; end
if ~exist('option','var'), option = []; end

%% Diffusion coefficient
time = cputime;  % record assembling time
if ~isfield(pde,'d'), pde.d = []; end
if isfield(pde,'d') && ~isempty(pde.d)
   if isnumeric(pde.d)
      K = pde.d;                   % d is an array
   else                            % d is a function
      center = (node(elem(:,1),:) + node(elem(:,2),:) + node(elem(:,3),:))/3;
      K = pde.d(center);  % take inverse sequencil.
   end
else
    K = [];
end

%% Data structure
elemold = elem;
[elem,bdFlag] = sortelem(elem,bdFlag);  % ascend ordering of elem
[elem2edge,edge] = dofedge(elem);
NT = size(elem,1); NE = size(edge,1);
[Dlambda,area,elemSign] = gradbasis(node,elem);

%% Assemble matrix
Nsigma = NE; Nu = NT; Ndof = Nsigma + Nu;

% M. Mass matrix for RT0 element
M = getmassmatvec(elem2edge,area,Dlambda,'RT0',K); % ascend ordering of loc edge

% B. divergence operator
B = icdmat(double(elem2edge),elemSign*[1 -1 1]); % inconsistency with the induced ordering

% C. zero matrix.
C = sparse(Nu,Nu);

A = [M B';B C];

%% Assemble right hand side.
fu = zeros(Nu,1);
if ~isfield(pde,'f') || (isreal(pde.f) && (pde.f==0))
    pde.f = [];
end
if ~isfield(option,'fquadorder')
    option.fquadorder = 2;   % default order
end
if ~isempty(pde.f)
	[lambda,w] = quadpts(option.fquadorder);
	nQuad = size(lambda,1);
    for p = 1:nQuad
		% quadrature points in the x-y coordinate
		pxy = lambda(p,1)*node(elem(:,1),:) ...
			+ lambda(p,2)*node(elem(:,2),:) ...
			+ lambda(p,3)*node(elem(:,3),:);
		fp = pde.f(pxy);
		fu = fu - fp*w(p);  % div u = -f;
    end
    fu = fu.*area;
end
clear fp
F((Nsigma+1):Ndof,1) = fu;

%% Boundary Conditions
if ~exist('bdFlag','var'), bdFlag = []; end
[AD,F,bigu,freeDof,isPureNeumannBC] = getbdRT0(F);
eqn = struct('M',AD(1:NE,1:NE),'B',AD(NE+1:end,1:NE),'C',AD(NE+1:end,NE+1:end),...
             'f',F(1:NE),'g',F(NE+1:end),'freeDof',freeDof,'A',AD);

%% Record assembling time
assembleTime = cputime - time;
if ~isfield(option,'printlevel'), option.printlevel = 1; end
if option.printlevel >= 2
    fprintf('Time to assemble matrix equation %4.2g s\n',assembleTime);
end

%% Solve the linear system.
% Set up solver type
if isempty(option) || ~isfield(option,'solver')    % no option.solver
    if Ndof <= 2e3  % Direct solver for small size systems
        option.solver = 'direct';
    else         % MGCG  solver for large size systems
        option.solver = 'tri';
    end
% elseif strcmp(option.solver,'mg')
%     option.solver = 'tripremixPoisson';
end
solver = option.solver;
% solve
switch lower(solver)
    case 'direct'
        t = cputime;
        bigu(freeDof) = AD(freeDof,freeDof)\F(freeDof);
        sigma = bigu(1:NE);
        u = bigu(NE+1:end);
        info = struct('solverTime',cputime - t,'itStep',1,'error',0,'flag',0,'stopErr',0);
    case 'none'
        sigma = zeros(NE,1); u = zeros(NT,1); info =[];
    case 'tri'
        [sigma,u,info] = tripremixPoisson(eqn.M,eqn.B,eqn.C,eqn.f,eqn.g,elemold);
    case 'uzawapcg'
        [sigma,u,info] = uzawapcg(eqn.M,eqn.B,eqn.C,eqn.f,eqn.g,elemold);
%     case 'mg'
%         option.freeEdge = freeEdge;
%         option.isPureNeumannBC = isPureNeumannBC;
%         [sigma0,u,info] = mgDarcy(eqn.M,eqn.B,eqn.f,eqn.g,elemold,option);
%         sigma = bigu(1:NE);
%         sigma(freeEdge) = sigma0;
end
if isPureNeumannBC % post process for u for pure Neumann boundary condition
    ubar = sum(u.*area)/sum(area);
    u = u - ubar;
end

%% Output information
info.assembleTime = assembleTime;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfunction getbdRT0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function [AD,F,bigu,freeDof,isPureNeumannBC] = getbdRT0(F)
    %% GETBDRT0 Boundary conditions for Poisson equation: RT0 element.
    %
    %  Created by Ming Wang. Improved the check of edgeSign by Long Chen.

    bigu = zeros(Ndof,1);

    %% No boundary conditions
    if ~isfield(pde,'g_D'), pde.g_D = []; end
    if ~isfield(pde,'g_N'), pde.g_N = []; end

    %% Set up bdFlag
    if isempty(bdFlag) % no bdFlag information
       if ~isempty(pde.g_N) % case: Neumann
           bdFlag = setboundary(node,elem,'Neumann');
       elseif ~isempty(pde.g_D) % case: Dirichlet
           bdFlag = setboundary(node,elem,'Dirichlet');
       end
    end

    %% Find Dirichlet and Neumann dofs
    edgeSign = ones(NE,1);
    isDirichlet = false(NE,1);
    isNeumann = false(NE,1);
    if ~isempty(bdFlag)
        isDirichlet(elem2edge(bdFlag(:)==1)) = true;
          isNeumann(elem2edge(bdFlag(:)==2)) = true;
        % Direction of boundary edges may not be the outwards normal direction
        % of the domain. edgeSign is introduced to record this inconsistency.
        edgeSign = ones(NE,1);
        idx = (bdFlag(:,1) ~= 0) & (elemSign == -1);% first edge is on boundary
        edgeSign(elem2edge(idx,1)) = -1;
        idx = (bdFlag(:,2) ~= 0) & (elemSign == 1); % second edge is on boundary
        edgeSign(elem2edge(idx,2)) = -1;
        idx = (bdFlag(:,3) ~= 0) & (elemSign == -1);% third edge is on boundary
        edgeSign(elem2edge(idx,3)) = -1;
    end
    Dirichlet = edge(isDirichlet,:);
    Neumann = edge(isNeumann,:);
    isBdDof = false(Ndof,1);
    isBdDof(isNeumann) = true;   % for mixed method, Neumann edges are fixed
    freeDof = find(~isBdDof);
%     isFreeEdge = true(NE,1);
%     isFreeEdge(isNeumann) = false;
%     freeEdge = find(isFreeEdge);

    %% Dirichlet boundary condition (Neumann BC in the mixed form)
    %   We need only modify the rhs on dof associated with Dirichlet
    %   boundary. Compute the int_e g_D \phi_e\cdot n ds on the boundary.
    if ~isempty(pde.g_D) && isnumeric(pde.g_D) && (pde.g_D==0)
        pde.g_D = [];
    end
    if ~isempty(pde.g_D) && any(isDirichlet)
        if ~isfield(option,'gNquadorder')
            option.gNquadorder = 2;   % default order exact for quadratic gN
        end
        [lambda,w] = quadpts1(option.gNquadorder);
        nQuad = size(lambda,1);
        % <\phi\cdot n, g_D> = 1/|e_{i,j}|\int e_{i,j} g_D ds
        for ip = 1:nQuad
        	pxy = lambda(ip,1)*node(Dirichlet(:,1),:)+...
                  lambda(ip,2)*node(Dirichlet(:,2),:);
            F(isDirichlet) = F(isDirichlet) + w(ip)*pde.g_D(pxy);
        end
        F(isDirichlet) = F(isDirichlet).*edgeSign(isDirichlet);
        % no edge length since the basis of sigma contains it.
    end

    %% Neumann boundary condition (Dirichlet BC in mixed form)
    % We compute the integral int_e g_N ds and assign to boundary dof
    if ~isempty(pde.g_N) && isnumeric(pde.g_N) && (pde.g_N==0)
        pde.g_N = [];
    end
    if ~isempty(pde.g_N) && any(isNeumann)
        % modify the rhs to include Dirichlet boundary condition
        ve = node(Neumann(:,1),:)-node(Neumann(:,2),:);
        edgeLength = sqrt(sum(ve.^2,2));
        % compute the integral int_e g_N ds
        if ~isfield(option,'gNquadorder')
            option.gNquadorder = 2;   % default order exact for quadratic gN
        end
        [lambda,w] = quadpts1(option.gNquadorder);
        nQuad = size(lambda,1);
        bigu(isNeumann) = 0;
        for ip = 1:nQuad
        	pxy = lambda(ip,1)*node(Neumann(:,1),:)+...
                  lambda(ip,2)*node(Neumann(:,2),:);
            bigu(isNeumann) = bigu(isNeumann) + w(ip)*pde.g_N(pxy).*edgeLength;
        end
        % correct the sign
        bigu(isNeumann) = bigu(isNeumann).*edgeSign(isNeumann);
        F = F - A*bigu;
        F(isNeumann) = bigu(isNeumann);
    end

    %% Pure Neumann boundary condition
    isPureNeumannBC = false;
    if ~any(isDirichlet) && any(isNeumann)
        freeDof = freeDof(1:end-1);  % eliminate the kernel by enforcing u(NT) = 0;
        isBdDof(end) = true;
        isPureNeumannBC = true;
        F(NE+1:end) = F(NE+1:end) - mean(F(NE+1:end)); % normalize
%         F(end) = 0;
    end

    %% Modify the matrix
    %  Build Neumann boundary condition(Dirichlet BC in mixed form) into the
    %  matrix AD by enforcing  |AD(bdNode,bdNode)=I,
    %  AD(bdNode,FreeNode)=0, AD(FreeNode,bdNode)=0|.
    if any(isBdDof)
       bdidx = zeros(Ndof,1);
       bdidx(isBdDof) = 1;
       Tbd = spdiags(bdidx,0,Ndof,Ndof);
       T = spdiags(1-bdidx,0,Ndof,Ndof);
       AD = T*A*T + Tbd;
    else
       AD = A;
    end
    end % end of getbdRT0
end

repo: ifem


text:
function y=PSNR(noisyImage,restoredImage)

noisyImage = double(noisyImage);
restoredImage = double(restoredImage);
% Compute the PSNR of two gray scale image
% Traditional progarmming using loops
% Class input : [0,1]
% july, 25 , 2012
% KHMOU Youssef

N=size(noisyImage);
if length(N)> 2
    error('Input must be grayscale image');
end
if size(noisyImage)~=size(restoredImage)
    error('The images must have the same size');
end

%if ~isa(noisyImage,'double')
%   noisyImage=double(noisyImage)./255.00;
%end
%if  ~isa(restoredImage,'double')
%    restoredImage=double(restoredImage)./255.00;
%end

% begin

d1=max(noisyImage(:));
d2=max(restoredImage(:));
d=max(d1,d2);

MSE=0;
for i=1:N(1)
    for j=1:N(2)
        if isnan(noisyImage(i,j)) || isinf(restoredImage(i,j))...
                || isnan(restoredImage(i,j)) || isinf(noisyImage(i,j))
            continue;
        end
        MSE=MSE+((abs(noisyImage(i,j)-restoredImage(i,j))).^2);
    end
end

MSE=MSE./(N(1)*N(2));

y=10*log10((d.^2) /MSE);

author: andrewssobral


text:
function [uxref, xref] = xsection(qmethod,xy,x_gal,yref,fig)
%xsection   plots/explores solution on horizontal cross-section
%   [uxref, xref] = xsection(qmethod,xy,x_gal,yref,fig);
%   input
%          qmethod    mixed method
%          xy         velocity nodal coordinate vector
%          x_gal   solution vector
%          yref       y-location of grid line
%          fig         figure number
%
%   IFISS function: DJS; 19 September 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
nvtx=length(xy(:,1));
kk=find(xy(:,2)==yref);
if isempty(kk), error('location yref is not a grid-line!'), end
uxref=x_gal(kk);  xref=xy(kk,1);
figure(fig)
plot(xref,uxref,'-k'), axis('square'), title('x-section of flow')
%%
%% compute volume of flow using appropriate quadrature
nny=length(xref); hy=xref(2)-xref(1);
if qmethod >1,
%% Simpson's rule
ww=ones(nny,1); ww(2:2:nny-1)=4; ww(3:2:nny-1)=2;
total=(hy/3)*sum(ww.*uxref);
else
%% Trapezium rule
ww=ones(nny,1); ww(2:2:nny-1)=2; ww(3:2:nny-1)=2;
total=(hy/2)*sum(ww.*uxref);
end
fprintf('\narea of x-section is %10.6e \n',total)
return


text:
function han = plot_vertical_line(x,color)
% ..
%    Edited by Tor, Feb 08, to plot multiple lines
% ..

if nargin < 2, color = 'k'; end

ylim = get(gca,'YLim');

for i = 1:length(x)
    han = plot([x(i) x(i)], ylim, color);

end

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function h = makemovie(matrix, cmap)
% h = makemovie(matrix, cmap)
h=figure;
% Resize the figure window to size of movie required.
% 3) Record the size of the plot window:
winsize = get(h,'Position');
% 4) Adjust size of this window to include the whole figure window (if you require the axes, title and axis labels
%      in the movie):
winsize(1:2) = [0 0];
% 5) Set the number of frames:
numframes=size(matrix,3);
% 6) Create the MATLAB movie matrix:
A=moviein(numframes,h,winsize);
% 7) Fix the features of the plot window (ensures each frame of the movie
% is the same size):
set(h,'NextPlot','replacechildren')
% 8) Within a loop, plot each picture and save to MATLAB movie matrix:
for i=1:numframes
   imagesc(matrix(:,:,i)); % plot command
   set (gca, 'DataAspectRatio',[1 1 1]);
   if exist ('cmap', 'var')
       colormap(cmap);
   end
   % add axis label, legends, titles, etc. in here
   A(:,i)=getframe(h,winsize);
end
% This procedure creates a movie stored in a special format that is only readable in MATLAB. The first thing you will want to do is to play the movie:
movie(h,A,30,3,winsize)

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB