Site icon R-bloggers

meanShiftR

[This article was first published on MeanMean, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
< svg style="display: none;">< defs id="MathJax_SVG_glyphs">< path stroke-width="1" id="MJMAIN-78" d="M201 0Q189 3 102 3Q26 3 17 0H11V46H25Q48 47 67 52T96 61T121 78T139 96T160 122T180 150L226 210L168 288Q159 301 149 315T133 336T122 351T113 363T107 370T100 376T94 379T88 381T80 383Q74 383 44 385H16V431H23Q59 429 126 429Q219 429 229 431H237V385Q201 381 201 369Q201 367 211 353T239 315T268 274L272 270L297 304Q329 345 329 358Q329 364 327 369T322 376T317 380T310 384L307 385H302V431H309Q324 428 408 428Q487 428 493 431H499V385H492Q443 385 411 368Q394 360 377 341T312 257L296 236L358 151Q424 61 429 57T446 50Q464 46 499 46H516V0H510H502Q494 1 482 1T457 2T432 2T414 3Q403 3 377 3T327 1L304 0H295V46H298Q309 46 320 51T331 63Q331 65 291 120L250 175Q249 174 219 133T185 88Q181 83 181 74Q181 63 188 55T206 46Q208 46 208 23V0H201Z">< path stroke-width="1" id="MJMAIN-A8" d="M95 612Q95 633 112 651T153 669T193 652T210 612Q210 588 194 571T152 554L127 560Q95 577 95 612ZM289 611Q289 634 304 649T335 668Q336 668 340 668T346 669Q369 669 386 652T404 612T387 572T346 554Q323 554 306 570T289 611Z">< path stroke-width="1" id="MJMAIN-2C7" d="M114 611L127 630L136 644Q138 644 193 612Q248 581 250 581L306 612Q361 644 363 644L385 611L318 562L249 513L114 611Z">< path stroke-width="1" id="MJMATHI-67" d="M311 43Q296 30 267 15T206 0Q143 0 105 45T66 160Q66 265 143 353T314 442Q361 442 401 394L404 398Q406 401 409 404T418 412T431 419T447 422Q461 422 470 413T480 394Q480 379 423 152T363 -80Q345 -134 286 -169T151 -205Q10 -205 10 -137Q10 -111 28 -91T74 -71Q89 -71 102 -80T116 -111Q116 -121 114 -130T107 -144T99 -154T92 -162L90 -164H91Q101 -167 151 -167Q189 -167 211 -155Q234 -144 254 -122T282 -75Q288 -56 298 -13Q311 35 311 43ZM384 328L380 339Q377 350 375 354T369 368T359 382T346 393T328 402T306 405Q262 405 221 352Q191 313 171 233T151 117Q151 38 213 38Q269 38 323 108L331 118L384 328Z">< path stroke-width="1" id="MJMAIN-221D" d="M56 124T56 216T107 375T238 442Q260 442 280 438T319 425T352 407T382 385T406 361T427 336T442 315T455 297T462 285L469 297Q555 442 679 442Q687 442 722 437V398H718Q710 400 694 400Q657 400 623 383T567 343T527 294T503 253T495 235Q495 231 520 192T554 143Q625 44 696 44Q717 44 719 46H722V-5Q695 -11 678 -11Q552 -11 457 141Q455 145 454 146L447 134Q362 -11 235 -11Q157 -11 107 56ZM93 213Q93 143 126 87T220 31Q258 31 292 48T349 88T389 137T413 178T421 196Q421 200 396 239T362 288Q322 345 288 366T213 387Q163 387 128 337T93 213Z">< path stroke-width="1" id="MJMATHI-4B" d="M285 628Q285 635 228 637Q205 637 198 638T191 647Q191 649 193 661Q199 681 203 682Q205 683 214 683H219Q260 681 355 681Q389 681 418 681T463 682T483 682Q500 682 500 674Q500 669 497 660Q496 658 496 654T495 648T493 644T490 641T486 639T479 638T470 637T456 637Q416 636 405 634T387 623L306 305Q307 305 490 449T678 597Q692 611 692 620Q692 635 667 637Q651 637 651 648Q651 650 654 662T659 677Q662 682 676 682Q680 682 711 681T791 680Q814 680 839 681T869 682Q889 682 889 672Q889 650 881 642Q878 637 862 637Q787 632 726 586Q710 576 656 534T556 455L509 418L518 396Q527 374 546 329T581 244Q656 67 661 61Q663 59 666 57Q680 47 717 46H738Q744 38 744 37T741 19Q737 6 731 0H720Q680 3 625 3Q503 3 488 0H478Q472 6 472 9T474 27Q478 40 480 43T491 46H494Q544 46 544 71Q544 75 517 141T485 216L427 354L359 301L291 248L268 155Q245 63 245 58Q245 51 253 49T303 46H334Q340 37 340 35Q340 19 333 5Q328 0 317 0Q314 0 280 1T180 2Q118 2 85 2T49 1Q31 1 31 11Q31 13 34 25Q38 41 42 43T65 46Q92 46 125 49Q139 52 144 61Q147 65 216 339T285 628Z">< path stroke-width="1" id="MJMAIN-2032" d="M79 43Q73 43 52 49T30 61Q30 68 85 293T146 528Q161 560 198 560Q218 560 240 545T262 501Q262 496 260 486Q259 479 173 263T84 45T79 43Z">< path stroke-width="1" id="MJMATHI-76" d="M173 380Q173 405 154 405Q130 405 104 376T61 287Q60 286 59 284T58 281T56 279T53 278T49 278T41 278H27Q21 284 21 287Q21 294 29 316T53 368T97 419T160 441Q202 441 225 417T249 361Q249 344 246 335Q246 329 231 291T200 202T182 113Q182 86 187 69Q200 26 250 26Q287 26 319 60T369 139T398 222T409 277Q409 300 401 317T383 343T365 361T357 383Q357 405 376 424T417 443Q436 443 451 425T467 367Q467 340 455 284T418 159T347 40T241 -11Q177 -11 139 22Q102 54 102 117Q102 148 110 181T151 298Q173 362 173 380Z">< path stroke-width="1" id="MJMAIN-28" d="M94 250Q94 319 104 381T127 488T164 576T202 643T244 695T277 729T302 750H315H319Q333 750 333 741Q333 738 316 720T275 667T226 581T184 443T167 250T184 58T225 -81T274 -167T316 -220T333 -241Q333 -250 318 -250H315H302L274 -226Q180 -141 137 -14T94 250Z">< path stroke-width="1" id="MJMAIN-30" d="M96 585Q152 666 249 666Q297 666 345 640T423 548Q460 465 460 320Q460 165 417 83Q397 41 362 16T301 -15T250 -22Q224 -22 198 -16T137 16T82 83Q39 165 39 320Q39 494 96 585ZM321 597Q291 629 250 629Q208 629 178 597Q153 571 145 525T137 333Q137 175 145 125T181 46Q209 16 250 16Q290 16 318 46Q347 76 354 130T362 333Q362 478 354 524T321 597Z">< path stroke-width="1" id="MJMAIN-29" d="M60 749L64 750Q69 750 74 750H86L114 726Q208 641 251 514T294 250Q294 182 284 119T261 12T224 -76T186 -143T145 -194T113 -227T90 -246Q87 -249 86 -250H74Q66 -250 63 -250T58 -247T55 -238Q56 -237 66 -225Q221 -64 221 250T66 725Q56 737 55 738Q55 746 60 749Z">< path stroke-width="1" id="MJMATHI-69" d="M184 600Q184 624 203 642T247 661Q265 661 277 649T290 619Q290 596 270 577T226 557Q211 557 198 567T184 600ZM21 287Q21 295 30 318T54 369T98 420T158 442Q197 442 223 419T250 357Q250 340 236 301T196 196T154 83Q149 61 149 51Q149 26 166 26Q175 26 185 29T208 43T235 78T260 137Q263 149 265 151T282 153Q302 153 302 143Q302 135 293 112T268 61T223 11T161 -11Q129 -11 102 10T74 74Q74 91 79 106T122 220Q160 321 166 341T173 380Q173 404 156 404H154Q124 404 99 371T61 287Q60 286 59 284T58 281T56 279T53 278T49 278T41 278H27Q21 284 21 287Z">< path stroke-width="1" id="MJMAIN-3D" d="M56 347Q56 360 70 367H707Q722 359 722 347Q722 336 708 328L390 327H72Q56 332 56 347ZM56 153Q56 168 72 173H708Q722 163 722 153Q722 140 707 133H70Q56 140 56 153Z">< path stroke-width="1" id="MJSZ1-2211" d="M61 748Q64 750 489 750H913L954 640Q965 609 976 579T993 533T999 516H979L959 517Q936 579 886 621T777 682Q724 700 655 705T436 710H319Q183 710 183 709Q186 706 348 484T511 259Q517 250 513 244L490 216Q466 188 420 134T330 27L149 -187Q149 -188 362 -188Q388 -188 436 -188T506 -189Q679 -189 778 -162T936 -43Q946 -27 959 6H999L913 -249L489 -250Q65 -250 62 -248Q56 -246 56 -239Q56 -234 118 -161Q186 -81 245 -11L428 206Q428 207 242 462L57 717L56 728Q56 744 61 748Z">< path stroke-width="1" id="MJMATHI-4E" d="M234 637Q231 637 226 637Q201 637 196 638T191 649Q191 676 202 682Q204 683 299 683Q376 683 387 683T401 677Q612 181 616 168L670 381Q723 592 723 606Q723 633 659 637Q635 637 635 648Q635 650 637 660Q641 676 643 679T653 683Q656 683 684 682T767 680Q817 680 843 681T873 682Q888 682 888 672Q888 650 880 642Q878 637 858 637Q787 633 769 597L620 7Q618 0 599 0Q585 0 582 2Q579 5 453 305L326 604L261 344Q196 88 196 79Q201 46 268 46H278Q284 41 284 38T282 19Q278 6 272 0H259Q228 2 151 2Q123 2 100 2T63 2T46 1Q31 1 31 10Q31 14 34 26T39 40Q41 46 62 46Q130 49 150 85Q154 91 221 362L289 634Q287 635 234 637Z">< path stroke-width="1" id="MJMATHI-6A" d="M297 596Q297 627 318 644T361 661Q378 661 389 651T403 623Q403 595 384 576T340 557Q322 557 310 567T297 596ZM288 376Q288 405 262 405Q240 405 220 393T185 362T161 325T144 293L137 279Q135 278 121 278H107Q101 284 101 286T105 299Q126 348 164 391T252 441Q253 441 260 441T272 442Q296 441 316 432Q341 418 354 401T367 348V332L318 133Q267 -67 264 -75Q246 -125 194 -164T75 -204Q25 -204 7 -183T-12 -137Q-12 -110 7 -91T53 -71Q70 -71 82 -81T95 -112Q95 -148 63 -167Q69 -168 77 -168Q111 -168 139 -140T182 -74L193 -32Q204 11 219 72T251 197T278 308T289 365Q289 372 288 376Z">< path stroke-width="1" id="MJMAIN-31" d="M213 578L200 573Q186 568 160 563T102 556H83V602H102Q149 604 189 617T245 641T273 663Q275 666 285 666Q294 666 302 660V361L303 61Q310 54 315 52T339 48T401 46H427V0H416Q395 3 257 3Q121 3 100 0H88V46H114Q136 46 152 46T177 47T193 50T201 52T207 57T213 61V578Z">< path stroke-width="1" id="MJMATHI-78" d="M52 289Q59 331 106 386T222 442Q257 442 286 424T329 379Q371 442 430 442Q467 442 494 420T522 361Q522 332 508 314T481 292T458 288Q439 288 427 299T415 328Q415 374 465 391Q454 404 425 404Q412 404 406 402Q368 386 350 336Q290 115 290 78Q290 50 306 38T341 26Q378 26 414 59T463 140Q466 150 469 151T485 153H489Q504 153 504 145Q504 144 502 134Q486 77 440 33T333 -11Q263 -11 227 52Q186 -10 133 -10H127Q78 -10 57 16T35 71Q35 103 54 123T99 143Q142 143 142 101Q142 81 130 66T107 46T94 41L91 40Q91 39 97 36T113 29T132 26Q168 26 194 71Q203 87 217 139T245 247T261 313Q266 340 266 352Q266 380 251 392T217 404Q177 404 142 372T93 290Q91 281 88 280T72 278H58Q52 284 52 289Z">< path stroke-width="1" id="MJMAIN-2212" d="M84 237T84 250T98 270H679Q694 262 694 250T679 230H98Q84 237 84 250Z">< path stroke-width="1" id="MJMATHI-54" d="M40 437Q21 437 21 445Q21 450 37 501T71 602L88 651Q93 669 101 677H569H659Q691 677 697 676T704 667Q704 661 687 553T668 444Q668 437 649 437Q640 437 637 437T631 442L629 445Q629 451 635 490T641 551Q641 586 628 604T573 629Q568 630 515 631Q469 631 457 630T439 622Q438 621 368 343T298 60Q298 48 386 46Q418 46 427 45T436 36Q436 31 433 22Q429 4 424 1L422 0Q419 0 415 0Q410 0 363 1T228 2Q99 2 64 0H49Q43 6 43 9T45 27Q49 40 55 46H83H94Q174 46 189 55Q190 56 191 56Q196 59 201 76T241 233Q258 301 269 344Q339 619 339 625Q339 630 310 630H279Q212 630 191 624Q146 614 121 583T67 467Q60 445 57 441T43 437H40Z">< path stroke-width="1" id="MJMATHI-48" d="M228 637Q194 637 192 641Q191 643 191 649Q191 673 202 682Q204 683 219 683Q260 681 355 681Q389 681 418 681T463 682T483 682Q499 682 499 672Q499 670 497 658Q492 641 487 638H485Q483 638 480 638T473 638T464 637T455 637Q416 636 405 634T387 623Q384 619 355 500Q348 474 340 442T328 395L324 380Q324 378 469 378H614L615 381Q615 384 646 504Q674 619 674 627T617 637Q594 637 587 639T580 648Q580 650 582 660Q586 677 588 679T604 682Q609 682 646 681T740 680Q802 680 835 681T871 682Q888 682 888 672Q888 645 876 638H874Q872 638 869 638T862 638T853 637T844 637Q805 636 794 634T776 623Q773 618 704 340T634 58Q634 51 638 51Q646 48 692 46H723Q729 38 729 37T726 19Q722 6 716 0H701Q664 2 567 2Q533 2 504 2T458 2T437 1Q420 1 420 10Q420 15 423 24Q428 43 433 45Q437 46 448 46H454Q481 46 514 49Q520 50 522 50T528 55T534 64T540 82T547 110T558 153Q565 181 569 198Q602 330 602 331T457 332H312L279 197Q245 63 245 58Q245 51 253 49T303 46H334Q340 38 340 37T337 19Q333 6 327 0H312Q275 2 178 2Q144 2 115 2T69 2T48 1Q31 1 31 10Q31 12 34 24Q39 43 44 45Q48 46 59 46H65Q92 46 125 49Q139 52 144 61Q147 65 216 339T285 628Q285 635 228 637Z">< path stroke-width="1" id="MJMAIN-32" d="M109 429Q82 429 66 447T50 491Q50 562 103 614T235 666Q326 666 387 610T449 465Q449 422 429 383T381 315T301 241Q265 210 201 149L142 93L218 92Q375 92 385 97Q392 99 409 186V189H449V186Q448 183 436 95T421 3V0H50V19V31Q50 38 56 46T86 81Q115 113 136 137Q145 147 170 174T204 211T233 244T261 278T284 308T305 340T320 369T333 401T340 431T343 464Q343 527 309 573T212 619Q179 619 154 602T119 569T109 550Q109 549 114 549Q132 549 151 535T170 489Q170 464 154 447T109 429Z">< path stroke-width="1" id="MJSZ1-28" d="M152 251Q152 646 388 850H416Q422 844 422 841Q422 837 403 816T357 753T302 649T255 482T236 250Q236 124 255 19T301 -147T356 -251T403 -315T422 -340Q422 -343 416 -349H388Q359 -325 332 -296T271 -213T212 -97T170 56T152 251Z">< path stroke-width="1" id="MJSZ1-29" d="M305 251Q305 -145 69 -349H56Q43 -349 39 -347T35 -338Q37 -333 60 -307T108 -239T160 -136T204 27T221 250T204 473T160 636T108 740T60 807T35 839Q35 850 50 850H56H69Q197 743 256 566Q305 425 305 251Z">< path stroke-width="1" id="MJMAIN-2E" d="M78 60Q78 84 95 102T138 120Q162 120 180 104T199 61Q199 36 182 18T139 0T96 17T78 60Z">< path stroke-width="1" id="MJMAIN-2208" d="M84 250Q84 372 166 450T360 539Q361 539 377 539T419 540T469 540H568Q583 532 583 520Q583 511 570 501L466 500Q355 499 329 494Q280 482 242 458T183 409T147 354T129 306T124 272V270H568Q583 262 583 250T568 230H124V228Q124 207 134 177T167 112T231 48T328 7Q355 1 466 0H570Q583 -10 583 -20Q583 -32 568 -40H471Q464 -40 446 -40T417 -41Q262 -41 172 45Q84 127 84 250Z">< path stroke-width="1" id="MJMAIN-7B" d="M434 -231Q434 -244 428 -250H410Q281 -250 230 -184Q225 -177 222 -172T217 -161T213 -148T211 -133T210 -111T209 -84T209 -47T209 0Q209 21 209 53Q208 142 204 153Q203 154 203 155Q189 191 153 211T82 231Q71 231 68 234T65 250T68 266T82 269Q116 269 152 289T203 345Q208 356 208 377T209 529V579Q209 634 215 656T244 698Q270 724 324 740Q361 748 377 749Q379 749 390 749T408 750H428Q434 744 434 732Q434 719 431 716Q429 713 415 713Q362 710 332 689T296 647Q291 634 291 499V417Q291 370 288 353T271 314Q240 271 184 255L170 250L184 245Q202 239 220 230T262 196T290 137Q291 131 291 1Q291 -134 296 -147Q306 -174 339 -192T415 -213Q429 -213 431 -216Q434 -219 434 -231Z">< path stroke-width="1" id="MJMAIN-2C" d="M78 35T78 60T94 103T137 121Q165 121 187 96T210 8Q210 -27 201 -60T180 -117T154 -158T130 -185T117 -194Q113 -194 104 -185T95 -172Q95 -168 106 -156T131 -126T157 -76T173 -3V9L172 8Q170 7 167 6T161 3T152 1T140 0Q113 0 96 17Z">< path stroke-width="1" id="MJMAIN-7D" d="M65 731Q65 745 68 747T88 750Q171 750 216 725T279 670Q288 649 289 635T291 501Q292 362 293 357Q306 312 345 291T417 269Q428 269 431 266T434 250T431 234T417 231Q380 231 345 210T298 157Q293 143 292 121T291 -28V-79Q291 -134 285 -156T256 -198Q202 -250 89 -250Q71 -250 68 -247T65 -230Q65 -224 65 -223T66 -218T69 -214T77 -213Q91 -213 108 -210T146 -200T183 -177T207 -139Q208 -134 209 3L210 139Q223 196 280 230Q315 247 330 250Q305 257 280 270Q225 304 212 352L210 362L209 498Q208 635 207 640Q195 680 154 696T77 713Q68 713 67 716T65 731Z">< path stroke-width="1" id="MJMAIN-7C" d="M139 -249H137Q125 -249 119 -235V251L120 737Q130 750 139 750Q152 750 159 735V-235Q151 -249 141 -249H139Z">< path stroke-width="1" id="MJMAIN-2F" d="M423 750Q432 750 438 744T444 730Q444 725 271 248T92 -240Q85 -250 75 -250Q68 -250 62 -245T56 -231Q56 -221 230 257T407 740Q411 750 423 750Z">< path stroke-width="1" id="MJMAIN-3C" d="M694 -11T694 -19T688 -33T678 -40Q671 -40 524 29T234 166L90 235Q83 240 83 250Q83 261 91 266Q664 540 678 540Q681 540 687 534T694 519T687 505Q686 504 417 376L151 250L417 124Q686 -4 687 -5Q694 -11 694 -19Z">< path stroke-width="1" id="MJMATHI-3B4" d="M195 609Q195 656 227 686T302 717Q319 716 351 709T407 697T433 690Q451 682 451 662Q451 644 438 628T403 612Q382 612 348 641T288 671T249 657T235 628Q235 584 334 463Q401 379 401 292Q401 169 340 80T205 -10H198Q127 -10 83 36T36 153Q36 286 151 382Q191 413 252 434Q252 435 245 449T230 481T214 521T201 566T195 609ZM112 130Q112 83 136 55T204 27Q233 27 256 51T291 111T309 178T316 232Q316 267 309 298T295 344T269 400L259 396Q215 381 183 342T137 256T118 179T112 130Z">

In this blog post, I will be introducing the meanShiftR package. meanShiftR is a rewrite of my original mean shift R package from 2013, based on the Fast Library for Approximate Nearest Neighbors (FLANN). The meanShiftR package is focused on providing to R users the most computationally efficient mean shift implementations available in the literature. This includes approximations to the mean shift algorithm through kernel truncations and approximate nearest-neighbor (ANN) approaches.

< !-- algorithm -->

The mean shift algorithm is a steepest ascent classification algorithm, where classification is performed by fixed point iteration to a local maxima of a kernel density estimate. This method is originally credited to (Fukunaga and Hostetler, 1975), but didn’t see wide-scale adoption until it was popularized by (Cheng, 1995). This algorithm has a fairly simple form of a convex combination of the support of the kernel density estimate, where the weight is calculated from the function < svg xmlns:xlink="http://www.w3.org/1999/xlink" width="6.993ex" height="2.676ex" style="vertical-align: -0.671ex;" viewBox="0 -863.1 3010.9 1152.1" role="img" focusable="false">< g stroke="currentColor" fill="currentColor" stroke-width="0" transform="matrix(1 0 0 -1 0 0)">< use x="0" y="0" xlink:href="#MJMATHI-67">< use x="758" y="0" xlink:href="#MJMAIN-221D">< g transform="translate(1814,0)">< use x="0" y="0" xlink:href="#MJMATHI-4B">< use transform="scale(0.707)" x="1274" y="513" xlink:href="#MJMAIN-2032">, the derivative of kernel < svg xmlns:xlink="http://www.w3.org/1999/xlink" width="2.066ex" height="2.176ex" style="vertical-align: -0.338ex;" viewBox="0 -791.3 889.5 936.9" role="img" focusable="false">< g stroke="currentColor" fill="currentColor" stroke-width="0" transform="matrix(1 0 0 -1 0 0)">< use x="0" y="0" xlink:href="#MJMATHI-4B">. In the computer science literature this function < svg xmlns:xlink="http://www.w3.org/1999/xlink" width="1.116ex" height="2.009ex" style="vertical-align: -0.671ex;" viewBox="0 -576.1 480.5 865.1" role="img" focusable="false">< g stroke="currentColor" fill="currentColor" stroke-width="0" transform="matrix(1 0 0 -1 0 0)">< use x="0" y="0" xlink:href="#MJMATHI-67"> is commonly refered to as the profile and < svg xmlns:xlink="http://www.w3.org/1999/xlink" width="2.066ex" height="2.176ex" style="vertical-align: -0.338ex;" viewBox="0 -791.3 889.5 936.9" role="img" focusable="false">< g stroke="currentColor" fill="currentColor" stroke-width="0" transform="matrix(1 0 0 -1 0 0)">< use x="0" y="0" xlink:href="#MJMATHI-4B"> is refered to as the shadow kernel.

The explicit formulation of the mean shift for < svg xmlns:xlink="http://www.w3.org/1999/xlink" width="3.461ex" height="2.676ex" style="vertical-align: -0.338ex;" viewBox="0 -1006.6 1490.2 1152.1" role="img" focusable="false">< g stroke="currentColor" fill="currentColor" stroke-width="0" transform="matrix(1 0 0 -1 0 0)">< use x="0" y="0" xlink:href="#MJMATHI-76">< g transform="translate(485,362)">< use transform="scale(0.707)" x="0" y="0" xlink:href="#MJMAIN-28">< use transform="scale(0.707)" x="389" y="0" xlink:href="#MJMAIN-30">< use transform="scale(0.707)" x="890" y="0" xlink:href="#MJMAIN-29"> (the point being classified) is,

< svg xmlns:xlink="http://www.w3.org/1999/xlink" width="49.73ex" height="8.176ex" style="vertical-align: -3.505ex;" viewBox="0 -2011.3 21411.3 3520.2" role="img" focusable="false">< g stroke="currentColor" fill="currentColor" stroke-width="0" transform="matrix(1 0 0 -1 0 0)">< use x="0" y="0" xlink:href="#MJMATHI-76">< g transform="translate(485,412)">< use transform="scale(0.707)" x="0" y="0" xlink:href="#MJMAIN-28">< use transform="scale(0.707)" x="389" y="0" xlink:href="#MJMATHI-69">< use transform="scale(0.707)" x="734" y="0" xlink:href="#MJMAIN-29">< use x="1658" y="0" xlink:href="#MJMAIN-3D">< g transform="translate(2436,0)">< g transform="translate(397,0)">< rect stroke="none" width="450" x="0" y="220">< g transform="translate(60,950)">< use x="0" y="0" xlink:href="#MJSZ1-2211">< use transform="scale(0.707)" x="1494" y="675" xlink:href="#MJMATHI-4E">< g transform="translate(1056,-287)">< use transform="scale(0.707)" x="0" y="0" xlink:href="#MJMATHI-6A">< use transform="scale(0.707)" x="412" y="0" xlink:href="#MJMAIN-3D">< use transform="scale(0.707)" x="1191" y="0" xlink:href="#MJMAIN-31">< g transform="translate(2519,0)">< use x="0" y="0" xlink:href="#MJMATHI-78">< use transform="scale(0.707)" x="809" y="-213" xlink:href="#MJMATHI-6A">< use x="3483" y="0" xlink:href="#MJMATHI-67">< g transform="translate(4130,0)">< use xlink:href="#MJSZ1-28">< use x="458" y="0" xlink:href="#MJMAIN-28">< g transform="translate(848,0)">< use x="0" y="0" xlink:href="#MJMATHI-76">< g transform="translate(485,362)">< use transform="scale(0.707)" x="0" y="0" xlink:href="#MJMAIN-28">< use transform="scale(0.707)" x="389" y="0" xlink:href="#MJMATHI-69">< use transform="scale(0.707)" x="734" y="0" xlink:href="#MJMAIN-2212">< use transform="scale(0.707)" x="1513" y="0" xlink:href="#MJMAIN-31">< use transform="scale(0.707)" x="2014" y="0" xlink:href="#MJMAIN-29">< use x="3355" y="0" xlink:href="#MJMAIN-2212">< g transform="translate(4355,0)">< use x="0" y="0" xlink:href="#MJMATHI-78">< use transform="scale(0.707)" x="809" y="-213" xlink:href="#MJMATHI-6A">< g transform="translate(5320,0)">< use x="0" y="0" xlink:href="#MJMAIN-29">< use transform="scale(0.707)" x="550" y="513" xlink:href="#MJMATHI-54">< g transform="translate(6307,0)">< use x="0" y="0" xlink:href="#MJMATHI-48">< g transform="translate(905,362)">< use transform="scale(0.707)" x="0" y="0" xlink:href="#MJMAIN-2212">< use transform="scale(0.707)" x="778" y="0" xlink:href="#MJMAIN-32">< use x="8217" y="0" xlink:href="#MJMAIN-28">< g transform="translate(8607,0)">< use x="0" y="0" xlink:href="#MJMATHI-76">< g transform="translate(485,362)">< use transform="scale(0.707)" x="0" y="0" xlink:href="#MJMAIN-28">< use transform="scale(0.707)" x="389" y="0" xlink:href="#MJMATHI-69">< use transform="scale(0.707)" x="734" y="0" xlink:href="#MJMAIN-2212">< use transform="scale(0.707)" x="1513" y="0" xlink:href="#MJMAIN-31">< use transform="scale(0.707)" x="2014" y="0" xlink:href="#MJMAIN-29">< use x="11114" y="0" xlink:href="#MJMAIN-2212">< g transform="translate(12115,0)">< use x="0" y="0" xlink:href="#MJMATHI-78">< use transform="scale(0.707)" x="809" y="-213" xlink:href="#MJMATHI-6A">< use x="13079" y="0" xlink:href="#MJMAIN-29">< use x="13469" y="-1" xlink:href="#MJSZ1-29">< g transform="translate(542,-982)">< use x="0" y="0" xlink:href="#MJSZ1-2211">< use transform="scale(0.707)" x="1494" y="675" xlink:href="#MJMATHI-4E">< g transform="translate(1056,-287)">< use transform="scale(0.707)" x="0" y="0" xlink:href="#MJMATHI-6A">< use transform="scale(0.707)" x="412" y="0" xlink:href="#MJMAIN-3D">< use transform="scale(0.707)" x="1191" y="0" xlink:href="#MJMAIN-31">< use x="2519" y="0" xlink:href="#MJMATHI-67">< g transform="translate(3166,0)">< use xlink:href="#MJSZ1-28">< use x="458" y="0" xlink:href="#MJMAIN-28">< g transform="translate(848,0)">< use x="0" y="0" xlink:href="#MJMATHI-76">< g transform="translate(485,288)">< use transform="scale(0.707)" x="0" y="0" xlink:href="#MJMAIN-28">< use transform="scale(0.707)" x="389" y="0" xlink:href="#MJMATHI-69">< use transform="scale(0.707)" x="734" y="0" xlink:href="#MJMAIN-2212">< use transform="scale(0.707)" x="1513" y="0" xlink:href="#MJMAIN-31">< use transform="scale(0.707)" x="2014" y="0" xlink:href="#MJMAIN-29">< use x="3355" y="0" xlink:href="#MJMAIN-2212">< g transform="translate(4355,0)">< use x="0" y="0" xlink:href="#MJMATHI-78">< use transform="scale(0.707)" x="809" y="-213" xlink:href="#MJMATHI-6A">< g transform="translate(5320,0)">< use x="0" y="0" xlink:href="#MJMAIN-29">< use transform="scale(0.707)" x="550" y="408" xlink:href="#MJMATHI-54">< g transform="translate(6307,0)">< use x="0" y="0" xlink:href="#MJMATHI-48">< g transform="translate(905,288)">< use transform="scale(0.707)" x="0" y="0" xlink:href="#MJMAIN-2212">< use transform="scale(0.707)" x="778" y="0" xlink:href="#MJMAIN-32">< use x="8217" y="0" xlink:href="#MJMAIN-28">< g transform="translate(8607,0)">< use x="0" y="0" xlink:href="#MJMATHI-76">< g transform="translate(485,288)">< use transform="scale(0.707)" x="0" y="0" xlink:href="#MJMAIN-28">< use transform="scale(0.707)" x="389" y="0" xlink:href="#MJMATHI-69">< use transform="scale(0.707)" x="734" y="0" xlink:href="#MJMAIN-2212">< use transform="scale(0.707)" x="1513" y="0" xlink:href="#MJMAIN-31">< use transform="scale(0.707)" x="2014" y="0" xlink:href="#MJMAIN-29">< use x="11114" y="0" xlink:href="#MJMAIN-2212">< g transform="translate(12115,0)">< use x="0" y="0" xlink:href="#MJMATHI-78">< use transform="scale(0.707)" x="809" y="-213" xlink:href="#MJMATHI-6A">< use x="13079" y="0" xlink:href="#MJMAIN-29">< use x="13469" y="-1" xlink:href="#MJSZ1-29">< use x="21132" y="0" xlink:href="#MJMAIN-2E">
Where, the support is defined by the points in < svg xmlns:xlink="http://www.w3.org/1999/xlink" width="15.691ex" height="2.843ex" style="vertical-align: -0.838ex;" viewBox="0 -863.1 6755.7 1223.9" role="img" focusable="false">< g stroke="currentColor" fill="currentColor" stroke-width="0" transform="matrix(1 0 0 -1 0 0)">< use x="0" y="0" xlink:href="#MJMATHI-78">< use transform="scale(0.707)" x="809" y="-213" xlink:href="#MJMATHI-69">< use x="1194" y="0" xlink:href="#MJMAIN-2208">< use x="2139" y="0" xlink:href="#MJMAIN-7B">< use x="2640" y="0" xlink:href="#MJMAIN-31">< use x="3140" y="0" xlink:href="#MJMAIN-2C">< use x="3586" y="0" xlink:href="#MJMAIN-2E">< use x="4031" y="0" xlink:href="#MJMAIN-2E">< use x="4476" y="0" xlink:href="#MJMAIN-2E">< use x="4921" y="0" xlink:href="#MJMAIN-2C">< use x="5366" y="0" xlink:href="#MJMATHI-4E">< use x="6255" y="0" xlink:href="#MJMAIN-7D"> and kernel bandwidth paramter is identified as < svg xmlns:xlink="http://www.w3.org/1999/xlink" width="2.064ex" height="2.176ex" style="vertical-align: -0.338ex;" viewBox="0 -791.3 888.5 936.9" role="img" focusable="false">< g stroke="currentColor" fill="currentColor" stroke-width="0" transform="matrix(1 0 0 -1 0 0)">< use x="0" y="0" xlink:href="#MJMATHI-48">. < svg xmlns:xlink="http://www.w3.org/1999/xlink" width="2.064ex" height="2.176ex" style="vertical-align: -0.338ex;" viewBox="0 -791.3 888.5 936.9" role="img" focusable="false">< g stroke="currentColor" fill="currentColor" stroke-width="0" transform="matrix(1 0 0 -1 0 0)">< use x="0" y="0" xlink:href="#MJMATHI-48"> is assumed to be of full rank and symmetric in this formula.

This algorithm identifies local maxima by updating < svg xmlns:xlink="http://www.w3.org/1999/xlink" width="3.739ex" height="2.843ex" style="vertical-align: -0.838ex;" viewBox="0 -863.1 1610 1223.9" role="img" focusable="false">< g stroke="currentColor" fill="currentColor" stroke-width="0" transform="matrix(1 0 0 -1 0 0)">< use x="0" y="0" xlink:href="#MJMATHI-76">< g transform="translate(485,0)">< use x="0" y="0" xlink:href="#MJMAIN-28">< use x="389" y="0" xlink:href="#MJMATHI-69">< use x="735" y="0" xlink:href="#MJMAIN-29"> at each iteration, starting with a set of initial points. Iteration continues until a fixed number of iterations is met or < svg xmlns:xlink="http://www.w3.org/1999/xlink" width="27.146ex" height="3.176ex" style="vertical-align: -0.838ex;" viewBox="0 -1006.6 11687.7 1367.4" role="img" focusable="false">< g stroke="currentColor" fill="currentColor" stroke-width="0" transform="matrix(1 0 0 -1 0 0)">< use x="0" y="0" xlink:href="#MJMAIN-7C">< use x="278" y="0" xlink:href="#MJMAIN-7C">< g transform="translate(557,0)">< use x="0" y="0" xlink:href="#MJMATHI-76">< g transform="translate(485,362)">< use transform="scale(0.707)" x="0" y="0" xlink:href="#MJMAIN-28">< use transform="scale(0.707)" x="389" y="0" xlink:href="#MJMATHI-69">< use transform="scale(0.707)" x="734" y="0" xlink:href="#MJMAIN-29">< use x="2159" y="0" xlink:href="#MJMAIN-2212">< g transform="translate(3160,0)">< use x="0" y="0" xlink:href="#MJMATHI-76">< g transform="translate(485,362)">< use transform="scale(0.707)" x="0" y="0" xlink:href="#MJMAIN-28">< use transform="scale(0.707)" x="389" y="0" xlink:href="#MJMATHI-69">< use transform="scale(0.707)" x="734" y="0" xlink:href="#MJMAIN-2212">< use transform="scale(0.707)" x="1513" y="0" xlink:href="#MJMAIN-31">< use transform="scale(0.707)" x="2014" y="0" xlink:href="#MJMAIN-29">< use x="5445" y="0" xlink:href="#MJMAIN-7C">< use x="5724" y="0" xlink:href="#MJMAIN-7C">< use x="6002" y="0" xlink:href="#MJMAIN-2F">< use x="6503" y="0" xlink:href="#MJMAIN-7C">< use x="6781" y="0" xlink:href="#MJMAIN-7C">< g transform="translate(7060,0)">< use x="0" y="0" xlink:href="#MJMATHI-76">< g transform="translate(485,362)">< use transform="scale(0.707)" x="0" y="0" xlink:href="#MJMAIN-28">< use transform="scale(0.707)" x="389" y="0" xlink:href="#MJMATHI-69">< use transform="scale(0.707)" x="734" y="0" xlink:href="#MJMAIN-2212">< use transform="scale(0.707)" x="1513" y="0" xlink:href="#MJMAIN-31">< use transform="scale(0.707)" x="2014" y="0" xlink:href="#MJMAIN-29">< use x="9345" y="0" xlink:href="#MJMAIN-7C">< use x="9623" y="0" xlink:href="#MJMAIN-7C">< use x="10179" y="0" xlink:href="#MJMAIN-3C">< use x="11236" y="0" xlink:href="#MJMATHI-3B4"> where < svg xmlns:xlink="http://www.w3.org/1999/xlink" width="1.049ex" height="2.343ex" style="vertical-align: -0.338ex;" viewBox="0 -863.1 451.5 1008.6" role="img" focusable="false">< g stroke="currentColor" fill="currentColor" stroke-width="0" transform="matrix(1 0 0 -1 0 0)">< use x="0" y="0" xlink:href="#MJMATHI-3B4"> is an acceptable tolerance.

< !-- implementations -->

R Implementations

On CRAN there are two other packages that perform the mean shift algorithm, MeanShift and LPCM. The MeanShift package provides multi-core acceleration, a large number of kernels, and support for blured mean shift. LPCM provides a fast and effective mean shift implementation, as well as plotting and diagnostic tools. LPCM is restricted to just Gaussian kernels and does not provide any multi-core acceleration.

The meanShiftR package presented here is considerably faster than the prior two packages, but currently lacks diagnostic tools and support for more kernels. Instead, the meanShiftR package provides a number of ways to speed mean shift. This includes a generalized Gaussian variant of the mean shift (Lisic, 2015) providing accelerated convergence at the expense of some numerical stability, and support for support truncation with nearest neighbor and approximate nearest neighbor searches. Furthermore, this mean shift implementaiton allows for mean shifting on a separate data set than the kernel support, a feature missing in the other two packages.

A quick speed comparison was done between the three packages, this was performed on a 13″ MacBook Pro with a dual-core hyper-threaded CPU. In this configuration the default number of cores used by MeanShiftR is four. A slight change had to be done to the MeanShift package to provide similar results, mainly the terminal condition had to be modified to be equivalent to meanShiftR and LPCM. This terminal condition is simply the relative L< svg xmlns:xlink="http://www.w3.org/1999/xlink" width="1.054ex" height="1.676ex" style="vertical-align: -0.671ex;" viewBox="0 -432.6 453.9 721.6" role="img" focusable="false">< g stroke="currentColor" fill="currentColor" stroke-width="0" transform="matrix(1 0 0 -1 0 0)">< use transform="scale(0.707)" x="0" y="-213" xlink:href="#MJMAIN-32"> distance.

library(meanShiftR)
library(LPCM)
library(MeanShift)

# set a seed to make this reproducible 
set.seed(100)

# set the number of iterations to test 
# (we will terminate well before this)
iter <- 1000

# set the number of points to simulate
n <- 500

# set the bandwidth
h <- c(0.5,0.5)

# create example data
x1 <- matrix( rnorm( n ),ncol=2)
x2 <- matrix( rnorm( n ),ncol=2) + 2
x <- rbind( x1, x2 )

#plot initial points
plot(x, col=rep(c('red','green'),each=n/2),
     cex=2, xlab='x',ylab='y',pch=20)

########### meanShiftR ###################
run.time <- proc.time()
result <- meanShift(
  x,
  x,
  algorithm="KDTREE",
  bandwidth=h,
  alpha=0,
  iterations = iter,
  parameters=c(10,100)
)
meanShiftR_kd_runtime <- (proc.time()-run.time)[3]

# assignment
meanShiftR_kd_assignment <- result$assignment

# value
meanShiftR_kd_value <- result$value


########### meanShiftR ###################
run.time <- proc.time()
result <- meanShift(
  x,
  x,
  bandwidth=h,
  alpha=0,
  iterations = iter
)
meanShiftR_runtime <- (proc.time()-run.time)[3]

# assignment
meanShiftR_assignment <- result$assignment

# value
meanShiftR_value <- result$value


########### LPCM ###################
runtime <- proc.time()
result <- ms(
            x,
            h=h,
            scaled=FALSE,
            iter=iter,
            plotms=-1)
LPCM_runtime <- (proc.time()-runtime)[3]

# assignment
LPCM_assignment <- result$cluster.label

# value
LPCM_value <- result$cluster.center[LPCM_assignment,]


########### MeanShift ###################
options(mc.cores = 4)
z <- t(x)
runtime <- proc.time()
result <- msClustering(
            X=z,
            h=h,
            kernel="gaussianKernel",
            tol.stop=1e-08,
            tol.epsilon=1e-04,
            multi.core=T)
MeanShift_runtime <- (proc.time()-runtime)[3]

MeanShift_assignment <- result$labels
MeanShift_value <- t(result$components[,result$labels])

# print 
plot(x, col=sapply(meanShiftR_assignment,function(x)c('red','green','blue')[x]),
   cex=1.5, xlab='x',ylab='y',pch=20)

result <- data.frame(
  runtime=c( meanShiftR_runtime,
                meanShiftR_kd_runtime,
                LPCM_runtime,
                MeanShift_runtime),
  maxDiff=c(max(abs(meanShiftR_value - LPCM_value)),
               max(abs(meanShiftR_kd_value - LPCM_value)),
               0,
               max(abs(MeanShift_value - LPCM_value))
               ),
  assignmentDiff=c(sum(meanShiftR_assignment != LPCM_assignment),
               sum(meanShiftR_kd_assignment != LPCM_assignment),
               0,
               sum(MeanShift_assignment != LPCM_assignment)
               )
  )

colnames(result) <- c('Run-Time',
                      'Maximum Absolute Difference',
                      'Label Disagreements')

rownames(result) <- c('meanShiftR',
                      'meanShiftR K-D Tree',
                      'LPCM ms',
                      'meanShift msClustering')

library(xtable)
print(xtable(result,digits=6,display=c('s','f','f','d')), type='html')
< !-- html table generated in R 3.3.1 by xtable 1.8-2 package -->

< !-- Mon Sep 19 10:21:51 2016 -->

Run-Time Maximum Absolute Difference Label Disagreements
meanShiftR 0.304000 0.000000 0
meanShiftR K-D Tree 4.871000 0.000000 0
LPCM ms 2.425000 0.000000 0
meanShift msClustering 43.168000 0.000000 0
Figure 1: Original image.
Figure 2: Classified image.

Here we can see some pretty good performance results for meanShiftR relative to the other R packages. Unfortunately the k-d tree based implementation wasn’t as fast as I would have liked, but its utility is really in tuncated support kernels. Support for these kernels will be added later.

It is unclear why the performance of MeanShift was so poor relative to the other approaches. To verify that it was not a function of the parallel implementation, the algorithm was run with one thread achieving approximately 50% of the performance with mc.cores=4.

Diagnostic plots are not yet included, but convergence path plots can be easily generated by looping over single iterations. An example is provided below, with results in Figure 3.

library(meanShiftR)
library(LPCM)
library(MeanShift)

# set a seed to make this reproducible 
set.seed(100)

# set the number of iterations to test 
# (we will terminate well before this)
iter <- 10

# set the number of points to simulate
n <- 500
m <- 20

# set the bandwidth
h <- c(0.5,0.5)

# create example data
x1 <- matrix( rnorm( n ),ncol=2)
x2 <- matrix( rnorm( n ),ncol=2) + 2
x <- rbind( x1, x2 )

# create some test data for diagnostic plots
y1 <- matrix( rnorm( m ) ,ncol=2)
y2 <- matrix( rnorm( m ),ncol=2) + 2
y <- rbind( y1, y2 )

plot(x, col=rep(c('salmon','greenyellow'),each=n/2),
   cex=1.5, xlab='x',ylab='y',pch=20)

points(y,col=rep(c('red','green'),each=m/2),
       cex=2,pch=19)

points(y,cex=2)

########### meanShiftR ###################

#initial run
result <- meanShift(
    y,
    x,
    algorithm="KDTREE",
    bandwidth=h,
    alpha=0,
    iterations = iter,
    parameters=c(10,100)
)

y0 <- rbind(y,result$value)

for( i in 2:iter) {

  result <- meanShift(
    result$value,
    x,
    algorithm="KDTREE",
    bandwidth=h,
    alpha=0,
    iterations = 1,
    parameters=c(10,100)
  )

  # concate on the result
  y0 <- rbind(y0,result$value)
}

# plot the paths

for( i in 1:m ) {
  pointIndex <- seq(from=0,to=(m*(iter-1)),by=m)+i
  points(y0[pointIndex,] , type='l',lwd=2)
}

Figure 3: Convergence paths.
< !-- plan -->

Future Plans

Due to issues with stack size limits in .C, I was unable to directly submit the original FLANN based package to CRAN, therefore a complete rewrite in C has been initiated. This 0.50 release, is the first release on-the-way to 1.00, where feature parity will be met with the original code. The release map for this package follows four more planned releases:

References

Cheng, Y. (1995). Mean shift, mode seeking, and clustering. IEEE transactions on pattern analysis and machine intelligence, 17(8), 790-799.

Fukunaga, K., & Hostetler, L. (1975). The estimation of the gradient of a density function, with applications in pattern recognition. IEEE Transactions on information theory, 21(1), 32-40.

Lisic, J. (2015). Parcel Level Agricultural Land Cover Prediction (Doctoral dissertation, George Mason University).

Wang, P., Lee, D., Gray, A. G., & Rehg, J. M. (2007, March). Fast Mean Shift with Accurate and Stable Convergence. In AISTATS (pp. 604-611).

Xiao, C., & Liu, M. (2010, September). Efficient Mean‐shift Clustering Using Gaussian KD‐Tree. In Computer Graphics Forum (Vol. 29, No. 7, pp. 2065-2073). Blackwell Publishing Ltd.

To leave a comment for the author, please follow the link and comment on their blog: MeanMean.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.