00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043 #include "sift.hpp"
00044 #include "sift-conv.tpp"
00045
00046 #include<algorithm>
00047 #include<iostream>
00048 #include<sstream>
00049 #include<cassert>
00050 #include<cstring>
00051
00052 using namespace VL ;
00053
00054
00055 namespace VL {
00056 namespace Detail {
00057
00058 int const expnTableSize = 256 ;
00059 VL::float_t const expnTableMax = VL::float_t(25.0) ;
00060 VL::float_t expnTable [ expnTableSize + 1 ] ;
00061
00062 struct buildExpnTable
00063 {
00064 buildExpnTable() {
00065 for(int k = 0 ; k < expnTableSize + 1 ; ++k) {
00066 expnTable[k] = std::exp( - VL::float_t(k) / expnTableSize * expnTableMax ) ;
00067 }
00068 }
00069 } _buildExpnTable ;
00070
00071 } }
00072
00073
00074 namespace VL {
00075 namespace Detail {
00076
00077
00078 class _cmnt {} cmnt ;
00079
00080
00081
00082
00083
00084
00085
00086
00087 std::istream&
00088 operator>>(std::istream& is, _cmnt& manip)
00089 {
00090 char c ;
00091 char b [1024] ;
00092 is>>c ;
00093 if( c != '#' )
00094 return is.putback(c) ;
00095 is.getline(b,1024) ;
00096 return is ;
00097 }
00098
00099 }
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113 std::ostream&
00114 insertPgm(std::ostream& os, pixel_t const* im, int width, int height)
00115 {
00116 os<< "P5" << "\n"
00117 << width << " "
00118 << height << "\n"
00119 << "255" << "\n" ;
00120 for(int y = 0 ; y < height ; ++y) {
00121 for(int x = 0 ; x < width ; ++x) {
00122 unsigned char v =
00123 (unsigned char)
00124 (std::max(std::min(*im++, 1.0f),0.f) * 255.0f) ;
00125 os << v ;
00126 }
00127 }
00128 return os ;
00129 }
00130
00131
00132 void
00133 createPgmBufferFromArray(int w, int h, pixel_t* d, PgmBuffer& buffer){
00134 buffer.width = w;
00135 buffer.height = h;
00136 buffer.data = new pixel_t[w * h];
00137
00138 for (int i = 0; i < w*h; i++){
00139 buffer.data[i] = d[i];
00140 }
00141 }
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164 std::istream&
00165 extractPgm(std::istream& in, PgmBuffer& buffer)
00166 {
00167 pixel_t* im_pt ;
00168 int width ;
00169 int height ;
00170 int maxval ;
00171
00172 char c ;
00173 in>>c ;
00174 if( c != 'P') VL_THROW("File is not in PGM format") ;
00175
00176 bool is_ascii ;
00177 in>>c ;
00178 switch( c ) {
00179 case '2' : is_ascii = true ; break ;
00180 case '5' : is_ascii = false ; break ;
00181 default : VL_THROW("File is not in PGM format") ;
00182 }
00183
00184 in >> Detail::cmnt
00185 >> width
00186 >> Detail::cmnt
00187 >> height
00188 >> Detail::cmnt
00189 >> maxval ;
00190
00191
00192 {char trash ; in.get(trash) ;}
00193
00194 if(maxval > 255)
00195 VL_THROW("Only <= 8-bit per channel PGM files are supported") ;
00196
00197 if(! in.good())
00198 VL_THROW("PGM header parsing error") ;
00199
00200 im_pt = new pixel_t [ width*height ];
00201
00202 try {
00203 if( is_ascii ) {
00204 pixel_t* start = im_pt ;
00205 pixel_t* end = start + width*height ;
00206 pixel_t norm = pixel_t( maxval ) ;
00207
00208 while( start != end ) {
00209 int i ;
00210 in >> i ;
00211 if( ! in.good() ) VL_THROW
00212 ("PGM parsing error file (width="<<width
00213 <<" height="<<height
00214 <<" maxval="<<maxval
00215 <<" at pixel="<<start-im_pt<<")") ;
00216 *start++ = pixel_t( i ) / norm ;
00217 }
00218 } else {
00219 std::streampos beg = in.tellg() ;
00220 char* newBuffer = new char [width*height] ;
00221 in.read(newBuffer, width*height) ;
00222 if( ! in.good() ) VL_THROW
00223 ("PGM parsing error file (width="<<width
00224 <<" height="<<height
00225 <<" maxval="<<maxval
00226 <<" at pixel="<<in.tellg()-beg<<")") ;
00227
00228 pixel_t* start = im_pt ;
00229 pixel_t* end = start + width*height ;
00230 uint8_t* src = reinterpret_cast<uint8_t*>(newBuffer) ;
00231 while( start != end ) *start++ = *src++ / 255.0f ;
00232 }
00233 } catch(...) {
00234 delete [] im_pt ;
00235 throw ;
00236 }
00237
00238 buffer.width = width ;
00239 buffer.height = height ;
00240 buffer.data = im_pt ;
00241
00242 return in ;
00243 }
00244
00245
00246
00247
00248
00249 namespace Detail {
00250
00251
00252
00253
00254
00255
00256
00257 void
00258 copy(pixel_t* dst, pixel_t const* src, int width, int height)
00259 {
00260 memcpy(dst, src, sizeof(pixel_t)*width*height) ;
00261 }
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273 void
00274 copyAndUpsampleRows
00275 (pixel_t* dst, pixel_t const* src, int width, int height)
00276 {
00277 for(int y = 0 ; y < height ; ++y) {
00278 pixel_t b, a ;
00279 b = a = *src++ ;
00280 for(int x = 0 ; x < width-1 ; ++x) {
00281 b = *src++ ;
00282 *dst = a ; dst += height ;
00283 *dst = (a+b)/2 ; dst += height ;
00284 a = b ;
00285 }
00286 *dst = b ; dst += height ;
00287 *dst = b ; dst += height ;
00288 dst += 1 - width * 2 * height ;
00289 }
00290 }
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305 void
00306 copyAndDownsample(pixel_t* dst, pixel_t const* src,
00307 int width, int height, int d)
00308 {
00309 for(int y = 0 ; y < height ; y+=d) {
00310 pixel_t const * srcrowp = src + y * width ;
00311 for(int x = 0 ; x < width - (d-1) ; x+=d) {
00312 *dst++ = *srcrowp ;
00313 srcrowp += d ;
00314 }
00315 }
00316 }
00317
00318 }
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334 void Sift::smooth (pixel_t* dst, pixel_t* newTemp,
00335 pixel_t const* src, int newWidth, int newHeight,
00336 VL::float_t s) {
00337
00338
00339 int W = int( ceil( VL::float_t(4.0) * s ) ) ;
00340 if( ! filter ) {
00341 filterReserved = 0 ;
00342 }
00343
00344 if( filterReserved < W ) {
00345 filterReserved = W ;
00346 if( filter ) delete [] filter ;
00347 filter = new pixel_t [ 2* filterReserved + 1 ] ;
00348 }
00349
00350
00351 for(int j = 0 ; j < 2*W+1 ; ++j)
00352 filter[j] = VL::pixel_t
00353 (std::exp
00354 (VL::float_t
00355 (-0.5 * (j-W) * (j-W) / (s*s) ))) ;
00356
00357
00358 normalize(filter, W) ;
00359
00360
00361 econvolve(newTemp, src, newWidth, newHeight, filter, W) ;
00362 econvolve(dst, newTemp, newHeight, newWidth, filter, W) ;
00363 }
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382 Sift::Sift(const pixel_t* _im_pt, int _width, int _height,
00383 VL::float_t _sigman,
00384 VL::float_t _sigma0,
00385 int _O, int __S,
00386 int _omin, int _smin, int _smax)
00387 : sigman( _sigman ),
00388 sigma0( _sigma0 ),
00389 sigmak(),
00390 O( _O ),
00391 S( __S ),
00392 omin( _omin ),
00393 smin( _smin ),
00394 smax( _smax ),
00395 width(),
00396 height(),
00397 magnif( 3.0f ),
00398 normalizeDescriptor( true ),
00399
00400 temp( NULL ),
00401 tempReserved(),
00402 tempIsGrad(),
00403 tempOctave(),
00404 octaves( NULL ),
00405 filter( NULL ),
00406 filterReserved(),
00407 keypoints()
00408 {
00409 process(_im_pt, _width, _height) ;
00410 }
00411
00412
00413
00414
00415 Sift::~Sift()
00416 {
00417 freeBuffers() ;
00418 }
00419
00420
00421
00422
00423 void
00424 Sift::
00425 prepareBuffers()
00426 {
00427
00428 int w = (omin >= 0) ? (width >> omin) : (width << -omin) ;
00429 int h = (omin >= 0) ? (height >> omin) : (height << -omin) ;
00430 int size = w*h* std::max
00431 ((smax - smin), 2*((smax+1) - (smin-2) +1)) ;
00432
00433 if( temp && tempReserved == size ) return ;
00434
00435 freeBuffers() ;
00436
00437
00438 temp = new pixel_t [ size ] ;
00439 tempReserved = size ;
00440 tempIsGrad = false ;
00441 tempOctave = 0 ;
00442
00443 octaves = new pixel_t* [ O ] ;
00444 for(int o = 0 ; o < O ; ++o) {
00445 octaves[o] = new pixel_t [ (smax - smin + 1) * w * h ] ;
00446 w >>= 1 ;
00447 h >>= 1 ;
00448 }
00449 }
00450
00451
00452
00453
00454
00455
00456
00457 void
00458 Sift::
00459 freeBuffers()
00460 {
00461 if( filter ) {
00462 delete [] filter ;
00463 }
00464 filter = 0 ;
00465
00466 if( octaves ) {
00467 for(int o = 0 ; o < O ; ++o) {
00468 delete [] octaves[ o ] ;
00469 }
00470 delete [] octaves ;
00471 }
00472 octaves = 0 ;
00473
00474 if( temp ) {
00475 delete [] temp ;
00476 }
00477 temp = 0 ;
00478 }
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498 Sift::Keypoint
00499 Sift::getKeypoint(VL::float_t x, VL::float_t y, VL::float_t sigma) const
00500 {
00501
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535
00536
00537
00538
00539
00540
00541
00542
00543
00544
00545
00546
00547
00548
00549
00550
00551
00552
00553
00554
00555
00556
00557
00558
00559 int o,ix,iy,is ;
00560 VL::float_t s,phi ;
00561
00562 phi = log2f(sigma/sigma0) ;
00563 o = fast_floor( phi - (VL::float_t(smin)+.5f)/S ) ;
00564 o = std::min(o, omin+O-1) ;
00565 o = std::max(o, omin ) ;
00566 s = S * (phi - o) ;
00567
00568 is = int(s + 0.5) ;
00569 is = std::min(is, smax - 2) ;
00570 is = std::max(is, smin + 1) ;
00571
00572 VL::float_t per = getOctaveSamplingPeriod(o) ;
00573 ix = int(x / per + 0.5) ;
00574 iy = int(y / per + 0.5) ;
00575
00576 Keypoint key ;
00577 key.o = o ;
00578
00579 key.ix = ix ;
00580 key.iy = iy ;
00581 key.is = is ;
00582
00583 key.x = x ;
00584 key.y = y ;
00585 key.s = s ;
00586
00587 key.sigma = sigma ;
00588
00589 return key ;
00590 }
00591
00592
00593
00594
00595
00596
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609 void
00610 Sift::
00611 process(const pixel_t* _im_pt, int _width, int _height)
00612 {
00613 using namespace Detail ;
00614
00615 width = _width ;
00616 height = _height ;
00617
00618 prepareBuffers() ;
00619
00620 VL::float_t sigmak1 = powf(2.0f, 1.0f / S) ;
00621 VL::float_t dsigma0 = sigma0 * std::sqrt (1.0f - 1.0f / (sigmak1*sigmak1) ) ;
00622
00623
00624
00625
00626 if( omin < 0 ) {
00627 copyAndUpsampleRows(temp, _im_pt, width, height ) ;
00628 copyAndUpsampleRows(octaves[0], temp, height, 2*width ) ;
00629
00630 for(int o = -1 ; o > omin ; --o) {
00631 copyAndUpsampleRows(temp, octaves[0], width << -o, height << -o) ;
00632 copyAndUpsampleRows(octaves[0], temp, height << -o, 2*(width << -o)) ; }
00633
00634 } else if( omin > 0 ) {
00635 copyAndDownsample(octaves[0], _im_pt, width, height, 1 << omin) ;
00636 } else {
00637 copy(octaves[0], _im_pt, width, height) ;
00638 }
00639
00640 {
00641 VL::float_t sa = sigma0 * powf(sigmak1, smin) ;
00642 VL::float_t sb = sigman / powf(2.0f, omin) ;
00643 if( sa > sb ) {
00644 VL::float_t sd = std::sqrt ( sa*sa - sb*sb ) ;
00645 smooth( octaves[0], temp, octaves[0],
00646 getOctaveWidth(omin),
00647 getOctaveHeight(omin),
00648 sd ) ;
00649 }
00650 }
00651
00652
00653
00654
00655 for(int o = omin ; o < omin+O ; ++o) {
00656
00657 if( o > omin ) {
00658 int sbest = std::min(smin + S, smax) ;
00659 copyAndDownsample(getLevel(o, smin ),
00660 getLevel(o-1, sbest),
00661 getOctaveWidth(o-1),
00662 getOctaveHeight(o-1), 2 ) ;
00663 VL::float_t sa = sigma0 * powf(sigmak1, smin ) ;
00664 VL::float_t sb = sigma0 * powf(sigmak1, sbest - S ) ;
00665 if(sa > sb ) {
00666 VL::float_t sd = std::sqrt ( sa*sa - sb*sb ) ;
00667 smooth( getLevel(o,0), temp, getLevel(o,0),
00668 getOctaveWidth(o), getOctaveHeight(o),
00669 sd ) ;
00670 }
00671 }
00672
00673
00674 for(int s = smin+1 ; s <= smax ; ++s) {
00675 VL::float_t sd = dsigma0 * powf(sigmak1, s) ;
00676 smooth( getLevel(o,s), temp, getLevel(o,s-1),
00677 getOctaveWidth(o), getOctaveHeight(o),
00678 sd ) ;
00679 }
00680 }
00681 }
00682
00683
00684
00685
00686
00687
00688
00689
00690
00691
00692
00693
00694
00695
00696
00697
00698
00699 void
00700 Sift::detectKeypoints(VL::float_t threshold, VL::float_t edgeThreshold)
00701 {
00702 keypoints.clear() ;
00703
00704 size_t nValidatedKeypoints = 0 ;
00705
00706
00707 for(int o = omin ; o < omin + O ; ++o) {
00708
00709 int const xo = 1 ;
00710 int const yo = getOctaveWidth(o) ;
00711 int const so = getOctaveWidth(o) * getOctaveHeight(o) ;
00712 int const ow = getOctaveWidth(o) ;
00713 int const oh = getOctaveHeight(o) ;
00714
00715 VL::float_t xperiod = getOctaveSamplingPeriod(o) ;
00716
00717
00718
00719
00720 pixel_t* dog = temp ;
00721 tempIsGrad = false ;
00722 {
00723 pixel_t* pt = dog ;
00724 for(int s = smin ; s <= smax-1 ; ++s) {
00725 pixel_t* srca = getLevel(o, s ) ;
00726 pixel_t* srcb = getLevel(o, s+1) ;
00727 pixel_t* enda = srcb ;
00728 while( srca != enda ) {
00729 *pt++ = *srcb++ - *srca++ ;
00730 }
00731 }
00732 }
00733
00734
00735
00736
00737 {
00738 pixel_t* pt = dog + xo + yo + so ;
00739 for(int s = smin+1 ; s <= smax-2 ; ++s) {
00740 for(int y = 1 ; y < oh - 1 ; ++y) {
00741 for(int x = 1 ; x < ow - 1 ; ++x) {
00742 pixel_t v = *pt ;
00743
00744
00745
00746 #define CHECK_NEIGHBORS(CMP,SGN) \
00747 ( v CMP ## = SGN 0.8 * threshold && \
00748 v CMP *(pt + xo) && \
00749 v CMP *(pt - xo) && \
00750 v CMP *(pt + so) && \
00751 v CMP *(pt - so) && \
00752 v CMP *(pt + yo) && \
00753 v CMP *(pt - yo) && \
00754 \
00755 v CMP *(pt + yo + xo) && \
00756 v CMP *(pt + yo - xo) && \
00757 v CMP *(pt - yo + xo) && \
00758 v CMP *(pt - yo - xo) && \
00759 \
00760 v CMP *(pt + xo + so) && \
00761 v CMP *(pt - xo + so) && \
00762 v CMP *(pt + yo + so) && \
00763 v CMP *(pt - yo + so) && \
00764 v CMP *(pt + yo + xo + so) && \
00765 v CMP *(pt + yo - xo + so) && \
00766 v CMP *(pt - yo + xo + so) && \
00767 v CMP *(pt - yo - xo + so) && \
00768 \
00769 v CMP *(pt + xo - so) && \
00770 v CMP *(pt - xo - so) && \
00771 v CMP *(pt + yo - so) && \
00772 v CMP *(pt - yo - so) && \
00773 v CMP *(pt + yo + xo - so) && \
00774 v CMP *(pt + yo - xo - so) && \
00775 v CMP *(pt - yo + xo - so) && \
00776 v CMP *(pt - yo - xo - so) )
00777
00778 if( CHECK_NEIGHBORS(>,+) || CHECK_NEIGHBORS(<,-) ) {
00779
00780 Keypoint k ;
00781 k.ix = x ;
00782 k.iy = y ;
00783 k.is = s ;
00784 keypoints.push_back(k) ;
00785 }
00786 pt += 1 ;
00787 }
00788 pt += 2 ;
00789 }
00790 pt += 2*yo ;
00791 }
00792 }
00793
00794
00795
00796
00797 {
00798 KeypointsIter siter ;
00799 KeypointsIter diter ;
00800
00801 for(diter = siter = keypointsBegin() + nValidatedKeypoints ;
00802 siter != keypointsEnd() ;
00803 ++siter) {
00804
00805 int x = int( siter->ix ) ;
00806 int y = int( siter->iy ) ;
00807 int s = int( siter->is ) ;
00808
00809 VL::float_t Dx=0,Dy=0,Ds=0,Dxx=0,Dyy=0,Dss=0,Dxy=0,Dxs=0,Dys=0 ;
00810 VL::float_t b [3] ;
00811 pixel_t* pt ;
00812 int dx = 0 ;
00813 int dy = 0 ;
00814
00815
00816 for(int iter = 0 ; iter < 5 ; ++iter) {
00817
00818 VL::float_t A[3*3] ;
00819
00820 x += dx ;
00821 y += dy ;
00822
00823 pt = dog
00824 + xo * x
00825 + yo * y
00826 + so * (s - smin) ;
00827
00828 #define at(dx,dy,ds) (*( pt + (dx)*xo + (dy)*yo + (ds)*so))
00829 #define Aat(i,j) (A[(i)+(j)*3])
00830
00831
00832 Dx = 0.5f * (at(+1,0,0) - at(-1,0,0)) ;
00833 Dy = 0.5f * (at(0,+1,0) - at(0,-1,0));
00834 Ds = 0.5f * (at(0,0,+1) - at(0,0,-1)) ;
00835
00836
00837 Dxx = (at(+1,0,0) + at(-1,0,0) - 2.0f * at(0,0,0)) ;
00838 Dyy = (at(0,+1,0) + at(0,-1,0) - 2.0f * at(0,0,0)) ;
00839 Dss = (at(0,0,+1) + at(0,0,-1) - 2.0f * at(0,0,0)) ;
00840
00841 Dxy = 0.25f * ( at(+1,+1,0) + at(-1,-1,0) - at(-1,+1,0) - at(+1,-1,0) ) ;
00842 Dxs = 0.25f * ( at(+1,0,+1) + at(-1,0,-1) - at(-1,0,+1) - at(+1,0,-1) ) ;
00843 Dys = 0.25f * ( at(0,+1,+1) + at(0,-1,-1) - at(0,-1,+1) - at(0,+1,-1) ) ;
00844
00845
00846 Aat(0,0) = Dxx ;
00847 Aat(1,1) = Dyy ;
00848 Aat(2,2) = Dss ;
00849 Aat(0,1) = Aat(1,0) = Dxy ;
00850 Aat(0,2) = Aat(2,0) = Dxs ;
00851 Aat(1,2) = Aat(2,1) = Dys ;
00852
00853 b[0] = - Dx ;
00854 b[1] = - Dy ;
00855 b[2] = - Ds ;
00856
00857
00858 for(int j = 0 ; j < 3 ; ++j) {
00859
00860
00861 VL::float_t maxa = 0 ;
00862 VL::float_t maxabsa = 0 ;
00863 int maxi = -1 ;
00864 int i ;
00865 for(i = j ; i < 3 ; ++i) {
00866 VL::float_t a = Aat(i,j) ;
00867 VL::float_t absa = fabsf( a ) ;
00868 if ( absa > maxabsa ) {
00869 maxa = a ;
00870 maxabsa = absa ;
00871 maxi = i ;
00872 }
00873 }
00874
00875
00876 if( maxabsa < 1e-10f ) {
00877 b[0] = 0 ;
00878 b[1] = 0 ;
00879 b[2] = 0 ;
00880 break ;
00881 }
00882
00883 i = maxi ;
00884
00885
00886
00887 for(int jj = j ; jj < 3 ; ++jj) {
00888 std::swap( Aat(j,jj) , Aat(i,jj) ) ;
00889 Aat(j,jj) /= maxa ;
00890 }
00891 std::swap( b[j], b[i] ) ;
00892 b[j] /= maxa ;
00893
00894
00895 for(int ii = j+1 ; ii < 3 ; ++ii) {
00896 VL::float_t x1 = Aat(ii,j) ;
00897 for(int jj = j ; jj < 3 ; ++jj) {
00898 Aat(ii,jj) -= x1 * Aat(j,jj) ;
00899 }
00900 b[ii] -= x1 * b[j] ;
00901 }
00902 }
00903
00904
00905 for(int i = 2 ; i > 0 ; --i) {
00906 VL::float_t x1 = b[i] ;
00907 for(int ii = i-1 ; ii >= 0 ; --ii) {
00908 b[ii] -= x1 * Aat(ii,i) ;
00909 }
00910 }
00911
00912
00913
00914
00915 dx= ((b[0] > 0.6 && x < ow-2) ? 1 : 0 )
00916 + ((b[0] < -0.6 && x > 1 ) ? -1 : 0 ) ;
00917
00918 dy= ((b[1] > 0.6 && y < oh-2) ? 1 : 0 )
00919 + ((b[1] < -0.6 && y > 1 ) ? -1 : 0 ) ;
00920
00921
00922
00923
00924
00925
00926
00927
00928 if( dx == 0 && dy == 0 ) break ;
00929 }
00930
00931
00932
00933
00934 {
00935 VL::float_t val = at(0,0,0) + 0.5f * (Dx * b[0] + Dy * b[1] + Ds * b[2]) ;
00936 VL::float_t score = (Dxx+Dyy)*(Dxx+Dyy) / (Dxx*Dyy - Dxy*Dxy) ;
00937 VL::float_t xn = x + b[0] ;
00938 VL::float_t yn = y + b[1] ;
00939 VL::float_t sn = s + b[2] ;
00940
00941 if(fast_abs(val) > threshold &&
00942 score < (edgeThreshold+1)*(edgeThreshold+1)/edgeThreshold &&
00943 score >= 0 &&
00944 fast_abs(b[0]) < 1.5 &&
00945 fast_abs(b[1]) < 1.5 &&
00946 fast_abs(b[2]) < 1.5 &&
00947 xn >= 0 &&
00948 xn <= ow-1 &&
00949 yn >= 0 &&
00950 yn <= oh-1 &&
00951 sn >= smin &&
00952 sn <= smax ) {
00953
00954 diter->o = o ;
00955
00956 diter->ix = x ;
00957 diter->iy = y ;
00958 diter->is = s ;
00959
00960 diter->x = xn * xperiod ;
00961 diter->y = yn * xperiod ;
00962 diter->s = sn ;
00963
00964 diter->sigma = getScaleFromIndex(o,sn) ;
00965
00966
00967
00968 ++diter ;
00969 }
00970 }
00971 }
00972
00973
00974 keypoints.resize( diter - keypoints.begin() ) ;
00975 nValidatedKeypoints = keypoints.size() ;
00976 }
00977
00978 }
00979 }
00980
00981
00982
00983
00984
00985
00986
00987
00988
00989
00990
00991
00992
00993
00994
00995
00996
00997
00998
00999
01000
01001 void
01002 Sift::prepareGrad(int o)
01003 {
01004 int const ow = getOctaveWidth(o) ;
01005 int const oh = getOctaveHeight(o) ;
01006 int const xo = 1 ;
01007 int const yo = ow ;
01008 int const so = oh*ow ;
01009
01010 if( ! tempIsGrad || tempOctave != o ) {
01011
01012
01013 for(int s = smin+1 ; s <= smax-2 ; ++s) {
01014 for(int y = 1 ; y < oh-1 ; ++y ) {
01015 pixel_t* src = getLevel(o, s) + xo + yo*y ;
01016 pixel_t* end = src + ow - 1 ;
01017 pixel_t* grad = 2 * (xo + yo*y + (s - smin -1)*so) + temp ;
01018 while(src != end) {
01019 VL::float_t Gx = 0.5f * ( *(src+xo) - *(src-xo) ) ;
01020 VL::float_t Gy = 0.5f * ( *(src+yo) - *(src-yo) ) ;
01021 VL::float_t m = fast_sqrt( Gx*Gx + Gy*Gy ) ;
01022 VL::float_t t = fast_mod_2pi( fast_atan2(Gy, Gx) + VL::float_t(2*M_PI) );
01023 *grad++ = pixel_t( m ) ;
01024 *grad++ = pixel_t( t ) ;
01025 ++src ;
01026 }
01027 }
01028 }
01029 }
01030
01031 tempIsGrad = true ;
01032 tempOctave = o ;
01033 }
01034
01035
01036
01037
01038
01039
01040
01041
01042
01043
01044
01045
01046
01047
01048
01049
01050
01051
01052
01053
01054
01055
01056
01057
01058 int
01059 Sift::computeKeypointOrientations(VL::float_t angles [4], Keypoint keypoint)
01060 {
01061 int const nbins = 36 ;
01062 VL::float_t const winFactor = 1.5f ;
01063 VL::float_t hist [nbins] ;
01064
01065
01066 int o = keypoint.o ;
01067 VL::float_t xperiod = getOctaveSamplingPeriod(o) ;
01068
01069
01070 const int ow = getOctaveWidth(o) ;
01071 const int oh = getOctaveHeight(o) ;
01072 const int xo = 2 ;
01073 const int yo = xo * ow ;
01074 const int so = yo * oh ;
01075
01076
01077 VL::float_t x = keypoint.x / xperiod ;
01078 VL::float_t y = keypoint.y / xperiod ;
01079 VL::float_t sigma = keypoint.sigma / xperiod ;
01080
01081
01082 int xi = ((int) (x+0.5)) ;
01083 int yi = ((int) (y+0.5)) ;
01084 int si = keypoint.is ;
01085
01086 VL::float_t const sigmaw = winFactor * sigma ;
01087 int W = (int) floor(3.0 * sigmaw) ;
01088
01089
01090 if(o < omin ||
01091 o >=omin+O ||
01092 xi < 0 ||
01093 xi > ow-1 ||
01094 yi < 0 ||
01095 yi > oh-1 ||
01096 si < smin+1 ||
01097 si > smax-2 ) {
01098 std::cerr<<"!"<<std::endl ;
01099 return 0 ;
01100 }
01101
01102
01103 prepareGrad(o) ;
01104
01105
01106 std::fill(hist, hist + nbins, 0) ;
01107
01108
01109 pixel_t* pt = temp + xi * xo + yi * yo + (si - smin -1) * so ;
01110
01111 #undef at
01112 #define at(dx,dy) (*(pt + (dx)*xo + (dy)*yo))
01113
01114 for(int ys = std::max(-W, 1-yi) ; ys <= std::min(+W, oh -2 -yi) ; ++ys) {
01115 for(int xs = std::max(-W, 1-xi) ; xs <= std::min(+W, ow -2 -xi) ; ++xs) {
01116
01117 VL::float_t dx = xi + xs - x;
01118 VL::float_t dy = yi + ys - y;
01119 VL::float_t r2 = dx*dx + dy*dy ;
01120
01121
01122 if(r2 >= W*W+0.5) continue ;
01123
01124 VL::float_t wgt = VL::fast_expn( r2 / (2*sigmaw*sigmaw) ) ;
01125 VL::float_t mod = *(pt + xs*xo + ys*yo) ;
01126 VL::float_t ang = *(pt + xs*xo + ys*yo + 1) ;
01127
01128
01129 int bin = (int) floor( nbins * ang / (2*M_PI) ) ;
01130 hist[bin] += mod * wgt ;
01131 }
01132 }
01133
01134
01135 #if defined VL_LOWE_STRICT
01136
01137
01138 for (int iter = 0; iter < 6; iter++) {
01139 VL::float_t prev = hist[nbins/2] ;
01140 for (int i = nbins/2-1; i >= -nbins/2 ; --i) {
01141 int const j = (i + nbins) % nbins ;
01142 int const jp = (i - 1 + nbins) % nbins ;
01143 VL::float_t newh = (prev + hist[j] + hist[jp]) / 3.0;
01144 prev = hist[j] ;
01145 hist[j] = newh ;
01146 }
01147 }
01148 #else
01149
01150 for (int iter = 0; iter < 6; iter++) {
01151 VL::float_t prev = hist[nbins-1] ;
01152 VL::float_t first = hist[0] ;
01153 int i ;
01154 for (i = 0; i < nbins - 1; i++) {
01155 VL::float_t newh = (prev + hist[i] + hist[(i+1) % nbins]) / 3.0f;
01156 prev = hist[i] ;
01157 hist[i] = newh ;
01158 }
01159 hist[i] = (prev + hist[i] + first)/3.0f ;
01160 }
01161 #endif
01162
01163
01164 VL::float_t maxh = * std::max_element(hist, hist + nbins) ;
01165
01166
01167 int nangles = 0 ;
01168 for(int i = 0 ; i < nbins ; ++i) {
01169 VL::float_t h0 = hist [i] ;
01170 VL::float_t hm = hist [(i-1+nbins) % nbins] ;
01171 VL::float_t hp = hist [(i+1+nbins) % nbins] ;
01172
01173
01174 if( h0 > 0.8*maxh && h0 > hm && h0 > hp ) {
01175
01176
01177
01178 VL::float_t di = -0.5f * (hp - hm) / (hp+hm-2*h0) ;
01179 VL::float_t th = 2*float(M_PI) * (i+di+0.5f) / nbins ;
01180 angles [ nangles++ ] = th ;
01181 if( nangles == 4 )
01182 goto enough_angles ;
01183 }
01184 }
01185 enough_angles:
01186 return nangles ;
01187 }
01188
01189
01190
01191
01192
01193 namespace Detail {
01194
01195
01196 void
01197 normalize_histogram(VL::float_t* L_begin, VL::float_t* L_end)
01198 {
01199 VL::float_t* L_iter ;
01200 VL::float_t norm = 0 ;
01201
01202 for(L_iter = L_begin; L_iter != L_end ; ++L_iter)
01203 norm += (*L_iter) * (*L_iter) ;
01204
01205 norm = fast_sqrt(norm) ;
01206
01207 for(L_iter = L_begin; L_iter != L_end ; ++L_iter)
01208 *L_iter /= (norm + std::numeric_limits<VL::float_t>::epsilon() ) ;
01209 }
01210
01211 }
01212
01213
01214
01215
01216
01217
01218
01219
01220
01221
01222
01223
01224
01225
01226
01227
01228
01229
01230 void
01231 Sift::computeKeypointDescriptor
01232 (VL::float_t* descr_pt,
01233 Keypoint keypoint,
01234 VL::float_t angle0)
01235 {
01236
01237
01238
01239
01240
01241
01242
01243
01244
01245
01246
01247
01248
01249
01250
01251
01252
01253 int o = keypoint.o ;
01254 VL::float_t xperiod = getOctaveSamplingPeriod(o) ;
01255
01256
01257 const int ow = getOctaveWidth(o) ;
01258 const int oh = getOctaveHeight(o) ;
01259 const int xo = 2 ;
01260 const int yo = xo * ow ;
01261 const int so = yo * oh ;
01262
01263
01264
01265
01266 VL::float_t x = keypoint.x / xperiod;
01267 VL::float_t y = keypoint.y / xperiod ;
01268 VL::float_t sigma = keypoint.sigma / xperiod ;
01269
01270
01271
01272 VL::float_t st0 = sinf( angle0 ) ;
01273 VL::float_t ct0 = cosf( angle0 ) ;
01274
01275
01276 int xi = ((int) (x+0.5)) ;
01277 int yi = ((int) (y+0.5)) ;
01278 int si = keypoint.is ;
01279
01280
01281 const int NBO = 8 ;
01282 const int NBP = 4 ;
01283 const VL::float_t SBP = magnif * sigma ;
01284 const int W = (int) floor (sqrt(2.0) * SBP * (NBP + 1) / 2.0 + 0.5) ;
01285
01286
01287
01288
01289 const int binto = 1 ;
01290 const int binyo = NBO * NBP ;
01291 const int binxo = NBO ;
01292
01293
01294 int bin ;
01295
01296
01297 if(o < omin ||
01298 o >=omin+O ||
01299 xi < 0 ||
01300 xi > ow-1 ||
01301 yi < 0 ||
01302 yi > oh-1 ||
01303 si < smin+1 ||
01304 si > smax-2 )
01305 return ;
01306
01307
01308 prepareGrad(o) ;
01309
01310 std::fill( descr_pt, descr_pt + NBO*NBP*NBP, 0 ) ;
01311
01312
01313
01314
01315 pixel_t const * pt = temp + xi*xo + yi*yo + (si - smin - 1)*so ;
01316 VL::float_t * dpt = descr_pt + (NBP/2) * binyo + (NBP/2) * binxo ;
01317
01318 #define atd(dbinx,dbiny,dbint) *(dpt + (dbint)*binto + (dbiny)*binyo + (dbinx)*binxo)
01319
01320
01321
01322
01323
01324 for(int dyi = std::max(-W, 1-yi) ; dyi <= std::min(+W, oh-2-yi) ; ++dyi) {
01325 for(int dxi = std::max(-W, 1-xi) ; dxi <= std::min(+W, ow-2-xi) ; ++dxi) {
01326
01327
01328 VL::float_t mod = *( pt + dxi*xo + dyi*yo + 0 ) ;
01329 VL::float_t angle = *( pt + dxi*xo + dyi*yo + 1 ) ;
01330 VL::float_t theta = fast_mod_2pi(-angle + angle0) ;
01331
01332
01333
01334
01335 VL::float_t dx = xi + dxi - x;
01336 VL::float_t dy = yi + dyi - y;
01337
01338
01339
01340 VL::float_t nx = ( ct0 * dx + st0 * dy) / SBP ;
01341 VL::float_t ny = (-st0 * dx + ct0 * dy) / SBP ;
01342 VL::float_t nt = NBO * theta / float(2*M_PI) ;
01343
01344
01345
01346
01347
01348
01349 VL::float_t const wsigma = NBP/2 ;
01350 VL::float_t win = VL::fast_expn((nx*nx + ny*ny)/(2.0f * wsigma * wsigma)) ;
01351
01352
01353
01354 int binx = fast_floor( nx - 0.5f ) ;
01355 int biny = fast_floor( ny - 0.5f ) ;
01356 int bint = fast_floor( nt ) ;
01357 VL::float_t rbinx = nx - (binx+0.5f) ;
01358 VL::float_t rbiny = ny - (biny+0.5f) ;
01359 VL::float_t rbint = nt - bint ;
01360 int dbinx ;
01361 int dbiny ;
01362 int dbint ;
01363
01364
01365 for(dbinx = 0 ; dbinx < 2 ; ++dbinx) {
01366 for(dbiny = 0 ; dbiny < 2 ; ++dbiny) {
01367 for(dbint = 0 ; dbint < 2 ; ++dbint) {
01368
01369 if( binx+dbinx >= -(NBP/2) &&
01370 binx+dbinx < (NBP/2) &&
01371 biny+dbiny >= -(NBP/2) &&
01372 biny+dbiny < (NBP/2) ) {
01373 VL::float_t weight = win
01374 * mod
01375 * fast_abs (1 - dbinx - rbinx)
01376 * fast_abs (1 - dbiny - rbiny)
01377 * fast_abs (1 - dbint - rbint) ;
01378
01379
01380
01381 atd(binx+dbinx, biny+dbiny, (bint+dbint) % NBO) += weight ;
01382 }
01383 }
01384 }
01385 }
01386 }
01387 }
01388
01389
01390 if( normalizeDescriptor ) {
01391
01392
01393 Detail::normalize_histogram(descr_pt, descr_pt + NBO*NBP*NBP) ;
01394
01395
01396 for(bin = 0; bin < NBO*NBP*NBP ; ++bin) {
01397 if (descr_pt[bin] > 0.2f) descr_pt[bin] = 0.2f;
01398 }
01399
01400
01401 Detail::normalize_histogram(descr_pt, descr_pt + NBO*NBP*NBP) ;
01402 }
01403
01404 }
01405
01406
01407 }