text:
plotAbsErr = true;

%% 500-1000

swf10 = load('swf_500_1000_10_dataset3.mat');
swf50 = load('swf_500_1000_50_dataset3.mat');
swf100 = load('swf_500_1000_100_dataset3.mat');
msckf5 = load('msckf_500_1000_min5_maxInf');
msckf10 = load('msckf_500_1000_min10_max50');
msckf20 = load('msckf_500_1000_min20_max100');
imu = load('imu_500_1000.mat');

load('../datasets/dataset3.mat');
kStart = 500;
kEnd = 1000;

transLim = [-0.5 0.5];
rotLim = [-0.25 0.25];
fontSize = 14;

if plotAbsErr
    imu.msckf_trans_err = abs(imu.msckf_trans_err);
    msckf5.msckf_trans_err = abs(msckf5.msckf_trans_err);
    msckf10.msckf_trans_err = abs(msckf10.msckf_trans_err);
    msckf20.msckf_trans_err = abs(msckf20.msckf_trans_err);
    swf10.swf_trans_err = abs(swf10.swf_trans_err);
    swf50.swf_trans_err = abs(swf50.swf_trans_err);
    swf100.swf_trans_err = abs(swf100.swf_trans_err);

    imu.msckf_rot_err = abs(imu.msckf_rot_err);
    msckf5.msckf_rot_err = abs(msckf5.msckf_rot_err);
    msckf10.msckf_rot_err = abs(msckf10.msckf_rot_err);
    msckf20.msckf_rot_err = abs(msckf20.msckf_rot_err);
    swf10.swf_rot_err = abs(swf10.swf_rot_err);
    swf50.swf_rot_err = abs(swf50.swf_rot_err);
    swf100.swf_rot_err = abs(swf100.swf_rot_err);

    transLim = transLim - transLim(1);
    rotLim = rotLim - rotLim(1);
end

figure
subplot(3,1,1)
plot(t(kStart:kEnd), imu.msckf_trans_err(1,:), '-k', 'LineWidth', 1.2)
hold on
plot(t(kStart:kEnd), msckf5.msckf_trans_err(1,:), '-b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), msckf10.msckf_trans_err(1,:), '-.b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), msckf20.msckf_trans_err(1,:), '--b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf10.swf_trans_err(1,:), '-g', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf50.swf_trans_err(1,:), '-.g', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf100.swf_trans_err(1,:), '--g', 'LineWidth', 1.2)

xlim([t(kStart) t(kEnd) ]);
ylim(transLim)
h_legend = legend('IMU Only','MSCKF 5-Inf', 'MSCKF 10-50', 'MSCKF 20-100', 'SWF 10', 'SWF 50', 'SWF 100','Location', 'northwest');
set(h_legend,'FontSize',10);

title(sprintf('Absolute Translational Error (Interval 1)'))
ylabel('\delta r_x [m]')
set(gca,'FontSize',fontSize)
grid on
grid minor
box on

subplot(3,1,2)
plot(t(kStart:kEnd), imu.msckf_trans_err(2,:), '-k', 'LineWidth', 1.2)
hold on
plot(t(kStart:kEnd), msckf5.msckf_trans_err(2,:), '-b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), msckf10.msckf_trans_err(2,:), '-.b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), msckf20.msckf_trans_err(2,:), '--b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf10.swf_trans_err(2,:), '-g', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf50.swf_trans_err(2,:), '-.g', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf100.swf_trans_err(2,:), '--g', 'LineWidth', 1.2)

xlim([t(kStart) t(kEnd) ]);
ylim(transLim)
ylabel('\delta r_y [m]')
set(gca,'FontSize',fontSize)
grid on
grid minor
box on

subplot(3,1,3)
plot(t(kStart:kEnd), imu.msckf_trans_err(3,:), '-k', 'LineWidth', 1.2)
hold on
plot(t(kStart:kEnd), msckf5.msckf_trans_err(3,:), '-b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), msckf10.msckf_trans_err(3,:), '-.b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), msckf20.msckf_trans_err(3,:), '--b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf10.swf_trans_err(3,:), '-g', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf50.swf_trans_err(3,:), '-.g', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf100.swf_trans_err(3,:), '--g', 'LineWidth', 1.2)

xlim([t(kStart) t(kEnd) ]);
ylim(transLim)
ylabel('\delta r_z [m]')
xlabel('t_k [s]')
set(gca,'FontSize',fontSize)
set(findall(gcf,'type','text'),'FontSize',fontSize)
grid on
grid minor
box on

filename = sprintf('6-Way-Comparison-500-1000-Trans.pdf');
export_fig(gcf, filename, '-transparent');

figure
subplot(3,1,1)
plot(t(kStart:kEnd), imu.msckf_rot_err(1,:), '-k', 'LineWidth', 1.2)
hold on
plot(t(kStart:kEnd), msckf5.msckf_rot_err(1,:), '-b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), msckf10.msckf_rot_err(1,:), '-.b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), msckf20.msckf_rot_err(1,:), '--b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf10.swf_rot_err(1,:), '-g', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf50.swf_rot_err(1,:), '-.g', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf100.swf_rot_err(1,:), '--g', 'LineWidth', 1.2)
xlim([t(kStart) t(kEnd) ]);
ylim(rotLim)
title(sprintf('Absolute Rotational Error (Interval 1)'))
ylabel('\delta\theta_x')
h_legend = legend('IMU Only','MSCKF 5-Inf', 'MSCKF 10-50', 'MSCKF 20-100', 'SWF 10', 'SWF 50', 'SWF 100','Location', 'northwest');
set(h_legend,'FontSize',10);
set(gca,'FontSize',fontSize)
grid on
grid minor
box on

subplot(3,1,2)
plot(t(kStart:kEnd), imu.msckf_rot_err(2,:), '-k', 'LineWidth', 1.2)
hold on
plot(t(kStart:kEnd), msckf5.msckf_rot_err(2,:), '-b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), msckf10.msckf_rot_err(2,:), '-.b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), msckf20.msckf_rot_err(2,:), '--b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf10.swf_rot_err(2,:), '-g', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf50.swf_rot_err(2,:), '-.g', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf100.swf_rot_err(2,:), '--g', 'LineWidth', 1.2)
xlim([t(kStart) t(kEnd) ]);
ylim(rotLim)
ylabel('\delta\theta_y')
set(gca,'FontSize',fontSize)
grid on
grid minor
box on

subplot(3,1,3)
plot(t(kStart:kEnd), imu.msckf_rot_err(3,:), '-k', 'LineWidth', 1.2)
hold on
plot(t(kStart:kEnd), msckf5.msckf_rot_err(3,:), '-b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), msckf10.msckf_rot_err(3,:), '-.b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), msckf20.msckf_rot_err(3,:), '--b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf10.swf_rot_err(3,:), '-g', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf50.swf_rot_err(3,:), '-.g', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf100.swf_rot_err(3,:), '--g', 'LineWidth', 1.2)
xlim([t(kStart) t(kEnd) ]);
ylim(rotLim)
ylabel('\delta\theta_z')
xlabel('t_k [s]')
grid on
grid minor
box on

set(gca,'FontSize',fontSize)
set(findall(gcf,'type','text'),'FontSize',fontSize)

filename = sprintf('6-Way-Comparison-500-1000-Rot.pdf');
export_fig(gcf, filename, '-transparent');

%% 1215-1715

swf10 = load('swf_1215_1715_10_dataset3.mat');
swf50 = load('swf_1215_1715_50_dataset3.mat');
swf100 = load('swf_1215_1715_100_dataset3.mat');
msckf5 = load('msckf_1215_1715_min5_maxInf');
msckf10 = load('msckf_1215_1715_min10_max50');
msckf20 = load('msckf_1215_1715_min20_max100');
imu = load('imu_1215_1715.mat');

plotAbsErr = true;

load('../datasets/dataset3.mat');
kStart = 1215;
kEnd = 1715;

transLim = [0 1];
rotLim = [0 1];
fontSize = 14;

if plotAbsErr
    imu.msckf_trans_err = abs(imu.msckf_trans_err);
    msckf5.msckf_trans_err = abs(msckf5.msckf_trans_err);
    msckf10.msckf_trans_err = abs(msckf10.msckf_trans_err);
    msckf20.msckf_trans_err = abs(msckf20.msckf_trans_err);
    swf10.swf_trans_err = abs(swf10.swf_trans_err);
    swf50.swf_trans_err = abs(swf50.swf_trans_err);
    swf100.swf_trans_err = abs(swf100.swf_trans_err);

    imu.msckf_rot_err = abs(imu.msckf_rot_err);
    msckf5.msckf_rot_err = abs(msckf5.msckf_rot_err);
    msckf10.msckf_rot_err = abs(msckf10.msckf_rot_err);
    msckf20.msckf_rot_err = abs(msckf20.msckf_rot_err);
    swf10.swf_rot_err = abs(swf10.swf_rot_err);
    swf50.swf_rot_err = abs(swf50.swf_rot_err);
    swf100.swf_rot_err = abs(swf100.swf_rot_err);

    transLim = transLim - transLim(1);
    rotLim = rotLim - rotLim(1);
end

figure
subplot(3,1,1)
plot(t(kStart:kEnd), imu.msckf_trans_err(1,:), '-k', 'LineWidth', 1.2)
hold on
plot(t(kStart:kEnd), msckf5.msckf_trans_err(1,:), '-b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), msckf10.msckf_trans_err(1,:), '-.b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), msckf20.msckf_trans_err(1,:), '--b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf10.swf_trans_err(1,:), '-g', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf50.swf_trans_err(1,:), '-.g', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf100.swf_trans_err(1,:), '--g', 'LineWidth', 1.2)

xlim([t(kStart) t(kEnd) ]);
ylim(transLim)
h_legend = legend('IMU Only','MSCKF 5-Inf', 'MSCKF 10-50', 'MSCKF 20-100', 'SWF 10', 'SWF 50', 'SWF 100','Location', 'northwest');
set(h_legend,'FontSize',10);

title(sprintf('Absolute Translational Error (Interval 2)'))
ylabel('\delta r_x [m]')
set(gca,'FontSize',fontSize)
grid on
grid minor
box on

subplot(3,1,2)
plot(t(kStart:kEnd), imu.msckf_trans_err(2,:), '-k', 'LineWidth', 1.2)
hold on
plot(t(kStart:kEnd), msckf5.msckf_trans_err(2,:), '-b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), msckf10.msckf_trans_err(2,:), '-.b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), msckf20.msckf_trans_err(2,:), '--b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf10.swf_trans_err(2,:), '-g', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf50.swf_trans_err(2,:), '-.g', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf100.swf_trans_err(2,:), '--g', 'LineWidth', 1.2)

xlim([t(kStart) t(kEnd) ]);
ylim(transLim)
ylabel('\delta r_y [m]')
set(gca,'FontSize',fontSize)
grid on
grid minor
box on

subplot(3,1,3)
plot(t(kStart:kEnd), imu.msckf_trans_err(3,:), '-k', 'LineWidth', 1.2)
hold on
plot(t(kStart:kEnd), msckf5.msckf_trans_err(3,:), '-b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), msckf10.msckf_trans_err(3,:), '-.b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), msckf20.msckf_trans_err(3,:), '--b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf10.swf_trans_err(3,:), '-g', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf50.swf_trans_err(3,:), '-.g', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf100.swf_trans_err(3,:), '--g', 'LineWidth', 1.2)

xlim([t(kStart) t(kEnd) ]);
ylim(transLim)
ylabel('\delta r_z [m]')
xlabel('t_k [s]')
set(gca,'FontSize',fontSize)
set(findall(gcf,'type','text'),'FontSize',fontSize)
grid on
grid minor
box on

filename = sprintf('6-Way-Comparison-1215-1715-Trans.pdf');
export_fig(gcf, filename, '-transparent');

figure
subplot(3,1,1)
plot(t(kStart:kEnd), imu.msckf_rot_err(1,:), '-k', 'LineWidth', 1.2)
hold on
plot(t(kStart:kEnd), msckf5.msckf_rot_err(1,:), '-b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), msckf10.msckf_rot_err(1,:), '-.b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), msckf20.msckf_rot_err(1,:), '--b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf10.swf_rot_err(1,:), '-g', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf50.swf_rot_err(1,:), '-.g', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf100.swf_rot_err(1,:), '--g', 'LineWidth', 1.2)
xlim([t(kStart) t(kEnd) ]);
ylim(rotLim)
title(sprintf('Absolute Rotational Error (Interval 2)'))
ylabel('\delta\theta_x')
h_legend = legend('IMU Only','MSCKF 5-Inf', 'MSCKF 10-50', 'MSCKF 20-100', 'SWF 10', 'SWF 50', 'SWF 100','Location', 'northwest');
set(h_legend,'FontSize',10);
set(gca,'FontSize',fontSize)
grid on
grid minor
box on

