text:
function plot_optfit(a,time,timef,bootloops,maepi)
    % function plot_optfit(a,time,timef,bootloops,maepi);
    % ---------------------------------------------------
    % Plots Ncum observed vs. Ncum modeled for specified time windows
    %
    % Input variables:
    % a         : earthquake catalog (Complete in magnitude!)
    % time      : learning period fo fit Omori parameters
    % timef     : forecast period
    % bootloops : Number of bootstraps
    % maepi     : mainshock
    %
    % Samuel Neukomm
    % last update: 01.03.04

    warning off;

    [fMc] = calc_Mc(a, 1, 0.1)+0.2;
    l = a.Magnitude >= fMc;
    a = a.subset(l);

    if size(a,2) == 9
        date_matlab = datenum(a.Date.Year,a.Date.Month,a.Date.Day,a.Date.Hour,a.Date.Minute,zeros(size(a,1),1));
        date_main = datenum(floor(maepi(3)),maepi(4),maepi(5),maepi(8),maepi(9),0);
    else
        date_matlab = datenum(a.Date.Year,a.Date.Month,a.Date.Day,a.Date.Hour,a.Date.Minute,a(:,10));
        date_main = datenum(floor(maepi(3)),maepi(4),maepi(5),maepi(8),maepi(9),maepi(10));
    end
    time_aftershock = date_matlab-date_main;
    % Select biggest aftershock earliest in time, but more than 1 day after mainshock
    fDay = 1;
    ft_c=fDay/365; % Time not considered to find biggest aftershock
    vSel = (a.Date > maepi(:,3)+ft_c & a.Date<= maepi(:,3)+time/365);
    mCat = a.subset(vSel);
    vSel = mCat(:,6) == max(mCat(:,6));
    vBigAf = mCat(vSel,:);
    if length(mCat(:,1)) > 1
        vSel = vBigAf(:,3) == min(vBigAf(:,3));
        vBigAf = vBigAf(vSel,:);
    end

    if size(a,2) == 9
        date_biga = datenum(floor(vBigAf(3)),vBigAf(4),vBigAf(5),vBigAf(8),vBigAf(9),0);
    else
        date_biga = datenum(floor(vBigAf(3)),vBigAf(4),vBigAf(5),vBigAf(8),vBigAf(9),vBigAf(10));
    end
    fT1 = date_biga - date_main; % Time of big aftershock

    % Aftershock times
    l = time_aftershock(:) > 0;
    tas = time_aftershock(l);
    eqcatalogue = a.subset(l);

    % time_as: Learning period
    l = tas <= time;
    time_as = tas(l);
    time_as = sort(time_as);
    if length(time_as) < 100 % at least 100 events in learning period
        return
    end

    % Times up to the forecast time
    lf = tas <= time+timef ;
    time_asf= [tas(lf) ];
    time_asf=sort(time_asf);

    figure_w_normalized_uicontrolunits('Numbertitle','off','Name','Forecast aftershock occurence')
    hold on

    n = length(time_as);
    loopout = []; nCnt = 0;

    for j = 1:bootloops
        clear newtas
        randnr = ceil(rand(n,1)*n);
        i = (1:n)';
        newtas(i,:) = time_as(randnr(i),:); % bootstrap sample

        % Calculate fits of different models
        mRes = [];
        % Modified Omori law (pck)
        nMod = 1; [pval1, pval2, cval1, cval2, kval1, kval2, fAIC, fL] = bruteforceloglike_a2(newtas,fT1,nMod);
        mRes = [mRes; nMod, pval1, pval2, cval1, cval2, kval1, kval2, fAIC, fL];
        % MOL with secondary aftershock (pckk)
        nMod = 2; [pval1, pval2, cval1, cval2, kval1, kval2, fAIC, fL] = bruteforceloglike_a2(newtas,fT1,nMod);
        mRes = [mRes; nMod, pval1, pval2, cval1, cval2, kval1, kval2, fAIC, fL];
        % MOL with secondary aftershock (ppckk)
        nMod = 3; [pval1, pval2, cval1, cval2, kval1, kval2, fAIC, fL] = bruteforceloglike_a2(newtas,fT1,nMod);
        mRes = [mRes; nMod, pval1, pval2, cval1, cval2, kval1, kval2, fAIC, fL];
        % MOL with secondary aftershock (ppcckk)
        nMod = 4; [pval1, pval2, cval1, cval2, kval1, kval2, fAIC, fL] = bruteforceloglike_a2(newtas,fT1,nMod);
        mRes = [mRes; nMod, pval1, pval2, cval1, cval2, kval1, kval2, fAIC, fL];

        % Select best fitting model by AIC
        vSel = (mRes(:,8)==min(mRes(:,8)));
        mRes = mRes(vSel,:);
        if length(mRes(:,1)) > 1
            vSel = (mRes(:,1)==min(mRes(:,1)));
            mRes = mRes(vSel,:);
        end
        % Model to use for bootstrapping as of lowest AIC to observed data
        nMod = mRes(1,1);
        pval1= mRes(1,2); pval2= mRes(1,3);
        cval1= mRes(1,4); cval2= mRes(1,5);
        kval1= mRes(1,6); kval2= mRes(1,7);

        % Calculate goodness of fit with KS-Test and RMS
        [H,P,KSSTAT,fRMS] = calc_llkstest_a2(newtas,fT1,pval1, pval2, cval1, cval2, kval1, kval2, nMod);
        if H == 1
            nCnt = nCnt + 1;
            if nCnt/bootloops >= 0.5
                return
            end

        else
            cumnr_model = [];

            if nMod == 1
                for i=1:length(time_asf)
                    if pval1 ~= 1
                        nrmod = kval1/(pval1-1)*(cval1^(1-pval1)-(time_asf(i)+cval1)^(1-pval1));
                    else
                        nrmod = kval1*log(time_asf(i)/cval1+1);
                    end
                    cumnr_model = [cumnr_model; nrmod];
                end % END of FOR on length(time_asf)

                if pval1 ~= 1
                    cm = kval1/(pval1-1)*(cval1^(1-pval1)-(max(time_asf)+cval1)^(1-pval1));
                    cl = kval1/(pval1-1)*(cval1^(1-pval1)-(max(time_as)+cval1)^(1-pval1));
                else
                    cm = kval1*log(max(time_asf)/cval1+1);
                    cl = kval1*log(max(time_as)/cval1+1);
                end
            else
                for i=1:length(time_asf)
                    if time_asf(i) <= fT1
                        if pval1 ~= 1
                            nrmod = kval1/(pval1-1)*(cval1^(1-pval1)-(time_asf(i)+cval1)^(1-pval1));
                        else
                            nrmod = kval1*log(time_asf(i)/cval1+1);
                        end
                        cumnr_model = [cumnr_model; nrmod];
                    else
                        if (pval1 ~= 1 & pval2 ~= 1)
                            nrmod = kval1/(pval1-1)*(cval1^(1-pval1)-(time_asf(i)+cval1)^(1-pval1))+ kval2/(pval2-1)*(cval2^(1-pval2)-(time_asf(i)-fT1+cval2)^(1-pval2));
                        elseif (pval1 ~= 1  &&  pval2 == 1)
                            nrmod = kval1/(pval1-1)*(cval1^(1-pval1)-(time_asf(i)+cval1)^(1-pval1))+ kval2*log((time_asf(i)-fT1)/cval2+1);
                        elseif (pval1 == 1  &&  pval2 ~= 1)
                            nrmod = kval1*log(time_asf(i)/cval1+1)+ kval2/(pval2-1)*(cval2^(1-pval2)-(time_asf(i)-fT1+cval2)^(1-pval2));
                        else
                            nrmod = kval1*log(time_asf(i)/cval1+1) + kval2*log((time_asf(i)-fT1)/cval2+1);
                        end
                        cumnr_model = [cumnr_model; nrmod];
                    end %END of IF on fT1
                end % End of FOR length(time_asf)

                if (pval1 ~= 1 & pval2 ~= 1)
                    cm = kval1/(pval1-1)*(cval1^(1-pval1)-(max(time_asf)+cval1)^(1-pval1)) + kval2/(pval2-1)*(cval2^(1-pval2)-(max(time_asf)-fT1+cval2)^(1-pval2));
                    cl = kval1/(pval1-1)*(cval1^(1-pval1)-(max(time_as)+cval1)^(1-pval1)) + kval2/(pval2-1)*(cval2^(1-pval2)-(max(time_as)-fT1+cval2)^(1-pval2));
                elseif (pval1 ~= 1  &&  pval2 == 1)
                    cm = kval1/(pval1-1)*(cval1^(1-pval1)-(max(time_asf)+cval1)^(1-pval1)) + kval2*log((max(time_asf)-fT1)/cval2+1);
                    cl = kval1/(pval1-1)*(cval1^(1-pval1)-(max(time_as)+cval1)^(1-pval1)) + kval2*log((max(time_as)-fT1)/cval2+1);
                elseif (pval1 == 1  &&  pval2 ~= 1)
                    cm = kval1*log(max(time_asf)/cval1+1) + kval2/(pval2-1)*(cval2^(1-pval2)-(max(time_asf)-fT1+cval2)^(1-pval2));
                    cl = kval1*log(max(time_as)/cval1+1) + kval2/(pval2-1)*(cval2^(1-pval2)-(max(time_as)-fT1+cval2)^(1-pval2));
                else
                    cm = kval1*log(max(time_asf)/cval1+1) + kval2*log((max(time_asf)-fT1)/cval2+1);
                    cl = kval1*log(max(time_as)/cval1+1) + kval2*log((max(time_as)-fT1)/cval2+1);
                end
            end % End of if on nMod
            loopout = [loopout; pval1, pval2, cval1, cval2, kval1, kval2, fAIC, fL, cm, cl, nMod];
            pfloop = plot(time_asf,cumnr_model,'color',[0.8 0.8 0.8]);
        end % End of KS-Test
    end

    % 2nd moment of bootstrap number of forecasted number of events
    fStdBst = calc_StdDev(loopout(:,9));

    % plot observed number in learning period
    cumnrf = (1:length(time_asf))';
    cumnr = (1:length(time_as))';
    p2 = plot(time_as,cumnr,'b','Linewidth',2,'Linestyle','--');

    % Plot observed events in forecast period from endpoint of modeled events in learning period
    vSel = time_asf >= max(time_as);
    vCumnr_forecast = cumnrf(vSel,:);
    vTime_forecast = time_asf(vSel,:);
    % Difference of modelled and observed number of events at time_as
    fDiff_timeas = mean(loopout(:,10))-cumnrf(length(time_as));
    vCumnr_forecast = vCumnr_forecast+fDiff_timeas;
    pf3 = plot(vTime_forecast, vCumnr_forecast,'m-.','Linewidth',2);

    xlim([0 max(time_asf)]);
    xlabel('Delay time [days]')
    ylabel('Cumulative number of aftershocks')

    % Plot standard deviation from bootstrap
    ps2=errorbar(max(time_asf),mean(loopout(:,9)),fStdBst,fStdBst);
    set(ps2,'Linewidth',2,'Color',[1 0 0])

    % calculate rate change
    nr_forecast = mean(loopout(:,9));
    nr_learn = mean(loopout(:,10));
    nummod = nr_forecast-nr_learn;
    l = time_asf <= time+timef & time_asf > time;
    numreal = sum(l);
    fRc_Bst = (numreal-nummod)/fStdBst;

    % Set line for learning period
    yy = get(gca,'ylim');
    plot([max(time_as) max(time_as)],[0 yy(2)],'k-.')
    string=['\sigma = ' num2str(fStdBst) '; RC = ' num2str(fRc_Bst) ];
    text(max(time_asf)*0.15,yy(2)*0.15,string,'FontSize',10);

    legend([p2 pf3 min(ps2)],'data','observed','\sigma (Bst)','location', 'NorthWest');

author: CelsoReyes


text:
function [h,parvar,Q,v] = sosd(X,ops)
%SOSD Returns sum-of-squares decomposition
%
% h = sosd(p)
%
% Example
%
%  x = sdpvar(1,1);sdisplay(sosd(x^4 + 1))

if isfield(X.extra,'sosid')
    % This is used internally when calling sosd on SOS constraint
    [h,parvar,Q,v] = yalmip('getsos',X.extra.sosid);
else
    F = sos(X);
    if nargin == 1
        solvesos(F,[],sdpsettings('verbose',0));
    else
        solvesos(F,[],ops);
    end
    [h,parvar,Q,v] = sosd(F);
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [VMap,NMap,tMap,CMap] = raycastingMEX(camRtC2W, castingRange)

% VMap and NMap are in the world coordinate

