text:
classdef Retina < handle
    %RETINA  A biological retina model for image spatio-temporal noise and luminance changes enhancement
    %
    % Class which allows the Gipsa/Listic Labs model to be used with OpenCV.
    %
    % This retina model allows spatio-temporal image processing (applied on
    % still images, video sequences).
    %
    % As a summary, these are the retina model properties:
    %
    % - It applies a spectral whithening (mid-frequency details enhancement)
    % - high frequency spatio-temporal noise reduction
    % - low frequency luminance to be reduced (luminance range compression)
    % - local logarithmic luminance compression allows details to be enhanced
    %   in low light conditions
    %
    % **USE**: this model can be used basically for spatio-temporal video
    % effects but also for:
    %
    % - using the `getParvo` method output matrix: texture analysis with
    %   enhanced signal to noise ratio and enhanced details robust against
    %   input images luminance ranges
    % - using the `getMagno` method output matrix: motion analysis also with
    %   the previously cited properties
    %
    % For more information, refer to the following papers [Benoit2010] and
    % [Herault2010].
    %
    % ## Retina class overview
    %
    % ### Introduction
    %
    % This class provides the main controls of the Gipsa/Listic labs human
    % retina model. This is a non separable spatio-temporal filter modelling
    % the two main retina information channels:
    %
    % - foveal vision for detailed color vision: the parvocellular pathway.
    % - peripheral vision for sensitive transient signals detection (motion
    %   and events): the magnocellular pathway.
    %
    % This model originates from Jeanny Herault work [Herault2010]. It has
    % been involved in Alexandre Benoit phd and his current research
    % [Benoit2010], [Benoit2014]. He currently maintains this module within
    % OpenCV. It includes the work of other Jeanny's phd student such as
    % [Chaix2007] and the log polar transformations of Barthelemy Durette
    % described in Jeanny's book.
    %
    % More into details here is an overview of the retina properties that are
    % implemented here:
    %
    % Regarding luminance and details enhancement:
    %
    % - local logarithmic luminance compression (at the entry point by
    %   photoreceptors and at the output by ganglion cells).
    % - spectral whitening at the Outer Plexiform Layer level (photoreceptors
    %   and horizontal cells spatio-temporal filtering).
    %
    % The former behavior compresses luminance range and allows very bright
    % areas and very dark ones to be visible on the same picture with lots of
    % details. The latter reduces low frequency luminance energy (mean
    % luminance) and enhances mid-frequencies (details). Applied all together,
    % retina well prepares visual signals prior high level analysis. Those
    % properties are really interesting with videos where light changes are
    % dramatically reduced with an interesting temporal consistency.
    %
    % Regarding noise filtering :
    %
    % - high frequency spatial and temporal noise is filtered out. Both
    %   outputs Parvo and Magno pathways benefit from this. Noise reduction
    %   benefits from the non separable spatio-temporal filtering.
    % - at the Parvo output, static textures are enhanced and noise is
    %   filtered (on videos, temporal noise is nicely removed). However, as
    %   human behaviors, moving textures are smoothed. Then, moving object
    %   details can be only enhanced  if the retina tracks it and keeps it
    %   static from its point of view.
    % - at Magno output, it allows a cleaner detection of events (motion,
    %   changes) with reduced noise errors even in difficult lighting
    %   conditions. As a compromise, the Magno output is a low spatial
    %   frequency signal and allows events' blobs to be reliably extracted
    %   (check the cv.TransientAreasSegmentationModule module for that).
    %
    % ### Use
    %
    % This model can be used as a preprocessing stage in the aim of:
    %
    % - performing texture analysis with enhanced signal to noise ratio and
    %   enhanced details which are robust against input images luminance
    %   ranges (check out the parvocellular retina channel output, by using
    %   the provided `getParvo` methods)
    % - performing motion analysis that is also taking advantage of the
    %   previously cited properties (check out the magnocellular retina
    %   channel output, by using the provided `getMagno` methods)
    % - general image/video sequence description using either one or both
    %   channels. An example of the use of Retina in a Bag of Words approach
    %   is given in [Benoit2014].
    %
    % Note:
    %
    % - For ease of use in computer vision applications, the two retina
    %   channels are applied on all the input images. This does not follow the
    %   real retina topology but it is practical from an image processing
    %   point of view. If retina mapping (foveal and parafoveal vision) is
    %   required, use the log sampling capabilities proposed within the class.
    % - Please do not hesitate to contribute by extending the retina
    %   description, code, use cases for complementary explanations and
    %   demonstrations.
    %
    % ### Use case illustrations
    %
    % #### Image preprocessing using the Parvocellular pathway (parvo retina output)
    %
    % As a preliminary presentation, let's start with a visual example. We
    % propose to apply the filter on a low quality color jpeg image with
    % backlight problems. Here is the considered input...
    % *"Well,i could see more with my eyes than what i captured with my camera..."*
    %
    % ![image](https://docs.opencv.org/3.3.1/retinaInput.jpg)
    %
    % Below, the retina foveal model applied on the entire image with default
    % parameters. Details are enforced whatever the local luminance is. Here
    % there contours are strongly enforced but the noise level is kept low.
    % Halo effects are voluntary visible with this configuration. See
    % parameters discussion below and increase `HorizontalCellsGain` near 1 to
    % remove them.
    %
    % ![image](https://docs.opencv.org/3.3.1/retinaOutput_default.jpg)
    %
    % Below, a second retina foveal model output applied on the entire image
    % with a parameters setup focused on naturalness perception.
    % *"Hey, i now recognize my cat, looking at the mountains at the end of the day!"*.
    % Here contours are enforced, luminance is corrected but halos are avoided
    % with this configuration. The backlight effect is corrected and highlight
    % details are still preserved. Then, even on a low quality jpeg image, if
    % some luminance's information remains, the retina is able to reconstruct
    % a proper visual signal. Such configuration is also useful for High
    % Dynamic Range (HDR) images compression to 8bit images as discussed in
    % [Benoit2010] and in the demonstration codes discussed below. As shown at
    % the end of the page, parameter changes from defaults are:
    %
    % - `HorizontalCellsGain = 0.3`
    % - `PhotoreceptorsLocalAdaptationSensitivity = GanglioncellsSensitivity = 0.89`.
    %
    % ![image](https://docs.opencv.org/3.3.1/retinaOutput_realistic.jpg)
    %
    % As observed in this preliminary demo, the retina can be settled up with
    % various parameters, by default, as shown on the figure above, the retina
    % strongly reduces mean luminance energy and enforces all details of the
    % visual scene. Luminance energy and halo effects can be modulated
    % (exaggerated to cancelled as shown on the two examples). In order to use
    % your own parameters, you can use at least one time the cv.Retina.write
    % method which will write a proper XML file with all default parameters.
    % Then, tweak it on your own and reload them at any time using method
    % cv.Retina.setup. These methods update a `RetinaParameters` member
    % structure that is described hereafter. XML parameters file samples are
    % shown at the end of the page.
    %
    % #### Tone mapping processing capability using the Parvocellular pathway (parvo retina output)
    %
    % This retina model naturally handles luminance range compression. Local
    % adaptation stages and spectral whitening contribute to luminance range
    % compression. In addition, high frequency noise that often corrupts tone
    % mapped images is removed at early stages of the process thus leading to
    % natural perception and noise free tone mapping.
    %
    % Compared to the demos shown above, setup differences are the following
    % ones:
    %
    % * load HDR images (OpenEXR format is supported by OpenCV) and cut
    %   histogram borders at ~5% and 95% to eliminate salt&pepper like pixel's
    %   corruption.
    % * apply retina with default parameters along with the following changes
    %   (generic parameters used for the presented illustrations of the
    %   section):
    %   * `HorizontalCellsGain=0.4` (the main change compared to the default
    %     configuration: it strongly reduces halo effects)
    %   * `PhotoreceptorsLocalAdaptationSensitivity=0.99` (a little higher
    %     than default value to enforce local adaptation)
    %   * `GanglionCellsSensitivity=0.95` (also slightly higher than default
    %     for local adaptation enforcement)
    % * get the parvo output using the cv.Retina.getParvo method.
    %
    % Have a look at the end of this page to see how to specify these
    % parameters in a configuration file.
    %
    % The following two illustrations show the effect of such configuration on
    % 2 image samples.
    %
    % - HDR image tone mapping example with generic parameters.
    %   Original image comes from [OpenEXR](http://openexr.com/) samples
    %   (`openexr-images-1.7.0/ScanLines/CandleGlass.exr`)
    %
    %   ![image](https://docs.opencv.org/3.3.1/HDRtoneMapping_candleSample.jpg)
    %
    % - HDR image tone mapping example with the same generic parameters.
    %   Original image comes from
    %   [memorial.exr](http://www.pauldebevec.com/Research/HDR/memorial.exr)
    %
    %   ![image](https://docs.opencv.org/3.3.1/HDRtoneMapping_memorialSample.jpg)
    %
    % #### Motion and event detection using the Magnocellular pathway (magno retina output)
    %
    % Spatio-temporal events can be easily detected using *magno* output of
    % the retina (use the cv.Retina.getMagno method). Its energy linearly
    % increases with motion speed.
    %
    % An event blob detector is proposed with the
    % cv.TransientAreasSegmentationModule class also provided in the
    % bioinspired module. The basic idea is to detect local energy drops with
    % regard of the neighborhood and then to apply a threshold. Such process
    % has been used in a bag of words description of videos on the TRECVid
    % challenge [Benoit2014] and only allows video frames description on
    % transient areas.
    %
    % We present here some illustrations of the retina outputs on some
    % examples taken from [CDNET](http://changedetection.net/) with RGB and
    % thermal videos.
    %
    % NOTE: here, we use the default retina setup that generates halos around
    % strong edges. Note that temporal constants allow a temporal effect to be
    % visible on moting objects (useful for still image illustrations of a
    % video). Halos can be removed by increasing retina Hcells gain while
    % temporal effects can be reduced by decreasing temporal constant values.
    % Also take into account that the two retina outputs are rescaled in range
    % `[0:255]` such that magno output can show a lot of "noise" when nothing
    % moves while drawing it. However, its energy remains low if you retrieve
    % it using cv.Retina.getMagnoRAW getter instead.
    %
    % - Retina processing on RGB image sequence: example from
    %   [CDNET](http://changedetection.net/) (baseline/PETS2006).
    %   Parvo enforces static signals but smooths moving persons since they do
    %   not remain static from its point of view.
    %   Magno channel highligths moving persons, observe the energy mapping on
    %   the one on top, partly behind a dark glass.
    %
    %   ![image](https://docs.opencv.org/3.3.1/VideoDemo_RGB_PETS2006.jpg)
    %
    % - Retina processing on gray levels image sequence: example from
    %   [CDNET](http://changedetection.net/) (thermal/park).
    %   On such grayscale images, parvo channel enforces contrasts while mango
    %   strongly reacts on moving pedestrians
    %
    %   ![image](https://docs.opencv.org/3.3.1/VideoDemo_thermal_park.jpg)
    %
    % ### Literature
    %
    % For more information, refer to the following papers:
    %
    % - Model description: [Benoit2010]
    % - Model use in a Bag of Words approach: [Benoit2014]
    % - Please have a look at the reference work of Jeanny Herault that you
    %   can read in his book: [Herault2010]
    %
    % This retina filter code includes the research contributions of
    % phd/research collegues from which code has been redrawn by the author:
    %
    % - take a look at the `retinacolor.hpp` module to discover Brice Chaix
    %   de Lavarene phD color mosaicing/demosaicing and his reference paper
    %   [Chaix2007]
    %
    % - take a look at `imagelogpolprojection.hpp` to discover retina spatial
    %   log sampling which originates from Barthelemy Durette phd with Jeanny
    %   Herault. A Retina / V1 cortex projection is also proposed and
    %   originates from Jeanny's discussions. More information in the above
    %   cited Jeanny Heraults's book.
    %
    % - Meylan et al. work on HDR tone mapping that is implemented as a
    %   specific method within the model [Meylan2007]
    %
    % ## Retina programming interfaces
    %
    % The proposed class allows the [Gipsa](http://www.gipsa-lab.inpg.fr)
    % (preliminary work) / [Listic](http://www.listic.univ-savoie.fr) labs
    % retina model to be used. It can be applied on still images, images
    % sequences and video sequences.
    %
    % ### Setting up Retina
    %
    % #### Managing the configuration file
    %
    % When using the cv.Retina.write and cv.Retina.load methods, you create or
    % load a XML file that stores Retina configuration.
    %
    % The default configuration is presented below.
    %
    % ```xml
    % <?xml version="1.0"?>
    % <opencv_storage>
    % <OPLandIPLparvo>
    %     <colorMode>1</colorMode>
    %     <normaliseOutput>1</normaliseOutput>
    %     <photoreceptorsLocalAdaptationSensitivity>7.5e-01</photoreceptorsLocalAdaptationSensitivity>
    %     <photoreceptorsTemporalConstant>9.0e-01</photoreceptorsTemporalConstant>
    %     <photoreceptorsSpatialConstant>5.7e-01</photoreceptorsSpatialConstant>
    %     <horizontalCellsGain>0.01</horizontalCellsGain>
    %     <hcellsTemporalConstant>0.5</hcellsTemporalConstant>
    %     <hcellsSpatialConstant>7.</hcellsSpatialConstant>
    %     <ganglionCellsSensitivity>7.5e-01</ganglionCellsSensitivity>
    % </OPLandIPLparvo>
    % <IPLmagno>
    %     <normaliseOutput>1</normaliseOutput>
    %     <parasolCells_beta>0.</parasolCells_beta>
    %     <parasolCells_tau>0.</parasolCells_tau>
    %     <parasolCells_k>7.</parasolCells_k>
    %     <amacrinCellsTemporalCutFrequency>2.0e+00</amacrinCellsTemporalCutFrequency>
    %     <V0CompressionParameter>9.5e-01</V0CompressionParameter>
    %     <localAdaptintegration_tau>0.</localAdaptintegration_tau>
    %     <localAdaptintegration_k>7.</localAdaptintegration_k>
    % </IPLmagno>
    % </opencv_storage>
    % ```
    %
    % Here are some words about all those parameters, tweak them as you wish
    % to amplify or moderate retina effects (contours enforcement, halos
    % effects, motion sensitivity, motion blurring, etc.)
    %
    % #### Basic parameters
    %
    % The simplest parameters are as follows:
    %
    % * __ColorMode__: let the retina process color information (if 1) or gray
    %   scale images (if 0). In that last case, only the first channels of the
    %   input will be processed.
    % * __NormaliseOutput__: each channel has such parameter: if the value is
    %   set to 1, then the considered channel's output is rescaled between 0
    %   and 255. Be aware at this case of the Magnocellular output level
    %   (motion/transient channel detection). Residual noise will also be
    %   rescaled!
    %
    % **Note**: using color requires color channels multiplexing/demultipexing
    % which also demands more processing. You can expect much faster processing
    % using gray levels: it would require around 30 product per pixel for all
    % of the retina processes and it has recently been parallelized for
    % multicore architectures.
    %
    % #### Photo-receptors parameters
    %
    % The following parameters act on the entry point of the retina
    % (photo-receptors) and has impact on all of the following processes.
    % These sensors are low pass spatio-temporal filters that smooth temporal
    % and spatial data and also adjust their sensitivity to local luminance,
    % thus, leads to improving details extraction and high frequency noise
    % canceling.
    %
    % * __PhotoreceptorsLocalAdaptationSensitivity__ between 0 and 1. Values
    %   close to 1 allow high luminance log compression's effect at the
    %   photo-receptors level. Values closer to 0 provide a more linear
    %   sensitivity. Increased alone, it can burn the *Parvo (details channel)*
    %   output image. If adjusted in collaboration with
    %   `GanglionCellsSensitivity`, images can be very contrasted whatever the
    %   local luminance there is at the cost of a naturalness decrease.
    % * __PhotoreceptorsTemporalConstant__ this setups the temporal constant
    %   of the low pass filter effect at the entry of the retina. High value
    %   leads to strong temporal smoothing effect: moving objects are blurred
    %   and can disappear while static object are favored. But when starting
    %   the retina processing, stable state is reached later.
    % * __PhotoreceptorsSpatialConstant__ specifies the spatial constant
    %   related to photo-receptors' lowpass filter's effect. Those parameters
    %   specify the minimum value of the spatial signal period allowed in what
    %   follows. Typically, this filter should cut high frequency noise. On
    %   the other hand, a 0 value cuts none of the noise while higher values
    %   start to cut high spatial frequencies, and progressively lower
    %   frequencies... Be aware to not go to high levels if you want to see
    %   some details of the input images! A good compromise for color images
    %   is a 0.53 value since such choice won't affect too much the color
    %   spectrum. Higher values would lead to gray and blurred output images.
    %
    % #### Horizontal cells parameters
    %
    % This parameter set tunes the neural network connected to the
    % photo-receptors, the horizontal cells. It modulates photo-receptors
    % sensitivity and completes the processing for final spectral whitening
    % (part of the spatial band pass effect thus favoring visual details
    % enhancement).
    %
    % * __HorizontalCellsGain__ here is a critical parameter! If you are not
    %   interested with the mean luminance and want just to focus on details
    %   enhancement, then, set this parameterto zero. However, if you want to
    %   keep some environment luminance's data, let some low spatial
    %   frequencies pass into the system and set a higher value (<1).
    % * __HCellsTemporalConstant__ similar to photo-receptors, this parameter
    %   acts on the temporal constant of a low pass temporal filter that
    %   smooths input data. Here, a high value generates a high retina after
    %   effect while a lower value makes the retina more reactive. This value
    %   should be lower than `PhotoreceptorsTemporalConstant` to limit strong
    %   retina after effects.
    % * __HCellsSpatialConstant__ is the spatial constant of these cells
    %   filter's low pass one. It specifies the lowest spatial frequency
    %   allowed in what follows. Visually, a high value leads to very low
    %   spatial frequencies processing and leads to salient halo effects.
    %   Lower values reduce this effect but has the limit of not go lower than
    %   the value of `PhotoreceptorsSpatialConstant`. Those 2 parameters
    %   actually specify the spatial band-pass of the retina.
    %
    % **NOTE**: Once the processing managed by the previous parameters is done,
    % input data is cleaned from noise and luminance is already partly
    % enhanced. The following parameters act on the last processing stages of
    % the two outing retina signals.
    %
    % #### Parvo (details channel) dedicated parameter
    %
    % * __GanglionCellsSensitivity__ specifies the strength of the final local
    %   adaptation occurring at the output of this details' dedicated channel.
    %   Parameter values remain between 0 and 1. Low value tend to give a
    %   linear response while higher values enforce the remaining low
    %   contrasted areas.
    %
    % **Note**: this parameter can correct eventual burned images by favoring
    % low energetic details of the visual scene, even in bright areas.
    %
    % #### IPL Magno (motion/transient channel) parameters
    %
    % Once image's information are cleaned, this channel acts as a high pass
    % temporal filter that selects only the signals related to transient
    % signals (events, motion, etc.). A low pass spatial filter smooths
    % extracted transient data while a final logarithmic compression enhances
    % low transient events thus enhancing event sensitivity.
    %
    % * __ParasolCellsBeta__ generally set to zero, can be considered as an
    %   amplifier gain at the entry point of this processing stage. Generally
    %   set to 0.
    % * __ParasolCellsTau__ the temporal smoothing effect that can be added
    % * __ParasolCellsK__ the spatial constant of the spatial filtering
    %   effect, set it at a high value to favor low spatial frequency signals
    %   that are lower subject for residual noise.
    % * __AmacrinCellsTemporalCutFrequency__ specifies the temporal constant
    %   of the high pass filter. High values let slow transient events to be
    %   selected.
    % * __V0CompressionParameter__ specifies the strength of the log
    %   compression. Similar behaviors to previous description but here
    %   enforces sensitivity of transient events.
    % * __LocalAdaptintegrationTau__ generally set to 0, has no real use
    %   actually in here.
    % * __LocalAdaptintegrationK__ specifies the size of the area on which
    %   local adaptation is performed. Low values lead to short range local
    %   adaptation (higher sensitivity to noise), high values secure log
    %   compression.
    %
    % ### Demos and experiments
    %
    % #### First time experiments
    %
    % Here are some code snippets to shortly show how to use Retina with
    % default parameters (with halo effects). Next section redirects to more
    % complete demos provided with the main retina class.
    %
    % Here is presented how to process a webcam stream with the following
    % steps:
    %
    % - load a frist input image to get its size
    % - allocate a retina instance with appropriate input size
    % - loop over grabbed frames:
    %   - grab a new frame
    %   - run on a frame
    %   - call the two output getters
    %   - display retina outputs
    %
    % See `retina_demo.m` MATLAB sample.
    %
    % #### More complete demos
    %
    % Note: Complementary to the following examples, have a look at the Retina
    % tutorial in the tutorial/contrib section for complementary explanations.
    %
    % ## References
    % [Benoit2010]:
    % > Benoit A., Caplier A., Durette B., Herault, J.; "Using Human Visual
    % > System Modeling For Bio-inspired Low Level Image Processing", Elsevier,
    % > Computer Vision and Image Understanding 114 (2010), pp. 758-773,
    % > [DOI](http://dx.doi.org/10.1016/j.cviu.2010.01.011)
    %
    % [Benoit2014]:
    % > Strat S.T., Benoit A., Lambert P.; "Retina enhanced bag of words
    % > descriptors for video classification". In Proceedings of the 22nd
    % > European Signal Processing Conference (EUSIPCO), 2014, pp. 1307-1311.
    %
    % [Herault2010]:
    % > Jeanny Herault; "Vision: Images, Signals and Neural Networks: Models
    % > of Neural Processing in Visual Perception" (Progress in Neural
    % > Processing), ISBN: 9814273686. WAPI (Tower ID): 113266891.
    %
    % [Chaix2007]:
    % > B. Chaix de Lavarene, D. Alleysson, B. Durette, J. Herault (2007).
    % > "Efficient demosaicing through recursive filtering". IEEE
    % > International Conference on Image Processing ICIP 2007.
    %
    % [Meylan2007]:
    % > Meylan L., Alleysson D., and Susstrunk S., "A Model of Retinal Local
    % > Adaptation for the Tone Mapping of Color Filter Array Images", Journal
    % > of Optical Society of America, A, Vol. 24, N. 9, September, 1st, 2007,
    % > pp. 2807-2816
    %
    % See also: cv.RetinaFastToneMapping, cv.TransientAreasSegmentationModule
    %

    properties (SetAccess = private)
        % Object ID
        id
    end

    methods
        function this = Retina(inputSize, varargin)
            %RETINA  Constructor from standardized interface to create a Retina instance
            %
            %     obj = cv.Retina(inputSize)
            %     obj = cv.Retina(inputSize, 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __inputSize__ the input frame size `[w,h]`.
            %
            % ## Options
            % * __ColorMode__ the chosen processing mode: with or without
            %   color processing. default true
            % * __ColorSamplingMethod__ specifies which kind of color sampling
            %   will be used, default 'Bayer'. One of:
            %   * __Random__ each pixel position is either R, G or B in a
            %     random choice
            %   * __Diagonal__ color sampling is RGBRGBRGB...,
            %     line 2 BRGBRGBRG..., line 3 GBRGBRGBR...
            %   * __Bayer__ standard bayer sampling
            % * __UseRetinaLogSampling__ activate retina log sampling. If true,
            %   the 2 following parameters can be used. default false
            % * __ReductionFactor__ only useful if param
            %   `UseRetinaLogSampling=true`, specifies the reduction factor of
            %   the output frame (as the center (fovea) is high resolution and
            %   corners can be underscaled, then a reduction of the output is
            %   allowed without precision leak). default 1.0
            % * __SamplingStrength__ only useful if param
            %   `UseRetinaLogSampling=true`, specifies the strength of the log
            %   scale that is applied. default 10.0
            %
            % See also: cv.Retina.run
            %
            this.id = Retina_(0, 'new', inputSize, varargin{:});
        end

        function delete(this)
            %DELETE  Destructor
            %
            %     obj.delete()
            %
            % See also: cv.Retina
            %
            if isempty(this.id), return; end
            Retina_(this.id, 'delete');
        end
    end

    %% Algorithm
    methods
        function clear(this)
            %CLEAR  Clears the algorithm state
            %
            %     obj.clear()
            %
            % See also: cv.Retina.empty, cv.Retina.load
            %
            Retina_(this.id, 'clear');
        end

        function b = empty(this)
            %EMPTY  Checks if detector object is empty
            %
            %     b = obj.empty()
            %
            % ## Output
            % * __b__ Returns true if the detector object is empty (e.g in the
            %   very beginning or after unsuccessful read).
            %
            % See also: cv.Retina.clear, cv.Retina.load
            %
            b = Retina_(this.id, 'empty');
        end

        function save(this, filename)
            %SAVE  Saves the algorithm parameters to a file
            %
            %     obj.save(filename)
            %
            % ## Input
            % * __filename__ Name of the file to save to.
            %
            % This method stores the algorithm parameters in the specified
            % XML or YAML file.
            %
            % See also: cv.Retina.load
            %
            Retina_(this.id, 'save', filename);
        end

        function load(this, fname_or_str, varargin)
            %LOAD  Loads algorithm from a file or a string
            %
            %     obj.load(fname)
            %     obj.load(str, 'FromString',true)
            %     obj.load(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __fname__ Name of the file to read.
            % * __str__ String containing the serialized model you want to
            %   load.
            %
            % ## Options
            % * __ObjName__ The optional name of the node to read (if empty,
            %   the first top-level node will be used). default empty
            % * __FromString__ Logical flag to indicate whether the input is a
            %   filename or a string containing the serialized model.
            %   default false
            %
            % This method reads algorithm parameters from the specified XML or
            % YAML file (either from disk or serialized string). The previous
            % algorithm state is discarded.
            %
            % See also: cv.Retina.save
            %
            Retina_(this.id, 'load', fname_or_str, varargin{:});
        end

        function name = getDefaultName(this)
            %GETDEFAULTNAME  Returns the algorithm string identifier
            %
            %     name = obj.getDefaultName()
            %
            % ## Output
            % * __name__ This string is used as top level XML/YML node tag
            %   when the object is saved to a file or string.
            %
            % See also: cv.Retina.save, cv.Retina.load
            %
            name = Retina_(this.id, 'getDefaultName');
        end
    end

    %% Retina
    methods
        function sz = getInputSize(this)
            %GETINPUTSIZE  Retreive retina input buffer size
            %
            %     sz = obj.getInputSize()
            %
            % ## Output
            % * __sz__ the retina input buffer size
            %
            % See also: cv.Retina.getOutputSize
            %
            sz = Retina_(this.id, 'getInputSize');
        end

        function sz = getOutputSize(this)
            %GETOUTPUTSIZE  Retreive retina output buffer size that can be different from the input if a spatial log transformation is applied
            %
            %     sz = obj.getOutputSize()
            %
            % ## Output
            % * __sz__ the retina output buffer size
            %
            % See also: cv.Retina.getInputSize
            %
            sz = Retina_(this.id, 'getOutputSize');
        end

        function setup(this, retinaParameterFile, varargin)
            %SETUP  Try to open an XML retina parameters file to adjust current retina instance setup
            %
            %     obj.setup(retinaParameterFile)
            %     obj.setup(retinaParameterFile, 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __retinaParameterFile__ the parameters filename.
            %
            % ## Options
            % * __ApplyDefaultSetupOnFailure__ set to true if an error must be
            %   thrown on error. default true
            %
            % If the xml file does not exist, then default setup is applied.
            % Warning: Exceptions are thrown if read XML file is not valid.
            %
            % You can retrieve the current parameters structure using the
            % method cv.Retina.getParameters and update it before running
            % method cv.Retina.setupParameters.
            %
            % See also: cv.Retina.setupParameters, cv.Retina.write
            %
            Retina_(this.id, 'setup', retinaParameterFile, varargin{:});
        end

        function setupParameters(this, varargin)
            %SETUPPARAMETERS  Pass retina parameters to adjust current retina instance setup
            %
            %     obj.setupParameters('OptionName',optionValue, ...)
            %
            % ## Options
            % * __OPLandIplParvo__ Outer Plexiform Layer (OPL) and Inner
            %   Plexiform Layer Parvocellular (IplParvo) parameters.
            % * __IplMagno__ Inner Plexiform Layer Magnocellular channel
            %   (IplMagno).
            %
            % ### `OPLandIplParvo` options
            % See cv.Retina.setupOPLandIPLParvoChannel options.
            %
            % ### `IplMagno` options
            % See cv.Retina.setupIPLMagnoChannel options.
            %
            % Specify retina model parameters structure, with the new target
            % configuration.
            %
            % For better clarity, check explenations on the comments of
            % methods `setupOPLandIPLParvoChannel` and `setupIPLMagnoChannel`.
            %
            % Here is the default configuration file of the retina module. It
            % gives results such as the first retina output shown on the top
            % of this page.
            %
            % ```xml
            % <?xml version="1.0"?>
            % <opencv_storage>
            % <OPLandIPLparvo>
            %     <colorMode>1</colorMode>
            %     <normaliseOutput>1</normaliseOutput>
            %     <photoreceptorsLocalAdaptationSensitivity>7.5e-01</photoreceptorsLocalAdaptationSensitivity>
            %     <photoreceptorsTemporalConstant>9.0e-01</photoreceptorsTemporalConstant>
            %     <photoreceptorsSpatialConstant>5.3e-01</photoreceptorsSpatialConstant>
            %     <horizontalCellsGain>0.01</horizontalCellsGain>
            %     <hcellsTemporalConstant>0.5</hcellsTemporalConstant>
            %     <hcellsSpatialConstant>7.</hcellsSpatialConstant>
            %     <ganglionCellsSensitivity>7.5e-01</ganglionCellsSensitivity>
            % </OPLandIPLparvo>
            % <IPLmagno>
            %     <normaliseOutput>1</normaliseOutput>
            %     <parasolCells_beta>0.</parasolCells_beta>
            %     <parasolCells_tau>0.</parasolCells_tau>
            %     <parasolCells_k>7.</parasolCells_k>
            %     <amacrinCellsTemporalCutFrequency>2.0e+00</amacrinCellsTemporalCutFrequency>
            %     <V0CompressionParameter>9.5e-01</V0CompressionParameter>
            %     <localAdaptintegration_tau>0.</localAdaptintegration_tau>
            %     <localAdaptintegration_k>7.</localAdaptintegration_k>
            % </IPLmagno>
            % </opencv_storage>
            % ```
            %
            % Here is the "realistic" setup used to obtain the second retina
            % output shown on the top of this page.
            %
            % ```xml
            % <?xml version="1.0"?>
            % <opencv_storage>
            % <OPLandIPLparvo>
            %     <colorMode>1</colorMode>
            %     <normaliseOutput>1</normaliseOutput>
            %     <photoreceptorsLocalAdaptationSensitivity>8.9e-01</photoreceptorsLocalAdaptationSensitivity>
            %     <photoreceptorsTemporalConstant>9.0e-01</photoreceptorsTemporalConstant>
            %     <photoreceptorsSpatialConstant>5.3e-01</photoreceptorsSpatialConstant>
            %     <horizontalCellsGain>0.3</horizontalCellsGain>
            %     <hcellsTemporalConstant>0.5</hcellsTemporalConstant>
            %     <hcellsSpatialConstant>7.</hcellsSpatialConstant>
            %     <ganglionCellsSensitivity>8.9e-01</ganglionCellsSensitivity>
            % </OPLandIPLparvo>
            % <IPLmagno>
            %     <normaliseOutput>1</normaliseOutput>
            %     <parasolCells_beta>0.</parasolCells_beta>
            %     <parasolCells_tau>0.</parasolCells_tau>
            %     <parasolCells_k>7.</parasolCells_k>
            %     <amacrinCellsTemporalCutFrequency>2.0e+00</amacrinCellsTemporalCutFrequency>
            %     <V0CompressionParameter>9.5e-01</V0CompressionParameter>
            %     <localAdaptintegration_tau>0.</localAdaptintegration_tau>
            %     <localAdaptintegration_k>7.</localAdaptintegration_k>
            % </IPLmagno>
            % </opencv_storage>
            % ```
            %
            % See also: cv.Retina.setup, cv.Retina.setupOPLandIPLParvoChannel,
            %  cv.Retina.setupIPLMagnoChannel
            %
            Retina_(this.id, 'setupParameters', varargin{:});
        end

        function setupOPLandIPLParvoChannel(this, varargin)
            %SETUPOPLANDIPLPARVOCHANNEL  Setup the OPL and IPL parvo channels (see biologocal model)
            %
            %     obj.setupOPLandIPLParvoChannel('OptionName',optionValue, ...)
            %
            % ## Options
            % * __ColorMode__ specifies if (true) color is processed of not
            %   (false) to then processing gray level image. default true
            % * __NormaliseOutput__ specifies if (true) output is rescaled
            %   between 0 and 255 of not (false). default true
            % * __PhotoreceptorsLocalAdaptationSensitivity__ the photoreceptors
            %   sensitivity renage is 0-1 (more log compression effect when
            %   value increases). default 0.7
            % * __PhotoreceptorsTemporalConstant__ the time constant of the
            %   first order low pass filter of the photoreceptors, use it to
            %   cut high temporal frequencies (noise or fast motion), unit is
            %   frames, typical value is 1 frame. default 0.5
            % * __PhotoreceptorsSpatialConstant__ the spatial constant of the
            %   first order low pass filter of the photoreceptors, use it to
            %   cut high spatial frequencies (noise or thick contours), unit
            %   is pixels, typical value is 1 pixel. default 0.53
            % * __HorizontalCellsGain__ gain of the horizontal cells network,
            %   if 0, then the mean value of the output is zero, if the
            %   parameter is near 1, then, the luminance is not filtered and
            %   is still reachable at the output, typicall value is 0.
            %   default 0.0
            % * __HCellsTemporalConstant__ the time constant of the first
            %   order low pass filter of the horizontal cells, use it to cut
            %   low temporal frequencies (local luminance variations), unit is
            %   frames, typical value is 1 frame, as the photoreceptors.
            %   default 1.0
            % * __HCellsSpatialConstant__ the spatial constant of the first
            %   order low pass filter of the horizontal cells, use it to cut
            %   low spatial frequencies (local luminance), unit is pixels,
            %   typical value is 5 pixel, this value is also used for local
            %   contrast computing when computing the local contrast
            %   adaptation at the ganglion cells level (Inner Plexiform Layer
            %   parvocellular channel model). default 7.0
            % * __GanglionCellsSensitivity__ the compression strengh of the
            %   ganglion cells local adaptation output, set a value between
            %   0.6 and 1 for best results, a high value increases more the
            %   low value sensitivity and the output saturates faster,
            %   recommended value: 0.7. default 0.7
            %
            % OPL is referred as Outer Plexiform Layer of the retina, it
            % allows the spatio-temporal filtering which withens the spectrum
            % and reduces spatio-temporal noise while attenuating global
            % luminance (low frequency energy) IPL parvo is the OPL next
            % processing stage, it refers to a part of the Inner Plexiform
            % layer of the retina, it allows high contours sensitivity in
            % foveal vision. See reference papers for more information. For
            % more information, please have a look at the paper [Benoit2010].
            %
            % See also: cv.Retina.setupIPLMagnoChannel
            %
            Retina_(this.id, 'setupOPLandIPLParvoChannel', varargin{:});
        end

        function setupIPLMagnoChannel(this, varargin)
            %SETUPIPLMAGNOCHANNEL  Set parameters values for the Inner Plexiform Layer (IPL) magnocellular channel
            %
            %     obj.setupIPLMagnoChannel('OptionName',optionValue, ...)
            %
            % ## Options
            % * __NormaliseOutput__ specifies if (true) output is rescaled
            %   between 0 and 255 of not (false). default true
            % * __ParasolCellsBeta__ the low pass filter gain used for local
            %   contrast adaptation at the IPL level of the retina (for
            %   ganglion cells local adaptation), typical value is 0.
            %   default 0.0
            % * __ParasolCellsTau__ the low pass filter time constant used for
            %   local contrast adaptation at the IPL level of the retina (for
            %   ganglion cells local adaptation), unit is frame, typical value
            %   is 0 (immediate response). default 0.0
            % * __ParasolCellsK__ the low pass filter spatial constant used
            %   for local contrast adaptation at the IPL level of the retina
            %   (for ganglion cells local adaptation), unit is pixels, typical
            %   value is 5. default 7.0
            % * __AmacrinCellsTemporalCutFrequency__ the time constant of the
            %   first order high pass fiter of the magnocellular way (motion
            %   information channel), unit is frames, typical value is 1.2.
            %   default 1.2
            % * __V0CompressionParameter__ the compression strengh of the
            %   ganglion cells local adaptation output, set a value between
            %   0.6 and 1 for best results, a high value increases more the
            %   low value sensitivity and the output saturates faster,
            %   recommended value: 0.95. default 0.95
            % * __LocalAdaptintegrationTau__ specifies the temporal constant
            %   of the low pas filter involved in the computation of the local
            %   "motion mean" for the local adaptation computation. default 0.0
            % * __LocalAdaptintegrationK__ specifies the spatial constant of
            %   the low pas filter involved in the computation of the local
            %   "motion mean" for the local adaptation computation. default 7.0
            %
            % This channel processes signals output from OPL processing stage
            % in peripheral vision, it allows motion information enhancement.
            % It is decorrelated from the details channel. See reference
            % papers for more details.
            %
            % See also: cv.Retina.setupOPLandIPLParvoChannel
            %
            Retina_(this.id, 'setupIPLMagnoChannel', varargin{:});
        end

        function params = getParameters(this)
            %GETPARAMETERS  Retrieve the current retina parameters values in a structure
            %
            %     params = obj.getParameters()
            %
            % ## Output
            % * __params__ the current parameters setup.
            %
            % See also: cv.Retina.setupParameters, cv.Retina.printSetup
            %
            params = Retina_(this.id, 'getParameters');
        end

        function str = printSetup(this)
            %PRINTSETUP  Outputs a string showing the used parameters setup
            %
            %     str = obj.printSetup()
            %
            % ## Output
            % * __str__ a string which contains formatted parameters
            %   information.
            %
            % See also: cv.Retina.getParameters
            %
            str = Retina_(this.id, 'printSetup');
        end

        function varargout = write(this, fs)
            %WRITE  Write xml/yml formatted parameters information
            %
            %     obj.write(fs)
            %     str = obj.write(fs)
            %
            % ## Input
            % * __fs__ the filename of the xml file that will be open and
            %   writen with formatted parameters information.
            %
            % ## Output
            % * __str__ optional output. If requested, the parameters are
            %   persisted to a string in memory instead of writing to disk.
            %
            % See also: cv.Retina.setup
            %
            [varargout{1:nargout}] = Retina_(this.id, 'write', fs);
        end

        function run(this, inputImage)
            %RUN  Method which allows retina to be applied on an input image
            %
            %     obj.run(inputImage)
            %
            % ## Input
            % * __inputImage__ the input image to be processed, can be gray
            %   level or BGR coded in any format (from 8bit to 16bits).
            %
            % After run, encapsulated retina module is ready to deliver its
            % outputs using dedicated acccessors, see `getParvo` and
            % `getMagno` methods.
            %
            % See also: cv.Retina.getParvo, cv.Retina.getMagno
            %
            if true
                %HACK: temp fix to get around an issue when OpenCL is enabled
                val = cv.Utils.useOptimized();
                cv.Utils.setUseOptimized(false);
                cObj = onCleanup(@() cv.Utils.setUseOptimized(val));
            end
            Retina_(this.id, 'run', inputImage);
        end

        function outputToneMappedImage = applyFastToneMapping(this, inputImage)
            %APPLYFASTTONEMAPPING  Method which processes an image in the aim to correct its luminance correct backlight problems, enhance details in shadows
            %
            %     outputToneMappedImage = obj.applyFastToneMapping(inputImage)
            %
            % ## Input
            % * __inputImage__ the input image to process (should be coded in
            %   float format `single`: 1/3/4-channels, the 4th channel won't
            %   be considered).
            %
            % ## Output
            % * __outputToneMappedImage__ the output 8bit/channel tone mapped
            %   image (`uint8` with 1/3-channels format).
            %
            % This method is designed to perform High Dynamic Range image tone
            % mapping (compress >8bit/pixel images to 8bit/pixel). This is a
            % simplified version of the Retina Parvocellular model (simplified
            % version of the `run`/`getParvo` methods call) since it does not
            % include the spatio-temporal filter modelling the Outer Plexiform
            % Layer of the retina that performs spectral whitening and many
            % other stuff. However, it works great for tone mapping and in a
            % faster way.
            %
            % Check the demos and experiments section to see examples and the
            % way to perform tone mapping using the original retina model and
            % the method.
            %
            % See also: cv.RetinaFastToneMapping
            %
            outputToneMappedImage = Retina_(this.id, 'applyFastToneMapping', inputImage);
        end

        function parvo = getParvo(this)
            %GETPARVO  Accessor of the details channel of the retina (models foveal vision)
            %
            %     parvo = obj.getParvo()
            %
            % ## Output
            % * __parvo__ the output buffer, format can be:
            %   * a matrix, this output is rescaled for standard 8bits image
            %     processing use in OpenCV
            %   * RAW methods actually return a 1D matrix (encoding is
            %     `R1, R2, ..., Rn, G1, G2, ..., Gn, B1, B2, ..., Bn`), this
            %     output is the original retina filter model output, without
            %     any quantification or rescaling.
            %
            % Warning, `getParvoRAW` methods return buffers that are not
            % rescaled within range [0;255] while the non RAW method allows a
            % normalized matrix to be retrieved.
            %
            % See also: cv.Retina.getParvoRAW
            %
            parvo = Retina_(this.id, 'getParvo');
        end

        function parvo = getParvoRAW(this)
            %GETPARVORAW  Accessor of the details channel of the retina (models foveal vision)
            %
            %     parvo = obj.getParvoRAW()
            %
            % ## Output
            % * __parvo__ the output buffer, format can be:
            %   * a matrix, this output is rescaled for standard 8bits image
            %     processing use in OpenCV
            %   * RAW methods actually return a 1D matrix (encoding is
            %     `R1, R2, ..., Rn, G1, G2, ..., Gn, B1, B2, ..., Bn`), this
            %     output is the original retina filter model output, without
            %     any quantification or rescaling.
            %
            % Warning, `getParvoRAW` methods return buffers that are not
            % rescaled within range [0;255] while the non RAW method allows a
            % normalized matrix to be retrieved.
            %
            % See also: cv.Retina.getParvo
            %
            parvo = Retina_(this.id, 'getParvoRAW');
        end

        function magno = getMagno(this)
            %GETMAGNO  Accessor of the motion channel of the retina (models peripheral vision)
            %
            %     magno = obj.getMagno()
            %
            % ## Output
            % * __magno__ the output buffer, format can be:
            %   * a matrix, this output is rescaled for standard 8bits image
            %     processing use in OpenCV
            %   * RAW methods actually return a 1D matrix (encoding is
            %     `M1, M2, ..., Mn`), this output is the original retina filter
            %     model output, without any quantification or rescaling.
            %
            % Warning, `getMagnoRAW` methods return buffers that are not
            % rescaled within range [0;255] while the non RAW method allows a
            % normalized matrix to be retrieved.
            %
            % See also: cv.Retina.getMagnoRAW
            %
            magno = Retina_(this.id, 'getMagno');
        end

        function magno = getMagnoRAW(this)
            %GETMAGNORAW  Accessor of the motion channel of the retina (models peripheral vision)
            %
            %     magno = obj.getMagnoRAW()
            %
            % ## Output
            % * __magno__ the output buffer, format can be:
            %   * a matrix, this output is rescaled for standard 8bits image
            %     processing use in OpenCV
            %   * RAW methods actually return a 1D matrix (encoding is
            %     `M1, M2, ..., Mn`), this output is the original retina filter
            %     model output, without any quantification or rescaling.
            %
            % Warning, `getMagnoRAW` methods return buffers that are not
            % rescaled within range [0;255] while the non RAW method allows a
            % normalized matrix to be retrieved.
            %
            % See also: cv.Retina.getMagno
            %
            magno = Retina_(this.id, 'getMagnoRAW');
        end

        function setColorSaturation(this, varargin)
            %SETCOLORSATURATION  Activate color saturation as the final step of the color demultiplexing process
            %
            %     obj.setColorSaturation('OptionName',optionValue, ...)
            %
            % ## Options
            % * __SaturateColors__ boolean that activates color saturation (if
            %   true) or desactivate (if false). default true
            % * __ColorSaturationValue__ the saturation factor: a simple
            %   factor applied on the chrominance buffers. default 4.0
            %
            % This saturation is a sigmoide function applied to each channel
            % of the demultiplexed image.
            %
            % See also: cv.Retina.run
            %
            Retina_(this.id, 'setColorSaturation', varargin{:});
        end

        function clearBuffers(this)
            %CLEARBUFFERS  Clears all retina buffers
            %
            %     obj.clearBuffers()
            %
            % (equivalent to opening the eyes after a long period of eye
            % close) whatchout the temporal transition occuring just after
            % this method call.
            %
            % See also: cv.Retina.run
            %
            Retina_(this.id, 'clearBuffers');
        end

        function activateMovingContoursProcessing(this, activate)
            %ACTIVATEMOVINGCONTOURSPROCESSING  Activate/desactivate the Magnocellular pathway processing (motion information extraction)
            %
            %     obj.activateMovingContoursProcessing(activate)
            %
            % ## Input
            % * __activate__ true if Magnocellular output should be activated,
            %   false if not. If activated, the Magnocellular output can be
            %   retrieved using the cv.Retina.getMagno method.
            %
            % By default, it is activated.
            %
            % See also: cv.Retina.activateContoursProcessing
            %
            Retina_(this.id, 'activateMovingContoursProcessing', activate);
        end

        function activateContoursProcessing(this, activate)
            %ACTIVATECONTOURSPROCESSING  Activate/desactivate the Parvocellular pathway processing (contours information extraction)
            %
            %     obj.activateContoursProcessing(activate)
            %
            % ## Input
            % * __activate__ true if Parvocellular (contours information
            %   extraction) output should be activated, false if not. If
            %   activated, the Parvocellular output can be retrieved using the
            %   cv.Retina.getParvo method.
            %
            % By default, it is activated.
            %
            % See also: cv.Retina.activateMovingContoursProcessing
            %
            Retina_(this.id, 'activateContoursProcessing', activate);
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function prob_test082 ( )

%*****************************************************************************80
%
%% TEST082 tests GEOMETRIC_MEAN, *_SAMPLE, *_VARIANCE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 January 2009
%
%  Author:
%
%    John Burkardt
%
  nsample = 1000;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'PROB_TEST082\n' );
  fprintf ( 1, '  For the Geometric PDF:\n' );
  fprintf ( 1, '  GEOMETRIC_MEAN computes the mean;\n' );
  fprintf ( 1, '  GEOMETRIC_SAMPLE samples;\n' );
  fprintf ( 1, '  GEOMETRIC_VARIANCE computes the variance.\n' );

  a = 0.25;

  check = geometric_check ( a );

  if ( ~check );
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PROB_TEST082 - Fatal error!\n' );
    fprintf ( 1, '  The parameters are not legal.\n' );
    return
  end

  mean = geometric_mean ( a );
  variance = geometric_variance ( a );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PDF parameter A =             %14f\n', a );
  fprintf ( 1, '  PDF mean =                    %14f\n', mean );
  fprintf ( 1, '  PDF variance =                %14f\n', variance );

  for i = 1 : nsample
    [ x(i), seed ] = geometric_sample ( a, seed );
  end

  mean = i4vec_mean ( nsample, x );
  variance = i4vec_variance ( nsample, x );
  xmax = max ( x(1:nsample) );
  xmin = min ( x(1:nsample) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Sample size =     %6d\n', nsample );
  fprintf ( 1, '  Sample mean =     %14f\n', mean );
  fprintf ( 1, '  Sample variance = %14f\n', variance );
  fprintf ( 1, '  Sample maximum =  %6d\n', xmax );
  fprintf ( 1, '  Sample minimum =  %6d\n', xmin );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function p = adjustable_CPD(CPD)
% ADJUSTABLE_CPD Does this CPD have any adjustable params? (gaussian)
% p = adjustable_CPD(CPD)

p = ~CPD.clamped_mean | ~CPD.clamped_cov | ~CPD.clamped_weights;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function prob = survProbStdModel(t,b,time)
% survProbStdModel: Computes survival probability using the standard model

prob = ones(size(t));
time0 = [0;time];
dtime = diff(time0);

for jdx=1:length(time)
   if jdx < length(time)
      tmpidx = t > time0(jdx) & t <= time0(jdx+1);
   else
      tmpidx = t > time0(jdx);
   end
   H = 0;
   if (jdx>1)
      H = dtime(1:jdx-1)'*b(1:jdx-1);
   end
   H = H + (t(tmpidx) - time0(jdx))*b(jdx);
   prob(tmpidx) = exp(-H);
end

author: Sable

repo: mcbench-benchmarks


text:
function i4mat_write ( output_filename, m, n, table )

%*****************************************************************************80
%
%% I4MAT_WRITE writes an I4MAT file.
%
%  Discussion:
%
%    An I4MAT is an array of I4's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 August 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the output filename.
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, integer TABLE(M,N), the points.
%
%    Input, logical HEADER, is TRUE if the header is to be included.
%

%
%  Open the file.
%
  output_unit = fopen ( output_filename, 'wt' );

  if ( output_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4MAT_WRITE - Error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'I4MAT_WRITE - Error!' );
    return;
  end
%
%  Write the data.
%
  for j = 1 : n
    for i = 1 : m
      fprintf ( output_unit, '  %12d', round ( table(i,j) ) );
    end
    fprintf ( output_unit, '\n' );
  end
%
%  Close the file.
%
  fclose ( output_unit );

  return
end


text:
function pass = test_bernoulli(~)
%TEST_BERNOULLI   Test that the Bernoulli method doesn't crash.
pass = doesNotCrash(@() cheb.bernoulli(4));

end

function pass = doesNotCrash(fn)
try
    fn();
    pass = true;
catch ME %#ok<NASGU>
    pass = false;
end
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%% IPOPT PARFOR testing
clc
clear all

%Objective
fun = @(x) 20 + x(1)^2 + x(2)^2 - 10*(cos(2*pi*x(1)) + cos(2*pi*x(2)));
%Constraints
lb = [5*pi;-20*pi];
ub = [20*pi;-4*pi];
%Setup Options
opts = optiset('solver','ipopt');
Opt = opti('fun',fun,'bounds',lb,ub,'opts',opts);

%Generate a series of starting points
n = 1000;
X0 = [linspace(lb(1),lb(end),n); linspace(ub(1),ub(end),n)];

%Preallocate solution vector
sols = zeros(2,n);
solp = zeros(2,n);

%% run serial (8.17s)
tic
for i = 1:n
    sols(:,i) = solve(Opt,X0(:,i));
end
toc

%%
parpool(4);

%% run parfor (2.6s)
tic
parfor i = 1:n
    solp(:,i) = solve(Opt,X0(:,i));
end
toc

%% check results
err = norm(sols-solp)

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function ui = mrViewROI(command, ui, varargin);
% ROI Editing and Display functions for mrViewer.
%
% ui = mrViewROI(command, ui, [options]);
%
% [The format is a bit of a break from the usual way of
% specifying the ui first, b/c the ui is usually omitted
% and the active one found w/ mrViewGet, but the command
% is always necessary]
%
% This function incorporates several small commands, which interact
% between the ROI tools and the mrViewer GUI. The commands are:
%   'new': create a new empty ROI
%
%   'add': add coords to the ROI, according to the current edit mode
%
%   'remove': remove coords from the ROI, according to the curr edit mode
%
%   'delete': delete an ROI, prompting with a dialog if one isn't specified
%
%   'draw': re-draw rois according to current view settings
%
%   'edit': edit selected ROI properties
%
%   'select': set selected ROI. If an extra 3rd argument is passed, which
%             is either an index into the ROI list, or a name of an ROI,
%             will select that ROI. Otherwise will select based on the
%             ROI select popup in the ROI panel.
%
%   'restrict': restrict ROI coordinates to within the overlays currently
%             being displayed. (Overlays in which the 'hide' option is
%             selected won't affect this.) mrView('restrict', ui, 'all')
%             will cause all ROIs to be restricted. Otherwise, the selected
%             ROI only will be restricted.
%
%   'load':   load a mrVista ROI into the viewer. The Session GUI must be
%             open for this.
%
%   'save':   save selected ROI in mrVista format. Will save as an Inplane
%             or Volume ROI, depending on the ROI's 'viewType' property.
%             (Edit the ROI to change this.)
%
%   'mesh':   If a mesh is opened, get the ROI points drawn on the mesh
%             as a new ROI. An optional mapping argument can be specified,
%             which determines how to grab the mesh ROI: it can be
%             'layer1', in which case only gray layer 1 nodes are drawn in the
%             new ROI, or 'all', in which case all layers are mapped. If
%             this option is omitted, defaults to using the 'layerMapMode'
%             pref in mrmPreferences.
%
%	'meshgrow': grow an ROI along the cortical surface, using the viewer's
%			  overlay settings to determine the boundaries of the ROI. The
%			  ROI will be a contiguous patch on the cortical surface
%			  bounded by the overlay conditions. (Note that if no
%			  thresholding is specified, the code will warn you and not
%			  make an ROI.) By default the ROI is grown starting from the
%			  mesh cursor. However, an alternate start coordinate can be
%			  specified as an optional input argument.
%
%   'disk':   If a gray-matter segmentation is installed, create a disk
%             ROI along the gray matter. ('disc' also works.) An optional
%             size argument can be specified, which determines the radius
%             of the disc in mm. If omitted, prompts the user.
%
%   'gray':   Create ROI containing all gray matter nodes for the selected
%             segmentation. (If a segmentation is loaded.)
%
% If called with no arguments, instead of returning a ui struct,
% mrViewROI will return the currently-selected ROI in the viewer.
%
% ras 08/05
if notDefined('ui'),        ui = mrViewGet;                         end
if ishandle(ui),            ui = get(ui, 'UserData');               end
if notDefined('command'),   ui = ui.rois(ui.settings.roi); return;  end

if ~isfield(ui,'display'), error('Need an active display');         end

% ras 03/07: I used to have a 'switch' command here, but it got pretty
% long, so I broke each command into its own subfunction to make finding a
% sub-command easier. This m-file has a lot of subfunctions! I'll probably
% move the dedicated ROI stuff (roi*Coords) to their own m-files soon...

% allow british spelling of the 'disk' command :)
if isequal(lower(command), 'disc'), command = 'disk'; end

% this switchyard will call the relevant sub-command functions
subfun = sprintf('mrView%s%sROI', upper(command(1)), lower(command(2:end)));
% try
	ui = feval(subfun, ui, varargin);
% catch
% 	error('Unknown command.')
% end

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

% /--------------------------------------------------------------------/ %
function ui = mrViewNewROI(ui, varargin);
%%%%%Create a new, empty ROI (and select it)
roi = roiCreate(ui.mr, []);
roi.name = sprintf('ROI%i',length(ui.rois)+1);
roi.lineHandles = []; % for drawing ROIs
roi.prevCoords = []; % for undo
if isempty(ui.rois), ui.rois=roi; else, ui.rois(end+1)=roi; end
ui.settings.roi = length(ui.rois); % select it

% select the ROI
mrViewROI('select', ui, length(ui.rois));
return
% /--------------------------------------------------------------------/ %

% /--------------------------------------------------------------------/ %
function ui = mrViewAddROI(ui, varargin)
%%%%%Add to current ROI, using selected method
% If no rois exist, make a new one
if isempty(ui.rois), ui = mrViewROI('new',ui);    end

% get new coords using selected method
method = ui.settings.roiEditMode;
switch method
	case 1, coords = roiRectCoords(ui);
	case 2, coords = roiCircleCoords(ui);
	case 3, coords = roiCubeCoords(ui);
	case 4, coords = roiSphereCoords(ui);
	case 5, coords = roiLineCoords(ui);
	case 6, coords = roiPointCoords(ui);
	case 7, coords = roiGrowCoords(ui);
end

% round and remove redundant
coords = round(coords);
coords = intersectCols(coords, coords);

% add to ROI
roi = ui.rois(ui.settings.roi);
roi.coords = [roi.coords coords];
roi.coords = intersectCols(roi.coords, roi.coords); % rm redundant
roi.modified = datestr(clock);
ui.rois(ui.settings.roi) = roi;

% show changes
mrViewROI('draw',ui);

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

% /--------------------------------------------------------------------/ %
function ui = mrViewRemoveROI(ui, varargin);
%%%%%Remove points from current ROI, using selected method
% If no rois exist, exit quietly
if isempty(ui.rois), return;  end

% get coords to remove using selected method
method = ui.settings.roiEditMode;
switch method
	case 1, coords = roiRectCoords(ui);
	case 2, coords = roiCircleCoords(ui);
	case 3, coords = roiCubeCoords(ui);
	case 4, coords = roiSphereCoords(ui);
	case 5, coords = roiLineCoords(ui);
	case 6, coords = roiPointCoords(ui);
	case 7, coords = roiGrowCoords(ui);
end

% remove coords from ROI
roi = ui.rois(ui.settings.roi);
[ignore oldInd] = setdiff(roi.coords', coords', 'rows');
roi.coords = roi.coords(:,sort(oldInd)); % preserve voxel order
roi.modified = datestr(clock);
ui.rois(ui.settings.roi) = roi;

% show changes
mrViewROI('draw', ui);

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

% /--------------------------------------------------------------------/ %
function ui = mrViewDeleteROI(ui, varargin);
%%%%%Delete an ROI from the UI
varargin = unNestCell(varargin);
% (optional arg is list of ROIs to delete)
if isempty(varargin)
	delRoi = ui.settings.roi;
elseif isequal(lower(varargin{1}),'all')
	delRoi = 1:length(ui.rois);
else
	delRoi = varargin{1};
end

% if we don't have any ROIs and try to delete, don't sweat it
if isempty(ui.rois), return; end

% indices of non-deleted ROIs
newInd = setdiff(1:length(ui.rois), delRoi);

% figure out what the selected ROI index will be after we delete
if ismember(ui.settings.roi, delRoi)
	% selected ROI was deleted -- select last one in list
	selRoi = length(newInd);
else
	% selected ROI was not deleted -- but its index may have
	% changed. Figure out the new index.
	selRoi = ui.rois(ui.settings.roi).name; % can select by name
end

% now, delete -- keep only non-deleted ROIs
ui.rois = ui.rois(newInd);

% set a new selected ROI
mrViewROI('select', ui, selRoi);

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

% /--------------------------------------------------------------------/ %
function ui = mrViewLoadROI(ui);
%%%%%Load an ROI
% mrVista 2, back-compatible version: load from GUI
sessionGUI_loadROI;
return
% /--------------------------------------------------------------------/ %

% /--------------------------------------------------------------------/ %
function ui = mrViewSaveROI(ui, varargin);
%%%%%Save an ROI
varargin = unNestCell(varargin);
% for now, we save to a back-compatible mrVista 1 format.
if ~isempty(varargin) & isequal(varargin{1}, 'all')
	roiList = 1:length(ui.rois);
else
	roiList = ui.settings.roi;
end

for r = roiList
	mrViewSave(ui, r, 'roi');
end

% update the sessionGUI to include this ROI
sessionGUI_selectROIType;

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

% /--------------------------------------------------------------------/ %
function ui = mrViewRestrictROI(ui, varargin);
%%%%%Restrict ROI(s) to within overlays
% (If no overlays, don't restrict -- would just delete ROI)
if isempty(ui.overlays), return; end

varargin = unNestCell(varargin);
if ~isempty(varargin) & isequal(lower(varargin{1}), 'all')
	roiList = 1:length(ui.rois);
else
	roiList = ui.settings.roi;
end

for r = roiList
	roi = ui.rois(r);
	hidden = [ui.overlays.hide];
	roi.coords = mrViewRestrict(ui, roi.coords, find(~hidden));
	roi.modified = datestr(clock);
	ui.rois(r) = roi;
end

% show changes
ui = mrViewROI('draw',ui);

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

% /--------------------------------------------------------------------/ %
function ui = mrViewEditROI(ui, varargin);
%%%%%Edit/view ROI traits, such as name/color
if isempty(ui.rois)
	myWarnDlg('No ROIs Created Yet.');
	return
end
roi = ui.rois(ui.settings.roi);
roi = roiEdit(roi);
ui.rois(ui.settings.roi) = roi;

% show changes (updating popup in this call)
try, ui = mrViewROI('draw', ui); end

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

% /--------------------------------------------------------------------/ %
function ui = mrViewDrawROI(ui, varargin);
%%%%%Draw rois on display
% delete any existing ROI lines
for ax = ui.display.axes
	delete(findobj('Parent', ax, 'Type', 'line', 'Tag', 'ROI'));
end

% get list of ROIs to draw, based on ROI view mode
if ui.settings.roiViewMode==1 | isempty(ui.rois)
	% hide ROIs, or no ROIs -- return quietly
	for r = 1:length(ui.rois)
		ui.rois(r).lineHandles = [];
	end
	set(ui.fig, 'UserData' ,ui);
	return
elseif ui.settings.roiViewMode==2 % show selected
	roiList = ui.settings.roi;
else,                              % show all
	roiList = 1:length(ui.rois);
end

% draw all selected ROIs on displayed images
for i = 1:length(ui.display.images)
	for r = roiList
		% find ROI points which may lie in image
		C = ui.rois(r).coords;
		if ui.settings.space==1
			% pixel space, use ROI coords directly
			if isempty(C), pts = [];
			else, pts = C(1:2,C(3,:)==ui.display.slices(i));
			end
		else
			% get ROI coords which correspond to imaged
			% region
			imgCoords = round(mrViewGet(ui,'DataCoords',i));
			[ignore ind] = intersectCols(imgCoords,C);
			imgSz = size(ui.display.images{i}); imgSz = imgSz(1:2);
			[rows cols] = ind2sub(imgSz,ind);
			AX = [get(ui.display.axes(i), 'XLim') ...
				  get(ui.display.axes(i), 'YLim')]; % will have to account for axis bounds
			pts = [rows+AX(3)-1 cols+AX(1)-1]';
		end

		% draw points, taking preferences into account
		if ~isempty(pts)
			prefs.color = ui.rois(r).color;
			prefs.axesHandle = ui.display.axes(i);
			prefs.method = cellfind({'perimeter' 'filled' 'patches'}, ...
									ui.rois(r).fillMode);
			if r==ui.settings.roi
				prefs.lineStyle = '-';
				prefs.lineWidth = 1;
				prefs.color = 'w';
			else
				prefs.lineStyle = ':';
				prefs.lineWidth = 0.2;
			end
			axes(ui.display.axes(i));
			ui.rois(r).lineHandles = outline(pts,prefs);
			set(ui.rois(r).lineHandles, 'Tag', 'ROI');
		end
	end
end

% update changes in UI
set(ui.fig, 'UserData', ui);

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

% /--------------------------------------------------------------------/ %
function ui = mrViewSelectROI(ui, varargin);
%%%%%Set selected ROI
% parse the ROI specification
varargin = unNestCell(varargin);
if isempty(varargin)  % get from selected ROI
	selRoi = get(ui.controls.roiSelect, 'Value');
elseif isnumeric(varargin{1}) % ROI index
	selRoi = varargin{1};
elseif ischar(varargin{1}) % ROI name
	selRoi = cellfind({ui.rois.name}, varargin{1});
else
	error('Invalid ROI specification.')
end

ui.settings.roi = selRoi;

% update selected ROI popup if it exists
ui = updateROIPopup(ui);

try,   mrViewROI('draw', ui);   end

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

% /--------------------------------------------------------------------/ %
function ui = mrViewMeshROI(ui, varargin);
%%%%%Get Mesh ROI
msh = mrViewGet(ui, 'CurMesh');
if isempty(msh), myErrorDlg('No Mesh Loaded.'); end
varargin = unNestCell(varargin);
if isempty(varargin)    % get map method via mrmPreferences
	prefs = mrmPreferences;
	mapMethod = lower(prefs.layerMapMode);
else
	mapMethod = lower(varargin{1});
end

seg = mrViewGet(ui, 'CurSegmentation');
roi = roiFromMesh(seg, mapMethod, ui.mr);
ui = mrViewSet(ui, 'NewROI', roi);
try,   mrViewROI('draw', ui);   end
% /--------------------------------------------------------------------/ %

% /--------------------------------------------------------------------/ %
function ui = mrViewMeshgrowROI(ui, varargin);
%%%%%Grow Mesh ROI along cortical surface
msh = mrViewGet(ui, 'CurMesh');
if isempty(msh), myErrorDlg('No Mesh Loaded.'); end

% get xform between the data and the I|P|R (mesh) coordinates
ii = cellfind( {ui.spaces.name}, 'I|P|R' );
xform = ui.spaces(ii).xform;

% get start coord (in base coordinate space)
varargin = unNestCell(varargin);
if isempty(varargin)
	% xform cursor position into base space
	startCoord = mrmGet(msh, 'CursorPosition');
	startCoord = coordsXform(inv(xform), startCoord);
else
	startCoord = varargin{1};
end

% xform the segmentation's nodes/edges into the base coordinate space, too
seg = mrViewGet(ui, 'segmentation');
seg.nodes([2 1 3],:) = coordsXform( inv(xform), seg.nodes([2 1 3],:) )';
seg.nodes(1:3,:) = round(seg.nodes(1:3,:));

% get gray matter node corresponding to the start coord
startCoord = floor(startCoord);
startNode = find( seg.nodes(2,:)==startCoord(1) & ...
				  seg.nodes(1,:)==startCoord(2) & ...
				  seg.nodes(3,:)==startCoord(3) );
if isempty(startNode)
	myWarnDlg( sprintf('Node not found for start coordinate %s.', ...
					num2str(startCoord)) );
	return
end

%% get data mask based on overlays, thresholds
mask3D = mrViewGet(ui, 'OverlayMask');
sz = size(mask3D);

% index mask according to the nodes
mask = logical( zeros(1, size(seg.nodes, 2)) );
allNodes = seg.nodes;
keep = find( allNodes(1,:) > 0 & allNodes(2,:) < sz(1) & ...
			allNodes(2,:) > 0 & allNodes(1,:) < sz(2) & ...
			allNodes(3,:) > 0 & allNodes(3,:) < sz(3) );
allNodes = allNodes(:,keep);
I = sub2ind( size(mask3D), allNodes(1,:), allNodes(2,:), allNodes(3,:) );
mask(I) = mask3D(I);
% clear mask3D

% does the start node pass the criterion? If not, we shouldn't make an ROI
% with it:
if mask(startNode)==0
	myWarnDlg('The selected start point doesn''t pass the threshold.');
	return
end

%% initialize an empty ROI
roi = roiCreate(ui.mr);
roi.name = ['Mesh ROI ' num2str(startCoord)];
roi.comments = sprintf(['Created by %s ''Meshgrow'' grown from seed %s'], mfilename, ...
						num2str(startCoord));

%% main part: "grow" coords in ROI (recursive)
h = msgbox('Growing Mesh ROI', mfilename);
iNodes = growMeshROICoords(startNode, seg, mask);
close(h);

roi.coords = coordsXform( inv(xform), seg.nodes([2 1 3],iNodes) );

%% add ROI to viewer
ui = mrViewSet(ui, 'AddAndSelectROI', roi);

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

% /-----------------------------------------------------------------/ %
function nodes = growMeshROICoords(nodes, seg, mask);
% recursive function designed to grow mesh ROI coordinates.
% The logic here is: the nodes passed in are all voxels that will
% definitely be included in the mesh ROI. The code looks at neighboring
% voxels (according to the gray matter nodes/edges), evalautes them for
% inclusion (mask==1), and selects those which pass the test.
% If there are any of these neighbors, the function is called again, now
% looking at the neighbors of these neighbors.

%% find neighbors of the input nodes
neighbors = [];
for n = nodes
	nEdges = seg.nodes(4,n);
	startEdge = seg.nodes(5,n);
	neighbors = [neighbors seg.edges(startEdge:startEdge+nEdges-1)];
end

%% which neighbors are contained in the mask?
neighbors = neighbors(mask(neighbors)==1);

% prevent the recursion from bouncing back and forth:
% some of these neighbors may already be included in nodes. Remove them;
% they're obviously okay.
% (If I didn't do this, then if A points to B, and B points to A, and both
% satisfy the criteria, the function would keep getting called, first with
% A as the neigbor, then B, then A, etc...)
neighbors = setdiff(neighbors, nodes);

%% do any neighbors pass the test?
% if no, we're done. If yes, we recurse: call this function on the existing
% nodes plus the neighbors:
if ~isempty(neighbors)
	nodes = growMeshROICoords([nodes neighbors], seg, mask);
end

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

% /--------------------------------------------------------------------/ %
function ui = mrViewDiskROI(ui, varargin);
%%%%%Create Disk ROI
varargin = unNestCell(varargin);
if isempty(varargin)    % get disk size via dialog
	q = {'Size of Gray Matter Disk?'};
	resp = inputdlg(q, 'Make Disk ROI', 1, {'3'});
	radius = str2num(resp{1});
else
	radius = varargin{1};
end

roi = roiGrayDisk(mrViewGet(ui, 'CurSegmentation'), ui.mr, radius);
ui = mrViewSet(ui, 'AddAndSelectROI', roi);
ui = mrViewRecenter(ui, 'roi');
try,   mrViewROI('draw', ui);   end
% /--------------------------------------------------------------------/ %

% /--------------------------------------------------------------------/ %
function ui = mrViewGrayROI(ui, varargin);
%%%%%Create Gray Matter ROI
seg = mrViewGet(ui, 'CurSegmentation');
roi = roiGrayAll(seg, ui.mr);
ui = mrViewSet(ui, 'newROI', roi);

try,   mrViewROI('draw', ui);   end

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

% /--------------------------------------------------------------------/ %
function coords = roiRectCoords(ui);
% Get points from the user specifying a rectangle on the mrViewer display,
% and compute a 3xN list of coordinates which that input represents.
figure(ui.fig);
tmp = get(ui.panels.display,'BackgroundColor');
set(ui.panels.display,'BackgroundColor','y');
hmsg = mrMessage('Select Corners of Rectangle');
[X Y] = ginput(2);
X = round(X); Y = round(Y);
set(ui.panels.display,'BackgroundColor',tmp);
close(hmsg);

% figure out what image was clicked on
imgNum = find(ui.display.axes==gca);
if isempty(imgNum)
    warning('Didn''t click on mrViewer display. No points added.')
    coords = [];
    return
end

% get range of rows (Y) and columns (X) corresponding to rectangle
% (the corner of the image isn't always (0,0), because of
% diff't coordinate transforms. So, account for the viewed range.)
AX = axis;
rows = round([min(Y):max(Y)]); %  - AX(3) + 1
cols = round([min(X):max(X)]); %  - AX(1) + 1

% get coordinates (in base MR data) of selected points
imgCoords = mrViewGet(ui,'DataCoords',imgNum); % whole image
imgSz = size(ui.display.images{imgNum}); imgSz = imgSz(1:2);
I = reshape(1:prod(imgSz),imgSz); % index matrix

% get indices that correspond to the rectange, and return those coords
ind = I(rows,cols);
coords = imgCoords(:,ind(:));

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

% /--------------------------------------------------------------------/ %
function coords = roiCircleCoords(ui);
% Get points from the user specifying a circle on the mrViewer display,
% and compute a 3xN list of coordinates which that input represents.
figure(ui.fig);
tmp = get(ui.panels.display,'BackgroundColor');
set(ui.panels.display,'BackgroundColor','y');
hmsg = mrMessage('Select center of circle, then a point on circumference');
AX = axis;
[X Y] = ginput(2);
X = round(X-AX(1)+1);
Y = round(Y-AX(3)+1);
close(hmsg);
set(ui.panels.display,'BackgroundColor',tmp);

% figure out what image was clicked on
imgNum = find(ui.display.axes==gca);
if isempty(imgNum)
    warning('Didn''t click on mrViewer display. No points added.')
    coords = [];
    return
end

% create a matrix R specifying the distance of each pixel in the
% image from the first selected point (center of circle):
imgSz = size(ui.display.images{imgNum});
[xx,yy] = meshgrid(1:imgSz(2),1:imgSz(1));
dX = xx-X(1); dY = yy-Y(1);
R = sqrt(dX.*dX + dY.*dY);

% the points in the circle are points where R is less than the
% radius (the distance between the first and second points):
radius = sqrt(diff(X)^2 + diff(Y)^2);
inCircle = find(R<radius);

% get coordinates (in base MR data) within the circle
allCoords = mrViewGet(ui,'DataCoords',imgNum); % whole image
coords = allCoords(:,inCircle);

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

% /--------------------------------------------------------------------/ %
function coords = roiCubeCoords(ui);
% Get points from the user specifying a cube on the mrViewer display,
% and compute a 3xN list of coordinates which that input represents.
coords = [];
myWarnDlg('Sorry, Not Implemented Yet');
return
% /--------------------------------------------------------------------/ %

% /--------------------------------------------------------------------/ %
function coords = roiSphereCoords(ui);
% Get points from the user specifying the center, and a point on the surface, of a sphere
% in the mrViewer display, and compute a 3xN list of coordinates which that input represents.
figure(ui.fig);
tmp = get(ui.panels.display,'BackgroundColor');
set(ui.panels.display,'BackgroundColor','y');
hmsg = mrMessage('Select center of sphere, then a point on circumference');
AX = axis;
[X Y] = ginput(2);
X = round(X-AX(1)+1);
Y = round(Y-AX(3)+1);
close(hmsg);
set(ui.panels.display,'BackgroundColor',tmp);

% figure out what image was clicked on
imgNum = find(ui.display.axes==gca);
if isempty(imgNum)
    warning('Didn''t click on mrViewer display. No points added.')
    coords = [];
    return
end

% create a matrix R specifying the distance of each pixel in the
% image from the first selected point (center of sphere):
% TODO: xform grid to current view space
[xx yy zz] = meshgrid(1:ui.mr.dims(2), 1:ui.mr.dims(1), 1:ui.mr.dims(3));
dX = xx-X(1); dY = yy-Y(1);  dZ = zz-mrViewGet(ui, 'curslice');
R = sqrt(dX.*dX + dY.*dY + dZ.*dZ);

% the points in the circle are points where R is less than the
% radius (the distance between the first and second points):
radius = sqrt(diff(X)^2 + diff(Y)^2);
inSphere = find(R<radius);

% get coordinates (in base MR data) within the circle
% (This is a temp fix, will only work in pixel space)
coords = [yy(inSphere)'; xx(inSphere)'; zz(inSphere)'];

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

% /--------------------------------------------------------------------/ %
function coords = roiLineCoords(ui);
% Get points from the user specifying a line on the mrViewer display,
% and compute a 3xN list of coordinates which that input represents.
figure(ui.fig);
tmp = get(ui.panels.display,'BackgroundColor');
set(ui.panels.display,'BackgroundColor','y');
hmsg = mrMessage('Select endpoints of line');
pts = round(ginput(2));
close(hmsg);
set(ui.panels.display,'BackgroundColor',tmp);

% figure out what image was clicked on
imgNum = find(ui.display.axes==gca);
if isempty(imgNum)
    warning('Didn''t click on mrViewer display. No points added.')
    coords = [];
    return
end

% find the 2D image coords that lie within the line specified by
% the two points
[cols rows] = findLinePoints(pts(1,:),pts(2,:));

% account for axis bounds
AX = axis;
cols = cols - AX(1) + 1;
rows = rows - AX(3) + 1;

% get coordinates (in base MR data) within the line
imgCoords = mrViewGet(ui,'DataCoords',imgNum); % whole image
imgSz = size(ui.display.images{imgNum}); imgSz=imgSz(1:2);
ind = sub2ind(imgSz,rows,cols);
coords = imgCoords(:,ind);

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

% /--------------------------------------------------------------------/ %
function coords = roiPointCoords(ui);
% Get points from the user specifying individual points on the mrViewer
% display, and compute a 3xN list of coordinates which that input
% represents (in the base MR data).
figure(ui.fig);
tmp = get(ui.panels.display,'BackgroundColor');
set(ui.panels.display,'BackgroundColor','y');
hmsg = mrMessage('Left click to select point, right click to quit.');
X = []; Y = []; button = 1;
while button==1
    [newX newY button] = ginput(1);
    X = [X; newX]; Y = [Y; newY];
end
X = round(X); Y = round(Y);
close(hmsg);
set(ui.panels.display,'BackgroundColor',tmp);

% figure out what image was clicked on
imgNum = find(ui.display.axes==gca);
if isempty(imgNum)
    warning('Didn''t click on mrViewer display. No points added.')
    coords = [];
    return
end

% get coordinates (in base MR data) for the selected points
allCoords = mrViewGet(ui,'DataCoords',imgNum); % whole image
imgSz = size(ui.display.images{imgNum}); imgSz=imgSz(1:2);
ind = sub2ind(imgSz,Y,X);
coords = allCoords(:,ind);

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

% /--------------------------------------------------------------------/ %
function coords = roiGrowCoords(ui);
% Get a point from the user from which to grow a 'blob' of
% interconnected points in 3D,  and return a 3xN list of coordinates
% for those points (in the base MR data).
figure(ui.fig);
tmp = get(ui.panels.display, 'BackgroundColor');
set(ui.panels.display, 'BackgroundColor', 'y');
hmsg = mrMessage('Select point from which to grow.');
pt = round(ginput(1));
close(hmsg);
set(ui.panels.display, 'BackgroundColor', tmp);

% figure out what image was clicked on
imgNum = find(ui.display.axes==gca);
if isempty(imgNum)
    warning('Didn''t click on mrViewer display. No points added.')
    coords = [];
    return
end

% find seed location in the data coords
% (take axis bounds into account)
imgCoords = mrViewGet(ui, 'DataCoords', imgNum);
imgSz = size(ui.display.images{imgNum}); imgSz=imgSz([1 2]);
AX = axis;
ind = sub2ind(imgSz, round(pt(2)-AX(3)+1), round(pt(1)-AX(1)+1));
% ind = sub2ind(imgSz, pt(2), pt(1));
seed = round(imgCoords(:,ind));

% First-pass method of growing a blob efficiently:
% (1) Create an initial box,  centered on seed
% (2) Restrict the box to coords which pass overlay thresholds
% (3) Find contiguous blob which contains seed
% (4) Test if box contains whole blob: if not,  repeat 1-4 iteratively
hmsg = msgbox('Growing ROI...');
boxSize = 10; stepSize = 10;
coords = [];
while 1
    % (1) create box (get coords specifying box in data)
    ymin = max(1, seed(1)-boxSize);
    ymax = min(ui.mr.dims(1), seed(1)+boxSize);
    xmin = max(1, seed(2)-boxSize);
    xmax = min(ui.mr.dims(2), seed(2)+boxSize);
    zmin = max(1, seed(3)-boxSize);
    zmax = min(ui.mr.dims(3), seed(3)+boxSize);
    [X Y Z] = meshgrid(xmin:xmax, ymin:ymax, zmin:zmax);
    boxCoords = [Y(:) X(:) Z(:)]'; clear X Y Z
    mask = logical(zeros(ui.mr.dims(1:3)));

    % (2) restrict box to overlays
    hidden = [ui.overlays.hide];
    boxCoords = mrViewRestrict(ui, boxCoords, find(~hidden));
    ok = sub2ind(size(mask), boxCoords(1,:), boxCoords(2,:), boxCoords(3,:));
    mask(round(ok)) = 1;

    % (3) Find contiguous blob which contains seed
    L = bwlabeln(mask,  6); % integer label matrix of 6-connected blobs
    seedLabel = L(seed(1), seed(2), seed(3));
    if seedLabel==0,  return;   end;  % no data so quit
    blob = (L==seedLabel); % binary matrix where blob is 1
    [i1,  i2,  i3] = ind2sub(size(mask), find(blob>0.5));
    coords = [i1 i2 i3]';

    % (4) Test if the box contains the entire blob
    % (Or is otherwise flush with the bounds of the data)
    inBoundsY=(all(i1>ymin)|ymin==1) & (all(i1<ymax)|ymax==ui.mr.dims(1));
    inBoundsX=(all(i2>xmin)|xmin==1) & (all(i2<xmax)|xmax==ui.mr.dims(2));
    inBoundsZ=(all(i3>zmin)|zmin==1) & (all(i3<zmax)|zmax==ui.mr.dims(3));
    if inBoundsX & inBoundsY & inBoundsZ
        break;
    else
        boxSize = boxSize + stepSize;
    end
end
close(hmsg);

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

% /--------------------------------------------------------------------/ %
function ui = updateROIPopup(ui);
% update selected ROI popup if it exists
if checkfields(ui, 'controls', 'roiSelect') & ishandle(ui.controls.roiSelect)
    str = get(ui.controls.roiSelect, 'String');
    if isempty(ui.rois)
        str = {'(No ROIs)'};
    else
        str = {ui.rois.name};
	end
	val = max(1, ui.settings.roi); % if sel==0, select the '(No ROIs)' string
    set(ui.controls.roiSelect, 'String', str, 'Value', val);
end
return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [ ft ] = gsp_translate_old(G, f, i)
%GSP_TRANSLATE Generalized translation of the signal f to the node i
%   Usage: ft = gsp_translate(G, f, i);
%
%   Input parameters
%       G   : Graph
%       f   : Signal (column)
%       i   : Indices of vertex (int)
%   Output parameters
%       ft  : translate signal
%
%   This function translate the column vector f onto the node i. If f*
%   is a matrix, the translation will be done to each column.
%
%

fhat=gsp_gft(G,f);
nt = size(f,2);

ft = sqrt(G.N)*gsp_igft(G,fhat .* ...
    repmat(transpose(G.U(i,:)),1,nt));

end

repo: gspbox


text:
function [W] = TW2W(TW)
% Convert power from terawatts to watts.
% Chad A. Greene 2012
W = TW*1e+12;


text:
%  Figure 10.38      Feedback Control of Dynamic Systems, 5e
%                        Franklin, Powell, Emami
%
% fig10_38.m is a script to create the symmetric
% root locus for the yaw damper control design
clf;
fpw =[-10.0000         0         0         0         0         0;
    0.0730   -0.0558   -0.9968    0.0802    0.0415         0     ;
   -4.7500    0.5980   -0.1150   -0.0318         0         0;
    1.5300   -3.0500    0.3880   -0.4650         0         0;
         0         0    0.0805    1.0000         0         0;
         0         0    1.0000         0         0   -0.3333];
 gpw =[10;
     0;
     0;
     0;
     0;
     0];
hpw = [0         0    1.0000         0         0   -0.3333];
jpw = 0;
% state space matrices for the symmetric root locus
a=[fpw, 0*fpw;
   -hpw'*hpw,  -fpw'];
b=[gpw;0*gpw];
c=[0*gpw' gpw'];
rlocus(a,b,c,0)
v=[-10 10 -1.5 1.5];
axis(v);
grid;
title('Fig. 10.38 Symmetric root locus of lateral dynamics')


text:
function uid = root(tree)
% XMLTREE/ROOT Root Method
% FORMAT uid = root(tree)
%
% tree   - XMLTree object
% uid    - UID of the root element of tree
%__________________________________________________________________________
%
% Return the uid of the root element of the tree.
%__________________________________________________________________________
% Copyright (C) 2002-2008  http://www.artefact.tk/

% Guillaume Flandin
% $Id: root.m 4460 2011-09-05 14:52:16Z guillaume $

% Actually root is necessarily the element whos UID is 1, by
% construction. However, xml_parser should return a tree with a ROOT
% element with as many children as needed but only ONE *element* child
% who would be the real root (and this method should return the UID of
% this element).

uid = 1;

% Update:
% xml_parser has been modified (not as explained above) to handle the
% case when several nodes are at the top level of the XML.
% Example: <!-- beginning --><root>blah blah</root><!-- end -->
% Now root is the first element node of the tree.

% Look for the first element in the XML Tree
for i=1:length(tree)
    if strcmp(get(tree,i,'type'),'element')
        uid = i;
        break
    end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% RES = reconSpyrLevs(PYR,INDICES,LOFILT,BFILTS,EDGES,LEVS,BANDS)
%
% Recursive function for reconstructing levels of a steerable pyramid
% representation.  This is called by reconSpyr, and is not usually
% called directly.

% Eero Simoncelli, 6/96.

function res = reconSpyrLevs(pyr,pind,lofilt,bfilts,edges,levs,bands);

nbands = size(bfilts,2);
lo_ind = nbands+1;
res_sz = pind(1,:);

% Assume square filters:
bfiltsz =  round(sqrt(size(bfilts,1)));

if any(levs > 1)

  if  (size(pind,1) > lo_ind)
    nres = reconSpyrLevs( pyr(1+sum(prod(pind(1:lo_ind-1,:)')):size(pyr,1)),  ...
	pind(lo_ind:size(pind,1),:), ...
	lofilt, bfilts, edges, levs-1, bands);
  else
    nres = pyrBand(pyr,pind,lo_ind); 	% lowpass subband
  end

  res = upConv(nres, lofilt, edges, [2 2], [1 1], res_sz);

else

  res = zeros(res_sz);

end

if any(levs == 1)
  ind = 1;
  for b = 1:nbands
    if any(bands == b)
      bfilt = reshape(bfilts(:,b), bfiltsz, bfiltsz);
      upConv(reshape(pyr(ind:ind+prod(res_sz)-1), res_sz(1), res_sz(2)), ...
	     bfilt, edges, [1 1], [1 1], res_sz, res); % Destructively modify res
    end
    ind = ind + prod(res_sz);
  end
end


text:
function [U,sm,X,V,W] = cgsvd(A,L)
%CGSVD Compact generalized SVD of a matrix pair in regularization problems.
%
% sm = cgsvd(A,L)
% [U,sm,X,V] = cgsvd(A,L) ,  sm = [sigma,mu]
% [U,sm,X,V,W] = cgsvd(A,L) ,  sm = [sigma,mu]
%
% Computes the generalized SVD of the matrix pair (A,L). The dimensions of
% A and L must be such that [A;L] does not have fewer rows than columns.
%
% If m >= n >= p then the GSVD has the form:
%    [ A ] = [ U  0 ]*[ diag(sigma)      0    ]*inv(X)
%    [ L ]   [ 0  V ] [      0       eye(n-p) ]
%                     [  diag(mu)        0    ]
% where
%    U  is  m-by-n ,    sigma  is  p-by-1
%    V  is  p-by-p ,    mu     is  p-by-1
%    X  is  n-by-n .
%
% Otherwise the GSVD has a more complicated form (see manual for details).
%
% A possible fifth output argument returns W = inv(X).

% Reference: C. F. Van Loan, "Computing the CS and the generalized
% singular value decomposition", Numer. Math. 46 (1985), 479-491.

% Per Christian Hansen, IMM, March 17, 2008.

% Initialization.
[m,n] = size(A); [p,n1] = size(L);
if (n1 ~= n)
  error('No. columns in A and L must be the same')
end
if (m+p < n)
  error('Dimensions must satisfy m+p >= n')
end

% Call Matlab's GSVD routine.
[U,V,W,C,S] = gsvd(full(A),full(L),0);

if (m >= n)
  % The overdetermined or square case.
  sm = [diag(C(1:p,1:p)),diag(S(1:p,1:p))];
  if (nargout < 2)
    U = sm;
  else
    % Full decomposition.
    X = inv(W');
  end
else
  % The underdetermined case.
  sm = [diag(C(1:m+p-n,n-m+1:p)),diag(S(n-m+1:p,n-m+1:p))];
  if (nargout < 2)
    U = sm;
  else
    % Full decomposition.
    X = inv(W');
    X = X(:,n-m+1:n);
  end
end

if (nargout==5), W = W'; end

repo: REST