Posted on

Photo Merge Source

This post contains full Python source code for Photo Merge. Source by Michal Lohnický.

<pre>#!/usr/bin/python

import cv

#######################################
# Class to hold helper functions
#######################################
class Helper():
image_set = 0
max_distance = 10000000
r_plane = None
g_plane = None
b_plane = None
hist1 = None
hist2 = None
is_init = False
color_bins = 2
color_range = [0, 255]

h_plane = None
s_plane = None
h_bins = 4#30
s_bins = 4#32
hist1_hsv = None
hist2_hsv = None

@staticmethod
def init(grid_size):
if Helper.is_init == True: return
Helper.is_init = True
Helper.r_plane = cv.CreateMat(grid_size, grid_size, cv.CV_8UC1)
Helper.g_plane = cv.CreateMat(grid_size, grid_size, cv.CV_8UC1)
Helper.b_plane = cv.CreateMat(grid_size, grid_size, cv.CV_8UC1)
color_bins = Helper.color_bins
color_range = Helper.color_range
Helper.hist1 = cv.CreateHist([color_bins, color_bins, color_bins], cv.CV_HIST_SPARSE, [color_range, color_range, color_range], 1)
Helper.hist2 = cv.CreateHist([color_bins, color_bins, color_bins], cv.CV_HIST_SPARSE, [color_range, color_range, color_range], 1)

Helper.h_plane = cv.CreateMat(grid_size, grid_size, cv.CV_8UC1)
Helper.s_plane = cv.CreateMat(grid_size, grid_size, cv.CV_8UC1)
h_bins = Helper.h_bins
s_bins = Helper.s_bins
ranges = [[0, 180],[0, 255]]
Helper.hist1_hsv = cv.CreateHist([h_bins, s_bins], cv.CV_HIST_ARRAY, ranges, 1)
Helper.hist2_hsv = cv.CreateHist([h_bins, s_bins], cv.CV_HIST_ARRAY, ranges, 1)
Helper.hsv_img = cv.CreateImage((grid_size, grid_size), 8, 3)

@staticmethod
def resize(im, new_width):
size = (new_width, int(im.height/float(im.width)*new_width) )
src_img_1_resized = cv.CreateImage(size, im.depth, im.nChannels)
cv.Resize(im, src_img_1_resized)
return src_img_1_resized

@staticmethod
def get_histogram_hsv(src, hist_used_ID=1):
# Convert to HSV
cv.CvtColor(src, Helper.hsv_img, cv.CV_BGR2HSV)
# Extract the H and S planes
cv.Split(Helper.hsv_img, Helper.h_plane, Helper.s_plane, None, None)
planes = [Helper.h_plane, Helper.s_plane]
if hist_used_ID==1:
cv.ClearHist(Helper.hist1_hsv)
cv.CalcHist([cv.GetImage(i) for i in planes], Helper.hist1_hsv)
return Helper.hist1_hsv
cv.ClearHist(Helper.hist2_hsv)
cv.CalcHist([cv.GetImage(i) for i in planes], Helper.hist2_hsv)
return Helper.hist2_hsv

@staticmethod
def get_histogram(src, hist_used_ID=1):
cv.Split(src, Helper.r_plane, Helper.g_plane, Helper.b_plane, None)
planes = [Helper.r_plane, Helper.g_plane, Helper.b_plane]
if hist_used_ID==1:
cv.ClearHist(Helper.hist1)
cv.CalcHist([cv.GetImage(i) for i in planes], Helper.hist1)
return Helper.hist1
cv.ClearHist(Helper.hist2)
cv.CalcHist([cv.GetImage(i) for i in planes], Helper.hist2)
return Helper.hist2

@staticmethod
def print_histogram(hist):
sum = 0
print "Printing histogram"
for r in range(0,Helper.color_bins):
for g in range(0,Helper.color_bins):
for b in range(0,Helper.color_bins):
intensity = cv.QueryHistValue_3D(hist,r,g,b)
sum+=intensity
if intensity>0.0000001:
print r,g,b,":", intensity
print "Overal histogram sum:",sum

@staticmethod
def print_histogram_double(hist1, hist2):
sum1 = 0
sum2 = 0
print "Printing histogram"
for r in range(0,Helper.color_bins):
for g in range(0,Helper.color_bins):
for b in range(0,Helper.color_bins):
intensity1 = cv.QueryHistValue_3D(hist1,r,g,b)
intensity2 = cv.QueryHistValue_3D(hist2,r,g,b)
sum1+=intensity1
sum2+=intensity2
if intensity1>0.0000001 or intensity2>0.0000001:
print r,g,b,":", intensity1, intensity2
print "Overal histogram sum:",sum1, sum2

