text:
function [position0,position,messages]= delineateP3D(heasig,w1,w2,w3,timenew,position0,position,intervalo,samp,ultimo_anot,messages)
%
% [position0,position,A,V,A2,V2]=% delineateP3D(heasig,w1,w2,w3,intreg,scale,timenew,position0,position,intervalo,samp,figuresoff,sig)
%
% Input Parameters:
%   heasig: struct vector with header information
%   w1: matrix with WT scales 1 to 5 for lead 1
%   w2: matrix with WT scales 1 to 5 for lead 2
%   w3: matrix with WT scales 1 to 5 for lead 1
%  intreg: interval to be considered in fitting the optimal direction
%  scale to use in delineation
%   timenew:  QRS times in inedexes refering the interval included in the current excerpt (borders excluded)
%   position0: struct vector with the detected points for 3D multilead  step1
%   position: struct vector with the detected points for 3D multilead  step2
%   position1: struct vector with the detected points for lead1 - to be replaced for multilead marks
%   samp: samples included in the current excerpt (borders excluded)
%   intervalo: numeration of the beats processed in this segment
%
%Output Parameters:
%   actualized parameters: position, position0
% A - points of fitted lines in step 1
% A2 - points of fitted lines in step 2
%  V- vectors director to the fitted line
%  V2- vectors director to the fitted line
%
% Rute Almeida: trabalho de doutoramento
% Last update: 07FEB2012 Rute Almeida
%
% MATLAB Version 6.5.0.180913a (R13)

if nargin<11
    messages = new_empty_mesg;
elseif nargin<10
    messages.errors=[messages.errors {'Fatal error in wavedet_3D\delineate3D: not enough inputs.'}];
    warning(char(messages.errors(end)))
    messages.errors_desc=[messages.errors_desc 'Mandatory inputs not defined.'];
    messages.status=0;
    return
end
if ~isfield(messages,'setup'), messages.setup.wavedet=[]; end
if ~isfield(messages,'errors'), messages.errors=[]; end
if ~isfield(messages,'errors_desc'), messages.errors_desc=[]; end
if ~isfield(messages,'warnings'), messages.warnings=[]; end
if isfield(messages,'status')
    if messages.status~=1
        messages.warnings=[messages.warnings {['Initial status=' num2str(messages.status) '. status changed to 1']}];
    end
end
messages.status=1;

global recursioncount OPT

%Threshols for onset and offset detection
%%%%%% Constants and Thresholds  !!!!!!!!!!!!!!!!!!!!!!!!!
if ~isfield(messages.setup.wavedet,'inivent_tol_P')
    messages.setup.wavedet.inivent_tol_P  =  0.34;
end
if ~isfield(messages.setup.wavedet,'finvent_tol_P')
    messages.setup.wavedet.finvent_tol_P  =  0.1;
end
if ~isfield(messages.setup.wavedet,'P_CSE_tol')
    messages.setup.wavedet.P_CSE_tol=10.2*2/1000; % based on CSE tolerance
end
if ~isfield(messages.setup.wavedet,'Pconvergence_crit')
    messages.setup.wavedet.Pconvergence_crit=0; % 1 maks are acepted as the same if they difer by Tconvergence_crit
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

A2=[];
if ~isempty(w3)
    V=NaN*ones(size(timenew,2),3);
    A=NaN*ones(size(timenew,2),3);
else
    V=NaN*ones(size(timenew,2),2);
    A=NaN*ones(size(timenew,2),2);
