Help! How can I create a normalized, background subtracted density map?

[Martianartist]

So, being new to python, I was very happy how short the code was to
create a density map from my 2D data. I shamelessly copied the
following code from a website and it works great (http://
pelican.rsvs.ulaval.ca/mediawiki/index.php/
Making_density_maps_using_Gnuplot) :
#!/usr/bin/python

from numpy import *

data_file='AGBstars'
histogram_file= 'AGBstars_hist'

# Here, angles are defined from [-180 to 180[ degrees.
# You could choose another domain. Points will be automatically
# wrapped inside that domain. This can be useful for symetrical
# side chains, for instance. Be careful to use an appropriate domain,
# otherwise the wrapping will produce meaningless data.
x_min, x_max, y_min, y_max = 275.0, 295.0, -37.0, -27.0

# Number of 2D regions in which the plot is divided.
x_resolution, y_resolution = 25, 25

def read_angles(line):
tokens = line.split()
x = float(tokens[1])
y = float(tokens[2])
while x < x_min:
x = x_max - (x_min - x)
while x >= x_max:
x = x_min + (x - x_max)
while y < y_min:
y = y_max - (y_min - y)
while y >= y_max:
y = y_min + (y - y_max)
return [x, y]

points = [read_angles(line) for line in open(data_file) if line[0] !
='#']
count = len(points)
histogram = zeros([x_resolution, y_resolution])
x_interval_length = (x_max - x_min) / x_resolution
y_interval_length = (y_max - y_min) / y_resolution
interval_surface = x_interval_length * y_interval_length
increment = 1000.0 / count / interval_surface

for i in points:
x = int((i[0] - x_min) / x_interval_length)
y = int((i[1] - y_min) / y_interval_length)
histogram[x,y] += increment

x_intervals = arange(x_min, x_max, (x_max - x_min) / x_resolution)
y_intervals = arange(y_min, y_max, (y_max - y_min) / y_resolution)

o = open(histogram_file, 'w')
for i, x in enumerate(x_intervals):
for j, y in enumerate(y_intervals):
o.write('%f %f %f \n' % (x, y, histogram[i,j]))
o.write('\n')
print histogram.max()

I have created several lovely density maps, but the problem is...
How do I subtract a background intensity distribution from my density
map? I created a density map with space coordinates and intensity,
but I want to normalize, and subtract the underlying distribution (one
histogram map file) from another blended map (a separate histogram
file which has background+something interesting). The problem is
first, I am not sure if I have really normalized each file as the
number of counts is very different, and second, even if each
individual file was normalized what I want is something like a
background subtracted image, which would hopefully show only 'true'
overdensities. Thanks