% DDA based ray casting

global raycastingDirectionC;
global voxel;
global tsdf_value;
global tsdf_color;

camCenterW = getCameraCenter(camRtC2W);

raycastingDirectionW = transformRTdir(raycastingDirectionC,camRtC2W);

camCenterWgrid = (camCenterW - voxel.range(:,1)) / voxel.unit + 1;

tMap = NaN(1,640*480);
NMap = NaN(3,640*480);
CMap = NaN(3,640*480);

nearPlane = 0.3;
farPlane = 8.0;

step = voxel.unit;
largestep = 0.75 * voxel.mu;

tPrev = 0;

for i=1:(640*480)
    % Ray-Box Intersection to get the range

    origin = camCenterW;
    direction = raycastingDirectionW(:,i);

    % intersect ray with a box
    % http://www.siggraph.org/education/materials/HyperGraph/raytrace/rtinter3.htm
    % compute intersection of ray with all six bbox planes
    invR = 1 ./ direction;
    tbot = invR .* (voxel.range(:,1) + repmat(voxel.unit,3,1) - origin);
    ttop = invR .* (voxel.range(:,2) - repmat(voxel.unit,3,1) - origin);

    % re-order intersections to find smallest and largest on each axis
    tmin = min(ttop, tbot);
    tmax = max(ttop, tbot);

    % find the largest tmin and the smallest tmax
    largest_tmin = max(max(tmin(1), tmin(2)), max(tmin(1), tmin(3)));
    smallest_tmax = min(min(tmax(1), tmax(2)), min(tmax(1), tmax(3)));

    % check against near and far plane
    tnear = max(largest_tmin, nearPlane);
    tfar = min(smallest_tmax, farPlane);

    if tnear < tfar
        t = tnear;
        % first walk with largesteps until we found a hit
        stepsize = largestep;
        %f_t = interpolateTrilineary(origin + direction * t);
        f_t = interpolateTrilineary(camCenterWgrid + direction * t / voxel.unit);

        f_tt = 0;
        if f_t > 0     % ups, if we were already in it, then don't render anything here
            while t < tfar
                %f_tt = interpolateTrilineary(origin + direction * t);
                f_tt = interpolateTrilineary(camCenterWgrid + direction * t / voxel.unit);
                if f_tt < 0                               % got it, jump out of inner loop
                    break;
                end
                if f_tt < 0.8                            % coming closer, reduce stepsize
                    stepsize = step;
                end
                f_t = f_tt;
                t = t + stepsize;
            end
            if f_tt < 0                               % got it, calculate accurate intersection
                tMap(i) = t + stepsize * f_tt / (f_t - f_tt);
            end
        end
    end
end

% normalize normal map
NMap = NMap ./ repmat(sqrt(sum(NMap.^2,1)),3,1);

% computer vertex map
VMap = repmat(camCenterW,1,640*480) + raycastingDirectionW .* (repmat(tMap,3,1));

%imagesc(reshape(tMap,480,640)); axis equal; axis tight

return;

%camCenterWcopy = repmat( (camCenterW - voxel.range(:,1) ) / voxel.unit + 1,1,640*480);
%startPoints = camCenterWcopy + raycastingDirectionW * (castingRange(1)/ voxel.unit);
%endPoints = camCenterWcopy + raycastingDirectionW * (castingRange(2)/ voxel.unit);

% http://stellar.mit.edu/S/course/6/fa10/6.837/courseMaterial/topics/topic1/lectureNotes/13_RayTracing-Acceleration/13_RayTracing-Acceleration.pdf

% we assume the camera must be inside. otherwise, it is an error.
% so we don't test the camera

raycastingDirectionWinv = raycastingDirectionW.^-1;

maxTx = max((voxel.size_grid(1)-2-camCenterWgrid(1))*raycastingDirectionWinv(1,:),(2-camCenterWgrid(1))*raycastingDirectionWinv(1,:));
maxTy = max((voxel.size_grid(2)-2-camCenterWgrid(2))*raycastingDirectionWinv(2,:),(2-camCenterWgrid(2))*raycastingDirectionWinv(2,:));
maxTz = max((voxel.size_grid(3)-2-camCenterWgrid(3))*raycastingDirectionWinv(3,:),(2-camCenterWgrid(3))*raycastingDirectionWinv(3,:));
maxT = min(maxTx, min(maxTy, maxTz));

castingRangeGrid = castingRange/voxel.unit;
maxT = min(maxT, castingRangeGrid(2));

% parallel setting
parallelBlock = matlabpool('size');
blockDim = (640*480)/parallelBlock;
if blockDim ~= round(blockDim)
    error('thread does not align');
end

backMove = voxel.mu_grid;

tMap = NaN(1,640*480);
NMap = NaN(3,640*480);
CMap = NaN(3,640*480);

initDis = 3/ voxel.unit;

% http://www.cse.yorku.ca/~amana/research/grid.pdf
prevT = initDis;

for i=1:(640*480)
    % Ray-Box Intersection to get the range

    %{
    raycast( const Volume volume, const uint2 pos, const Matrix4 view, const float nearPlane, const float farPlane, const float step, const float largestep){
        const float3 origin = view.get_translation();
        const float3 direction = rotate(view, make_float3(pos.x, pos.y, 1.f));

        // intersect ray with a box
        // http://www.siggraph.org/education/materials/HyperGraph/raytrace/rtinter3.htm
        // compute intersection of ray with all six bbox planes
        const float3 invR = make_float3(1.0f) / direction;
        const float3 tbot = -1 * invR * origin;
        const float3 ttop = invR * (volume.dim - origin);

        // re-order intersections to find smallest and largest on each axis
        const float3 tmin = fminf(ttop, tbot);
        const float3 tmax = fmaxf(ttop, tbot);

        // find the largest tmin and the smallest tmax
        const float largest_tmin = fmaxf(fmaxf(tmin.x, tmin.y), fmaxf(tmin.x, tmin.z));
        const float smallest_tmax = fminf(fminf(tmax.x, tmax.y), fminf(tmax.x, tmax.z));

        // check against near and far plane
        const float tnear = fmaxf(largest_tmin, nearPlane);
        const float tfar = fminf(smallest_tmax, farPlane);

        if(tnear < tfar) {
            // first walk with largesteps until we found a hit
            float t = tnear;
            float stepsize = largestep;
            float f_t = volume.interp(origin + direction * t);
            float f_tt = 0;
            if( f_t > 0){     // ups, if we were already in it, then don't render anything here
                for(; t < tfar; t += stepsize){
                    f_tt = volume.interp(origin + direction * t);
                    if(f_tt < 0)                               // got it, jump out of inner loop
                        break;
                    if(f_tt < 0.8f)                            // coming closer, reduce stepsize
                        stepsize = step;
                    f_t = f_tt;
                }
                if(f_tt < 0){                               // got it, calculate accurate intersection
                    t = t + stepsize * f_tt / (f_t - f_tt);
                    return make_float4(origin + direction * t, t);
                }
            }
        }
        return make_float4(0);
    }
    %}

    tMin = castingRangeGrid(1);
    tMax = maxT(i);

    if tMin<tMax % max range > min range

        tStart = max(castingRangeGrid(1),min(tMax,prevT-backMove));

        rayDir = raycastingDirectionW(:,i);

        startPoint = camCenterWgrid + rayDir*tStart;

        % The initialization phase begins by identifying the voxel in which the ray origin, ?u, is found.
        % If the ray origin is outside the grid, we find the point in which the ray enters the grid and take the adjacent voxel.
        % The integer variables X and Y are initialized to the starting voxel coordinates.
        X = round(startPoint(1));
        Y = round(startPoint(2));
        Z = round(startPoint(3));
        valCurrentVoxel = tsdf_value(X,Y,Z);
        %valCurrentVoxel = interpValue(startPoint);

        tOptimal = NaN;

        if valCurrentVoxel == 0
            % lucky, you got the surface immediately!
            %tMap(i) = tStart;
            prevT = tStart;

            tOptimal = tStart;
        else
            if valCurrentVoxel>0
                localDir = rayDir;
                tMax = tMax-tStart;
                lookforSign = -1;
            else
                localDir = -rayDir;
                tMax = tStart-tMin;
                lookforSign = +1;
            end

            % In addition, the variables stepX and stepY are initialized to either 1 or -1 indicating
            % whether X and Y are incremented or decremented as the ray crosses voxel boundaries
            % (this is determined by the sign of the x and y components of ?v).
            stepX = sign(localDir(1));
            stepY = sign(localDir(2));
            stepZ = sign(localDir(3));

            % Next, we determine the value of t at which the ray crosses the ?rst vertical voxel boundary
            % and store it in variable tMaxX. We perform a similar computation in y and store the result in tMaxY.
            % The minimum of these two values will indicate how much we can travel along the ray and still remain in the current voxel.

            if localDir(1)>0
                tMaxX = (X+0.5 - startPoint(1))/localDir(1);
            elseif localDir(1)<0
                tMaxX = (X-0.5 - startPoint(1))/localDir(1);
            else
                tMaxX = Inf;
            end
            if localDir(2)>0
                tMaxY = (Y+0.5 - startPoint(2))/localDir(2);
            elseif localDir(2)<0
                tMaxY = (Y-0.5 - startPoint(2))/localDir(2);
            else
                tMaxY = Inf;
            end
            if localDir(3)>0
                tMaxZ = (Z+0.5 - startPoint(3))/localDir(3);
            elseif localDir(3)<0
                tMaxZ = (Z-0.5 - startPoint(3))/localDir(3);
            else
                tMaxZ = Inf;
            end

            if min(min(tMaxX,tMaxY),tMaxZ)<0
                error('tStart<0');
            end

            % Finally, we compute tDeltaX and tDeltaY.
            % tDeltaX indicates how far along the ray we must move (in units of t) for the horizontal component of such a movement to equal the width of a voxel.
            % Similarly,we store in tDeltaY the amount of movement along the ray which has a vertical component equal to the height of a voxel.
            tDeltaX = 1/abs(localDir(1));
            tDeltaY = 1/abs(localDir(2));
            tDeltaZ = 1/abs(localDir(3));

            t = 0;
            while t<tMax
                if tMaxX < tMaxY
                    if tMaxX < tMaxZ
                        X= X + stepX;
                        tMaxX= tMaxX + tDeltaX;
                    else
                        Z= Z + stepZ;
                        tMaxZ= tMaxZ + tDeltaZ;
                    end
                else
                    if tMaxY < tMaxZ
                        Y= Y + stepY;
                        tMaxY= tMaxY + tDeltaY;
                    else
                        Z= Z + stepZ;
                        tMaxZ= tMaxZ + tDeltaZ;
                    end
                end
                prevT = t;
                t = min(min(tMaxX,tMaxY),tMaxZ);

                if tsdf_value(X,Y,Z) * lookforSign > 0

                    if (tsdf_value(X,Y,Z)==-1 && valCurrentVoxel==+1) || (tsdf_value(X,Y,Z)==+1 && valCurrentVoxel==-1)
                        tOptimal = NaN;
                    else

                        % found the zero crossing points!

                        % simplest
                        % tOptimal = t;

                        % simple average
                        tOptimal = (prevT + t)/2;

                        %{
                        % buggy: Ftdt == Ft

                        if lookforSign>0
                            prevTswap = prevT;
                            prevT = t;
                            t= prevTswap;
                        end
                        XYZt   = startPoint + localDir*t;
                        XYZtdt = startPoint + localDir*prevT;

                        Ft = interpolateTrilineary(XYZt(1),XYZt(2),XYZt(3));
                        Ftdt = interpolateTrilineary(XYZtdt(1),XYZtdt(2),XYZtdt(3));

                        tOptimal = t - abs(t-prevT) * Ft / (Ftdt - Ft);
                        if tOptimal>1000
                            error('tOptimal>1000');
                        end
                        %}

                        %fprintf('i=%d: t=%f tOptimal=%f Ft=%f Ftdt=%f\n',i,t,t - abs(t-prevT) * Ft / (Ftdt - Ft), Ft, Ftdt);
                        tOptimal = tStart - tOptimal * lookforSign;

                        prevT = tOptimal;

                    end
                    break;
                end

                valCurrentVoxel = tsdf_value(X,Y,Z);
            end

        end
        if ~isnan(tOptimal)
            tMap(i) = tOptimal;

            % compute normal map

            XYZgrid = camCenterWgrid + rayDir*tOptimal;

            NMap(1,i) = interpolateTrilineary(XYZgrid(1)+1,XYZgrid(2),XYZgrid(3))-interpolateTrilineary(XYZgrid(1)-1,XYZgrid(2),XYZgrid(3));
            NMap(2,i) = interpolateTrilineary(XYZgrid(1),XYZgrid(2)+1,XYZgrid(3))-interpolateTrilineary(XYZgrid(1),XYZgrid(2)-1,XYZgrid(3));
            NMap(3,i) = interpolateTrilineary(XYZgrid(1),XYZgrid(2),XYZgrid(3)+1)-interpolateTrilineary(XYZgrid(1),XYZgrid(2),XYZgrid(3)-1);

            if ~isempty(tsdf_color)
                CMap(:,i) = interpolateTrilinearyColor(XYZgrid(1),XYZgrid(2),XYZgrid(3));
            end
        end
    end