mode_list = {
(cv.IPL_DEPTH_8U, 3) : "RGBA",
(cv.IPL_DEPTH_8U, 3) : "RGB",
(cv.IPL_DEPTH_8U, 1) : "L",
(cv.IPL_DEPTH_32F, 1) : "F"
}
mode_list_r = {
"RGBA" : (cv.IPL_DEPTH_8U, 3),
"RGB" : (cv.IPL_DEPTH_8U, 3),
"L" : (cv.IPL_DEPTH_8U, 1),
"F" : (cv.IPL_DEPTH_32F, 1)
}
@staticmethod
def convertToPIL(cv_im):
from PIL import Image
mode = Helper.mode_list[(cv_im.depth, cv_im.nChannels)]
return Image.fromstring(mode, cv.GetSize(cv_im), cv_im.tostring()), mode

@staticmethod
def convertFromPIL(pi):
mode = Helper.mode_list_r[pi.mode]
cv_im = cv.CreateImageHeader(pi.size, mode[0], mode[1])
cv.SetData(cv_im, pi.tostring(), pi.size[0]*mode[1])
return cv_im

@staticmethod
def reduce_colors(cv_im, colors_count):
from PIL import Image
pi, mode = Helper.convertToPIL(cv_im)
pi = pi.convert("RGB").convert("P", palette=Image.ADAPTIVE, colors=colors_count)
pi = pi.convert(mode)
return Helper.convertFromPIL(pi)

@staticmethod
def pasteMask(im1,im2,mask):
pi1, mode = Helper.convertToPIL(im1)
pi2, mode = Helper.convertToPIL(im2)
pimask, mode = Helper.convertToPIL(mask)
pi1.paste(pi2,(0,0),pimask)
return Helper.convertFromPIL(pi1)

#######################################
# Structure to hold matched values
#######################################
class SURFFeatureMatch():
keypoint_img_1 = None
descriptor_img_1 = None
keypoint_img_2 = None
descriptor_img_2 = None
distance = Helper.max_distance

#######################################
# Match feature points of 2 images
#######################################
class SURFFeaturesMatcher():

#=========================
# Match feature points of 2 images (main function)
#=========================
def findPairs(self, img_proc_1, img_proc_2):
matches = []
for i in range(0,len(img_proc_1.keypoints)):
featureMatch = SURFFeatureMatch()
featureMatch.keypoint_img_1 = img_proc_1.keypoints[i]
featureMatch.descriptor_img_1 = img_proc_1.descriptors[i]
#find nearest feature point
for j in range(0,len(img_proc_2.keypoints)):
if featureMatch.keypoint_img_1[1] == img_proc_2.keypoints[j][1]: #laplacians has to be same
if abs(featureMatch.keypoint_img_1[3] - img_proc_2.keypoints[j][3])<45: #similar orientation
if self.realSqrDistance(featureMatch.keypoint_img_1[0], img_proc_2.keypoints[j][0])<4000:
dist = self.descriptorsDistance(featureMatch.descriptor_img_1, img_proc_2.descriptors[j], featureMatch.distance)
if dist current_best: break
return total_cost

#=========================
# Equlide distance TODO consider image size
#=========================
def realSqrDistance(self, p1, p2):
return (p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1])

#######################################
# Image processing class
#######################################
class ImageProcessor:
keypoints = None
descriptors = None
src_img_orig = None
src_img = None
src_img_colored = None
src_img_path = ""
matched_features = None
binded_img_proc = None
top_feature_points_count = 500
grid_size = 10
grid_tamplateMatching_offset = 8

def __init__(self, path):
self.homography = None
self.src_img_path = path
self.src_img = Helper.resize(cv.LoadImage(self.src_img_path, cv.CV_LOAD_IMAGE_GRAYSCALE), 800)
self.src_img_colored = Helper.resize(cv.LoadImage(self.src_img_path, cv.CV_LOAD_IMAGE_COLOR), 800)
self.src_img_colored = Helper.reduce_colors(self.src_img_colored,255)
self.src_img_orig = Helper.resize(cv.LoadImage(self.src_img_path, cv.CV_LOAD_IMAGE_COLOR), 800)
Helper.init(self.grid_size)

def extractSURFFeatures(self):
try:
self.border_w = int(self.src_img.width * 0.1)
self.border_h = int(self.src_img.height * 0.1)
cv.SetImageROI(self.src_img, ( self.border_w, self.border_h, self.src_img.width-2*self.border_w, self.src_img.height-2*self.border_h ) )
(self.keypoints, self.descriptors) = cv.ExtractSURF(self.src_img, None, cv.CreateMemStorage(), (1, 500, 3, 4))
print self.src_img_path+":\n"+" keypoints: "+str(len(self.keypoints))+"\n descriptors: "+str(len(self.descriptors))
cv.ResetImageROI(self.src_img)
except Exception, e:
print e