end
V2=[];
intervalo=intervalo(1):intervalo(end);
for i=1:size(timenew,2)

    if ~isnan(position.QRSon(i+intervalo(1)-1))
        qrson = position.QRSon(i+intervalo(1)-1)-samp(1)+1;
    else
        qrson = timenew(i) - messages.setup.wavedet.finvent_tol_P*messages.setup.wavedet.freq;
    end

    ultimo_anot = ultimo_anot - samp(1)+1;
    % last anotation of the previous segment can overlap with the present one
    if (i==1 && ultimo_anot>=1),   % Protection in case it is the first beat
        begwin = max([1 qrson-round(messages.setup.wavedet.inivent_tol_P*messages.setup.wavedet.freq)+1 ultimo_anot+1]);
    else
        begwin = max(1, qrson-round(messages.setup.wavedet.inivent_tol_P*messages.setup.wavedet.freq)+1);
    end
    endwin = qrson- round(messages.setup.wavedet.finvent_tol_P*messages.setup.wavedet.freq)+1;

    if begwin>=endwin % empty search window: ultimo_anot should be too close to qrson
        position.Ptipoon(intervalo(i))=NaN;
        position.Ppicon(intervalo(i))=NaN;
        position.PscalePon(intervalo(i))=NaN;
        position.Pon(intervalo(i))=NaN;
        recursioncount.Pon(intervalo(i))=NaN;
        position.contadorPon(intervalo(i))=NaN;
        position.PscalePoff(intervalo(i))=NaN;
        position.Ppicoff(intervalo(i))=NaN;
        position.Poff(intervalo(i))=NaN;
        position.Ptipooff(intervalo(i))=NaN;
        recursioncount.Poff(intervalo(i))=NaN;
        position.contadorPoff(intervalo(i))=NaN; %Rute
        position.P(intervalo(i))=NaN;
    else

        scale=4;
        if ~isempty(w3)
            pontos=[w1(begwin:min(endwin,length(w1)),scale) w2(begwin:min(endwin,length(w2)),scale) w3(begwin:min(endwin,length(w3)),scale)];
        else
            pontos=[w1(begwin:min(endwin,length(w1)),scale) w2(begwin:min(endwin,length(w2)),scale)];
        end
        weight=ones(size(pontos,1),1);
        [a,v] = optimline(pontos, weight,OPT);
        A(i,:)=a;
        V(i,:)=v;
        %WT projected: from the previous QRS complex to the following QRS complex
        newleadbeat=[];
        if ~isempty(w3)
            if i>1,
                interval=timenew(i-1):timenew(i);
            else                        % if first beat of the segment
                interval=1:timenew(i);
            end

            newleadbeat(interval,3)=(w1(interval,3).*v(1)+w2(interval,3)*v(2)+w3(interval,3)*v(3))./norm(v); %#ok<AGROW>
            newleadbeat(interval,4)=(w1(interval,4).*v(1)+w2(interval,4)*v(2)+w3(interval,4)*v(3))./norm(v); %#ok<AGROW>
            newleadbeat(interval,5)=(w1(interval,5).*v(1)+w2(interval,5)*v(2)+w3(interval,5)*v(3))./norm(v); %#ok<AGROW>
            %signew(interval)=(sig(interval,1).*v(1)+sig(interval,2)*v(2)+sig(interval,3)*v(3))./norm(v);
        else
            if i>1,
                interval=timenew(i-1):timenew(i);
            else                        % if first beat of the segment
                interval=1:timenew(i);
            end
            % newleadbeat=NaN*ones(length(interval),5);
            newleadbeat(:,3)=(w1(interval,3).*v(1)+w2(timenew(i):timenew(i+1),3)*v(2))./norm(v);
            newleadbeat(:,4)=(w1(interval,4).*v(1)+w2(timenew(i):timenew(i+1),4)*v(2))./norm(v);
            newleadbeat(:,5)=(w1(interval,5).*v(1)+w2(timenew(i):timenew(i+1),5)*v(2))./norm(v);
            %signew=(sig(timenew(i):timenew(i+1),1).*v(1)+sig(timenew(i):timenew(i+1),2)*v(2))./norm(v);
        end
        if i>1,
            pos.qrs=position.qrs([i-1 i]);
            pos.QRSon=position.QRSon([i-1 i]);
            pos.QRSoff=position.QRSoff([i-1 i]);
            pos.Toff=position.Toff([i-1 i]);
            [pos,picon,picoff]= pwavef(heasig,samp,timenew([i i]),position,newleadbeat,[intervalo(i) intervalo(i)],ultimo_anot,messages);
        else
            pos.qrs=position.qrs(i);
            pos.QRSon=position.QRSon(i);
            pos.QRSoff=position.QRSoff(i);
            pos.Toff=position.Toff(i);
            [pos,picon,picoff]= pwavef(heasig,samp,timenew(i),position,newleadbeat,[intervalo(i) intervalo(i)] ,ultimo_anot,messages);
        end

        position0.Pscale(intervalo(i))=pos.Pscale(intervalo(i));
        position0.Pon(intervalo(i))=pos.Pon(intervalo(i));
        position0.Poff(intervalo(i))=pos.Poff(intervalo(i));
        position0.Ptipo(intervalo(i))=pos.Ptipo(intervalo(i));
        position0.direccionPonopt(:,intervalo(i))=NaN*ones(size(intervalo(i)));
        %position0.P(intervalo(i))=pos.P(intervalo(i));
        %position0.Pprima(intervalo(i))=pos.Pprima(intervalo(i));

        %%%%% P wave onset
        if ~isnan(picon)

            begaux=max([(position0.Pon(intervalo(i))-samp(1)+1-round(messages.setup.wavedet.P_CSE_tol*messages.setup.wavedet.freq)) 1]);
            if ~isempty(w3)
                pontos2on=[w1(begaux:picon,scale) w2(begaux:picon,scale) w3(begaux:picon,scale)];
            else
                pontos2on=[w1(begaux:picon,scale) w2(begaux:picon,scale)];
            end
            if size(pontos2on,1)<=1 % unable to do another step: empty search window
                position.contadorPon(intervalo(i))=-1;
                position.PscalePon(intervalo(i))=position0.Pscale(intervalo(i));
                position.Pon(intervalo(i))=position0.Pon(intervalo(i));
                position.Ppicon=picon+samp(1)-1;
                %position.P(intervalo(i))= position0.P(intervalo(i));
                %position.Pprima(intervalo(i))=position0.Pprima(intervalo(i));
                position.direccionPonopt(intervalo(i),:)=V(end,:)';
                position.Ptipoon(intervalo(i))=position0.Ptipo(intervalo(i));
            else
                weight2=ones(size(pontos2on,1),1);
                [a,v] = optimline(pontos2on, weight2,OPT);
                A2=[A2;a]; %#ok<AGROW>
                V2=[V2;v]; %#ok<AGROW>
                newleadbeat2on = [];
                if ~isempty(w3)
                    newleadbeat2on(interval,3)=(w1(interval,3).*v(1)+w2(interval,3)*v(2)+w3(interval,3)*v(3))./norm(v);%#ok<AGROW>
                    newleadbeat2on(interval,4)=(w1(interval,4).*v(1)+w2(interval,4)*v(2)+w3(interval,4)*v(3))./norm(v);%#ok<AGROW>
                    newleadbeat2on(interval,5)=(w1(interval,5).*v(1)+w2(interval,5)*v(2)+w3(interval,5)*v(3))./norm(v);%#ok<AGROW>
                    % signew2on(interval)=(sig(interval,1).*v(1)+sig(interval,2)*v(2)+sig(interval,3)*v(3))./norm(v);
                else
                    newleadbeat2on(interval,3)=(w1(interval,3).*v(1)+w2(timenew(i):timenew(i+1),3)*v(2))./norm(v);%#ok<AGROW>
                    newleadbeat2on(interval,4)=(w1(interval,4).*v(1)+w2(timenew(i):timenew(i+1),4)*v(2))./norm(v);%#ok<AGROW>
                    newleadbeat2on(interval,5)=(w1(interval,5).*v(1)+w2(timenew(i):timenew(i+1),5)*v(2))./norm(v);%#ok<AGROW>
                    %signew2on(interval)=(sig(timenew(i):timenew(i+1),1).*v(1)+sig(timenew(i):timenew(i+1),2)*v(2))./norm(v);
                end

                if i>1,
                    [pos2,picon2,picoff2]= pwavef(heasig,samp,timenew([i  i]),pos,newleadbeat2on,[intervalo(i) intervalo(i)],ultimo_anot,messages); %#ok<NASGU>
                else
                    [pos2,picon2,picoff2]= pwavef(heasig,samp,timenew(i),pos,newleadbeat2on,[intervalo(i) intervalo(i)] ,ultimo_anot,messages); %#ok<NASGU>
                end
                piconall=[picon picon2];
                if isnan(picon2)
                    amppiconall=[abs(newleadbeat(picon,scale)) NaN];
                else
                    amppiconall=[abs(newleadbeat(picon,scale)) abs(newleadbeat2on(picon2,scale))];
                end
                if (amppiconall(end)<amppiconall(end-1)) || isnan(amppiconall(end))% second step worse than first step
                    position.contadorPon(intervalo(i))=-1;
                    position.Ppicon(intervalo(i))=piconall(end-1)+samp(1)-1;
                    position.PscalePon(intervalo(i))=position0.Pscale(intervalo(i));
                    position.Pon(intervalo(i))=position0.Pon(intervalo(i));
                    position.direccionPonopt(:,intervalo(i))=V(end,:);
                    position.Ptipoon(intervalo(i))=position0.Ptipo(intervalo(i));
                elseif pos2.Pon(intervalo(i))==position0.Pon(intervalo(i)) % first and second steps give same mark
                    position.contadorPon(intervalo(i))=0;
                    position.Ppicon(intervalo(i))=piconall(end)+samp(1)-1;
                    position.PscalePon(intervalo(i))=pos2.Pscale(intervalo(i));
                    position.Pon(intervalo(i))=pos2.Pon(intervalo(i));
                    position.direccionPonopt(:,intervalo(i))=V2(end,:);
                    position.Ptipoon(intervalo(i))=pos2.Ptipo(intervalo(i));
                else % more steps
                    newPon=[position0.Pon(intervalo(i)) pos2.Pon(intervalo(i))];
                    position.contadorPon(intervalo(i))=0;
                    posnewon.Pscale(intervalo(i))=NaN;
                    posnewon.Pon(intervalo(i))=NaN;
                    posnewon.Ptipo(intervalo(i))=NaN;
                    position.Ppicon(intervalo(i))=piconall(end)+samp(1)-1;
                    while (~isnan(piconall(end)) &&  abs(newPon(end)-newPon(end-1))>messages.setup.wavedet.Pconvergence_crit) && (amppiconall(end)>amppiconall(end-1))
                        position.contadorPon(intervalo(i))=position.contadorPon(intervalo(i))+1;
                        begaux=max([(position0.Pon(intervalo(i))-samp(1)+1-round(messages.setup.wavedet.P_CSE_tol*messages.setup.wavedet.freq)) 1]);
                        if ~isempty(w3)
                            pontosnewon=[w1( begaux:piconall(end),scale) w2( begaux:piconall(end),scale) w3( begaux:piconall(end),scale)];
                        else
                            pontosnewon=[w1(begaux:piconall(end),scale) w2( begaux:piconall(end),scale)];
                        end
                        if size(pontosnewon,1)<=1 % unable to do another step: empty search window
                            posnewon.Pon(intervalo(i))=newPon(end-1); % previous step
                            posnewon.Pscale(intervalo(i))=pos2.Pscale(intervalo(i));
                            posnewon.Ptipo(intervalo(i))=pos2.Ptipo(intervalo(i));
                            position.contadorPon(intervalo(i))=position.contadorPon(intervalo(i))-1;
                            piconall=[piconall NaN]; %#ok<AGROW>
                            amppiconall=[amppiconall NaN]; %#ok<AGROW>
                        else
                            weight2=ones(size(pontosnewon,1),1);
                            [a,v] = optimline(pontosnewon, weight2,OPT);
                            A2=[A2;a]; %#ok<AGROW>
                            V2=[V2;v]; %#ok<AGROW>
                            newleadbeatnewon = [];
                            if ~isempty(w3)
                                % newleadbeatnewon=NaN*ones(interval(end),5);
                                newleadbeatnewon(interval,3)=(w1(interval,3).*v(1)+w2(interval,3)*v(2)+w3(interval,3)*v(3))./norm(v);%#ok<AGROW>
                                newleadbeatnewon(interval,4)=(w1(interval,4).*v(1)+w2(interval,4)*v(2)+w3(interval,4)*v(3))./norm(v);%#ok<AGROW>
                                newleadbeatnewon(interval,5)=(w1(interval,5).*v(1)+w2(interval,5)*v(2)+w3(interval,5)*v(3))./norm(v);%#ok<AGROW>
                                %signewon(interval)=(sig(interval,1).*v(1)+sig(interval,2)*v(2)+sig(interval,3)*v(3))./norm(v);
                            else
                                % newleadbeatnewon=NaN*ones(interval(end),5);
                                newleadbeatnewon(interval,3)=(w1(interval,3).*v(1)+w2(timenew(i):timenew(i+1),3)*v(2))./norm(v);%#ok<AGROW>
                                newleadbeatnewon(interval,4)=(w1(interval,4).*v(1)+w2(timenew(i):timenew(i+1),4)*v(2))./norm(v);%#ok<AGROW>
                                newleadbeatnewon(interval,5)=(w1(interval,5).*v(1)+w2(timenew(i):timenew(i+1),5)*v(2))./norm(v);%#ok<AGROW>
                                %signewon(interval)=(sig(timenew(i):timenew(i+1),1).*v(1)+sig(timenew(i):timenew(i+1),2)*v(2))./norm(v);
                            end
                            if i>1,
                                [posnewon,piconnew,picoffnew]= pwavef(heasig,samp,timenew([i  i]),pos,newleadbeatnewon,[intervalo(i) intervalo(i)],ultimo_anot,messages); %#ok<NASGU>
                            else
                                [posnewon,piconnew,picoffnew]= pwavef(heasig,samp,timenew(i),pos,newleadbeatnewon,[intervalo(i) intervalo(i)] ,ultimo_anot,messages); %#ok<NASGU>
                            end
                            newPon=[newPon posnewon.Pon(intervalo(i))]; %#ok<AGROW>
                            if isnan(piconnew)
                                amppiconall=[amppiconall NaN]; %#ok<AGROW>
                            else
                                amppiconall=[amppiconall abs(newleadbeatnewon(piconnew,scale))]; %#ok<AGROW>
                            end
                            piconall=[piconall piconnew]; %#ok<AGROW>

                            if isnan(posnewon.Pon(intervalo(i))) || (amppiconall(end)<amppiconall(end-1)) % lead got worse
                                posnewon.Pon(intervalo(i))=newPon(end-1); % previous step
                                position.Ppicon(intervalo(i))=piconall(end-1)+samp(1)-1;
                                posnewon.Pscale(intervalo(i))=pos2.Pscale(intervalo(i));
                                posnewon.Ptipo(intervalo(i))=pos2.Ptipo(intervalo(i));
                                position.contadorPon(intervalo(i))=position.contadorPon(intervalo(i))-1;
                            else
                                position.direccionPonopt(:,intervalo(i))=V2(end,:);
                                position.Ppicon(intervalo(i))=piconall(end)+samp(1)-1;
                            end
                            if  prod(piconall(1:end-1)-piconall(end))==0  % ciclo infinito no teste 14
                                piconall(end)=NaN;
                            end

                        end % if size(pontosnewon,1)>1
                    end %while
                    position.PscalePon(intervalo(i))=posnewon.Pscale(intervalo(i));
                    position.Pon(intervalo(i))=posnewon.Pon(intervalo(i));
                    position.Ptipoon(intervalo(i))=posnewon.Ptipo(intervalo(i));
                end % more steps
            end   %size(pontosnewon,1)>1
        else % isnan(picon)
            position.Ptipoon(intervalo(i))=NaN;
            position.Ppicon(intervalo(i))=NaN;
            position.PscalePon(intervalo(i))=NaN;
            position.Pon(intervalo(i))=NaN;
            recursioncount.Pon(intervalo(i))=NaN;
            position.contadorPon(intervalo(i))=NaN; %Rute
        end
        %%%%%% P wave onset
        %%%%% P wave end
        if ~isnan(picoff)

            endaux=min([(position0.Poff(intervalo(i))-samp(1)+1+round(messages.setup.wavedet.P_CSE_tol*messages.setup.wavedet.freq)) length(w1)]);
            if ~isempty(w3)
                pontos2off=[w1(picoff:endaux,scale) w2(picoff:endaux,scale) w3(picoff:endaux,scale)];
            else
                pontos2off=[w1(picoff:endaux,scale) w2(picoff:endaux,scale)];
            end
            if size(pontos2off,1)<=1 % unable to do another step: empty search window
                position.contadorPoff(intervalo(i))=-1;
                position.PscalePoff(intervalo(i))=position0.Pscale(intervalo(i));
                position.Poff(intervalo(i))=position0.Poff(intervalo(i));
                position.Ppicoff=picoff+samp(1)-1;
                %position.P(intervalo(i))= position0.P(intervalo(i));
                %position.Pprima(intervalo(i))=position0.Pprima(intervalo(i));
                position.direccionPoffopt(:,intervalo(i))=V(end,:);
                position.Ptipooff(intervalo(i))=position0.Ptipo(intervalo(i));
            else
                weight2=ones(size(pontos2off,1),1);
                [a,v] = optimline(pontos2off, weight2,OPT);
                A2=[A2;a]; %#ok<AGROW>
                V2=[V2;v]; %#ok<AGROW>
                newleadbeat2off = [];
                if ~isempty(w3)
                    %   newleadbeat2off=NaN*ones(interval(end),5);
                    newleadbeat2off(interval,3)=(w1(interval,3).*v(1)+w2(interval,3)*v(2)+w3(interval,3)*v(3))./norm(v);%#ok<AGROW>
                    newleadbeat2off(interval,4)=(w1(interval,4).*v(1)+w2(interval,4)*v(2)+w3(interval,4)*v(3))./norm(v);%#ok<AGROW>
                    newleadbeat2off(interval,5)=(w1(interval,5).*v(1)+w2(interval,5)*v(2)+w3(interval,5)*v(3))./norm(v);%#ok<AGROW>
                    %signew2off(interval)=(sig(interval,1).*v(1)+sig(interval,2)*v(2)+sig(interval,3)*v(3))./norm(v);
                else
                    % newleadbeat2off=NaN*ones(interval(end),5);
                    newleadbeat2off(interval,3)=(w1(interval,3).*v(1)+w2(timenew(i):timenew(i+1),3)*v(2))./norm(v);%#ok<AGROW>
                    newleadbeat2off(interval,4)=(w1(interval,4).*v(1)+w2(timenew(i):timenew(i+1),4)*v(2))./norm(v);%#ok<AGROW>
                    newleadbeat2off(interval,5)=(w1(interval,5).*v(1)+w2(timenew(i):timenew(i+1),5)*v(2))./norm(v);%#ok<AGROW>
                    %signew2off(interval)=(sig(timenew(i):timenew(i+1),1).*v(1)+sig(timenew(i):timenew(i+1),2)*v(2))./norm(v);
                end

                if i>1,
                    [pos2off,picon2off,picoff2]= pwavef(heasig,samp,timenew([i  i]),pos,newleadbeat2off,[intervalo(i) intervalo(i)],ultimo_anot,messages); %#ok<ASGLU>
                else
                    [pos2off,picon2off,picoff2]= pwavef(heasig,samp,timenew(i),pos,newleadbeat2off,[intervalo(i) intervalo(i)] ,ultimo_anot,messages); %#ok<ASGLU>
                end

                picoffall=[picoff picoff2];
                if isnan(picoff2)
                    amppicoffall=[abs(newleadbeat(picoff,scale)) NaN];
                else
                    amppicoffall=[abs(newleadbeat(picoff,scale)) abs(newleadbeat2off(picoff2,scale))];
                end
                if (amppicoffall(end)<amppicoffall(end-1))  || isnan(amppicoffall(end)) % second step worse than first step
                    position.contadorPoff(intervalo(i))=-1;
                    position.Ppicoff(intervalo(i))=picoffall(end-1)+samp(1)-1;
                    position.PscalePoff(intervalo(i))=position0.Pscale(intervalo(i));
                    position.Poff(intervalo(i))=position0.Poff(intervalo(i));
                    %position.P(intervalo(i))= position0.P(intervalo(i));
                    %position.Pprima(intervalo(i))=position0.Pprima(intervalo(i));
                    position.direccionPoffopt(:,intervalo(i))=V(end,:);
                    position.Ptipooff(intervalo(i))=position0.Ptipo(intervalo(i));
                elseif pos2off.Poff(intervalo(i))==position0.Poff(intervalo(i)) % first and second steps give same mark
                    position.contadorPoff(intervalo(i))=0;
                    position.Ppicoff(intervalo(i))=picoffall(end)+samp(1)-1;
                    position.PscalePoff(intervalo(i))=pos2off.Pscale(intervalo(i));
                    position.Poff(intervalo(i))=pos2off.Poff(intervalo(i));
                    %position.P(intervalo(i))= pos2off.P(intervalo(i));
                    %position.Pprima(intervalo(i))=pos2off.Pprima(intervalo(i));
                    position.direccionPoffopt(:,intervalo(i))=V2(end,:);
                    position.Ptipooff(intervalo(i))=pos2off.Ptipo(intervalo(i));
                else % more steps
                    newPoff=[position0.Poff(intervalo(i)) pos2off.Poff(intervalo(i))];
                    position.contadorPoff(intervalo(i))=0;
                    position.Ppicoff(intervalo(i))=picoffall(end)+samp(1)-1;
                    posnewoff.Pscale(intervalo(i))=NaN;
                    posnewoff.Poff(intervalo(i))=NaN;
                    posnewoff.Ptipo(intervalo(i))=NaN;
                    while (~isnan(picoffall(end)) &&  abs(newPoff(end)-newPoff(end-1))>messages.setup.wavedet.Pconvergence_crit ) && (amppicoffall(end)>amppicoffall(end-1))
                        position.contadorPoff(intervalo(i))=position.contadorPoff(intervalo(i))+1;
                        endaux=min([(position0.Poff(intervalo(i))-samp(1)+1+round(messages.setup.wavedet.P_CSE_tol*messages.setup.wavedet.freq)) length(w1)]);
                        if ~isempty(w3)
                            pontosnewoff=[w1(picoffall(end):endaux,scale) w2(picoffall(end):endaux,scale) w3(picoffall(end):endaux,scale)];
                        else
                            pontosnewoff=[w1(picoffall(end):endaux,scale) w2(picoffall(end):endaux,scale)];
                        end
                        if size(pontosnewoff,1)<=1 % unable to do another step: empty search window
                            posnewoff.Poff(intervalo(i))=newPoff(end-1); % previous step
                            posnewoff.Pscale(intervalo(i))=pos2off.Pscale(intervalo(i));
                            posnewoff.Ptipo(intervalo(i))=pos2off.Ptipo(intervalo(i));
                            %posnew.P=pos2off.P;
                            %posnew.Pprima=pos2off.Pprima;
                            position.contadorPoff(intervalo(i))=position.contadorPoff(intervalo(i))-1;
                            picoffall=[picoffall NaN]; %#ok<AGROW>
                            amppicoffall=[amppicoffall NaN]; %#ok<AGROW>
                        else
                            weight2=ones(size(pontosnewoff,1),1);
                            [a,v] = optimline(pontosnewoff, weight2,OPT);
                            A2=[A2;a]; %#ok<AGROW>
                            V2=[V2;v]; %#ok<AGROW>
                            newleadbeatnewoff = [];
                            if ~isempty(w3)
                                %newleadbeatnewoff=NaN*ones(interval(end),5);
                                newleadbeatnewoff(interval,3)=(w1(interval,3).*v(1)+w2(interval,3)*v(2)+w3(interval,3)*v(3))./norm(v);%#ok<AGROW>
                                newleadbeatnewoff(interval,4)=(w1(interval,4).*v(1)+w2(interval,4)*v(2)+w3(interval,4)*v(3))./norm(v);%#ok<AGROW>
                                newleadbeatnewoff(interval,5)=(w1(interval,5).*v(1)+w2(interval,5)*v(2)+w3(interval,5)*v(3))./norm(v);%#ok<AGROW>
                                %signewoff(interval)=(sig(interval,1).*v(1)+sig(interval,2)*v(2)+sig(interval,3)*v(3))./norm(v);
                            else
                                %newleadbeatnewoff=NaN*ones(interval(end),5);
                                newleadbeatnewoff(interval,3)=(w1(interval,3).*v(1)+w2(timenew(i):timenew(i+1),3)*v(2))./norm(v);%#ok<AGROW>
                                newleadbeatnewoff(interval,4)=(w1(interval,4).*v(1)+w2(timenew(i):timenew(i+1),4)*v(2))./norm(v);%#ok<AGROW>
                                newleadbeatnewoff(interval,5)=(w1(interval,5).*v(1)+w2(timenew(i):timenew(i+1),5)*v(2))./norm(v);%#ok<AGROW>
                                %signewoff(interval)=(sig(timenew(i):timenew(i+1),1).*v(1)+sig(timenew(i):timenew(i+1),2)*v(2))./norm(v);
                            end
                            if i>1,
                                [posnewoff,piconnew,picoffnew]= pwavef(heasig,samp,timenew([i  i]),pos,newleadbeatnewoff,[intervalo(i) intervalo(i)],ultimo_anot,messages); %#ok<ASGLU>
                            else
                                [posnewoff,piconnew,picoffnew]= pwavef(heasig,samp,timenew(i),pos,newleadbeatnewoff,[intervalo(i) intervalo(i)] ,ultimo_anot,messages); %#ok<ASGLU>
                            end
                            newPoff=[newPoff posnewoff.Poff(intervalo(i))]; %#ok<AGROW>
                            if isnan(picoffnew)
                                amppicoffall=[amppicoffall NaN]; %#ok<AGROW>
                            else
                                amppicoffall=[amppicoffall abs(newleadbeatnewoff(picoffnew,scale))]; %#ok<AGROW>
                            end
                            picoffall=[picoffall picoffnew]; %#ok<AGROW>

                            if isnan(posnewoff.Poff(intervalo(i))) || (amppicoffall(end)<amppicoffall(end-1)) % lead got worse
                                position.Ppicoff(intervalo(i))=picoffall(end-1)+samp(1)-1;
                                posnewoff.Poff(intervalo(i))=newPoff(end-1); % previous step
                                posnewoff.Pscale(intervalo(i))=pos2off.Pscale(intervalo(i));
                                posnewoff.Ptipo(intervalo(i))=pos2off.Ptipo(intervalo(i));
                                %posnewoff.P=pos2.P;
                                %posnewoff.Pprima=pos2.Pprima;
                                position.contadorPoff(intervalo(i))=position.contadorPoff(intervalo(i))-1;
                            else
                                position.direccionPoffopt(:,intervalo(i))=V2(end,:);
                                position.Ppicoff(intervalo(i))=picoffall(end)+samp(1)-1;
                            end

                            if  prod(picoffall(1:end-1)-picoffall(end))==0
                                picoffall(end)=NaN;
                            end
                        end % if size(pontosnewoff,1)>1
                    end %while
                    position.PscalePoff(intervalo(i))=posnewoff.Pscale(intervalo(i));
                    position.Poff(intervalo(i))=posnewoff.Poff(intervalo(i));
                    position.Ptipooff(intervalo(i))=posnewoff.Ptipo(intervalo(i));
                    % position.P(intervalo(i))= posnewoff.P(intervalo(i));
                    % position.Pprima(intervalo(i))=posnewoff.Pprima(intervalo(i));
                end % more steps
            end   %size(pontosnewoff,1)>1
        else % isnan(picoff)
            position.PscalePoff(intervalo(i))=NaN;
            position.Ppicoff(intervalo(i))=NaN;
            position.Poff(intervalo(i))=NaN;
            position.Ptipooff(intervalo(i))=NaN;
            %position.P(intervalo(i))=NaN;
            %position.Pprima(intervalo(i))=NaN;
            recursioncount.Poff(intervalo(i))=NaN;
            position.contadorPoff(intervalo(i))=NaN; %Rute
        end
        %%%%%% P wave end

        %P peak
        if ~isnan(position.Ppicon(intervalo(i))) && ~isnan(position.Ppicoff(intervalo(i))) && position.Ppicon(intervalo(i))<position.Ppicoff(intervalo(i))
            pontosP=[w1((position.Ppicon(intervalo(i))-samp(1)+1):(position.Ppicoff(intervalo(i))-samp(1)+1),scale) w2((position.Ppicon(intervalo(i))-samp(1)+1):(position.Ppicoff(intervalo(i))-samp(1)+1),scale) w3((position.Ppicon(intervalo(i))-samp(1)+1):(position.Ppicoff(intervalo(i))-samp(1)+1),scale)];
            [m1,m2]=min((pontosP(:,1)).^2+(pontosP(:,2)).^2+(pontosP(:,3)).^2); %#ok<ASGLU> % main P peak
            position.P(intervalo(i))=m2+position.Ppicon(intervalo(i));
        else
            position.P(intervalo(i))=position0.P(intervalo(i));
        end
        %%
        %writennot protection: at least one sample between peak and border
        if (position.P(intervalo(i))-position.Pon(intervalo(i)))<1
            position.P(intervalo(i))=NaN;
        end
        if (position.Poff(intervalo(i))-position.P(intervalo(i)))<1
            position.P(intervalo(i))=NaN;
        end
        clear pos pos2 pontos pontos2on pontos2off posnew pos2off posnewoff
    end