end

% normalize normal map
NMap = NMap ./ repmat(sqrt(sum(NMap.^2,1)),3,1);

% computer vertex map
VMap = repmat(camCenterW,1,640*480) + raycastingDirectionW .* (repmat(tMap*voxel.unit,3,1));

%imagesc(reshape(tMap,480,640)); axis equal; axis tight

author: jianxiongxiao

repo: ProfXkit


text:
function PlotCSD(csd,varargin)

%PlotCSD - Plot current source density.
%
%  USAGE
%
%    PlotCSD(csd,<options>)
%
%    csd            current source density (see CSD)
%    <options>      optional list of property-value pairs (see table below)
%
%    =========================================================================
%     Properties    Values
%    -------------------------------------------------------------------------
%     'lfp'         local field potential
%     'scale'       scale (arbitrary units) for LFP traces (default = 1)
%     'cutoffs'     cutoff values (default = [-M M] where M is the maximum
%                   amplitude of the CSD)
%    =========================================================================
%
%  SEE
%
%    See CSD.

% Copyright (C) 2008-2011 by Michaël Zugaro
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.

% Default values
lfp = [];
scale = 1;
cutoffs = [];

% Check number of parameters
if nargin < 1 | mod(length(varargin),2) ~= 0,
  error('Incorrect number of parameters (type ''help PlotCSD'' for details).');
end

% Parse parameter list
for i = 1:2:length(varargin),
	if ~ischar(varargin{i}),
		error(['Parameter ' num2str(i+1) ' is not a property (type ''help PlotCSD'' for details).']);
	end
	switch(lower(varargin{i})),
		case 'lfp',
			lfp = varargin{i+1};
			if ~isdmatrix(lfp) | size(lfp,1) ~= size(csd,1) | size(lfp,2) ~= size(csd,2)+2,
				error('Incorrect size for parameter ''lfp'' (type ''help PlotCSD'' for details).');
			end
		case 'scale',
			scale = varargin{i+1};
			if ~isdscalar(scale,'>0'),
				error('Incorrect value for property ''scale'' (type ''help PlotCSD'' for details).');
			end
		case 'cutoffs',
			cutoffs = varargin{i+1};
			if ~isdvector(cutoffs,'#2','<'),
				error('Incorrect value for property ''cutoffs'' (type ''help PlotCSD'' for details).');
			end
		otherwise,
			error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help PlotCSD'' for details).']);
	end
end

d = csd(:,2:end);
d = interp2(d);
d = d(1:2:size(d,1),:);

pcolor(csd(:,1),1:size(d,2),flipud(transpose(d)));
shading interp;
if ~isempty(cutoffs),
	set(gca,'clim',cutoffs);
end

if ~isempty(lfp),
	hold on;
	y = lfp(:,2:end);
	y = y - repmat(mean(y),size(y,1),1);
	y = y / max(max(abs(y)))*scale;
	n = size(y,2);
	for i = 1:n,
		plot(lfp(:,1),y(:,i)+(n-i)*2-1,'k');
    end
    ylim([-2 (n-1)*2]);
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function pass = test_diag(pref)
% Test for CHEBFUN/DIAG()

if ( nargin < 1 )
    pref = chebfunpref();
end

%% Test a few things
% No breakpoints
x = chebfun(@(x) x, [-1 6]);
tol = 1e-15;
f = sin(x);
g = cos(4*x);

D = diag(f);
err = norm(D*g -f.*g);
pass(1) = ( err < tol );

%% Breakpoints in g
x = chebfun(@(x) x, [-1 4 6]);
g = cos(4*x);
err = norm(D*g -f.*g);
pass(2) = ( err < tol );

%% Breakpoints in f and g
x = chebfun(@(x) x, [-1 4 6]);
f = sin(x);
g = cos(4*x);
err = norm(D*g -f.*g);
pass(3) = ( err < 10*tol );

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function I = eye(disc)
%EYE   Identity operator for VALSDISCRETIZATION discretization.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

n = disc.dimension;
I = eye(sum(n));

end


text:
function geometry_test056 ( )