subplot(3,1,2)
plot(t(kStart:kEnd), imu.msckf_rot_err(2,:), '-k', 'LineWidth', 1.2)
hold on
plot(t(kStart:kEnd), msckf5.msckf_rot_err(2,:), '-b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), msckf10.msckf_rot_err(2,:), '-.b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), msckf20.msckf_rot_err(2,:), '--b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf10.swf_rot_err(2,:), '-g', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf50.swf_rot_err(2,:), '-.g', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf100.swf_rot_err(2,:), '--g', 'LineWidth', 1.2)
xlim([t(kStart) t(kEnd) ]);
ylim(rotLim)
ylabel('\delta\theta_y')
set(gca,'FontSize',fontSize)
grid on
grid minor
box on

subplot(3,1,3)
plot(t(kStart:kEnd), imu.msckf_rot_err(3,:), '-k', 'LineWidth', 1.2)
hold on
plot(t(kStart:kEnd), msckf5.msckf_rot_err(3,:), '-b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), msckf10.msckf_rot_err(3,:), '-.b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), msckf20.msckf_rot_err(3,:), '--b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf10.swf_rot_err(3,:), '-g', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf50.swf_rot_err(3,:), '-.g', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf100.swf_rot_err(3,:), '--g', 'LineWidth', 1.2)
xlim([t(kStart) t(kEnd) ]);
ylim(rotLim)
ylabel('\delta\theta_z')
xlabel('t_k [s]')
grid on
grid minor
box on

set(gca,'FontSize',fontSize)
set(findall(gcf,'type','text'),'FontSize',fontSize)

filename = sprintf('6-Way-Comparison-1215-1715-Rot.pdf');
export_fig(gcf, filename, '-transparent');

%% Less Noisy Dataset

msckf20 = load('msckf_1215_1715_min20_max100_20lessnoisy.mat');
%msckf40 = load('msckf_1215_1715_min20_max100_40lessnoisy.mat');
msckf60 = load('msckf_1215_1715_min20_max100_60lessnoisy.mat');
%msckf80 = load('msckf_1215_1715_min20_max100_60lessnoisy.mat');
msckf100 = load('msckf_1215_1715_min20_max100_100lessnoisy.mat');

imu = load('imu_1215_1715.mat');

load('../datasets/dataset3.mat');
kStart = 1215;
kEnd = 1715;
plotAbsErr = true;
transLim = [0 1];
rotLim = [0 1];
fontSize = 14;

if plotAbsErr
    imu.msckf_trans_err = abs(imu.msckf_trans_err);
    msckf20.msckf_trans_err = abs(msckf20.msckf_trans_err);
    msckf100.msckf_trans_err = abs(msckf100.msckf_trans_err);
    msckf60.msckf_trans_err = abs(msckf60.msckf_trans_err);
%     swf20.swf_trans_err = abs(swf20.swf_trans_err);
%     swf40.swf_trans_err = abs(swf40.swf_trans_err);
%     swf60.swf_trans_err = abs(swf60.swf_trans_err);

    imu.msckf_rot_err = abs(imu.msckf_rot_err);
    msckf20.msckf_rot_err = abs(msckf20.msckf_rot_err);
    msckf100.msckf_rot_err = abs(msckf100.msckf_rot_err);
    msckf60.msckf_rot_err = abs(msckf60.msckf_rot_err);
%     swf20.swf_rot_err = abs(swf20.swf_rot_err);
%     swf40.swf_rot_err = abs(swf40.swf_rot_err);
%     swf60.swf_rot_err = abs(swf60.swf_rot_err);

    transLim = transLim - transLim(1);
    rotLim = rotLim - rotLim(1);
end

figure
subplot(3,1,1)
plot(t(kStart:kEnd), imu.msckf_trans_err(1,:), '-k', 'LineWidth', 1.2)
hold on
plot(t(kStart:kEnd), msckf20.msckf_trans_err(1,:), '-b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), msckf60.msckf_trans_err(1,:), '--b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), msckf100.msckf_trans_err(1,:), '-.b', 'LineWidth', 1.2)
xlim([t(kStart) t(kEnd) ]);
ylim(transLim)
h_legend = legend('IMU Only','MSCKF 20', 'MSCKF 60', 'MSCKF 100', 'Location', 'northwest');
set(h_legend,'FontSize',10);

title(sprintf('Feature Density Comparison: Absolute Translational Error'))
ylabel('\delta r_x [m]')
set(gca,'FontSize',fontSize)
grid on
grid minor
box on

subplot(3,1,2)
plot(t(kStart:kEnd), imu.msckf_trans_err(2,:), '-k', 'LineWidth', 1.2)
hold on
plot(t(kStart:kEnd), msckf20.msckf_trans_err(2,:), '-b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), msckf60.msckf_trans_err(2,:), '--b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), msckf100.msckf_trans_err(2,:), '-.b', 'LineWidth', 1.2)

xlim([t(kStart) t(kEnd) ]);
ylim(transLim)
ylabel('\delta r_y [m]')
set(gca,'FontSize',fontSize)
grid on
grid minor
box on

subplot(3,1,3)
plot(t(kStart:kEnd), imu.msckf_trans_err(3,:), '-k', 'LineWidth', 1.2)
hold on
plot(t(kStart:kEnd), msckf20.msckf_trans_err(3,:), '-b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), msckf60.msckf_trans_err(3,:), '--b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), msckf100.msckf_trans_err(3,:), '-.b', 'LineWidth', 1.2)

xlim([t(kStart) t(kEnd) ]);
ylim(transLim)
ylabel('\delta r_z [m]')
xlabel('t_k [s]')
set(gca,'FontSize',fontSize)
set(findall(gcf,'type','text'),'FontSize',fontSize)
grid on
grid minor
box on

filename = sprintf('6-Way-Comparison-Feat-Density-1215-1715-Trans.pdf');
export_fig(gcf, filename, '-transparent');

figure
subplot(3,1,1)
plot(t(kStart:kEnd), imu.msckf_rot_err(1,:), '-k', 'LineWidth', 1.2)
hold on
plot(t(kStart:kEnd), msckf20.msckf_rot_err(1,:), '-b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), msckf60.msckf_rot_err(1,:), '--b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), msckf100.msckf_rot_err(1,:), '-.b', 'LineWidth', 1.2)
xlim([t(kStart) t(kEnd) ]);
ylim(rotLim)
title(sprintf('Feature Density Comparison: Absolute Rotational Error'))
ylabel('\delta\theta_x')
h_legend = legend('IMU Only','MSCKF 20', 'MSCKF 60', 'MSCKF 100', 'Location', 'northwest');
set(h_legend,'FontSize',10);
set(gca,'FontSize',fontSize)
grid on
grid minor
box on

subplot(3,1,2)
plot(t(kStart:kEnd), imu.msckf_rot_err(2,:), '-k', 'LineWidth', 1.2)
hold on
plot(t(kStart:kEnd), msckf20.msckf_rot_err(2,:), '-b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), msckf60.msckf_rot_err(2,:), '--b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), msckf100.msckf_rot_err(2,:), '-.b', 'LineWidth', 1.2)
xlim([t(kStart) t(kEnd)]);
ylim(rotLim)
ylabel('\delta\theta_y')
set(gca,'FontSize',fontSize)
grid on
grid minor
box on

subplot(3,1,3)
plot(t(kStart:kEnd), imu.msckf_rot_err(3,:), '-k', 'LineWidth', 1.2)
hold on
plot(t(kStart:kEnd), msckf20.msckf_rot_err(3,:), '-b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), msckf60.msckf_rot_err(3,:), '--b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), msckf100.msckf_rot_err(3,:), '-.b', 'LineWidth', 1.2)
xlim([t(kStart) t(kEnd) ]);
ylim(rotLim)
ylabel('\delta\theta_z')
xlabel('t_k [s]')
grid on
grid minor
box on

set(gca,'FontSize',fontSize)
set(findall(gcf,'type','text'),'FontSize',fontSize)

filename = sprintf('6-Way-Comparison-Feat-Density-1215-1715-Rot.pdf');
export_fig(gcf, filename, '-transparent');

%% RMSE

msckf20 = load('msckf_1215_1715_min20_max100_20lessnoisy.mat');
msckf40 = load('msckf_1215_1715_min20_max100_40lessnoisy.mat');
msckf60 = load('msckf_1215_1715_min20_max100_60lessnoisy.mat');
msckf100 = load('msckf_1215_1715_min20_max100_100lessnoisy.mat');
swf40 = load('swf_1215_1715_25_40lessnoisy.mat');
swf60 = load('swf_1215_1715_25_60lessnoisy.mat');
swf100 = load('swf_1215_1715_25_100lessnoisy.mat');

imu = load('imu_1215_1715.mat');

load('../datasets/dataset3.mat');
kStart = 1215;
kEnd = 1715;
plotAbsErr = true;
transLim = [0 0.7];
rotLim = [0 0.5];
fontSize = 14;

figure
set(gcf, 'Position', [100 100 600 200]);
plot(t(kStart:kEnd), sqrt(sum(imu.msckf_trans_err.^2, 1)/3), '-k', 'LineWidth', 1.2)
hold on
plot(t(kStart:kEnd), sqrt(sum(msckf40.msckf_trans_err.^2, 1)/3), '-b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), sqrt(sum(msckf60.msckf_trans_err.^2, 1)/3), '--b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), sqrt(sum(msckf100.msckf_trans_err.^2, 1)/3), '-.b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), sqrt(sum(swf40.swf_trans_err.^2, 1)/3), '-g', 'LineWidth', 1.2)
plot(t(kStart:kEnd), sqrt(sum(swf60.swf_trans_err.^2, 1)/3), '--g', 'LineWidth', 1.2)
plot(t(kStart:kEnd), sqrt(sum(swf100.swf_trans_err.^2, 1)/3), '-.g', 'LineWidth', 1.2)

xlim([t(kStart) t(kEnd) ]);
ylim(transLim)
h_legend = legend('IMU Only','MSCKF 40', 'MSCKF 60', 'MSCKF 100', 'SWF 40','SWF 60','SWF 100', 'Location', 'northwest');
set(h_legend,'FontSize',10);
title(sprintf('Feature Density Comparison: Position RMSE'))
xlabel('t_k [s]');
ylabel('RMSE [m]')
set(gca,'FontSize',fontSize)
grid on
grid minor
box on
filename = sprintf('RMSE-Comparison-Feat-Density-1215-1715-Trans.pdf');
export_fig(gcf, filename, '-transparent');

figure
set(gcf, 'Position', [100 100 600 200]);
plot(t(kStart:kEnd), sqrt(sum(imu.msckf_rot_err.^2, 1)/3), '-k', 'LineWidth', 1.2)
hold on
plot(t(kStart:kEnd), sqrt(sum(msckf40.msckf_rot_err.^2, 1)/3), '-b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), sqrt(sum(msckf60.msckf_rot_err.^2, 1)/3), '--b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), sqrt(sum(msckf100.msckf_rot_err.^2, 1)/3), '-.b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), sqrt(sum(swf40.swf_rot_err.^2, 1)/3), '-g', 'LineWidth', 1.2)
plot(t(kStart:kEnd), sqrt(sum(swf60.swf_rot_err.^2, 1)/3), '--g', 'LineWidth', 1.2)
plot(t(kStart:kEnd), sqrt(sum(swf100.swf_rot_err.^2, 1)/3), '-.g', 'LineWidth', 1.2)

xlim([t(kStart) t(kEnd) ]);
ylim([0 0.25])
title(sprintf('Feature Density Comparison: Rotation RMSE'))
ylabel('RMSE')
xlabel('t_k [s]');
h_legend = legend('IMU Only','MSCKF 40', 'MSCKF 60', 'MSCKF 100', 'SWF 40','SWF 60','SWF 100', 'Location', 'northwest');
set(h_legend,'FontSize',10);
set(gca,'FontSize',fontSize)
grid on
grid minor
box on
set(gca,'FontSize',fontSize)
set(findall(gcf,'type','text'),'FontSize',fontSize)

filename = sprintf('RMSE-Comparison-Feat-Density-1215-1715-Rot.pdf');
export_fig(gcf, filename, '-transparent');

%% RMSE Noisy

swf40 = load('swf_1215_1715_10_40noisy.mat');
swf60 = load('swf_1215_1715_10_60noisy.mat');
swf100 = load('swf_1215_1715_10_100noisy.mat');
imu = load('imu_1215_1715.mat');

load('../datasets/dataset3.mat');
kStart = 1215;
kEnd = 1715;
plotAbsErr = true;
transLim = [0 0.7];
rotLim = [0 0.5];
fontSize = 14;

figure
set(gcf, 'Position', [100 100 600 200]);
plot(t(kStart:kEnd), sqrt(sum(imu.msckf_trans_err.^2, 1)/3), '-k', 'LineWidth', 1.2)
hold on
plot(t(kStart:kEnd), sqrt(sum(swf40.swf_trans_err.^2, 1)/3), '-g', 'LineWidth', 1.2)
plot(t(kStart:kEnd), sqrt(sum(swf60.swf_trans_err.^2, 1)/3), '--g', 'LineWidth', 1.2)
plot(t(kStart:kEnd), sqrt(sum(swf100.swf_trans_err.^2, 1)/3), '-.g', 'LineWidth', 1.2)

