text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Reconstruction of N-th order nonuniformly sampled bandlimited signals
% using digital filter banks
% Authors: S. K. Sindhi, K. M. M. Prabhu
%%%%%%%%%%%%%%%%%%%%%%%  with Knab window   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LF=13;w_c=0.85; BETA=4 Proposed/Prendergast and BETA=3 ITAMI gives best %
%%%%%%%%%%%%%%%%%%%%%%%  with kaiser_mine1  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LF=13;w_c=0.85; BETA=3 Proposed/Prendergast and BETA=3 ITAMI gives best %

% clear all;
% close all;
% clc;

display('Nth (2) order Reconstruction');

N = 2;                  % Nth order nonuniform sampling
TQ = 1;Fs = 1/TQ;                 % Nyquist Period

T = [1.5*TQ 3*TQ];      % Decimation Periods
K = 0.5*lcm(2*T(1), 2*T(2))/TQ; % number of samples in recurrent period
capT = K*TQ; % the full sampling period - of all samplers
M = capT./T;
ML = 400; % number of slices
w_c = 0.85;
NS = 100;  % Number of Sinusoids

LF = lcm(M(1),M(2))*2*K+1;      % length of LF should be Multiple of LCM{M(p)}*2*K
n = -(LF-1)/2:1:(LF-1)/2;
Hd = firpm(LF-1,[0 w_c],[0 w_c*pi],'differentiator');
delayV = (LF-1)/2;

k = -(K-1):1:(K-1);
m = (0:1:(2*K-1))';
F = exp(1i*(pi/K).*kron(m,k));
std = [1e-6 1e-5 1e-4 1e-3 1e-2 1e-1];%5*1e-1];

serP = zeros(size(std));
serE = zeros(size(std));
serI = zeros(size(std));
serV = zeros(size(std));
serPr = zeros(size(std));
serJ = zeros(size(std));
serNO = zeros(size(std));

MCruns = 25;
MCruns1 = 25;

for tt = 1:length(std)
aa = 0;
display(tt);
for rrr = 1:MCruns1
taus = [0 1+std(tt)*randn]*TQ;
if or(taus(2)==1.5*TQ,taus(1)==taus(2))
    aa = aa+1;
    continue;
end
tausI = sort([taus(1) taus(2) T(1)+taus(1)]);
tauI = zeros(K,ML);
for p = 1:K
    tauI(p,:) = tausI(p)+(0:ML-1)*capT;
end;
a = zeros(1,K);
for p = 1:K
    a(p) = 1;
    for q = 1:K
        if q ~= p
                a(p) = a(p)/sin(pi*(tausI(p)-tausI(q))/capT);
        end;
    end;
end;
% c = sin(pi*(tausI(0+1))/capT);
% s = cos(pi*(tausI(0+1))/capT);
% b(1,1) = 0.5*(c+1i*s);
% c = sin(pi*(tausI(1+1))/capT);
% s = cos(pi*(tausI(1+1))/capT);
% b(1,2) = 0.5*(c+1i*s);
% b(3,:) = conj(b(1,:));
% b(2,:) = 0;

b = zeros(2*K-1,K);
c = -0.5*cos(pi*(tausI(1+1)+tausI(2+1))/capT);
s = -0.5*sin(pi*(tausI(1+1)+tausI(2+1))/capT);
b(1,1) = 0.5*(c+1i*s);
c = -0.5*cos(pi*(tausI(0+1)+tausI(2+1))/capT);
s = -0.5*sin(pi*(tausI(0+1)+tausI(2+1))/capT);
b(1,2) = 0.5*(c+1i*s);
c = -0.5*cos(pi*(tausI(0+1)+tausI(1+1))/capT);
s = -0.5*sin(pi*(tausI(0+1)+tausI(1+1))/capT);
b(1,3) = 0.5*(c+1i*s);
b(5,:) = conj(b(1,:));
b(2,:) = 0;
b(4,:) = conj(b(2,:));
c = 0.5*cos(pi*(tausI(2+1)-tausI(1+1))/capT);
b(3,1) = c;
c = 0.5*cos(pi*(tausI(2+1)-tausI(0+1))/capT);
b(3,2) = c;
c = 0.5*cos(pi*(tausI(1+1)-tausI(0+1))/capT);
b(3,3) = c;

r = tausI-TQ*(0:K-1);
rr = r(mod((0:ML*K-1),K)+1);

tauPr = zeros(K,ML);
for p = 1:K
    tauPr(p,:) = -tausI(p)+(0:ML-1)*capT;