end % for

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function [ out ] = opt_stackHistory(in, c )
%OPT_STACKHISTORY Summary of this function goes here
%   Detailed explanation goes here
if isstruct(in)
    c = strsplit(c,'\');
    if isfield(in,'stack')
        in.stack{end+1}=c{end};
    end
else
    return;
end
out=in;
end

author: PatternRecognition

repo: OpenBMI


text:
function g = dnetGradient(params, model)

% DNETGRADIENT Density Network gradient wrapper.
% FORMAT
% DESC is a wrapper function for the gradient of the negative log
% likelihood of an Density Network model with respect to the latent postions
% and parameters.
% ARG params : vector of parameters and latent postions where the
% gradient is to be evaluated.
% ARG model : the model structure into which the latent positions
% and the parameters will be placed.
% RETURN g : the gradient of the negative log likelihood with
% respect to the latent positions and the parameters at the given
% point.
%
% SEEALSO : dnetLogLikeGradients, dnetExpandParam
%
% COPYRIGHT : Neil D. Lawrence, 2008

% MLTOOLS

model = dnetExpandParam(model, params);
g = - dnetLogLikeGradients(model);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


function [x_out] = applyBeamforming(arraySignal, dMic, fs, DOA, typeBF, flagGSC)

nChannels = size(arraySignal,1);
x = arraySignal;  % input signal, matrix  [nChannels * time]

%=========initialize parameters for spectral processing ===================%
len=floor(20 * fs / 1000); % Frame size in samples
if rem(len,2)==1,
    len=len+1;
end;
overlap = 0.50; % window overlap in percent of frame size
len1 = floor(len * overlap);
len2 = len - len1;

nFrames = floor(size(x, 2)/len2) - floor(len/len2);
xfinal = zeros(1, nFrames*len2+len1);

win = hamming(len)';  % define window
win2D = repmat(win, nChannels, 1);

nFFT = len;
f = linspace(0, fs, nFFT);% freqeuncy vector
%==========================================================================%

%============================Initialize covariance matirx==================%
noise_mean=zeros(nChannels, nFFT);
noiseConv = zeros(nChannels, nChannels, nFFT);
j=1;
for k=1:10   %use fisrt 10 segments to  estimate noise spectrum
    noise_spec = fft(win2D .* x(:, j:j+len-1), nFFT, 2);
    noise_mean = noise_mean + noise_spec;
    for freqidx = 1:nFFT
        noiseConv(:,:,freqidx) = noiseConv(:,:,freqidx) + (noise_spec(:, freqidx)' * noise_spec(:, freqidx));
    end
    j = j + len2;
end
noiseConv = noiseConv / 10;
%==========================================================================%

%================= generate beamformer ====================================%
hw = zeros(nChannels, nFFT);
if strcmp(typeBF, 'DSB')
    hw = dsb(nChannels, dMic, DOA, f);        % delay-and-sum filter
elseif  strcmp(typeBF, 'MVDR')
    hw = mvdrFF(nChannels, dMic, DOA, f, noiseConv);
end
%==========================================================================%

%================= speech enhancement per frame============================%
k = 1;
for n = 1 : nFrames

    xk = win2D .* x(:, k:k+len-1);
    spec = fft(xk, nFFT, 2);

    %vad or speech presence probability
    %power covariance matrix
    %updatebeamformer
    %rtf = rtf(, noiseConv)
    %hw = mvdr(rtf, noiseConv);

    if flagGSC == false
        spec_bf = sum(hw .* spec, 1);
    else
        [spec_bf, ak] = lmsGSC(X(:, 1:(L/2)+1), hw, 0, B, ak, mu);
    end

    xi_w = ifft( spec_bf, nFFT);
    xi_w = real(xi_w);

   % Overlap and add
    xfinal( k:k+len-1) =  xfinal(k:k+len-1) + win.* xi_w;

    k = k + len2;

end
%==========================================================================%

%===============================Output=====================================%
x_out = xfinal;
%==========================================================================%

end


text:
function [uOutput] = import_KJRelm(nFunction, sFilename)

% Filter function switchyard
if nFunction == FilterOp.getDescription
  uOutput = 'RELM - Kagan & Jackson PDE-based catalog';
elseif nFunction == FilterOp.importCatalog
  % Read formated data
  mData = textread(sFilename, '%s', 'delimiter', '\n', 'whitespace', '');
  % Create empty catalog
  uOutput = zeros(length(mData), 9);
  % Loop thru all lines of catalog and convert them
  for i = 1:length(mData)
    if rem(i,100) == 0 ; disp([ num2str(i) ' of ' num2str(length(mData)) ' events processed ']); end
    try
      uOutput(i,1) = str2num(mData{i}(33:40));  % Longitude
      uOutput(i,2) = str2num(mData{i}(25:31));  % Latitude
      uOutput(i,3) = str2num(mData{i}(7:10));   % Year
      uOutput(i,4) = str2num(mData{i}(12:13));  % Month
      uOutput(i,5) = str2num(mData{i}(15:16));  % Day
      uOutput(i,6) = str2num(mData{i}(59:62));  % Magnitude
      uOutput(i,7) = str2num(mData{i}(53:57));  % Depth
      uOutput(i,8) = str2num(mData{i}(19:20));  % Hour
      uOutput(i,9) = str2num(mData{i}(22:23));  % Minute
      %Create decimal year
      uOutput(i,3) = decyear([uOutput(i,3) uOutput(i,4) uOutput(i,5) uOutput(i,8) uOutput(i,9)]);
    catch
      msg.dbfprintf('Import: Problem in line %d of %s. Line ignored.\n',i, sFilename);
      uOutput(i,:)=nan;
    end
  end
  l = isnan(uOutput(:,1));
  uOutput(l,:) = [];
end

author: CelsoReyes


text:
%   y: dimensionanlity-reduced data
%	eigVector: eigen-vector obtained in kPCA
%   X: data matrix
%   para: parameter of Gaussian kernel
%	z: pre-image of y

%   Copyright by Quan Wang, 2011/05/10
%   Please cite: Quan Wang. Kernel Principal Component Analysis and its
%   Applications in Face Recognition and Active Shape Models.
%   arXiv:1207.3538 [cs.CV], 2012.

function z=kPCA_PreImage(y,eigVector,X,para)

iter=1000;
N=size(X,1);
d=max(size(y));

gamma=zeros(1,N);
for i=1:N
    gamma(i)=eigVector(i,1:d)*y;
end

z=mean(X)'; % initialization
fprintf('\nReconstruction: \n');
for count=1:iter
    fprintf('%d ', count);
    if mod(count,10)==0
        fprintf('\n');
    end

    pre_z=z;
    xx=bsxfun(@minus,X',z);
    xx=xx.^2;
    xx=-sum(xx)/(2*para.^2);
    xx=exp(xx);
    xx=xx.^gamma;

    z=xx*X/sum(xx);
    z=z';
    if norm(pre_z-z)/norm(z)<0.0001
        break;
    end
end
fprintf('\n');

author: Sable

repo: mcbench-benchmarks


text:
function [m]=qtt_mg_interp(d)
%The simplest multigrid operator in QTT as a TT-matrix
%   [M]=QTT_MG_INTERP(D) Computes the simplest MG-operator in the
%   QTT-format on 2^D grid. Its stencil is
% 1,	0
% 0.5,	0.5
% 0,	1
% 0,	0.5
%
%
% TT-Toolbox 2.2, 2009-2012
%
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------

prol0 = [1 0; 0.5 0.5; 0 1; 0 0.5];
h0 = [0 0; 0 0; 0 0; 0.5 0];
I = eye(2);
I12 = [0 1; 0 0];
I21 = [0 0; 1 0];

m = cell(d,1);

m{1}=zeros(4,2,2);
m{1}(:,:,1)=prol0; m{1}(:,:,2)=h0;
for i=2:d-1
    m{i}=zeros(2,2,2,2);
    m{i}(:,:,1,1)=I;
    m{i}(:,:,2,1)=I12; m{i}(:,:,2,2)=I21;
end;
m{d} = zeros(2,2,2);
m{d}(:,:,1)=I; m{d}(:,:,2)=I12;

m = tt_matrix(m);

end

author: oseledets

repo: TT-Toolbox


text:
function varargout = drawPolyhedron(nodes, faces, varargin)
%DRAWPOLYHEDRON Draw polyhedron defined by vertices and faces
%
%   drawPolyhedron(NODES, FACES)
%   Draws the polyhedron defined by vertices NODES and the faces FACES.
%   NODES is a NV-by-3 array containing coordinates of vertices, and FACES
%   is either a NF-by3 or NF-by-4 array containing indices of vertices of
%   the triangular or rectangular faces.
%   FACES can also be a cell array, in the content of each cell is an array
%   of indices to the nodes of the current face. Faces can have different
%   number of vertices.
%
%   H = drawPolyhedron(...);
%   Also returns a handle to the created patche.
%
%   Example:
%   [n f] = createSoccerBall;
%   drawPolyhedron(n, f);
%
%   See also:
%   polyhedra, drawMesh, drawPolygon
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 10/02/2005.
%

%   HISTORY
%   07/11/2005 update doc.
%   04/01/2007 typo
%   18/01/2007 add support for 2D polyhedra ("nodes" is N-by-2 array), and
%       make 'cnodes' a list of points instead of a list of indices
%   14/08/2007 add comment, add support for NaN in faces (complex polygons)
%   14/09/2007 rename as drawPolyhedron
%   16/10/2008 better support for colors
%   27/07/2010 copy to 'drawMesh'

%% Initialisations

% process input arguments
switch length(varargin)
    case 0
        % default color is red
        varargin = {'facecolor', [1 0 0]};
    case 1
        % use argument as color for faces
        varargin = {'facecolor', varargin{1}};
    otherwise
        % otherwise do nothing
end

% overwrites on current figure
hold on;

% if nodes are 2D points, add a z=0 coordinate
if size(nodes, 2) == 2
    nodes(1,3) = 0;
end

%% main loop : for each face

if iscell(faces)
    % array FACES is a cell array
    h = zeros(length(faces(:)), 1);

    for f = 1:length(faces(:))
        % get nodes of the cell
        face = faces{f};

        if sum(isnan(face))~=0
            % Special processing in case of multiple polygonal face.
            % each polygonal loop is separated by a NaN.

            % find indices of loops breaks
            inds = find(isnan(face));

            % replace NaNs by index of first vertex of each polygon
            face(inds(2:end))   = face(inds(1:end-1)+1);
            face(inds(1))       = face(1);
            face(length(face)+1)= face(inds(end)+1);
        end

        % draw current face
        cnodes  = nodes(face, :);
        h(f)    = patch(cnodes(:, 1), cnodes(:, 2), cnodes(:, 3), [1 0 0]);
    end

else
    % array FACES is a NC*NV indices array, with NV : number of vertices of
    % each face, and NC number of faces
    h = zeros(size(faces, 1), 1);
    for f = 1:size(faces, 1)
        % get nodes of the cell
        cnodes = nodes(faces(f,:)', :);
        h(f) = patch(cnodes(:, 1), cnodes(:, 2), cnodes(:, 3), [1 0 0]);
    end
end

% set up drawing options
if ~isempty(varargin)
    set(h, varargin{:});
end

% format output parameters
if nargout > 0
    varargout = {h};
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function [xpos, u] = firstzero(s)

%tstoolbox/@signal/firstzero
%   Syntax:
%     * [xpos, unit] = firstzero(s)
%
%   Give information about first zero of scalar signal s, using linear
%   interpolation.
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt
narginchk(1,1);

if ndim(s) ~= 1
	help(mfilename)
	return
end

d = data(s);
ze = find(d==0);
if isempty(ze)		% no exact zero found, looking for zero crossings
	tmp = d(2:end) .* d(1:end-1);
	ze = find(tmp < 0);
	if isempty(ze)
		error('no zero crossing found')
	else
		A = getaxis(s, 1);
		a = abs(d(ze(1)));
		b = abs(d(ze(1)+1));
		if (a+b) > 0
			xpos = first(A) + (ze(1)-1 + a/(a+b) )*delta(A);
		else
			xpos = first(A) + (ze(1)-1)*delta(A);
		end
		u = unit(A);
	% 	disp(['first zero/zero crossing at x = ' num2str(xpos) ' ' label(A)]);
	end
else		% use first sample that matches zero exactly
	A = getaxis(s, 1);
	xpos = first(A) + (ze(1)-1)*delta(A);
	u = unit(A);
end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
% exampleTextureCalc1.m
%
% Example script for texture calculation
%
% APA, 05/23/2016

global planC

%% EXAMPLE 1: Patch-wise texture
scanNum     = 1;
structNum   = 3;
descript    = 'CTV texture';
patchUnit   = 'vox'; % or 'cm'
patchSizeV  = [1 1 1];
category    = 1; % Haralick texture
dirctn      = 1; % 1: 3d neighbors , 2: 2d neighbors
numGrLevels = 16; % 32, 64, 256 etc..
energyFlg = 1; % or 0
entropyFlg = 1; % or 0
sumAvgFlg = 1; % or 0
homogFlg = 1; % or 0
contrastFlg = 1; % or 0
corrFlg = 1; % or 0
clustShadFlg = 1; % or 0
clustPromFlg = 1; % or 0
haralCorrFlg = 1; % or 0
flagsV = [energyFlg, entropyFlg, sumAvgFlg, corrFlg, homogFlg, ...
    contrastFlg, clustShadFlg, clustPromFlg, haralCorrFlg];
planC = createTextureMaps(scanNum,structNum,descript,...
    patchUnit,patchSizeV,category,dirctn,numGrLevels,flagsV,planC);

%% EXAMPLE 2: Texture for the entire structure
global planC
indexS = planC{end};
scanNum     = 1;
structNum   = 4;
numGrLevels = 16;
dirctn      = 1; % 2: 2d neighbors
cooccurType = 1; % 2: build separate cooccurrence for each direction

% Quantize the volume of interest
[rasterSegments, planC, isError]    = getRasterSegments(structNum,planC);
[mask3M, uniqueSlices]              = rasterToMask(rasterSegments, scanNum, planC);
scanArray3M                         = getScanArray(planC{indexS.scan}(scanNum));
SUVvals3M                           = mask3M.*double(scanArray3M(:,:,uniqueSlices));
[minr, maxr, minc, maxc, mins, maxs]= compute_boundingbox(mask3M);
maskBoundingBox3M                   = mask3M(minr:maxr,minc:maxc,mins:maxs);
volToEval                           = SUVvals3M(minr:maxr,minc:maxc,mins:maxs);
volToEval(maskBoundingBox3M==0)     = NaN;
quantizedM = imquantize_cerr(volToEval,numGrLevels);

% Buiild cooccurrence matrix
offsetsM = getOffsets(dirctn);
cooccurM = calcCooccur(quantizedM, offsetsM, numGrLevels, cooccurType);

% Reduce cooccurrence matrix to scalar features
glcmFlagS.energy = 1;
glcmFlagS.jointEntropy = 1;
glcmFlagS.jointMax = 1;
glcmFlagS.jointAvg = 1;
glcmFlagS.jointVar = 1;
glcmFlagS.contrast = 1;
glcmFlagS.invDiffMoment = 1;
glcmFlagS.sumAvg = 1;
glcmFlagS.corr = 1;
glcmFlagS.clustShade = 1;
glcmFlagS.clustProm = 1;
glcmFlagS.haralickCorr = 1;
glcmFlagS.invDiffMomNorm = 1;
glcmFlagS.invDiff = 1;
glcmFlagS.invDiffNorm = 1;
glcmFlagS.invVar = 1;
glcmFlagS.dissimilarity = 1;
glcmFlagS.diffEntropy = 1;
glcmFlagS.diffVar = 1;
glcmFlagS.diffAvg = 1;
glcmFlagS.sumVar = 1;
glcmFlagS.sumEntropy = 1;
glcmFlagS.clustTendency = 1;
glcmFlagS.autoCorr = 1;
glcmFlagS.invDiffMomNorm = 1;
glcmFlagS.firstInfCorr = 1;
glcmFlagS.secondInfCorr = 1;

featureS = cooccurToScalarFeatures(cooccurM, glcmFlagS);

%% Dominant orientation
global planC
indexS = planC{end};

scanNum     = 1;
structNum   = 6;
patchSizeV  = [3 3 3];

[rasterSegments, planC, isError]    = getRasterSegments(structNum,planC);
[mask3M, uniqueSlices]              = rasterToMask(rasterSegments, scanNum, planC);
scanArray3M                         = getScanArray(planC{indexS.scan}(scanNum));

SUVvals3M                           = mask3M.*double(scanArray3M(:,:,uniqueSlices));
[minr, maxr, minc, maxc, mins, maxs]= compute_boundingbox(mask3M);
maskBoundingBox3M                   = mask3M(minr:maxr,minc:maxc,mins:maxs);
volToEval                           = SUVvals3M(minr:maxr,minc:maxc,mins:maxs);
volToEval(maskBoundingBox3M==0)     = NaN;

% volToEval = scanArray3M; % for ITK comparison

position = [400 400 300 50];
waitFig = figure('name','Creating Texture Maps','numbertitle','off',...
            'MenuBar','none','ToolBar','none','position',position);
waitAx = axes('parent',waitFig,'position',[0.1 0.3 0.8 0.4],...
    'nextplot','add','XTick',[],'YTick',[],'yLim',[0 1],'xLim',[0 1]);
waitH = patch([0 0 0 0], [0 1 1 0], [0.1 0.9 0.1],...
    'parent', waitAx);

domOrient3M = calcDominantOrientation(volToEval, patchSizeV, waitH);

vol3M = zeros(size(volToEval));
vol3M(:,30:60,5) = 1;
vol3M(30:62,:,5) = 1;
dom3M = calcDominantOrientation(vol3M, patchSizeV, waitH);

close(waitFig)

%% Law's texture
global planC
indexS = planC{end};

scanNum     = 1;
structNum   = 1;
patchSizeV  = [3 3 3];

[rasterSegments, planC, isError]    = getRasterSegments(structNum,planC);
[mask3M, uniqueSlices]              = rasterToMask(rasterSegments, scanNum, planC);
scanArray3M                         = getScanArray(planC{indexS.scan}(scanNum));

SUVvals3M                           = mask3M.*double(scanArray3M(:,:,uniqueSlices));
[minr, maxr, minc, maxc, mins, maxs]= compute_boundingbox(mask3M);
maskBoundingBox3M                   = mask3M(minr:maxr,minc:maxc,mins:maxs);
volToEval                           = SUVvals3M(minr:maxr,minc:maxc,mins:maxs);
%volToEval(maskBoundingBox3M==0)     = NaN;
volToEval(maskBoundingBox3M==0)     = nanmean(volToEval(:));
meanVol = nanmean(volToEval(:));
if exist('padarray.m','file')
    paddedVolM = padarray(volToEval,[5 5 5],meanVol,'both');
else
    paddedVolM = padarray_oct(volToEval,[5 5 5],meanVol,'both');
end
lawsMasksS = getLawsMasks();

fieldNamesC = fieldnames(lawsMasksS);
numFeatures = length(fieldNamesC);
featuresM = zeros(sum(maskBoundingBox3M(:)),numFeatures);
for i = 1:numFeatures
    disp(i)
    text3M = convn(paddedVolM,lawsMasksS.(fieldNamesC{i}),'same');
    text3M = text3M(6:end-5,6:end-5,6:end-5);
    featuresM(:,i) = text3M(maskBoundingBox3M);
end

% Haralick textures

featuresM = zeros(sum(maskBoundingBox3M(:)),0);
for patchSiz = 1:3

%patchSizeV  = [2 2 2];
patchSizeV  = [patchSiz patchSiz 0];

%numGrLevels = 16;
for numGrLevels = [8 16 32]
offsetsM = getOffsets(2);

[rasterSegments, planC, isError]    = getRasterSegments(structNum,planC);
[mask3M, uniqueSlices]              = rasterToMask(rasterSegments, scanNum, planC);
scanArray3M                         = getScanArray(planC{indexS.scan}(scanNum));

SUVvals3M                           = mask3M.*double(scanArray3M(:,:,uniqueSlices));
[minr, maxr, minc, maxc, mins, maxs]= compute_boundingbox(mask3M);
maskBoundingBox3M                   = mask3M(minr:maxr,minc:maxc,mins:maxs);
volToEval                           = SUVvals3M(minr:maxr,minc:maxc,mins:maxs);
volToEval(maskBoundingBox3M==0)     = NaN;

% Haralick texture
energyFlg = 1; % or 0
entropyFlg = 1; % or 0
sumAvgFlg = 1; % or 0
homogFlg = 1; % or 0
contrastFlg = 1; % or 0
corrFlg = 1; % or 0
clustShadFlg = 1; % or 0
clustPromFlg = 1; % or 0
haralCorrFlg = 1; % or 0
flagsV = [energyFlg, entropyFlg, sumAvgFlg, corrFlg, homogFlg, ...
    contrastFlg, clustShadFlg, clustPromFlg, haralCorrFlg];
waitH = NaN;
[energy3M,entropy3M,sumAvg3M,corr3M,invDiffMom3M,contrast3M, ...
    clustShade3M,clustPromin3M,haralCorr3M] = textureByPatchCombineCooccur(volToEval,...
    numGrLevels,patchSizeV,offsetsM,flagsV,waitH);

featuresM(:,end+1) = energy3M(maskBoundingBox3M);
featuresM(:,end+1) = entropy3M(maskBoundingBox3M);
featuresM(:,end+1) = sumAvg3M(maskBoundingBox3M);
featuresM(:,end+1) = corr3M(maskBoundingBox3M);
featuresM(:,end+1) = invDiffMom3M(maskBoundingBox3M);
featuresM(:,end+1) = contrast3M(maskBoundingBox3M);
featuresM(:,end+1) = clustShade3M(maskBoundingBox3M);
featuresM(:,end+1) = clustPromin3M(maskBoundingBox3M);
featuresM(:,end+1) = haralCorr3M(maskBoundingBox3M);

end
end

[coeff,score,latVar] = pca(featuresM,'NumComponents',20);
figure, plot(cumsum(latVar)./sum(latVar)*100,'linewidth',2)
xlabel('Number of components','fontsize',20)
ylabel('Explained variance','fontsize',20)
set(gca,'fontsize',20)

figure,
for i = 1:4
    comp1M = zeros(size(maskBoundingBox3M));
    comp1M(maskBoundingBox3M) = score(:,3);
    comp1M = volToEval;
    subplot(2,2,i), imagesc(comp1M(:,:,i)), title(['slice: ',num2str(i)])
    axis equal, colormap('gray')
    axis off
end

comp1M = NaN*ones(size(maskBoundingBox3M));
comp1M(maskBoundingBox3M) = score(:,2);
figure, hist(comp1M(:),30)
title('Component 2','fontsize',20)

%% Neighborhood Gray Tone Difference Matrix (NGTDM)
global planC
indexS = planC{end};

% L:\Ziad\HUH_from_Maria\ASTRO_AAPM\Included_after_AAPM\mat_cropped\HUH2

scanNum     = 1;
structNum   = 27;
patchRadiusV  = [1 1 1];
numGrLevels = 16;

[rasterSegments, planC, isError]    = getRasterSegments(structNum,planC);
[mask3M, uniqueSlices]              = rasterToMask(rasterSegments, scanNum, planC);
scanArray3M                         = getScanArray(planC{indexS.scan}(scanNum));

SUVvals3M                           = mask3M.*double(scanArray3M(:,:,uniqueSlices));
[minr, maxr, minc, maxc, mins, maxs]= compute_boundingbox(mask3M);
maskBoundingBox3M                   = mask3M(minr:maxr,minc:maxc,mins:maxs);
volToEval                           = SUVvals3M(minr:maxr,minc:maxc,mins:maxs);
volToEval(maskBoundingBox3M==0)     = NaN;
quantizedM = imquantize_cerr(volToEval,numGrLevels);

hWait = NaN;
[s,p] = calcNGTDM(quantizedM,patchRadiusV,numGrLevels,hWait);

if exist('padarray.m','file')
    paddedM = padarray(quantizedM,[1 1 1],0,'both');
    paddedMaskM = padarray(maskBoundingBox3M,[1 1 1],0,'both');
else
    paddedM = padarray_oct(quantizedM,[1 1 1],0,'both');
    paddedMaskM = padarray_oct(maskBoundingBox3M,[1 1 1],0,'both');
end
[rV,cV,sV] = find3d(paddedMaskM);
mask_rcs = [rV(:),cV(:),sV(:)];
[NGTDM,vox_occurances_NGD26] = compute_3D_NGTDM_full_vol(paddedM,numGrLevels,mask_rcs);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function value = p34_param ( action, name, value )

%*****************************************************************************80
%
%% P34_PARAM defaults, sets, or gets the parameters for problem 34.
%
%  Discussion:
%
%    In typical use, the user calls
%      p34_param ( 'DEFAULT' )
%    Then, once the variables have been given default values, any
%    value can be retrieved, as for example by:
%      value = p34_param ( 'GET', 'G' )
%    or set to a new value, as in:
%      p34_param ( 'SET', 'G', 1 )
%
%    G is the gravitational force;
%    K is the damping coefficient;
%    L is the length of the string;
%    M is the mass of the pendulum bob.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 March 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string ACTION:
%    'DEFAULT' sets the variables to default values.
%    'SET' sets the variable chosen by NAME to have the value VALUE.
%    'GET' gets the variable chosen by NAME.
%
%    Input, string NAME, the name of the variable, which should
%    be 'G', 'K', 'L', 'M'.  NAME is ignored if the ACTION is 'DEFAULT'.
%
%    Input, real VALUE, the value to be assigned the variable NAME,
%    if the value of ACTION was 'SET'.  Otherwise VALUE is ignored.
%
%    Output, real VALUE, depends on the input value of ACTION:
%    'DEFAULT', returns the number of variables;
%    'SET', returns the input value of VALUE;
%    'GET', returns the value of the variable NAME.
%
  persistent g;
  persistent k;
  persistent l;
  persistent m;

  if ( s_eqi ( action, 'DEFAULT' ) )

    g = 32.0;
    k = 1.0;
    l = 1.0;
    m = 1.0;

    value = 4;

  elseif ( s_eqi ( action, 'SET' ) )

    if ( s_eqi ( name, 'G' ) )
      g = value;
    elseif ( s_eqi ( name, 'K' ) )
      k = value;
    elseif ( s_eqi ( name, 'L' ) )
      l = value;
    elseif ( s_eqi ( name, 'M' ) )
      m = value;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P34_PARAM - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized variable name: "%s".\n', name );
      error ( 'P34_PARAM - Fatal error!' );
    end

  elseif ( s_eqi ( action, 'GET' ) )

    if ( s_eqi ( name, 'G' ) )
      value = g;
    elseif ( s_eqi ( name, 'K' ) )
      value = k;
    elseif ( s_eqi ( name, 'L' ) )
      value = l;
    elseif ( s_eqi ( name, 'M' ) )
      value = m;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P34_PARAM - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized variable name: "%s".\n', name );
      error ( 'P34_PARAM - Fatal error!' );
    end

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'P34_PARAM - Fatal error!\n' );
    fprintf ( 1, '  Unrecognized action: "%s".\n', action );
    error ( 'P34_PARAM - Fatal error!' );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [MDP] = spm_MDP_check(MDP)
% MDP structure checking
% FORMAT [MDP] = spm_MDP_check(MDP)
%
% MDP.V(T - 1,P,F)      - P allowable policies of T moves over F factors
% or
% MDP.U(1,P,F)          - P allowable actions at each move
% MDP.T                 - number of outcomes
%
% MDP.A{G}(O,N1,...,NF) - likelihood of O outcomes given hidden states
% MDP.B{F}(NF,NF,MF)    - transitions among hidden under MF control states
% MDP.C{G}(O,T)         - prior preferences over O outcomes in modality G
% MDP.D{F}(NF,1)        - prior probabilities over initial states
%
% MDP.a{G}              - concentration parameters for A
% MDP.b{F}              - concentration parameters for B
% MDP.d{F}              - concentration parameters for D
%
% optional:
% MDP.s(F,T)            - vector of true states - for each hidden factor
% MDP.o(G,T)            - vector of outcome     - for each outcome modality
% MDP.u(F,T - 1)        - vector of action      - for each hidden factor
% MDP.w(1,T)            - vector of precisions
%
% if C or D are not specified, they will be set to default values (of no
% preferences and uniform priors over initial steps).  If there are no
% policies, it will be assumed that I = 1 and all policies (for each
% marginal hidden state) are allowed.
%__________________________________________________________________________
% Copyright (C) 2005 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_MDP_check.m 7766 2020-01-05 21:37:39Z karl $

% deal with a sequence of trials
%==========================================================================

% if there are multiple structures check each separately
%--------------------------------------------------------------------------
if numel(MDP) > 1
    for m = 1:size(MDP,1)
        for i = 1:size(MDP,2)
            mdp(m,i) = spm_MDP_check(MDP(m,i));
        end
    end
    MDP   = mdp;
    return
end

% fill in (posterior or  process) likelihood and priors
%--------------------------------------------------------------------------
if ~isfield(MDP,'A'), MDP.A = MDP.a; end
if ~isfield(MDP,'B'), MDP.B = MDP.b; end

% check format of likelihood and priors
%--------------------------------------------------------------------------
if ~iscell(MDP.A), MDP.A = {full(MDP.A)}; end
if ~iscell(MDP.B), MDP.B = {full(MDP.B)}; end

if isfield(MDP,'a'), if ~iscell(MDP.a), MDP.a = {full(MDP.a)}; end; end
if isfield(MDP,'b'), if ~iscell(MDP.b), MDP.b = {full(MDP.b)}; end; end

% check dimensions and orders
%==========================================================================

% numbers of transitions, policies and states
%--------------------------------------------------------------------------
Nf  = numel(MDP.B);                 % number of hidden state factors
Ng  = numel(MDP.A);                 % number of outcome factors
for f = 1:Nf
    Nu(f)    = size(MDP.B{f},3);    % number of hidden controls
    Ns(f)    = size(MDP.B{f},1);    % number of hidden states
    MDP.B{f} = double(MDP.B{f});
end
for g = 1:Ng
    No(g)    = size(MDP.A{g},1);    % number of outcomes
    MDP.A{g} = double(MDP.A{g});
end

% check policy specification (create default moving policy U, if necessary)
% V = V(Nt,Np,Nf)
% U = U(Np,Nf)
%--------------------------------------------------------------------------
if isfield(MDP,'U')
    if size(MDP.U,1) == 1 && size(MDP.U,3) == Nf
        MDP.U = shiftdim(MDP.U,1);
    end
end
try
    V(1,:,:) = MDP.U;               % allowable actions (1,Np)
catch
    try
        V = MDP.V;                  % allowable policies (T - 1,Np)
    catch

        % allowable (moving) policies using all allowable actions
        %------------------------------------------------------------------
        for f = 1:Nf
            u = 1;
            for i = 1:Nf
                if i == f
                    u = kron(1:Nu(i),u);
                else
                    u = kron(ones(1,Nu(i)),u);
                end
            end
            MDP.U(:,f)  = u;
        end
        V(1,:,:) = MDP.U;
    end
end
MDP.V = V;

% check policy specification
%--------------------------------------------------------------------------
if Nf ~= size(V,3) && size(V,3) > 1
    error('please ensure V(:,:,1:Nf) is consistent with MDP.B{1:Nf}')
end

% check preferences
%--------------------------------------------------------------------------
if ~isfield(MDP,'C')
    for g = 1:Ng
        MDP.C{g} = zeros(No(g),1);
    end
end
for g = 1:Ng
    if iscell(MDP.C)
        if isvector(MDP.C{g})
            MDP.C{g} = spm_vec(MDP.C{g});
        end
        if No(g) ~= size(MDP.C{g},1)
            error(['please ensure A{' num2str(g) '} and C{' num2str(g) '} are consistent'])
        end
    end
end

% check iinitial states
%--------------------------------------------------------------------------
if ~isfield(MDP,'D')
    for f = 1:Nf
        MDP.D{f} = ones(Ns(f),1);
    end
end
if Nf  ~= numel(MDP.D)
    error('please ensure V(:,:,1:Nf) is consistent with MDP.D{1:Nf}')
end

% check iinitial states and internal consistency
%--------------------------------------------------------------------------
if Nf  ~= numel(MDP.D)
    error('please ensure V(:,:,1:Nf) is consistent with MDP.D{1:Nf}')
end
for f = 1:Nf
    if Ns(f) ~= size(MDP.D{f},1)
        error(['please ensure B{' num2str(f) '} and D{' num2str(f) '} are consistent'])
    end
    if size(V,3) > 1
        if Nu(f) < max(spm_vec(V(:,:,f)))
            error(['please check V(:,:,' num2str(f) ') or U(:,:,' num2str(f) ')'])
        end
    end
    for g = 1:Ng
        Na  = size(MDP.A{g});
        if ~all(Na(2:end) == Ns)
            error(['please ensure A{' num2str(g) '} and D{' num2str(f) '} are consistent'])
        end
    end
end

% check probability matrices are properly specified
%--------------------------------------------------------------------------
for f = 1:Nf
    if ~all(spm_vec(any(MDP.B{f},1)))
        error(['please check B{' num2str(f) '} for missing entries'])
    end
end
for g = 1:Ng
    if ~all(spm_vec(any(MDP.A{g},1)))
        error(['please check A{' num2str(g) '} for missing entries'])
    end
end

% check initial states
%--------------------------------------------------------------------------
if isfield(MDP,'s')
    if size(MDP.s,1) > Nf
        error('please specify an initial state MDP.s for %i factors',Nf)
    end
    f  = max(MDP.s,[],2)';
    if any(f > Ns(1:numel(f)))
        error('please ensure initial states MDP.s are consistent with MDP.B')
    end
end

% check outcomes if specified
%--------------------------------------------------------------------------
if isfield(MDP,'o')
    if numel(MDP.o)
        if size(MDP.o,1) ~= Ng
            error('please specify an outcomes MDP.o for %i modalities',Ng)
        end
        if any(max(MDP.o,[],2) > No(:))
            error('please ensure outcomes MDP.o are consistent with MDP.A')
        end
    end
end

% check (primary link array if necessary)
%--------------------------------------------------------------------------
if isfield(MDP,'link')

    % cardinality of subordinate level
    %----------------------------------------------------------------------
    nf    = numel(MDP.MDP(1).B);               % number of hidden factors
    for f = 1:nf
        ns(f)    = size(MDP.MDP(1).B{f},1);    % number of hidden states
    end

    % check the size of link
    %----------------------------------------------------------------------
    if ~all(size(MDP.link) == [nf,Ng]);
        error('please check the size of link {%i,%i}',nf,Ng)
    end

    % convert matrix to cell array if necessary
    %----------------------------------------------------------------------
    if isnumeric(MDP.link)
        link  = cell(nf,Ng);
        for f = 1:size(MDP.link,1)
            for g = 1:size(MDP.link,2)
                if MDP.link(f,g)
                    link{f,g} = spm_speye(ns(f),No(g),0);
                end
            end
        end
        MDP.link = link;
    end

    % check sizes of cell array
    %----------------------------------------------------------------------
    for f = 1:size(MDP.link,1)
        for g = 1:size(MDP.link,2)
            if ~isempty(MDP.link{f,g})
                if ~all(size(MDP.link{f,g}) == [ns(f),No(g)]);
                    error('please check link{%i,%i}',f,g)
                end
            end
        end
    end

end

% Empirical prior preferences
%--------------------------------------------------------------------------
if isfield(MDP,'linkC')
    if isnumeric(MDP.linkC)
        linkC  = cell(numel(MDP.MDP.C),Ng);
        for f = 1:size(MDP.linkC,1)
            for g = 1:size(MDP.linkC,2)
                if MDP.linkC(f,g)
                    linkC{f,g} = spm_speye(size(MDP.MDP.C{f},1),No(g),0);
                end
            end
        end
        MDP.linkC = linkC;
    end
end

% Empirical priors over policies
%--------------------------------------------------------------------------
if isfield(MDP,'linkE')
    if isnumeric(MDP.linkE)
        linkE  = cell(1,Ng);
        for g = 1:size(MDP.linkE,2)
            if MDP.linkE(g)
                linkE{g} = spm_speye(size(MDP.MDP.E,1),No(g),0);
            end
        end
        MDP.linkE = linkE;
    end
end

% check factors and outcome modalities have proper labels
%--------------------------------------------------------------------------
for i = 1:Nf

    % name of factors
    %----------------------------------------------------------------------
    try
        MDP.label.factor(i);
    catch
        try
            MDP.label.factor{i} = MDP.Bname{i};
        catch
            MDP.label.factor{i} = sprintf('factor %i',i);
        end
    end

    % name of levels of each factor
    %----------------------------------------------------------------------
    for j = 1:Ns(i)
        try
            MDP.label.name{i}(j);
        catch
            try
                MDP.label.name{i}{j} = MDP.Sname{i}{j};
            catch
                MDP.label.name{i}{j} = sprintf('state %i(%i)',j,i);
            end
        end
    end

    % name of actions under each factor
    %----------------------------------------------------------------------
    for j = 1:Nu(i)
        try
            MDP.label.action{i}(j);
        catch
            MDP.label.action{i}{j} = sprintf('act %i(%i)',j,i);
        end
    end
end

% name of outcomes under each modality
%--------------------------------------------------------------------------
for i = 1:Ng
    try
        MDP.label.modality(i);
    catch
        try
            MDP.label.modality{i} = MDP.Bname{i};
        catch
            MDP.label.modality{i} = sprintf('modality %i',i);
        end
    end
    for j = 1:No(i)
        try
            MDP.label.outcome{i}(j);
        catch
            try
                MDP.label.outcome{i}{j} = MDP.Oname{i}{j};
            catch
                MDP.label.outcome{i}{j} = sprintf('outcome %i(%i)',j,i);
            end
        end
    end
end

% check names are specified properly
%--------------------------------------------------------------------------
if isfield(MDP,'Aname')
    if numel(MDP.Aname) ~= Ng
        error('please specify an MDP.Aname for each modality')
    end
else
    MDP.Aname = MDP.label.modality;
end
if isfield(MDP,'Bname')
    if numel(MDP.Bname) ~= Nf
        error('please specify an MDP.Bname for each factor')
    end
else
    MDP.Bname = MDP.label.factor;
end

author: spm


text:
function check_roi_extraction()

% test ROI extraction
mask_image = fmri_data(which('atlas_labels_combined.img'));
dat = mask_image;
wh_region = dat.dat;

% create timeseries data -- add data "after" 1st image
nimgs = 50;

n = unique(wh_region);
ts = linspace(-10, 10, nimgs);
v = size(dat.dat, 1);
dat.dat = zeros(v, nimgs);

for i = 1:length(n)
    my_indx = wh_region == n(i);

    dat.dat(my_indx, :) = repmat(ts .* sqrt(i), sum(my_indx), 1);
end

fprintf('\nDATA GENERATED\n\n');

cl1 = extract_roi_averages(dat, mask_image, 'unique_mask_values');

% extracted average values for all the regions in mask_image.
all_reg1 = cat(2, cl1(:).dat);

dat.fullpath = fullfile(pwd, 'test_image.img');
write(dat);
fprintf('\nSAVED DATA TO %s\n\n', dat.fullpath);

% reload
fprintf('\nLOADING DATA FROM %s\n\n', dat.fullpath);
dat = fmri_data(dat.fullpath, which('gray_matter_mask.img'));

% region obj -> extract_data
% test extract_roi_averages
cl2 = extract_roi_averages(dat, mask_image, 'unique_mask_values');

% extracted average values for all the regions in mask_image.
all_reg2 = cat(2, cl2(:).dat);

% uniq_reg = sort(unique(mask_image.dat));

% compare roi averages of initially data generated and data loaded from file
if is_equal(all_reg1, all_reg2, 0.001)
    fprintf('\nExtracted ROI averages are equal with an error threshold of 0.001\nPASS\n');
else
    fprintf('\nExtracted ROI averages are not equal with an error threshold of 0.001\nFAIL\n');
end

% ------------------------------------------------------
% INLINE FUNCTION
% ------------------------------------------------------

    function res = is_equal(A, B, error_threshold)
        err = abs((A-B)./B) <= error_threshold;

        res = all(err(:));
    end % is_equal

end  % main function

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function p=classprob(t,j)
%CLASSPROB Class probabilities for tree nodes.
%   P=CLASSPROB(T) returns an N-by-M array P of class probabilities for the
%   nodes in the classification tree T, where N is the number of nodes
%   and M is the number of classes.  For any node number K, the class
%   probabilities P(K,:) are the estimated probabilities for each class
%   for a point satisfying the conditions for node K.
%
%   P=CLASSPROB(T,J) takes an array J of node numbers and returns the
%   class probabilities for the specified nodes.
%
%   See also CLASSREGTREE, CLASSREGTREE/NUMNODES.

%   Copyright 2006-2007 The MathWorks, Inc.
%   $Revision: 1.1.6.3 $  $Date: 2007/02/15 21:48:02 $

if isequal(t.method,'regression')
    error('stats:classregtree:classprob:NotClassification',...
          'The CLASSPROB method is not available for regression trees.');
end
if nargin>=2 && ~validatenodes(t,j)
    error('stats:classregtree:classprob:InvalidNode',...
          'J must be an array of node numbers or a logical array of the proper size.');
end

if nargin<2
    p = t.classprob;
else
    p = t.classprob(j,:);
end

author: Cloud-CV

repo: object-proposals


text:
classdef TTeMPS_op_laplace
% TTeMPS_op_laplace
%
%   A MATLAB class for representing and manipulating
%   Laplace-like operators in the TT/MPS operator format,
%
%   Laplace-like operators are of the form
%         L \otimes I \otimes I \otimes I
%       + I \otimes L \otimes I \otimes I
%       + ...
%       + I \otimes I \otimes I \otimes L
%
%   with e.g. L being the discrete 1D-Laplacian matrix.

%   TTeMPS Toolbox.
%   Michael Steinlechner, 2013-2016
%   Questions and contact: michael.steinlechner@epfl.ch
%   BSD 2-clause license, see LICENSE.txt

properties( SetAccess = public, GetAccess = public )

    L0
    U           % core tensors as 4D doubles
    rank
    order
    size_row
    size_col
    V_L
    E_L

end

% Set methods for Cores
methods

    function A = set.U( A, U_);

        A.U = U_;
        A = update_properties( A );

    end

    function A = update_properties( A );

        A.rank = [1,  2*ones(1, length(A.U)-1), 1];  % the TT rank is always two for such Laplace-like tensors
        size_col_ = cellfun( @(y) size(y,1), A.U);
        A.size_col = size_col_ ./ (A.rank(1:end-1).*A.rank(2:end));
        A.size_row = cellfun( @(y) size(y,2), A.U);
        A.order = length( A.size_row );

    end

end

methods( Access = public )

    function A = TTeMPS_op_laplace(varargin)
    %LAPLACE Construct a tensor in TT/MPS format and return TTeMPS_op_laplace object.
    %
    %   A = TTEMPS_OP_LAPLACE( L, D ) creates the D-dimensional Laplace-like
    %       TT/MPS operator using the supplied matrix L.
    %
    %
        if nargin == 1
            % debug constructor
            A.U = varargin{1};
            A = update_properties( A );
            A.L0 = [];
            A.V_L = [];
            A.E_L = [];
        % Default constructor
        elseif (nargin == 2)

            % only one matrix passed
            L = varargin{1};
            A.L0 = L;
            d = varargin{2};

            A.V_L = [];
            A.E_L = [];

            [m,n] = size( L );
            E = speye( m, n );
            a_1 = sparse( 1, 1, 1, 2, 1 );
            a_mid = sparse( 2, 1, 1, 4, 1 );
            a_end = sparse( 2, 1, 1, 2, 1 );
            b_1 = sparse( 2, 1, 1, 2, 1 );
            b_mid = sparse( [1;4], [1;1], [1;1], 4, 1 );
            b_end = sparse( 1, 1, 1, 2, 1 );

            A.U = cell( 1, d );
            A.U{1} = kron( L, a_1 ) + kron( E, b_1 );
            A_mid = kron( L, a_mid ) + kron( E, b_mid );
            for i=2:d-1
                A.U{i} = A_mid;
            end
            A.U{d} = kron( L, a_end ) + kron( E, b_end );

            A = update_properties( A );

        else
            error('Invalid number of arguments.')
        end
    end

    % Other public functions
    y = apply( A, x, idx );
    A = mtimes( B, A );
    res = contract( A, x, y, idx );

    disp( A, name );
    display( A );

    B = TTeMPS_op_laplace_to_TTeMPS_op( A );
    B = TTeMPS_op_laplace_to_TT_matrix( A );
    expB = constr_precond_inner( A, X, mu );

    function A = initialize_precond( A )
        [A.V_L, A.E_L] = eig(full(A.L0));
        A.E_L = diag( A.E_L );
    end
end

end

author: NicolasBoumal

repo: manopt


text:
function obj_write ( output_file_name, node_num, face_num, normal_num, ...
  order_max, node_xyz, face_order, face_node, normal_vector, vertex_normal )

%*****************************************************************************80
%
%% OBJ_WRITE writes graphics information to an Alias OBJ file.
%
%  Discussion:
%
%    If no normal vectors are supplied (NORMAL_NUM <= 0) then
%    a simple format is used for the "F" records.  Otherwise,
%    the "v//vn" format is used.
%
%  Example:
%
%    #  no_normals.obj
%
%    g Group002
%
%    v -3.269770 -39.572201 0.876128
%    v -3.263720 -39.507999 2.160890
%    ...
%    v 0.000000 -9.988540 0.000000
%
%    f 8 9 11 10
%    f 12 13 15 14
%    ...
%    f 788 806 774
%
%    #  normals_supplied.obj
%
%    g Group001
%
%    v -3.269770 -39.572201 0.876128
%    v -3.263720 -39.507999 2.160890
%    ...
%    v 0.000000 -9.988540 0.000000
%
%    vn 0.0 1.0 0.0
%    vn 1.0 0.0 0.0
%    ...
%    vn 0.0 0.0 1.0
%
%    f 8//1 9//2 11//3 10//4
%    f 12//5 13//6 15//7 14//8
%    ...
%    f 788//800 806//803 774//807
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILE_NAME, the name of the output file.
%
%    Input, integer NODE_NUM, the number of points.
%
%    Input, integer FACE_NUM, the number of faces.
%
%    Input, integer NORMAL_NUM, the number of normal vectors.
%
%    Input, integer ORDER_MAX, the maximum number of vertices per face.
%
%    Input, real NODE_XYZ(3,NODE_NUM), the coordinates of points.
%
%    Input, integer FACE_ORDER(FACE_NUM), the number of vertices per face.
%
%    Input, integer FACE_NODE(ORDER_MAX,FACE_NUM), the nodes making faces.
%
%    Input, real NORMAL_VECTOR(3,NORMAL_NUM), normal vectors.
%
%    Input, integer VERTEX_NORMAL(ORDER_MAX,FACE_NUM), the indices of normal
%    vectors per vertex.
%
  output_file_unit = fopen ( output_file_name, 'wt' );

  if ( output_file_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'OBJ_WRITE - Fatal error!\n' );
    fprintf ( 1, '  Could not open the output file "%s".\n', ...
      output_file_name );
    return
  end

  text_num = 0;

  fprintf ( output_file_unit, '# %s\n', output_file_name );
  fprintf ( output_file_unit, '#  created by OBJ_WRITE.\n' );
  fprintf ( output_file_unit, '\n' );
  fprintf ( output_file_unit, 'g Group001\n' );

  text_num = text_num + 4;
%
%  V: vertex coordinates.
%  For some reason, a fourth "coordinate" may be recommended.
%  What is its meaning?
%
  if ( 0 < node_num )
    fprintf ( output_file_unit, '\n' );
    text_num = text_num + 1;
  end

  w = 1.0;
  for node = 1 : node_num
    fprintf ( output_file_unit, 'v  %f  %f  %f  %f\n',  node_xyz(1:3,node), w );
    text_num = text_num + 1;
  end
%
%  VN: normal vectors.
%
  if ( 0 < normal_num )

    fprintf ( output_file_unit, '\n' );
    text_num = text_num + 1;

    for normal = 1 : normal_num
      fprintf ( output_file_unit, 'vn  %f  %f  %f\n', normal_vector(1:3,normal) );
      text_num = text_num + 1;
    end

  end
%
%  F: Faces, specified as a list of triples, one triple for each vertex:
%    vertex index/vertex texture index/vertex normal index
%
  if ( 0 < face_num )
    fprintf ( output_file_unit, '\n' );
    text_num = text_num + 1;
  end

  for face = 1 : face_num

    text = 'f';

    if ( normal_num <= 0 )

      for vertex = 1 : face_order(face)
        fprintf ( output_file_unit, '  %d', face_node(vertex,face) );
      end
      fprintf ( output_file_unit, '\n' );
      text_num = text_num + 1;

    else

      for vertex = 1 : face_order(face)
        fprintf ( output_file_unit, '  %d//%d', ...
          face_node(vertex,face), vertex_normal(vertex,face) );
      end
      fprintf ( output_file_unit, '\n' );
      text_num = text_num + 1;

    end

  end

  fclose ( output_file_unit );
%
%  Report.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'OBJ_WRITE:\n' );
  fprintf ( 1, '  Wrote %d text lines to "%s".\n', ...
    text_num, output_file_name );

  return
end


text:
%% Homework 11
% We wish to solve the following PDE using the approximate factorization
% method:
%
% $\frac{\partial \phi }{\partial t}+U(x,y)\frac{\partial \phi }{\partial
% x}+V(x,y)\frac{\partial \phi }{\partial y}=\alpha \left (
% \frac{\partial^2 \alpha }{\partial x^2}-\frac{\partial^2 \phi }{\partial
% y^2} \right)$
%
% This is 2D convection-diffusion equation

close all
clear all

%% Grid

ly = [-1,1];                lx =[0,10];
dy = 0.2;                   dx = 0.5;
y  = ly(1):dy:ly(2);        x  = lx(1):dx:lx(2); % Physical grid
q  = length(y);             p  = length(x);
phi = zeros(q,p); % Numerical grid

t_end = 1;      dt = 0.1;   n = (t_end-0)/dt;

%% Boundary Conditions

phi(:,1) = 1;
phi(1,:) = 1;

%% Constants & coheficients
% assume: $\alpha = const$.

alpha = 1;

% 1. assume: $u(x,y)=x*y% and $v(x,y)=x/(y+0.1)$,
% 2. assume: $u(x,y)=1$ and $v(x,y)=1$.

u = zeros(q,p);
for j = 1:q
    for i = 1:p
        u(j,i) = x(i)*y(j);
        % u(j,i) = 1
    end
end

v = zeros(q,p);
for j = 1:q
    for i = 1:p
        v(j,i) = (1/(y(j)+0.1))*x(i);
        % v(j,i) = 1
    end
end

%% Formulate Base matrices
% Dy
% Constants
a = dt/2*(alpha/dy^2 - v/(2*dy));
b = (1 - dt/dy^2);
c = dt/2*(alpha/dy^2 + v/(2*dy));

Qy = zeros(q-2,q-2,p);
for i = 2:p-1
    e = zeros(q-2,3);
    e(:,3) =  a(1:q-2,i);
    e(:,2) =  b*ones(q-2,1); % But the last No. in de Diag must be modif.:
    e(q-2,2) =  a(q-1,i)-b(q-1,i);
    e(:,1) =  c(3:q,i);
    Qy(:,:,i) = spdiags(e,-1:1,q-2,q-2);
end

% Dx
% Constants
% f = dt/2*(alpha/dx^2 - u/(2*dx));
% g = (1 - dt/dx^2);
% h = dt/2*(alpha/dx^2 - u/(2*dx));
%
% e1 = a(1,3:p-1);
% e2 = b*ones(1,p-2);
% e3 = c(1,2:p-2);
% Qx = spdiags([e3 e2 e1],[-1 0 1],p-2,p-2);

%% Compute $f$ for time n step

%% Boundary Conditions

repo: Aero-matlab


text:
function wNew = response_apply(w,filterObj,sourceType,source)

%RESPONSE_APPLY Deconvolve an instrument response from a waveform object.
% W = RESPONSE_APPLY(W,FILTEROBJ,SOURCETYPE,SOURCE) applies an instrument
% response to waveform W, filtered on the filter range specified by
% filterobject FILTEROBJ. The response is read from SOURCE which is of type
% SOURCETYPE. SOURCETYPE can be:
%   polezero:  A polezero structure. See HELP RESPONSE_GET_FROM_POLEZERO
%   antelope:  An Antelope database. See HELP RESPONSE_GET_FROM_DB
%  structure:  A response structure. See HELP RESPONSE_GET_FROM_DB
%
%  Example uses
%     w = response_apply(w,filterObj,'polezero',polezero)
%     w = response_apply(w,filterObj,'antelope',dbName)
%     w = response_apply(w,filterObj,'structure',response)
%
%
% see also response_get_from_db, response_get_from_polezero

% Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
%         (borrowing liberally from codes written by M. Haney)
% $Date$
% $Revision$

% CHECK INPUTS
if ~ischar(sourceType) || size(sourceType,1)~=1
    error('SOURCETYPE must be the same for all waveforms');
end
if strncmpi(sourceType,'antelope',3) && size(sourceType,1)~=1
    error('All responses must be sourced from a single database');
end

% PREPARE INDIVIDUAL WAVEFORMS

filterObj = repmat(filterObj,size(w));
source = repmat({source},size(w));
for n = 1:numel(w)
    wNew(n) = response_apply_one(w(n),filterObj(n),sourceType,source(n));
end

wNew = reshape(wNew,size(w));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PROCESS INDIVIDUAL WAVEFORMS

function w = response_apply_one(w,filterObj,sourceType,source)

if numel(w)~=1
    error('Only one waveform should be passed to this subroutine.');
end

% % CHECK THE SOURCE OF THE RESPONSE
% if isa(source,'char')                 % antelope db
%     disp(['Attempting to retreive response information from database: ' source ' ...']);
% elseif  isa(source,'struct')          % polezero structure
%     disp('Using explicit poles/zeros ...');
% else
%     error('response source argument not recognized');
% end

% PREPARE VARIABLE SPACE
filterBand = get(filterObj,'CUTOFF');
period = get(w,'PERIOD');
nyquist = get(w,'NYQ');
rawData = double(w);
dataLength = numel(rawData);
rawData = reshape(rawData,1,dataLength);

% PREPARE TRACE DATA
% Create a zero padded Tukey taper
%    Taper size is determined by the high pass filter frequency
%    Tapered waveform is zeros padded on either end to triple the trace length
% Glenn 20160519: could the follow commands be simplified using tukeywin?
taperFullWidth = round(0.5/(filterBand(1)*period))*2; % guarantees an even number
taperAmp = hanning(taperFullWidth)';
taperAmp = [taperAmp(1:(taperFullWidth/2)) ones(1,(dataLength-taperFullWidth)) taperAmp(((taperFullWidth/2)+1):taperFullWidth)];
taperAmp = [ zeros(1,dataLength) taperAmp zeros(1,dataLength) ];
rawData = [ zeros(1,dataLength) rawData zeros(1,dataLength) ].*taperAmp;
dataLength = numel(rawData);

% CREATE FULL INVERSE RESPONSE VECTORS
% Seperate routines for even and odd trace length
% Full responses are reconstituted from positive only side
% ws and respInv are positive frequencies only
% wsFull and respInvFull have negative frequencies as well
if floor(dataLength/2) == ceil(dataLength/2)
    halfDataLength=dataLength/2;
    ws = 2*pi*[0:(halfDataLength)]*(1/(halfDataLength))*nyquist;
    [respInv,response] = prep_inverse_response(w,ws,sourceType,source);
    wsFull = [ -ws(halfDataLength+1:-1:2) ws(1:halfDataLength) ];
    respInvFull = ([ real(respInv(end:-1:2)) real(respInv(1:end-1)) ] + 1i*[ -imag(respInv(end:-1:2)) imag(respInv(1:end-1)) ]);
else
    halfDataLength = (dataLength+1)/2;
    ws = 2*pi * [0:(halfDataLength-1)] * (2/dataLength) * nyquist;
    [respInv,response] = prep_inverse_response(w,ws,sourceType,source);
    wsFull = [ -ws(end:-1:2) ws(1:end) ];
    respInvFull = ([ real(respInv(end:-1:2)) real(respInv(1:end)) ] + 1i*[ -imag(respInv(end:-1:2)) imag(respInv(1:end)) ]);
end

% APPLY THE RESPONSE
newData = real(ifft(ifftshift(fftshift(fft(rawData)).*(respInvFull))));

% BANDPASS FILTER THE DATA
[z q] = butter(get(filterObj,'POLES'),[(filterBand(1)/nyquist) (filterBand(2)/nyquist)]);
% Glenn 20160519 could the follow 3 lines be replaced using filtfilt?
newData = filter(z,q,newData);
newData = filter(z,q,fliplr(newData));
newData = fliplr(newData);

% REMOVE TRACE PADDING
newData = newData( (dataLength/3)+1 : 2*(dataLength/3) );

% RETURN INSTURMENT CORRECTED WAVEFORM
w = set(w,'DATA',newData);
w = addhistory(w,'A frequency response was deconvolved using RESPONSE_APPLY. See RESPONSE field.');
w = addhistory(w,'Butterworth filter was doubly applied (FILTFILT). See FILTER field.');
w = addfield(w,'RESPONSE',response);
w = addfield(w,'FILTER',filterObj);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SET UP THE INVERSE RESPONSE FUNCTIONS
% could facilitate calls to a different response loader

function [respInv,response] = prep_inverse_response(w,ws,sourceType,source)

switch lower(sourceType)

    case {'antelope'}
        sta = get(w,'STATION');
        chan = get(w,'CHANNEL');
        time = get(w,'START_MATLAB');
        response = response_get_from_db(sta,chan,time,(ws/(2*pi)),source);
        respInv = 1./response.values;
        respInv(1) = 0;
        respInv = reshape(respInv,1,numel(respInv));

    case {'polezero'}
        response = response_get_from_polezero(ws,source);
        respInv = 1./response.values;
        respInv(1) = 0;
        respInv = reshape(respInv,1,numel(respInv));

    case {'structure'}
        disp('Structure is not yet implemented')
    otherwise
        error('response source argument not recognized');

end

author: geoscience-community-codes

repo: GISMO


text:
function conv2out_mem()
    global config mem;
    idx = reshape(1:config.output_size(1) * config.output_size(2), config.output_size(1), config.output_size(2));
    idx = im2col(idx, config.kernel_size(size(config.kernel_size, 1), :));
    %mem.gen_out_matrix = config.NEW_MEM(idx);
    mem.output = config.NEW_MEM(zeros(config.output_size(1), config.output_size(2), config.output_size(3), config.batch_size));

    counts = ones(size(idx));
    counts = reshape(accumarray(idx(:), counts(:)), config.output_size(1), config.output_size(2));
    mem.one_over_add_counts = config.NEW_MEM(1 ./ counts);
    mem.gen_out_matrix = config.NEW_MEM(zeros(size(idx, 1)*config.output_size(3), size(idx, 2)));
    for m = 1:config.output_size(3)%config.chs
        mem.gen_out_matrix((m-1)*size(idx, 1)+1:m*size(idx, 1), :) = idx + ((m-1)*max(max(idx)));
    end

    mem.gen_out_matrix = mem.gen_out_matrix(:);
    if config.batch_size > 1
        h = size(mem.gen_out_matrix, 1);
        mem.gen_out_matrix = repmat(mem.gen_out_matrix, [config.batch_size 2]);
        for m = 1:config.batch_size
            mem.gen_out_matrix((m-1)*h+1:m*h, 2) = m;
        end
    end
end

author: jimmy-ren


text:
function [lat, lon, gam, k] = utm_inv(zone, northp, x, y)
%UTM_INV  Forward transverse Mercator projection
%
%   [LAT, LON] = UTM_INV(ZONE, NORTHP, X, Y)
%   [LAT, LON, GAM, K] = UTM_INV(ZONE, NORTHP, X, Y)
%
%   performs the inverse universal transverse Mercator projection of points
%   (X,Y) to (LAT,LON) using ZONE and NORTHP.  X and Y can be scalars or
%   arrays of equal size.  ZONE should be an integer in [1,60] and NORTHP
%   is a logical indicating whether the transformation should use the false
%   northing for the northern (NORTHP = true) or southern (NORTHP = false)
%   hemisphere.  The forward projection is given by UTM_FWD.
%
%   GAM and K give metric properties of the projection at (LAT,LON); GAM is
%   the meridian convergence at the point and K is the scale.
%
%   LAT, LON, GAM are in degrees.  The projected coordinates X, Y are in
%   meters.  K is dimensionless.
%
%   This implementation for the UTM projection is based on the series
%   method described in
%
%     C. F. F. Karney, Transverse Mercator with an accuracy of a few
%     nanometers, J. Geodesy 85(8), 475-485 (Aug. 2011);
%     Addenda: http://geographiclib.sf.net/tm-addenda.html
%
%   This extends the series given by Krueger (1912) to sixth order in the
%   flattening.  This is a substantially better series than that used by
%   the MATLAB mapping toolbox.  In particular the errors in the projection
%   are less than 5 nanometers withing 3900 km of the central meridian (and
%   less than 1 mm within 7600 km of the central meridian).  The mapping
%   can be continued accurately over the poles to the opposite meridian.
%
%   This routine depends on the MATLAB File Exchange package "Geodesics on
%   an ellipsoid of revolution":
%
%     http://www.mathworks.com/matlabcentral/fileexchange/39108
%
%   See also GEODPROJ, UTM_FWD, TRANMERC_INV.

% Copyright (c) Charles Karney (2012) <charles@karney.com>.
%
% This file was distributed with GeographicLib 1.29.

  if nargin < 4, error('Too few input arguments'), end
  lon0 = -183 + 6 * zone; lat0 = 0;
  fe = 500e3; fn = cvmgt(0, 10000e3, logical(northp)); k0 = 0.9996;
  x = (x - fe) / k0; y = (y - fn) / k0;
  [lat, lon, gam, k] = tranmerc_inv(lat0, lon0, x, y);
  k = k * k0;
end


text:
function formulas = printRxnFormula(model, varargin)
% Prints the reaction formulas for a list of reactions
%
% Reactions that have an upperbound <= 0 and lowerbound < 0 will have
% its directionality reversed unless directionFlag = false.
%
% USAGE:
%
%    formulas = printRxnFormula(model, varargin)
%
% INPUTS:
%    model:             COBRA model structure
%
% OPTIONAL INPUTS:
%
%    varargin:           Optional Inputs provided as 'ParameterName', Value
%                        pairs. the following parameternames are available:
%
%                       * rxnAbbrList:       Cell array of rxnIDs to be printed (Default = print all reactions)
%                       * printFlag:        Print formulas or just return them (Default = true)
%                       * lineChangeFlag:    Append a line change at the end of each line
%                                            (Default = true)
%                       * metNameFlag:       Print full met names instead of abbreviations
%                                            (Default = false)
%                       * fid:               Optional file identifier for printing in files
%                                            (default 1, i.e. stdout)
%                       * directionFlag:     Checks directionality of reaction. See Note.
%                                            (Default = false)
%                       * gprFlag:           Print gene protein reaction association
%                                            (Default = false)
%                       * proteinFlag:       Print the protein names associated with the genes in the
%                                            GPRs associated with the reactions. (Default = false)
%                       * printBounds:       Print the upper and lower Bounds of the reaction (Default = false)
%
%
% OUTPUT:
%    formulas:          Cell array containing formulas of specified reactions
%
% EXAMPLES:
%    1) print only ATPM and TKT1:
%       printRxnFormula(model,'rxnAbbrList',{'ATPM','TKT1'});
%    2) print the reactions of the model with the metabolite Names instead of ids.
%       printRxnFormula(model, 'metNameFlag',true);
%    3) print all reactions with the metabolite names to a given fileID
%       printRxnFormula(model, 'metNameFlag',true, 'fid', fileID);
%
% .. Authors:
%       - Markus Herrgard 11/17/05
%       - Ronan Fleming 04/30/08 (altered code since findRxnIDs used
%                                abbreviations not names of reactions)
%       - Jeff Ortn 10/11/09 (added metNameFlag option)
%       - Richard Que 03/10/10 (added lb < 0 requirement for reversing
%                              directionality)
%       - Ronan Fleming 21/11/14 (printing gpr optional)
%       - Thomas Pfau 15/12/14 (corrected line end)
%       - Ronan Fleming 16/07/16 (directionality not flipped by default
%                                anymore)
%       - Thomas Pfau May 2017 - Changed to Parameter value pair input