xlim([t(kStart) t(kEnd) ]);
ylim(transLim)
h_legend = legend('IMU Only','MSCKF 40', 'MSCKF 60', 'MSCKF 100', 'SWF 40','SWF 60','SWF 100', 'Location', 'northwest');
set(h_legend,'FontSize',10);
title(sprintf('Feature Density Comparison: Position RMSE'))
xlabel('t_k [s]');
ylabel('RMSE [m]')
set(gca,'FontSize',fontSize)
grid on
grid minor
box on
filename = sprintf('RMSE-Comparison-Feat-Density-1215-1715-Noisy-Trans.pdf');
export_fig(gcf, filename, '-transparent');

figure
set(gcf, 'Position', [100 100 600 200]);
plot(t(kStart:kEnd), sqrt(sum(imu.msckf_rot_err.^2, 1)/3), '-k', 'LineWidth', 1.2)
hold on
plot(t(kStart:kEnd), sqrt(sum(swf40.swf_rot_err.^2, 1)/3), '-g', 'LineWidth', 1.2)
plot(t(kStart:kEnd), sqrt(sum(swf60.swf_rot_err.^2, 1)/3), '--g', 'LineWidth', 1.2)
plot(t(kStart:kEnd), sqrt(sum(swf100.swf_rot_err.^2, 1)/3), '-.g', 'LineWidth', 1.2)

xlim([t(kStart) t(kEnd) ]);
ylim([0 0.25])
title(sprintf('Feature Density Comparison: Rotation RMSE'))
ylabel('RMSE')
xlabel('t_k [s]');
h_legend = legend('IMU Only','MSCKF 40', 'MSCKF 60', 'MSCKF 100', 'SWF 40','SWF 60','SWF 100', 'Location', 'northwest');
set(h_legend,'FontSize',10);
set(gca,'FontSize',fontSize)
grid on
grid minor
box on
set(gca,'FontSize',fontSize)
set(findall(gcf,'type','text'),'FontSize',fontSize)

filename = sprintf('RMSE-Comparison-Feat-Density-1215-1715-Noisy-Rot.pdf');
export_fig(gcf, filename, '-transparent');

%% RMSE 500-1000 or 1215_1715
% Run relevent calculcateRMSE_1215_1715.m file first

load('../datasets/dataset3.mat');
kStart = 1215;
kEnd = 1715;
plotAbsErr = true;
transLim = [0 0.7];
rotLim = [0 0.5];
fontSize = 8;

figure
set(gcf, 'Position', [100 100 600 200]);
plot(t(kStart:kEnd), imu_trans, '-k', 'LineWidth', 1.2)
hold on
plot(t(kStart:kEnd), msckf5_inf_trans, '-b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), msckf10_50_trans, '--b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), msckf20_100_trans, '-.b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf10_trans, '-g', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf50_trans, '--g', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf100_trans, '-.g', 'LineWidth', 1.2)

xlim([t(kStart) t(kEnd) ]);
ylim(transLim)
h_legend = legend('IMU Only','MSCKF 5-Inf', 'MSCKF 10-20', 'MSCKF 20-100', 'SWF 10','SWF 50','SWF 100', 'Location', 'northwest');
set(h_legend,'FontSize',fontSize);
title(sprintf('Translation RMSE (Interval 2)'))
xlabel('t_k [s]');
ylabel('RMSE [m]')
set(gca,'FontSize',fontSize)
grid on
grid minor
box on
filename = sprintf('RMSE-Comparison-1215-1715-Trans.pdf');
export_fig(gcf, filename, '-transparent');

figure
set(gcf, 'Position', [100 100 600 200]);
plot(t(kStart:kEnd), imu_rot, '-k', 'LineWidth', 1.2)
hold on
plot(t(kStart:kEnd), msckf5_inf_rot, '-b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), msckf10_50_rot, '--b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), msckf20_100_rot, '-.b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf10_rot, '-g', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf50_rot, '--g', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf100_rot, '-.g', 'LineWidth', 1.2)

xlim([t(kStart) t(kEnd) ]);
ylim(rotLim)
title(sprintf('Rotation RMSE (Interval 2)'))
ylabel('RMSE')
xlabel('t_k [s]');
h_legend = legend('IMU Only','MSCKF 40', 'MSCKF 60', 'MSCKF 100', 'SWF 40','SWF 60','SWF 100', 'Location', 'northwest');
set(h_legend,'FontSize',fontSize);
set(gca,'FontSize',fontSize)
grid on
grid minor
box on
set(gca,'FontSize',fontSize)
set(findall(gcf,'type','text'),'FontSize',fontSize)

filename = sprintf('RMSE-Comparison-1215-1715-Rot.pdf');
export_fig(gcf, filename, '-transparent');

%% RMSE 500-1000 or 1215_1715
% Run calculcateRMSE_500_1000.m file first

load('../datasets/dataset3.mat');
kStart = 500;
kEnd = 1000;
plotAbsErr = true;
transLim = [0 0.7];
rotLim = [0 0.5];
fontSize = 8;

figure
set(gcf, 'Position', [100 100 600 200]);
plot(t(kStart:kEnd), imu_trans, '-k', 'LineWidth', 1.2)
hold on
plot(t(kStart:kEnd), msckf5_inf_trans, '-b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), msckf10_50_trans, '--b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), msckf20_100_trans, '-.b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf10_trans, '-g', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf50_trans, '--g', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf100_trans, '-.g', 'LineWidth', 1.2)

xlim([t(kStart) t(kEnd) ]);
ylim(transLim)
h_legend = legend('IMU Only','MSCKF 5-Inf', 'MSCKF 10-20', 'MSCKF 20-100', 'SWF 10','SWF 50','SWF 100', 'Location', 'northwest');
set(h_legend,'FontSize',fontSize);
title(sprintf('Translation RMSE (Interval 1)'))
xlabel('t_k [s]');
ylabel('RMSE [m]')
set(gca,'FontSize',fontSize)
grid on
grid minor
box on
filename = sprintf('RMSE-Comparison-500-1000-Trans.pdf');
export_fig(gcf, filename, '-transparent');

figure
set(gcf, 'Position', [100 100 600 200]);
plot(t(kStart:kEnd), imu_rot, '-k', 'LineWidth', 1.2)
hold on
plot(t(kStart:kEnd), msckf5_inf_rot, '-b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), msckf10_50_rot, '--b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), msckf20_100_rot, '-.b', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf10_rot, '-g', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf50_rot, '--g', 'LineWidth', 1.2)
plot(t(kStart:kEnd), swf100_rot, '-.g', 'LineWidth', 1.2)

xlim([t(kStart) t(kEnd) ]);
ylim(rotLim)
title(sprintf('Rotation RMSE (Interval 1)'))
ylabel('RMSE')
xlabel('t_k [s]');
h_legend = legend('IMU Only','MSCKF 40', 'MSCKF 60', 'MSCKF 100', 'SWF 40','SWF 60','SWF 100', 'Location', 'northwest');
set(h_legend,'FontSize',fontSize);
set(gca,'FontSize',fontSize)
grid on
grid minor
box on
set(gca,'FontSize',fontSize)
set(findall(gcf,'type','text'),'FontSize',fontSize)

filename = sprintf('RMSE-Comparison-500-1000-Rot.pdf');
export_fig(gcf, filename, '-transparent');

repo: MSCKF


text:
function [Jul1,Jul2]=Cal2UTC(year,month,day,hour,minute,second)
%%CAL2UTC Convert dates in terms of the Gregorian calendar in years,
%         months, days, hours, minutes and seconds with the time in
%         universal coordinated time (UTC) to a two-part pseudo-Julian
%         date in UTC.
%
%INPUTS: year A matrix of integer  years in the Gregorian calendar under
%             UTC time.
%       month A matrix of integer months in the Gregorian calendar under
%             UTC time. 1<=month<=12
%         day A matrix of integer days in the Gregorian calendar under UTC
%             time. Days count from 1.
%        hour A matrix of integer hours under the Gregorian calendar. UTC
%             time 0<=hour<=23
%      minute A matrix of integer minutes in the Gregorian calendar under
%             UTC time. 0<=minute<=59.
%      second A matrix of double floating point seconds in the Gregorian
%             calendar under UTC time. This is >=0 and normally less than
%             60, but can be a value less than 61 or 59 at the right hour
%             on a day with a leap second.
%
%Since the default format of numbers in Matlab is double-precision
%floating point numbers and one must go to lengths to use values that are
%differently formatted, it is assumed that all inputs, including the
%integer values, are double-precision floats.
%
%OUTPUTS: Jul1, Jul2 Matrices of the time as pseudo-Julian dates in UTC
%                    where each row/column corresponds to the values in
%                    the same row/column of the input matrices. Jul1
%                    corresponds to an integer number of days starting at
%                    midnight, which means that it ends in 0.5, since
%                    Julian days start at noon. Jul2 is the fraction of a
%                    day (midnight-to-midnight) after that.
%
%This is a mex wrapper for the function iauDtf2d in the International
%Astronomical Union's (IAU) Standard's of Fundamental Astronomy library.
%The UTC date is only pseudo-Julian, because there is not a fixed number
%of seconds in a Julian day. The convention used in the IAU standard is
%that the Julian day matches the UTC day regardless of whether the UTC day
%is 86399, 86400 or 86401 SI seconds (depending on the presence of leap
%seconds).
%
%The algorithm can be compiled for use in Matlab  using the
%CompileCLibraries function.
%
%The algorithm is run in Matlab using the command format
%[Jul1,Jul2]=Cal2UTC(year,month,day,hour,minute,second);
%
%October 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

error('This function is only implemented as a mexed C or C++ function. Please run CompileCLibraries.m to compile the function for use.')

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function Average = LMaverage(D, objectname, HOMEIMAGES, object_size, average_size)
%
% Average = LMaverage(D, objectname, HOMEIMAGES)

% Parameters:
if nargin<4
    object_size = [256 256]; % scale normalized
end
if nargin<5
    average_size = object_size*4; % scale normalized
end

b = [ceil((average_size(1)- object_size(1))/2) ceil((average_size(2)- object_size(2))/2)];

L = 128;

disp('selecting objects to average')
if isempty(D)
    D = LMdatabase(HOMEIMAGES, HOMEIMAGES);
    for i = 1:length(D(1:10))
        img = LMimread(D, i, HOMEIMAGES);
        disp(i)
        bb = labelimage(img); % [xmin ymin width height]
        if size(bb,1) >0
            for n = 1:size(bb,1)
                min(bb(n,3), bb(n,4))

                D(i).annotation.object(n).name = objectname;
                D(i).annotation.object(n).crop = '0';

                x = [bb(n,1) bb(n,1)+bb(n,3) bb(n,1)+bb(n,3) bb(n,1)];
                y = [bb(n,2) bb(n,2) bb(n,2)+bb(n,4) bb(n,2)+bb(n,4)];

                D(i).annotation.object(n).polygon.x = x;
                D(i).annotation.object(n).polygon.y = y;

            end
        end
        close
    end
end

D = LMquery(D, 'object.name', objectname,'exact');

disp('removing small and cropped objects from the averaging')
D = addsmallobjectlabel(D, object_size(1)/4, object_size(2)/4);
D = LMquery(D, 'object.name', '-smallobject');
D = LMquery(D, 'object.occluded', 'no');

% Align all images (scale and translate) and compute averages
Average = zeros([average_size 3], 'single');

[nrows, ncols, cc] = size(Average);
Counts  = zeros([nrows ncols], 'single');

%[x,y] = meshgrid(0:ncols-1, 0:nrows-1);
%x = x - ncols/2;
%y = y - nrows/2;

