BVB Source Codes

jsfeat Show jsfeat.js Source code

Return Download jsfeat: download jsfeat.js Source code - Download jsfeat Source code - Type:.js
  1. /**
  2.  * @author Eugene Zatepyakin / http://inspirit.ru/
  3.  */
  4.  
  5. // namespace ?
  6. var jsfeat = jsfeat || { REVISION: 'ALPHA' };
  7. /**
  8.  * @author Eugene Zatepyakin / http://inspirit.ru/
  9.  */
  10.  
  11. (function(global) {
  12.     "use strict";
  13.     //
  14.  
  15.     // CONSTANTS
  16.     var EPSILON = 0.0000001192092896;
  17.     var FLT_MIN = 1E-37;
  18.  
  19.     // implementation from CCV project
  20.     // currently working only with u8,s32,f32
  21.     var U8_t = 0x0100,
  22.         S32_t = 0x0200,
  23.         F32_t = 0x0400,
  24.         S64_t = 0x0800,
  25.         F64_t = 0x1000;
  26.  
  27.     var C1_t = 0x01,
  28.         C2_t = 0x02,
  29.         C3_t = 0x03,
  30.         C4_t = 0x04;
  31.  
  32.     var _data_type_size = new Int32Array([ -1, 1, 4, -1, 4, -1, -1, -1, 8, -1, -1, -1, -1, -1, -1, -1, 8 ]);
  33.  
  34.     var get_data_type = (function () {
  35.         return function(type) {
  36.             return (type & 0xFF00);
  37.         }
  38.     })();
  39.  
  40.     var get_channel = (function () {
  41.         return function(type) {
  42.             return (type & 0xFF);
  43.         }
  44.     })();
  45.  
  46.     var get_data_type_size = (function () {
  47.         return function(type) {
  48.             return _data_type_size[(type & 0xFF00) >> 8];
  49.         }
  50.     })();
  51.  
  52.     // color conversion
  53.     var COLOR_RGBA2GRAY = 0;
  54.     var COLOR_RGB2GRAY = 1;
  55.     var COLOR_BGRA2GRAY = 2;
  56.     var COLOR_BGR2GRAY = 3;
  57.  
  58.     // box blur option
  59.     var BOX_BLUR_NOSCALE = 0x01;
  60.     // svd options
  61.     var SVD_U_T = 0x01;
  62.     var SVD_V_T = 0x02;
  63.  
  64.     var data_t = (function () {
  65.         function data_t(size_in_bytes, buffer) {
  66.             // we need align size to multiple of 8
  67.             this.size = ((size_in_bytes + 7) | 0) & -8;
  68.             if (typeof buffer === "undefined") {
  69.                 this.buffer = new ArrayBuffer(this.size);
  70.             } else {
  71.                 this.buffer = buffer;
  72.                 this.size = buffer.length;
  73.             }
  74.             this.u8 = new Uint8Array(this.buffer);
  75.             this.i32 = new Int32Array(this.buffer);
  76.             this.f32 = new Float32Array(this.buffer);
  77.             this.f64 = new Float64Array(this.buffer);
  78.         }
  79.         return data_t;
  80.     })();
  81.  
  82.     var matrix_t = (function () {
  83.         // columns, rows, data_type
  84.         function matrix_t(c, r, data_type, data_buffer) {
  85.             this.type = get_data_type(data_type)|0;
  86.             this.channel = get_channel(data_type)|0;
  87.             this.cols = c|0;
  88.             this.rows = r|0;
  89.             if (typeof data_buffer === "undefined") {
  90.                 this.allocate();
  91.             } else {
  92.                 this.buffer = data_buffer;
  93.                 // data user asked for
  94.                 this.data = this.type&U8_t ? this.buffer.u8 : (this.type&S32_t ? this.buffer.i32 : (this.type&F32_t ? this.buffer.f32 : this.buffer.f64));
  95.             }
  96.         }
  97.         matrix_t.prototype.allocate = function() {
  98.             // clear references
  99.             delete this.data;
  100.             delete this.buffer;
  101.             //
  102.             this.buffer = new data_t((this.cols * get_data_type_size(this.type) * this.channel) * this.rows);
  103.             this.data = this.type&U8_t ? this.buffer.u8 : (this.type&S32_t ? this.buffer.i32 : (this.type&F32_t ? this.buffer.f32 : this.buffer.f64));
  104.         }
  105.         matrix_t.prototype.copy_to = function(other) {
  106.             var od = other.data, td = this.data;
  107.             var i = 0, n = (this.cols*this.rows*this.channel)|0;
  108.             for(; i < n-4; i+=4) {
  109.                 od[i] = td[i];
  110.                 od[i+1] = td[i+1];
  111.                 od[i+2] = td[i+2];
  112.                 od[i+3] = td[i+3];
  113.             }
  114.             for(; i < n; ++i) {
  115.                 od[i] = td[i];
  116.             }
  117.         }
  118.         matrix_t.prototype.resize = function(c, r, ch) {
  119.             if (typeof ch === "undefined") { ch = this.channel; }
  120.             // relocate buffer only if new size doesnt fit
  121.             var new_size = (c * get_data_type_size(this.type) * ch) * r;
  122.             if(new_size > this.buffer.size) {
  123.                 this.cols = c;
  124.                 this.rows = r;
  125.                 this.channel = ch;
  126.                 this.allocate();
  127.             } else {
  128.                 this.cols = c;
  129.                 this.rows = r;
  130.                 this.channel = ch;
  131.             }
  132.         }
  133.  
  134.         return matrix_t;
  135.     })();
  136.  
  137.     var pyramid_t = (function () {
  138.  
  139.         function pyramid_t(levels) {
  140.             this.levels = levels|0;
  141.             this.data = new Array(levels);
  142.             this.pyrdown = jsfeat.imgproc.pyrdown;
  143.         }
  144.  
  145.         pyramid_t.prototype.allocate = function(start_w, start_h, data_type) {
  146.             var i = this.levels;
  147.             while(--i >= 0) {
  148.                 this.data[i] = new matrix_t(start_w >> i, start_h >> i, data_type);
  149.             }
  150.         }
  151.  
  152.         pyramid_t.prototype.build = function(input, skip_first_level) {
  153.             if (typeof skip_first_level === "undefined") { skip_first_level = true; }
  154.             // just copy data to first level
  155.             var i = 2, a = input, b = this.data[0];
  156.             if(!skip_first_level) {
  157.                 var j=input.cols*input.rows;
  158.                 while(--j >= 0) {
  159.                     b.data[j] = input.data[j];
  160.                 }
  161.             }
  162.             b = this.data[1];
  163.             this.pyrdown(a, b);
  164.             for(; i < this.levels; ++i) {
  165.                 a = b;
  166.                 b = this.data[i];
  167.                 this.pyrdown(a, b);
  168.             }
  169.         }
  170.  
  171.         return pyramid_t;
  172.     })();
  173.  
  174.     var keypoint_t = (function () {
  175.         function keypoint_t(x,y,score,level,angle) {
  176.             if (typeof x === "undefined") { x=0; }
  177.             if (typeof y === "undefined") { y=0; }
  178.             if (typeof score === "undefined") { score=0; }
  179.             if (typeof level === "undefined") { level=0; }
  180.             if (typeof angle === "undefined") { angle=-1.0; }
  181.  
  182.             this.x = x;
  183.             this.y = y;
  184.             this.score = score;
  185.             this.level = level;
  186.             this.angle = angle;
  187.         }
  188.         return keypoint_t;
  189.     })();
  190.  
  191.  
  192.     // data types
  193.     global.U8_t = U8_t;
  194.     global.S32_t = S32_t;
  195.     global.F32_t = F32_t;
  196.     global.S64_t = S64_t;
  197.     global.F64_t = F64_t;
  198.     // data channels
  199.     global.C1_t = C1_t;
  200.     global.C2_t = C2_t;
  201.     global.C3_t = C3_t;
  202.     global.C4_t = C4_t;
  203.  
  204.     // popular formats
  205.     global.U8C1_t = U8_t | C1_t;
  206.     global.U8C3_t = U8_t | C3_t;
  207.     global.U8C4_t = U8_t | C4_t;
  208.  
  209.     global.F32C1_t = F32_t | C1_t;
  210.     global.F32C2_t = F32_t | C2_t;
  211.     global.S32C1_t = S32_t | C1_t;
  212.     global.S32C2_t = S32_t | C2_t;
  213.  
  214.     // constants
  215.     global.EPSILON = EPSILON;
  216.     global.FLT_MIN = FLT_MIN;
  217.  
  218.     // color convert
  219.     global.COLOR_RGBA2GRAY = COLOR_RGBA2GRAY;
  220.     global.COLOR_RGB2GRAY = COLOR_RGB2GRAY;
  221.     global.COLOR_BGRA2GRAY = COLOR_BGRA2GRAY;
  222.     global.COLOR_BGR2GRAY = COLOR_BGR2GRAY;
  223.  
  224.     // options
  225.     global.BOX_BLUR_NOSCALE = BOX_BLUR_NOSCALE;
  226.     global.SVD_U_T = SVD_U_T;
  227.     global.SVD_V_T = SVD_V_T;
  228.  
  229.     global.get_data_type = get_data_type;
  230.     global.get_channel = get_channel;
  231.     global.get_data_type_size = get_data_type_size;
  232.  
  233.     global.data_t = data_t;
  234.     global.matrix_t = matrix_t;
  235.     global.pyramid_t = pyramid_t;
  236.     global.keypoint_t = keypoint_t;
  237.  
  238. })(jsfeat);
  239. /**
  240.  * @author Eugene Zatepyakin / http://inspirit.ru/
  241.  */
  242.  
  243. (function(global) {
  244.     "use strict";
  245.     //
  246.  
  247.     var cache = (function() {
  248.  
  249.         // very primitive array cache, still need testing if it helps
  250.         // of course V8 has its own powerful cache sys but i'm not sure
  251.         // it caches several multichannel 640x480 buffer creations each frame
  252.  
  253.         var _pool_node_t = (function () {
  254.             function _pool_node_t(size_in_bytes) {
  255.                 this.next = null;
  256.                 this.data = new jsfeat.data_t(size_in_bytes);
  257.                 this.size = this.data.size;
  258.                 this.buffer = this.data.buffer;
  259.                 this.u8 = this.data.u8;
  260.                 this.i32 = this.data.i32;
  261.                 this.f32 = this.data.f32;
  262.                 this.f64 = this.data.f64;
  263.             }
  264.             _pool_node_t.prototype.resize = function(size_in_bytes) {
  265.                 delete this.data;
  266.                 this.data = new jsfeat.data_t(size_in_bytes);
  267.                 this.size = this.data.size;
  268.                 this.buffer = this.data.buffer;
  269.                 this.u8 = this.data.u8;
  270.                 this.i32 = this.data.i32;
  271.                 this.f32 = this.data.f32;
  272.                 this.f64 = this.data.f64;
  273.             }
  274.             return _pool_node_t;
  275.         })();
  276.  
  277.         var _pool_head, _pool_tail;
  278.         var _pool_size = 0;
  279.  
  280.         return {
  281.  
  282.             allocate: function(capacity, data_size) {
  283.                 _pool_head = _pool_tail = new _pool_node_t(data_size);
  284.                 for (var i = 0; i < capacity; ++i) {
  285.                     var node = new _pool_node_t(data_size);
  286.                     _pool_tail = _pool_tail.next = node;
  287.  
  288.                     _pool_size++;
  289.                 }
  290.             },
  291.  
  292.             get_buffer: function(size_in_bytes) {
  293.                 // assume we have enough free nodes
  294.                 var node = _pool_head;
  295.                 _pool_head = _pool_head.next;
  296.                 _pool_size--;
  297.  
  298.                 if(size_in_bytes > node.size) {
  299.                     node.resize(size_in_bytes);
  300.                 }
  301.  
  302.                 return node;
  303.             },
  304.  
  305.             put_buffer: function(node) {
  306.                 _pool_tail = _pool_tail.next = node;
  307.                 _pool_size++;
  308.             }
  309.         };
  310.     })();
  311.  
  312.     global.cache = cache;
  313.     // for now we dont need more than 30 buffers
  314.     // if having cache sys really helps we can add auto extending sys
  315.     cache.allocate(30, 640*4);
  316.  
  317. })(jsfeat);
  318. /**
  319.  * @author Eugene Zatepyakin / http://inspirit.ru/
  320.  */
  321.  
  322. (function(global) {
  323.     "use strict";
  324.     //
  325.  
  326.     var math = (function() {
  327.  
  328.         var qsort_stack = new Int32Array(48*2);
  329.  
  330.         return {
  331.             get_gaussian_kernel: function(size, sigma, kernel, data_type) {
  332.                 var i=0,x=0.0,t=0.0,sigma_x=0.0,scale_2x=0.0;
  333.                 var sum = 0.0;
  334.                 var kern_node = jsfeat.cache.get_buffer(size<<2);
  335.                 var _kernel = kern_node.f32;//new Float32Array(size);
  336.  
  337.                 if((size&1) == 1 && size <= 7 && sigma <= 0) {
  338.                     switch(size>>1) {
  339.                         case 0:
  340.                         _kernel[0] = 1.0;
  341.                         sum = 1.0;
  342.                         break;
  343.                         case 1:
  344.                         _kernel[0] = 0.25, _kernel[1] = 0.5, _kernel[2] = 0.25;
  345.                         sum = 0.25+0.5+0.25;
  346.                         break;
  347.                         case 2:
  348.                         _kernel[0] = 0.0625, _kernel[1] = 0.25, _kernel[2] = 0.375,
  349.                         _kernel[3] = 0.25, _kernel[4] = 0.0625;
  350.                         sum = 0.0625+0.25+0.375+0.25+0.0625;
  351.                         break;
  352.                         case 3:
  353.                         _kernel[0] = 0.03125, _kernel[1] = 0.109375, _kernel[2] = 0.21875,
  354.                         _kernel[3] = 0.28125, _kernel[4] = 0.21875, _kernel[5] = 0.109375, _kernel[6] = 0.03125;
  355.                         sum = 0.03125+0.109375+0.21875+0.28125+0.21875+0.109375+0.03125;
  356.                         break;
  357.                     }
  358.                 } else {
  359.                     sigma_x = sigma > 0 ? sigma : ((size-1)*0.5 - 1.0)*0.3 + 0.8;
  360.                     scale_2x = -0.5/(sigma_x*sigma_x);
  361.  
  362.                     for( ; i < size; ++i )
  363.                     {
  364.                         x = i - (size-1)*0.5;
  365.                         t = Math.exp(scale_2x*x*x);
  366.  
  367.                         _kernel[i] = t;
  368.                         sum += t;
  369.                     }
  370.                 }
  371.  
  372.                 if(data_type & jsfeat.U8_t) {
  373.                     // int based kernel
  374.                     sum = 256.0/sum;
  375.                     for (i = 0; i < size; ++i) {
  376.                         kernel[i] = (_kernel[i] * sum + 0.5)|0;
  377.                     }
  378.                 } else {
  379.                     // classic kernel
  380.                     sum = 1.0/sum;
  381.                     for (i = 0; i < size; ++i) {
  382.                         kernel[i] = _kernel[i] * sum;
  383.                     }
  384.                 }
  385.  
  386.                 jsfeat.cache.put_buffer(kern_node);
  387.             },
  388.  
  389.             // model is 3x3 matrix_t
  390.             perspective_4point_transform: function(model, src_x0, src_y0, dst_x0, dst_y0,
  391.                                                         src_x1, src_y1, dst_x1, dst_y1,
  392.                                                         src_x2, src_y2, dst_x2, dst_y2,
  393.                                                         src_x3, src_y3, dst_x3, dst_y3) {
  394.                 var t1 = src_x0;
  395.                 var t2 = src_x2;
  396.                 var t4 = src_y1;
  397.                 var t5 = t1 * t2 * t4;
  398.                 var t6 = src_y3;
  399.                 var t7 = t1 * t6;
  400.                 var t8 = t2 * t7;
  401.                 var t9 = src_y2;
  402.                 var t10 = t1 * t9;
  403.                 var t11 = src_x1;
  404.                 var t14 = src_y0;
  405.                 var t15 = src_x3;
  406.                 var t16 = t14 * t15;
  407.                 var t18 = t16 * t11;
  408.                 var t20 = t15 * t11 * t9;
  409.                 var t21 = t15 * t4;
  410.                 var t24 = t15 * t9;
  411.                 var t25 = t2 * t4;
  412.                 var t26 = t6 * t2;
  413.                 var t27 = t6 * t11;
  414.                 var t28 = t9 * t11;
  415.                 var t30 = 1.0 / (t21-t24 - t25 + t26 - t27 + t28);
  416.                 var t32 = t1 * t15;
  417.                 var t35 = t14 * t11;
  418.                 var t41 = t4 * t1;
  419.                 var t42 = t6 * t41;
  420.                 var t43 = t14 * t2;
  421.                 var t46 = t16 * t9;
  422.                 var t48 = t14 * t9 * t11;
  423.                 var t51 = t4 * t6 * t2;
  424.                 var t55 = t6 * t14;
  425.                 var Hr0 = -(t8-t5 + t10 * t11 - t11 * t7 - t16 * t2 + t18 - t20 + t21 * t2) * t30;
  426.                 var Hr1 = (t5 - t8 - t32 * t4 + t32 * t9 + t18 - t2 * t35 + t27 * t2 - t20) * t30;
  427.                 var Hr2 = t1;
  428.                 var Hr3 = (-t9 * t7 + t42 + t43 * t4 - t16 * t4 + t46 - t48 + t27 * t9 - t51) * t30;
  429.                 var Hr4 = (-t42 + t41 * t9 - t55 * t2 + t46 - t48 + t55 * t11 + t51 - t21 * t9) * t30;
  430.                 var Hr5 = t14;
  431.                 var Hr6 = (-t10 + t41 + t43 - t35 + t24 - t21 - t26 + t27) * t30;
  432.                 var Hr7 = (-t7 + t10 + t16 - t43 + t27 - t28 - t21 + t25) * t30;
  433.                
  434.                 t1 = dst_x0;
  435.                 t2 = dst_x2;
  436.                 t4 = dst_y1;
  437.                 t5 = t1 * t2 * t4;
  438.                 t6 = dst_y3;
  439.                 t7 = t1 * t6;
  440.                 t8 = t2 * t7;
  441.                 t9 = dst_y2;
  442.                 t10 = t1 * t9;
  443.                 t11 = dst_x1;
  444.                 t14 = dst_y0;
  445.                 t15 = dst_x3;
  446.                 t16 = t14 * t15;
  447.                 t18 = t16 * t11;
  448.                 t20 = t15 * t11 * t9;
  449.                 t21 = t15 * t4;
  450.                 t24 = t15 * t9;
  451.                 t25 = t2 * t4;
  452.                 t26 = t6 * t2;
  453.                 t27 = t6 * t11;
  454.                 t28 = t9 * t11;
  455.                 t30 = 1.0 / (t21-t24 - t25 + t26 - t27 + t28);
  456.                 t32 = t1 * t15;
  457.                 t35 = t14 * t11;
  458.                 t41 = t4 * t1;
  459.                 t42 = t6 * t41;
  460.                 t43 = t14 * t2;
  461.                 t46 = t16 * t9;
  462.                 t48 = t14 * t9 * t11;
  463.                 t51 = t4 * t6 * t2;
  464.                 t55 = t6 * t14;
  465.                 var Hl0 = -(t8-t5 + t10 * t11 - t11 * t7 - t16 * t2 + t18 - t20 + t21 * t2) * t30;
  466.                 var Hl1 = (t5 - t8 - t32 * t4 + t32 * t9 + t18 - t2 * t35 + t27 * t2 - t20) * t30;
  467.                 var Hl2 = t1;
  468.                 var Hl3 = (-t9 * t7 + t42 + t43 * t4 - t16 * t4 + t46 - t48 + t27 * t9 - t51) * t30;
  469.                 var Hl4 = (-t42 + t41 * t9 - t55 * t2 + t46 - t48 + t55 * t11 + t51 - t21 * t9) * t30;
  470.                 var Hl5 = t14;
  471.                 var Hl6 = (-t10 + t41 + t43 - t35 + t24 - t21 - t26 + t27) * t30;
  472.                 var Hl7 = (-t7 + t10 + t16 - t43 + t27 - t28 - t21 + t25) * t30;
  473.  
  474.                 // the following code computes R = Hl * inverse Hr
  475.                 t2 = Hr4-Hr7*Hr5;
  476.                 t4 = Hr0*Hr4;
  477.                 t5 = Hr0*Hr5;
  478.                 t7 = Hr3*Hr1;
  479.                 t8 = Hr2*Hr3;
  480.                 t10 = Hr1*Hr6;
  481.                 var t12 = Hr2*Hr6;
  482.                 t15 = 1.0 / (t4-t5*Hr7-t7+t8*Hr7+t10*Hr5-t12*Hr4);
  483.                 t18 = -Hr3+Hr5*Hr6;
  484.                 var t23 = -Hr3*Hr7+Hr4*Hr6;
  485.                 t28 = -Hr1+Hr2*Hr7;
  486.                 var t31 = Hr0-t12;
  487.                 t35 = Hr0*Hr7-t10;
  488.                 t41 = -Hr1*Hr5+Hr2*Hr4;
  489.                 var t44 = t5-t8;
  490.                 var t47 = t4-t7;
  491.                 t48 = t2*t15;
  492.                 var t49 = t28*t15;
  493.                 var t50 = t41*t15;
  494.                 var mat = model.data;
  495.                 mat[0] = Hl0*t48+Hl1*(t18*t15)-Hl2*(t23*t15);
  496.                 mat[1] = Hl0*t49+Hl1*(t31*t15)-Hl2*(t35*t15);
  497.                 mat[2] = -Hl0*t50-Hl1*(t44*t15)+Hl2*(t47*t15);
  498.                 mat[3] = Hl3*t48+Hl4*(t18*t15)-Hl5*(t23*t15);
  499.                 mat[4] = Hl3*t49+Hl4*(t31*t15)-Hl5*(t35*t15);
  500.                 mat[5] = -Hl3*t50-Hl4*(t44*t15)+Hl5*(t47*t15);
  501.                 mat[6] = Hl6*t48+Hl7*(t18*t15)-t23*t15;
  502.                 mat[7] = Hl6*t49+Hl7*(t31*t15)-t35*t15;
  503.                 mat[8] = -Hl6*t50-Hl7*(t44*t15)+t47*t15;
  504.             },
  505.  
  506.             // The current implementation was derived from *BSD system qsort():
  507.             // Copyright (c) 1992, 1993
  508.             // The Regents of the University of California.  All rights reserved.
  509.             qsort: function(array, low, high, cmp) {
  510.                 var isort_thresh = 7;
  511.                 var t,ta,tb,tc;
  512.                 var sp = 0,left=0,right=0,i=0,n=0,m=0,ptr=0,ptr2=0,d=0;
  513.                 var left0=0,left1=0,right0=0,right1=0,pivot=0,a=0,b=0,c=0,swap_cnt=0;
  514.  
  515.                 var stack = qsort_stack;
  516.  
  517.                 if( (high-low+1) <= 1 ) return;
  518.  
  519.                 stack[0] = low;
  520.                 stack[1] = high;
  521.  
  522.                 while( sp >= 0 ) {
  523.                
  524.                     left = stack[sp<<1];
  525.                     right = stack[(sp<<1)+1];
  526.                     sp--;
  527.  
  528.                     for(;;) {
  529.                         n = (right - left) + 1;
  530.  
  531.                         if( n <= isort_thresh ) {
  532.                         //insert_sort:
  533.                             for( ptr = left + 1; ptr <= right; ptr++ ) {
  534.                                 for( ptr2 = ptr; ptr2 > left && cmp(array[ptr2],array[ptr2-1]); ptr2--) {
  535.                                     t = array[ptr2];
  536.                                     array[ptr2] = array[ptr2-1];
  537.                                     array[ptr2-1] = t;
  538.                                 }
  539.                             }
  540.                             break;
  541.                         } else {
  542.                             swap_cnt = 0;
  543.  
  544.                             left0 = left;
  545.                             right0 = right;
  546.                             pivot = left + (n>>1);
  547.  
  548.                             if( n > 40 ) {
  549.                                 d = n >> 3;
  550.                                 a = left, b = left + d, c = left + (d<<1);
  551.                                 ta = array[a],tb = array[b],tc = array[c];
  552.                                 left = cmp(ta, tb) ? (cmp(tb, tc) ? b : (cmp(ta, tc) ? c : a))
  553.                                                   : (cmp(tc, tb) ? b : (cmp(ta, tc) ? a : c));
  554.  
  555.                                 a = pivot - d, b = pivot, c = pivot + d;
  556.                                 ta = array[a],tb = array[b],tc = array[c];
  557.                                 pivot = cmp(ta, tb) ? (cmp(tb, tc) ? b : (cmp(ta, tc) ? c : a))
  558.                                                   : (cmp(tc, tb) ? b : (cmp(ta, tc) ? a : c));
  559.  
  560.                                 a = right - (d<<1), b = right - d, c = right;
  561.                                 ta = array[a],tb = array[b],tc = array[c];
  562.                                 right = cmp(ta, tb) ? (cmp(tb, tc) ? b : (cmp(ta, tc) ? c : a))
  563.                                                   : (cmp(tc, tb) ? b : (cmp(ta, tc) ? a : c));
  564.                             }
  565.  
  566.                             a = left, b = pivot, c = right;
  567.                             ta = array[a],tb = array[b],tc = array[c];
  568.                             pivot = cmp(ta, tb) ? (cmp(tb, tc) ? b : (cmp(ta, tc) ? c : a))  
  569.                                                : (cmp(tc, tb) ? b : (cmp(ta, tc) ? a : c));
  570.                             if( pivot != left0 ) {
  571.                                 t = array[pivot];
  572.                                 array[pivot] = array[left0];
  573.                                 array[left0] = t;
  574.                                 pivot = left0;
  575.                             }
  576.                             left = left1 = left0 + 1;
  577.                             right = right1 = right0;
  578.  
  579.                             ta = array[pivot];
  580.                             for(;;) {
  581.                                 while( left <= right && !cmp(ta, array[left]) ) {
  582.                                     if( !cmp(array[left], ta) ) {
  583.                                         if( left > left1 ) {
  584.                                             t = array[left1];
  585.                                             array[left1] = array[left];
  586.                                             array[left] = t;
  587.                                         }
  588.                                         swap_cnt = 1;
  589.                                         left1++;
  590.                                     }
  591.                                     left++;
  592.                                 }
  593.  
  594.                                 while( left <= right && !cmp(array[right], ta) ) {
  595.                                     if( !cmp(ta, array[right]) ) {
  596.                                         if( right < right1 ) {
  597.                                             t = array[right1];
  598.                                             array[right1] = array[right];
  599.                                             array[right] = t;
  600.                                         }
  601.                                         swap_cnt = 1;
  602.                                         right1--;
  603.                                     }
  604.                                     right--;
  605.                                 }
  606.  
  607.                                 if( left > right ) break;
  608.                                
  609.                                 t = array[left];
  610.                                 array[left] = array[right];
  611.                                 array[right] = t;
  612.                                 swap_cnt = 1;
  613.                                 left++;
  614.                                 right--;
  615.                             }
  616.  
  617.                             if( swap_cnt == 0 ) {
  618.                                 left = left0, right = right0;
  619.                                 //goto insert_sort;
  620.                                 for( ptr = left + 1; ptr <= right; ptr++ ) {
  621.                                     for( ptr2 = ptr; ptr2 > left && cmp(array[ptr2],array[ptr2-1]); ptr2--) {
  622.                                         t = array[ptr2];
  623.                                         array[ptr2] = array[ptr2-1];
  624.                                         array[ptr2-1] = t;
  625.                                     }
  626.                                 }
  627.                                 break;
  628.                             }
  629.  
  630.                             n = Math.min( (left1 - left0), (left - left1) );
  631.                             m = (left-n)|0;
  632.                             for( i = 0; i < n; ++i,++m ) {
  633.                                 t = array[left0+i];
  634.                                 array[left0+i] = array[m];
  635.                                 array[m] = t;
  636.                             }
  637.  
  638.                             n = Math.min( (right0 - right1), (right1 - right) );
  639.                             m = (right0-n+1)|0;
  640.                             for( i = 0; i < n; ++i,++m ) {
  641.                                 t = array[left+i];
  642.                                 array[left+i] = array[m];
  643.                                 array[m] = t;
  644.                             }
  645.                             n = (left - left1);
  646.                             m = (right1 - right);
  647.                             if( n > 1 ) {
  648.                                 if( m > 1 ) {
  649.                                     if( n > m ) {
  650.                                         ++sp;
  651.                                         stack[sp<<1] = left0;
  652.                                         stack[(sp<<1)+1] = left0 + n - 1;
  653.                                         left = right0 - m + 1, right = right0;
  654.                                     } else {
  655.                                         ++sp;
  656.                                         stack[sp<<1] = right0 - m + 1;
  657.                                         stack[(sp<<1)+1] = right0;
  658.                                         left = left0, right = left0 + n - 1;
  659.                                     }
  660.                                 } else {
  661.                                     left = left0, right = left0 + n - 1;
  662.                                 }
  663.                             }
  664.                             else if( m > 1 )
  665.                                 left = right0 - m + 1, right = right0;
  666.                             else
  667.                                 break;
  668.                         }
  669.                     }
  670.                 }
  671.             },
  672.  
  673.             median: function(array, low, high) {
  674.                 var w;
  675.                 var middle=0,ll=0,hh=0,median=(low+high)>>1;
  676.                 for (;;) {
  677.                     if (high <= low) return array[median];
  678.                     if (high == (low + 1)) {
  679.                         if (array[low] > array[high]) {
  680.                             w = array[low];
  681.                             array[low] = array[high];
  682.                             array[high] = w;
  683.                         }
  684.                         return array[median];
  685.                     }
  686.                     middle = ((low + high) >> 1);
  687.                     if (array[middle] > array[high]) {
  688.                         w = array[middle];
  689.                         array[middle] = array[high];
  690.                         array[high] = w;
  691.                     }
  692.                     if (array[low] > array[high]) {
  693.                         w = array[low];
  694.                         array[low] = array[high];
  695.                         array[high] = w;
  696.                     }
  697.                     if (array[middle] > array[low]) {
  698.                         w = array[middle];
  699.                         array[middle] = array[low];
  700.                         array[low] = w;
  701.                     }
  702.                     ll = (low + 1);
  703.                     w = array[middle];
  704.                     array[middle] = array[ll];
  705.                     array[ll] = w;
  706.                     hh = high;
  707.                     for (;;) {
  708.                         do ++ll; while (array[low] > array[ll]);
  709.                         do --hh; while (array[hh] > array[low]);
  710.                         if (hh < ll) break;
  711.                         w = array[ll];
  712.                         array[ll] = array[hh];
  713.                         array[hh] = w;
  714.                     }
  715.                     w = array[low];
  716.                     array[low] = array[hh];
  717.                     array[hh] = w;
  718.                     if (hh <= median)
  719.                         low = ll;
  720.                     else if (hh >= median)
  721.                         high = (hh - 1);
  722.                 }
  723.                 return 0;
  724.             }
  725.         };
  726.  
  727.     })();
  728.  
  729.     global.math = math;
  730.  
  731. })(jsfeat);
  732. /**
  733.  * @author Eugene Zatepyakin / http://inspirit.ru/
  734.  *
  735.  */
  736.  
  737. (function(global) {
  738.     "use strict";
  739.     //
  740.  
  741.     var matmath = (function() {
  742.        
  743.         return {
  744.             identity: function(M, value) {
  745.                 if (typeof value === "undefined") { value=1; }
  746.                 var src=M.data;
  747.                 var rows=M.rows, cols=M.cols, cols_1=(cols+1)|0;
  748.                 var len = rows * cols;
  749.                 var k = len;
  750.                 while(--len >= 0) src[len] = 0.0;
  751.                 len = k;
  752.                 k = 0;
  753.                 while(k < len)  {
  754.                     src[k] = value;
  755.                     k = k + cols_1;
  756.                 }
  757.             },
  758.  
  759.             transpose: function(At, A) {
  760.                 var i=0,j=0,nrows=A.rows,ncols=A.cols;
  761.                 var Ai=0,Ati=0,pAt=0;
  762.                 var ad=A.data,atd=At.data;
  763.  
  764.                 for (; i < nrows; Ati += 1, Ai += ncols, i++) {
  765.                     pAt = Ati;
  766.                     for (j = 0; j < ncols; pAt += nrows, j++) atd[pAt] = ad[Ai+j];
  767.                 }
  768.             },
  769.  
  770.             // C = A * B
  771.             multiply: function(C, A, B) {
  772.                 var i=0,j=0,k=0;
  773.                 var Ap=0,pA=0,pB=0,p_B=0,Cp=0;
  774.                 var ncols=A.cols,nrows=A.rows,mcols=B.cols;
  775.                 var ad=A.data,bd=B.data,cd=C.data;
  776.                 var sum=0.0;
  777.  
  778.                 for (; i < nrows; Ap += ncols, i++) {
  779.                     for (p_B = 0, j = 0; j < mcols; Cp++, p_B++, j++) {
  780.                         pB = p_B;
  781.                         pA = Ap;
  782.                         sum = 0.0;
  783.                         for (k = 0; k < ncols; pA++, pB += mcols, k++) {
  784.                             sum += ad[pA] * bd[pB];
  785.                         }
  786.                         cd[Cp] = sum;
  787.                     }
  788.                 }
  789.             },
  790.  
  791.             // C = A * B'
  792.             multiply_ABt: function(C, A, B) {
  793.                 var i=0,j=0,k=0;
  794.                 var Ap=0,pA=0,pB=0,Cp=0;
  795.                 var ncols=A.cols,nrows=A.rows,mrows=B.rows;
  796.                 var ad=A.data,bd=B.data,cd=C.data;
  797.                 var sum=0.0;
  798.  
  799.                 for (; i < nrows; Ap += ncols, i++) {
  800.                     for (pB = 0, j = 0; j < mrows; Cp++, j++) {
  801.                         pA = Ap;
  802.                         sum = 0.0;
  803.                         for (k = 0; k < ncols; pA++, pB++, k++) {
  804.                             sum += ad[pA] * bd[pB];
  805.                         }
  806.                         cd[Cp] = sum;
  807.                     }
  808.                 }
  809.             },
  810.  
  811.             // C = A' * B
  812.             multiply_AtB: function(C, A, B) {
  813.                 var i=0,j=0,k=0;
  814.                 var Ap=0,pA=0,pB=0,p_B=0,Cp=0;
  815.                 var ncols=A.cols,nrows=A.rows,mcols=B.cols;
  816.                 var ad=A.data,bd=B.data,cd=C.data;
  817.                 var sum=0.0;
  818.  
  819.                 for (; i < ncols; Ap++, i++) {
  820.                     for (p_B = 0, j = 0; j < mcols; Cp++, p_B++, j++) {
  821.                         pB = p_B;
  822.                         pA = Ap;
  823.                         sum = 0.0;
  824.                         for (k = 0; k < nrows; pA += ncols, pB += mcols, k++) {
  825.                             sum += ad[pA] * bd[pB];
  826.                         }
  827.                         cd[Cp] = sum;
  828.                     }
  829.                 }
  830.             },
  831.  
  832.             // C = A * A'
  833.             multiply_AAt: function(C, A) {
  834.                 var i=0,j=0,k=0;
  835.                 var pCdiag=0,p_A=0,pA=0,pB=0,pC=0,pCt=0;
  836.                 var ncols=A.cols,nrows=A.rows;
  837.                 var ad=A.data,cd=C.data;
  838.                 var sum=0.0;
  839.  
  840.                 for (; i < nrows; pCdiag += nrows + 1, p_A = pA, i++) {
  841.                     pC = pCdiag;
  842.                     pCt = pCdiag;
  843.                     pB = p_A;
  844.                     for (j = i; j < nrows; pC++, pCt += nrows, j++) {
  845.                         pA = p_A;
  846.                         sum = 0.0;
  847.                         for (k = 0; k < ncols; k++) {
  848.                             sum += ad[pA++] * ad[pB++];
  849.                         }
  850.                         cd[pC] = sum
  851.                         cd[pCt] = sum;
  852.                     }
  853.                 }
  854.             },
  855.  
  856.             // C = A' * A
  857.             multiply_AtA: function(C, A) {
  858.                 var i=0,j=0,k=0;
  859.                 var p_A=0,pA=0,pB=0,p_C=0,pC=0,p_CC=0;
  860.                 var ncols=A.cols,nrows=A.rows;
  861.                 var ad=A.data,cd=C.data;
  862.                 var sum=0.0;
  863.  
  864.                 for (; i < ncols; p_C += ncols, i++) {
  865.                     p_A = i;
  866.                     p_CC = p_C + i;
  867.                     pC = p_CC;
  868.                     for (j = i; j < ncols; pC++, p_CC += ncols, j++) {
  869.                         pA = p_A;
  870.                         pB = j;
  871.                         sum = 0.0;
  872.                         for (k = 0; k < nrows; pA += ncols, pB += ncols, k++) {
  873.                             sum += ad[pA] * ad[pB];
  874.                         }
  875.                         cd[pC] = sum
  876.                         cd[p_CC] = sum;
  877.                     }
  878.                 }
  879.             },
  880.  
  881.             // various small matrix operations
  882.             identity_3x3: function(M, value) {
  883.                 if (typeof value === "undefined") { value=1; }
  884.                 var dt=M.data;
  885.                 dt[0] = dt[4] = dt[8] = value;
  886.                 dt[1] = dt[2] = dt[3] = 0;
  887.                 dt[5] = dt[6] = dt[7] = 0;
  888.             },
  889.  
  890.             invert_3x3: function(from, to) {
  891.                 var A = from.data, invA = to.data;
  892.                 var t1 = A[4];
  893.                 var t2 = A[8];
  894.                 var t4 = A[5];
  895.                 var t5 = A[7];
  896.                 var t8 = A[0];
  897.  
  898.                 var t9 = t8*t1;
  899.                 var t11 = t8*t4;
  900.                 var t13 = A[3];
  901.                 var t14 = A[1];
  902.                 var t15 = t13*t14;
  903.                 var t17 = A[2];
  904.                 var t18 = t13*t17;
  905.                 var t20 = A[6];
  906.                 var t21 = t20*t14;
  907.                 var t23 = t20*t17;
  908.                 var t26 = 1.0/(t9*t2-t11*t5-t15*t2+t18*t5+t21*t4-t23*t1);
  909.                 invA[0] = (t1*t2-t4*t5)*t26;
  910.                 invA[1] = -(t14*t2-t17*t5)*t26;
  911.                 invA[2] = -(-t14*t4+t17*t1)*t26;
  912.                 invA[3] = -(t13*t2-t4*t20)*t26;
  913.                 invA[4] = (t8*t2-t23)*t26;
  914.                 invA[5] = -(t11-t18)*t26;
  915.                 invA[6] = -(-t13*t5+t1*t20)*t26;
  916.                 invA[7] = -(t8*t5-t21)*t26;
  917.                 invA[8] = (t9-t15)*t26;
  918.             },
  919.             // C = A * B
  920.             multiply_3x3: function(C, A, B) {
  921.                 var Cd=C.data, Ad=A.data, Bd=B.data;
  922.                 var m1_0 = Ad[0], m1_1 = Ad[1], m1_2 = Ad[2];
  923.                 var m1_3 = Ad[3], m1_4 = Ad[4], m1_5 = Ad[5];
  924.                 var m1_6 = Ad[6], m1_7 = Ad[7], m1_8 = Ad[8];
  925.  
  926.                 var m2_0 = Bd[0], m2_1 = Bd[1], m2_2 = Bd[2];
  927.                 var m2_3 = Bd[3], m2_4 = Bd[4], m2_5 = Bd[5];
  928.                 var m2_6 = Bd[6], m2_7 = Bd[7], m2_8 = Bd[8];
  929.  
  930.                 Cd[0] = m1_0 * m2_0 + m1_1 * m2_3 + m1_2 * m2_6;
  931.                 Cd[1] = m1_0 * m2_1 + m1_1 * m2_4 + m1_2 * m2_7;
  932.                 Cd[2] = m1_0 * m2_2 + m1_1 * m2_5 + m1_2 * m2_8;
  933.                 Cd[3] = m1_3 * m2_0 + m1_4 * m2_3 + m1_5 * m2_6;
  934.                 Cd[4] = m1_3 * m2_1 + m1_4 * m2_4 + m1_5 * m2_7;
  935.                 Cd[5] = m1_3 * m2_2 + m1_4 * m2_5 + m1_5 * m2_8;
  936.                 Cd[6] = m1_6 * m2_0 + m1_7 * m2_3 + m1_8 * m2_6;
  937.                 Cd[7] = m1_6 * m2_1 + m1_7 * m2_4 + m1_8 * m2_7;
  938.                 Cd[8] = m1_6 * m2_2 + m1_7 * m2_5 + m1_8 * m2_8;
  939.             },
  940.  
  941.             mat3x3_determinant: function(M) {
  942.                 var md=M.data;
  943.                 return  md[0] * md[4] * md[8] -
  944.                         md[0] * md[5] * md[7] -
  945.                         md[3] * md[1] * md[8] +
  946.                         md[3] * md[2] * md[7] +
  947.                         md[6] * md[1] * md[5] -
  948.                         md[6] * md[2] * md[4];
  949.             },
  950.  
  951.             determinant_3x3: function(M11, M12, M13,
  952.                                       M21, M22, M23,
  953.                                       M31, M32, M33) {
  954.                 return  M11 * M22 * M33 - M11 * M23 * M32 -
  955.                           M21 * M12 * M33 + M21 * M13 * M32 +
  956.                           M31 * M12 * M23 - M31 * M13 * M22;
  957.             }
  958.         };
  959.  
  960.     })();
  961.  
  962.     global.matmath = matmath;
  963.  
  964. })(jsfeat);
  965. /**
  966.  * @author Eugene Zatepyakin / http://inspirit.ru/
  967.  *
  968.  */
  969.  
  970. (function(global) {
  971.     "use strict";
  972.     //
  973.  
  974.     var linalg = (function() {
  975.  
  976.         var swap = function(A, i0, i1, t) {
  977.             t = A[i0];
  978.             A[i0] = A[i1];
  979.             A[i1] = t;
  980.         }
  981.  
  982.         var hypot = function(a, b) {
  983.             a = Math.abs(a);
  984.             b = Math.abs(b);
  985.             if( a > b ) {
  986.                 b /= a;
  987.                 return a*Math.sqrt(1.0 + b*b);
  988.             }
  989.             if( b > 0 ) {
  990.                 a /= b;
  991.                 return b*Math.sqrt(1.0 + a*a);
  992.             }
  993.             return 0.0;
  994.         }
  995.  
  996.         var JacobiImpl = function(A, astep, W, V, vstep, n) {
  997.             var eps = jsfeat.EPSILON;
  998.             var i=0,j=0,k=0,m=0,l=0,idx=0,_in=0,_in2=0;
  999.             var iters=0,max_iter=n*n*30;
  1000.             var mv=0.0,val=0.0,p=0.0,y=0.0,t=0.0,s=0.0,c=0.0,a0=0.0,b0=0.0;
  1001.  
  1002.             var indR_buff = jsfeat.cache.get_buffer(n<<2);
  1003.             var indC_buff = jsfeat.cache.get_buffer(n<<2);
  1004.             var indR = indR_buff.i32;
  1005.             var indC = indC_buff.i32;
  1006.  
  1007.             if(V) {
  1008.                 for(; i < n; i++) {
  1009.                     k = i*vstep;
  1010.                     for(j = 0; j < n; j++) {
  1011.                         V[k + j] = 0.0;
  1012.                     }
  1013.                     V[k + i] = 1.0;
  1014.                 }
  1015.             }
  1016.  
  1017.             for(k = 0; k < n; k++) {
  1018.                 W[k] = A[(astep + 1)*k];
  1019.                 if(k < n - 1) {
  1020.                     for(m = k+1, mv = Math.abs(A[astep*k + m]), i = k+2; i < n; i++) {
  1021.                         val = Math.abs(A[astep*k+i]);
  1022.                         if(mv < val)
  1023.                             mv = val, m = i;
  1024.                     }
  1025.                     indR[k] = m;
  1026.                 }
  1027.                 if(k > 0) {
  1028.                     for(m = 0, mv = Math.abs(A[k]), i = 1; i < k; i++) {
  1029.                         val = Math.abs(A[astep*i+k]);
  1030.                         if(mv < val)
  1031.                             mv = val, m = i;
  1032.                     }
  1033.                     indC[k] = m;
  1034.                 }
  1035.             }
  1036.  
  1037.             if(n > 1) for( ; iters < max_iter; iters++) {
  1038.                 // find index (k,l) of pivot p
  1039.                 for(k = 0, mv = Math.abs(A[indR[0]]), i = 1; i < n-1; i++) {
  1040.                     val = Math.abs(A[astep*i + indR[i]]);
  1041.                     if( mv < val )
  1042.                         mv = val, k = i;
  1043.                 }
  1044.                 l = indR[k];
  1045.                 for(i = 1; i < n; i++) {
  1046.                     val = Math.abs(A[astep*indC[i] + i]);
  1047.                     if( mv < val )
  1048.                         mv = val, k = indC[i], l = i;
  1049.                 }
  1050.                
  1051.                 p = A[astep*k + l];
  1052.  
  1053.                 if(Math.abs(p) <= eps) break;
  1054.  
  1055.                 y = (W[l] - W[k])*0.5;
  1056.                 t = Math.abs(y) + hypot(p, y);
  1057.                 s = hypot(p, t);
  1058.                 c = t/s;
  1059.                 s = p/s; t = (p/t)*p;
  1060.                 if(y < 0)
  1061.                     s = -s, t = -t;
  1062.                 A[astep*k + l] = 0;
  1063.                
  1064.                 W[k] -= t;
  1065.                 W[l] += t;
  1066.                
  1067.                 // rotate rows and columns k and l
  1068.                 for (i = 0; i < k; i++) {
  1069.                     _in = (astep * i + k);
  1070.                     _in2 = (astep * i + l);
  1071.                     a0 = A[_in];
  1072.                     b0 = A[_in2];
  1073.                     A[_in] = a0 * c - b0 * s;
  1074.                     A[_in2] = a0 * s + b0 * c;
  1075.                 }
  1076.                 for (i = (k + 1); i < l; i++) {
  1077.                     _in = (astep * k + i);
  1078.                     _in2 = (astep * i + l);
  1079.                     a0 = A[_in];
  1080.                     b0 = A[_in2];
  1081.                     A[_in] = a0 * c - b0 * s;
  1082.                     A[_in2] = a0 * s + b0 * c;
  1083.                 }
  1084.                 i = l + 1;
  1085.                 _in = (astep * k + i);
  1086.                 _in2 = (astep * l + i);
  1087.                 for (; i < n; i++, _in++, _in2++) {
  1088.                     a0 = A[_in];
  1089.                     b0 = A[_in2];
  1090.                     A[_in] = a0 * c - b0 * s;
  1091.                     A[_in2] = a0 * s + b0 * c;
  1092.                 }
  1093.                
  1094.                 // rotate eigenvectors
  1095.                 if (V) {
  1096.                     _in = vstep * k;
  1097.                     _in2 = vstep * l;
  1098.                     for (i = 0; i < n; i++, _in++, _in2++) {
  1099.                         a0 = V[_in];
  1100.                         b0 = V[_in2];
  1101.                         V[_in] = a0 * c - b0 * s;
  1102.                         V[_in2] = a0 * s + b0 * c;
  1103.                     }
  1104.                 }
  1105.                
  1106.                 for(j = 0; j < 2; j++) {
  1107.                     idx = j == 0 ? k : l;
  1108.                     if(idx < n - 1) {
  1109.                         for(m = idx+1, mv = Math.abs(A[astep*idx + m]), i = idx+2; i < n; i++) {
  1110.                             val = Math.abs(A[astep*idx+i]);
  1111.                             if( mv < val )
  1112.                                 mv = val, m = i;
  1113.                         }
  1114.                         indR[idx] = m;
  1115.                     }
  1116.                     if(idx > 0) {
  1117.                         for(m = 0, mv = Math.abs(A[idx]), i = 1; i < idx; i++) {
  1118.                             val = Math.abs(A[astep*i+idx]);
  1119.                             if( mv < val )
  1120.                                 mv = val, m = i;
  1121.                         }
  1122.                         indC[idx] = m;
  1123.                     }
  1124.                 }
  1125.             }
  1126.  
  1127.             // sort eigenvalues & eigenvectors
  1128.             for(k = 0; k < n-1; k++) {
  1129.                 m = k;
  1130.                 for(i = k+1; i < n; i++) {
  1131.                     if(W[m] < W[i])
  1132.                         m = i;
  1133.                 }
  1134.                 if(k != m) {
  1135.                     swap(W, m, k, mv);
  1136.                     if(V) {
  1137.                         for(i = 0; i < n; i++) {
  1138.                             swap(V, vstep*m + i, vstep*k + i, mv);
  1139.                         }
  1140.                     }
  1141.                 }
  1142.             }
  1143.  
  1144.  
  1145.             jsfeat.cache.put_buffer(indR_buff);
  1146.             jsfeat.cache.put_buffer(indC_buff);
  1147.         }
  1148.  
  1149.         var JacobiSVDImpl = function(At, astep, _W, Vt, vstep, m, n, n1) {
  1150.             var eps = jsfeat.EPSILON * 2.0;
  1151.             var minval = jsfeat.FLT_MIN;
  1152.             var i=0,j=0,k=0,iter=0,max_iter=Math.max(m, 30);
  1153.             var Ai=0,Aj=0,Vi=0,Vj=0,changed=0;
  1154.             var c=0.0, s=0.0, t=0.0;
  1155.             var t0=0.0,t1=0.0,sd=0.0,beta=0.0,gamma=0.0,delta=0.0,a=0.0,p=0.0,b=0.0;
  1156.             var seed = 0x1234;
  1157.             var val=0.0,val0=0.0,asum=0.0;
  1158.  
  1159.             var W_buff = jsfeat.cache.get_buffer(n<<3);
  1160.             var W = W_buff.f64;
  1161.            
  1162.             for(; i < n; i++) {
  1163.                 for(k = 0, sd = 0; k < m; k++) {
  1164.                     t = At[i*astep + k];
  1165.                     sd += t*t;
  1166.                 }
  1167.                 W[i] = sd;
  1168.                
  1169.                 if(Vt) {
  1170.                     for(k = 0; k < n; k++) {
  1171.                         Vt[i*vstep + k] = 0;
  1172.                     }
  1173.                     Vt[i*vstep + i] = 1;
  1174.                 }
  1175.             }
  1176.            
  1177.             for(; iter < max_iter; iter++) {
  1178.                 changed = 0;
  1179.                
  1180.                 for(i = 0; i < n-1; i++) {
  1181.                     for(j = i+1; j < n; j++) {
  1182.                         Ai = (i*astep)|0, Aj = (j*astep)|0;
  1183.                         a = W[i], p = 0, b = W[j];
  1184.                        
  1185.                         k = 2;
  1186.                         p += At[Ai]*At[Aj];
  1187.                         p += At[Ai+1]*At[Aj+1];
  1188.  
  1189.                         for(; k < m; k++)
  1190.                             p += At[Ai+k]*At[Aj+k];
  1191.                        
  1192.                         if(Math.abs(p) <= eps*Math.sqrt(a*b)) continue;
  1193.                        
  1194.                         p *= 2.0;
  1195.                         beta = a - b, gamma = hypot(p, beta);
  1196.                         if( beta < 0 ) {
  1197.                             delta = (gamma - beta)*0.5;
  1198.                             s = Math.sqrt(delta/gamma);
  1199.                             c = (p/(gamma*s*2.0));
  1200.                         } else {
  1201.                             c = Math.sqrt((gamma + beta)/(gamma*2.0));
  1202.                             s = (p/(gamma*c*2.0));
  1203.                         }
  1204.                        
  1205.                         a=0.0, b=0.0;
  1206.                        
  1207.                         k = 2; // unroll
  1208.                         t0 = c*At[Ai] + s*At[Aj];
  1209.                         t1 = -s*At[Ai] + c*At[Aj];
  1210.                         At[Ai] = t0; At[Aj] = t1;
  1211.                         a += t0*t0; b += t1*t1;
  1212.  
  1213.                         t0 = c*At[Ai+1] + s*At[Aj+1];
  1214.                         t1 = -s*At[Ai+1] + c*At[Aj+1];
  1215.                         At[Ai+1] = t0; At[Aj+1] = t1;
  1216.                         a += t0*t0; b += t1*t1;
  1217.  
  1218.                         for( ; k < m; k++ )
  1219.                         {
  1220.                             t0 = c*At[Ai+k] + s*At[Aj+k];
  1221.                             t1 = -s*At[Ai+k] + c*At[Aj+k];
  1222.                             At[Ai+k] = t0; At[Aj+k] = t1;
  1223.                            
  1224.                             a += t0*t0; b += t1*t1;
  1225.                         }
  1226.                        
  1227.                         W[i] = a; W[j] = b;
  1228.                        
  1229.                         changed = 1;
  1230.                        
  1231.                         if(Vt) {
  1232.                             Vi = (i*vstep)|0, Vj = (j*vstep)|0;
  1233.  
  1234.                             k = 2;
  1235.                             t0 = c*Vt[Vi] + s*Vt[Vj];
  1236.                             t1 = -s*Vt[Vi] + c*Vt[Vj];
  1237.                             Vt[Vi] = t0; Vt[Vj] = t1;
  1238.  
  1239.                             t0 = c*Vt[Vi+1] + s*Vt[Vj+1];
  1240.                             t1 = -s*Vt[Vi+1] + c*Vt[Vj+1];
  1241.                             Vt[Vi+1] = t0; Vt[Vj+1] = t1;
  1242.  
  1243.                             for(; k < n; k++) {
  1244.                                 t0 = c*Vt[Vi+k] + s*Vt[Vj+k];
  1245.                                 t1 = -s*Vt[Vi+k] + c*Vt[Vj+k];
  1246.                                 Vt[Vi+k] = t0; Vt[Vj+k] = t1;
  1247.                             }
  1248.                         }
  1249.                     }
  1250.                 }
  1251.                 if(changed == 0) break;
  1252.             }
  1253.            
  1254.             for(i = 0; i < n; i++) {
  1255.                 for(k = 0, sd = 0; k < m; k++) {
  1256.                     t = At[i*astep + k];
  1257.                     sd += t*t;
  1258.                 }
  1259.                 W[i] = Math.sqrt(sd);
  1260.             }
  1261.            
  1262.             for(i = 0; i < n-1; i++) {
  1263.                 j = i;
  1264.                 for(k = i+1; k < n; k++) {
  1265.                     if(W[j] < W[k])
  1266.                         j = k;
  1267.                 }
  1268.                 if(i != j) {
  1269.                     swap(W, i, j, sd);
  1270.                     if(Vt) {
  1271.                         for(k = 0; k < m; k++) {
  1272.                             swap(At, i*astep + k, j*astep + k, t);
  1273.                         }
  1274.                        
  1275.                         for(k = 0; k < n; k++) {
  1276.                             swap(Vt, i*vstep + k, j*vstep + k, t);
  1277.                         }
  1278.                     }
  1279.                 }
  1280.             }
  1281.            
  1282.             for(i = 0; i < n; i++) {
  1283.                 _W[i] = W[i];
  1284.             }
  1285.            
  1286.             if(!Vt) {
  1287.                 jsfeat.cache.put_buffer(W_buff);
  1288.                 return;
  1289.             }
  1290.  
  1291.             for(i = 0; i < n1; i++) {
  1292.  
  1293.                 sd = i < n ? W[i] : 0;
  1294.                
  1295.                 while(sd <= minval) {
  1296.                     // if we got a zero singular value, then in order to get the corresponding left singular vector
  1297.                     // we generate a random vector, project it to the previously computed left singular vectors,
  1298.                     // subtract the projection and normalize the difference.
  1299.                     val0 = (1.0/m);
  1300.                     for(k = 0; k < m; k++) {
  1301.                         seed = (seed * 214013 + 2531011);
  1302.                         val = (((seed >> 16) & 0x7fff) & 256) != 0 ? val0 : -val0;
  1303.                         At[i*astep + k] = val;
  1304.                     }
  1305.                     for(iter = 0; iter < 2; iter++) {
  1306.                         for(j = 0; j < i; j++) {
  1307.                             sd = 0;
  1308.                             for(k = 0; k < m; k++) {
  1309.                                 sd += At[i*astep + k]*At[j*astep + k];
  1310.                             }
  1311.                             asum = 0.0;
  1312.                             for(k = 0; k < m; k++) {
  1313.                                 t = (At[i*astep + k] - sd*At[j*astep + k]);
  1314.                                 At[i*astep + k] = t;
  1315.                                 asum += Math.abs(t);
  1316.                             }
  1317.                             asum = asum ? 1.0/asum : 0;
  1318.                             for(k = 0; k < m; k++) {
  1319.                                 At[i*astep + k] *= asum;
  1320.                             }
  1321.                         }
  1322.                     }
  1323.                     sd = 0;
  1324.                     for(k = 0; k < m; k++) {
  1325.                         t = At[i*astep + k];
  1326.                         sd += t*t;
  1327.                     }
  1328.                     sd = Math.sqrt(sd);
  1329.                 }
  1330.                
  1331.                 s = (1.0/sd);
  1332.                 for(k = 0; k < m; k++) {
  1333.                     At[i*astep + k] *= s;
  1334.                 }
  1335.             }
  1336.  
  1337.             jsfeat.cache.put_buffer(W_buff);
  1338.         }
  1339.        
  1340.         return {
  1341.  
  1342.             lu_solve: function(A, B) {
  1343.                 var i=0,j=0,k=0,p=1,astep=A.cols;
  1344.                 var ad=A.data, bd=B.data;
  1345.                 var t,alpha,d,s;
  1346.  
  1347.                 for(i = 0; i < astep; i++) {
  1348.                     k = i;                    
  1349.                     for(j = i+1; j < astep; j++) {
  1350.                         if(Math.abs(ad[j*astep + i]) > Math.abs(ad[k*astep+i])) {
  1351.                             k = j;
  1352.                         }
  1353.                     }
  1354.                    
  1355.                     if(Math.abs(ad[k*astep+i]) < jsfeat.EPSILON) {
  1356.                         return 0; // FAILED
  1357.                     }
  1358.                    
  1359.                     if(k != i) {
  1360.                         for(j = i; j < astep; j++ ) {
  1361.                             swap(ad, i*astep+j, k*astep+j, t);
  1362.                         }
  1363.                        
  1364.                         swap(bd, i, k, t);
  1365.                         p = -p;
  1366.                     }
  1367.                    
  1368.                     d = -1.0/ad[i*astep+i];
  1369.                    
  1370.                     for(j = i+1; j < astep; j++) {
  1371.                         alpha = ad[j*astep+i]*d;
  1372.                        
  1373.                         for(k = i+1; k < astep; k++) {
  1374.                             ad[j*astep+k] += alpha*ad[i*astep+k];
  1375.                         }
  1376.                        
  1377.                         bd[j] += alpha*bd[i];
  1378.                     }
  1379.                    
  1380.                     ad[i*astep+i] = -d;
  1381.                 }
  1382.                
  1383.                 for(i = astep-1; i >= 0; i--) {
  1384.                     s = bd[i];
  1385.                     for(k = i+1; k < astep; k++) {
  1386.                         s -= ad[i*astep+k]*bd[k];
  1387.                     }
  1388.                     bd[i] = s*ad[i*astep+i];
  1389.                 }
  1390.  
  1391.                 return 1; // OK
  1392.             },
  1393.  
  1394.             cholesky_solve: function(A, B) {
  1395.                 var col=0,row=0,col2=0,cs=0,rs=0,i=0,j=0;
  1396.                 var size = A.cols;
  1397.                 var ad=A.data, bd=B.data;
  1398.                 var val,inv_diag;
  1399.  
  1400.                 for (col = 0; col < size; col++) {
  1401.                     inv_diag = 1.0;
  1402.                     cs = (col * size);
  1403.                     rs = cs;
  1404.                     for (row = col; row < size; row++)
  1405.                     {
  1406.                         // correct for the parts of cholesky already computed
  1407.                         val = ad[(rs+col)];
  1408.                         for (col2 = 0; col2 < col; col2++) {
  1409.                             val -= ad[(col2*size+col)] * ad[(rs+col2)];
  1410.                         }
  1411.                         if (row == col) {
  1412.                             // this is the diagonal element so don't divide
  1413.                             ad[(rs+col)] = val;
  1414.                             if(val == 0) {
  1415.                                 return 0;
  1416.                             }
  1417.                             inv_diag = 1.0 / val;
  1418.                         } else {
  1419.                             // cache the value without division in the upper half
  1420.                             ad[(cs+row)] = val;
  1421.                             // divide my the diagonal element for all others
  1422.                             ad[(rs+col)] = val * inv_diag;
  1423.                         }
  1424.                         rs = (rs + size);
  1425.                     }
  1426.                 }
  1427.  
  1428.                 // first backsub through L
  1429.                 cs = 0;
  1430.                 for (i = 0; i < size; i++) {
  1431.                     val = bd[i];
  1432.                     for (j = 0; j < i; j++) {
  1433.                         val -= ad[(cs+j)] * bd[j];
  1434.                     }
  1435.                     bd[i] = val;
  1436.                     cs = (cs + size);
  1437.                 }
  1438.                 // backsub through diagonal
  1439.                 cs = 0;
  1440.                 for (i = 0; i < size; i++) {
  1441.                     bd[i] /= ad[(cs + i)];
  1442.                     cs = (cs + size);
  1443.                 }
  1444.                 // backsub through L Transpose
  1445.                 i = (size-1);
  1446.                 for (; i >= 0; i--) {
  1447.                     val = bd[i];
  1448.                     j = (i + 1);
  1449.                     cs = (j * size);
  1450.                     for (; j < size; j++) {
  1451.                         val -= ad[(cs + i)] * bd[j];
  1452.                         cs = (cs + size);
  1453.                     }
  1454.                     bd[i] = val;
  1455.                 }
  1456.  
  1457.                 return 1;
  1458.             },
  1459.  
  1460.             svd_decompose: function(A, W, U, V, options) {
  1461.                 if (typeof options === "undefined") { options = 0; };
  1462.                 var at=0,i=0,j=0,_m=A.rows,_n=A.cols,m=_m,n=_n;
  1463.                 var dt = A.type | jsfeat.C1_t; // we only work with single channel
  1464.  
  1465.                 if(m < n) {
  1466.                     at = 1;
  1467.                     i = m;
  1468.                     m = n;
  1469.                     n = i;
  1470.                 }
  1471.  
  1472.                 var a_buff = jsfeat.cache.get_buffer((m*m)<<3);
  1473.                 var w_buff = jsfeat.cache.get_buffer(n<<3);
  1474.                 var v_buff = jsfeat.cache.get_buffer((n*n)<<3);
  1475.  
  1476.                 var a_mt = new jsfeat.matrix_t(m, m, dt, a_buff.data);
  1477.                 var w_mt = new jsfeat.matrix_t(1, n, dt, w_buff.data);
  1478.                 var v_mt = new jsfeat.matrix_t(n, n, dt, v_buff.data);
  1479.  
  1480.                 if(at == 0) {
  1481.                     // transpose
  1482.                     jsfeat.matmath.transpose(a_mt, A);
  1483.                 } else {
  1484.                     for(i = 0; i < _n*_m; i++) {
  1485.                         a_mt.data[i] = A.data[i];
  1486.                     }
  1487.                     for(; i < n*m; i++) {
  1488.                         a_mt.data[i] = 0;
  1489.                     }
  1490.                 }
  1491.  
  1492.                 JacobiSVDImpl(a_mt.data, m, w_mt.data, v_mt.data, n, m, n, m);
  1493.  
  1494.                 if(W) {
  1495.                     for(i=0; i < n; i++) {
  1496.                         W.data[i] = w_mt.data[i];
  1497.                     }
  1498.                     for(; i < _n; i++) {
  1499.                         W.data[i] = 0;
  1500.                     }
  1501.                 }
  1502.  
  1503.                 if (at == 0) {
  1504.                     if(U && (options & jsfeat.SVD_U_T)) {
  1505.                         i = m*m;
  1506.                         while(--i >= 0) {
  1507.                             U.data[i] = a_mt.data[i];
  1508.                         }
  1509.                     } else if(U) {
  1510.                         jsfeat.matmath.transpose(U, a_mt);
  1511.                     }
  1512.  
  1513.                     if(V && (options & jsfeat.SVD_V_T)) {
  1514.                         i = n*n;
  1515.                         while(--i >= 0) {
  1516.                             V.data[i] = v_mt.data[i];
  1517.                         }
  1518.                     } else if(V) {
  1519.                         jsfeat.matmath.transpose(V, v_mt);
  1520.                     }
  1521.                 } else {
  1522.                     if(U && (options & jsfeat.SVD_U_T)) {
  1523.                         i = n*n;
  1524.                         while(--i >= 0) {
  1525.                             U.data[i] = v_mt.data[i];
  1526.                         }
  1527.                     } else if(U) {
  1528.                         jsfeat.matmath.transpose(U, v_mt);
  1529.                     }
  1530.  
  1531.                     if(V && (options & jsfeat.SVD_V_T)) {
  1532.                         i = m*m;
  1533.                         while(--i >= 0) {
  1534.                             V.data[i] = a_mt.data[i];
  1535.                         }
  1536.                     } else if(V) {
  1537.                         jsfeat.matmath.transpose(V, a_mt);
  1538.                     }
  1539.                 }
  1540.  
  1541.                 jsfeat.cache.put_buffer(a_buff);
  1542.                 jsfeat.cache.put_buffer(w_buff);
  1543.                 jsfeat.cache.put_buffer(v_buff);
  1544.  
  1545.             },
  1546.  
  1547.             svd_solve: function(A, X, B) {
  1548.                 var i=0,j=0,k=0;
  1549.                 var pu=0,pv=0;
  1550.                 var nrows=A.rows,ncols=A.cols;
  1551.                 var sum=0.0,xsum=0.0,tol=0.0;
  1552.                 var dt = A.type | jsfeat.C1_t;
  1553.  
  1554.                 var u_buff = jsfeat.cache.get_buffer((nrows*nrows)<<3);
  1555.                 var w_buff = jsfeat.cache.get_buffer(ncols<<3);
  1556.                 var v_buff = jsfeat.cache.get_buffer((ncols*ncols)<<3);
  1557.  
  1558.                 var u_mt = new jsfeat.matrix_t(nrows, nrows, dt, u_buff.data);
  1559.                 var w_mt = new jsfeat.matrix_t(1, ncols, dt, w_buff.data);
  1560.                 var v_mt = new jsfeat.matrix_t(ncols, ncols, dt, v_buff.data);
  1561.  
  1562.                 var bd = B.data, ud = u_mt.data, wd = w_mt.data, vd = v_mt.data;
  1563.  
  1564.                 this.svd_decompose(A, w_mt, u_mt, v_mt, 0);
  1565.  
  1566.                 tol = jsfeat.EPSILON * wd[0] * ncols;
  1567.  
  1568.                 for (; i < ncols; i++, pv += ncols) {
  1569.                     xsum = 0.0;
  1570.                     for(j = 0; j < ncols; j++) {
  1571.                         if(wd[j] > tol) {
  1572.                             for(k = 0, sum = 0.0, pu = 0; k < nrows; k++, pu += ncols) {
  1573.                                 sum += ud[pu + j] * bd[k];
  1574.                             }
  1575.                             xsum += sum * vd[pv + j] / wd[j];
  1576.                         }
  1577.                     }
  1578.                     X.data[i] = xsum;
  1579.                 }
  1580.  
  1581.                 jsfeat.cache.put_buffer(u_buff);
  1582.                 jsfeat.cache.put_buffer(w_buff);
  1583.                 jsfeat.cache.put_buffer(v_buff);
  1584.             },
  1585.  
  1586.             svd_invert: function(Ai, A) {
  1587.                 var i=0,j=0,k=0;
  1588.                 var pu=0,pv=0,pa=0;
  1589.                 var nrows=A.rows,ncols=A.cols;
  1590.                 var sum=0.0,tol=0.0;
  1591.                 var dt = A.type | jsfeat.C1_t;
  1592.  
  1593.                 var u_buff = jsfeat.cache.get_buffer((nrows*nrows)<<3);
  1594.                 var w_buff = jsfeat.cache.get_buffer(ncols<<3);
  1595.                 var v_buff = jsfeat.cache.get_buffer((ncols*ncols)<<3);
  1596.  
  1597.                 var u_mt = new jsfeat.matrix_t(nrows, nrows, dt, u_buff.data);
  1598.                 var w_mt = new jsfeat.matrix_t(1, ncols, dt, w_buff.data);
  1599.                 var v_mt = new jsfeat.matrix_t(ncols, ncols, dt, v_buff.data);
  1600.  
  1601.                 var id = Ai.data, ud = u_mt.data, wd = w_mt.data, vd = v_mt.data;
  1602.  
  1603.                 this.svd_decompose(A, w_mt, u_mt, v_mt, 0);
  1604.  
  1605.                 tol = jsfeat.EPSILON * wd[0] * ncols;
  1606.  
  1607.                 for (; i < ncols; i++, pv += ncols) {
  1608.                     for (j = 0, pu = 0; j < nrows; j++, pa++) {
  1609.                         for (k = 0, sum = 0.0; k < ncols; k++, pu++) {
  1610.                             if (wd[k] > tol) sum += vd[pv + k] * ud[pu] / wd[k];
  1611.                         }
  1612.                         id[pa] = sum;
  1613.                     }
  1614.                 }
  1615.  
  1616.                 jsfeat.cache.put_buffer(u_buff);
  1617.                 jsfeat.cache.put_buffer(w_buff);
  1618.                 jsfeat.cache.put_buffer(v_buff);
  1619.             },
  1620.  
  1621.             eigenVV: function(A, vects, vals) {
  1622.                 var n=A.cols,i=n*n;
  1623.                 var dt = A.type | jsfeat.C1_t;
  1624.  
  1625.                 var a_buff = jsfeat.cache.get_buffer((n*n)<<3);
  1626.                 var w_buff = jsfeat.cache.get_buffer(n<<3);
  1627.                 var a_mt = new jsfeat.matrix_t(n, n, dt, a_buff.data);
  1628.                 var w_mt = new jsfeat.matrix_t(1, n, dt, w_buff.data);
  1629.  
  1630.                 while(--i >= 0) {
  1631.                     a_mt.data[i] = A.data[i];
  1632.                 }
  1633.  
  1634.                 JacobiImpl(a_mt.data, n, w_mt.data, vects ? vects.data : null, n, n);
  1635.  
  1636.                 if(vals) {
  1637.                     while(--n >= 0) {
  1638.                         vals.data[n] = w_mt.data[n];
  1639.                     }
  1640.                 }
  1641.  
  1642.                 jsfeat.cache.put_buffer(a_buff);
  1643.                 jsfeat.cache.put_buffer(w_buff);
  1644.             }
  1645.  
  1646.         };
  1647.  
  1648.     })();
  1649.  
  1650.     global.linalg = linalg;
  1651.  
  1652. })(jsfeat);
  1653. /**
  1654.  * @author Eugene Zatepyakin / http://inspirit.ru/
  1655.  *
  1656.  */
  1657.  
  1658. (function(global) {
  1659.     "use strict";
  1660.     //
  1661.  
  1662.     var motion_model = (function() {
  1663.  
  1664.         var sqr = function(x) {
  1665.                 return x*x;
  1666.         }
  1667.  
  1668.         // does isotropic normalization
  1669.         var iso_normalize_points = function(from, to, T0, T1, count) {
  1670.                         var i=0;
  1671.                     var cx0=0.0, cy0=0.0, d0=0.0, s0=0.0;
  1672.                     var cx1=0.0, cy1=0.0, d1=0.0, s1=0.0;
  1673.                     var dx=0.0,dy=0.0;
  1674.  
  1675.                     for (; i < count; ++i) {
  1676.                         cx0 += from[i].x;
  1677.                         cy0 += from[i].y;
  1678.                         cx1 += to[i].x;
  1679.                         cy1 += to[i].y;
  1680.                     }
  1681.  
  1682.                     cx0 /= count; cy0 /= count;
  1683.                     cx1 /= count; cy1 /= count;
  1684.  
  1685.                     for (i = 0; i < count; ++i) {
  1686.                         dx = from[i].x - cx0;
  1687.                         dy = from[i].y - cy0;
  1688.                         d0 += Math.sqrt(dx*dx + dy*dy);
  1689.                         dx = to[i].x - cx1;
  1690.                         dy = to[i].y - cy1;
  1691.                         d1 += Math.sqrt(dx*dx + dy*dy);
  1692.                     }
  1693.  
  1694.                     d0 /= count; d1 /= count;
  1695.  
  1696.                     s0 = Math.SQRT2 / d0; s1 = Math.SQRT2 / d1;
  1697.  
  1698.                     T0[0] = T0[4] = s0;
  1699.                     T0[2] = -cx0*s0;
  1700.                     T0[5] = -cy0*s0;
  1701.                     T0[1] = T0[3] = T0[6] = T0[7] = 0.0;
  1702.                     T0[8] = 1.0;
  1703.  
  1704.                     T1[0] = T1[4] = s1;
  1705.                     T1[2] = -cx1*s1;
  1706.                     T1[5] = -cy1*s1;
  1707.                     T1[1] = T1[3] = T1[6] = T1[7] = 0.0;
  1708.                     T1[8] = 1.0;
  1709.                 }
  1710.  
  1711.                 var have_collinear_points = function(points, count) {
  1712.                     var j=0,k=0,i=(count-1)|0;
  1713.                     var dx1=0.0,dy1=0.0,dx2=0.0,dy2=0.0;
  1714.  
  1715.                     // check that the i-th selected point does not belong
  1716.                     // to a line connecting some previously selected points
  1717.                     for(; j < i; ++j) {
  1718.                         dx1 = points[j].x - points[i].x;
  1719.                         dy1 = points[j].y - points[i].y;
  1720.                         for(k = 0; k < j; ++k) {
  1721.                             dx2 = points[k].x - points[i].x;
  1722.                             dy2 = points[k].y - points[i].y;
  1723.                             if( Math.abs(dx2*dy1 - dy2*dx1) <= jsfeat.EPSILON*(Math.abs(dx1) + Math.abs(dy1) + Math.abs(dx2) + Math.abs(dy2)))
  1724.                                 return true;
  1725.                         }
  1726.                     }
  1727.                     return false;
  1728.                 }
  1729.  
  1730.                 var T0 = new jsfeat.matrix_t(3, 3, jsfeat.F32_t|jsfeat.C1_t);
  1731.         var T1 = new jsfeat.matrix_t(3, 3, jsfeat.F32_t|jsfeat.C1_t);
  1732.         var AtA = new jsfeat.matrix_t(6, 6, jsfeat.F32_t|jsfeat.C1_t);
  1733.         var AtB = new jsfeat.matrix_t(6, 1, jsfeat.F32_t|jsfeat.C1_t);
  1734.        
  1735.         var affine2d = (function () {
  1736.  
  1737.                 function affine2d() {
  1738.                         // empty constructor
  1739.                 }
  1740.  
  1741.                 affine2d.prototype.run = function(from, to, model, count) {
  1742.                         var i=0,j=0;
  1743.                         var dt=model.type|jsfeat.C1_t;
  1744.                         var md=model.data, t0d=T0.data, t1d=T1.data;
  1745.                         var pt0,pt1,px=0.0,py=0.0;
  1746.  
  1747.                     iso_normalize_points(from, to, t0d, t1d, count);
  1748.  
  1749.                     var a_buff = jsfeat.cache.get_buffer((2*count*6)<<3);
  1750.                 var b_buff = jsfeat.cache.get_buffer((2*count)<<3);
  1751.  
  1752.                 var a_mt = new jsfeat.matrix_t(6, 2*count, dt, a_buff.data);
  1753.                 var b_mt = new jsfeat.matrix_t(1, 2*count, dt, b_buff.data);
  1754.                 var ad=a_mt.data, bd=b_mt.data;
  1755.  
  1756.                             for (; i < count; ++i) {
  1757.                                 pt0 = from[i];
  1758.                                 pt1 = to[i];
  1759.  
  1760.                                 px = t0d[0]*pt0.x + t0d[1]*pt0.y + t0d[2];
  1761.                                 py = t0d[3]*pt0.x + t0d[4]*pt0.y + t0d[5];
  1762.  
  1763.                                 j = i*2*6;
  1764.                                 ad[j]=px, ad[j+1]=py, ad[j+2]=1.0, ad[j+3]=0.0, ad[j+4]=0.0, ad[j+5]=0.0;
  1765.  
  1766.                                 j += 6;
  1767.                                 ad[j]=0.0, ad[j+1]=0.0, ad[j+2]=0.0, ad[j+3]=px, ad[j+4]=py, ad[j+5]=1.0;
  1768.  
  1769.                                 bd[i<<1] = t1d[0]*pt1.x + t1d[1]*pt1.y + t1d[2];
  1770.                                 bd[(i<<1)+1] = t1d[3]*pt1.x + t1d[4]*pt1.y + t1d[5];
  1771.                             }
  1772.  
  1773.                             jsfeat.matmath.multiply_AtA(AtA, a_mt);
  1774.                             jsfeat.matmath.multiply_AtB(AtB, a_mt, b_mt);
  1775.  
  1776.                             jsfeat.linalg.lu_solve(AtA, AtB);
  1777.  
  1778.                             md[0] = AtB.data[0], md[1]=AtB.data[1], md[2]=AtB.data[2];
  1779.                             md[3] = AtB.data[3], md[4]=AtB.data[4], md[5]=AtB.data[5];
  1780.                             md[6] = 0.0, md[7] = 0.0, md[8] = 1.0; // fill last row
  1781.  
  1782.                             // denormalize
  1783.                             jsfeat.matmath.invert_3x3(T1, T1);
  1784.                             jsfeat.matmath.multiply_3x3(model, T1, model);
  1785.                             jsfeat.matmath.multiply_3x3(model, model, T0);
  1786.  
  1787.                             // free buffer
  1788.                             jsfeat.cache.put_buffer(a_buff);
  1789.                             jsfeat.cache.put_buffer(b_buff);
  1790.  
  1791.                             return 1;
  1792.                 }
  1793.  
  1794.                 affine2d.prototype.error = function(from, to, model, err, count) {
  1795.                         var i=0;
  1796.                         var pt0,pt1;
  1797.                         var m=model.data;
  1798.  
  1799.                             for (; i < count; ++i) {
  1800.                                 pt0 = from[i];
  1801.                                 pt1 = to[i];
  1802.  
  1803.                                 err[i] = sqr(pt1.x - m[0]*pt0.x - m[1]*pt0.y - m[2]) +
  1804.                                          sqr(pt1.y - m[3]*pt0.x - m[4]*pt0.y - m[5]);
  1805.                             }
  1806.                 }
  1807.  
  1808.                 affine2d.prototype.check_subset = function(from, to, count) {
  1809.                     return true; // all good
  1810.                 }
  1811.  
  1812.                 return affine2d;
  1813.             })();
  1814.  
  1815.             var mLtL = new jsfeat.matrix_t(9, 9, jsfeat.F32_t|jsfeat.C1_t);
  1816.             var Evec = new jsfeat.matrix_t(9, 9, jsfeat.F32_t|jsfeat.C1_t);
  1817.  
  1818.             var homography2d = (function () {
  1819.  
  1820.                 function homography2d() {
  1821.                         // empty constructor
  1822.                         //this.T0 = new jsfeat.matrix_t(3, 3, jsfeat.F32_t|jsfeat.C1_t);
  1823.                         //this.T1 = new jsfeat.matrix_t(3, 3, jsfeat.F32_t|jsfeat.C1_t);
  1824.                         //this.mLtL = new jsfeat.matrix_t(9, 9, jsfeat.F32_t|jsfeat.C1_t);
  1825.                         //this.Evec = new jsfeat.matrix_t(9, 9, jsfeat.F32_t|jsfeat.C1_t);
  1826.                 }
  1827.  
  1828.                 homography2d.prototype.run = function(from, to, model, count) {
  1829.                         var i=0,j=0;
  1830.                         var md=model.data, t0d=T0.data, t1d=T1.data;
  1831.                         var LtL=mLtL.data, evd=Evec.data;
  1832.                         var x=0.0,y=0.0,X=0.0,Y=0.0;
  1833.  
  1834.                             // norm
  1835.                                 var smx=0.0, smy=0.0, cmx=0.0, cmy=0.0, sMx=0.0, sMy=0.0, cMx=0.0, cMy=0.0;
  1836.  
  1837.                                 for(; i < count; ++i) {
  1838.                                     cmx += to[i].x;
  1839.                                     cmy += to[i].y;
  1840.                                     cMx += from[i].x;
  1841.                                     cMy += from[i].y;
  1842.                                 }
  1843.  
  1844.                             cmx /= count; cmy /= count;
  1845.                             cMx /= count; cMy /= count;
  1846.  
  1847.                             for(i = 0; i < count; ++i)
  1848.                             {
  1849.                                     smx += Math.abs(to[i].x - cmx);
  1850.                                     smy += Math.abs(to[i].y - cmy);
  1851.                                     sMx += Math.abs(from[i].x - cMx);
  1852.                                     sMy += Math.abs(from[i].y - cMy);
  1853.                                 }
  1854.  
  1855.                             if( Math.abs(smx) < jsfeat.EPSILON
  1856.                                 || Math.abs(smy) < jsfeat.EPSILON
  1857.                                 || Math.abs(sMx) < jsfeat.EPSILON
  1858.                                 || Math.abs(sMy) < jsfeat.EPSILON ) return 0;
  1859.  
  1860.                             smx = count/smx; smy = count/smy;
  1861.                             sMx = count/sMx; sMy = count/sMy;
  1862.  
  1863.                             t0d[0] = sMx;       t0d[1] = 0;     t0d[2] = -cMx*sMx;
  1864.                             t0d[3] = 0;         t0d[4] = sMy;   t0d[5] = -cMy*sMy;
  1865.                             t0d[6] = 0;         t0d[7] = 0;     t0d[8] = 1;
  1866.  
  1867.                                 t1d[0] = 1.0/smx;       t1d[1] = 0;             t1d[2] = cmx;
  1868.                                 t1d[3] = 0;             t1d[4] = 1.0/smy;       t1d[5] = cmy;
  1869.                                 t1d[6] = 0;             t1d[7] = 0;             t1d[8] = 1;
  1870.                                 //
  1871.  
  1872.                                 // construct system
  1873.                                 i = 81;
  1874.                                 while(--i >= 0) {
  1875.                                         LtL[i] = 0.0;
  1876.                                 }
  1877.                                 for(i = 0; i < count; ++i) {
  1878.                                         x = (to[i].x - cmx) * smx;
  1879.                                         y = (to[i].y - cmy) * smy;
  1880.                                         X = (from[i].x - cMx) * sMx;
  1881.                                         Y = (from[i].y - cMy) * sMy;
  1882.  
  1883.                                         LtL[0] += X*X;
  1884.                                         LtL[1] += X*Y;
  1885.                                         LtL[2] += X;
  1886.  
  1887.                                         LtL[6] += X*-x*X;
  1888.                                         LtL[7] += X*-x*Y;
  1889.                                         LtL[8] += X*-x;
  1890.                                         LtL[10] += Y*Y;
  1891.                                         LtL[11] += Y;
  1892.  
  1893.                                         LtL[15] += Y*-x*X;
  1894.                                         LtL[16] += Y*-x*Y;
  1895.                                         LtL[17] += Y*-x;
  1896.                                         LtL[20] += 1.0;
  1897.  
  1898.                                         LtL[24] += -x*X;
  1899.                                         LtL[25] += -x*Y;
  1900.                                         LtL[26] += -x;
  1901.                                         LtL[30] += X*X;
  1902.                                         LtL[31] += X*Y;
  1903.                                         LtL[32] += X;
  1904.                                         LtL[33] += X*-y*X;
  1905.                                         LtL[34] += X*-y*Y;
  1906.                                         LtL[35] += X*-y;
  1907.                                         LtL[40] += Y*Y;
  1908.                                         LtL[41] += Y;
  1909.                                         LtL[42] += Y*-y*X;
  1910.                                         LtL[43] += Y*-y*Y;
  1911.                                         LtL[44] += Y*-y;
  1912.                                         LtL[50] += 1.0;
  1913.                                         LtL[51] += -y*X;
  1914.                                         LtL[52] += -y*Y;
  1915.                                         LtL[53] += -y;
  1916.                                         LtL[60] += -x*X*-x*X + -y*X*-y*X;
  1917.                                         LtL[61] += -x*X*-x*Y + -y*X*-y*Y;
  1918.                                         LtL[62] += -x*X*-x + -y*X*-y;
  1919.                                         LtL[70] += -x*Y*-x*Y + -y*Y*-y*Y;
  1920.                                         LtL[71] += -x*Y*-x + -y*Y*-y;
  1921.                                         LtL[80] += -x*-x + -y*-y;
  1922.                                 }
  1923.                                 //
  1924.  
  1925.                                 // symmetry
  1926.                             for(i = 0; i < 9; ++i) {
  1927.                                 for(j = 0; j < i; ++j)
  1928.                                     LtL[i*9+j] = LtL[j*9+i];
  1929.                             }
  1930.  
  1931.                                 jsfeat.linalg.eigenVV(mLtL, Evec);
  1932.  
  1933.                                 md[0]=evd[72], md[1]=evd[73], md[2]=evd[74];
  1934.                             md[3]=evd[75], md[4]=evd[76], md[5]=evd[77];
  1935.                             md[6]=evd[78], md[7]=evd[79], md[8]=evd[80];
  1936.  
  1937.                                 // denormalize
  1938.                             jsfeat.matmath.multiply_3x3(model, T1, model);
  1939.                             jsfeat.matmath.multiply_3x3(model, model, T0);
  1940.  
  1941.                             // set bottom right to 1.0
  1942.                             x = 1.0/md[8];
  1943.                             md[0] *= x; md[1] *= x; md[2] *= x;
  1944.                             md[3] *= x; md[4] *= x; md[5] *= x;
  1945.                             md[6] *= x; md[7] *= x; md[8] = 1.0;
  1946.  
  1947.                             return 1;
  1948.                 }
  1949.  
  1950.                 homography2d.prototype.error = function(from, to, model, err, count) {
  1951.                         var i=0;
  1952.                         var pt0,pt1,ww=0.0,dx=0.0,dy=0.0;
  1953.                         var m=model.data;
  1954.  
  1955.                             for (; i < count; ++i) {
  1956.                                 pt0 = from[i];
  1957.                                 pt1 = to[i];
  1958.  
  1959.                                 ww = 1.0/(m[6]*pt0.x + m[7]*pt0.y + 1.0);
  1960.                                 dx = (m[0]*pt0.x + m[1]*pt0.y + m[2])*ww - pt1.x;
  1961.                                 dy = (m[3]*pt0.x + m[4]*pt0.y + m[5])*ww - pt1.y;
  1962.                                 err[i] = (dx*dx + dy*dy);
  1963.                             }
  1964.                 }
  1965.  
  1966.                 homography2d.prototype.check_subset = function(from, to, count) {
  1967.                         // seems to reject good subsets actually
  1968.                         //if( have_collinear_points(from, count) || have_collinear_points(to, count) ) {
  1969.                                 //return false;
  1970.                         //}
  1971.                         if( count == 4 ) {
  1972.                                 var negative = 0;
  1973.  
  1974.                                 var fp0=from[0],fp1=from[1],fp2=from[2],fp3=from[3];
  1975.                                 var tp0=to[0],tp1=to[1],tp2=to[2],tp3=to[3];
  1976.  
  1977.                                 // set1
  1978.                                 var A11=fp0.x, A12=fp0.y, A13=1.0;
  1979.                                 var A21=fp1.x, A22=fp1.y, A23=1.0;
  1980.                                 var A31=fp2.x, A32=fp2.y, A33=1.0;
  1981.  
  1982.                                 var B11=tp0.x, B12=tp0.y, B13=1.0;
  1983.                                 var B21=tp1.x, B22=tp1.y, B23=1.0;
  1984.                                 var B31=tp2.x, B32=tp2.y, B33=1.0;
  1985.  
  1986.                                 var detA = jsfeat.matmath.determinant_3x3(A11,A12,A13, A21,A22,A23, A31,A32,A33);
  1987.                                         var detB = jsfeat.matmath.determinant_3x3(B11,B12,B13, B21,B22,B23, B31,B32,B33);
  1988.  
  1989.                                         if(detA*detB < 0) negative++;
  1990.  
  1991.                                         // set2
  1992.                                         A11=fp1.x, A12=fp1.y;
  1993.                                 A21=fp2.x, A22=fp2.y;
  1994.                                 A31=fp3.x, A32=fp3.y;
  1995.  
  1996.                                 B11=tp1.x, B12=tp1.y;
  1997.                                 B21=tp2.x, B22=tp2.y;
  1998.                                 B31=tp3.x, B32=tp3.y;
  1999.  
  2000.                                 detA = jsfeat.matmath.determinant_3x3(A11,A12,A13, A21,A22,A23, A31,A32,A33);
  2001.                                         detB = jsfeat.matmath.determinant_3x3(B11,B12,B13, B21,B22,B23, B31,B32,B33);
  2002.  
  2003.                                         if(detA*detB < 0) negative++;
  2004.  
  2005.                                         // set3
  2006.                                         A11=fp0.x, A12=fp0.y;
  2007.                                 A21=fp2.x, A22=fp2.y;
  2008.                                 A31=fp3.x, A32=fp3.y;
  2009.  
  2010.                                 B11=tp0.x, B12=tp0.y;
  2011.                                 B21=tp2.x, B22=tp2.y;
  2012.                                 B31=tp3.x, B32=tp3.y;
  2013.  
  2014.                                 detA = jsfeat.matmath.determinant_3x3(A11,A12,A13, A21,A22,A23, A31,A32,A33);
  2015.                                         detB = jsfeat.matmath.determinant_3x3(B11,B12,B13, B21,B22,B23, B31,B32,B33);
  2016.  
  2017.                                         if(detA*detB < 0) negative++;
  2018.  
  2019.                                         // set4
  2020.                                         A11=fp0.x, A12=fp0.y;
  2021.                                 A21=fp1.x, A22=fp1.y;
  2022.                                 A31=fp3.x, A32=fp3.y;
  2023.  
  2024.                                 B11=tp0.x, B12=tp0.y;
  2025.                                 B21=tp1.x, B22=tp1.y;
  2026.                                 B31=tp3.x, B32=tp3.y;
  2027.  
  2028.                                 detA = jsfeat.matmath.determinant_3x3(A11,A12,A13, A21,A22,A23, A31,A32,A33);
  2029.                                         detB = jsfeat.matmath.determinant_3x3(B11,B12,B13, B21,B22,B23, B31,B32,B33);
  2030.  
  2031.                                         if(detA*detB < 0) negative++;
  2032.  
  2033.                                 if(negative != 0 && negative != 4) {
  2034.                                         return false;
  2035.                                 }
  2036.                             }
  2037.                     return true; // all good
  2038.                 }
  2039.  
  2040.                 return homography2d;
  2041.             })();
  2042.  
  2043.             return {
  2044.  
  2045.                 affine2d:affine2d,
  2046.                 homography2d:homography2d
  2047.  
  2048.         };
  2049.  
  2050.     })();
  2051.  
  2052.     var ransac_params_t = (function () {
  2053.         function ransac_params_t(size, thresh, eps, prob) {
  2054.             if (typeof size === "undefined") { size=0; }
  2055.             if (typeof thresh === "undefined") { thresh=0.5; }
  2056.             if (typeof eps === "undefined") { eps=0.5; }
  2057.             if (typeof prob === "undefined") { prob=0.99; }
  2058.  
  2059.             this.size = size;
  2060.             this.thresh = thresh;
  2061.             this.eps = eps;
  2062.             this.prob = prob;
  2063.         };
  2064.         ransac_params_t.prototype.update_iters = function(_eps, max_iters) {
  2065.                 var num = Math.log(1 - this.prob);
  2066.                 var denom = Math.log(1 - Math.pow(1 - _eps, this.size));
  2067.                 return (denom >= 0 || -num >= max_iters*(-denom) ? max_iters : Math.round(num/denom))|0;
  2068.         };
  2069.         return ransac_params_t;
  2070.     })();
  2071.  
  2072.     var motion_estimator = (function() {
  2073.  
  2074.         var get_subset = function(kernel, from, to, need_cnt, max_cnt, from_sub, to_sub) {
  2075.                 var max_try = 1000;
  2076.                 var indices = [];
  2077.                     var i=0, j=0, ssiter=0, idx_i=0, ok=false;
  2078.                     for(; ssiter < max_try; ++ssiter)  {
  2079.                         i = 0;
  2080.                         for (; i < need_cnt && ssiter < max_try;) {
  2081.                             ok = false;
  2082.                             idx_i = 0;
  2083.                             while (!ok) {
  2084.                                 ok = true;
  2085.                                 idx_i = indices[i] = Math.floor(Math.random() * max_cnt)|0;
  2086.                                 for (j = 0; j < i; ++j) {
  2087.                                     if (idx_i == indices[j])
  2088.                                     { ok = false; break; }
  2089.                                 }
  2090.                             }
  2091.                             from_sub[i] = from[idx_i];
  2092.                             to_sub[i] = to[idx_i];
  2093.                             if( !kernel.check_subset( from_sub, to_sub, i+1 ) ) {
  2094.                                 ssiter++;
  2095.                                 continue;
  2096.                             }
  2097.                             ++i;
  2098.                         }
  2099.                         break;
  2100.                     }
  2101.  
  2102.                     return (i == need_cnt && ssiter < max_try);
  2103.         }
  2104.  
  2105.         var find_inliers = function(kernel, model, from, to, count, thresh, err, mask) {
  2106.                 var numinliers = 0, i=0, f=0;
  2107.                 var t = thresh*thresh;
  2108.  
  2109.                 kernel.error(from, to, model, err, count);
  2110.  
  2111.                     for(; i < count; ++i) {
  2112.                         f = err[i] <= t;
  2113.                         mask[i] = f;
  2114.                         numinliers += f;
  2115.                     }
  2116.                     return numinliers;
  2117.         }
  2118.  
  2119.         return {
  2120.  
  2121.                 ransac: function(params, kernel, from, to, count, model, mask, max_iters) {
  2122.                         if (typeof max_iters === "undefined") { max_iters=1000; }
  2123.  
  2124.                         if(count < params.size) return false;
  2125.  
  2126.                         var model_points = params.size;
  2127.                             var niters = max_iters, iter=0;
  2128.                             var result = false;
  2129.  
  2130.                             var subset0 = [];
  2131.                             var subset1 = [];
  2132.                             var found = false;
  2133.  
  2134.                             var mc=model.cols,mr=model.rows;
  2135.                 var dt = model.type | jsfeat.C1_t;
  2136.  
  2137.                             var m_buff = jsfeat.cache.get_buffer((mc*mr)<<3);
  2138.                             var ms_buff = jsfeat.cache.get_buffer(count);
  2139.                             var err_buff = jsfeat.cache.get_buffer(count<<2);
  2140.                             var M = new jsfeat.matrix_t(mc, mr, dt, m_buff.data);
  2141.                             var curr_mask = new jsfeat.matrix_t(count, 1, jsfeat.U8C1_t, ms_buff.data);
  2142.  
  2143.                             var inliers_max = -1, numinliers=0;
  2144.                             var nmodels = 0;
  2145.  
  2146.                             var err = err_buff.f32;
  2147.  
  2148.                             // special case
  2149.                             if(count == model_points) {
  2150.                                 if(kernel.run(from, to, M, count) <= 0) {
  2151.                                         jsfeat.cache.put_buffer(m_buff);
  2152.                                         jsfeat.cache.put_buffer(ms_buff);
  2153.                                         jsfeat.cache.put_buffer(err_buff);
  2154.                                         return false;
  2155.                                 }
  2156.  
  2157.                                 M.copy_to(model);
  2158.                                 if(mask) {
  2159.                                         while(--count >= 0) {
  2160.                                                 mask.data[count] = 1;
  2161.                                         }
  2162.                                 }
  2163.                                 jsfeat.cache.put_buffer(m_buff);
  2164.                                 jsfeat.cache.put_buffer(ms_buff);
  2165.                                 jsfeat.cache.put_buffer(err_buff);
  2166.                                 return true;
  2167.                             }
  2168.  
  2169.                             for (; iter < niters; ++iter) {
  2170.                                 // generate subset
  2171.                                 found = get_subset(kernel, from, to, model_points, count, subset0, subset1);
  2172.                                 if(!found) {
  2173.                                     if(iter == 0) {
  2174.                                         jsfeat.cache.put_buffer(m_buff);
  2175.                                         jsfeat.cache.put_buffer(ms_buff);
  2176.                                         jsfeat.cache.put_buffer(err_buff);
  2177.                                         return false;
  2178.                                     }
  2179.                                     break;
  2180.                                 }
  2181.  
  2182.                                 nmodels = kernel.run( subset0, subset1, M, model_points );
  2183.                                 if(nmodels <= 0)
  2184.                                     continue;
  2185.  
  2186.                                 // TODO handle multimodel output
  2187.  
  2188.                                 numinliers = find_inliers(kernel, M, from, to, count, params.thresh, err, curr_mask.data);
  2189.  
  2190.                                 if( numinliers > Math.max(inliers_max, model_points-1) ) {
  2191.                                     M.copy_to(model);
  2192.                                     inliers_max = numinliers;
  2193.                                     if(mask) curr_mask.copy_to(mask);
  2194.                                     niters = params.update_iters((count - numinliers)/count, niters);
  2195.                                     result = true;
  2196.                                 }
  2197.                             }
  2198.  
  2199.                             jsfeat.cache.put_buffer(m_buff);
  2200.                             jsfeat.cache.put_buffer(ms_buff);
  2201.                             jsfeat.cache.put_buffer(err_buff);
  2202.  
  2203.                             return result;
  2204.                 },
  2205.  
  2206.                 lmeds: function(params, kernel, from, to, count, model, mask, max_iters) {
  2207.                         if (typeof max_iters === "undefined") { max_iters=1000; }
  2208.  
  2209.                         if(count < params.size) return false;
  2210.  
  2211.                         var model_points = params.size;
  2212.                             var niters = max_iters, iter=0;
  2213.                             var result = false;
  2214.  
  2215.                             var subset0 = [];
  2216.                             var subset1 = [];
  2217.                             var found = false;
  2218.  
  2219.                             var mc=model.cols,mr=model.rows;
  2220.                 var dt = model.type | jsfeat.C1_t;
  2221.  
  2222.                             var m_buff = jsfeat.cache.get_buffer((mc*mr)<<3);
  2223.                             var ms_buff = jsfeat.cache.get_buffer(count);