optionalParameters = {'rxnAbbrList','printFlag', 'lineChangeFlag', 'metNameFlag', 'fid', 'directionFlag', 'gprFlag', 'proteinFlag','printBounds'};
if (numel(varargin) > 0 && (~ischar(varargin{1}) || ~any(ismember(varargin{1},optionalParameters))))
    %We have an old style thing....
    %Now, we need to check, whether this is a formula, or a complex setup
        tempargin = cell(1,2*(numel(varargin)));

        for i = 1:numel(varargin)

                tempargin{2*(i-1)+1} = optionalParameters{i};
                tempargin{2*(i-1)+2} = varargin{i};
        end
        varargin = tempargin;

end

if ~isfield(model,'lb')
    model.lb=ones(size(model.S,2),1)*-inf;
end
if ~isfield(model,'ub')
    model.lb=ones(size(model.S,2),1)*inf;
end

parser = inputParser();
parser.addRequired('model',@isstruct) % we only check, whether its a struct, no details for speed
parser.addParamValue('rxnAbbrList',model.rxns,@(x) iscell(x) || ischar(x))
parser.addParamValue('printFlag',true,@(x) isnumeric(x) || islogical(x))
parser.addParamValue('lineChangeFlag',true,@(x) isnumeric(x) || islogical(x));
parser.addParamValue('metNameFlag',false,@(x) isnumeric(x) || islogical(x));
parser.addParamValue('fid',1, @isnumeric);
parser.addParamValue('directionFlag',false,@(x) isnumeric(x) || islogical(x));
parser.addParamValue('gprFlag',false,@(x) isnumeric(x) || islogical(x));
parser.addParamValue('proteinFlag',false,@(x) isnumeric(x) || islogical(x));
parser.addParamValue('printBounds',false,@(x) isnumeric(x) || islogical(x));