%*****************************************************************************80
%
%% TEST056 tests PLANE_IMP_LINE_PAR_INT_3D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 April 2009
%
%  Author:
%
%    John Burkardt
%
  dim_num = 3;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST056\n' );
  fprintf ( 1, '  PLANE_IMP_LINE_PAR_INT_3D finds the \n' );
  fprintf ( 1, '    intersection of an implicit plane and\n' );
  fprintf ( 1, '    a parametric line, in 3D.\n' );

  a = 1.0;
  b = -2.0;
  c = -3.0;
  d = 6.0;

  f = 2.0;
  g = 1.0;
  h = 5.0;
  x0 = 3.0;
  y0 = 0.0;
  z0 = -7.0;

  [ intersect, p ] = plane_imp_line_par_int_3d ( a, b, c, d, x0, y0, z0, f, g, h );

  if ( intersect )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  The plane and line intersect at \n' );
    fprintf ( 1, '  %f  %f  %f\n', p(1:dim_num) );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, '  The plane and the line do not intersect.\n' );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Expected answer:\n' );
  fprintf ( 1, '    The plane and line intersect at \n' );
  fprintf ( 1, '    7, 2, 3.\n' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Test file for @chebtech/any.m.

function pass = test_any(pref)

if ( nargin < 1 )
    pref = chebtech.techPref();
end

for n = 1:2
    if ( n == 1 )
        testclass = chebtech1();
    else
        testclass = chebtech2();
    end

    % Check behavior for any() down columns.
    pass(n, 1) = ~any(testclass);

    f = testclass.make(@(x) [sin(x) 0*x cos(x)]);
    pass(n, 2) = isequal(any(f), [1 0 1]);

    % Check behavior for any() across rows.
    g = any(f, 2);
    pass(n, 3) = isequal(g.coeffs, 1);

    f = testclass.make(@(x) [0*x 0*x]);
    g = any(f, 2);
    pass(n, 4) = isequal(g.coeffs, 0);
end

end


text:
filename = 'Cantilever';
ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance'};
weights = [1];
constraint = {'volumeConstraint'};
constraint_case = 'EQUALITY';
optimizerUnconstrained = 'SLERP';
optimizer = 'AlternatingPrimalDual';
incrementFactor = 1.05;
designVariable = 'LevelSet';
filterType = 'PDE';
line_search_initiator = 'INCREASING LAST STEP';

nsteps = 1;
Vfrac_final = 0.3;
optimality_final =1e-3;
constr_final = 1e-5;

Vfrac_initial = 1;
optimality_initial = 1e-3;
constr_initial = 1e-3;
Perimeter_target = 5;

TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

plotting = false;
printing = false;
printing_physics = false;
monitoring = false;
monitoring_interval = 10;
maxiter = 15;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [img, scale, im_scaled_size] = prepare_img_blob(img, scale, mean_pix, max_size)
%
% This file is part of the code that implements the following paper:
% Title      : "LocNet: Improving Localization Accuracy for Object Detection"
% Authors    : Spyros Gidaris, Nikos Komodakis
% Institution: Universite Paris Est, Ecole des Ponts ParisTech
% ArXiv link : http://arxiv.org/abs/1511.07763
% code       : https://github.com/gidariss/LocNet
%
% Part of the code in this file comes from the Faster R-CNN code and the matlab
% re-implementation of Fast-RCNN (https://github.com/ShaoqingRen/faster_rcnn)
%
% AUTORIGHTS
% --------------------------------------------------------
% Copyright (c) 2016 Spyros Gidaris
%
% Title     : "LocNet: Improving Localization Accuracy for Object Detection"
% ArXiv link: http://arxiv.org/abs/1511.07763
% Licensed under The MIT License [see LICENSE for details]
% ---------------------------------------------------------

if numel(mean_pix) == 1, mean_pix = [mean_pix, mean_pix, mean_pix]; end
assert(numel(mean_pix)==3);

im_height = size(img, 1);
im_width = size(img, 2);
% get the image size of the scaled image
im_scaled_size = prepare_img_blob_size([im_height, im_width], scale, max_size);
scale = min(im_scaled_size);

% scale image
if (scale <= 224)
    img = imresize(img, [im_scaled_size(1), im_scaled_size(2)], 'bilinear');
else
    img = imresize(img, [im_scaled_size(1), im_scaled_size(2)], 'bilinear', 'antialiasing', false);
end

img = single(img); % transform to single precision
img = img(:,:,[3 2 1]); % transform from RGB -> BGR (necessary for Caffe)
% subtruct the mean pixel from the image
img = bsxfun(@minus, img, reshape(mean_pix, [1, 1, 3]));
img = permute(img, [2 1 3]); % permute from [Height x Width x 3] to [Width x Height x 3] (necessary for Caffe)
end

author: gidariss

repo: LocNet

save_path: github-repos/MATLAB/gidariss-LocNet


text:
%% AKAZE and ORB planar tracking
%
% In this demo, we will compare AKAZE and ORB local features by using them to
% find matches between video frames and track object movements.
%
% The algorithm is as follows: First, we detect and describe keypoints on the
% first frame inside a manually set object boundaries. For every next frame:
%
% * Detect and describe keypoints
% * Match them using bruteforce matcher
% * Estimate homography transformation using RANSAC
% * Filter inliers from all the matches
% * Apply homography transformation to the bounding box to find the object
% * Draw bounding box and inliers, compute inlier ratio as evaluation metric
%
% To do the tracking we need a video and object position on the first frame.
% You can download an example video and data from
% <https://docs.google.com/file/d/0B72G7D4snftJandBb0taLVJHMFk here>.
%
% To run the sample, you specify an input (video file or camera id), then
% select a bounding box with the mouse, to start tracking.
%
% Example video: <http://www.youtube.com/watch?v=pzVbhxx6aog>.
%
% Sources:
%
% * <https://docs.opencv.org/3.3.1/dc/d16/tutorial_akaze_tracking.html>
% * <https://github.com/opencv/opencv/blob/3.3.1/samples/cpp/tutorial_code/features2D/AKAZE_tracking/planar_tracking.cpp>
% * <https://github.com/opencv/opencv/blob/3.3.1/samples/python/plane_tracker.py>
%

function varargout = planar_tracker_demo(fname)
    % initialize app state and options
    app = defaultOptions();

    % load video
    if nargin < 1
        fname = fullfile(mexopencv.root(), 'test', 'blais.mp4');
        win = [136 0 366 433];     % book
        %win = [135 165 285 175];  % face
        if exist(fname, 'file') ~= 2
            error('Missing "blais.mp4" video, download it from here: %s', ...
                'https://docs.google.com/file/d/0B72G7D4snftJandBb0taLVJHMFk')
        end
    else
        win = [];
    end
    vid = cv.VideoCapture(fname);
    assert(vid.isOpened(), 'Failed to open video');

    % get first frame
    frame = vid.read();
    assert(~isempty(frame), 'Failed to read frames');

    % prompt user for an object ROI
    if isempty(win)
        win = selectROI(frame)
        assert(~isempty(win), 'No object specified');
    end

    % store first frame and bounding box points
    app.sz = size(frame);
    app.frame0 = frame;
    app.bb0 = bsxfun(@plus, win(1:2), [0 0; win(3) 0; win(3:4); 0 win(4)]);

    % create feature detectors and descriptor matchers objects
    %names = {'AKAZE', 'KAZE', 'ORB', 'BRISK', 'SIFT', 'SURF'};
    names = {'AKAZE', 'ORB'};
    f = createTrackers(app, names);

    % detect keypoints and compute descriptors in first frame (ROI region)
    for i=1:numel(f)
        if true && i>1 && strcmp(f(i).name, 'ORB')
            % to ensure detectors locate roughly same number of keypoints
            f(i).detector.MaxFeatures = f(1).stats(1);
        end
        [f(i).kp0, f(i).desc0, f(i).stats] = processFirstFrame(f(i), app);
    end
    stats = cat(1, f.stats);

    % prepare plot
    h = createUI(f, app);
    if nargout > 0, varargout{1} = h; end

    % main loop
    counter = 0;
    while ishghandle(h(1).img)
        % get next frame
        counter = counter + 1;
        frame = vid.read();
        if isempty(frame), break; end
        if rem(counter, app.show_every_frame) > 0, continue; end

        for i=1:numel(f)
            % feature matching + homography to track object
            [out, bb, f(i).stats] = processFrame(f(i), app, frame);

            % check quality of match
            if all(f(i).stats(3:4) >= [app.min_inliers app.min_inliers_ratio])
                lin_opts = {'LineWidth',2, 'LineStyle','-'};
            else
                % uncertain: not enough matches and/or too many outliers
                lin_opts = {'LineWidth',0.5, 'LineStyle',':'};
            end

            % show results
            set(h(i).txt, 'String',printStats(f(i).stats));
            set(h(i).lin, 'XData',bb([1:end 1],1) + app.sz(2), ...
                'YData',bb([1:end 1],2), lin_opts{:});
            set(h(i).img, 'CData',out);
        end
        stats = stats + cat(1, f.stats);
        drawnow;
    end
    vid.release();

    % show average stats
    stats = stats ./ (counter / app.show_every_frame);
    stats(:,1:3) = round(stats(:,1:3));
    for i=1:numel(f)
        disp(['-- ' f(i).name ' --']);
        cellfun(@disp, printStats(stats(i,:)));
    end
end

%% Helper functions

function app = defaultOptions()
    %DEFAULTOPTIONS  Create default options structure

    app = struct();
    app.akaze_thresh = 3e-4;     % AKAZE threshold set to locate ~1000 keypoints
    app.ransac_thresh = 2.5;     % RANSAC inlier threshold
    app.match_ratio = 0.8;       % nearest-neighbour matching ratio
    app.min_inliers = 100;       % minimal number of inliers to draw bounding box
    app.min_inliers_ratio = 0.2; % minimal inliers ratio to draw bounding box
    app.draw_top_matches = 100;  % draw the top 50 matches only
    app.show_every_frame = 3;    % display every 3 frames to maintain high fps
end

function win = selectROI(img)
    %SELECTROI  Interactively select region with mouse

    hImg = imshow(img);
    axis on
    hRect = imrect();
    win = wait(hRect);
    close(ancestor(hImg, 'Figure'));
    if ~isempty(win)
        win(1:2) = win(1:2) - 1;
    end
end

function f = createTrackers(app, names)
    %CREATETRACKERS  Create detectors and matchers objects

    f = struct();
    for i=1:numel(names)
        % create detector object
        f(i).name = names{i};
        switch names{i}
            case 'AKAZE'
                f(i).detector = cv.AKAZE('Threshold',app.akaze_thresh);
            case 'KAZE'
                f(i).detector = cv.KAZE('Threshold',app.akaze_thresh);
            case 'ORB'
                f(i).detector = cv.ORB('MaxFeatures',1000);
            case 'BRISK'
                f(i).detector = cv.BRISK();
            case 'SIFT'
                f(i).detector = cv.SIFT();
            case 'SURF'
                f(i).detector = cv.SURF();
            otherwise
                error('Unexpected detector type')
        end

        % create matcher object
        if true
            if false
                opts = {'CrossCheck',true};
            else
                opts = {};
            end
            f(i).matcher = cv.DescriptorMatcher('BFMatcher', ...
                'NormType',f(i).detector.defaultNorm(), opts{:});
        else
            f(i).matcher = cv.DescriptorMatcher('FlannBasedMatcher', ...
                'Index',...
                {'LSH', 'TableNumber',6, 'KeySize',12, 'MultiProbeLevel',1});
        end
    end
end

function h = createUI(f, app)
    %CREATEUI  Creates the UI

    x = app.bb0([1:end 1], 1);
    y = app.bb0([1:end 1], 2);
    lin_opts = {'Color','b', 'LineWidth',2, 'LineStyle','-'};

    % plot images, objects, and stats
    h = struct();
    figure('Position',get(0, 'DefaultFigurePosition') .* [0.4 0.2 1.5 1.7])
    for i=1:numel(f)
        subplot(numel(f),1,i)
        h(i).img = imshow([app.frame0, app.frame0]);
        line(x, y, lin_opts{:})
        h(i).lin = line(x + app.sz(2), y, lin_opts{:});
        h(i).txt = text(2*app.sz(2), app.sz(1), printStats(f(i).stats), ...
            'Color','y', 'HorizontalAlignment','right', 'VerticalAlignment','bottom');
        title(sprintf('%s, %d keypoints', f(i).name, f(i).stats(1)))
    end
end

function str = printStats(stats)
    %PRINTSTATS  Pretty print matching statistics

    str = {
        sprintf('Keypoints: %d', stats(1))
        sprintf('Matches: %d', stats(2))
        sprintf('Inliers: %d', stats(3))
        sprintf('Inliers Ratio: %.0f%%', stats(4)*100)
        sprintf('FPS: %.2f', stats(5))
    };
end

function [kp, desc, stats] = processFirstFrame(f, app)
    %PROCESSFIRSTFRAME Process first frame

    % init stats
    stats = [0, 0, 0, 0, 0];

    % create ROI mask
    mask = zeros(app.sz(1:2), 'uint8');
    mask = cv.fillPoly(mask, app.bb0, 'Color',255);

    tm = cv.TickMeter();
    tm.start();

    % detect and compute features of first frame in ROI region
    [kp, desc] = f.detector.detectAndCompute(app.frame0, 'Mask',mask);
    assert(~isempty(kp), 'No keypoints detected in first frame');
    stats(1) = numel(kp);

    tm.stop();
    stats(5) = 1 / tm.TimeSec;
end

function [out, bb1, stats] = processFrame(f, app, frame1)
    %PROCESSFRAME  Process subsequent frames

    % initialize output
    out = [app.frame0, frame1];
    bb1 = nan(size(app.bb0));
    stats = [0 0 0 0 0];

    tm = cv.TickMeter();
    tm.start();

    % detect and compute features in current frame
    [kp1, desc1] = f.detector.detectAndCompute(frame1);
    stats(1) = numel(kp1);
    if numel(kp1) < 2, return; end

    % match against first frame features
    if true
        % 2-NN matching with ratio test (if closest match is MATCH_RATIO
        % closer than the second closest one, then its a good match)
        matches = f.matcher.knnMatch(f.desc0, desc1, 2);
        idx = cellfun(@(m) numel(m) == 2 && ...
            (m(1).distance < app.match_ratio * m(2).distance), matches);
        matches = cellfun(@(m) m(1), matches(idx));
    elseif false
        matches = f.matcher.knnMatch(f.desc0, desc1, 1);
        matches = [matches{:}];
    else
        matches = f.matcher.match(f.desc0, desc1);
    end
    stats(2) = numel(matches);
    if numel(matches) < 4, return; end

    % estimate homography using RANSAC
    pt0 = cat(1, f.kp0([matches.queryIdx] + 1).pt);
    pt1 = cat(1, kp1([matches.trainIdx] + 1).pt);
    [H, inliers] = cv.findHomography(pt0, pt1, ...
        'Method','Ransac', 'RansacReprojThreshold',app.ransac_thresh);
    inliers = logical(inliers);
    stats(3) = nnz(inliers);
    stats(4) = stats(3) / stats(2);
    if isempty(H), return; end

    tm.stop();
    stats(5) = 1 / tm.TimeSec;

    % project object bounding box using homography to locate it in new frame
    bb1 = cv.perspectiveTransform(app.bb0, H);

    % draw matches
    if app.draw_top_matches > 0
        % only draw top inlier matches (sorted by distance) to reduce clutter
        [~,idx] = sort([matches.distance]);
        idx(~inliers(idx)) = [];
        inliers(idx(app.draw_top_matches+1:end)) = false;
    end
    out = cv.drawMatches(app.frame0, f.kp0, frame1, kp1, matches, ...
        'MatchesMask',inliers, ...
        'MatchColor',[255 0 0], 'SinglePointColor',[255 0 255]);
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [deformation_field, deformed_image] = PTKSolveForFluidRegistration(image_to_transform, reference_image, reporting, affine_initial_matrix)
    % PTKSolveForFluidRegistration. Computes the deformation field to register one
    %     image segmentation to another based on solving a fluid registration.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    if nargin < 4
        affine_initial_matrix = [];
    end

    if isempty(affine_initial_matrix)
        affine_initial_matrix = GetInitialTransformation(image_to_transform, reference_image, reporting);
    end

    % For the fluid registration, the images must be the same size and have the
    % same voxel size
    isotropic_reference_image = MimImageUtilities.MakeImageApproximatelyIsotropic(reference_image, 'PTK smoothed binary');

    % The fluid registration will be performed after the rigid registration has
    % been performed
    transformed_image = PTKRegisterImageAffine(image_to_transform, isotropic_reference_image, affine_initial_matrix, '*nearest', reporting);

    deformation_field = PTKSolveMatchedImagesForFluidRegistration(transformed_image, isotropic_reference_image, reporting);

    deformation_field = MimImageCoordinateUtilities.AdjustDeformationFieldForInitialAffineTransformation(deformation_field, affine_initial_matrix);

    deformed_image = PTKRegisterImageFluid(image_to_transform, deformation_field, '*nearest', reporting);
end

function affine_initial_matrix = GetInitialTransformation(image_to_transform, reference_image, reporting)
    % Initially perform a centroid registration
    [affine_initial_matrix, ~] = PTKRegisterCentroid(image_to_transform, reference_image, reporting);

%    % Alternative: perform an initial rigid registration
%     [affine_initial_matrix, transformed_matrix] = PTKSolveForRigidRegistration(image_to_transform, reference_image, reporting);
end

author: tomdoel

repo: pulmonarytoolkit


text:
clear all; close all; clc;
addpath('./toolbox');

% load feature trajectories
load('./data/running.mat');

% define some variables
dtrng			= 20:0.5:40;
arng			= 0.8:0.02:1.2;
projmodel = 'affine';
a_con			= [];						% set to a value to constrain a

% define noise parameters
nnoise		= 21;
nits			= 20;
noisevec	= linspace(0,5,nnoise);
dts				= zeros(nits,nnoise);

% offset W1 by 30 frames
W1		= W1(31:end-30,:);

for n = 1:nnoise
	for it = 1:nits
		W1n = W1 + noisevec(n)*randn(size(W1));
		W2n	= W2 + noisevec(n)*randn(size(W2));

		% run synchronization script
		[a,dt] = sync(W1n,W2n,a_con,projmodel);

		dts(it,n)	= dt;
	end
end

% this took a while so save to a temporary file just in case
save('temp.mat','dts');

mn	= mean(dts,1);
sd	= 3*std(dts,[],1);

figure;
	plot([0,5],[30,30],'r--',...
				noisevec,dts,'b.');
	axis('equal',[minmax(noisevec),28,32]);
	xlabel('\sigma_n (pixels)'); ylabel('Recovered offset (frames)');

author: Sable

repo: mcbench-benchmarks


text:
function [steeringModel, ok] = promptForSteeringModelType(selectedEnum)
    [listboxStr, enums] = SteerModelTypeEnum.getListBoxStr();

    if(not(isempty(selectedEnum)))
        initialValue = find(selectedEnum == enums, 1, 'first');

        if(isempty(initialValue))
            initialValue = 1;
        end
    else
        initialValue = 1;
    end

    out = AppDesignerGUIOutput();
    listdlgARH_App('ListString',listboxStr, ...
                    'Name','Select Steering Model Type', ...
                    'PromptString',{'Select steering model type:'}, ...
                    'SelectionMode','single', ...
                    'ListSize',[300 300], ...
                    'InitialValue', initialValue, ...
                    'out',out);
    Selection = out.output{1};
    ok = out.output{2};

    if(ok)
        enum = SteerModelTypeEnum.getEnumForListboxStr(listboxStr{Selection});

        switch enum
            case SteerModelTypeEnum.PolyAngles
                steeringModel = RollPitchYawPolySteeringModel.getDefaultSteeringModel();

            case SteerModelTypeEnum.QuaterionInterp
                steeringModel = GenericQuatInterpSteeringModel.getDefaultSteeringModel();

            case SteerModelTypeEnum.LinearTangentAngles
                steeringModel = GenericLinearTangentSteeringModel.getDefaultSteeringModel();

            case SteerModelTypeEnum.SumOfSinesAngles
                steeringModel = GenericSumOfSinesSteeringModel.getDefaultSteeringModel();

            case SteerModelTypeEnum.SelectableModelAngles
                steeringModel = GenericSelectableSteeringModel.getDefaultSteeringModel();

            otherwise
                error('Unknown steering model type: %s', enum.name);
        end
    else
        steeringModel = RollPitchYawPolySteeringModel.getDefaultSteeringModel();
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
% PQ_TOP queries for the topmost element of the priority queue (not removing it)
%
% SYNTAX
% [idx, cost] = pq_top(pq)
%
% INPUT PARAMETERS
%   pq: a pointer to the priority queue
%
% OUTPUT PARAMETERS
%   idx: the index of the topmost element
%   cost: the cost of the topmost element
%
% DESCRIPTION
% Queries the topmost element from a priority queue returning its
% index and associated cost.
%
% See also:
% PQ_DEMO, PQ_CREATE, PQ_PUSH, PQ_POP, PQ_SIZE, PQ_TOP, PQ_DELETE
%
% References:
% Gormen, T.H. and Leiserson, C.E. and Rivest, R.L., "introduction to
% algorithms", 1990, MIT Press/McGraw-Hill, Chapter 6.

% Copyright (c) 2008 Andrea Tagliasacchi
% All Rights Reserved
% email: andrea.tagliasacchi@gmail.com
% $Revision: 1.0$  Created on: May 22, 2009

author: Cloud-CV

repo: object-proposals


text:
function y = power(x,d)
%POWER (overloaded)

% Vectorize x if d is vector
if numel(x)==1 & (numel(d)>1)
    x = x.*ones(size(d));
end
% Vectorize if x is a vector
if numel(d)==1 & (numel(x)>1)
    d = d.*ones(size(x));
end
if ~isequal(size(d),size(x))
    error('Dimension mismatch in power');
end

% Reuse code
if numel(x)==1 && numel(d)==1
    y = mpower(x,d);
    return
end

if isa(d,'sdpvar')
    % Call helper which vectorizes the elements
    y = powerinternalhelper(d,x);
    if isa(y,'sdpvar')
        y.extra.createTime = definecreationtime;
    end
    return
end

% Sanity Check
if prod(size(d))>1
    if any(size(d)~=size(x))
        error('Matrix dimensions must agree.');
    end
else
    d = ones(x.dim(1),x.dim(2))*d;
end

% Trivial cases
if isnumeric(d)
    if all(all(d==0))
        y = ones(x.dim(1),x.dim(2));
        return
    end
    if all(all(d==1))
        y = x;
        return
    end
    if isnan(d)
        disp('You have NaNs in model (learn to debug)')
        error('NaN power makes no sense.');
    end
end

% Fractional, negative or different powers are
% treated less efficiently using simple code.
fractional = any(any((ceil(d)-d>0)));
negative = any(any(d<0));
different = ~all(all(d==d(1)));
if fractional | negative | different
    if x.dim(1)>1 | x.dim(2)>1
        if isequal(x.basis,[spalloc(prod(x.dim),1,0) speye(prod(x.dim))]) & all(d==d(1))
            % Simple case x.^d
            y = vectorizedUnitPower(x,d);
            y.extra.createTime = definecreationtime;
            return
        end
        [n,m] = size(x);
        y = [];
        for i = 1:n % FIX : Vectorize!
            if m == 1
                temp = extsubsref(x,i,1).^d(i,1);
            else
                temp = [];
                for j = 1:m
                    temp = [temp extsubsref(x,i,j).^d(i,j)];
                end
            end
            y = [y;temp];
        end
        y.extra.createTime = definecreationtime;
        return
    else
        base = getbase(x);
        if isequal(base,[0 1])
            mt = yalmip('monomtable');
            var = getvariables(x);
            previous_var = find((mt(:,var)==d)  & (sum(mt~=0,2)==1));
            if isempty(previous_var)
                mt(end+1,:) = mt(getvariables(x),:)*d;
                yalmip('setmonomtable',mt);
                y = recover(size(mt,1));
            else
                y = recover(previous_var);
            end
        elseif (size(base,2) == 2) & base(1)==0
            % Something like a*t^-d
            y = base(2)^d*recover(getvariables(x))^d;
        else
            error('Only unit scalars can have negative or non-integer powers.');
        end
    end
    y.extra.createTime = definecreationtime;
    return
end

if isequal(x.basis,[spalloc(prod(x.dim),1,0) speye(prod(x.dim))]) & all(d==d(1))
     % Simple case x.^d
     y = vectorizedUnitPower(x,d);
     y.extra.createTime = definecreationtime;
     return
 end

% Back to scalar power...
d = d(1,1);
if x.dim(1)>1 | x.dim(2)>1
    switch d
        case 0
            y = 1;
        case 1
            y = x;
        otherwise
            y = x.*power(x,d-1);
    end
else
    error('This should not appen. Report bug (power does not use mpower)')
end

function y = vectorizedUnitPower(x,d)
d = d(1);
[mt,variabletype,hashM,hash] = yalmip('monomtable');
var = getvariables(x);

usedmt = mt(var,:);
newmt = usedmt*d;
hashV = newmt*hash;
if ~any(ismember(hashV,hashM))
    variabletype = [variabletype newvariabletypegen(newmt)];
    y = size(mt,1) + (1:length(var));
    mt = [mt;newmt];
else
    y = [];
    allnewmt = [];
    newvariables = 0;
    keep = zeros(size(usedmt,1),1);
    for i = 1:length(hashV)
        previous_var = find(abs(hashM - hashV(i)) < 1e-20);
        if isempty(previous_var)
            newmt =  usedmt(i,:)*d;
            variabletype = [variabletype newvariabletypegen(newmt)];
            newvariables = newvariables + 1;
            keep(i) = 1;
            y = [y size(mt,1)+newvariables];
        else
            y = [y previous_var];
        end
    end
    mt = [mt;usedmt(find(keep),:)*d];
end
yalmip('setmonomtable',mt,variabletype);
y = reshape(recover(y),x.dim);


text:
function [ u, v ] = uv_spiral ( n, x, y, c )

%*****************************************************************************80
%
%% UV_SPIRAL computes a spiral velocity vector field.
%
%  Discussion:
%
%    Note that the continuous velocity field (U,V)(X,Y) that is discretely
%    sampled here satisfies the homogeneous continuity equation, that is,
%    it has zero divergence.  In other words:
%
%      dU/dX + dV/dY = 0.
%
%    This is by construction, since we have
%
%      U(X,Y) =  10 * d/dY ( PHI(X) * PHI(Y) )
%      V(X,Y) = -10 * d/dX ( PHI(X) * PHI(Y) )
%
%    which guarantees zero divergence.
%
%    The underlying function PHI is defined by
%
%      PHI(Z) = ( 1 - cos ( C * pi * Z ) ) * ( 1 - Z )^2
%
%    where C is a parameter.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 January 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real X(N), Y(N), the coordinates of the evaluation points.
%
%    Input, real C, a parameter, typically between 0 and 2 * PI.
%
%    Output, real U(N), V(N), the velocity components.
%
  u = zeros ( n, 1 );
  v = zeros ( n, 1 );

  u =   10.0 * ( 1.0 - cos ( c * pi * x ) ) ...
           .* ( 1.0 - x ).^2 ...
           .* ( ...
               c * pi * sin ( c * pi * y ) .* ( 1.0 - y ).^2 ...
             - ( 1.0 - cos ( c * pi * y ) ) ...
             .* 2.0 .* ( 1.0 - y ) ...
             );

  v = - 10.0 * ( 1.0 - cos ( c * pi * y ) ) ...
           .* ( 1.0 - y ).^2 ...
           .* ( ...
               c * pi * sin ( c * pi * x ) .* ( 1.0 - x ).^2 ...
             - ( 1.0 - cos ( c * pi * x ) ) ...
             .* 2.0 .* ( 1.0 - x ) ...
             );

  return
end


text:
clear all

dt = 0.2;
t  = 0:dt:10;

Nsamples = length(t);

Avgsaved = zeros(Nsamples, 1);
Xmsaved  = zeros(Nsamples, 1);

for k=1:Nsamples
  xm  = GetVolt();
  avg = AvgFilter(xm);

  Avgsaved(k) = avg;
  Xmsaved(k)  = xm;
end

figure
plot(t, Xmsaved, 'r:*')
hold on
plot(t, Avgsaved, 'o-')

author: philbooks

repo: Kalman-Filter-for-Beginners


text:
classdef FMS_Cmd < Simulink.IntEnumType
    enumeration
        None(0)
        PreArm(1000)
        Arm(1001)
        Disarm(1002)
        Takeoff(1003)
        Land(1004)
        Return(1005)
        Pause(1006)
        Continue(1007)
    end

    methods (Static)

        function defaultValue = getDefaultValue()
            % GETDEFAULTVALUE  Returns the default enumerated value.
            %   If this method is not defined, the first enumeration is used.
            defaultValue = FMS_Cmd.None;
        end

        function dScope = getDataScope()
            % GETDATASCOPE  Specifies whether the data type definition should be imported from,
            %               or exported to, a header file during code generation.
            dScope = 'Auto';
        end

        function desc = getDescription()
            % GETDESCRIPTION  Returns a description of the enumeration.
            desc = 'enumeration of FMS command';
        end

        function fileName = getHeaderFile()
            % GETHEADERFILE  Returns path to header file if non-empty.
            fileName = '';
        end

        function flag = addClassNameToEnumNames()
            % ADDCLASSNAMETOENUMNAMES  Indicate whether code generator applies the class name as a prefix
            %                          to the enumeration.
            flag = true;
        end

    end
end

author: Firmament-Autopilot

repo: FMT-Model


text:
function point_total_num = sparse_grid_mixed_size_total ( dim_num, level_max, rule )

%*****************************************************************************80
%
%% SPARSE_GRID_MIXED_SIZE sizes a sparse grid, counting duplicate points.
%
%  Discussion:
%
%    The sparse grid is the logical sum of product grids with total LEVEL
%    between LEVEL_MIN and LEVEL_MAX.
%
%    In some cases, the same point may occur in different product grids
%    used to form the sparse grid.
%
%    This routine counts the total number of points used to construct the sparse
%    grid; if the same point occurs several times, each occurrence is added
%    to the sum.
%
%    This computation is useful in order to be able to allocate enough
%    space for the full set of points, before they are compressed by removing
%    duplicates.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 March 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Fabio Nobile, Raul Tempone, Clayton Webster,
%    A Sparse Grid Stochastic Collocation Method for Partial Differential
%    Equations with Random Input Data,
%    SIAM Journal on Numerical Analysis,
%    Volume 46, Number 5, 2008, pages 2309-2345.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer LEVEL_MAX, the maximum value of LEVEL.
%
%    Input, integer RULE(DIM_NUM), the rule in each dimension.
%     1, "CC",  Clenshaw Curtis, Closed Fully Nested rule.
%     2, "F2",  Fejer Type 2, Open Fully Nested rule.
%     3, "GP",  Gauss Patterson, Open Fully Nested rule.
%     4, "GL",  Gauss Legendre, Open Weakly Nested rule.
%     5, "GH",  Gauss Hermite, Open Weakly Nested rule.
%     6, "GGH", Generalized Gauss Hermite, Open Weakly Nested rule.
%     7, "LG",  Gauss Laguerre, Open Non Nested rule.
%     8, "GLG", Generalized Gauss Laguerre, Open Non Nested rule.
%     9, "GJ",  Gauss Jacobi, Open Non Nested rule.
%    10, "GW",  Golub Welsch, (presumed) Open Non Nested rule.
%    11, "CC_SE", Clenshaw Curtis Slow Exponential, Closed Fully Nested rule.
%    12, "F2_SE", Fejer Type 2 Slow Exponential, Closed Fully Nested rule.
%    13, "GP_SE", Gauss Patterson Slow Exponential, Closed Fully Nested rule.
%    14, "CC_ME", Clenshaw Curtis Moderate Exponential, Closed Fully Nested rule.
%    15, "F2_ME", Fejer Type 2 Moderate Exponential, Closed Fully Nested rule.
%    16, "GP_ME", Gauss Patterson Moderate Exponential, Closed Fully Nested rule.
%    17, "CCN", Clenshaw Curtis Nested, Linear, Closed Fully Nested rule.
%
%    Output, integer POINT_TOTAL_NUM, the total number of points in the grid.
%

%
%  Special case.
%
  if ( level_max == 0 )
    point_total_num = 1;
    return
  end

  point_total_num = 0;
%
%  The outer loop generates values of LEVEL.
%
  level_min = max ( 0, level_max + 1 - dim_num );

  for level = level_min : level_max
%
%  The middle loop generates a GRID,
%  based on the next partition that adds up to LEVEL.
%
    level_1d = [];
    more_grids = 0;
    h = 0;
    t = 0;

    while ( 1 )

      [ level_1d, more_grids, h, t ] = comp_next ( level, dim_num, level_1d, ...
        more_grids, h, t );

      order_1d = level_to_order_default ( dim_num, level_1d, rule );

      point_total_num = point_total_num + prod ( order_1d(1:dim_num) );

      if ( ~more_grids )
        break
      end

    end

  end

  return
end


text:
function manoptADhelp()
% Automatic differentation (AD) in Manopt requires the following:
%   Matlab version R2021a or later.
%   Deep Learning Toolbox version 14.2 or later.
%
% First, read the documentation of manoptAD by typing:
%
%   help manoptAD
%
% Examples using AD can be found in the folder:
%
%   /manopt/autodiff/basic_examples_AD
%
% More sophisticated examples can also be found in
%
%   /examples
%
% by reading the comments inside the scripts.
%
% The comments here provide further information about how to use AD. These
% comments are necessary because of certain limitations both of the AD
% capabilities of the Deep Learning Toolbox and of Manopt's ability to
% handle certain delicate geometries, e.g., fixed-rank matrices, with AD.
%
% The basic usage of AD in Manopt goes as follows:
%
%   problem.M = ...; %call a factory to get a manifold structure
%   problem.cost = @(x) ...; % implement your cost function
%   problem = manoptAD(problem); % ask AD to figure out gradient + Hessian
%   ...; % call a solver on problem, e.g., x = trustregions(problem);
%
% The current implementation of AD is based on Matlab's Deep Learning
% Toolbox. The latter relies on dlarray to represent data (arrays).
% While this works very well, unfortunately, there are a few limitations.
% For example, certain functions do not support dlarray.
% See the following official website for a list of compatbile functions:
%
%   https://ch.mathworks.com/help/deeplearning/ug/list-of-functions-with-dlarray-support.html.
%
% Moreover, dlarray only supports complex variables starting with
% Matlab R2021b.
%
% To handle complex numbers in R2021a and earlier, and also to handle
% functions that are not supported by dlarray at this moment, Manopt
% provides a limited collection of backup functions which are compatible
% with AD. Explicitly, the left column below lists commonly used functions
% which are not supported by dlarray at this time, and the right column
% lists corresponding replacement functions that can be used:
%
%     trace                 ctrace
%     diag                  cdiag
%     triu                  ctriu
%     tril                  ctril
%     norm(..., 'fro')      cnormfro
%     norm(..., 'fro')^2    cnormsqfro
%     multiscale            cmultiscale
%
% All the other multi*** functions in Manopt support AD.
%
% Moreover, bsxfun is not supported. The user may have to translate it
% into repmat and point-wise expressions.
% Concatenating arrays along the third or higher dimension such as in
% cat(A, B, 3+) is not supported for AD.
% The matrix functions sqrtm, logm, expm, eig, svd, det, cumsum,  movsum,
% cumprod, \, inv, mod, rem, vecnorm, bandwidth are not supported.
%
% For Matlab R2021a or earlier, in addition to the functions mentioned
% above, dot is not supported.
% Element-wise multiplication can be replaced by cdottimes.
%
% To deal with complex variables in R2021a and earlier, Manopt converts
% complex arrays into a structure with fields real and imag.
% See mat2dl_complex and dl2mat_complex for more details. In this case,
% the user should try using the basic functions in the folder
%
%   /manopt/autodiff/functions_AD
%
% when defining the cost function. An alternative way is to define one's
% own basic functions. These functions should accept both numeric arrays
% and structures with fields real and imag.
% See cprod and complex_example_AD as examples.
%
% See also: manoptAD complex_example_AD

% This file is part of Manopt: www.manopt.org.
% Original author: Xiaowen Jiang, Aug. 31, 2021.
% Contributors: Nicolas Boumal
% Change log:

    fprintf('The file manoptADhelp is just for documentation: run ''help manoptADhelp''.\n');

end

author: NicolasBoumal

repo: manopt


text:
function r83row_print_part ( n, a, max_print, title )

%*****************************************************************************80
%
%% R83ROW_PRINT_PART prints "part" of an R83ROW.
%
%  Discussion:
%
%    An R83ROW is a (3,N) array of R8's.
%
%    The user specifies MAX_PRINT, the maximum number of lines to print.
%
%    If N, the size of the vector, is no more than MAX_PRINT, then
%    the entire vector is printed, one entry per line.
%
%    Otherwise, if possible, the first MAX_PRINT-2 entries are printed,
%    followed by a line of periods suggesting an omission,
%    and the last entry.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 September 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries of the vector.
%
%    Input, real A(3,N), the vector to be printed.
%
%    Input, integer MAX_PRINT, the maximum number of lines
%    to print.
%
%    Input, string TITLE, a title.
%
  if ( max_print <= 0 )
    return
  end

  if ( n <= 0 )
    return
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '%s\n', title );
  fprintf ( 1, '\n' );

  if ( n <= max_print )

    for i = 1 : n
      fprintf ( 1, '  %4d  %14g  %14g  %14g\n', i, a(1:3,i) );
    end

  elseif ( 3 <= max_print )

    for i = 1 : max_print-2
      fprintf ( 1, '  %4d  %14g  %14g  %14g\n', i, a(1:3,i) );
    end
    fprintf ( 1, '  ....  ..............  ..............  ..............\n' );
    i = n;
    fprintf ( 1, '  %4d  %14g  %14g  %14g\n', i, a(1:3,i) );

  else

    for i = 1 : max_print - 1
      fprintf ( 1, '  %4d  %14g  %14g  %14g\n', i, a(1:3,i) );
    end
    i = max_print;
    fprintf ( 1, '  %4d  %14g  %14g  %14g  ...more entries...\n', i, a(1:3,i) );

  end

  return
end


text:
function area = sphere01_area ( )

%*****************************************************************************80
%
%% SPHERE01_AREA returns the area of the unit sphere.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 January 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real AREA, the area of the unit sphere.
%
  r = 1.0;
  area = 4.0 * pi * r * r;

  return
end


text:
function V = lumScotopic(imgXYZ)
%
%       V = lumScotopic(imgXYZ)
%
%       This function calculates the scotopic luminance
%
%       input:
%           img: an image in the XYZ color space
%
%       output:
%           V: scotopic luminance approximation by Larson, Rushmeier and Piatko 1997
%
%     Copyright (C) 2011-14  Francesco Banterle
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%

check3Color(imgXYZ);

eps = 1e-6;
t = (imgXYZ(:,:,2) + imgXYZ(:,:,3)) ./ (imgXYZ(:,:,1) + eps);
V = imgXYZ(:,:,2) .* (1.33 * (1.0 + t) - 1.68);

end

author: banterle


text:
classdef Tx < adi.AD9364.Tx
    % adi.FMComms4.Tx Transmit data from the FMComms4 evaluation board
    %   The adi.FMComms4.Tx System object is a signal sink that can
    %   transmit complex data from the FMComms4.
    %
    %   tx = adi.FMComms4.Tx;
    %   tx = adi.FMComms4.Tx('uri','192.168.2.1');
    %
    %   Product Page
    %
    %   See also adi.AD9364.Tx
    methods
        %% Constructor
        function obj = Tx(varargin)
            % Returns the matlabshared.libiio.base object
            coder.allowpcode('plain');
            obj = obj@adi.AD9364.Tx(varargin{:});
        end
    end

end

author: analogdevicesinc


text: plot(p,h,'LineWidth',2,'Color','b'),grid,ylabel('Altitude: h [m]'),xlabel('Pressure: p [Pa]')


text:
% prettyPolorScatter
%--------------------------------------------------------------------------
%
% plot with projection an area of the world
%
% POSSIBLE SINTAXES:
%   scatterHandler = prettyPolorScatter(dataScatter, phi, lambda);
%
%   prettyPolorScatter(dataScatter, phi, lambda, shape);
%   prettyPolorScatter(dataScatter, phi, lambda, projection);
%   prettyPolorScatter(dataScatter, phi, lambda, lineCol);
%
%   prettyPolorScatter(dataScatter, phi, lambda, phiMin, phiMax, lambdaMin, lambdaMax);
%   prettyPolorScatter(dataScatter, phi, lambda, phiMin, phiMax, lambdaMin, lambdaMax, shape);
%   prettyPolorScatter(dataScatter, phi, lambda, phiMin, phiMax, lambdaMin, lambdaMax, projection);
%   prettyPolorScatter(dataScatter, phi, lambda, phiMin, phiMax, lambdaMin, lambdaMax, lineCol);
%
%   prettyPolorScatter(dataScatter, phi, lambda, phiMin, phiMax, lambdaMin, lambdaMax, shape, projection);
%   prettyPolorScatter(dataScatter, phi, lambda, phiMin, phiMax, lambdaMin, lambdaMax, shape, lineCol);
%   prettyPolorScatter(dataScatter, phi, lambda, phiMin, phiMax, lambdaMin, lambdaMax, projection, shape);
%   prettyPolorScatter(dataScatter, phi, lambda, phiMin, phiMax, lambdaMin, lambdaMax, projection, lineCol);
%
%   prettyPolorScatter(dataScatter, phi, lambda, phiMin, phiMax, lambdaMin, lambdaMax, projection, shape, lineCol);
%   prettyPolorScatter(dataScatter, phi, lambda, phiMin, phiMax, lambdaMin, lambdaMax, shape, projection, lineCol);
%
% EXAMPLE:
%   prettyPolorScatter(dataScatter, phiMin, phiMax, lambdaMin, lambdaMax, 'Miller Cylindrical');
%
% INPUT:
%   dataScatter     array containing the data value (z = colour)
%   phi             array [degree]
%   lambda          array [degree]
%   phiMin          minimum latitude    [degree]
%   phiMax          maximum latitude    [degree]
%   lambdaMin       minimum longitude   [degree]
%   lambdaMax       maximum longitude   [degree]
%   projection      type of projection to be used "standard" is the default
%   shape           shapefile to load as coast (or country) contour
%                       - coast         only coasts coarse
%                       - 50m           1:50000000 scale country contours
%                       - 30m           1:30000000 scale country contours
%                       - 10m           1:10000000 scale country contours
%   lineCol         [1 1 1] array of RGB component to draw the contour lines
%
% DEFAULT VALUES:
%    projection = 'Lambert'
%
% AVAILABLE PROJECTION:
%    * Lambert
%      Stereographic
%      Orthographic
%      Azimuthal Equal-area
%      Azimuthal Equidistant
%      Gnomonic
%      Satellite
%      Albers Equal-Area Conic
%      Lambert Conformal Conic
%      Mercator
%    * Miller Cylindrical
%    * Equidistant Cylindrical (world dataScatter)
%      Oblique Mercator
%      Transverse Mercator
%      Sinusoidal
%      Gall-Peters
%      Hammer-Aitoff
%      Mollweide
%      Robinson
%    * UTM
%
% SEE ALSO:
%   dataScatterPlot, dataScatterPlot3D, scatter
%
% REQUIREMENTS:
%   M_Map: http://www.eos.ubc.ca/~rich/dataScatter.html
%   shape files with contours
%
% VERSION: 2.1
%
% CREDITS:
%   http://www.eos.ubc.ca/~rich/dataScatter.html
%
%   Andrea Gatti
%   DIIAR - Politecnico di Milano
%   2013-12-19
%

function scatterHandler = prettyPolorScatter(dataScatter, phi, lambda, phiMin, phiMax, lambdaMin, lambdaMax, projection, shape, lineCol)

%shape = 'coast';
%shape = '10m';
%shape = '30m';
%shape = '50m';

% lineCol = [0 0 0];

dataScatter = dataScatter(:);
phi = phi(:);
lambda = lambda(:);

limitsOk = false;

% Manage opening a new figure;
tohold = false;
if length(findall(0,'Type','figure'))>=1
    if ishold
        tohold = true;
    else
        figure;
    end
end

switch (nargin)
    case 3
        shape = 'coast';
        lineCol = [0 0 0];
        if (ischar(phi))
            projection = 'Miller Cylindrical';
            if (sum(strcmp(phi,[{'coast'},{'10m'},{'30m'},{'50m'}])))
                shape = phi;
                if (ischar(lambda))
                    projection = lambda;                              % prettyPolorScatter(dataScatter, shape, projection);
                else
                    lineCol = lambda;                                 % prettyPolorScatter(dataScatter, shape, lineCol);
                end
            else
                projection = phi;
                if (ischar(lambda))
                    shape = lambda;                                   % prettyPolorScatter(dataScatter, projection, shape);
                else
                    lineCol = lambda;                                 % prettyPolorScatter(dataScatter, projection, lineCol);
                end
            end
            phiMin = 90;
            phiMax = -90;
            lambdaMin = -180;
            lambdaMax = 180;

            deltaPhi = (phiMax-phiMin)/size(dataScatter,1);
            deltaLambda = (lambdaMax-lambdaMin)/size(dataScatter,2);

            phi    = (phiMin + deltaPhi/2 : deltaPhi : phiMax - deltaPhi/2)';
            lambda = (lambdaMin + deltaLambda/2 :  deltaLambda :  lambdaMax  - deltaLambda/2)';
        else                                                              % prettyPolorScatter(dataScatter, phi, lambda);
            projection = 'Miller Cylindrical';
            phiMin = max(phi);
            phiMax = min(phi);
            lambdaMin = min(lambda);
            lambdaMax = max(lambda);
        end
    case 4
        shape = 'coast';
        lineCol = [0 0 0];
        projection = 'Miller Cylindrical';
        if (ischar(phiMin))
            if (sum(strcmp(phiMin,[{'coast'},{'10m'},{'30m'},{'50m'}])))  % prettyPolorScatter(dataScatter, phi, lambda, shape);
                shape = phiMin;
            else                                                          % prettyPolorScatter(dataScatter, phi, lambda, projection);
                projection = phiMin;
            end
        elseif (length(phiMin) == 3)                                      % prettyPolorScatter(dataScatter, phi, lambda, lineCol);
            lineCol = phiMin;
        end

        phiMin = max(phi);
        phiMax = min(phi);
        lambdaMin = min(lambda);
        lambdaMax = max(lambda);
    case 5
        shape = 'coast';
        lineCol = [0 0 0];
        projection = 'Miller Cylindrical';
        if (ischar(phiMin))
            if (sum(strcmp(phiMin,[{'coast'},{'10m'},{'30m'},{'50m'}])))
                shape = phiMin;
                if (ischar(phiMax))
                    projection = phiMax;                                  % prettyPolorScatter(dataScatter, phiMin, phiMax, shape, projection);
                else
                    lineCol = phiMax;                                     % prettyPolorScatter(dataScatter, phiMin, phiMax, shape, lineCol);
                end
            else
                projection = phiMin;
                if (ischar(phiMax))
                    shape = phiMax;                                       % prettyPolorScatter(dataScatter, phiMin, phiMax, projection, shape);
                else
                    lineCol = phiMax;                                     % prettyPolorScatter(dataScatter, phiMin, phiMax, projection, lineCol);
                end
            end
            phiMin = max(phi);
            phiMax = min(phi);
            lambdaMin = min(lambda);
            lambdaMax = max(lambda);
        else                                                             %  prettyPolorScatter(dataScatter, phiMin, phiMax, lambdaMin, lambdaMax);
            limitsOk = true;
            lambdaMin = phiMin;
            lambdaMax = phiMax;
            phiMin = phi;
            phiMax = lambda;

            if (phiMin < phiMax)
                tmp = phiMin;
                phiMin = phiMax;
                phiMax = tmp;
            end

            deltaPhi = (phiMax-phiMin)/size(dataScatter,1);
            deltaLambda = (lambdaMax-lambdaMin)/size(dataScatter,2);

            phi    = (phiMin + deltaPhi/2 : deltaPhi : phiMax - deltaPhi/2)';
            lambda = (lambdaMin + deltaLambda/2 :  deltaLambda :  lambdaMax  - deltaLambda/2)';
        end
    case 6
        shape = 'coast';
        lineCol = [0 0 0];
        limitsOk = true;
        projection = 'Lambert';
        if (ischar(lambdaMin))
            if (sum(strcmp(lambdaMin,[{'coast'},{'10m'},{'30m'},{'50m'}])))
                shape = lambdaMin;                                        % prettyPolorScatter(dataScatter, phiMin, phiMax, lambdaMin, lambdaMax, shape);
            else
                projection = lambdaMin;
            end
        elseif (length(lambdaMin) == 3)
            lineCol = lambdaMin;                                          % prettyPolorScatter(dataScatter, phiMin, phiMax, lambdaMin, lambdaMax, lineCol);
        end

        lambdaMax = phiMax;
        lambdaMin = phiMin;
        phiMin = phi;
        phiMax = lambda;

        if (phiMin < phiMax)
            tmp = phiMin;
            phiMin = phiMax;
            phiMax = tmp;
        end

        deltaPhi = (phiMax-phiMin)/size(dataScatter,1);
        deltaLambda = (lambdaMax-lambdaMin)/size(dataScatter,2);

        phi    = (phiMin + deltaPhi/2 : deltaPhi : phiMax - deltaPhi/2)';
        lambda = (lambdaMin + deltaLambda/2 :  deltaLambda :  lambdaMax  - deltaLambda/2)';
    case 7
        shape = 'coast';
        lineCol = [0 0 0];
        limitsOk = true;
        projection = 'Lambert';
        if (ischar(lambdaMin))
            if (sum(strcmp(lambdaMin,[{'coast'},{'10m'},{'30m'},{'50m'}])))
                shape = lambdaMin;
                if (ischar(lambdaMax))
                    projection = lambdaMax;                               % prettyPolorScatter(dataScatter, phiMin, phiMax, lambdaMin, lambdaMax, shape, projection);
                else
                    lineCol = lambdaMax;                                  % prettyPolorScatter(dataScatter, phiMin, phiMax, lambdaMin, lambdaMax, shape, lineCol);
                end
            else
                projection = lambdaMin;
                if (ischar(lambdaMax))
                    shape = lambdaMax;                                    % prettyPolorScatter(dataScatter, phiMin, phiMax, lambdaMin, lambdaMax, projection, shape);
                else
                    lineCol = lambdaMax;                                  % prettyPolorScatter(dataScatter, phiMin, phiMax, lambdaMin, lambdaMax, projection, lineCol);
                end
            end

            lambdaMin = phiMin;
            lambdaMax = phiMax;
            phiMin = phi;
            phiMax = lambda;

            if (phiMin < phiMax)
                tmp = phiMin;
                phiMin = phiMax;
                phiMax = tmp;
            end

            deltaPhi = (phiMax-phiMin)/size(dataScatter,1);
            deltaLambda = (lambdaMax-lambdaMin)/size(dataScatter,2);

            phi    = (phiMin + deltaPhi/2 : deltaPhi : phiMax - deltaPhi/2)';
            lambda = (lambdaMin + deltaLambda/2 :  deltaLambda :  lambdaMax  - deltaLambda/2)';
        else
            projection = 'lambert';                                       % prettyPolorScatter(dataScatter, phi, lambda, phiMin, phiMax, lambdaMin, lambdaMax);
        end
    case 8
        shape = 'coast';
        lineCol = [0 0 0];
        limitsOk = true;
        if (ischar(projection))
            if (sum(strcmp(projection,[{'coast'},{'10m'},{'30m'},{'50m'}])))
                shape = projection;                                       % prettyPolorScatter(dataScatter, phi, lambda, phiMin, phiMax, lambdaMin, lambdaMax, shape);
                projection = 'Lambert';
            else
                % prettyPolorScatter(dataScatter, phi, lambda, phiMin, phiMax, lambdaMin, lambdaMax, projection);
            end
        elseif (length(projection) == 3)
            lineCol = projection;                                         % prettyPolorScatter(dataScatter, phi, lambda, phiMin, phiMax, lambdaMin, lambdaMax, lineCol);
            projection = 'Lambert';
        end
    case 9
        lineCol = [0 0 0];
        limitsOk = true;
        if (ischar(projection))
            if (sum(strcmp(projection,[{'coast'},{'10m'},{'30m'},{'50m'}])))
                tmp = shape;
                shape = projection;
                if (ischar(tmp))
                    projection = tmp;                                     % prettyPolorScatter(dataScatter, phi, lambda, phiMin, phiMax, lambdaMin, lambdaMax, shape, projection);
                else
                    lineCol = tmp;                                        % prettyPolorScatter(dataScatter, phi, lambda, phiMin, phiMax, lambdaMin, lambdaMax, shape, lineCol);
                    projection = 'UTM';
                end
            else
                if (ischar(shape))
                    % prettyPolorScatter(dataScatter, phi, lambda, phiMin, phiMax, lambdaMin, lambdaMax, projection, shape);
                else
                    lineCol = shape;                                      % prettyPolorScatter(dataScatter, phi, lambda, phiMin, phiMax, lambdaMin, lambdaMax, projection, lineCol);
                    shape = 'coast';
                end
            end
        end
    case 10                                                               % prettyPolorScatter(dataScatter, phi, lambda, phiMin, phiMax, lambdaMin, lambdaMax, projection, shape, lineCol)
        limitsOk = true;
        if (sum(strcmp(projection,[{'coast'},{'10m'},{'30m'},{'50m'}])))   % prettyPolorScatter(dataScatter, phi, lambda, phiMin, phiMax, lambdaMin, lambdaMax, shape, projection, lineCol)
            tmp = shape;
            shape = projection;
            projection = tmp;
        end
end

if (phiMin < phiMax)
    tmp = phiMin;
    phiMin = phiMax;
    phiMax = tmp;
end

lambdaTmp = sort(lambda);
[val, idMax] = max(diff(lambdaTmp));
if (sum(diff(lambdaTmp) == val) == 1) && val > 10
    lambdaTmp(1:idMax) = lambdaTmp(1:idMax)+360;
    if ~limitsOk
        lambdaMax = lambdaTmp(idMax);
        lambdaMin = lambdaTmp(idMax+1);
    end
end

if(lambdaMax<lambdaMin)
    lambdaMax = lambdaMax+360;
end

if (~limitsOk)
    lambdaMin = lambdaMin-3;
    lambdaMax = lambdaMax+3;
    phiMin = min(phiMin+3,90);
    phiMax = max(phiMax-3,-90);
end

projection = 'Transverse Mercator';

% setup the projection
if ~tohold
    if (strcmpi(projection,{'lambert'}) && abs(phiMax==-phiMin))
        projection='Miller Cylindrical';
    end

    if (sum(strcmpi(projection,[{'lambert'},{'UTM'},{'Sinusoidal'},{'Transverse Mercator'},{'Oblique Mercator'},{'Miller Cylindrical'}])))
        m_proj(projection,'long',[lambdaMin lambdaMax],'lat',[phiMax phiMin]);
    else
        m_proj(projection);
    end

    % Printing projection
    fprintf('Using projection: %s\n', projection);

    if sum(diff(lambda)<-200)
        lambda(lambda<0)=lambda(lambda<0)+360;
    end
    % plot the dataScatter
    m_pcolor([lambdaMin lambdaMax],[phiMin phiMax], nan(2));
    % set the light
    shading flat;

end

ids = (lambda<lambdaMin) & (lambda < 0);
lambda(ids) = lambda(ids)+360;

ids = (lambda>lambdaMax);
lambda(ids) = lambda(ids)-360;

hold on;
[xlocal,ylocal] = m_ll2xy(lambda(:),phi(:));
scatterHandler = scatter(xlocal,ylocal,30,dataScatter,'filled'); % <========================= SCATTER function is here

% read shapefile
if (~strcmp(shape,'coast'))
    if (strcmp(shape,'10m'))
        M=m_shaperead('countries_10m');
    elseif (strcmp(shape,'30m'))
        M=m_shaperead('countries_30m');
    else
        M=m_shaperead('countries_50m');
    end
    [xMin,yMin] = m_ll2xy(lambdaMin,phiMin);
    [xMax,yMax] = m_ll2xy(lambdaMax,phiMax);
    for k=1:length(M.ncst)
        lamC = M.ncst{k}(:,1);
        ids = lamC < lambdaMin;
        lamC(ids) = lamC(ids) + 360;
        phiC = M.ncst{k}(:,2);
        [x,y] = m_ll2xy(lamC,phiC);
        if sum(~isnan(x))>1
            x(find(abs(diff(x))>=abs(xMax-xMin)*0.90)+1) = nan; % Romove lines that occupy more than th 90% of the plot
            line(x,y,'color', lineCol);
        end
    end;
else
    m_coast('line','color', lineCol);
end

m_grid('box','fancy','tickdir','in');
drawnow;
if (phiMax >= 0)
    h=get(gca,'Children'); hf = findobj(h(1:end),'Tag','m_grid_fancybox1');
    delete(h(18));
    delete(hf([1:5 8]));
else
    h=get(gca,'Children');
    hf = findobj(h(1:end),'Tag','m_grid_fancybox1');
    %delete(h(18));
    delete(hf([1:4 6 7]));
end

colorbar;

if tohold
    hold on;
else
    hold off;
end

author: goGPS-Project


% =========================================================================
% This code is part of the Matlab-based toolbox
% LagLDDDM - A Lagrangian Gauss--Newton--Krylov Solver for Mass- and
%                        Intensity-Preserving Diffeomorphic Image Registration
%
% For details and license info see
% - https://github.com/C4IR/FAIR.m/tree/master/add-ons/LagLDDMM
%
% 2D Multilevel Mass-Preserving LDDMM Example using stationary velocity
% field and diffusion regularizer as described in Sec. 4 of the paper:
%
% @article{MangRuthotto2017,
%   Title = {A {L}agrangian {G}auss--{N}ewton--{K}rylov solver for mass- and intensity-preserving diffeomorphic image registration},
%   Year = {2017},
%   Journal = {SIAM Journal on Scientific Computing},
%   Author = {A. Mang, L. Ruthotto},
% }
%
% =========================================================================
close all; clc; clear all;
setup3DmouseData;
alpha = [50 0];
lvl = 4;
pad = 0;
mV = @(m) m;
nt = 0;
N  = 2;
minLevel = 4;
maxLevel = 6;
imgModel('reset','imgModel','linearInterMex')
% 1) setup grid for velocities (padded)
omegaV = omega; omegaV(1:2:end) = omegaV(1:2:end)-pad;  omegaV(2:2:end) = omega(2:2:end)+pad;

