In [301]:
# first of all, let us define some parameters
#
# size of the images
sx = 400
sy = 400
# number of sources on each image
nsrc1 = 3
nsrc2 = nsrc1
# typical error radius (in pixels)
rerr1 = 40
rerr2 = rerr1
In [302]:
# generate coordinate pairs for each image
def generate_positions(npts,img_shape):
"""
Generate 'npts' points uniformly across 'image_shape'.
Args:
npts : number of points to generate
img_shape : (y,x) shape where to generate points
Returns:
Pair_Coordinates_List : list of (y,x) tuples
"""
import numpy as np
_sy,_sx = img_shape
assert _sy>=5 and _sx>=5 # because I want
indy = np.random.random_integers(0,_sy-1,npts)
indx = np.random.random_integers(0,_sx-1,npts)
_inds = zip(indy,indx)
return _inds
# "sources 1"
coords1 = generate_positions(nsrc1,(sy,sx))
assert isinstance(coords1,list) and len(coords1) is nsrc1
In [303]:
# create the positions table
def create_positions_table(coords,err_radius):
"""
"""
tab = {}
for i,oo in enumerate(coords):
tab[i] = [i,oo[1],oo[0],err_radius]
return tab
# table for "sources 1"
tab1 = create_positions_table(coords1,rerr1)
In [304]:
# create and draw each source on black(null) images
def draw_image_sources(tab_positions,img_shape):
"""
"""
from math import ceil
from PIL import Image,ImageDraw
assert(isinstance(img_shape,tuple) and len(img_shape) is 2)
size = img_shape[::-1]
# Modification to accomplish color codes ---
#mode = 'L'
mode = 'RGB'
# ---
color = "black"
img = Image.new(mode,size,color)
assert(len(tab_positions)>=1)
dictColorId = {}
for i,src in tab_positions.items():
assert isinstance(src,list) and src is tab_positions[i]
assert len(src)>=4, "length of table raw %d is %d" % (i,len(src))
assert i==src[0]
draw = ImageDraw.Draw(img)
x = src[1]
assert 0<=x and x<size[0], "coordinate x is %d" % x
y = src[2]
assert 0<=y and y<size[1], "coordinate y is %d" % y
r = src[3]
assert r<size[0]/2 and r<size[1]/2
box = (x-r,y-r,x+r,y+r)
# Modification to accomplish color codes ---
#fill=255
_R = int(255 - ( int(x/256) + int(y/256)*(1 + ceil(size[0]/256)) )) #TODO: restrict total size of image to avoid _R<=0
_G = x%256
_B = y%256
fill = (_R,_G,_B)
# ---
dictColorId[str(fill)] = i
print i,fill,src
draw.ellipse(box,fill=fill)
del draw,box,x,y,r
return img,dictColorId
img1,cor2id1 = draw_image_sources(tab1,(sy,sx))
img1.show()
In [305]:
# do the same steps for "sources 2"
coords2 = generate_positions(nsrc2,(sy,sx))
assert isinstance(coords2,list) and len(coords2) is nsrc2
tab2 = create_positions_table(coords2,rerr2)
img2,cor2id2 = draw_image_sources(tab2,(sy,sx))
img2.show()
In [306]:
# intersect the images
def intersect_images(img1,img2):
"""
"""
from PIL import Image
assert isinstance(img1,Image.Image) and isinstance(img2,Image.Image)
import numpy as np
im1 = np.asarray(img1,dtype=bool)
im2 = np.asarray(img2,dtype=bool)
assert im1.shape == im2.shape
imout = np.zeros(im1.shape,dtype=bool)
np.logical_and(im1,im2,imout)
imgout = Image.fromarray(np.uint8(imout)*255)
del im1,im2
return imgout#,imout
#imgs_pack = intersect_images(img1,img2)
#imgx = imgs_pack[0]
imgx = intersect_images(img1,img2)
imgx.show()
In [307]:
# save the figs
img1.save('sources1.png')
img2.save('sources2.png')
imgx.save('sourcesx.png')
In [308]:
# At this point we have the matching of the (two) catalogs. So far, we (roughly) know how many objects had matched.
# Now we need to know which objects match each other. We need now iterate through the matched regions and query
# which objects are there.
In [309]:
# Now we want to label each of the intersections so that we can iterate each of them...
# OBS: to label an image, it needs to be a two dimensional array.
# So we need to convert the match (RGB,3D) image to a matrix (2D).
def rgb2bool(pil_image):
"""
"""
from PIL.Image import Image
import numpy as np
assert isinstance(pil_image,Image)
# there are two ways for doing it
if True:
# 1) converting still on PIL and then to a numpy.array
imgL = pil_image.convert('L')
imC = np.uint8(imgL)
imC = imC.astype(bool)
else:
# 2) transform to a numpy.array and then collapse the RGB layers
imB = np.asarray(pil_image,dtype=np.uint8)
imC = np.sum(imB,axis=2)
imC = imC.astype(bool)
return imC
imxbin = rgb2bool(imgx)
assert len(imxbin.shape)==2
# and then, label it
from scipy import ndimage as ndi
regionsx,numlabels = ndi.label(imxbin)
In [310]:
# (let us transform the color images back from PIL to Numpy)
import numpy as np
def get_colors(indices,image):
# query image-1
pxcor = image[indices]
# now we want to remove duplicates
pcl = [ list(raw) for raw in pxcor ]
pcl.sort(key = lambda l:tuple(l))
pxcor = np.array(pcl)
_diff = np.diff(pxcor,axis=0)
_mask = np.ones(len(pxcor),'bool')
_mask[1:] = (_diff!=0).any(axis=1)
return pxcor[_mask]
colors = []
cnt1 = 0
for rgb in pcl:
cnt1 = cnt1+1
print "rgb in pcl: %d" % cnt1
if len(colors) == 0:
print "0: %s" % rgb
colors.append(rgb)
continue
cnt2 = 0
for t in colors:
cnt2 = cnt2+1
print "t in colors: %d " % cnt2
#print "t in colors: %s ; len(colors): %d" % (str(t),len(colors))
if rgb[0]==t[0] and rgb[1]==t[1] and rgb[2]==t[2]:
print "t==rgb: %s , %s" % (t,rgb)
else:
print "t!=rgb: %s, %s" % (t,rgb)
colors.append(tuple(rgb))
return colors[:]
im1 = np.asarray(img1)
im2 = np.asarray(img2)
In [311]:
# iterate over the regions querying each image for the respective objects
tojoinTab1 = {}
tojoinTab2 = {}
for i in xrange(1,numlabels+1):
# get the index for label-i
indx = np.where(regionsx==i)
colors1 = get_colors(indx,im1)
tojoinTab1[i] = colors1
colors2 = get_colors(indx,im2)
tojoinTab2[i] = colors2
def print_join_tables(dictTab1,cor2id1,dictTab2,cor2id2):
"""
"""
print "%-15s \t\t %-15s" % ("Image-1:","Image-2:")
for i in xrange(1,numlabels+1):
print "---label: %d" % i
l1 = dictTab1[i].tolist()
l2 = dictTab2[i].tolist()
_n = max(len(l1),len(l2))
for ni in xrange(_n):
if len(l1):
rgb1 = tuple(l1.pop())
id1 = cor2id1[str(rgb1)]
else:
rgb1 = ''
id1 = ''
if len(l2):
rgb2 = tuple(l2.pop())
id2 = cor2id2[str(rgb2)]
else:
rgb2 = ''
id2 = ''
print "%-15s [%s] \t %-15s [%s]" % (rgb1,id1,rgb2,id2)
print_join_tables(tojoinTab1,cor2id1,tojoinTab2,cor2id2)
In [312]:
# let us join the tables
def get_coordinates(ID,tab):
"""
"""
src = tab[ID]
assert ID==src[0]
x = src[1]
y = src[2]
return (x,y)
match_table = []
for i in xrange(1,numlabels+1):
l1 = tojoinTab1[i].tolist()
l2 = tojoinTab2[i].tolist()
_n = max(len(l1),len(l2))
for ni in xrange(_n):
if len(l1):
rgb1 = tuple(l1.pop())
id1 = cor2id1[str(rgb1)]
if len(l2):
rgb2 = tuple(l2.pop())
id2 = cor2id2[str(rgb2)]
coords1 = get_coordinates(id1,tab1)
coords2 = get_coordinates(id2,tab2)
match_raw = (i,id1,coords1,id2,coords2)
match_table.append(match_raw)
print "match_ID \t tab1_ID tab1_coords \t tab2_ID tab2_coords"
for mraw in match_table:
print "%d \t\t %d %s \t\t %d %s" % (mraw[0],mraw[1],mraw[2],mraw[3],mraw[4])
_x1,_y1 = mraw[2]
_x2,_y2 = mraw[4]
regionsx[_y1,_x1] = mraw[0]
regionsx[_y2,_x2] = mraw[0]