parser.parse(model,varargin{:})

model = parser.Results.model;
rxnAbbrList = parser.Results.rxnAbbrList;
printFlag = parser.Results.printFlag;
lineChangeFlag = parser.Results.lineChangeFlag;
metNameFlag = parser.Results.metNameFlag;
fid = parser.Results.fid;
directionFlag = parser.Results.directionFlag;
gprFlag = parser.Results.gprFlag;
proteinFlag = parser.Results.proteinFlag;
printBounds = parser.Results.printBounds;

if proteinFlag && ~isfield(model,'proteins')
    %If no proteins field is present, we will use the genes field.
    fprintf('Proteins requested, but no proteins Field exists in the model, using genes instead!\n');
    model.proteins = model.genes;
end

if metNameFlag && ~isfield(model,'metNames')
    %if we want to print the metNames, but they don't exist, just use the mets instead.
    fprintf('metNames requested, but no metNames Field exists in the model, using mets instead!\n');
    model.metNames = model.mets;
end

if (~iscell(rxnAbbrList))
    if (strcmp(rxnAbbrList, 'all'))
        rxnAbbrList = model.rxns;
    else
        rxnAbbrTmp = rxnAbbrList;
        clear rxnAbbrList;
        rxnAbbrList{1} = rxnAbbrTmp;
    end
