pyvision.point.DetectorDOG

1 # PyVision License 2 # 3 # Copyright (c) 2006-2008 David S. Bolme 4 # All rights reserved. 5 # 6 # Redistribution and use in source and binary forms, with or without 7 # modification, are permitted provided that the following conditions 8 # are met: 9 # 10 # 1. Redistributions of source code must retain the above copyright 11 # notice, this list of conditions and the following disclaimer. 12 # 13 # 2. Redistributions in binary form must reproduce the above copyright 14 # notice, this list of conditions and the following disclaimer in the 15 # documentation and/or other materials provided with the distribution. 16 # 17 # 3. Neither name of copyright holders nor the names of its contributors 18 # may be used to endorse or promote products derived from this software 19 # without specific prior written permission. 20 # 21 # 22 # THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 23 # ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 24 # LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR 25 # A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR 26 # CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 27 # EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 28 # PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 29 # PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 30 # LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 31 # NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 32 # SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 33 34 ''' 35 Purpose: 36 Find regions of interest using the multiscale differance of 37 gaussians opereator. 38 39 Inspired by: 40 M. Brown and D. Lowe. Invariant features from interest point groups. 41 British Machine Vision Conference. 2002. 42 43 Contributed by: 44 David Bolme 2008. 45 ''' 46 47 48 from pyvision.point.DetectorROI import DetectorROI 49 import unittest 50 import os.path 51 import pyvision as pv 52 #from pyvision.types.img import Image 53 #from pyvision.types.Point import Point 54 from scipy.ndimage.interpolation import zoom 55 from scipy.ndimage.filters import gaussian_filter,maximum_filter,minimum_filter 56 from scipy.linalg import lstsq 57 58 #from pyvision.analysis.ImageLog import ImageLog 59 from numpy import sqrt,array,nonzero 60 #from numpy.lib.polynomial import polyfit 61 62 63 EXTREMA_SIZE=3 64 DEFAULT_SIGMA=sqrt(2) # from Lowe99 65 DEFAULT_TYPE='d' 66 DOG_SCALES = 2 67 MAX_EXTREMA = 1024 68

69 -class DetectorDOG(DetectorROI):

70 - def __init__(self, sigma=DEFAULT_SIGMA, scales=DOG_SCALES, min_size=20, min_contrast=0.03, max_curvature_ratio=10, **kwargs):

71 ''' 72 min_size - Image pyramid terminates when one of the image demensions reaches this size. 73 ''' 74 DetectorROI.__init__(self,**kwargs) 75 76 self.min_size = min_size 77 self.scales = scales 78 self.sigma = sigma 79 self.min_contrast = min_contrast 80 self.max_curvature_ratio = max_curvature_ratio

81 82 83

84 - def _detect(self,im):

85 86 mat = im.asMatrix2D() 87 88 levels = [] 89 while mat.shape[0] > self.min_size and mat.shape[1] > self.min_size: 90 levels.append(mat) 91 mat = zoom(mat,0.5) 92 93 gaussians = [] 94 k = 2.0**(1.0/self.scales) 95 for level in levels: 96 gs = [] 97 sigma = self.sigma 98 for _ in range(self.scales+3): 99 g = gaussian_filter(level,sigma) 100 gs.append(g) 101 sigma = k*sigma 102 gaussians.append(gs) 103 104 dogs = [] 105 for gs in gaussians: 106 ds = [] 107 for i in range(len(gs)-1): 108 d = gs[i]-gs[i+1] 109 ds.append(d) 110 ds = array(ds) 111 dogs.append(ds) 112 #dogs = array(dogs,'d') 113 114 #points = [] 115 sigma = 2*k*self.sigma # approx 95% bounds 116 extrema = [] 117 scale = 1 118 for ds in dogs: 119 # find extrema 120 mins = minimum_filter(ds,(3,3,3)) 121 maxs = maximum_filter(ds,(3,3,3)) 122 123 # Find the extrema but not on the edges of the image 124 minima = nonzero(mins[1:-1,1:-1,1:-1] == ds[1:-1,1:-1,1:-1]) 125 maxima = nonzero(maxs[1:-1,1:-1,1:-1] == ds[1:-1,1:-1,1:-1]) 126 127 for i in range(len(minima[0])): 128 # Correct for removing the edges in the previous step 129 s = minima[0][i]+1 130 x = minima[1][i]+1 131 y = minima[2][i]+1 132 133 # Get a 3 by 3 block 134 block = ds[s-1:s+2,x-1:x+2,y-1:y+2] 135 params = TaylorFit(block) 136 ts,tx,ty = TaylorSubpixel(params) 137 td = TaylorApprox(params,ts,tx,ty) 138 ts -= 1.0 139 tx -= 1.0 140 ty -= 1.0 141 142 # Only select extrema with high contrast 143 if abs(td) < self.min_contrast: continue 144 145 # Chech the ratios of the principal curvatures (see Lowe 2004 Sec 4.1: 146 Dxx=2*params[1] 147 Dyy=2*params[2] 148 Dxy=params[3] 149 TrH = Dxx+Dyy 150 DetH = Dxx*Dyy-Dxy*Dxy 151 152 r = self.max_curvature_ratio 153 if DetH < 0: continue # discard because curvatures have different signs 154 if r*TrH > DetH*(r+1)*(r+1): continue # Ratio of curvatuers is greater than R 155 156 if not (-1.0 < tx and tx < 1.0 and -1.0 < ty and ty < 1.0): continue 157 158 extrema.append([-td,scale*(x+tx),scale*(y+ty),(k**((s+ts)-1))*sigma]) 159 160 for i in range(len(maxima[0])): 161 s = maxima[0][i]+1 162 x = maxima[1][i]+1 163 y = maxima[2][i]+1 164 165 # Get a 3 by 3 block 166 block = ds[s-1:s+2,x-1:x+2,y-1:y+2] 167 params = TaylorFit(block) 168 ts,tx,ty = TaylorSubpixel(params) 169 td = TaylorApprox(params,ts,tx,ty) 170 ts -= 1.0 171 tx -= 1.0 172 ty -= 1.0 173 174 # Only select extrema with high contrast 175 if abs(td) < self.min_contrast: continue 176 177 # Chech the ratios of the principal curvatures (see Lowe 2004 Sec 4.1: 178 Dxx=2*params[1] 179 Dyy=2*params[2] 180 Dxy=params[3] 181 TrH = Dxx+Dyy 182 DetH = Dxx*Dyy-Dxy*Dxy 183 184 r = self.max_curvature_ratio 185 if DetH < 0: continue # discard because curvatures have different signs 186 if r*TrH > DetH*(r+1)*(r+1): continue # Ratio of curvatuers is greater than R 187 188 if not (-1.0 < tx and tx < 1.0 and -1.0 < ty and ty < 1.0): continue 189 190 extrema.append([td,scale*(x+tx),scale*(y+ty),(k**((s+ts)-1))*sigma]) 191 192 sigma = (k**2.0)*sigma 193 scale *= 2 194 195 return extrema

