*Deprecated*
Normalized Cross-Correlation Script
By TC
For usage see Normalized Cross-Correlation
Header=''' # USAGE: python norm_xcorr.py [-option] outfile image1 image2 # # -o Use this followed by 'outfile' to # specify a unique output destination. # Default is corrOut.png # # -h Use this to print usage directions # to StdOut. (Prints the header) # # -v Use this to only output current # version number and exit. By default # the version will be sent to StdOut # at the beginning of each use. # ######################################################## ''' ##~VERSION NOTES~## # 1.0 - first release # 1.1 - added dx,dy print out for template shift # - added sub-sample of template for equal size images # - few minor style changes # 1.2 - added pos. axis display to the output # - removed secondary "normalization" # - added printout of max correlation # - fixed output to remove blank 4th plot # 1.2.1 - removed option print redundancies # 1.3 - fixed 0.5 pixel dx, dy bias # - added colorbar to plot # - scaled correlation plot between 0.2 and 1 # - fixed data print out to remove [] # - changed marker color and type for visibility # - fixed autoscaled ax3 ## ##~FILE DEPENDENCIES~## # User specified: # - image files, the two images to correlate (small first) # - output image file, plotted figure showing matched # location and both images # # Required: # N/A ## ####################~INITIALIZE~#################### import re import sys import numpy as np from scipy.ndimage import convolve from scipy.fftpack import fftn, ifftn from matplotlib import pyplot as plt version = '1.3' opt = sys.argv efile = 'corrOut.png' if len(opt) == 1: print('You forgot to specify the image files to compare!') im1,im2 = raw_input('Please enter the image filenames now: \n').split() elif len(opt) == 3 and opt[1][0] != '-': im1 = opt[1] im2 = opt[2] elif opt[1][0] == '-': if opt[1][1] == 'v': sys.exit('Version: '+version) elif opt[1][1] == 'h': sys.exit(Header) elif opt[1][1] == 'o': efile = opt[2] im1 = opt[3] im2 = opt[4] else: sys.exit('Try again :(') print('norm_xcorr.py version: '+version) print('Image to search for: '+im1) print('Image to search in: '+im2) print('Output figure file: '+efile+'\n') ####################~OBJECT~#################### class normX(object): def __init__(self,img): self.img = img def __call__(self,a): if a.ndim != self.img.ndim: raise Exception('Search area must have the same '\ 'dimensions as the template') return norm_xcorr(self.img,a) ####################~MODULES~#################### def readPGM(file): f=open(file, 'r') all = f.read() try: header, width, height, maxval = re.search( b"(^P5\s(?:\s*#.*[\r\n])*" b"(\d+)\s(?:\s*#.*[\r\n])*" b"(\d+)\s(?:\s*#.*[\r\n])*" b"(\d+)\s(?:\s*#.*[\r\n]\s)*)", all).groups() except AttributeError: raise ValueError("Not a raw PGM file: '%s'" % file) return np.frombuffer(all, dtype='u1' if int(maxval) < 256 else '< u2', count=int(width)*int(height), offset=len(header) ).reshape((int(height), int(width))) def norm_xcorr(t,a): if t.size <= 2: raise Exception('Image is too small.') std_t,mean_t = np.std(t),np.mean(t) if std_t == 0: raise Exception('The image is blank.') t = np.float64(t) a = np.float64(a) outdim = np.array([a.shape[i]+t.shape[i]-1 for i in xrange(a.ndim)]) spattime, ffttime = get_time(t,a,outdim) if spattime < ffttime: method = 'spatial' else: method = 'fourier' if method == 'fourier': af = fftn(a,shape=outdim) tf = fftn(nflip(t),shape=outdim) xcorr = np.real(ifftn(tf*af)) else: xcorr = convolve(a,t,mode='constant',cval=0) ls_a = lsum(a,t.shape) ls2_a = lsum(a**2,t.shape) xcorr = padArray(xcorr,ls_a.shape) ls_diff = ls2_a-(ls_a**2)/t.size ls_diff = np.where(ls_diff < 0,0,ls_diff) sigma_a = np.sqrt(ls_diff) sigma_t = np.sqrt(t.size-1.)*std_t den = sigma_t*sigma_a num = (xcorr - ls_a*mean_t) tol = np.sqrt(np.finfo(den.dtype).eps) nxcorr = np.where(den < tol,0,num/den) # nxcorr = np.where((np.abs(nxcorr)-1.) > np.sqrt(np.finfo(nxcorr.dtype).eps),0,nxcorr) # nxcorr = nxcorr/np.abs(nxcorr).max() nxcorr = padArray(nxcorr,a.shape) return nxcorr def lsum(a,tsh): a = pad(a,tsh) def shiftdiff(a,tsh,shiftdim): ind1 = [slice(None,None),]*a.ndim ind2 = [slice(None,None),]*a.ndim ind1[shiftdim] = slice(tsh[shiftdim],a.shape[shiftdim]-1) ind2[shiftdim] = slice(0,a.shape[shiftdim]-tsh[shiftdim]-1) return a[ind1] - a[ind2] for i in xrange(a.ndim): a = np.cumsum(a,i) a = shiftdiff(a,tsh,i) return a def get_time(t,a,outdim): k_conv = 1.21667E-09 k_fft = 2.65125E-08 convtime = k_conv*(t.size*a.size) ffttime = 3*k_fft*(np.prod(outdim)*np.log(np.prod(outdim))) return convtime,ffttime def pad(a,sh=None,padval=0): if sh == None: sh = np.ones(a.ndim) elif np.isscalar(sh): sh = (sh,)*a.ndim padsize = [a.shape[i]+2*sh[i] for i in xrange(a.ndim)] b = np.ones(padsize,a.dtype)*padval ind = [slice(np.floor(sh[i]),a.shape[i]+np.floor(sh[i])) for i in xrange(a.ndim)] b[ind] = a return b def padArray(a,target,padval=0): b = np.ones(target,a.dtype)*padval aind = [slice(None,None)]*a.ndim bind = [slice(None,None)]*a.ndim for i in xrange(a.ndim): if a.shape[i] > target[i]: diff = (a.shape[i]-target[i])/2. aind[i] = slice(np.floor(diff),a.shape[i]-np.ceil(diff)) elif a.shape[i] < target[i]: diff = (target[i]-a.shape[i])/2. bind[i] = slice(np.floor(diff),target[i]-np.ceil(diff)) b[bind] = a[aind] return b def nflip(a): ind = (slice(None,None,-1),)*a.ndim return a[ind] ####################~READ DATA~#################### temp = readPGM(im1) #temp = temp[20:120,50:150] #for debugging temp=reg reg = readPGM(im2) if temp.size >= reg.size: print('ERROR: Image must be smaller than search area.') print('Template size: '+str(temp.shape)) ans = raw_input('Would you like to sub-sample search template? (y/n) \n') while ans != 'y' and ans != 'n': ans = raw_input('Would you like to sub-sample search template? (y/n) \n') if ans == 'y': x1,x2,y1,y2 = raw_input('Enter px/ln box (min 50px x 50px) to use as x1 x2 y1 y2: \n').split() x1 = int(x1) x2 = int(x2) y1 = int(y1) y2 = int(y2) temp = temp[y1:y2,x1:x2] if ans == 'n': sys.exit('Choose a smaller image.') A = normX(temp) ncc = A(reg) nccloc = np.nonzero(ncc == ncc.max()) x = int(nccloc[1]) y = int(nccloc[0]) expx = int(float(reg.shape[1])/2) expy = int(float(reg.shape[0])/2) dx = x-expx dy = y-expy ####################~OUTPUT~#################### axes = '\n.----> +x \n| \n| \nv \n+y' print('Found match in search area at px/ln: '+str(x)+' '+str(y)) print('Max correlation value: '+str(ncc.max())) print('Template moved dx = '+str(dx)+' dy = '+str(dy)+ \ ' pixels from center of search area. '+axes) fig = plt.figure() ax1 = plt.subplot(2,2,1) ax1.imshow(reg,plt.cm.gray,interpolation='nearest') ax1.set_title('Search Image') ax2 = plt.subplot(2,2,2) ax2.imshow(temp,plt.cm.gray,interpolation='nearest') ax2.set_title('Search Template') ax3 = plt.subplot(2,2,3) ax3.hold(True) im = ax3.imshow(ncc,vmin=0.2,vmax=1,interpolation='nearest') ax3.plot(x,y,'rs',ms=6,fillstyle='none') ax3.set_title('Normalized Cross-Correlation') plt.xlim(0,reg.shape[1]) plt.ylim(reg.shape[0],0) cax = fig.add_axes([0.9, 0.1, 0.03, 0.8]) fig.colorbar(im,cax=cax) plt.savefig(efile)