end

formulas = cell(size(rxnAbbrList));

for i = 1:length(rxnAbbrList)

    rxnAbbr = rxnAbbrList{i};

    rxnID = findRxnIDs(model, rxnAbbr);

    if (rxnID > 0)

        Srxn = full(model.S(:, rxnID));

        if directionFlag && (isfield(model, 'ub') && model.ub(rxnID) <= 0) && (isfield(model, 'lb') && model.lb(rxnID) < 0)
            Srxn = -Srxn;
        end

        Sprod = (Srxn(Srxn > 0));
        if metNameFlag
            prodMets = model.metNames(Srxn > 0);
        else
            prodMets = model.mets(Srxn > 0);
        end

        Sreact = (Srxn(Srxn < 0));
        if metNameFlag
            reactMets = model.metNames(Srxn < 0);
        else
            reactMets = model.mets(Srxn < 0);
        end

        formulaStr = '';
        for j = 1:length(reactMets)
            if (j > 1)
                formulaStr = sprintf('%s+ ', formulaStr);
            end
            if (abs(Sreact(j)) ~= 1)
                formulaStr = sprintf('%s%g %s ', formulaStr, abs(Sreact(j)), reactMets{j});
            else
                formulaStr = sprintf('%s%s ', formulaStr, reactMets{j});
            end
        end

        if (model.lb(rxnID) < 0)
            formulaStr = sprintf('%s <=> ', formulaStr);
        else
            formulaStr = sprintf('%s -> ', formulaStr);
        end

        for j = 1:length(prodMets)
            if (j > 1)
                formulaStr = sprintf('%s+ ', formulaStr);
            end
            if (Sprod(j) ~= 1)
                formulaStr = sprintf('%s%g %s ', formulaStr, abs(Sprod(j)), prodMets{j});
            else
                formulaStr = sprintf('%s%s ', formulaStr, prodMets{j});
            end
        end
        formulas{i} = formulaStr;
        if (printFlag)
            formulaStr = regexprep(formulaStr, ' <=> ', '\t<=>\t');
            formulaStr = regexprep(formulaStr, ' -> ', '\t->\t');
            fprintf(fid, '%s\t%s', rxnAbbr, formulaStr);
            if gprFlag
                if (isfield(model, 'grRules'))
                    fprintf(fid, '\t%s', model.grRules{rxnID});
                elseif (isfield(model,'rules'))
                    rule = regexprep(model.rules{rxnID},'|','or');
                    rule = regexprep(rule,'&','and');
                    rule = regexprep(rule,'x\((?<id>[0-9]+)\)','${model.genes{num2str($1)}}');
                    fprintf(fid, '\t%s', rule);
                end
            end
            if proteinFlag
                if (isfield(model, 'rules'))
                    rule = regexprep(model.rules{rxnID},'|','or');
                    rule = regexprep(rule,'&','and');
                    rule = regexprep(rule,'x\((?<id>[0-9]+)\)','${model.proteins{num2str($1)}}');
                    fprintf(fid, '\t%s', rule);
                end
            end
            if printBounds
                fprintf('\tlb:%f\tub:%f',model.lb(rxnID),model.ub(rxnID));
            end
        end
    else
        if (printFlag)
            fprintf(fid, 'not in model');
        end
        formulaStr = 'NA';
        formulas{i} = formulaStr;
    end
    if (lineChangeFlag && printFlag)
        fprintf(fid, '\n');
    end