end;
H = zeros(K,LF);
for i=1:K
%     H(i,:) = sinc(n-tausI(i)).*conv(sinc(n-tausI(i)),kaiser(LF,10).','same');
%     H(i,:) = sinc(n-tausI(i)).*knab(LF,3,-tausI(i)).';
    H(i,:) = sinc(n-tausI(i)).*kaiser_mine1(LF,3,-tausI(i));
end;
H = H(:,1:end-1);
EP = reshape(H.',K,length(H(1,:))/K,K);
capE = [];
for k=1:K
    temp = [];
    for i=1:K
        temp=[temp,toeplitz([EP(i,1,k),zeros(1,(length(H(1,:))/K)-1)],[EP(i,:,k),zeros(1,(length(H(1,:))/K)-1)])];
    end
    capE = [capE;temp];
end
d = ceil(size(capE,2)/(2*K));
P = kron(eye(K),[zeros(1,d),1,zeros(1,(size(capE,2)/K)-d-1)]);
R = P/capE;
% size(capE)
% size(zeros(LF-1,2*LF-2-K))
% size(R)
% size(zeros(K,LF-1))
% size(P)
% size(zeros(K,2*LF-2-K))
R = upsample(R.',K).';
rt = size(R,2)/K;
FR = zeros(K,rt);
for j=1:K
    for i=1:K
        temp = filter([zeros(1,K-i),1],1,R(i,(j-1)*rt+1:j*rt));
        FR(j,:) = FR(j,:)+temp;
    end
end

r = r.';
w_o = w_c*pi*TQ;
hJ = zeros(K,LF);
C = zeros(1,LF);
Nt = (LF-1)/2;
for i = 1:K
    C = -2*sin(w_o*(n-r(1+(mod(i-1-n,K)))'))./(pi*(n-r(1+(mod(i-1-n,K)))'));
    C(isnan(C)==1)=-2*w_o/pi;
    C = C.';
    S = zeros(LF,LF);
    for k = 1:LF
        S(k,:) = sin(w_o*(-Nt+k-1-r(1+(mod(i-1-(-Nt+k-1),K)))-(n-r(1+(mod(i-1-n,K)))')))./(pi*(-Nt+k-1-r(1+(mod(i-1-(-Nt+k-1),K)))-(n-r(1+(mod(i-1-n,K)))')));
    end;
    S(isnan(S)==1)=w_o/pi;
    hJ(i,:) = -0.5*S\C;
end;

for pp = 1:MCruns;
Frq = rand(1,NS)*w_c/2;
Amp = rand(1,NS)/(sqrt(NS)*2);
Phi = rand(1,NS)*2*pi;
input = zeros(1,ML*K);
for k = 1:NS
  input = input + Amp(k)*sin(2*pi*Frq(k)*(0:ML*K-1)*TQ+Phi(k));
end;
xp = zeros(N,ML*K);
for p = 1:N
    tau = taus(p)+(0:ML*M(p)-1)*T(p);
    x1 = zeros(1,ML*M(p));
    for k = 1:NS
        x1 = x1 + Amp(k)*sin(2*pi*Frq(k)*tau+Phi(k));
    end;

    m = (0:1:M(p)-1)'; lemda = 0:1:M(p)-1;
    W = exp(1i*(2*pi/M(p)).*kron(m,lemda)); % m*lemda

    bb = zeros(2*(K-M(p))+1,M(p));
    aaa = ones(1,M(p));
    for l = 1:M(p)
        for q = 1:N
            if q ~= p
                    aaa(l) = aaa(l)/sin(pi*M(q)*(taus(p)-taus(q)+(l-1)*T(p))/capT);
                    c = sin(pi*M(q)*(-taus(q)+(l-1)*T(p))/capT);
                    s = cos(pi*M(q)*(taus(q)-(l-1)*T(p))/capT);
                    bb(2*M(q)+1,l) = 0.5*(c-1i*s);
                    bb(1,l) = conj(bb(2*M(q)+1,l));
            end;
        end;
    end;
    A = diag(aaa); % display(A);
    B = bb; % display(B);

    y1 = upsample(x1,K);
    y1 = reshape(y1,M(p),length(y1)/M(p));

    if M(p)>1
        y1(2:end,:) = flipud(y1(2:end,:));
        for i = 1:M(p)-1
            y1(i+1,:) = filter([0,1],1,y1(i+1,:));
        end;
    end;

    dim = K-M(p); w = -dim:1:dim;

    xlemda = zeros(M(p),size(y1,2));
    for lemda = 0:M(p)-1

        rP = (lemda/M(p))+(0:1:(2*K-1))';
        Fshift = exp(1i*(pi/K).*kron(rP,w));   % r*w

%         h = sinc((n*TQ/T(p))+(lemda/K)-(taus(p)/T(p))).*conv(sinc((lemda/M(p))-(taus(p)/TQ)), kaiser(LF,3),'same');
        h = sinc((n*TQ/T(p))+(lemda/K)-(taus(p)/T(p))).*kaiser_mine1(LF,3,(lemda/M(p))-(taus(p)/TQ));
%        h = sinc((n*TQ/T(p))+(lemda/K)-(taus(p)/T(p))).*knab(LF,3,(lemda/M(p))-(taus(p)/TQ)).';
        y2 = Fshift*B*A*W*W(:,lemda+1)*y1(lemda+1,:);
        for i=1:2*K
            h1 = upsample(downsample(h,2*K,i-1),2*K);
            y2(i,:) = filter(h1,1,y2(i,:));
            y2(i,:) = filter([zeros(1,i-1),1],1,y2(i,:));%,zeros(1,2*K-i)
        end;
        xlemda(lemda+1,:) = sum(y2,1);
    end;
    xp(p,:) = sum(xlemda,1)/M(p);
    clear bb;
end;
y = real(sum(xp,1));
delayP = (length(h)-1)/2;
y = y(1+delayP:end);
x = input(1:end-delayP);
y = y(160:end-60);
x = x(160:end-60);
serP(tt) = serP(tt)+20*log10(norm(x,2)/norm(y-x,2));

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % Filterbank Reconstruction of Bandlimited Signals from Nonuniform and
% % Generalized Samples
% % Authors: Y C Eldar and A V Oppenheim
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

y = zeros(N,ML*K);
for p = 1:N
    tau = taus(p)+(0:ML*M(p)-1)*T(p);
    x1 = zeros(1,ML*M(p));
    for k = 1:NS
        x1 = x1 + Amp(k)*sin(2*pi*Frq(k)*tau+Phi(k));
    end;
    y1 = upsample(x1,K);

    LFE = M(p)*lcm(M(1),M(2))*2*K+1;      % length of LF should be Multiple of LCM{M(p)}*2*K
    nE = -(LFE-1)/2:1:(LFE-1)/2;
    h = sinc((nE/K)-(taus(p)/T(p))).*kaiser_mine1(LFE,3,-K*(taus(p)/T(p)));

%%%%%%%%%%%%%%%%%%%%%%%%%%% Implementation 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%
%     bb = zeros(2*(K-M(p))+1,M(p));
%     aaa = ones(1,M(p));
%     for l = 1:M(p)
%         for q = 1:N
%             if q ~= p
%                 aaa(l) = aaa(l)/sin(pi*(taus(p)-taus(q)+(l-1)*T(p))/T(q));
%                 c = sin(pi*M(q)*(-taus(q)+(l-1)*T(p))/capT);
%                 s = cos(pi*M(q)*(taus(q)-(l-1)*T(p))/capT);
%                 bb(2*M(q)+1,l) = 0.5*(c-1i*s);
%                 bb(1,l) = conj(bb(2*M(q)+1,l));
%             end;
%         end;
%     end;
%     AE = diag(aaa);
%     BE = bb.';
%
%     dim = K-M(p); w = -dim:1:dim;
%     FE = exp(1i*(pi/(K*M(p))).*kron(nE,w'));
%
%     E1E = exp(1i*(2*pi/M(p)).*kron((0:M(p)-1),(0:M(p)-1)'))/M(p);
%
%     E2E = exp(1i*(2*pi/M(p)).*kron(nE,(0:M(p)-1)'));
%
%     temp = E1E*AE*BE*FE;
%     bbn = h.*sum(E2E.*temp,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%% Implementation 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%
%     bb = zeros(2*(K-M(p))+1,M(p));
%     aaa = ones(1,M(p));
%     bbvl = zeros(M(p),LFE);
%     for l = 1:M(p)
%         for q = 1:N
%             if q ~= p
%                 aaa(l) = aaa(l)*sin(pi*(taus(p)-taus(q)+(l-1)*T(p))/T(q));
%                 c = sin(pi*M(q)*(-taus(q)+(l-1)*T(p))/capT);
%                 s = cos(pi*M(q)*(taus(q)-(l-1)*T(p))/capT);
%                 bb(2*M(q)+1,l) = 0.5*(c-1i*s);
%                 bb(1,l) = conj(bb(2*M(q)+1,l));
%             end;
%         end;
%         bbv = zeros(2*(K-M(p))+1,LFE);
%         for v = -(K-M(p)):K-M(p)
%             bbv(K-M(p)+1+v,:) = bb(K-M(p)+1+v,l)*exp(1i*(pi/(K*M(p)))*v*nE);
%         end;
%         bbvl(l,:) = sum(bbv,1)/aaa(l);
%     end;
%
%     bbb = zeros(M(p),LFE);
%     bbn = zeros(M(p),LFE);
%     for m = 1:M(p)
%         for l = 1:M(p)
%             bbb(l,:) = bbvl(l,:)*exp(1i*(2*pi/M(p))*(m-1)*(l-1));
%         end
%         bbb = sum(bbb,1);
%         bbn(m,:) = bbb.*exp(1i*(2*pi/M(p))*(m-1)*nE);
%     end;
%     bbn = sum(bbn,1);
%     bbn = bbn.*h/M(p);

%%%%%%%%%%%%%%%%%%%%%%%%%%% Implementation 3 %%%%%%%%%%%%%%%%%%%%%%%%%%%
    aaa = ones(1,M(p));
    bb = ones(M(p),LFE);
    for l = 1:M(p)
        for q = 1:N
            if q ~= p
                aaa(l) = aaa(l)*sin(pi*(taus(p)-taus(q)+(l-1)*T(p))/T(q));
                bb(l,:) = bb(l,:).*sin(pi*((nE*TQ/M(p))-taus(q)+(l-1)*T(p))/T(q));
            end;
        end;
        bb(l,:) = bb(l,:)/aaa(l);
    end;

    bbb = zeros(M(p),LFE);
    bbn = zeros(M(p),LFE);
    for m = 1:M(p)
        for l = 1:M(p)
            bbb(l,:) = bb(l,:)*exp(1i*(2*pi/M(p))*(m-1)*(l-1));
        end
        bbb = sum(bbb,1);
        bbn(m,:) = bbb.*exp(1i*(2*pi/M(p))*(m-1)*nE);
    end;
    bbn = sum(bbn,1);
    bbn = bbn.*h/M(p);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    y1 = conv(y1,bbn);
    delay = (length(h)-1)/2;
    y(p,:) = y1(1+delay:M(p):end-delay);
end;
y = sum(real(y),1);
x = input;
y = y(160:end-60);
x = x(160:end-60);
serE(tt) = serE(tt)+20*log10(norm(x,2)/norm(y-x,2));

% % figure();
% subplot(2,1,1);
% plot(([x' y']));
% title('input / output signals');
% xlabel('sample');
% ylabel('signal value');
% grid on;
% subplot(2,1,2);
% plot((x'-y'));
% xlabel('time (sample)');
% ylabel('error value');
% grid on;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% A realization of Digital Filter Banks for Reconstruction of Uniformly
% sampled signals from nonuniform samples
% Authors: Itami, Watanabe, Nishihara
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

x11 = zeros(K,ML);
for k = 1:NS
    x11 = x11 + Amp(k)*sin(2*pi*Frq(k)*tauI+Phi(k));
end;
y1=upsample(x11.',K).';

y = zeros(K,size(y1,2));
for r = 1:K
    y2 = F*b(:,r)*a(r)*y1(r,:);
%     h = sinc((n/K)-tausI(r)/capT).*conv(sinc(-tausI(r)/TQ), kaiser(LF,3),'same');
    h = sinc((n/K)-tausI(r)/capT).*kaiser_mine1(LF,3,-tausI(r)/TQ);
%     h = sinc((n/K)-tausI(r)/capT).*knab(LF,3,-tausI(r)/TQ).';
    for i=1:2*K
            h1 = upsample(downsample(h,2*K,i-1),2*K);
            y2(i,:) = filter(h1,1,y2(i,:));
            y2(i,:) = filter([zeros(1,i-1),1],1,y2(i,:));%,zeros(1,2*K-i)
    end;
    y(r,:) = sum(y2,1);
end;
y = (sum(y,1));
delayI = (length(n)-1)/2;
x=input(1:end-delayI);
y=y(1+delayI:end);
y = y(160:end-60);
x = x(160:end-60);
serI(tt) = serI(tt)+20*log10(norm(x,2)/norm(y-x,2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Without any reconstruction, the SNR value calculation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wo_reconst = reshape(x11,1,[]);
serNO(tt) = serNO(tt)+20*log10(norm(input(160:end-60),2)/...
            norm(wo_reconst(160:end-60)-input(160:end-60),2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Reconstruction of Nonuniformly Sampled Band-Limited Signals
% Using a Differentiator-Multiplier Cascade
% Authors: Stefan Tertinek and Christian Vogel
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

x1 = reshape(x11,1,K*size(x1,2));
% x1 = x11;% LF = (0:11); Fs = 1;
% Differentiator Design
% figure();
% NFFT = 2^nextpow2(length(Hd)); % Next power of 2 from length of Hd
% HD = fftshift(fft(Hd,NFFT))/length(Hd);
% f = Fs*linspace(-1,1,NFFT);
% % Plot double-sided amplitude spectrum.
% plot(f,2*abs(HD(1:NFFT)))
% title('Double-Sided Amplitude Spectrum of Hd(n)')
% xlabel('Frequency (Hz)')
% ylabel('|HD(f)|')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
y1 = filter(Hd,1,x1);
x1 = filter([zeros(1,delayV),1],1,x1);
r2 = filter([zeros(1,delayV),1],1,rr);
e = y1.*r2;
y1 = x1-e;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
y2 = filter(Hd,1,y1);
temp = filter([zeros(1,delayV),1],1,y2);
x1 = filter([zeros(1,2*delayV),1],1,x1);
r2 = filter([zeros(1,2*delayV),1],1,r2);
e1 = temp.*r2;
y2 = filter(Hd,1,y2);
e2 = 0.5*y2.*r2.^2;
y2 = x1-e1-e2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
y3 = filter(Hd,1,y2);
temp = filter([zeros(1,2*delayV),1],1,y3);
x1 = filter([zeros(1,3*delayV),1],1,x1);
r2 = filter([zeros(1,3*delayV),1],1,r2);
e1 = temp.*r2;
y3 = filter(Hd,1,y3);
temp = filter([zeros(1,delayV),1],1,y3);
e2 = 0.5*temp.*r2.^2;
y3 = filter(Hd,1,y3);
e3 = (y3.*r2.^3)/6;
y3 = x1-e1-e2-e3;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
y = real(y3(1+6*delayV:end));
x = input(1:end-6*delayV);
y = y(160:end-60);
x = x(160:end-60);
serV(tt) = serV(tt)+20*log10(norm(x,2)/norm(y-x,2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Reconstruction of Band-Limited Periodic Nonuniformly Sampled Signals
% Through Multirate Filter Banks
% Ryan S Prendergast, Bernard C Levy, Paul J Hurst
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

x1 = zeros(K,ML);
for k = 1:NS
    x1 = x1 + Amp(k)*sin(2*pi*Frq(k)*tauPr+Phi(k));
end;

yb = upsample(x1.',K).';
for i = 1:K
    yb(i,:) = filter(FR(i,:),1,yb(i,:));
end;
y = sum(yb,1);
delayPr = (size(FR,2))/2+K-1;
x=input(1:end-delayPr);
y=y(1+delayPr:end);
y = y(160:end-60);
x = x(160:end-60);
serPr(tt) = serPr(tt)+20*log10(norm(x,2)/norm(y-x,2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Reconstruction of Periodically Nonuniformly Sampled Bandlimited Signals
% Using Time-Varying FIR Filters
% Authors: H. Johansson and Per Lowenborg
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

x1 = reshape(x11,1,size(x11,2)*K);
y1 = zeros(K,length(x1));
for j=1:K
    y1(j,:) = filter(hJ(j,:),1,x1);
    y1(j,:) = upsample(downsample(y1(j,:),K,j-1),K)/K;
    y1(j,:) = filter([zeros(1,j-1),1],1,y1(j,:));
end;
y = K*0.25*sum(y1,1);
delayJ = (size(hJ,2)-1)/2;
y = real(y(1+delayJ:end));
x = input(1:end-delayJ);
y = y(160:end-60);
x = x(160:end-60);
serJ(tt) = serJ(tt)+20*log10(norm(x,2)/norm(y-x,2));
end
end
end
serP = serP/(MCruns*(MCruns1-aa));
serE = serE/(MCruns*(MCruns1-aa));
serI = serI/(MCruns*(MCruns1-aa));
serV = serV/(MCruns*(MCruns1-aa));
serPr = serPr/(MCruns*(MCruns1-aa));
serJ = serJ/(MCruns*(MCruns1-aa));
serNO = serNO/(MCruns*(MCruns1-aa));

figure();hold all;
% plot(std,serJ,'kp-','LineWidth',2);
% plot(std,serPr,'ko-','LineWidth',2);
% plot(std,serV,'ks-','LineWidth',2);
% plot(std,serP,'kd-','LineWidth',2);
% plot(std,serI,'k>-','LineWidth',2);
% plot(std,serE,'k+-','LineWidth',2);
% plot(std,serNO,'k+-','LineWidth',2);
plot(std,serJ);
plot(std,serPr);
plot(std,serV);
plot(std,serP);
plot(std,serI);
plot(std,serE);
plot(std,serNO);
% legend('Johansson','Prendergast','Tertinek','Proposed','Itami','Eldar');
legend('Johansson','Prendergast','Tertinek','Proposed','Itami','Eldar', 'W/O Reconst');
xlabel('Standard Deviation (\sigma)','fontsize',14,'fontweight','b');
ylabel('SNR in dB','fontsize',14,'fontweight','b');
grid on;box on;
set(gca,'fontsize',14,'fontweight','b')

% display('Proposed Method');
% display(sprintf('Delay imposed due to reconstruction system = %d', delayP));
% display(sprintf('Length of prototype filter = %d\n', length(n)));
%
% display('Itami Method');
% display(sprintf('Delay imposed due to reconstruction system = %d', delayI));
% display(sprintf('Length of prototype filter = %d', length(n)));
% display(sprintf('Length of fractional delay filter = %d\n', length(n)));
%
% display('Vogel Method');
% display(sprintf('Delay imposed due to reconstruction system = %d', 6*delayV));
% display(sprintf('Length of differentiator = %d\n', LF));
%
% display('Prendergast Method');
% display(sprintf('Delay imposed due to reconstruction system = %d', delayPr));
% display(sprintf('Length of prototype filter = %d', size(FR,2)));
% display(sprintf('Length of fractional delay filter = %d\n', length(n)));
%
% display('Johansson Method');
% display(sprintf('Delay imposed due to reconstruction system = %d', delayJ));
% display(sprintf('Length of prototype filter = %d\n', size(h,2)));

% figure();
% subplot(2,1,1);
% plot(([x' y']));
% title('input / output signals');
% xlabel('sample');
% ylabel('signal value');
% grid on;
% subplot(2,1,2);
% plot((x'-y'));
% xlabel('time (sample)');
% ylabel('error value');
% grid on;

author: yueyuzhao

repo: gyrophone

save_path: github-repos/MATLAB/yueyuzhao-gyrophone


text:
%STEREOCALIBRATE  Calibrates the stereo camera
%
%     S = cv.stereoCalibrate(objectPoints, imagePoints1, imagePoints2, imageSize)
%     [...] = cv.stereoCalibrate(..., 'OptionName', optionValue, ...)
%
% ## Input
% * __objectPoints__ A cell array of cells of calibration pattern points in
%   the calibration pattern coordinate space  `{{[x,y,z], ..}, ...}`.
% * __imagePoints1__ A cell array of cells of the projections of calibration
%   pattern points `{{[x,y], ..}, ...}`, observed by the first camera.
% * __imagePoints2__ A cell array of cells of the projections of calibration
%   pattern points `{{[x,y], ..}, ...}`, observed by the second camera.
% * __imageSize__ Size of the image used only to initialize the intrinsic
%   camera matrix `[w,h]`.
%
% ## Output
% * __S__ scalar struct having the following fields:
%   * __cameraMatrix1__ output first camera matrix
%     `A = [fx1 0 cx1; 0 fy1 cy1; 0 0 1]`.
%   * __distCoeffs1__ output vector of distortion coefficients
%     `[k1,k2,p1,p2,k3,k4,k5,k6,s1,s2,s3,s4,taux,tauy]` of 4, 5, 8, 12
%     or 14 elements. The output vector length depends on the options.
%   * __cameraMatrix2__ output second camera matrix
%     `A = [fx2 0 cx2; 0 fy2 cy2; 0 0 1]`. The parameter is similar to
%     `cameraMatrix1`.
%   * __distCoeffs2__ output lens distortion coefficients for the second
%     camera. The parameter is similar to `distCoeffs1`.
%   * __R__ output 3x3 rotation matrix between the 1st and the 2nd camera
%     coordinate systems.
%   * __T__ output 3x1 translation vector between the coordinate systems of
%     the cameras.
%   * __E__ output 3x3 essential matrix.
%   * __F__ output 3x3 fundamental matrix.
%   * __reprojErr__ output final re-projection error (scalar).
% * __perViewErrors__ Output matrix of the RMS re-projection error estimated
%   for each pattern view.
%
% ## Options
% * __CameraMatrix1__, __CameraMatrix2__ Initial camera matrices. If any of
%   `UseIntrinsicGuess`, `FixAspectRatio`, `FixIntrinsic` (default), or
%   `FixFocalLength` are specified, some or all of the matrix components must
%   be initialized. See the flags description for details.
% * __DistCoeffs1__, __DistCoeffs2__ Initial lens distortion coefficients.
% * __R__, __T__ Initial extrinsic parameters `R` and `T`.
%   If `UseExtrinsicGuess` is set, these must be initialized. Not set by
%   default.
% * __FixIntrinsic__ Fix `cameraMatrix1`,`cameraMatrix2` and `distCoeffs1`,
%   `distCoeffs2` so that only `R`, `T`, `E`, and `F` matrices are estimated.
%   default true.
% * __UseIntrinsicGuess__ Optimize some or all of the intrinsic parameters
%   according to the specified flags. Initial values are provided by the user.
%   default false.
% * __UseExtrinsicGuess__ `R`, `T` contain valid initial values that are
%   optimized further. Otherwise `R`, `T` are initialized to the median value
%   of the pattern views (each dimension separately). default false
% * __FixPrincipalPoint__ Fix the principal points during the optimization.
%   default false.
% * __FixFocalLength__ Fix `fx1`,`fx2` and `fy1`,`fy2`. default false.
% * __FixAspectRatio__ Optimize `fy1`,`fy2` and fix the ratio `fx1/fy1`,
%   `fx2/fy2`. default false.
% * __SameFocalLength__ Enforce same `fx1=fx2` and `fy1=fy2`. default false.
% * __ZeroTangentDist__ Tangential distortion coefficients for each camera are
%   set to zeros and stay fixed. default false.
% * __FixTangentDist__ The tangential distortion coefficients are not
%   changed during the optimization. If `UseIntrinsicGuess` is set, the
%   coefficient from the supplied `DistCoeffs` matrix is used. Otherwise, it
%   is set to 0. default false.
% * __FixK1__, ..., __FixK6__ The corresponding radial distortion coefficient
%   is not changed during the optimization. If `UseIntrinsicGuess` is set, the
%   coefficient from the supplied `DistCoeffs` matrix is used. Otherwise, it
%   is set to 0. default false.
% * __RationalModel__ Coefficients `k4`, `k5`, and `k6` are enabled. To
%   provide the backward compatibility, this extra flag should be explicitly
%   specified to make the calibration function use the rational model and
%   return 8 coefficients. If the flag is not set, the function computes and
%   returns only 5 distortion coefficients. default false.
%   (`RationalModel` as false implies `FixK4`,`FixK5`,`FixK6` as true).
% * __ThinPrismModel__ Coefficients `s1`, `s2`, `s3` and `s4` are enabled. To
%   provide the backward compatibility, this extra flag should be explicitly
%   specified to make the calibration function use the thin prism model and
%   return 12 coefficients. If the flag is not set, the function computes and
%   returns only 5 distortion coefficients. default false.
%   (`ThinPrismModel` as false implies `FixS1S2S3S4` as true).
% * __FixS1S2S3S4__ The thin prism distortion coefficients are not changed
%   during the optimization. If `UseIntrinsicGuess` is set, the coefficient
%   from the supplied `DistCoeffs` matrix is used. Otherwise, it is set to 0.
%   default false.
% * __TiltedModel__ Coefficients `tauX` and `tauY` are enabled. To provide the
%   backward compatibility, this extra flag should be explicitly specified to
%   make the calibration function use the tilted sensor model and return 14
%   coefficients. If the flag is not set, the function computes and returns
%   only 5 distortion coefficients. default false.
%   (`TiltedModel` as false implies `FixTauXTauY` as true).
% * __FixTauXTauY__ The coefficients of the tilted sensor model are not
%   changed during the optimization. If `UseIntrinsicGuess` is set, the
%   coefficient from the supplied `DistCoeffs` matrix is used. Otherwise, it
%   is set to 0. default false.
% * __UseLU__ Use LU instead of SVD decomposition for solving. Much faster but
%   potentially less precise. default false.
% * __UseQR__ Use QR instead of SVD decomposition for solving. Faster but
%   potentially less precise. default false.
% * __Criteria__ Termination criteria for the iterative optimization algorithm.
%   default `struct('type','Count+EPS', 'maxCount',30, 'epsilon',1e-6)`
%
% The function estimates transformation between two cameras making a stereo
% pair. If you have a stereo camera where the relative position and
% orientation of two cameras is fixed, and if you computed poses of an object
% relative to the first camera and to the second camera, `(R1,T1)` and
% `(R2,T2)`, respectively (this can be done with cv.solvePnP), then those
% poses definitely relate to each other. This means that, given `(R1,T1)`, it
% should be possible to compute `(R2,T2)`. You only need to know the position
% and orientation of the second camera relative to the first camera. This is
% what the described function does. It computes `(R,T)` so that:
%
%     R2 = R * R1
%     T2 = R * T1 + T
%
% Optionally, it computes the essential matrix `E`:
%
%     E = [ 0 -T2  T1;
%          T2   0 -T0;
%         -T1  T0   0] * R
%
% where `Ti` are components of the translation vector `T`: `T = [T0,T1,T2]'`.
% And the function can also compute the fundamental matrix `F`:
%
%     F = inv(cameraMatrix2)' * E * inv(cameraMatrix1)
%
% Besides the stereo-related information, the function can also perform a full
% calibration of each of two cameras. However, due to the high dimensionality
% of the parameter space and noise in the input data, the function can diverge
% from the correct solution. If the intrinsic parameters can be estimated with
% high accuracy for each of the cameras individually (for example, using
% cv.calibrateCamera), you are recommended to do so and then pass
% `FixIntrinsic` flag to the function along with the computed intrinsic
% parameters. Otherwise, if all the parameters are estimated at once, it makes
% sense to restrict some parameters, for example, pass `SameFocalLength` and
% `ZeroTangentDist` flags, which is usually a reasonable assumption.
%
% Similarly to cv.calibrateCamera, the function minimizes the total
% re-projection error for all the points in all the available views from both
% cameras. The function returns the final value of the re-projection error.
%
% See also: cv.calibrateCamera, cv.solvePnP, cv.stereoRectify,
%  estimateCameraParameters, stereoParameters, extrinsics
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function test_suite = test_regression1

%   Run specific demo, save values and compare the results to the expected.
%   Works for both xUnit Test Framework package by Steve Eddins and for
%   the built-in Unit Testing Framework (as of Matlab version 2013b).
%
%   See also
%     TEST_ALL, DEMO_REGRESSION1
%
% Copyright (c) 2011-2012 Ville Tolvanen
% Copyright (c) 2014 Tuomas Sivula

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

  % Check which package to use
  if exist('initTestSuite', 'file')
    initTestSuite_path = which('initTestSuite');
  else
    initTestSuite_path = '';
  end
  if ~isempty(initTestSuite_path) && ...
     exist([initTestSuite_path(1:end-15) 'runtests'], 'file')
    % xUnit package
    initTestSuite;
  else
    % Built-in package
    % Use all functions except the @setup
    tests = localfunctions;
    tests = tests(~cellfun(@(x)strcmp(func2str(x),'setup'),tests));
    test_suite = functiontests(tests);
  end
end

% -------------
%     Tests
% -------------

function testRunDemo(testCase)
  % Run the correspondin demo and save the values. Note this test has to
  % be run at lest once before the other test may succeed.
  rundemo(getName(), { ...
    'K', 'C', 'w', ...
    {'Eft_map', @(x)x(1:50)}, {'Varft_map', @(x)x(1:50)}, ...
    {'Eft_ia', @(x)x(1:50)}, {'Varft_ia', @(x)x(1:50)}, ...
    {'Eft_mc', @(x)x(1:50)}, {'Varft_mc', @(x)x(1:50)} })
end

function testCovarianceMatrices(testCase)
  % Test diagonal with 1e-7 tolerance (because of varying jitterSigma2)
  verifyVarsEqual(testCase, getName(), {'K','C'}, ...
    @(x)diag(x), 'AbsTol', 1e-7)
  % Test non-diagonal elements
  verifyVarsEqual(testCase, getName(), {'K','C'}, ...
    @(x)x(~eye(size(x))), 'AbsTol', eps)
end

function testOptimizedParameter(testCase)
  verifyVarsEqual(testCase, getName(), {'w'}, ...
    'RelTolElement', 0.05)
end

function testPredictedMeanVarianceGrid(testCase)
  verifyVarsEqual(testCase, getName(), ...
    {'Eft_map','Varft_map'}, ...
    'RelTolElement', 0.05, 'RelTolRange', 0.01)
end

function testPredictedMeanVarianceMC(testCase)
  verifyVarsEqual(testCase, getName(), {'Eft_mc'}, ...
    'RelTolElement', 0.1, 'RelTolRange', 0.02)
  verifyVarsEqual(testCase, getName(), {'Varft_mc'}, ...
    'AbsTol', 0.3)
end

function testPredictedMeanVarianceIA(testCase)
  verifyVarsEqual(testCase, getName(), ...
    {'Eft_ia','Varft_ia'}, ...
    'RelTolElement', 0.1, 'RelTolRange', 0.02)
end

% ------------------------
%     Helper functions
% ------------------------

function testCase = setup
  % Helper function to suply empty array into variable testCase as an
  % argument for each test function, if using xUnit package. Not to be
  % used with built-in test framework.
  testCase = [];
end

function name = getName
  % Helperfunction that returns the name of the demo, e.g. 'binomial1'.
  name = mfilename;
  name = name(6:end);
end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function ttf = sin_test()

d = 10;
r = 2;
n = 2;
N = 20000;
x = 2*pi*rand(d,N);
y = sin(sum(x));
tt_rank = [1 ; r*ones(d-1,1) ; 1];
n_basis = n*ones(d,1);
basis_ps = cumsum([1 ; n_basis*N]);
basis_cr = zeros(basis_ps(d+1) - basis_ps(1), 1);
for dim = 1: d
    t = zeros(n, N);
    t(1,:) = sin(x(dim,:));
    t(2,:) = cos(x(dim,:));
	basis_cr(basis_ps(dim):basis_ps(dim+1)-1) = reshape(t, [n*N 1]);
end
coeff = reg_als(basis_cr, y, tt_rank, n_basis);
ttf = tt_function(@fun, coeff);

end

function fx = fun(~, j, x)
% i-th dimension, j-th basis function at x

if j == 1
    fx = sin(x);
elseif j == 2
    fx = cos(x);
end

end

author: oseledets

repo: TT-Toolbox


text:
% The COBRAToolbox: testDetectDeadEnds.m
%
% Purpose:
%     - Tests whether all dead ends are correctly determined from a given
%       model
%
% Authors:
%     - Thomas Pfau, April 2017
%     - CI integration: Laurent Heirendt

% save the current path
currentDir = pwd;

% The testmodel used is structured as follows:
%
%   <=> A -> F ---> H <=> J ->
%        \     /          ^
%         -> E            |
%            | \          L
%            v   -> I     |
%   <=> B -> D       \    v
%        \             -> K <=>
%         -> G <=>
%   <=> C
% Thus, C, D, and L are dead ends, and C is also an exchanged metabolite
% and should thus be detected by removeExternalMets set to true.

fileDir = fileparts(which('testDetectDeadEnds'));
cd(fileDir);

% load the test models
load('DeadEndTestModel','model')

%Determine all dead end metabolites
[mets] = detectDeadEnds(model);
assert(all(ismember(model.mets(mets),{'D','L','C'})));

%When detectDeadEnds is changed according to Ronans suggestion, we need to test
%multiple solvers.
solverPkgs = {'gurobi', 'tomlab_cplex', 'glpk'};

for k = 1:length(solverPkgs)

    % set the solver
    solverOK = changeCobraSolver(solverPkgs{k}, 'LP', 0);

    if solverOK == 1

        %Only determine those, which are not involved in an exchange reaction
        [mets] = detectDeadEnds(model,1);
        assert(all(ismember(model.mets(mets),{'D','L'})));
    end
end
%Clean up after test
clear mets
clear model

cd(currentDir)

author: opencobra

repo: cobratoolbox


text:
classdef SubcellsMesher_Boundary < SubcellsMesher

    properties (GetAccess = protected, SetAccess = private)
        nCellNodes
    end

    methods (Access = protected, Abstract)

        computeFacetsConnectivities(obj)

    end

    methods (Access = protected)

        function computeCoordinates(obj)
            obj.coord_iso = obj.cutPoints.ISO;
            obj.coord_global = obj.cutPoints.GLOBAL;
        end

        function computeConnectivities(obj)
            obj.computeFacetsConnectivities();
        end

        function computeLevelSet(obj)
            obj.levelSet = obj.boundary_levelSet;
        end

    end

    methods

        function nCellNodes = get.nCellNodes(obj)
            nCellNodes = size(obj.posNodes,1);
        end

    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function settings = meshSettingsList(msh)
% Load a list of settings from the MeshSettings.mat file for a given mesh,
% updating any GUI listboxes as well.
%
% USAGE:
%	settings = meshSettingsList([msh=cur mesh of cur view]);
%
%
% ras, 07/07.
if notDefined('msh')
	msh = viewGet(getSelectedVolume, 'CurMesh');
end

% if there is no mesh, then clear the settings and return
if isempty(msh)
    settings = [];
    return
end

meshDir = fileparts(msh.path);
settingsFile = fullfile(meshDir, 'MeshSettings.mat');

if ~exist(settingsFile, 'file')
	fprintf('[%s]: No Settings File Found: %s.\n', mfilename, settingsFile);
	settings = [];
	return
end

load(settingsFile, 'settings');

%% check for GUI listboxes; update them if found
allH = findobj('Tag', 'MeshSettingsList');
for h = allH(:)'
	if isequal( get(h, 'Type'), 'uicontrol' )
		%% uicontrol listbox: set list name
		% make sure the selected entry (value) is within the range
		% of available settings
		n = get(h, 'Value');
		n = min(n, length(settings)); %#ok<*NODEF>

		set(h, 'String', {settings.name}, 'Value', n);

	elseif isequal( get(h, 'Type'), 'uimenu' )
		%% uimenu (gray menu)
		vw = getSelectedGray;

		% if several of these menus are found, find the one belonging to
		% the current view
		if length(h) > 1
			h = vw.ui.menus.meshSettings;
		end

		% if previous submenus exist, get rid of 'em
		a = get(h, 'Children');
		delete(a(1:end-4));

		% add a series of menu items to retrieve each stored setting
		for ii = 1:length(settings)
			if ii==1, sep = 'on';	else sep = 'off'; 	end

			cb = sprintf('meshRetrieveSettings(viewGet(%s, ''CurMesh''), ''%s''); ', ...
						 vw.name, settings(ii).name);

			uimenu(h, 'Label', settings(ii).name, 'Separator', sep, ...
				   'Callback', cb);
		end
	end
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function out = conwaylaws(nhood)
%CONWAYLAWS Applies Conway's genetic laws to a single pixel.
%   OUT = CONWAYLAWS(NHOOD) applies Conway's genetic laws to a single
%   pixel and its 3-by-3 neighborhood, NHOOD.
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

num_neighbors = sum(nhood(:)) - nhood(2, 2);
if nhood(2, 2) == 1
   if num_neighbors <= 1
      out = 0; % Pixel dies from isolation.
   elseif num_neighbors >= 4
      out = 0; % Pixel dies from overpopulation.
   else
      out = 1; % Pixel survives.
   end
else
   if num_neighbors == 3
      out = 1; % Birth pixel.
   else
      out = 0; % Pixel remains empty.
   end
end

author: dipum

repo: dipum-toolbox

save_path: github-repos/MATLAB/dipum-dipum-toolbox


text:
function F = binary_implies_linearequality(f,X,M,m,eps);

if nargin < 3 || isempty(M)
    [M,m,infbound] = derivebounds(f);
    if infbound
        warning('You have unbounded variables in an implication leading to a lousy big-M relaxation.');
    end
end

if all(m==0)
    lhs = m;
elseif all(m==1)
    lhs = 1-X;
else
    lhs = m.*(1-X);
end
if all(M==0)
    rhs = M;
elseif all(M==1)
    rhs = 1-X;
else
    rhs = M.*(1-X);
end
F = [f <= rhs, -f <= -lhs];

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
% Simple script to read a video file and track a single red marker
% in 2D using the hue and saturation values. The largest red blob is
% detected and the mean co-ordinate position of this blob is tracked. This
% code does not have any tracking built-in: it treats each frame
% idependently.
%
% This assumes that the video is a constant frame rate video. The resulting
% co-ordinate time series are stored in the variables: t, Xc and Yc.

fname = 'Capture 1 (5-10-2011 3-29 PM).wmv';
% Download from http://www.youtube.com/watch?v=Thnv0IQqXFE
% for a sample video on which this program was tested.

hmfr = video.MultimediaFileReader(fname);
hp  = video.VideoPlayer;
hpI = video.VideoPlayer;

hcsc = video.ColorSpaceConverter;
hcsc.Conversion = 'RGB to HSV';

fi = 0;
fps = 30;
refColor = [0.8 0.8];
cD = 0.05;
start_t = 0;
Xc = [];
Yc = [];
t  = [];
while ~isDone(hmfr)
    fi = fi + 1;
    tf = fi/fps;
    videoFrame = step(hmfr);
    if tf < start_t
        continue;
    end
    vI =  step(hcsc, videoFrame);
    [xc, yc, vI] = FrameByColorDistance(vI, refColor, cD);
    Xc = [Xc; xc]; %#ok<AGROW>
    Yc = [Yc; yc]; %#ok<AGROW>
    t  = [t; tf];  %#ok<AGROW>
    DisplayFrame(hp, videoFrame, fi, fps);
    DisplayFrame(hpI, floor(vI*255), fi, fps);
end
release(hp);
release(hmfr);

author: Sable

repo: mcbench-benchmarks


text:
%% Find a bipolar colormap, linear in grayscale, close to a template map.

%%
% Linearity in gray is important for a colormap to look reasonable when
% printed in grayscale ("black and white"), but note that precise linearity
% depends on the particular (non-unique) choice for conversion from RGB
% values to grayscale brightness/luminance, which would ideally be
% (printing-) device-dependent in this context, and using ICC color
% profiles, etc. In practice, we use a simpler conversion, like rgb2gray,
% (see code comments in colormap_optimization for further details).

%%
% We explore a range of template colormaps, but note that others could be
% better still... In particular, a colormap "dusk" has recently been added
% to real2rgb (see below), which I have not had time to investigate...
% Another area for future work would be a CIELAB-based investigation, see:
%%
% <http://www.mathworks.com/matlabcentral/fileexchange/11037-lab-color-scale>
%%
% <http://en.wikipedia.org/wiki/Lab_color_space>

%% CMRmap from Carey Rappaport
% http://www.mathworks.com/matlabcentral/fileexchange/2662
cmr = [
    0         0         0
    0.1500    0.1500    0.5000
    0.3000    0.1500    0.7500
    0.6000    0.2000    0.5000
    1.0000    0.2500    0.1500
    0.9000    0.5000    0
    0.9000    0.7500    0.1000
    0.9000    0.9000    0.5000
    1.0000    1.0000    1.0000
    ];
cmr = colormap_optimization(cmr); display(cmr)
colormap_visualization(cmr, 1)

%%
% This is nice, but I would argue it is perceptually asymmetric in that
% there is a relatively smooth transition in the positive half, while the
% negative half has a more distinct mauve band around -0.2 to -0.3.

%% Thermal from Oliver Woodford's real2rgb
% http://www.mathworks.com/matlabcentral/fileexchange/23342
therm = thermal(inf); display(therm)
therm = colormap_optimization(therm); display(therm)
colormap_visualization(therm, 1)

%%
% I think this is a really nice bipolar colormap. My only minor quibble
% is that the off-white and off-black ends are not very appealing (to me).
% Note that they are pure white and pure black in Oliver's original thermal
% scheme, but I don't like that, since I want hardcopy to be distinguished
% from any black lines or text labels and from a white background.

%% Based on recommendations in Brewer (1996), but with gray central color
% http://www.ingentaconnect.com/content/maney/caj/1996/00000033/00000002/art00002
brew1 = [
    0.2500         0    0.3333 % purple
    0.5000    0.5000    0.5000 % grey
    1.0000    1.0000    0.1667 % yellow
    ];
brew1 = colormap_optimization(brew1); display(brew1)
colormap_visualization(brew1, 1)

%%
% Although academically well-motivated, I find the single-color transitions
% away from the origin don't seem to make as good use of color to aid
% visualisation compared to some of the other maps here.

%% Based on my adaption of Brewer's (1996) recommendations
% with a neutral central colour, and two colours on each side, which I feel
% better discriminates between values within each of the two halves.
% Brewer (1996) considers issues such as colorblindness, etc. She
% recommends (separately) blue-red and purple-yellow schemes (but not
% blue-purple or red-yellow). The following seems like a reasonable
% compromise given that Brewer's table 2 shows no ideal paths through four
% colors (only a few ideal color-pairs, which sadly cannot be linked up).
brew2 = [
    0.2500         0    0.3333 % purple
    0.0000    0.2500    1.0000 % blue
    0.5000    0.5000    0.5000 % grey
    1.0000    0.2500    0.3333 % red
    1.0000    1.0000    0.1667 % yellow
    ];
brew2 = colormap_optimization(brew2); display(brew2)
brew2fine = colormap_visualization(brew2, 1);

%%
% I am probably biased, but I think the two-color transitions here are
% helpful, and there is reasonable (though not perfect) perceptual symmetry
% between the orange-yellow and blue-purple transitions (if anything, the
% purple transitions a little too sharply near the end, around -0.95). It
% could undoubtedly be improved further, but seems good enough to me...

%% Check linearity in grayscale luminance of raw and interpolated maps
colormap_visualization(brew2, 2)
colormap_visualization(brew2fine, 2)

%%
% See also: colormap_visualization, colormap_optimization, bipolar
%%
% Example:
%  map = bipolar(m, 0.5); % (gives the optimized brew2 that we get here)
%%
% Copyright 2009 Ged Ridgway at gmail com


text:
function TR=IcosahedronMesh
% Name speaks for itself.

% Get the vertex coordinates
t=(1+sqrt(5))/2; % golden ratio
x=[0 1 t];
s=[1 1 1; 1 1 -1; 1 -1 -1; 1 -1 1];
x=repmat(x,[4 1]).*s;
x=[x;circshift(x,[0 -1]);circshift(x,[0 -2])];
x_L2=sqrt(sum(x.^2,2));
x=bsxfun(@rdivide,x,x_L2);

% Triangulate the points
Tri = fliplr(convhulln(x));
TR=TriRep(Tri,x);

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
%  COVARY - For vectors, covary(X) returns the variance of X.
%             For matrices, covary(X)is a row vector containing the
%             variance of each column of X.
%
% Notes:
%   covary(X) normalizes by N-1 where N is the sequence length.
%   This makes covary(X) the best unbiased estimate of the
%   covariance if X are from a normal distribution.
%   Does not require the Matlab Signal Processing Toolbox
%
% Author: Scott Makeig, SCCN/INC/UCSD, La Jolla, 2000

% Copyright (C) 2000 Scott Makeig, SCCN/INC/UCSD, scott@sccn.ucsd.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

% 01-25-02 reformated help & license -ad

function covout = covary(data)

data = data - mean(mean(data));
if size(data,1) == 1
    data = data';   % make column vector
end
covout = sum(data.*data)/(size(data,1)-1);

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
classdef SegmentationLossPixel < dagnn.Loss
    % SegmentationLossPixel
    %
    % Similar to dagnn.SegmentationLoss, but also sets pixel weights.
    %
    % Inputs: scoresMap, labels, classWeights
    % Outputs: loss
    %
    % Note: All weights can be empty, which means they are ignored.
    % Note: If you use this for weakly supervised, the loss output will be
    % wrong (divided by instances, not images)
    %
    % Copyright by Holger Caesar, 2016

    properties (Transient)
        instanceWeights
    end

    methods
        function outputs = forward(obj, inputs, params) %#ok<INUSD>

            % Get inputs
            assert(numel(inputs) == 3);
            scoresMap = inputs{1};
            labels = inputs{2};
            classWeights = inputs{3};

            % Check inputs
            assert(~isempty(scoresMap));
            assert(~isempty(labels));

            % Compute invMass
            mass = sum(sum(labels > 0, 2), 1); % Removed the +1
            invMass = zeros(size(mass));
            nonEmpty = mass ~= 0;
            invMass(nonEmpty) = 1 ./ mass(nonEmpty);

            % Compute pixelWeights
            if isempty(classWeights)
                pixelWeights = [];
            else
                classWeightsPad = [0; classWeights(:)];

                %%% Pixel weighting
                pixelWeights = classWeightsPad(labels + 1);

                % Make sure mass of the image does not change
                curMasses = sum(sum(pixelWeights, 1), 2);
                divisor = curMasses ./ mass;
                valid = mass ~= 0 && curMasses ~= 0;
                if any(valid)
                    pixelWeights(:, :, :, valid) = bsxfun(@rdivide, pixelWeights(:, :, :, valid), divisor(valid));
                end

                % Checks
                pixelWeightsSum = sum(sum(pixelWeights, 1), 2);
                assert(all(abs(pixelWeightsSum - mass) < 1e-6) | pixelWeightsSum == 0);
            end;

            % Combine mass invMass and pixelWeights in instanceWeights
            obj.instanceWeights = invMass;
            if ~isempty(pixelWeights)
                obj.instanceWeights = bsxfun(@times, obj.instanceWeights, pixelWeights);
            end

            % Checks
            if ~isempty(obj.instanceWeights)
                assert(~any(isnan(obj.instanceWeights(:))))
            end

            % Compute loss
            loss = vl_nnloss(scoresMap, labels, [], ...
                'loss', obj.loss, ...
                'instanceWeights', obj.instanceWeights);

            assert(gather(~isnan(loss) && ~isinf(loss)));
            outputs{1} = loss;
            n = obj.numAveraged;
            m = n + size(scoresMap, 4);
            obj.average = (n * obj.average + double(gather(outputs{1}))) / m;
            obj.numAveraged = m;
        end

        function [derInputs, derParams] = backward(obj, inputs, params, derOutputs) %#ok<INUSL>

            % Get inputs
            scoresMap = inputs{1};
            labels = inputs{2};

            derInputs{1} = vl_nnloss(scoresMap, labels, derOutputs{1}, ...
                'loss', obj.loss, ...
                'instanceWeights', obj.instanceWeights);
            derInputs{2} = [];
            derInputs{3} = [];
            derParams = {};
        end

        function obj = SegmentationLossPixel(varargin)
            obj.load(varargin);
        end

        function forwardAdvanced(obj, layer)
            % Modification: Overrides standard forward pass to avoid giving up when any of
            % the inputs is empty.

            in = layer.inputIndexes;
            out = layer.outputIndexes;
            par = layer.paramIndexes;
            net = obj.net;
            inputs = {net.vars(in).value};

            % clear inputs if not needed anymore
            for v = in
                net.numPendingVarRefs(v) = net.numPendingVarRefs(v) - 1;
                if net.numPendingVarRefs(v) == 0
                    if ~net.vars(v).precious && ~net.computingDerivative && net.conserveMemory
                        net.vars(v).value = [];
                    end
                end
            end

            % call the simplified interface
            outputs = obj.forward(inputs, {net.params(par).value});
            for oi = 1:numel(out)
                net.vars(out(oi)).value = outputs{oi};
            end
        end
    end
end

author: nightrome

repo: matconvnet-calvin


text:
function [recall, precision, th, averagePrecision] = precisionRecall(confidence, testClass, col)
% Plot precision-recall curve
%    [recall, precision, th, averagePrecision] = precisionRecall(confidence, testClass, col)
%
% Input
%    confidence: score of the classifier
%    testclass:  1 = inlier, 0 = outlier
%    col: color (it will make a plot)
%
% Output:
%    recall
%    precision
%    th: threhold needed for each point in the precision-recall curve
%    area: area of the precision-recall curve
%
% Definition of precision-recall:
% Assuming that:
%    * RET is the set of all items the system has retrieved for a specific inquiry;
%    * REL is the set of relevant items for a specific inquiry;
%    * RETREL is the set of the retrieved relevant items, i.e. RETREL = RET REL.
% then precision and recall measures are obtained as follows:
%    precision = RETREL / RET
%    recall = RETREL / REL

confidence = double(confidence);

S = rand('state');
rand('state',0);
confidence = confidence + rand(size(confidence))*10^(-10);
rand('state',S)

% [th, j] = sort(confidence); th = th(:);
% %th = th(fix(linspace(1, length(th), 150))); % here the number of points is hardcoded to be 150.
% th = th(2:end-1);
%
% relevant = sum(testClass == 1);
% for t=1:length(th)
%     j = find(confidence > th(t));
%     retrieved(t) = length(j);
%     retrievedrelevant(t)  = sum(testClass(j) == 1);
% end
%
% precision = 100*retrievedrelevant ./ retrieved;
% recall    = 100*retrievedrelevant / relevant;

% Compute Precision-Recall
[S,j] = sort(-confidence); % retrieved elements are the ones above or equal to the threshold
C = testClass(j);
n = length(C);

REL    = sum(testClass);
if n>0
    RETREL = cumsum(C);
    RET    = 1:n;
else
    RETREL = 0;
    RET    = 1;
end

precision = 100*RETREL ./ RET;
recall    = 100*RETREL  / REL;
th = -S;

% compute average precision (from PASCAL source code)
ap=0;
T = linspace(0,100,101);
for t=T % why so few bins?
    p=max(precision(recall>=t));
    if isempty(p)
        p=0;
    end
    ap=ap+p/length(T);
end
averagePrecision = ap;

% Visualization
if nargin == 3
    plot(recall, precision, [col '-']); axis([0 100 0 100])

    grid on
    ylabel('Precision')
    xlabel('Recall')
    axis('square')
end

author: CSAILVision

repo: LabelMeToolbox


text:
classdef prtClassRvm < prtClass
    % prtClassRvm  Relevance vector machine classifier
    %
    %   CLASSIFIER = prtClassRvm returns a relevance vector machine classifier
    %
    %   CLASSIFIER = prtClassRvm(PROPERTY1, VALUE1, ...) constructs a
    %   prtClassRvm object CLASSIFIER with properties as specified by
    %   PROPERTY/VALUE pairs.
    %
    %   A prtClassRvm object inherits all properties from the abstract class
    %   prtClass. In addition is has the following properties:
    %
    %   kernels                - A cell array of prtKernel objects specifying
    %                            the kernels to use
    %   verbosePlot            - Flag indicating whether or not to plot during
    %                            training
    %   verboseText            - Flag indicating whether or not to output
    %                            verbose updates during training
    %   learningMaxIterations  - The maximum number of iterations
    %
    %   A prtClassRvm also has the following read-only properties:
    %
    %   learningConverged  - Flag indicating if the training converged
    %   beta               - The regression weights, estimated during training
    %   sparseBeta         - The sparse regression weights, estimated during
    %                        training
    %   sparseKernels      - The sparse regression kernels, estimated during
    %                        training
    %
    %   For information on relevance vector machines, please
    %   refer to the following URL:
    %
    %   http://en.wikipedia.org/wiki/Relevance_vector_machine
    %
    %   By default, prtClassRvm uses the Laplacian approximation as found
    %   in the paper:
    %
    %   Michael E. Tipping. 2001. Sparse bayesian learning and the
    %   relevance vector machine. J. Mach. Learn. Res. 1 (September 2001),
    %
    %   The code is based on the algorithm in:
    %
    %   Herbrich, Learning Kernel Classifiers, The MIT Press, 2002
    %   http://www.learning-kernel-classifiers.org/
    %
    %   A prtClassRvm object inherits the TRAIN, RUN, CROSSVALIDATE and
    %   KFOLDS methods from prtAction. It also inherits the PLOT method
    %   from prtClass.
    %
    %   Example:
    %
    %   TestDataSet = prtDataGenUnimodal;      % Create some test and
    %   TrainingDataSet = prtDataGenUnimodal;  % training data
    %   classifier = prtClassRvm;              % Create a classifier
    %   classifier = classifier.train(TrainingDataSet);    % Train
    %   classified = run(classifier, TestDataSet);         % Test
    %   % Plot the results
    %   subplot(2,1,1);
    %   classifier.plot;
    %   subplot(2,1,2);
    %   [pf,pd] = prtScoreRoc(classified,TestDataSet);
    %   h = plot(pf,pd,'linewidth',3);
    %   title('ROC'); xlabel('Pf'); ylabel('Pd');
    %
    %   % Example 2, using a different kernel
    %
    %   TestDataSet = prtDataGenUnimodal;      % Create some test and
    %   TrainingDataSet = prtDataGenUnimodal;  % training data
    %   classifier = prtClassRvm;              % Create a classifier
    %
    %   % Create a prtKernelSet object with a different pair of
    %   % prtKernels and assign them to the classifier
    %   kernSet = prtKernelDirect & prtKernelRbf;
    %   classifier.kernels = kernSet;
    %
    %   classifier = classifier.train(TrainingDataSet);    % Train
    %   classified = run(classifier, TestDataSet);         % Test
    %   % Plot
    %   subplot(2,1,1);
    %   classifier.plot;
    %   subplot(2,1,2);
    %   [pf,pd] = prtScoreRoc(classified,TestDataSet);
    %   h = plot(pf,pd,'linewidth',3);
    %   title('ROC'); xlabel('Pf'); ylabel('Pd');
    %
    %   See also prtClass, prtClassLogisticDiscriminant, prtClassBagging,
    %   prtClassMap, prtClassCap, prtClassBinaryToMaryOneVsAll, prtClassDlrt,
    %   prtClassPlsda, prtClassFld, prtClassRvmFigueiredo, prtClassRvmSequential, prtClassGlrt,  prtClass

    properties (SetAccess=private)
        name = 'Relevance Vector Machine'  % Relevance Vector Machine
        nameAbbreviation = 'RVM'           % RVM
        isNativeMary = false;  % False
    end

    properties
        kernels = prtKernelDc & prtKernelRbfNdimensionScale;  % The kernels to be used

        verboseText = false;  % Whether or not to display text during training
        verbosePlot = false;  % Whether or not to plot during training

        learningMaxIterations = 1000;       % The maximum number of iterations
        learningConvergedTolerance = 1e-5;  % Learning tolerance;
        % at iteration i, if ||if \theta_{i}-\theta_{i-1}|| / length(theta)
        % < learningConvergedTolerance, learning has converged

        learningRelevantTolerance = 1e-5;   %Tolerance below which a kernel is marked as irrelevant and removed
        % Tolerance on \theta; if \theta is < learningConvergedTolerance, the kernel is irrelevant and can be ignored
    end

    % Estimated Parameters
    properties (GetAccess = public, SetAccess = protected)
        beta = [];    % Regression weights
        sparseBeta = [];  % Sparse Beta
        sparseKernels = {};  % Sparse Kernel array
        learningConverged = false;   % Flag indicating whether or not training convereged
    end

    properties
    end

    methods

        function Obj = prtClassRvm(varargin)
            Obj = prtUtilAssignStringValuePairs(Obj,varargin{:});
        end

        function Obj = set.learningMaxIterations(Obj,val)
            if ~prtUtilIsPositiveInteger(val)
                error('prt:prtClassRvm:learningMaxIterations','learningMaxIterations must be a positive integer');
            end
            Obj.learningMaxIterations = val;
        end

        function Obj = set.learningConvergedTolerance(Obj,val)
            if ~prtUtilIsPositiveScalar(val)
                error('prt:prtClassRvm:learningConvergedTolerance','learningConvergedTolerance must be a positive scalar');
            end
            Obj.learningConvergedTolerance = val;
        end

        function Obj = set.learningRelevantTolerance(Obj,val)
            if ~prtUtilIsPositiveScalar(val)
                error('prt:prtClassRvm:learningRelevantTolerance','learningRelevantTolerance must be a positive scalar');
            end
            Obj.learningRelevantTolerance = val;
        end

        function Obj = set.kernels(Obj,val)
            assert(numel(val)==1 &&  isa(val,'prtKernel'),'prt:prtClassRvm:kernels','kernels must be a prtKernel');

            Obj.kernels = val;
        end

        function Obj = set.verbosePlot(Obj,val)
            assert(isscalar(val) && (islogical(val) || prtUtilIsPositiveInteger(val)),'prt:prtClassRvm:verbosePlot','verbosePlot must be a logical value or a positive integer');
            Obj.verbosePlot = val;
        end

        function Obj = set.verboseText(Obj,val)
            assert(isscalar(val) && islogical(val),'prt:prtClassRvm:verboseText','verboseText must be a logical value, but value provided is a %s',class(val));
            Obj.verboseText = val;
        end

        function varargout = plot(Obj)
            % plot - Plot output confidence of the prtClassRvm object
            %
            %   CLASS.plot plots the output confidence of the prtClassRvm
            %   object. The dimensionality of the dataset must be 3 or
            %   less, and verboseStorage must be true.

            HandleStructure = plot@prtClass(Obj);

            holdState = get(gca,'nextPlot');
            hold on;
            Obj.sparseKernels.plot;
            set(gca, 'nextPlot', holdState);

            varargout = {};
            if nargout > 0
                varargout = {HandleStructure};
            end
        end

    end

    methods (Access=protected, Hidden = true)

        function Obj = trainAction(Obj,DataSet)
            %Rvm = trainAction(Rvm,DataSet) (Private; see prtClass\train)
            %   Implements Jefferey's prior based training of a relevance
            %   vector machine.  The Rvm output from this function contains
            %   fields "sparseBeta" and "sparseKernels"
            %

            assert(DataSet.isBinary,'prtClassRvm Requires a binary dataset');

            warningState = warning;
            warning off MATLAB:nearlySingularMatrix

            %Note: do not assume that getTargets returns a double array or
            %values "0" and "1", instead use this:
            y = Obj.getMinusOneOneTargets(DataSet);
            y(y==-1) = 0;

            localKernels = Obj.kernels.train(DataSet);
            gram = localKernels.run_OutputDoubleArray(DataSet);

            theta = ones(size(gram,2),1);
            Obj.beta = zeros(size(theta));
            deltaThetaNorm = ones(Obj.learningMaxIterations,1)*nan;

            if Obj.verboseText
                fprintf('RVM training with %d possible vectors.\n', size(gram,2));
            end

            for iteration = 1:Obj.learningMaxIterations

                %%%%
                %%See: Herbrich: Learning Kernel Classifiers, Algorithm 7, Page 328
                %%%%

                %check tolerance for basis removal
                cRelevant = theta > Obj.learningRelevantTolerance;

                Obj.beta(~cRelevant) = 0;

                cGram = gram(:,cRelevant);
                cTheta = theta(cRelevant);
                cThetaInv = diag(1./cTheta);

                if isempty(cGram)
                    error('prt:prtClassRvm:noRelevantVectors','No relevant vectors were retained; this indicates a kernel that is scaled improperly with regards to the classification problem.  Please try choosing a different kernel or modifying the kernel parameters');
                end
                [newBeta, SigmaInvChol] = prtUtilPenalizedIrls(y,cGram,Obj.beta(cRelevant),cThetaInv);

                Obj.beta(cRelevant) = newBeta;

                SigmaChol = inv(SigmaInvChol);
                sigma = SigmaChol*SigmaChol'; %#ok<MINV>

                zeta = ones(size(diag(cThetaInv))) - (1./cTheta).*diag(sigma);

                previousTheta= theta;
                theta(cRelevant) = Obj.beta(cRelevant).^2./zeta;

                deltaThetaNorm(iteration) = norm(previousTheta-theta)./length(theta);

                if ~mod(iteration,Obj.verbosePlot)
                    if DataSet.nFeatures == 2
                        Obj.verboseIterationPlot(DataSet,cRelevant);
                    elseif iteration == 1
                        warning('prt:prtClassRvm','Learning iteration plot can only be produced for training Datasets with 2 features');
                    end
                end

                if deltaThetaNorm(iteration) < Obj.learningConvergedTolerance && iteration > 1
                    % Converged

                    Obj.learningConverged = true;

                    if Obj.verboseText
                        fprintf('Convergence reached. Exiting...\n\n');
                    end

                    break;
                end

                if Obj.verboseText
                    fprintf('\t Iteration %d: %d RV''s, Convergence tolerance: %g \n',iteration, sum(cRelevant), deltaThetaNorm(iteration));
                end

            end

            if Obj.verboseText && iteration == Obj.learningMaxIterations
                fprintf('Exiting...Convergence not reached before the maximum allowed iterations was reached.\n\n');
            end

            % Make sparse represenation
            Obj.sparseBeta = Obj.beta(cRelevant,1);
            Obj.sparseKernels = localKernels.retainKernelDimensions(cRelevant);

            % Very bad training
            if isempty(find(cRelevant,1));
                warning('prt:prtClassRvm:NoRelevantFeatures','No relevant features were found during training.');
            end

            % Reset warning
            warning(warningState);

        end

        function DataSet = runAction(Obj,DataSet)

            if isempty(Obj.sparseBeta)
                DataSet = DataSet.setObservations(nan(DataSet.nObservations,DataSet.nFeatures));
                return
            end

            n = DataSet.nObservations;

            largestMatrixSize = prtOptionsGet('prtOptionsComputation','largestMatrixSize');

            memChunkSize = max(floor(largestMatrixSize/length(Obj.sparseBeta)),1);

            OutputMat = zeros(n,1);
            for i = 1:memChunkSize:n
                cI = i:min(i+memChunkSize,n);
                cDataSet = prtDataSetClass(DataSet.X(cI,:));

                gram = Obj.sparseKernels.run(cDataSet);

                OutputMat(cI) = prtRvUtilNormCdf(gram.getObservations*Obj.sparseBeta);
            end

            DataSet.X = OutputMat;
        end
    end

    methods (Access=protected, Hidden = true)

        function y = getMinusOneOneTargets(Obj, DataSet) %#ok<MANU>
            yMat = double(DataSet.getTargetsAsBinaryMatrix());
            y = nan(size(yMat,1),1);
            y(yMat(:,1) == 1) = -1;
            y(yMat(:,2) == 1) = 1;
        end

        function G = regularizeGramInnerProduct(Obj, gram)
            nBasis = size(gram,2);

            sigmaSquared = 1e-6;

            %Check to make sure the problem is well-posed.  This can be fixed either
            %with changes to kernels, or by regularization
            G = gram'*gram;
            while rcond(G) < 1e-6
                if sigmaSquared == eps && Obj.verboseText
                    %warning('prt:prtClassRvm:illConditionedG','RVM initial G matrix ill-conditioned; regularizing diagonal of G to resolve; this can be modified by changing kernel parameters\n');
                    fprintf('\n\tRegularizing Gram matrix...\n');
                end
                G = (sigmaSquared*eye(nBasis) + gram'*gram);
                sigmaSquared = sigmaSquared*2;
            end

        end

        function verboseIterationPlot(Obj,DataSet,relevantIndices)
            DsSummary = DataSet.summarize;

            [linGrid, gridSize,xx,yy] = prtPlotUtilGenerateGrid(DsSummary.lowerBounds, DsSummary.upperBounds, Obj.plotOptions); %#ok<ASGLU>

            localKernels = Obj.kernels.train(DataSet);
            cKernels = localKernels.retainKernelDimensions(relevantIndices);
            cPhiDataSet = cKernels.run(prtDataSetClass([xx(:),yy(:)]));
            cPhi = cPhiDataSet.getObservations;

            confMap = reshape(prtRvUtilNormCdf(cPhi*Obj.beta(relevantIndices)),gridSize);
            imagesc(xx(1,:),yy(:,1),confMap,[0,1])
            colormap(Obj.plotOptions.twoClassColorMapFunction());
            axis xy
            hold on
            plot(DataSet);
            cKernels.plot();
            hold off;

            set(gcf,'color',[1 1 1]);
            drawnow;
        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function stack = tiffread(filename, indices)
% TIFFREAD  Reads 8,16,32 bits uncompressed grayscale and (some) color tiff
% files, as well as stacks or multiple tiff images, for example those
% produced by metamorph, Zeiss LSM or NIH-image
%
% tiffread, version 2.91 Nov 1, 2010
%
% stack = tiffread;
% stack = tiffread(filename);
% stack = tiffread(filename, indices);
%
% Reads 8,16,32 bits uncompressed grayscale and (some) color tiff files,
% as well as stacks or multiple tiff images, for example those produced
% by metamorph, Zeiss LSM or NIH-image.
%
% The function can be called with a file name in the current directory,
% or without argument, in which case it pops up a file opening dialog
% to allow for a manual selection of the file.
% If the stacks contains multiples images, reading can be restricted by
% specifying the indices of the desired images (eg. 1:5), or just one index (eg. 2).
%
% The returned value 'stack' is a vector struct containing the images
% and their meta-data. The length of the vector is the number of images.
% The image pixels values are stored in a field .data, which is a simple
% matrix for gray-scale images, or a cell-array of matrices for color images.
%
% The pixels values are returned in their native (usually integer) format,
% and must be converted to be used in most matlab functions.
%
% Example:
% im = tiffread('spindle.stk');
% imshow( double(im(5).data) );
%
% Only a fraction of the TIFF standard is supported, but you may extend support
% by modifying this file. If you do so, please return your modification to us,
% such that the added functionality can be redistributed to everyone.

% ------------------------------------------------------------------------------
%
% If you would like to acknowledge tiffread2.m in a publication,
% please cite the article for which the macro was first written:
%
% Dynamic Concentration of Motors in Microtubule Arrays
% Francois Nedelec, Thomas Surrey and A.C. Maggs
% Physical Review Letters 86: 3192-3195; 2001.
% DOI:	10.1103/PhysRevLett.86.3192
%
% Thank you!
%
%
% Francois Nedelec
% nedelec -at- embl.de
% Cell Biology and Biophysics, EMBL; Meyerhofstrasse 1; 69117 Heidelberg; Germany
% http://www.embl.org
% http://www.cytosim.org
%
%
% ------------------------------------------------------------------------------
%
% Copyright (C) 1999-2010 Francois Nedelec,
% with contributions from:
%   Kendra Burbank for the waitbar
%   Hidenao Iwai for the code to read floating point images,
%   Stephen Lang to be more compliant with PlanarConfiguration
%   Jan-Ulrich Kreft for Zeiss LSM support
%   Elias Beauchanp and David Kolin for additional Metamorph support
%   Jean-Pierre Ghobril for requesting that image indices may be specified
%   Urs Utzinger for the better handling of color images, and LSM meta-data
%   O. Scott Sands for support of GeoTIFF tags
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, version 3 of the License.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details:
% <http://www.gnu.org/licenses/>.
%
% ------------------------------------------------------------------------------
%
% Ramon Casero <rcasero@gmail.com>: Minor edits

%Optimization: join adjacent TIF strips: this results in faster reads
consolidateStrips = 1;

%without argument, we ask the user to choose a file:
if nargin < 1
    [filename, pathname] = uigetfile('*.tif;*.stk;*.lsm', 'select image file');
    filename = [ pathname, filename ];
end

if (nargin<=1);  indices = 1:10000; end

% not all valid tiff tags have been included, as they are really a lot...
% if needed, tags can easily be added to this code
% See the official list of tags:
% http://partners.adobe.com/asn/developer/pdfs/tn/TIFF6.pdf
%
% the structure IMG is returned to the user, while TIF is not.
% so tags usefull to the user should be stored as fields in IMG, while
% those used only internally can be stored in TIF.

global TIF;
TIF = [];

%counters for the number of images read and skipped
img_skip  = 0;
img_read  = 1;
hWaitbar  = [];

%% set defaults values :
TIF.SampleFormat     = 1;
TIF.SamplesPerPixel  = 1;
TIF.BOS              = 'ieee-le';          %byte order string

if  isempty(findstr(filename,'.'))
    filename = [filename,'.tif'];
end

TIF.file = fopen(filename,'r','l');
if TIF.file == -1
    stkname = strrep(filename, '.tif', '.stk');
    TIF.file = fopen(stkname,'r','l');
    if TIF.file == -1
        error(['File "',filename,'" not found.']);
    else
        filename = stkname;
    end
end
[s, m] = fileattrib(filename);

% obtain the full file path:
filename = m.Name;

% find the file size in bytes:
% m = dir(filename);
% filesize = m.bytes;

%% read header
% read byte order: II = little endian, MM = big endian
byte_order = fread(TIF.file, 2, '*char');
if ( strcmp(byte_order', 'II') )
    TIF.BOS = 'ieee-le';                                % Intel little-endian format
elseif ( strcmp(byte_order','MM') )
    TIF.BOS = 'ieee-be';
else
    error('This is not a TIFF file (no MM or II).');
end

%% ---- read in a number which identifies file as TIFF format
tiff_id = fread(TIF.file,1,'uint16', TIF.BOS);
if (tiff_id ~= 42)
    error('This is not a TIFF file (missing 42).');
end

%% ---- read the byte offset for the first image file directory (IFD)
TIF.img_pos = fread(TIF.file, 1, 'uint32', TIF.BOS);

while  TIF.img_pos ~= 0

    clear IMG;
    IMG.filename = filename;
    % move in the file to the first IFD
    status = fseek(TIF.file, TIF.img_pos, -1);
    if status == -1
        error('invalid file offset (error on fseek)');
    end

    %disp(strcat('reading img at pos :',num2str(TIF.img_pos)));

    %read in the number of IFD entries
    num_entries = fread(TIF.file,1,'uint16', TIF.BOS);
    %disp(strcat('num_entries =', num2str(num_entries)));

    %read and process each IFD entry
    for i = 1:num_entries

        % save the current position in the file
        file_pos  = ftell(TIF.file);

        % read entry tag
        TIF.entry_tag = fread(TIF.file, 1, 'uint16', TIF.BOS);
        % read entry
        entry = readIFDentry;
        %disp(strcat('reading entry <',num2str(TIF.entry_tag),'>'));

        switch TIF.entry_tag
            case 254
                TIF.NewSubfiletype = entry.val;
            case 256         % image width - number of column
                IMG.width          = entry.val;
            case 257         % image height - number of row
                IMG.height         = entry.val;
                TIF.ImageLength    = entry.val;
            case 258         % BitsPerSample per sample
                TIF.BitsPerSample  = entry.val;
                TIF.BytesPerSample = TIF.BitsPerSample / 8;
                IMG.bits           = TIF.BitsPerSample(1);
                %fprintf('BitsPerSample %i %i %i\n', entry.val);
            case 259         % compression
                if ( entry.val ~= 1 )
                    error(['Compression format ', num2str(entry.val),' not supported.']);
                end
            case 262         % photometric interpretation
                TIF.PhotometricInterpretation = entry.val;
                if ( TIF.PhotometricInterpretation == 3 )
                    warning('tiffread2:LookUp', 'Ignoring TIFF look-up table');
                end
            case 269
                IMG.document_name  = entry.val;
            case 270         % comments:
                IMG.info           = entry.val;
            case 271
                IMG.make           = entry.val;
            case 273         % strip offset
                TIF.StripOffsets   = entry.val;
                TIF.StripNumber    = entry.cnt;
                %fprintf('StripNumber = %i, size(StripOffsets) = %i %i\n', TIF.StripNumber, size(TIF.StripOffsets));
            case 277         % sample_per pixel
                TIF.SamplesPerPixel  = entry.val;
                %fprintf('Color image: sample_per_pixel=%i\n',  TIF.SamplesPerPixel);
            case 278         % rows per strip
                TIF.RowsPerStrip   = entry.val;
            case 279         % strip byte counts - number of bytes in each strip after any compressio
                TIF.StripByteCounts= entry.val;
            case 282         % X resolution
                IMG.x_resolution   = entry.val;
            case 283         % Y resolution
                IMG.y_resolution   = entry.val;
            case 284         %planar configuration describe the order of RGB
                TIF.PlanarConfiguration = entry.val;
            case 296         % resolution unit
                IMG.resolution_unit= entry.val;
            case 305         % software
                IMG.software       = entry.val;
            case 306         % datetime
                IMG.datetime       = entry.val;
            case 315
                IMG.artist         = entry.val;
            case 317        %predictor for compression
                if (entry.val ~= 1); error('unsuported predictor value'); end
            case 320         % color map
                IMG.cmap           = entry.val;
                IMG.colors         = entry.cnt/3;
            case 339
                TIF.SampleFormat   = entry.val;
            case 33550       % GeoTIFF ModelPixelScaleTag
                IMG.ModelPixelScaleTag    = entry.val;
            case 33628       %metamorph specific data
                IMG.MM_private1    = entry.val;
            case 33629       %this tag identify the image as a Metamorph stack!
                TIF.MM_stack       = entry.val;
                TIF.MM_stackCnt    = entry.cnt;
            case 33630       %metamorph stack data: wavelength
                TIF.MM_wavelength  = entry.val;
            case 33631       %metamorph stack data: gain/background?
                TIF.MM_private2    = entry.val;
            case 33922       % GeoTIFF ModelTiePointTag
                IMG.ModelTiePointTag    = entry.val;
            case 34412       % Zeiss LSM data
                LSM_info           = entry.val;
            case 34735       % GeoTIFF GeoKeyDirectory
                IMG.GeoKeyDirTag       = entry.val;
            case 34737       % GeoTIFF GeoASCIIParameters
                IMG.GeoASCII       = entry.val;
            case 42113       % GeoTIFF GDAL_NODATA
                IMG.GDAL_NODATA    = entry.val;
            otherwise
                fprintf( 'Ignored TIFF entry with tag %i (cnt %i)\n', TIF.entry_tag, entry.cnt);
        end

        % calculate bounding box  if you've got the stuff
        if isfield(IMG, 'ModelPixelScaleTag') && isfield(IMG, 'ModelTiePointTag') && isfield(IMG, 'height')&& isfield(IMG, 'width'),
            IMG.North=IMG.ModelTiePointTag(5)-IMG.ModelPixelScaleTag(2)*IMG.ModelTiePointTag(2);
            IMG.South=IMG.North-IMG.height*IMG.ModelPixelScaleTag(2);
            IMG.West=IMG.ModelTiePointTag(4)+IMG.ModelPixelScaleTag(1)*IMG.ModelTiePointTag(1);
            IMG.East=IMG.West+IMG.width*IMG.ModelPixelScaleTag(1);
        end

        % move to next IFD entry in the file
        status = fseek(TIF.file, file_pos+12, -1);
        if status == -1
            error('invalid file offset (error on fseek)');
        end
    end

    %Planar configuration is not fully supported
    %Per tiff spec 6.0 PlanarConfiguration irrelevent if SamplesPerPixel==1
    %Contributed by Stephen Lang
    if (TIF.SamplesPerPixel ~= 1) && ( ~isfield(TIF, 'PlanarConfiguration') || TIF.PlanarConfiguration == 1 )
        error('PlanarConfiguration = 1 is not supported');
    end

    %total number of bytes per image:
    PlaneBytesCnt = IMG.width * IMG.height * TIF.BytesPerSample;

    %% try to consolidate the TIFF strips if possible

    if consolidateStrips
        %Try to consolidate the strips into a single one to speed-up reading:
        BytesCnt = TIF.StripByteCounts(1);

        if BytesCnt < PlaneBytesCnt

            ConsolidateCnt = 1;
            %Count how many Strip are needed to produce a plane
            while TIF.StripOffsets(1) + BytesCnt == TIF.StripOffsets(ConsolidateCnt+1)
                ConsolidateCnt = ConsolidateCnt + 1;
                BytesCnt = BytesCnt + TIF.StripByteCounts(ConsolidateCnt);
                if ( BytesCnt >= PlaneBytesCnt ); break; end
            end

            %Consolidate the Strips
            if ( BytesCnt <= PlaneBytesCnt(1) ) && ( ConsolidateCnt > 1 )
                %fprintf('Consolidating %i stripes out of %i', ConsolidateCnt, TIF.StripNumber);
                TIF.StripByteCounts = [BytesCnt; TIF.StripByteCounts(ConsolidateCnt+1:TIF.StripNumber ) ];
                TIF.StripOffsets = TIF.StripOffsets( [1 , ConsolidateCnt+1:TIF.StripNumber] );
                TIF.StripNumber  = 1 + TIF.StripNumber - ConsolidateCnt;
            end
        end
    end

    %% read the next IFD address:
    TIF.img_pos = fread(TIF.file, 1, 'uint32', TIF.BOS);
    %if (TIF.img_pos) disp(['next ifd at', num2str(TIF.img_pos)]); end

    if isfield( TIF, 'MM_stack' )

        sel = ( indices <= TIF.MM_stackCnt );
        indices = indices(sel);

        if numel(indices) > 1
            hWaitbar = waitbar(0,'Reading images...','Name','TiffRead');
        end

        %this loop reads metamorph stacks:
        for ii = indices

            TIF.StripCnt = 1;
            offset = PlaneBytesCnt * (ii-1);

            %read the image channels
            for c = 1:TIF.SamplesPerPixel
                IMG.data{c} = read_plane(offset, IMG.width, IMG.height, c);
            end

            % print a text timer on the main window, or update the waitbar
            % fprintf('img_read %i img_skip %i\n', img_read, img_skip);
            if ~isempty( hWaitbar )
                waitbar(img_read/numel(indices), hWaitbar);
            end

            [ IMG.MM_stack, IMG.MM_wavelength, IMG.MM_private2 ] = splitMetamorph(ii);

            stack(img_read) = IMG;
            img_read = img_read + 1;

        end
        break;

    else

        %this part reads a normal TIFF stack:

        read_img = any( img_skip+img_read == indices );
        if exist('stack','var')
            if IMG.width ~= stack(1).width || IMG.height ~= stack(1).height
                %setting read_img=0 will skip dissimilar images:
                %comment-out the line below to allow dissimilar stacks
                read_img = 0;
            end
        end

        if read_img
            TIF.StripCnt = 1;
            %read the image channels
            for c = 1:TIF.SamplesPerPixel
                IMG.data{c} = read_plane(0, IMG.width, IMG.height, c);
            end

            try
                stack(img_read) = IMG;  % = orderfields(IMG);
                img_read = img_read + 1;
            catch
                fprintf('Tiffread skipped dissimilar image %i\n', img_read+img_skip);
                img_skip = img_skip + 1;
            end

            if  all( img_skip+img_read > indices )
                break;
            end

        else
            img_skip = img_skip + 1;
        end

    end
end

%% remove the cell structure if there is always only one channel
flat = 1;
for i = 1:numel(stack)
    if numel(stack(i).data) ~= 1
        flat = 0;
        break;
    end
end

if flat
    for i = 1:numel(stack)
        stack(i).data = stack(i).data{1};
    end
end

%% distribute the MetaMorph info
if isfield(TIF, 'MM_stack') && isfield(IMG, 'info') && ~isempty(IMG.info)
    MM = parseMetamorphInfo(IMG.info, TIF.MM_stackCnt);
    for i = 1:numel(stack)
        stack(i).MM = MM(i);
    end
end

%% duplicate the LSM info
if exist('LSM_info', 'var')
    for i = 1:numel(stack)
        stack(i).lsm = LSM_info;
    end
end

%% return

if ~ exist('stack', 'var')
    stack = [];
end

%clean-up
fclose(TIF.file);
if ~isempty( hWaitbar )
    delete( hWaitbar );
end

end

%% ===========================================================================

function plane = read_plane(offset, width, height, plane_nb)

global TIF;

%return an empty array if the sample format has zero bits
if ( TIF.BitsPerSample(plane_nb) == 0 )
    plane=[];
    return;
end

%fprintf('reading plane %i size %i %i\n', plane_nb, width, height);

%determine the type needed to store the pixel values:
switch( TIF.SampleFormat )
    case 1
        classname = sprintf('uint%i', TIF.BitsPerSample(plane_nb));
    case 2
        classname = sprintf('int%i', TIF.BitsPerSample(plane_nb));
    case 3
        if ( TIF.BitsPerSample(plane_nb) == 32 )
            classname = 'single';
        else
            classname = 'double';
        end
    otherwise
        error('unsuported TIFF sample format %i', TIF.SampleFormat);
end

% Preallocate a matrix to hold the sample data:
try
    plane = zeros(width, height, classname);
catch
    %compatibility with older matlab versions:
    eval(['plane = ', classname, '(zeros(width, height));']);
end

% Read the strips and concatenate them:
line = 1;
while ( TIF.StripCnt <= TIF.StripNumber )

    strip = read_strip(offset, width, plane_nb, TIF.StripCnt, classname);
    TIF.StripCnt = TIF.StripCnt + 1;

    % copy the strip onto the data
    plane(:, line:(line+size(strip,2)-1)) = strip;

    line = line + size(strip,2);
    if ( line > height )
        break;
    end

end

% Extract valid part of data if needed
if ~all(size(plane) == [width height]),
    plane = plane(1:width, 1:height);
    warning('tiffread2:Crop','Cropping data: found more bytes than needed');
end

% transpose the image (otherwise display is rotated in matlab)
plane = plane';

end

%% ================== sub-functions to read a strip ===================

function strip = read_strip(offset, width, plane_nb, stripCnt, classname)

global TIF;

%fprintf('reading strip at position %i\n',TIF.StripOffsets(stripCnt) + offset);
StripLength = TIF.StripByteCounts(stripCnt) ./ TIF.BytesPerSample(plane_nb);

%fprintf( 'reading strip %i\n', stripCnt);
status = fseek(TIF.file, TIF.StripOffsets(stripCnt) + offset, 'bof');
if status == -1
    error('invalid file offset (error on fseek)');
end

bytes = fread( TIF.file, StripLength, classname, TIF.BOS );

if any( length(bytes) ~= StripLength )
    error('End of file reached unexpectedly.');
end

strip = reshape(bytes, width, StripLength / width);

end

%% ==================sub-functions that reads an IFD entry:===================

function [nbBytes, matlabType] = convertType(tiffType)
switch (tiffType)
    case 1
        nbBytes=1;
        matlabType='uint8';
    case 2
        nbBytes=1;
        matlabType='uchar';
    case 3
        nbBytes=2;
        matlabType='uint16';
    case 4
        nbBytes=4;
        matlabType='uint32';
    case 5
        nbBytes=8;
        matlabType='uint32';
    case 7
        nbBytes=1;
        matlabType='uchar';
    case 11
        nbBytes=4;
        matlabType='float32';
    case 12
        nbBytes=8;
        matlabType='float64';
    otherwise
        error('tiff type %i not supported', tiffType)
end
end

%% ==================sub-functions that reads an IFD entry:===================

function  entry = readIFDentry()

global TIF;
entry.tiffType = fread(TIF.file, 1, 'uint16', TIF.BOS);
entry.cnt      = fread(TIF.file, 1, 'uint32', TIF.BOS);
%disp(['tiffType =', num2str(entry.tiffType),', cnt = ',num2str(entry.cnt)]);

[ entry.nbBytes, entry.matlabType ] = convertType(entry.tiffType);

if entry.nbBytes * entry.cnt > 4
    %next field contains an offset:
    offset = fread(TIF.file, 1, 'uint32', TIF.BOS);
    %disp(strcat('offset = ', num2str(offset)));
    status = fseek(TIF.file, offset, -1);
    if status == -1
        error('invalid file offset (error on fseek)');
    end

end

if TIF.entry_tag == 33629   % metamorph 'rationals'
    entry.val = fread(TIF.file, 6*entry.cnt, entry.matlabType, TIF.BOS);
elseif TIF.entry_tag == 34412  %TIF_CZ_LSMINFO
    entry.val = readLSMinfo;
else
    if entry.tiffType == 5
        entry.val = fread(TIF.file, 2*entry.cnt, entry.matlabType, TIF.BOS);
    else
        entry.val = fread(TIF.file, entry.cnt, entry.matlabType, TIF.BOS);
    end
end

if ( entry.tiffType == 2 );
    entry.val = char(entry.val');
end

end

%% =============distribute the metamorph infos to each frame:
function [MMstack, MMwavelength, MMprivate2] = splitMetamorph(imgCnt)

global TIF;

MMstack = [];
MMwavelength = [];
MMprivate2 = [];

if TIF.MM_stackCnt == 1
    return;
end

left  = imgCnt - 1;

if isfield( TIF, 'MM_stack' )
    S = length(TIF.MM_stack) / TIF.MM_stackCnt;
    MMstack = TIF.MM_stack(S*left+1:S*left+S);
end

if isfield( TIF, 'MM_wavelength' )
    S = length(TIF.MM_wavelength) / TIF.MM_stackCnt;
    MMwavelength = TIF.MM_wavelength(S*left+1:S*left+S);
end

if isfield( TIF, 'MM_private2' )
    S = length(TIF.MM_private2) / TIF.MM_stackCnt;
    MMprivate2 = TIF.MM_private2(S*left+1:S*left+S);
end

end

%% %%  Parse the Metamorph camera info tag into respective fields
% EVBR 2/7/2005, FJN Dec. 2007
function mm = parseMetamorphInfo(info, cnt)

info   = regexprep(info, '\r\n|\o0', '\n');
parse  = textscan(info, '%s %s', 'Delimiter', ':');
tokens = parse{1};
values = parse{2};

first = char(tokens(1,1));

k = 0;
mm = struct('Exposure', zeros(cnt,1));
for i=1:size(tokens,1)
    tok = char(tokens(i,1));
    val = char(values(i,1));
    %fprintf( '"%s" : "%s"\n', tok, val);
    if strcmp(tok, first)
        k = k + 1;
    end
    if strcmp(tok, 'Exposure')
        [v, c, e, pos] = sscanf(val, '%i');
        unit = val(pos:length(val));
        %return the exposure in milli-seconds
        switch( unit )
            case 'ms'
                mm(k).Exposure = v;
            case 's'
                mm(k).Exposure = v * 1000;
            otherwise
                warning('tiffread2:Unit', ['Exposure unit "',unit,'" not recognized']);
                mm(k).Exposure = v;
        end
    else
        switch tok
            case 'Binning'
                % Binning: 1 x 1 -> [1 1]
                mm(k).Binning = sscanf(val, '%d x %d')';
            case 'Region'
                mm(k).Region = sscanf(val, '%d x %d, offset at (%d, %d)')';
            otherwise
                field  = regexprep(tok, ' ', '');
                if strcmp(val, 'Off')
                    eval(['mm(k).',field,'=0;']);
                elseif strcmp(val, 'On')
                    eval(['mm(k).',field,'=1;']);
                elseif isstrprop(val,'digit')
                    eval(['mm(k).',field,'=str2num(val)'';']);
                else
                    eval(['mm(k).',field,'=val;']);
                end
        end
    end
end

end

%% ==============partial-parse of LSM info:

function R = readLSMinfo()

% Read part of the LSM info table version 2
% this provides only very partial information, since the offset indicate that
% additional data is stored in the file
global TIF;

R.MagicNumber            = sprintf('0x%09X',fread(TIF.file, 1, 'uint32', TIF.BOS));
StructureSize          = fread(TIF.file, 1, 'uint32', TIF.BOS);
R.DimensionX             = fread(TIF.file, 1, 'uint32', TIF.BOS);
R.DimensionY             = fread(TIF.file, 1, 'uint32', TIF.BOS);
R.DimensionZ             = fread(TIF.file, 1, 'uint32', TIF.BOS);
R.DimensionChannels      = fread(TIF.file, 1, 'uint32', TIF.BOS);
R.DimensionTime          = fread(TIF.file, 1, 'uint32', TIF.BOS);
R.IntensityDataType      = fread(TIF.file, 1, 'uint32', TIF.BOS);
R.ThumbnailX             = fread(TIF.file, 1, 'uint32', TIF.BOS);
R.ThumbnailY             = fread(TIF.file, 1, 'uint32', TIF.BOS);
R.VoxelSizeX             = fread(TIF.file, 1, 'float64', TIF.BOS);
R.VoxelSizeY             = fread(TIF.file, 1, 'float64', TIF.BOS);
R.VoxelSizeZ             = fread(TIF.file, 1, 'float64', TIF.BOS);
R.OriginX                = fread(TIF.file, 1, 'float64', TIF.BOS);
R.OriginY                = fread(TIF.file, 1, 'float64', TIF.BOS);
R.OriginZ                = fread(TIF.file, 1, 'float64', TIF.BOS);
R.ScanType               = fread(TIF.file, 1, 'uint16', TIF.BOS);
R.SpectralScan           = fread(TIF.file, 1, 'uint16', TIF.BOS);
R.DataType               = fread(TIF.file, 1, 'uint32', TIF.BOS);
OffsetVectorOverlay    = fread(TIF.file, 1, 'uint32', TIF.BOS);
OffsetInputLut         = fread(TIF.file, 1, 'uint32', TIF.BOS);
OffsetOutputLut        = fread(TIF.file, 1, 'uint32', TIF.BOS);
OffsetChannelColors    = fread(TIF.file, 1, 'uint32', TIF.BOS);
R.TimeInterval           = fread(TIF.file, 1, 'float64', TIF.BOS);
OffsetChannelDataTypes = fread(TIF.file, 1, 'uint32', TIF.BOS);
OffsetScanInformation  = fread(TIF.file, 1, 'uint32', TIF.BOS);
OffsetKsData           = fread(TIF.file, 1, 'uint32', TIF.BOS);
OffsetTimeStamps       = fread(TIF.file, 1, 'uint32', TIF.BOS);
OffsetEventList        = fread(TIF.file, 1, 'uint32', TIF.BOS);
OffsetRoi              = fread(TIF.file, 1, 'uint32', TIF.BOS);
OffsetBleachRoi        = fread(TIF.file, 1, 'uint32', TIF.BOS);
OffsetNextRecording    = fread(TIF.file, 1, 'uint32', TIF.BOS);

% There are more information stored in this table, which is not read here

%read real acquisition times:
if ( OffsetTimeStamps > 0 )

    status = fseek(TIF.file, OffsetTimeStamps, -1);
    if status == -1
        error('error on fseek');
    end

    StructureSize          = fread(TIF.file, 1, 'int32', TIF.BOS);
    NumberTimeStamps       = fread(TIF.file, 1, 'int32', TIF.BOS);
    for i=1:NumberTimeStamps
        R.TimeStamp(i)       = fread(TIF.file, 1, 'float64', TIF.BOS);
    end

    %calculate elapsed time from first acquisition:
    R.TimeOffset = R.TimeStamp - R.TimeStamp(1);

end

end