function res = BR_unfold(step)
%

vals = [0:step:1];
num_vals = length(vals);
num_pts = num_vals^4;

V = [0 0 1 1; 1 1 0 1; 1 1 1 1];
pb = [0.7;0.7;1];

for c = 1:8
    res(c).pts4D = zeros(num_pts,4);
    res(c).pts2D = zeros(num_pts,2);
end

V = [0 0 0 0;...
     1 0 0 0;...
     1 0 1 0;...
     0 0 1 0;...
     0 1 0 0;...
     1 1 0 0;...
     1 1 1 0;...
     0 1 1 0;...
     0 0 0 1;...
     1 0 0 1;...
     1 0 1 1;...
     0 0 1 1;...
     0 1 0 1;...
     1 1 0 1;...
     1 1 1 1;...
     0 1 1 1];
 
 cube_indexes = ...
     [ 1, 2, 3, 4, 5, 6, 7, 8;...
       9,10,11,12,13,14,15,16;...
       3, 4,11,12, 7, 8,15,16;...
       9,10, 1, 2,13,14, 5, 6;...
       1, 4, 9,12, 5, 8,13,16;...
       2, 3, 6, 7,10,11,14,15;...
       5, 6, 7, 8,13,14,15,16;...
       1, 2, 3, 4, 9,10,11,12];
 
cubes = zeros(8,4,8);

for c = 1:8
    for v = 1:8
        cubes(v,:,c) = V(cube_indexes(c,v),:);
    end
end

for c = 1:8
    pts = cubes(:,:,c);
    m = mean(pts);
    indexes1 = find(abs(m-0.5)<0.1);
    indexes2 = find(abs(m-0.5)>0.1);
    pt = zeros(1,4);
    pt(indexes2) = m(indexes2);
    index = 0;
    for ind1 = 1:num_vals
        pt(indexes1(1)) = vals(ind1);
        for ind2 = 1:num_vals
            pt(indexes1(2)) = vals(ind2);
            for ind3 = 1:num_vals
                pt(indexes1(3)) = vals(ind3);
                pa = pt(3) + pt(4);
                pb = pt(2) + pt(4);
                index = index + 1;
                res(c).pts4D(index,:) = pt;
                res(c).pts2D(index,:) = [pa,pb];
            end
        end
    end
end

tch = 0;

