function I1 = syn_Img(r, ps)

range = -r:r;
[cx, cy, cz] = meshgrid(range, range, range);
c1 = zeros(2*r+1,2*r+1,2*r+1);
c1(cx.^2 + cy.^2 + cz.^2 < r^2) = 255;
I1 = zeros(ps);

sx = [8 36 91 45 50 76 85 67 51 21];
sy = [23 12 56 43 87 63 89 51 61 92];
sz = [8 10 9 14 17 10 11 15 13 10 11];

for kk = 1:10
% 	set = [randi([r+1,ps(1)-r]) randi([r+1,ps(2)-r]) randi([r+1,ps(3)-r])];
% 	I1(set(1)-r:set(1)+r, set(2)-r:set(2)+r, set(3)-r:set(3)+r)= c1;
	I1(sx(kk)-r:sx(kk)+r, sy(kk)-r:sy(kk)+r, sz(kk)-r:sz(kk)+r)= c1;
end

% -------------------------------------------------------------------