def findSURFFeaturesPairs(self,img_proc):
if self.keypoints == None: self.extractSURFFeatures()
if img_proc.keypoints == None: img_proc.extractSURFFeatures()
self.binded_img_proc = img_proc
featureMatcher = SURFFeaturesMatcher()
self.matched_features = featureMatcher.findPairs(self,img_proc)

def getImportantFeaturePoints(self, top_points_count):
ret_1 = []
ret_2 = []
for featureMatch in self.matched_features[:top_points_count]:
point_1 = (int(featureMatch.keypoint_img_1[0][0]+self.border_w), int(featureMatch.keypoint_img_1[0][1]+self.border_h))
point_2 = (int(featureMatch.keypoint_img_2[0][0]+self.border_w), int(featureMatch.keypoint_img_2[0][1]+self.border_h))
ret_1.append(point_1)
ret_2.append(point_2)
return ret_1, ret_2

def getHomography(self):
if self.homography!=None: return self.homography
srcKeypoints, dstKeypoints = self.getImportantFeaturePoints(self.top_feature_points_count)
dim=len(srcKeypoints)
srcMat=cv.CreateMat(dim,2,cv.CV_32FC1)
dstMat=cv.CreateMat(dim,2,cv.CV_32FC1)
for i in range(1,dim):
srcMat[i,0]=srcKeypoints[i][0]
srcMat[i,1]=srcKeypoints[i][1]
dstMat[i,0]=dstKeypoints[i][0]
dstMat[i,1]=dstKeypoints[i][1]
h=cv.CreateMat(3,3,cv.CV_64F)
cv.FindHomography(dstMat,srcMat, h,cv.CV_RANSAC,5)
self.homography = h
return h

def warpImage(self):
src_img_warped = cv.CloneImage(self.binded_img_proc.src_img_colored)
src_img_orig_warped = cv.CloneImage(self.binded_img_proc.src_img_orig)
homo = self.getHomography()
cv.WarpPerspective(self.binded_img_proc.src_img_colored, src_img_warped, homo)
cv.WarpPerspective(self.binded_img_proc.src_img_orig, src_img_orig_warped, homo)
return src_img_warped, src_img_orig_warped

def drawSURFFeatures(self):
window_name = "SURF "+self.src_img_path
cv.NamedWindow(window_name, 1)
for ((x, y), laplacian, size, dir, hessian) in self.keypoints:
radio = size*1.2/9.*2
color = (255, 0, 0)
if radio < 3:
radio = 2
color = (0, 255, 0)
cv.Circle(self.src_img_colored, (int(self.border_w+x),int(self.border_h+y)), int(radio), (0,255,0))
cv.ShowImage(window_name, self.src_img_colored)

def drawMatchedPoints(self):
#merge photos
merged_size = (self.src_img.width, self.src_img.height + self.binded_img_proc.src_img.height)
correspond = cv.CreateImage( merged_size, self.src_img.depth, self.src_img.nChannels)
cv.SetImageROI( correspond, ( 0, 0, self.src_img.width, self.src_img.height ) )
cv.Copy( self.src_img, correspond )
cv.SetImageROI( correspond, ( 0, self.src_img.height, correspond.width, correspond.height ) )
cv.Copy( self.binded_img_proc.src_img, correspond )
cv.ResetImageROI( correspond )
#create window
window_name = "Matched features"
cv.NamedWindow(window_name, 1)
#draw lines
points_1, points_2 = self.getImportantFeaturePoints(self.top_feature_points_count)
for i in range(0,len(points_1)):
point_2 = (int(points_2[i][0]), int(points_2[i][1]+ self.src_img.height))
cv.Line( correspond, points_1[i], point_2, 128 )
#show image
cv.ShowImage(window_name, correspond)

