from pylab import *

def ringsprojfn(shape, rcen, rs, drs, shell=False):
    ''' Draw the projection of rings.

        shape - shape of image to return (row, col format)
        r0 - the center in (row, col) format
        rs - starting radii
        drs - thicknesses of radii
    '''
    n = len(rs)

    x = np.arange(shape[0])
    y = np.arange(shape[1])
    # X is rows, Y is columns
    X,Y = np.meshgrid(x,y,indexing='ij')
    r = np.sqrt((X-rcen[0])**2 + (Y-rcen[1])**2)

    img = np.zeros(shape)

    if shell:
        for i in range(n):
            img += ((1-((r-rs[i])/(drs[i]*.5))**2) > 0)
    else:
        for i in range(n):
            img += np.sqrt(np.maximum(1-((r-rs[i])/(drs[i]*.5))**2,0))

    return img

rs = [10, 100,200]
drs = [5, 20, 40]
shape = 1000,1000
rcen = 500,500
img1 = ringsprojfn(shape, rcen, rs, drs)
img2 = ringsprojfn(shape, rcen, rs, drs,shell=True)

figure(3);clf();
imshow(img1)

figure(4);clf();
imshow(img2)