% 2) setup regularizer (and decide for stationary or nonstationary velocity)
regularizer('reset','regularizer','mfDiffusionCC','alpha',alpha,'nt',nt,'HessianShift',1e-2); % nonstationary velocity

%%
plots = 1;
NPIRpara    = optPara('NPIR-GN');
NPIRpara.maxIter = 40;
NPIRpara.scheme = @GaussNewtonLDDMM;

[vc,~,wc,his] = MLLDDMM(ML,'minLevel',minLevel,'maxLevel',maxLevel,'omegaV',omegaV,...
    'mV',mV,'N',N,'parametric',0,'plots',plots,'NPIRpara',NPIRpara,'NPIRobj',@MPLDDMMobjFctn);
switch regularizer
    case {'mfDiffusionST','mfCurvatureST'}
        yInv = getTrafoFromInstationaryVelocityRK4(vc,getNodalGrid(omega,m),'omega',omegaV,'m',m,'N',N,'nt',nt,'tspan',[1 0]);
    case {'mfDiffusionCC','mfCurvatureCC'}
        yInv = getTrafoFromVelocityRK4(vc,getNodalGrid(omega,m),'omega',omegaV,'m',m,'N',N,'tspan',[1 0]);
        nt = 0;
end

%%
yc = getTrafoFromVelocityRK4(vc,getNodalGrid(omega,m),'omega',omegaV,'m',m,'nt',nt,'tspan',[1,0],'N',N);
Jac = geometry(yInv,m,'Jac','omega',omega);
Topt = linearInterMex(dataT,omega,center(yInv,m)) .* Jac;
D0  = distance(dataT(:),dataR(:),omega,m);
DOpt = distance(Topt(:),dataR(:),omega,m);