figure
for n = 1:length(D)
    n
    clear img Tmp
    img = LMimread(D, n, HOMEIMAGES);

    %img = single(img);
    %img = img - min(img(:));
    %img = 256*img / max(img(:));
    %img = uint8(253*img / max(img(:))+2);

    if size(img,3)>1
        for k = 1:length(D(n).annotation.object);
            [imgCrop, ~, ~, ~, valid] = LMobjectnormalizedcrop(img, D(n).annotation, k, b, object_size(1), object_size(2));

            if valid(1,1) < 0.5 || valid(end,end)<0.5
                valid = 2*single(valid)-1;
                w = hamming(L); w= w/sum(w);
                valid2 = conv2(valid, w, 'same');
                valid2 = conv2(valid2, w', 'same');
                valid = max(0,valid2.*(valid>0));
            end

            imgCrop(isnan(imgCrop))=0;
            imgCrop = single(imgCrop);
            imgCrop(:,:,1) = imgCrop(:,:,1).*valid;
            imgCrop(:,:,2) = imgCrop(:,:,2).*valid;
            imgCrop(:,:,3) = imgCrop(:,:,3).*valid;

            Counts = Counts + valid;
            Average = Average + imgCrop;

            if 1
                Tmp = Average(1:4:end,1:4:end,:) ./ repmat(Counts(1:4:end,1:4:end)+.000000001, [1 1 3]);
                Tmp = Tmp - min(Tmp(:));
                Tmp = Tmp / max(Tmp(:))*256;

                imshow(uint8(Tmp))
                title(n)
                drawnow
            end
        end
    end
end

Average = Average ./ repmat(Counts+.00000001, [1 1 3]);
Average = Average - min(Average(:));
Average = Average / max(Average(:))*65535;

%average = average-prctile(average(:), 3);
%average = 255*average/prctile(average(:), 97);

author: CSAILVision

repo: LabelMeToolbox


text:
function  [isnr, y_hat] = Demo_IDDBM3D(experiment_number, test_image_name)
% ------------------------------------------------------------------------------------------
%
%     Demo software for BM3D-frame based image deblurring
%               Public release ver. 0.8 (beta) (June 03, 2011)
%
% ------------------------------------------------------------------------------------------
%
%  This function implements the IDDBM3D image deblurring algorithm proposed in:
%
%  [1] A.Danielyan, V. Katkovnik, and K. Egiazarian, "BM3D frames and
%   variational image deblurring," submitted to IEEE TIP, May 2011
%
% ------------------------------------------------------------------------------------------
%
% authors:               Aram Danielyan
%                        Vladimir Katkovnik
%
% web page:              http://www.cs.tut.fi/~foi/GCF-BM3D/
%
% contact:               firstname.lastname@tut.fi
%
% ------------------------------------------------------------------------------------------
% Copyright (c) 2011 Tampere University of Technology.
% All rights reserved.
% This work should be used for nonprofit purposes only.
% ------------------------------------------------------------------------------------------
%
% Disclaimer
% ----------
%
% Any unauthorized use of these routines for industrial or profit-oriented activities is
% expressively prohibited. By downloading and/or using any of these files, you implicitly
% agree to all the terms of the TUT limited license (included in the file Legal_Notice.txt).
% ------------------------------------------------------------------------------------------

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  FUNCTION INTERFACE:
%
%  [psnr, y_hat] = Demo_IDDBM3D(experiment_number, test_image_name)
%
%  INPUT:
%   1) experiment_number: 1 -> PSF 1, sigma^2 = 2
%                         2 -> PSF 1, sigma^2 = 8
%                         3 -> PSF 2, sigma^2 = 0.308
%                         4 -> PSF 3, sigma^2 = 49
%                         5 -> PSF 4, sigma^2 = 4
%                         6 -> PSF 5, sigma^2 = 64
%                         7-13 -> experiments 7-13 are not described in [1].
%                         see this file for the blur and noise parameters.
%   2) test_image_name:   a valid filename of a grayscale test image
%
%  OUTPUT:
%   1) isnr           the output improvement in SNR, dB
%   2) y_hat:         the restored image
%
%  ! The function can work without any of the input arguments,
%   in which case, the internal default ones are used !
%
%   To run this demo functions within the BM3D package should be accessible to Matlab
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
addpath('../')

if ~exist('experiment_number','var'), experiment_number=3; end
if ~exist('test_image_name','var'), test_image_name='Cameraman256.png'; end

filename=test_image_name;

if 1 %
    initType = 'bm3ddeb'; %use output of the BM3DDEB to initialize the algorithm
else
	initType = 'zeros'; %use zero image to initialize the algorithm
end

matchType = 'bm3ddeb'; %build groups using output of the BM3DDEB algorithm
numIt = 200;

fprintf('Experiment number: %d\n', experiment_number);
fprintf('Image: %s\n', filename);

%% ------- Generating bservation ---------------------------------------------
disp('--- Generating observation ----');
y=im2double(imread(filename));

[yN,xN]=size(y);

switch experiment_number
    case 1
        sigma=sqrt(2)/255;
        for x1=-7:7; for x2=-7:7; h(x1+8,x2+8)=1/(x1^2+x2^2+1); end, end; h=h./sum(h(:));
    case 2
        sigma=sqrt(8)/255;
        s1=0; for a1=-7:7; s1=s1+1; s2=0; for a2=-7:7; s2=s2+1; h(s1,s2)=1/(a1^2+a2^2+1); end, end;  h=h./sum(h(:));
    case 3
        BSNR=40;
        sigma=-1; % if "sigma=-1", then the value of sigma depends on the BSNR
        h=ones(9); h=h./sum(h(:));
    case 4
        sigma=7/255;
        h=[1 4 6 4 1]'*[1 4 6 4 1]; h=h./sum(h(:));  % PSF
    case 5
        sigma=2/255;
        h=fspecial('gaussian', 25, 1.6);
    case 6
        sigma=8/255;
        h=fspecial('gaussian', 25, .4);
    %extra experiments
    case 7
        BSNR=30;
        sigma=-1;
        h=ones(9); h=h./sum(h(:));
    case 8
        BSNR=20;
        sigma=-1;
        h=ones(9); h=h./sum(h(:));
    case 9
        BSNR=40;
        sigma=-1;
        h=fspecial('gaussian', 25, 1.6);
    case 10
        BSNR=20;
        sigma=-1;
        h=fspecial('gaussian', 25, 1.6);
    case 11
        BSNR=15;
        sigma=-1;
        h=fspecial('gaussian', 25, 1.6);
    case 12
        BSNR=40;
        sigma=-1; % if "sigma=-1", then the value of sigma depends on the BSNR
        h=ones(19); h=h./sum(h(:));
    case 13
        BSNR=25;
        sigma=-1; % if "sigma=-1", then the value of sigma depends on the BSNR
        h=ones(19); h=h./sum(h(:));
end

y_blur = imfilter(y, h, 'circular'); % performs blurring (by circular convolution)

if sigma == -1;   %% check whether to use BSNR in order to define value of sigma
    sigma=sqrt(norm(y_blur(:)-mean(y_blur(:)),2)^2 /(yN*xN*10^(BSNR/10)));
    %     Xv% compute sigma from the desired BSNR
end

%%%% Create a blurred and noisy observation
randn('seed',0);
z = y_blur + sigma*randn(yN, xN);

bsnr=10*log10(norm(y_blur(:)-mean(y_blur(:)),2)^2 /sigma^2/yN/xN);
psnr_z =PSNR(y,z,1,0);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fprintf('Observation BSNR: %4.2f, PSNR: %4.2f\n', bsnr, psnr_z);

%% ----- Computing initial estimate ---------------------
disp('--- Computing initial estimate  ----');

[dummy, y_hat_RI,y_hat_RWI,zRI] = BM3DDEB_init(experiment_number, y, z, h, sigma);

switch lower(initType)
    case 'zeros'
        y_hat_init=zeros(size(z));
    case 'zri'
        y_hat_init=zRI;
    case 'ri'
        y_hat_init=y_hat_RI;
    case 'bm3ddeb'
        y_hat_init=y_hat_RWI;

end

switch lower(matchType)
    case 'z'
        match_im = z;
    case 'y'
        match_im = y;
    case 'zri'
        match_im = zRI;
    case 'ri'
        match_im = y_hat_RI;
    case 'bm3ddeb'
        match_im = y_hat_RWI;
end

psnr_init = PSNR(y, y_hat_init,1,0);

fprintf('Initialization method: %s\n', initType);
fprintf('Initial estimate ISNR: %4.2f, PSNR: %4.2f\n', psnr_init-psnr_z, psnr_init);

%% ------- Core algorithm ---------------------
%------ Description of the parameters of the IDDBM3D function ----------
%y - true image (use [] if true image is unavaliable)
%z - observed
%h - blurring PSF
%y_hat_init - initial estimate y_0
%match_im - image used to constuct groups and calculate weights g_r
%sigma - standard deviation of the noise
%threshType = 'h'; %use 's' for soft thresholding
%numIt - number of iterations
%gamma - regularization parameter see [1]
%tau - regularization parameter see [1] (thresholding level)
%xi - regularization parameter see [1], it is always set to 1 in this implementation
%showFigure - set to True to display figure with current estimate
%--------------------------------------------------------------------

threshType = 'h';
showFigure = true;

switch threshType
    case {'s'}
        gamma_tau_xi_inits= [
            0.0004509 0.70 1;%1
            0.0006803 0.78 1;%2
            0.0003485 0.65 1;%3
            0.0005259 0.72 1;%4
            0.0005327 0.82 1;%5
            7.632e-05 0.25 1;%6
            0.0005818 0.81 1;%7
            0.001149  1.18 1;%8
            0.0004155 0.74 1;%9
            0.0005591 0.74 1;%10
            0.0007989 0.82 1;%11
            0.0006702 0.75 1;%12
            0.001931  1.83 1;%13
        ];
    case {'h'}
        gamma_tau_xi_inits= [
            0.00051   3.13 1;%1
            0.0006004 2.75 1;%2
            0.0004573 2.91 1;%3
            0.0005959 2.82 1;%4
            0.0006018 3.63 1;%5
            0.0001726 2.24 1;%6
            0.00062   2.98 1;%7
            0.001047  3.80 1;%8
            0.0005125 3.00 1;%9
            0.0005685 2.80 1;%10
            0.0005716 2.75 1;%11
            0.0005938 2.55 1;%12
            0.001602  4.16 1;%13
        ];
end

gamma = gamma_tau_xi_inits(experiment_number,1);
tau   = gamma_tau_xi_inits(experiment_number,2)/255*2.7;
xi    = gamma_tau_xi_inits(experiment_number,3);

disp('-------- Start ----------');
fprintf('Number of iterations to perform: %d\n', numIt);
fprintf('Thresholding type: %s\n', threshType);

y_hat = IDDBM3D(y, h, z, y_hat_init, match_im, sigma, threshType, numIt, gamma, tau, xi, showFigure);
psnr = PSNR(y,y_hat,1,0);
isnr = psnr-psnr_z;

disp('-------- Results --------');
fprintf('Final estimate ISNR: %4.2f, PSNR: %4.2f\n', isnr, psnr);
return;

end

function PSNRdb = PSNR(x, y, maxval, borders)
    if ~exist('borders', 'var'), borders = 0; end
    if ~exist('maxval', 'var'), maxval = 255; end

    xx=borders+1:size(x,1)-borders;
    yy=borders+1:size(x,2)-borders;

    PSNRdb = zeros(1,size(x,3));
    for fr=1:size(x,3)
        err = x(xx,yy,fr) - y(xx,yy,fr);
        PSNRdb(fr) = 10 * log10((maxval^2)/mean2(err.^2));
    end
end

author: xialeiliu

repo: RankIQA

save_path: github-repos/MATLAB/xialeiliu-RankIQA


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

function [xR, residuum, gradnorm, cost, times] = RiemannPrecondSteep( L, F, X0, Lh, Ph, opts )
% L is the actual operator
%    needs a apply(L, X) interface but can be anything
% Lh is the operator that represents the (inexact) Euclidean Hessian
%    needs a apply(L, X) interface but can be anything
% Ph is the preconditioner for Lh; should be a TTeMPS_op_laplace operator
%
% When L is Laplace+perturbation, taking Lh and Ph both the Laplacian works
% well.

t_start = tic();
% set default opts
if ~exist( 'opts', 'var');       opts = struct();     end
if ~isfield( opts, 'maxiter');   opts.maxiter = 500;  end
if ~isfield( opts, 'tol');       opts.tol = 1e-16;     end
if ~isfield( opts, 'safe_norm');       opts.safe_norm = false;     end

d = X0.order;
n = X0.size;

[xL, xR, G] = gauge_matrices( X0 );

%xL = orthogonalize(X, d);
%xR = orthogonalize(X, 1);

cost = zeros(opts.maxiter, 1);
residuum = zeros(opts.maxiter, 1);
gradnorm = zeros(opts.maxiter, 1);
times = zeros(opts.maxiter, 1);
normRHS = norm(F);

for i = 1:opts.maxiter

    g = euclid_grad( L, xR, F );

    cost(i) = cost_function_res( xR, g );
    residuum(i) = norm(g, opts.safe_norm) / normRHS;
    times(i) = toc(t_start);

    % test for stopping criterion
    if abs(residuum(i)) < opts.tol
        sprintf( 'Current residual: %g', residuum(i))
        residuum = residuum(1:i);
        cost = cost(1:i);
        times = times(1:i);
        sprintf( 'RiemannLinsolve CONVERGED after %i iterations', i )
        break
    end

    grad = TTeMPS_tangent_orth( xL, xR, g );
    gradnorm(i) = norm( grad );

    sprintf('steepest descent step %i', i)
    %P_grad = solvePrecond( L, P, grad, xL, xR, opts );
    P_grad = solvePrecond_noSaddle( Lh, Ph, grad, xL, xR, opts, G );

    %check_precond_laplace(P, grad, P_grad)
    %eta = -P_grad;

    %line search
    alpha = linesearch_linearized( L, P_grad, g )
    %alpha = linesearch_linearized2( L, P_grad, grad )
    %alpha = -1;
    xR = tangentAdd(  P_grad, alpha, true );

    [xL, G] = left_orth_with_gauge( xR );
    %xL = orthogonalize( X, d );
    %xR = orthogonalize( X, 1 );