def findMid(self, histogram):
#normalize
offset = len(histogram)/4
hist = []
for i in range(offset,len(histogram)-offset):
hist.append((histogram[i-2]*0.3+histogram[i-1]*0.7+histogram[i]+histogram[i+1]*0.7+histogram[i+2]*0.3)/3)
maxi = max(hist)
for i in range(0,len(hist)): hist[i]/=maxi
#find maxes
max1_val = 0
max1_ind = 0
max2_val = 0
max2_ind = 0
ind1 = 0
ind2 = len(hist)-1
while ind1max1_val:
max1_val = hist[ind1]
max1_ind = ind1
if hist[ind2]>max2_val:
max2_val = hist[ind2]
max2_ind = ind2
ind1+=1
ind2-=1
print max1_ind, max2_ind
#find min
mid = len(hist)/2
min_ind = 0
min_val = 10000000
for i in range(max1_ind,max2_ind):
new_val = hist[i]#*((abs(mid-i)/float(mid))*0.3+0.7)
print i, ":", hist[i]
if new_val=mid_ind:
cv.Rectangle(mask2, (0, 0), (grid_size, grid_size), 0, thickness=-1)
cv.ResetImageROI(mask1)
cv.ResetImageROI(mask2)
element_width_half = grid_size*2
element = cv.CreateStructuringElementEx(element_width_half*2+1, element_width_half*2+1, element_width_half, element_width_half, cv.CV_SHAPE_RECT)
cv.Erode(mask2, mask2, element, 1)
cv.Dilate(mask2, mask2, element, 1)
cv.Erode(mask2, mask2, element, 1)
cv.Erode(mask1, mask1, element, 1)
cv.Dilate(mask1, mask1, element, 1)
cv.Erode(mask1, mask1, element, 1)
#finalize
cv.Smooth(mask1, mask1, cv.CV_BLUR, grid_size*2, grid_size*2)
cv.Smooth(mask2, mask2, cv.CV_BLUR, grid_size*2, grid_size*2)
mask1 = Helper.resize(mask1, img_orig1.width)
mask2 = Helper.resize(mask2, img_orig2.width)
img_res1 = Helper.pasteMask(img_orig1,img_orig2,mask1)
img_res2 = Helper.pasteMask(img_orig1,img_orig2,mask2)
img_proc_1.crop(img_res1)
cv.SaveImage("res1"+str(Helper.image_set)+".png",img_res1)
cv.SaveImage("res2"+str(Helper.image_set)+".png",img_res2)

def crop(self,im):
roz = [[0, im.width, 0, im.width],[0, 0, im.height, im.height]]
p = []
for i in range(0,4):
p.append(cv.CreateMat(1,3,cv.CV_64FC1))
p[i][0,0]=roz[0][i]
p[i][0,1]=roz[1][i]
p[i][0,2]=1
cv.MatMul(p[i],self.homography,p[i])
p[i][0,0]/=p[i][0,2]
p[i][0,1]/=p[i][0,2]
cropped = cv.CreateImage( (im.width, im.height), im.depth, im.nChannel)
src_region = cvGetSubRect(image, opencv.cvRect(left, top, new_width, new_height) )
cvCopy(src_region, cropped)

import sys
if __name__ == "__main__":
Helper.image_set = str(sys.argv[1])
image_set = Helper.image_set
img_proc_1 = ImageProcessor(".\\test_imgs\\"+str(image_set).zfill(2)+"_01.JPG")
img_proc_2 = ImageProcessor(".\\test_imgs\\"+str(image_set).zfill(2)+"_02.JPG")
img_final = cv.CreateImage( (img_proc_1.src_img_colored.width, img_proc_1.src_img_colored.height), img_proc_1.src_img_colored.depth, img_proc_1.src_img_colored.nChannels)
img_histograms = cv.CreateImage( (img_proc_1.src_img_colored.width, img_proc_1.src_img_colored.height), 8, 1)

img_proc_1.extractSURFFeatures()
img_proc_2.extractSURFFeatures()
img_proc_1.findSURFFeaturesPairs(img_proc_2)
img_proc_1.drawMatchedPoints()
src_img_warped, src_img_orig_warped = img_proc_1.warpImage()

histograms_diffs = []

templateMatching_offset = 8
templateMatching_result = cv.CreateImage((templateMatching_offset+1, templateMatching_offset+1), cv.IPL_DEPTH_32F, 1)
grid_size = ImageProcessor.grid_size
act_width = templateMatching_offset/2
while act_width+grid_size(maxi-mini)*0.5:
grid_mask_cols.append(1)
frame_thickness = 1
if first and act_width<400 and act_height<230 and act_width>300 and act_height>180:
hist1 = Helper.get_histogram(img_proc_1.src_img_colored, hist_used_ID=1)
hist2 = Helper.get_histogram(src_img_warped, hist_used_ID=2)
frame_thickness = -1
first = False
cv.Copy( img_proc_1.src_img_colored, img_final )
frame_color = 128
cv.Rectangle(img_histograms, (1, 1), (grid_size-1, grid_size-1), 0)
else:
grid_mask_cols.append(0)
cv.Copy( src_img_warped, img_final )
grid_mask.append(grid_mask_cols)

cv.ResetImageROI(img_proc_1.src_img_colored)
cv.ResetImageROI(src_img_warped)
cv.ResetImageROI(img_final)
cv.ResetImageROI(img_histograms)

img_proc_1.combinePictures(img_proc_1.src_img_orig, src_img_orig_warped, img_proc_1.src_img_colored, src_img_warped, grid_mask)

cv.ShowImage("tmp",img_proc_1.src_img_colored)
cv.ShowImage("tmp2",src_img_warped)
cv.ShowImage("tmp1",img_histograms)
cv.ShowImage("FINAL",img_final)
cv.WaitKey()
cv.DestroyAllWindows()</pre>