fig = figure(); clf;
fig.Name = sprintf('Results for %s',mfilename);

subplot(2,3,1);
viewImage(dataR,omega,m);
title('reference');

subplot(2,3,4);
viewImage(dataT,omega,m);
title('template');

subplot(2,3,2);
viewImage(Topt,omega,m);
title('T(yc)')

subplot(2,3,3);
imgmontage(dataT(:)-dataR(:),omega,m);
title('init. residual, SSD=100%');

subplot(2,3,5);
imgmontage(Jac,omega,m);
title(sprintf('Jac, min=%1.2f max=%1.2f',min(Jac(:)),max(Jac(:))));

subplot(2,3,6);
imgmontage(Topt(:)-dataR(:),omega,m);
title(sprintf('opt residual, SSD=%1.2f%%',100*DOpt/D0));


text:
function [N] = kN2N(kN)
% Convert force from kilonewtons to newtons.
% Chad A. Greene 2012
N = kN*1000;


text:
function tf = isVarInLaunchVehicle(var, lvdData)
    lv = lvdData.launchVehicle;
    [lvVars,~] = lv.getActiveOptVars();

    tf = false;
    for(i=1:length(lvVars))
        lvVar = lvVars(i);

        if(var == lvVar)
            tf = true;
            break;
        end
    end
