In = struct;
In.nEventReps = [89,43,56,21]';
In.ContrastMats = {eye(4);[0.5,0.5,-0.5,-0.5];[0.5,-0.5,0.5,-0.5]};
In.ContrastNames = {'OverallEffect';'Main1';'Main2'};
In.nNulls = 99;
In.ITI = 5.5;
In.TR = 2.65;
In.nScans = round((((sum(In.nEventReps,1) + In.nNulls) * In.ITI) + 2) / In.TR);
In.nIterations = 1000;
[Out] = SamsOptseq(In);
scatter3([Out.Main1]',[Out.Main2]',[Out.OverallEffect]',20,zeros(size([Out.OverallEffect]',1),3),'filled');
Betas = regress([Out.OverallEffect]',[ones(size(Out,1),1),[Out.Main1]',[Out.Main2]']);
[X1,X2]=ndgrid(min([Out.Main1]):std([Out.Main1])/10:max([Out.Main1]),min([Out.Main2]):std([Out.Main2])/10:max([Out.Main2]));
Y = Betas(1,1) + Betas(2,1)*X1 + Betas(3,1)*X2;
hold on;
mesh(X1,X2,Y);