end

% %pass out a character string if only one reaction in the abbreviation list
% if length(rxnAbbrList)==1
%     formulas=formulas{1};
% end

author: opencobra

repo: cobratoolbox


text:
function [data, opt] = GrayRawPixelExtractor(im, tmpl, opt)
sz = opt.FeatureExtractor.tmplsize;
if (ndims(im) == 3)
    im = rgb2gray(im);
end
data.tmpl = tmpl;

features = zeros(prod(sz)*size(im, 3), size(tmpl, 1));

% pad
minW = min(round(tmpl(:, 1) - tmpl(:, 3) / 2)) - 1;
maxW = max(round(tmpl(:, 1) + tmpl(:, 3) / 2)) + 1;
minH = min(round(tmpl(:, 2) - tmpl(:, 4) / 2)) - 1;
maxH = max(round(tmpl(:, 2) + tmpl(:, 4) / 2)) + 1;
[h, w, c] = size(im);
if (minW < 1)
    im_new = zeros(h, w + abs(minW) + 1, c);
    im_new(:, abs(minW) + 2:end, :) = im;
    im = im_new;
    tmpl(:, 1) = tmpl(:, 1) + abs(minW) + 1;
end
if (maxW > w)
    im_new = zeros(h, size(im, 2) + maxW - w, c);
    im_new(:, 1:size(im, 2), :) = im;
    im = im_new;
end
if (minH < 1)
    im_new = zeros(h + abs(minH) + 1, size(im, 2), c);
    im_new(abs(minH) + 2:end, :, :) = im;
    im = im_new;
    tmpl(:, 2) = tmpl(:, 2) + abs(minH) + 1;
end
if (maxH > h)
    im_new = zeros(size(im, 1) + maxH - h, size(im, 2), c);
    im_new(1:size(im, 1), :, :) = im;
    im = im_new;
end

for i = 1:size(tmpl, 1)
    midW    = tmpl(i, 1);
    midH    = tmpl(i, 2);
    w       = tmpl(i, 3);
    h       = tmpl(i, 4);

    tempIm = im(round(midH-h/2) : round(midH+h/2),...
                round(midW-w/2) : round(midW+w/2), :);
%     tempIm = imresize(tempIm, sz);
    tempIm = mexResize(tempIm, sz, 'auto');
    features(:, i) = tempIm(:);
    if (norm(features(:, i)) > 1e-6)
        features(:, i) = features(:, i) / norm(features(:, i));
    end
end
% features = features - 0.5;

data.feat = features;

author: flyers

repo: drone-tracking