end


text:
function [coarse_approximations,prediction_errors]=gsp_pyramid_analysis_old(signal,Gs,num_levels,varargin)
%GSP_PYRAMID_ANALYSIS Compute the graph pyramid transform coefficients
%   Usage:  [coarse_approximations,prediction_errors]=gsp_pyramid_analysis(signal,Gs,num_levels);
%           [coarse_approximations,prediction_errors]=gsp_pyramid_analysis(signal,Gs,num_levels,param);
%
%   Input parameters:
%         signal                    : Graph signal to analyze.
%         Gs                        : A multiresolution sequence of graph structures, including the idx parameters tracking the subsampling pattern.
%         num_levels                : Number of levels in the pyramid transform.
%   Output parameters:
%         coarse_approximations     : Cell array with the coarse approximations at each level.
%         prediction_errors         : Cell array with the prediction errors at each level.
%   Additional parameters:
%         param.use_exact           : To use exact graph spectral filtering instead of the Chebyshev approximation.
%         param.order               : Degree of the Chebyshev approximation (default=30).
%         param.regularize_epsilon  : Interpolation parameter.
%         param.h_filters           : A cell array of graph spectral filters. If just one filter is included, it is used at every level of the pyramid.
%
%   'gsp__pyramid_analysis(signal,Gs,num_levels)' computes
%   the graph pyramid transform coefficients of a signal $f$.
%
%   See also:
%
%   Demos:
%
%   References:

