ps = [128, 128, 128];
n = 50;
I1 = zeros(ps);
for i = 1:n
	
	r1 = randi([4 10]);
	r2 = r1;
	range1 = -r1:r1;
	range2 = -r2:r2;
	[cx, cy, cz] = meshgrid(range1, range1, range2);
	c1 = zeros(2*r1+1,2*r1+1,2*r2+1);
	c1(cx.^2/r1^2 + cy.^2/r1^2 + cz.^2/r2^2 <= 1) = 1;
% 	c2 = zeros(2*r+1,2*r+1,2*r+1);
% 	c2(cx.^2 + cy.^2 + cz.^2 < (r-2)^2) = 230;
% 	c = c1 - c2;
	% Get random coordinates in the big image where we will place the upper left corner of the small circle.
	dx = ps(1) - 2 * r1;
	dy = ps(2) - 2 * r1;
	dz = ps(3) - 2 * r2;
	sx = dx * rand(n, 1);
	sy = dy * rand(n, 1);
	sz = dz * rand(n, 1);

	% Find the square in the big image where we're going to add a small circle.
	x1 = int16(sx(i));
	y1 = int16(sy(i));
	z1 = int16(sz(i));
	x2 = int16(x1 + 2*r1);
	y2 = int16(y1 + 2*r1);
	z2 = int16(z1 + 2*r2);

	% Add in one small circle to the existing big image
	I1(x1:x2, y1:y2, z1:z2) = I1(x1:x2, y1:y2, z1:z2) + c1;
	I = single(I1);
% 	I1(I1==0) = 150;
end
I (I>1) = 1;
fid = fopen('128-128-128/s50-r4-10.vol', 'w');
fwrite(fid, I1, '*single');
fclose(fid);
% -------------------------------------------------------------------