end

end

function res = cost_function( L, X, F )
res = 0.5*innerprod( X, apply(L, X) ) - innerprod( X, F );
end
function res = cost_function_res( X, res )
res = 0.5*innerprod( X, res );
end

function res = euclid_grad( L, X, F )
res = apply(L, X) - F;
end

function alpha = linesearch_linearized( L, xi, g )
eta = tangent_to_TTeMPS( xi );
alpha = -innerprod( eta, g );
alpha = alpha / innerprod( eta, apply(L, eta) );
end

function alpha = linesearch_linearized2( L, xi, grad )
alpha = -innerprod( xi, grad );
eta = tangent_to_TTeMPS( xi );
alpha = alpha / innerprod( eta, apply(L, eta) );
end

author: NicolasBoumal

repo: manopt


text:
function [c,Ls] = unsdgtreal(f,g,a,M)
%UNSDGTREAL  Uniform non-stationary Discrete Gabor transform
%   Usage:  c=unsdgtreal(f,g,a,M);
%           [c,Ls]=unsdgtreal(f,g,a,M);
%
%   Input parameters:
%         f     : Input signal.
%         g     : Cell array of window functions.
%         a     : Vector of time positions of windows.
%         M     : Vector of numbers of frequency channels.
%   Output parameters:
%         c     : Cell array of coefficients.
%         Ls    : Length of input signal.
%
%   `unsdgtreal(f,g,a,M)` computes the non-stationary Gabor coefficients of the
%   input signal *f*. The signal *f* can be a multichannel signal, given in
%   the form of a 2D matrix of size $Ls \times W$, with *Ls* the signal
%   length and *W* the number of signal channels.
%
%   As opposed to |nsdgt| only the coefficients of the positive frequencies
%   of the output are returned. `unsdgtreal` will refuse to work for complex
%   valued input signals.
%
%   The non-stationary Gabor theory extends standard Gabor theory by
%   enabling the evolution of the window over time. It is therefore
%   necessary to specify a set of windows instead of a single window.  This
%   is done by using a cell array for *g*. In this cell array, the n'th
%   element `g{n}` is a row vector specifying the n'th window. The
%   uniformity means that the number of channels is not allowed to vary over
%   time.
%
%   The resulting coefficients is stored as a $M/2+1 \times N \times W$
%   array. `c(m,n,l)` is thus the value of the coefficient for time index *n*,
%   frequency index *m* and signal channel *l*.
%
%   The variable *a* contains the distance in samples between two
%   consecutive blocks of coefficients. The variable *M* contains the
%   number of channels for each block of coefficients. Both *a* and *M* are
%   vectors of integers.
%
%   The variables *g*, *a* and *M* must have the same length, and the result *c*
%   will also have the same length.
%
%   The time positions of the coefficients blocks can be obtained by the
%   following code. A value of 0 correspond to the first sample of the
%   signal::
%
%     timepos = cumsum(a)-a(1);
%
%   `[c,Ls]=unsdgtreal(f,g,a,M)` additionally returns the length *Ls* of the input
%   signal *f*. This is handy for reconstruction::
%
%     [c,Ls]=unsdgtreal(f,g,a,M);
%     fr=insdgtreal(c,gd,a,Ls);
%
%   will reconstruct the signal *f* no matter what the length of *f* is,
%   provided that *gd* are dual windows of *g*.
%
%   Notes:
%   ------
%
%   `unsdgtreal` uses circular border conditions, that is to say that the signal is
%   considered as periodic for windows overlapping the beginning or the
%   end of the signal.
%
%   The phaselocking convention used in `unsdgtreal` is different from the
%   convention used in the |dgt| function. `unsdgtreal` results are phaselocked (a
%   phase reference moving with the window is used), whereas |dgt| results are
%   not phaselocked (a fixed phase reference corresponding to time 0 of the
%   signal is used). See the help on |phaselock| for more details on
%   phaselocking conventions.
%
%   See also:  nsdgt, insdgtreal, nsgabdual, nsgabtight, phaselock
%
%   Demos:  demo_nsdgt
%
%   References: ltfatnote018

%   AUTHOR : Florent Jaillet
%   TESTING: TEST_NSDGTREAL
%   REFERENCE:

if ~isnumeric(a)
  error('%s: a must be numeric.',upper(mfilename));
end;

if ~isnumeric(M)
  error('%s: M must be numeric.',upper(mfilename));
end;

L=sum(a);

[f,Ls,W,wasrow,remembershape]=comp_sigreshape_pre(f,'UNSDGTREAL',0);
f=postpad(f,L);

[g,info]=nsgabwin(g,a,M);

if ~info.isuniform
    error('%s: M must be a scalar or a constant vector.',upper(mfilename));
end;
M=M(1);

timepos=cumsum(a)-a(1);

N=length(a); % Number of time positions

M2=floor(M/2)+1;
c=zeros(M2,N,W,assert_classname(f,g{1})); % Initialisation of the result

for ii=1:N
  shift=floor(length(g{ii})/2);
  temp=zeros(M,W,assert_classname(f,g{1}));

  % Windowing of the signal.
  % Possible improvements: The following could be computed faster by
  % explicitely computing the indexes instead of using modulo and the
  % repmat is not needed if the number of signal channels W=1 (but the time
  % difference when removing it whould be really small)
  temp(1:length(g{ii}))=f(mod((1:length(g{ii}))+timepos(ii)-shift-1,L)+1,:).*...
    repmat(conj(circshift(g{ii},shift)),1,W);

  temp=circshift(temp,-shift);
  if M<length(g{ii})
    % Fft size is smaller than window length, some aliasing is needed
    x=floor(length(g{ii})/M);
    y=length(g{ii})-x*M;
    % Possible improvements: the following could probably be computed
    % faster using matrix manipulation (reshape, sum...)
    temp1=temp;
    temp=zeros(M,size(temp,2),assert_classname(f,g{1}));
    for jj=0:x-1
      temp=temp+temp1(jj*M+(1:M),:);
    end
    temp(1:y,:)=temp(1:y,:)+temp1(x*M+(1:y),:);
  end

  % FFT of the windowed signal
  c(:,ii,:) = reshape(fftreal(temp),M2,1,W);
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function varargout = integral(varargin)
%INTEGRAL   Complete definite integral of CHEBFUN2.
%
%   I = INTEGRAL(F), returns the definite integral of a CHEBFUN2. Integrated
%   over its domain of definition.
%
%   I = INTEGRAL(F, g), returns the integral of a CHEBFUN2 along the curve
%   defined by the complex-valued CHEBFUN g.
%
% See also INTEGRAL2, SUM2, QUAD2D.

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