196 197

198 -def TaylorFit(block):

199 ''' 200 a*s*s + b*x*x + c*y*y + d*x*y + e*s + f*x + g*y + h 201 ''' 202 A = [] 203 b = [] 204 for s in range(3): 205 for x in range(3): 206 for y in range(3): 207 z = block[s,x,y] 208 row = [s*s,x*x,y*y,x*y,s,x,y,1] 209 A.append(row) 210 b.append([z]) 211 params,_,_,s = lstsq(A,b) 212 return params.flatten()

213

214 -def TaylorSubpixel(params):

215 ''' 216 a*s*s + b*x*x + c*y*y + d*x*y + e*s + f*x + g*y + h 217 ''' 218 a,b,c,_,e,f,g,_ = params 219 s = -e/(2*a) 220 x = -f/(2*b) 221 y = -g/(2*c) 222 return s,x,y

223

224 -def TaylorApprox(params,s,x,y):

225 ''' 226 a*s*s + b*x*x + c*y*y + d*x*y + e*s + f*x + g*y + h 227 ''' 228 a,b,c,d,e,f,g,h = params 229 return a*s*s + b*x*x + c*y*y + d*x*y + e*s + f*x + g*y + h

230 231 232

233 -class _DetectorDOGTestCase(unittest.TestCase):

234 - def setUp(self):

235 pass

236 237

238 - def testDetectorDOG1(self):

239 detector = DetectorDOG(selector='best',n=100) 240 filename = os.path.join(pv.__path__[0],'data','nonface','NONFACE_1.jpg') 241 im = pv.Image(filename,bw_annotate=True) 242 243 points = detector.detect(im) 244 #print len(points) 245 for _,pt,radius in points: 246 im.annotateCircle(pt,radius) 247 #im.show() 248 self.assertEquals(len(points),100)

249

250 - def testDetectorDOG2(self):

251 detector = DetectorDOG(selector='best') 252 filename = os.path.join(pv.__path__[0],'data','nonface','NONFACE_19.jpg') 253 im = pv.Image(filename,bw_annotate=True) 254 255 points = detector.detect(im) 256 for _,pt,_ in points: 257 im.annotatePoint(pt) 258 259 self.assertEquals(len(points),250)

260

261 - def testDetectorDOG3(self):

262 detector = DetectorDOG() 263 filename = os.path.join(pv.__path__[0],'data','nonface','NONFACE_22.jpg') 264 im = pv.Image(filename,bw_annotate=True) 265 266 points = detector.detect(im) 267 for _,pt,_ in points: 268 im.annotatePoint(pt) 269 270 self.assertEquals(len(points),326)

271

272 - def testDetectorDOG4(self):

273 detector = DetectorDOG() 274 filename = os.path.join(pv.__path__[0],'data','nonface','NONFACE_37.jpg') 275 im = pv.Image(filename,bw_annotate=True) 276 277 points = detector.detect(im) 278 for _,pt,_ in points: 279 im.annotatePoint(pt) 280 281 self.assertEquals(len(points),295)

282

283 - def testDetectorDOG5(self):

284 detector = DetectorDOG(selector='best') 285 filename = os.path.join(pv.__path__[0],'data','nonface','NONFACE_37.jpg') 286 im = pv.Image(filename,bw_annotate=True) 287 288 points = detector.detect(im) 289 for _,pt,_ in points: 290 im.annotatePoint(pt) 291 292 self.assertEquals(len(points),250)

293

294 - def testDetectorDOG6(self):

295 detector = DetectorDOG(selector='all') 296 filename = os.path.join(pv.__path__[0],'data','nonface','NONFACE_37.jpg') 297 im = pv.Image(filename,bw_annotate=True) 298 299 points = detector.detect(im) 300 for _,pt,_ in points: 301 im.annotatePoint(pt) 302 303 self.assertEquals(len(points),561)

304

Source Code for Module pyvision.point.DetectorDOG