%   AUTHOR : David I Shuman.
%   TESTING:
%   REFERENCE:

% Read input parameters and check that inputs have the correct sizes
if nargin>3
    param=varargin{1};
else
    param=0;
end

if length(signal) ~= Gs{1}.N
    error('The signal to analyze should have the same dimension as the first graph');
end

if num_levels >= length(Gs)
    error('Not enough graphs provided to compute that many levels of the graph Laplacian pyramid');
end

if ~isfield(param,'h_filters')
    h_filters=cell(num_levels,1);
    for i=1:num_levels
        h_filters{i}=@(x) .5./(.5+x);
    end
elseif length(param.h_filters)==1
    h_filters=cell(num_levels,1);
    for i=1:num_levels
        h_filters{i}=param.h_filters;
    end
elseif length(param.h_filters)==num_levels
    h_filters=param.h_filters;
else
    error('param.h_filters should be a cell array of length 1 or num_levels');
end

% Compute the pyramid transform
coarse_approximations=cell(num_levels+1,1);
coarse_approximations{1}=signal;
prediction_errors=cell(num_levels,1);

for i=1:num_levels
    [coarse_approximations{i+1},prediction_errors{i}]=gsp_pyramid_analysis_single_interpolation_old(coarse_approximations{i},Gs{i},Gs{i+1}.idx,h_filters{i},param);
end

repo: gspbox


text:
function [curr_segments] = generate_mincut_segments(gp_obj)
% Top level function which performs parametric min-cut, once the pairwise
% and unary capacities have been computed. The main job of this function is
% to call the right function according to which method needs to be used for
% parametric min-cut, and then collate the results once that function
% returns. gp_obj is a GraphProcessor object
%
% @authors:     Ahmad Humayun
% @contact:     ahumayun@cc.gatech.edu
% @affiliation: Georgia Institute of Technology
% @date:        Fall 2013 - Summer 2014

    t_pmc_full = tic;

    fprintf('\tComputing min st cut for all unaries ... ');

    pmc_maxflow_method = gp_obj.segm_params.pmc_maxflow_method;

    if ~isempty(strfind(pmc_maxflow_method, 'nodynamic')) || ...
       ~isempty(strfind(pmc_maxflow_method, 'kohli')) || ...
       ~isempty(strfind(pmc_maxflow_method, 'multiseed'))

        num_graph_types = length(gp_obj.graph_sets_per_method);

        start_ind = cumsum([1 gp_obj.graph_sets_per_method]);
        sols_to_unary_mapping = repmat({[]}, 1, num_graph_types);

        graphcut_params = [];
        if ~isempty(gp_obj.graph_sub_methods_cut_param)
            graphcut_params = gp_obj.graph_sub_methods_cut_param;
        end

        t_pmc_total_all = tic;

        [partitions, lambdas, seed_mapping, mincut_vals, t_pmc_all, ...
         cut_info] = gp_obj.multiseed_param_min_st_cut(...
                        gp_obj.segm_params.pmc_maxflow_method, ...
                        start_ind(1:end-1), graphcut_params);

        t_pmc_total_all = toc(t_pmc_total_all);
        t_pmc_total_all = repmat(t_pmc_total_all / num_graph_types, ...
                                 num_graph_types, 1);

        segs_meta_info.lambdas = lambdas;
        segs_meta_info.mincut_vals = mincut_vals;
        segs_meta_info.extra_cut_info = cut_info;

        curr_segments.cut_segs = partitions;
        curr_segments.segs_meta_info.sols_to_unary_mapping = seed_mapping;

    elseif strcmp(pmc_maxflow_method, 'hochbaum')
        num_graphs = size(gp_obj.graph_unaries_all.lambda_s,2);

        t_pmc_all = nan(num_graphs, 1);
        t_pmc_total_all = nan(num_graphs, 1);
        num_segs = nan(num_graphs, 1);

        segs_meta_info = repmat(struct, 1, num_graphs);
        partitions = repmat({[]}, 1, num_graphs);

        pmc_parallel_loop = tic;
        parfor(graph_idx = 1:num_graphs, 8)
%         for graph_idx = 1:num_graphs
            t_pmc_total = tic;

            [part, lambdas, mincut_vals, t_pmc] = ...
                gp_obj.parametric_min_st_cut(graph_idx);

            partitions{graph_idx} = ~part;
            segs_meta_info(graph_idx).lambdas = lambdas;
            segs_meta_info(graph_idx).mincut_vals = mincut_vals;
            segs_meta_info(graph_idx).extra_cut_info = [];

            num_segs(graph_idx) = size(part,2);
            t_pmc_all(graph_idx) = t_pmc;
            t_pmc_total_all(graph_idx) = toc(t_pmc_total);
        end
        pmc_parallel_loop = toc(pmc_parallel_loop);

        submethod_inds = ...
            duplicateElems(1:length(gp_obj.graph_sets_per_method), ...
                           gp_obj.graph_sets_per_method);
        parallel_fact = pmc_parallel_loop / sum(t_pmc_total_all);
        t_pmc_all = accumarray(submethod_inds(:), t_pmc_all) .* parallel_fact;
        t_pmc_total_all = accumarray(submethod_inds(:), t_pmc_total_all) ...
                          .* parallel_fact;

        curr_cut_segs = cell2mat(partitions);

        curr_segments.cut_segs = curr_cut_segs(1:end-2,:);
        curr_segments.segs_meta_info.sols_to_unary_mapping = ...
                        duplicateElems(1:length(num_segs), num_segs);
    else
        error('GraphProcessor:generate_mincut_segments', ...
              'Invalid min cut method: %s', pmc_maxflow_method);
    end

    % if debug, print segments overlayed over the image
    graph_num = length(gp_obj.seg_obj.num_segs.init_segs) + 1;
    diagnostic_methods('gen_masks_overlay', gp_obj.seg_obj, ...
                       curr_segments, ...
                       fullfile(sprintf('cuts%02d',graph_num), ...
                                '%03d_%04d.png'));

    total_segs = size(curr_segments.cut_segs,2);

    gp_obj.seg_obj.num_segs.init_segs = ...
        [gp_obj.seg_obj.num_segs.init_segs, total_segs];

    curr_segments.segs_meta_info.lambdas = [segs_meta_info.lambdas];
    curr_segments.segs_meta_info.mincut_vals = [segs_meta_info.mincut_vals];
    curr_segments.segs_meta_info.extra_cut_info = [segs_meta_info.extra_cut_info];
    curr_segments.segs_meta_info.seg_mapping_final_to_orig = 1:total_segs;

    time_util(gp_obj.seg_obj, 'pmc_time', t_pmc_full, 1, 1);
    fprintf('\t\t produced %d segments\n', ...
            gp_obj.seg_obj.num_segs.init_segs(end));

    gp_obj.seg_obj.timings.pmc_parallel_cut_time = ...
        [gp_obj.seg_obj.timings.pmc_parallel_cut_time, t_pmc_all'];
    gp_obj.seg_obj.timings.pmc_parallel_overhead_time = ...
        [gp_obj.seg_obj.timings.pmc_parallel_overhead_time, ...
         t_pmc_total_all' - t_pmc_all'];
end