[varargout{1:nargout}] = integral@separableApprox(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% function tutorial_phantom_elekta_lcmv(tutorial_dir)
% TUTORIAL_PHANTOM_ELEKTA_LCMV: Script that runs the tests for the Elekta phantom (Mosher's LCMV validation).
%
% CORRESPONDING ONLINE TUTORIAL:
%     https://neuroimage.usc.edu/brainstorm/Tutorials/PhantomElekta
%
% INPUTS:
%     tutorial_dir: Directory where the sample_phantom.zip file has been unzipped

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Author: Ken Taylor, John C Mosher, Francois Tadel, 2016

%% ========= Parameters ===========
% The important criteria for CAPON is the data included in
% the data covariance calculation. We will test various sets of data for
% use in this calculation

DipNdx = [1:4:32]; % targeted dipoles for testing
% for smaller set used for data covariance calculations
RefDip = [5 13 17 25]; % widely separated and shallow in each quad
% RefDip = [4 12 20 28]; % widely separated and deep in each quad

DipNdx = sort(DipNdx); % other brainstorm routines end up sorting the index, so ensure

% pick one:
%SUBARRAY = 'MEG'; % all sensors
SUBARRAY = 'MEG GRAD';
%SUBARRAY = 'MEG MAG';

% Set the noise and data covariance intervals
Baseline = [-0.1 -0.01];
DataWindow = [0 .150];
sampling_rate = 100; % phantom data is collected at 1000 fs, can resample here
% epochs will be collected from beginning of baseline to end of data window

%% ===== FILES TO IMPORT =====

if ~exist('tutorial_dir','var'),
    tutorial_dir = '/Users/mosher'; % default
end

% Does the sample folder exist
if ~exist(sprintf('%s',fullfile(tutorial_dir,'sample_phantom_elekta')),'dir'),
    error(sprintf('Missing ''sample_phantom_elekta'' folder in the tutorial directory ''%s''',tutorial_dir));
end
% Build the path of the files to import
FifFile{1} = fullfile(tutorial_dir, 'sample_phantom_elekta', 'kojak_all_2000nAm_pp_no_chpi_no_ms_raw.fif');
FifFile{2} = fullfile(tutorial_dir, 'sample_phantom_elekta', 'kojak_all_200nAm_pp_no_chpi_no_ms_raw.fif');
FifFile{3} = fullfile(tutorial_dir, 'sample_phantom_elekta', 'kojak_all_20nAm_pp_no_chpi_no_ms_raw.fif');
% Select which file is processed with this script: 1=2000nAm, 2=200nAm, 3=20nAm
iFile = 2;
% Check if the folder contains the required files
if ~file_exist(FifFile{iFile})
    error(['The folder ' tutorial_dir ' does not contain the data from the file sample_phantom_elekta.zip.']);
end

%%  ===== CREATE PROTOCOL =====
% The protocol name has to be a valid folder name (no spaces, no weird characters...)
ProtocolName = sprintf('TutorialPhantomElekta_LCMV_%.0f',iFile);

% Subject name
SubjectName = 'Kojak';

% Start brainstorm without the GUI
if ~brainstorm('status')
    brainstorm nogui
end
% Delete existing protocol

iProtocol = bst_get('Protocol', ProtocolName);
% If protocol exists
YorN = 'n'; % default to delete an existing protocol

if false, % only if we care to delete
    if ~isempty(iProtocol)
        % Set as current protocol
        gui_brainstorm('SetCurrentProtocol', iProtocol)

        YorN = input(sprintf('Protocol %s already exists, delete and rebuild (y/n): ',ProtocolName),'s');

    end
end

if isempty(iProtocol) || strncmp(YorN,'y',1),
    % build protocol

    if strncmp(YorN,'y',1), % first delete it
        gui_brainstorm('DeleteProtocol', ProtocolName);
    end

    % Create new protocol
    gui_brainstorm('CreateProtocol', ProtocolName, 0, 0);

    % ===== ANATOMY =====

    % Generate the phantom anatomy
    DipTrueFile = generate_phantom_elekta(SubjectName);

else
    % or whatever the true name is from previous run
    DipTrueFile = fullfile(SubjectName,'TrueDipoles','dipoles_160816_1755.mat');

end

% Start a new report
% bst_report('Start');

%% ===== LINK CONTINUOUS FILES =====
% Process: Create link to raw files
sFilesKojak = bst_process('CallProcess', 'process_import_data_raw', [], [], ...
    'subjectname',  SubjectName, ...
    'datafile',     {FifFile{iFile}, 'FIF'}, ...
    'channelalign', 0);

%% ===== READ EVENTS =====
% Process: Read from channel
bst_process('CallProcess', 'process_evt_read', sFilesKojak, [], ...
    'stimchan',  'STI201', ...
    'trackmode', 1, ...  % Value: detect the changes of channel value
    'zero',      0);
% Process: Delete spurious other events unrelated to dipoles
bst_process('CallProcess', 'process_evt_delete', sFilesKojak, [], ...
    'eventname', '256, 768, 1792, 3840, 4096, 6144, 7168, 7680, 7936');
% Process: Rename events to have a leading zero, for proper sorting
for i = 1:9
    bst_process('CallProcess', 'process_evt_rename', sFilesKojak, [], ...
        'src',  sprintf('%.0f',i), ...
        'dest', sprintf('%02.0f',i));
end
% Delete the first event of the first category (there is always an artifact)
LinkFile = file_fullpath(sFilesKojak(1).FileName);
LinkMat = load(LinkFile, 'F');
if ~isempty(LinkMat.F.events) && ~isempty(LinkMat.F.events(1).times)
    LinkMat.F.events(1).times(1)   = [];
    LinkMat.F.events(1).epochs(1)  = [];
    if ~isempty(LinkMat.F.events(1).channels)
        LinkMat.F.events(1).channels(1)= [];
    end
    if ~isempty(LinkMat.F.events(1).notes)
        LinkMat.F.events(1).notes(1) = [];
    end
end
bst_save(LinkFile, LinkMat, 'v6', 1);

%%  ===== IMPORT EVENTS =====
% Process: Import MEG/EEG: Events
%          DC offset correction: [-100ms,-10ms]
%          Resample: 100Hz
sFilesEpochs = bst_process('CallProcess', 'process_import_data_event', sFilesKojak, [], ...
    'subjectname', SubjectName, ...
    'condition',   '', ...
    'eventname',   ['01' sprintf(',%02.0f',2:32)], ...
    'timewindow',  [], ...
    'epochtime',   [Baseline(1), DataWindow(2)], ...
    'createcond',  0, ...
    'ignoreshort', 1, ...
    'usectfcomp',  0, ...
    'usessp',      0, ...
    'freq',        sampling_rate, ...
    'baseline',    Baseline);

%%  Process: Average: By trial group (folder average)
sAvgKojak = bst_process('CallProcess', 'process_average', sFilesEpochs, [], ...
    'avgtype',    5, ...  % By trial group (folder average)
    'avg_func',   1, ...  % Arithmetic average:  mean(x)
    'weighted',   0, ...
    'keepevents', 0);

%% Process: Compute Head Model
bst_process('CallProcess', 'process_headmodel', sAvgKojak, [], ...
    'Comment',     '', ...
    'sourcespace', 2, ...  % MRI volume
    'meg',         3, ...  % Overlapping spheres
    'volumegrid',  struct(...
    'Method',        'isotropic', ...
    'nLayers',       17, ...
    'Reduction',     3, ...
    'nVerticesInit', 4000, ...
    'Resolution',    0.0025, ...
    'FileName',      []));

%% ===== SOURCE ESTIMATION =====

%% First, we use the prestim, to establish a baseline with GLS method

fprintf('\n GLS METHOD Calculating full noise prestim from all %.0f epochs.\n\n',length(sFilesEpochs))

% Process: Compute covariance (noise or data)
bst_process('CallProcess', 'process_noisecov', sFilesEpochs, [], ...
    'baseline',       Baseline, ...
    'datatimewindow', DataWindow, ...
    'sensortypes',    'MEG', ...
    'target',         1, ...  % Noise covariance     (covariance over baseline time window)
    'dcoffset',       1, ...  % Block by block, to avoid effects of slow shifts in data
    'identity',       0, ...
    'copycond',       0, ...
    'copysubj',       0, ...
    'replacefile',    1);  % Replace

% Process: Compute covariance (noise or data)
bst_process('CallProcess', 'process_noisecov', sFilesEpochs, [], ...
    'baseline',       Baseline, ...
    'datatimewindow', DataWindow, ...
    'sensortypes',    'MEG', ...
    'target',         2, ...  % Data covariance      (covariance over data time window)
    'dcoffset',       1, ...  % Block by block, to avoid effects of slow shifts in data
    'identity',       0, ...
    'copycond',       0, ...
    'copysubj',       0, ...
    'replacefile',    1);  % Replace

% and now run the source estimation and display the results

tutorial_phantom_elekta_lcmv_src(sAvgKojak,DipNdx,'gls',SUBARRAY,DipTrueFile)

%% Now the proper LCMV with same noise and data covariances
% and now run the source estimation and display the results

fprintf('\nResults using all data as covariance\n')

tutorial_phantom_elekta_lcmv_src(sAvgKojak,DipNdx,'lcmv',SUBARRAY,DipTrueFile)

%% Testing Brainstorm's LCMV vs pseudoLCMV(GLS) routine

fprintf('\n GLS USED TO GENERATE LCMV Calculating the Noise as the full data covariance from all %.0f epochs.\n\n',length(sFilesEpochs))

% Process: Compute covariance (noise or data)
bst_process('CallProcess', 'process_noisecov', sFilesEpochs, [], ...
    'baseline',       DataWindow, ...
    'datatimewindow', DataWindow, ...
    'sensortypes',    'MEG', ...
    'target',         1, ...  % Noise covariance     (covariance over baseline time window)
    'dcoffset',       1, ...  % Block by block, to avoid effects of slow shifts in data
    'identity',       0, ...
    'copycond',       0, ...
    'copysubj',       0, ...
    'replacefile',    1);  % Replace

% Process: Compute covariance (noise or data)
bst_process('CallProcess', 'process_noisecov', sFilesEpochs, [], ...
    'baseline',       DataWindow, ...
    'datatimewindow', DataWindow, ...
    'sensortypes',    'MEG', ...
    'target',         2, ...  % Data covariance      (covariance over data time window)
    'dcoffset',       1, ...  % Block by block, to avoid effects of slow shifts in data
    'identity',       0, ...
    'copycond',       0, ...
    'copysubj',       0, ...
    'replacefile',    1);  % Replace

% and now run the source estimation and display the results

fprintf('\nResults using  data as covariances and running as GLS:\n')

tutorial_phantom_elekta_lcmv_src(sAvgKojak,DipNdx,'gls',SUBARRAY,DipTrueFile)

%% Testing Brainstorm's LCMV vs pseudoLCMV(GLS) routine

fprintf('\n LCMV USED FOR GLS Calculating the Data Covariance as the  noise covariance from all %.0f epochs.\n\n',length(sFilesEpochs))

% Process: Compute covariance (noise or data)
bst_process('CallProcess', 'process_noisecov', sFilesEpochs, [], ...
    'baseline',       Baseline, ...
    'datatimewindow', Baseline, ...
    'sensortypes',    'MEG', ...
    'target',         1, ...  % Noise covariance     (covariance over baseline time window)
    'dcoffset',       1, ...  % Block by block, to avoid effects of slow shifts in data
    'identity',       0, ...
    'copycond',       0, ...
    'copysubj',       0, ...
    'replacefile',    1);  % Replace

% Process: Compute covariance (noise or data)
bst_process('CallProcess', 'process_noisecov', sFilesEpochs, [], ...
    'baseline',       Baseline, ...
    'datatimewindow', Baseline, ...
    'sensortypes',    'MEG', ...
    'target',         2, ...  % Data covariance      (covariance over data time window)
    'dcoffset',       1, ...  % Block by block, to avoid effects of slow shifts in data
    'identity',       0, ...
    'copycond',       0, ...
    'copysubj',       0, ...
    'replacefile',    1);  % Replace

% and now run the source estimation and display the results

fprintf('\nResults using prestim as covariance and running as LCMV:\n')

tutorial_phantom_elekta_lcmv_src(sAvgKojak,DipNdx,'lcmv',SUBARRAY,DipTrueFile)

%% DISABLED: Now estimate the noise covariance from just the selected reference dipoles, to establish again a baseline

if false, % if desired, not really different from full baseline
    % want just the epochs of the Reference Dipole
    EpochNdx = []; % find all of the corresponding epochs
    for i = 1:length(RefDip),
        EpochNdx = [EpochNdx find(strncmp(sprintf('%02.0f',RefDip(i)),{sFilesEpochs.Comment},2))];
    end

    % Set the baseline to the prestim

    fprintf('\nCalculating the prestim baseline from just dipole %.0f for %.0f epochs\n\n',RefDip,length(EpochNdx))

 % Process: Compute covariance (noise or data)
bst_process('CallProcess', 'process_noisecov', sFilesEpochs, [], ...
    'baseline',       Baseline, ...
    'datatimewindow', DataWindow, ...
    'sensortypes',    'MEG', ...
    'target',         1, ...  % Noise covariance     (covariance over baseline time window)
    'dcoffset',       1, ...  % Block by block, to avoid effects of slow shifts in data
    'identity',       0, ...
    'copycond',       0, ...
    'copysubj',       0, ...
    'replacefile',    1);  % Replace

% Process: Compute covariance (noise or data)
bst_process('CallProcess', 'process_noisecov', sFilesEpochs, [], ...
    'baseline',       Baseline, ...
    'datatimewindow', DataWindow, ...
    'sensortypes',    'MEG', ...
    'target',         2, ...  % Data covariance      (covariance over data time window)
    'dcoffset',       1, ...  % Block by block, to avoid effects of slow shifts in data
    'identity',       0, ...
    'copycond',       0, ...
    'copysubj',       0, ...
    'replacefile',    1);  % Replace

    fprintf('\nBaseline results a few dipole prestims as noise covariance\n')

    tutorial_phantom_elekta_lcmv_src(sAvgKojak,DipNdx,'gls',SUBARRAY,DipTrueFile)

end

%% Now estimate the data covariance from just the reference dipoles

EpochNdx = []; % find all of the corresponding epochs
for i = 1:length(RefDip),
    EpochNdx = [EpochNdx find(strncmp(sprintf('%02.0f',RefDip(i)),{sFilesEpochs.Comment},2))];
end

fprintf('\nReference Dipole set: %.0f',RefDip(1)), fprintf(', %.0f',RefDip(2:end)), fprintf('\n')
fprintf('Calculating the data covariance from %.0f epochs\n\n',length(EpochNdx))

% Process: Compute covariance (noise or data)
bst_process('CallProcess', 'process_noisecov', sFilesEpochs(EpochNdx), [], ...
    'baseline',       Baseline, ...
    'datatimewindow', DataWindow, ...
    'sensortypes',    'MEG', ...
    'target',         1, ...  % Noise covariance     (covariance over baseline time window)
    'dcoffset',       1, ...  % Block by block, to avoid effects of slow shifts in data
    'identity',       0, ...
    'copycond',       0, ...
    'copysubj',       0, ...
    'replacefile',    1);  % Replace

% Process: Compute covariance (noise or data)
bst_process('CallProcess', 'process_noisecov', sFilesEpochs(EpochNdx), [], ...
    'baseline',       Baseline, ...
    'datatimewindow', DataWindow, ...
    'sensortypes',    'MEG', ...
    'target',         2, ...  % Data covariance      (covariance over data time window)
    'dcoffset',       1, ...  % Block by block, to avoid effects of slow shifts in data
    'identity',       0, ...
    'copycond',       0, ...
    'copysubj',       0, ...
    'replacefile',    1);  % Replace

fprintf('\nLCMV results a few dipoles for their own data covariance\n\n')

tutorial_phantom_elekta_lcmv_src(sAvgKojak,DipNdx,'lcmv',SUBARRAY,DipTrueFile)

%% Now sum the reference set for a new data set, to show beamformer operation

% use sAvgKojak from above

% which files have the reference dipoles
FileNdx = [];
for ii = 1:length(RefDip),
    STR = sprintf('data_%02.0f',RefDip(ii)); % search string
    for i = 1:length(sAvgKojak),
        if ~isempty(strfind(sAvgKojak(i).FileName,STR)),
            FileNdx = [FileNdx i];
        end
    end
end

sFilesRef = sAvgKojak(FileNdx);

% % Input files
% sFiles = {...
%     'Kojak/kojak_all_200nAm_pp_no_chpi_no_ms_raw/data_05_average_160817_1419.mat', ...
%     'Kojak/kojak_all_200nAm_pp_no_chpi_no_ms_raw/data_13_average_160817_1419.mat', ...
%     'Kojak/kojak_all_200nAm_pp_no_chpi_no_ms_raw/data_17_average_160817_1419.mat', ...
%     'Kojak/kojak_all_200nAm_pp_no_chpi_no_ms_raw/data_25_average_160817_1419.mat'};

% Process: Scale each one by the total number in preparation of "mean"
sFilesRef = bst_process('CallProcess', 'process_scale', sFilesRef, [], ...
    'factor',      length(RefDip), ...
    'sensortypes', '', ...
    'overwrite',   0);

% Process: Average: By folder (subject average). Each signal maintains its
% strength
sFilesRefAvg = bst_process('CallProcess', 'process_average', sFilesRef, [], ...
    'avgtype',    3, ...  % By folder (subject average)
    'avg_func',   1, ...  % Arithmetic average:  mean(x)
    'weighted',   0, ...
    'keepevents', 0);

%%
% Save and display report
%ReportFile = bst_report('Save', []);
%bst_report('Open', ReportFile);

author: brainstorm-tools


text:
function [ ave_recall, recall,  thresholds] = compute_ave_recall_of_bbox( bbox_pred, bbox_gt, thresholds )

if ~exist('thresholds','var') || isempty(thresholds)
    thresholds = 0.5:0.05:0.95;
end

if iscell(bbox_pred)
    assert(iscell(bbox_gt));
    assert(length(bbox_pred) == length(bbox_gt));
    num_imgs = length(bbox_pred);
    overlap  = cell(num_imgs,1);
    for i = 1:num_imgs
        bbox_pred_this = bbox_pred{i};
        bbox_gt_this   = bbox_gt{i};
        num_bbox_gt    = size(bbox_gt_this,1);
        num_bbox_pred  = size(bbox_pred_this,1);
        overlap{i}      = zeros(num_bbox_gt,1,'single');

        if num_bbox_gt && num_bbox_pred
            [overlap{i}, ~] = closest_candidates(bbox_gt_this(:,1:4), bbox_pred_this(:,1:4));
            overlap{i} = single(overlap{i});
%             overlap{i} = max(boxoverlap(bbox_pred_this(:,1:4), bbox_gt_this(:,1:4)),[],1)';
        end
    end
    overlap = cell2mat(overlap);
else
    assert(size(bbox_pred,1) == size(bbox_gt,1));
    assert(size(bbox_pred,2) == size(bbox_gt,2));
    assert(size(bbox_pred,2) == 4);
    overlap = boxoverlap(bbox_pred, bbox_gt, true);
    assert(size(bbox_pred,1) == size(overlap,1));
    assert(size(overlap,2)==1);
end

[thresholds, recall, ave_recall] = compute_average_recall(overlap);

end

function [overlap, recall, AR] = compute_average_recall(unsorted_overlaps)
all_overlaps = sort(unsorted_overlaps(:)', 'ascend');
num_pos = numel(all_overlaps);
dx = 0.001;

overlap = 0:dx:1;
overlap(end) = 1;
recall = zeros(length(overlap), 1);
for i = 1:length(overlap)
recall(i) = sum(all_overlaps >= overlap(i)) / (num_pos+eps);
end

good_recall = recall(overlap >= 0.5);
AR = 2 * dx * trapz(good_recall);

if num_pos == 0
    AR = 0;
end
end

function [best_overlap,best_boxes] = closest_candidates(gt_boxes, candidates)
% do a matching between gt_boxes and candidates

  num_gt_boxes = size(gt_boxes, 1);
  num_candidates = size(candidates, 1);

%   iou_matrix = zeros(num_gt_boxes, num_candidates);
  iou_matrix = boxoverlap(candidates, gt_boxes)';
%   for i = 1:num_gt_boxes
%     iou = boxoverlap(gt_boxes(i,:), candidates);
%     iou_matrix(i,:) = iou';
%   end

  best_overlap = zeros(num_gt_boxes, 1);
  best_boxes = -ones(num_gt_boxes, 4);

  [best_overlap,best_boxes] = greedy_matching(iou_matrix, gt_boxes, candidates);
end

function [best_overlap,best_boxes] = greedy_matching(iou_matrix, gt_boxes, candidates)
  [n, m] = size(iou_matrix);
  assert(n == size(gt_boxes, 1));
  assert(m == size(candidates, 1));
  if n > m
    gt_matching = greedy_matching_rowwise(iou_matrix');
    candidate_matching = (1:m)';
  else
    gt_matching = (1:n)';
    candidate_matching = greedy_matching_rowwise(iou_matrix);
  end

  best_overlap = zeros(n, 1);
  best_boxes = zeros(n, 4);
  for pair_idx = 1:numel(gt_matching)
    gt_idx = gt_matching(pair_idx);
    candidate_idx = candidate_matching(pair_idx);

    best_overlap(gt_idx) = iou_matrix(gt_idx, candidate_idx);
    best_boxes(gt_idx,:) = candidates(candidate_idx, :);
  end
end

function [matching, objective] = greedy_matching_rowwise(iou_matrix)
  assert(size(iou_matrix, 1) <= size(iou_matrix, 2));
  n = size(iou_matrix, 1);
  matching = zeros(n, 1);
  objective = 0;
  for i = 1:n
    % find max element int matrix
    [max_per_row, max_col_per_row] = max(iou_matrix, [], 2);
    [max_iou,row] = max(max_per_row);
    if max_iou == -inf
      break
    end

    objective = objective + max_iou;
    col = max_col_per_row(row);
    matching(row) = col;
    iou_matrix(row,:) = -inf;
    iou_matrix(:,col) = -inf;
  end
end

author: gidariss

repo: mrcnn-object-detection


text:
function [ region_image ] = average_feature_region(im, region_size)
%compute cell wise averages, where a cell is a region_size*region_sized
%region in the image. Input can be uint8t, single or double matrices
%of arbitrary dimension

% region_image = zeros(floor(size(im,1)/region_size),...
%                     floor(size(im,2)/region_size),...
%                     size(im,3), 'single');

region_area = region_size.^2;

if isa(im,'double') || isa(im,'single')
    maxval = 1.0;
elseif isa(im,'unit8')
    maxval = 255;
end;

% for n = 1:size(im,3);

  %compute the integral image
  iImage = integralVecImage(im);
  %region indices
  i1 = (region_size:region_size:size(im,1)) + 1;
  i2 = (region_size:region_size:size(im,2)) + 1;
  %sum over region, divided by number of elements, and normalize to [0,1]
  %range if integer image
  region_image = (iImage(i1,i2,:,:) - iImage(i1,i2-region_size,:,:) - iImage(i1-region_size,i2,:,:) + iImage(i1-region_size,i2-region_size,:,:)) ./ (region_area * maxval);
% end;

end

author: flyers

repo: drone-tracking


text:
function [ptheta] = tapas_sem_mixed_validate_ptheta(ptheta)
%% Make sute that the input is consistent.
%
% Input
%
% Output
%

% aponteeduardo@gmail.com
% copyright (C) 2017
%

if ~isfield(ptheta, 'x')
    error('tapas:sem:mixed:ptheta', 'Design matrix ptheta.m is not specified');
end

if ~isfield(ptheta, 'jm')
    error('tapas:sem:mixed:ptheta', ...
        'Reduction matrix ptheta.jm is not specified');
end

if ~isfield(ptheta, 'mixed')
    error('tapas:sem:mixed:ptheta', ...
        'Random effects ptheta.mixed is not specified');
end

% Now make sure that the every thing is consistent

jm1, jm2 = size(ptheta.jm);
x1, x2 = size(ptheta.x);
m1, m2 = size(ptheta.mixed);

test_matrix(ptheta.jm, jm1, jm2, 'jm');
test_matrix(ptheta.x, jm2, x2, 'x');
test_matrix(ptheta.mixed, jm2, m2, 'mixed');

ptheta.x = logical(ptheta.x);

end

function test_matrix(amatrix, d1, d2, name)
% Make standard test on a matrix

if ~isnumeric(amatrix)
    error(sprintf('tapas:sem:mixed:ptheta:%s', name), ...
        '%s not numeric', name);

end

if ~all(size(amatrix) == [d1, d2])
    error(sprintf('tapas:sem:mixed:ptheta:%s', name), ...
        '%s dimension are %d, %d, instead %d, %d', ...
        name, d1, d2, size(amatrix, 1), size(amatrix, 2));
end

end

author: translationalneuromodeling

repo: tapas


text:
% std_ersp() - Compute ERSP and/or ITC transforms for ICA components
%              or data channels of a dataset. Save results into Matlab
%              float files. When these output files already exist, loads
%              the ERSP/ITC information from them unless the requested
%              flag specifies differently. If so, a query window
%              pops up.
%
% Function description:
%              The function returns the masked (as per the requested alpha)
%              mean ERSP or ITC for the selected dataset ICA components or
%              data channels in the requested frequency range and time window
%              (the two are dependent). Frequencies are equally log spaced.
%              Options specify component numbers, desired frequency range,
%              time window length, frequency resolution, significance level,
%              and wavelet cycles. See >> help newtimef and >> timef details
%
%              Two Matlab files are saved (for ERSP and ITC). These contain
%              the ERSP|ITC image, plus the transform parameters
%              used to compute them. Saves the computed dataset mean images
%              in dataset-name files with extensions '.icaersp' and '.icaitc'
%              for ICA components or '.datersp', '.datitc' for data channels.
%              If the ERSPs/ITCs were previously saved into these files,
%              and the same set of ERSP/ITC parameters are used, the values
%              are not recomputed, but the information is read from these
%              files. Vectors of frequencies and latencies for the ERSP/ITC
%              images are returned separately. Returned 'EEG.etc' fields
%              are modified with pointers to the output float files and some
%              information about them.
% Usage:
%              >> [X times logfreqs ] = std_ersp(EEG, 'key', 'val', ...);
% Inputs:
%   EEG          - an EEG dataset structure.
%
% Optional inputs:
%   'components' - [numeric vector] components in the EEG structure for which
%                  ERSP and ITC data will be computed {default|[]: all
%                  components if no 'channels' are specified (see below)}
%   'channels'   - [numeric vector or cell array of channel labels] channels
%                  in the EEG structure for which ERSP and ITC will be computed
%                  {default|[]: no channels}
%   'freqs'      - [minHz maxHz] the ERSP/ITC frequency range to compute
%                  and return. {default: 3 to EEG sampling rate divided by 3}
%   'timelimits' - [minms maxms] time window (in ms) to compute.
%                  {default: whole input epoch}.
%   'timewindow' - [minms maxms] time window (in ms) to plot.
%                  {default: all output latencies}
%   'cycles'     - [wavecycles (factor)]. If 0 -> DFT (constant window length
%                  across frequencies).
%                  If >0 -> the number of cycles in each analysis wavelet.
%                  If [wavecycles factor], wavelet cycles increase with
%                  frequency, beginning at wavecyles. (0 < factor < 1)
%                  factor = 0 -> fixed epoch length (DFT, as in FFT).
%                  factor = 1 -> no increase (standard wavelets)
%                  {default: [0]}
%   'padratio'   - (power of 2). Multiply the number of output frequencies
%                  by dividing their frequency spacing through 0-padding.
%                  Output frequency spacing is (low_freq/padratio).
%   'alpha'      - If in (0, 1), compute two-tailed permutation-based
%                  probability thresholds and use these to mask the output
%                  ERSP/ITC images {default: NaN}
%   'type'       - ['ersp'|'itc'] though both ERSP and ITC images are computed
%                  and saved to disk, only this transform is returned to the
%                  command line (see first output, X, below) {default: 'ersp'}
%   'savetrials' - ['on'|'off'] Save single-trial time-freq. decompositions in
%                  a file with extension '.dattimef' (channels) or '.icatimef'
%                  (components). {default: 'off'}
%   'powbase'    - [ncomps,nfreqs] optional input matrix giving baseline power
%                  spectra (not dB power, see >> help timef).
%                  For use in repeated calls to timef() using the same baseine
%                  {default|[] -> none; data windows centered before 0 latency}
%   'recompute'  - ['on'|'off'] 'on' forces recomputation of both ERSP and ITC.
%                  {default: 'off'}
%
% Other optional inputs:
%   This function will take any of the newtimef() optional inputs (for instance
%   to compute log-space frequencies)...
%
% Outputs:
%   X         - the masked log ERSP/ITC of the requested ICA components/channels
%               in the selected frequency and time range.
%   times     - vector of time points for which the ERSPs/ITCs were computed.
%   logfreqs  - vector of (equally log spaced) frequencies (in Hz) at which the
%               log ERSP/ITC was evaluated.
%
% Files written or modified:
%              [dataset_filename].icaersp   <-- saved component ERSPs
%              [dataset_filename].icaitc    <-- saved component ITCs
%  OR for channels
%              [dataset_filename].datersp   <-- saved channel ERSPs
%              [dataset_filename].datitc    <-- saved channel ITCs
% Example:
%            % Create mean ERSP and ITC images on disk for all comps from
%            % dataset EEG use three-cycle wavelets (at 3 Hz) to more than
%            % three-cycle wavelets at 50 Hz. Use probability masking at
%            % p < 0.01 with padratio 4. See >> timef details.
%            % Return the (equally log-freq spaced, probability-masked) ERSP.
%            >> [Xersp, times, logfreqs] = std_ersp(EEG, ...
%                       'type', 'ersp', 'freqs', [3 50], ...
%                                 'cycles', [3 0.5], 'alpha', 0.01);
%
% See also: timef(), std_itc(), std_erp(), std_spec(), std_topo(), std_preclust()
%
% Authors: Arnaud Delorme, Hilit Serby, SCCN, INC, UCSD, January, 2005-

% Copyright (C) Hilit Serby, SCCN, INC, UCSD, October 11, 2004, hilit@sccn.ucsd.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function [X, times, freqs, parameters] = std_ersp(EEG, varargin)

if nargin < 1
    help std_ersp;
    return;
end;

options = {};
if length(varargin) > 1
    if ~isstr(varargin{1})
        if length(varargin) > 0, options = { options{:} 'components' varargin{1} }; end;
        if length(varargin) > 1, options = { options{:} 'freqs'      varargin{2} }; end;
        if length(varargin) > 2, options = { options{:} 'timewindow' varargin{3} }; end;
        if length(varargin) > 3, options = { options{:} 'cycles'     varargin{4} }; end;
        if length(varargin) > 4, options = { options{:} 'padratio'   varargin{5} }; end;
        if length(varargin) > 5, options = { options{:} 'alpha'      varargin{6} }; end;
        if length(varargin) > 6, options = { options{:} 'type'       varargin{7} }; end;
        if length(varargin) > 7, options = { options{:} 'powbase'    varargin{8} }; end;
    else
        options = varargin;
    end;
end;

[g timefargs] = finputcheck(options, { ...
                        'components'    'integer'     []      [];
                        'channels'      { 'cell' 'integer' }  { [] [] }     {};
                        'outputfile'    'string'      []      '';
                        'powbase'       'real'        []      [];
                        'trialindices' { 'integer' 'cell' } []         [];
                        'savetrials'    'string'      { 'on' 'off' }      'off';
                        'plot'          'string'      { 'on' 'off' }      'off';
                        'recompute'     'string'      { 'on' 'off' }      'off';
                        'getparams'     'string'      { 'on' 'off' }      'off';
                        'timewindow'    'real'        []      [];
                        'fileout'       'string'  []         '';
                        'timelimits'    'real'        []      [EEG(1).xmin EEG(1).xmax]*1000;
                        'cycles'        'real'        []      [3 .5];
                        'padratio'      'real'        []      1;
                        'freqs'         'real'        []      [0 EEG(1).srate/2];
                        'rmcomps'       'cell'        []      cell(1,length(EEG));
                        'interp'        'struct'      { }     struct([]);
                        'freqscale'     'string'      []      'log';
                        'alpha'         'real'        []      NaN;
                        'type'          'string'      { 'ersp' 'itc' 'both' 'ersp&itc' }  'both'}, 'std_ersp', 'ignore');
if isstr(g), error(g); end;
if isempty(g.trialindices), g.trialindices = cell(length(EEG)); end;
if ~iscell(g.trialindices), g.trialindices = { g.trialindices }; end;

% checking input parameters
% -------------------------
if isempty(g.components) & isempty(g.channels)
    if isempty(EEG(1).icaweights)
        error('EEG.icaweights not found');
    end
    g.components = 1:size(EEG(1).icaweights,1);
    disp('Computing ERSP with default values for all components of the dataset');
end

% select ICA components or data channels
% --------------------------------------
if isempty(g.fileout), g.fileout = fullfile(EEG(1).filepath, EEG(1).filename(1:end-4)); end;
if ~isempty(g.components)
    g.indices = g.components;
    prefix = 'comp';
    filenameersp   = [ g.fileout '.icaersp'  ];
    filenameitc    = [ g.fileout '.icaitc'   ];
    filenametrials = [ g.fileout '.icatimef' ];
    if ~isempty(g.channels)
        error('Cannot compute ERSP/ITC for components and channels at the same time');
    end;
elseif ~isempty(g.channels)
    if iscell(g.channels)
        if ~isempty(g.interp)
            g.indices = eeg_chaninds(g.interp, g.channels, 0);
        else
            g.indices = eeg_chaninds(EEG(1), g.channels, 0);
            for ind = 2:length(EEG)
                if ~isequal(eeg_chaninds(EEG(ind), g.channels, 0), g.indices)
                    error([ 'Channel information must be consistant when ' 10 'several datasets are merged for a specific design' ]);
                end;
            end;
        end;
    else
        g.indices = g.channels;
    end;
    prefix = 'chan';
    filenameersp   = [ g.fileout '.datersp'  ];
    filenameitc    = [ g.fileout '.datitc'   ];
    filenametrials = [ g.fileout '.dattimef' ];
end;

powbaseexist = 1; % used also later
if isempty(g.powbase) | isnan(g.powbase)
    powbaseexist = 0;
    g.powbase = NaN*ones(length(g.indices),1);  % default for timef()
end;
if size(g.powbase,1) ~= length(g.indices)
    error('powbase should be of size (ncomps,nfreqs)');
end

% Check if ERSP/ITC information found in datasets and if fits requested parameters
% ----------------------------------------------------------------------------
if exist( filenameersp ) & strcmpi(g.recompute, 'off')
    fprintf('Use existing file for ERSP: %s\n', filenameersp);
    return;
end;
%    tmpersp  = load( '-mat', filenameersp, 'parameters'); % AND IT SHOULD BE USED HERE TOO - ARNO
%	params   = struct(tmpersp.parameters{:});
%    if ~isequal(params.cycles, g.cycles)                   ...
%            | (g.padratio ~= params.padratio) ...
%            | ( (g.alpha~= params.alpha) & ~( isnan(g.alpha) & isnan(params.alpha)) )
%        % if not computed with the requested parameters, recompute ERSP/ITC
%        % i.e., continue
%    else
%        disp('File ERSP/ITC data already present, computed with the same parameters: no need to recompute...');
%        return; % no need to compute ERSP/ITC
%    end
%end;

% Compute ERSP parameters
% -----------------------
parameters = { 'cycles', g.cycles, 'padratio', g.padratio, ...
               'alpha', g.alpha, 'freqscale', g.freqscale, timefargs{:} };
defaultlowfreq = 3;
[time_range] = compute_ersp_times(g.cycles,  EEG(1).srate, ...
                                 [EEG(1).xmin EEG(1).xmax]*1000 , defaultlowfreq, g.padratio);
if time_range(1) < time_range(2) && g.freqs(1) == 0
     g.freqs(1) = defaultlowfreq; % for backward compatibility
end
parameters = { parameters{:} 'freqs' g.freqs };
if strcmpi(g.plot, 'off')
    parameters = { parameters{:} 'plotersp', 'off', 'plotitc', 'off', 'plotphase', 'off' };
end;
if powbaseexist & time_range(1) >= 0
    parameters{end+1} = 'baseboot';
    parameters{end+1} = 0;
    fprintf('No pre-0 baseline spectral estimates: Using whole epoch for timef() "baseboot"\n');
end

% return parameters
% -----------------
if strcmpi(g.getparams, 'on')
    X = []; times = []; freqs = [];
    if strcmpi(g.savetrials, 'on')
        parameters = { parameters{:} 'savetrials', g.savetrials };
    end;
    return;
end;

% No usable ERSP/ITC information available
% ---------------------------------
% tmpdata = [];
% for index = 1:length(EEG)
%     if isstr(EEG(index).data)
%         TMP = eeg_checkset( EEG(index), 'loaddata' );  % load EEG.data and EEG.icaact
%     else
%         TMP = EEG;
%     end
%     if ~isempty(g.components)
%         if isempty(TMP.icaact)                      % make icaact if necessary
%             TMP.icaact = (TMP.icaweights*TMP.icasphere)* ...
%                           reshape(TMP.data(TMP.icachansind,:,:), [ length(TMP.icachansind) size(TMP.data,2)*size(TMP.data,3) ]);
%         end;
%         tmpdata    = reshape(TMP.icaact, [ size(TMP.icaact,1) size(TMP.data,2) size(TMP.data,3) ]);
%         tmpdata    = tmpdata(g.indices, :,:);
%     else
%         if isempty(tmpdata)
%             tmpdata = TMP.data(g.indices,:,:);
%         else
%             tmpdata(:,:,end+1:end+size(TMP.data,3)) = TMP.data(g.indices,:,:);
%         end;
%     end;
% end;

options = {};
if ~isempty(g.rmcomps), options = { options{:} 'rmcomps' g.rmcomps }; end;
if ~isempty(g.interp),  options = { options{:} 'interp' g.interp }; end;
X       = [];
for dat = 1:length(EEG)
    if strcmpi(prefix, 'comp')
        tmpdata = eeg_getdatact(EEG(dat), 'component', [1:size(EEG(dat).icaweights,1)], 'trialindices', g.trialindices{dat} );
    else
        EEG(dat).data = eeg_getdatact(EEG(dat), 'channel', [1:EEG(dat).nbchan], 'rmcomps', g.rmcomps{dat}, 'trialindices', g.trialindices{dat});
        EEG(dat).trials = size(EEG(dat).data,3);
        EEG(dat).event  = [];
        EEG(dat).epoch  = [];
        if ~isempty(g.interp),
            TMPEEG = eeg_interp(EEG(dat), g.interp, 'spherical');
            tmpdata = TMPEEG.data;
        else
            tmpdata = EEG(dat).data;
        end;
    end;
    if isempty(X), X = tmpdata;
    else
        if size(X,1) ~= size(tmpdata,1), error('Datasets to be concatenated do not have the same number of channels'); end;
        if size(X,2) ~= size(tmpdata,2), error('Datasets to be concatenated do not have the same number of time points'); end;
        X(:,:,end+1:end+size(tmpdata,3)) = tmpdata; % concatenating trials
    end;
end;

% frame range
% -----------
pointrange1 = round(max((g.timelimits(1)/1000-EEG(1).xmin)*EEG(1).srate, 1));
pointrange2 = round(min(((g.timelimits(2)+1000/EEG(1).srate)/1000-EEG(1).xmin)*EEG(1).srate, EEG(1).pnts));
pointrange = [pointrange1:pointrange2];

% Compute ERSP & ITC
% ------------------
all_ersp   = [];
all_trials = [];
all_itc    = [];
for k = 1:length(g.indices)  % for each (specified) component
    if k>size(X,1), break; end; % happens for components
    if powbaseexist
        tmpparams = parameters;
        tmpparams{end+1} = 'powbase';
        tmpparams{end+1} = g.powbase(k,:);
    else
        tmpparams = parameters;
    end;

    % Run timef() to get ERSP
    % ------------------------
    timefdata  = reshape(X(k,pointrange,:), 1, length(pointrange)*size(X,3));
    if strcmpi(g.plot, 'on'), figure; end;
    [logersp,logitc,logbase,times,logfreqs,logeboot,logiboot,alltfX] ...
          = newtimef( timefdata, length(pointrange), g.timelimits, EEG(1).srate, tmpparams{2:end});
    %figure; newtimef( TMP.data(32,:), EEG.pnts, [EEG.xmin EEG.xmax]*1000, EEG.srate, cycles, 'freqs', freqs);
    %figure; newtimef( timefdata, length(pointrange), g.timelimits, EEG.srate, cycles, 'freqs', freqs);
    if strcmpi(g.plot, 'on'), return; end;

    all_ersp = setfield( all_ersp, [ prefix int2str(g.indices(k)) '_ersp'     ], single(logersp ));
    all_ersp = setfield( all_ersp, [ prefix int2str(g.indices(k)) '_erspbase' ], single(logbase ));
    all_ersp = setfield( all_ersp, [ prefix int2str(g.indices(k)) '_erspboot' ], single(logeboot));
    all_itc  = setfield( all_itc , [ prefix int2str(g.indices(k)) '_itc'      ], single(logitc  ));
    all_itc  = setfield( all_itc , [ prefix int2str(g.indices(k)) '_itcboot'  ], single(logiboot));

    if strcmpi(g.savetrials, 'on')
        all_trials = setfield( all_trials, [ prefix int2str(g.indices(k)) '_timef'     ], single( alltfX ));
    end;
end

% Save ERSP into file
% -------------------
all_ersp.freqs      = logfreqs;
all_ersp.times      = times;
all_ersp.datatype   = 'ERSP';
all_itc.freqs       = logfreqs;
all_itc.times       = times;
all_itc.parameters  = parameters;
all_itc.datatype    = 'ITC';
all_trials.freqs    = logfreqs;
all_trials.times    = times;
all_trials.parameters = { options{:} parameters{:} };
all_trials.datatype   = 'TIMEF';

if powbaseexist
    all_ersp.parameters = { parameters{:}, 'baseline', g.powbase };
else
    all_ersp.parameters = parameters;
end;
if ~isempty(g.channels)
    if ~isempty(g.interp)
        all_ersp.chanlabels   = { g.interp(g.indices).labels };
        all_itc.chanlabels    = { g.interp(g.indices).labels };
        all_trials.chanlabels = { g.interp(g.indices).labels };
    elseif ~isempty(EEG(1).chanlocs)
        tmpchanlocs = EEG(1).chanlocs;
        all_ersp.chanlabels   = { tmpchanlocs(g.indices).labels };
        all_itc.chanlabels    = { tmpchanlocs(g.indices).labels };
        all_trials.chanlabels = { tmpchanlocs(g.indices).labels };
    end;
end;

if strcmpi(g.type, 'both') | strcmpi(g.type, 'ersp') | strcmpi(g.type, 'ersp&itc')
    std_savedat( filenameersp, all_ersp);
end;
if strcmpi(g.type, 'both') | strcmpi(g.type, 'itc') | strcmpi(g.type, 'ersp&itc')
    std_savedat( filenameitc , all_itc );
end;
if strcmpi(g.savetrials, 'on')
    std_savedat( filenametrials , all_trials );
end;

author: PatternRecognition

repo: OpenBMI


text:
function [J, grad] = linearRegCostFunction(X, y, theta, lambda)
%LINEARREGCOSTFUNCTION Compute cost and gradient for regularized linear
%regression with multiple variables
%   [J, grad] = LINEARREGCOSTFUNCTION(X, y, theta, lambda) computes the
%   cost of using theta as the parameter for linear regression to fit the
%   data points in X and y. Returns the cost in J and the gradient in grad

% Initialize some useful values
m = length(y); % number of training examples

% You need to return the following variables correctly
J = 0;
grad = zeros(size(theta));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost and gradient of regularized linear
%               regression for a particular choice of theta.
%
%               You should set J to the cost and grad to the gradient.
%
tempTheta = theta;
tempTheta(1) = 0;

J = (1 / (2*m) ) * sum(((X * theta)-y).^2) + (lambda / (2 * m))*sum(tempTheta.^2);
temp = X * theta;
error = temp - y;
grad = (1 / m) * (X' * error) + (lambda/m)*tempTheta;

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

grad = grad(:);

end

author: AvaisP


text:
function title = p19_title ( )

%*****************************************************************************80
%
%% P19_TITLE returns the title of problem P19.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, the title of the problem.
%
  title = 'Problem 19, Enright and Pryce #D4';

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function y = In
y = cvxin;

repo: REST