def bg_compensate(img, sigma, splinepoints, scale):
    '''Reads file, subtracts background. Returns [compensated image, background].'''
    
    from PIL import Image
    import pylab
    from matplotlib.image import pil_to_array
    from centrosome.filter import canny
    import matplotlib

    img = Image.open(img)
    if img.mode=='I;16':
        # 16-bit image
        # deal with the endianness explicitly... I'm not sure
        # why PIL doesn't get this right.
        imgdata = np.fromstring(img.tostring(),np.uint8)
        imgdata.shape=(int(imgdata.shape[0]/2),2)
        imgdata = imgdata.astype(np.uint16)
        hi,lo = (0,1) if img.tag.prefix == 'MM' else (1,0)
        imgdata = imgdata[:,hi]*256 + imgdata[:,lo]
        img_size = list(img.size)
        img_size.reverse()
        new_img = imgdata.reshape(img_size)
        # The magic # for maximum sample value is 281
        if 281 in img.tag:
            img = new_img.astype(np.float32) / img.tag[281][0]
        elif np.max(new_img) < 4096:
            img = new_img.astype(np.float32) / 4095.
        else:
            img = new_img.astype(np.float32) / 65535.
    else:
        img = pil_to_array(img)
    
    pylab.subplot(1,3,1).imshow(img, cmap=matplotlib.cm.Greys_r)
    pylab.show()
    
    if len(img.shape)>2:
        raise ValueError('Image must be grayscale')

## Create mask that will fix problem when image has black areas outside of well
    edges = canny(img, np.ones(img.shape, bool), 2, .1, .3)
    ci = np.cumsum(edges, 0)
    cj = np.cumsum(edges, 1)
    i,j = np.mgrid[0:img.shape[0], 0:img.shape[1]]
    mask = ci > 0
    mask = mask & (cj > 0)
    mask[1:,:] &= (ci[0:-1,:] < ci[-1,j[0:-1,:]])
    mask[:,1:] &= (cj[:,0:-1] < cj[i[:,0:-1],-1])
    
    import time
    t0 = time.clock()
    bg = backgr(img, mask, MODE_AUTO, sigma, splinepoints=splinepoints, scale=scale)
    print("Executed in %f sec" % (time.clock() - t0))
    bg[~mask] = img[~mask]

    pylab.subplot(1,3,2).imshow(img - bg, cmap=matplotlib.cm.Greys_r)
    pylab.subplot(1,3,3).imshow(bg, cmap=matplotlib.cm.Greys_r)
    pylab.show()