libcm is a C development framework with an emphasis on audio signal processing applications.
Você não pode selecionar mais de 25 tópicos Os tópicos devem começar com uma letra ou um número, podem incluir traços ('-') e podem ter até 35 caracteres.

cmVectOpsRICode.h 29KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218
  1. #ifdef cmVectOpsRICode_h
  2. VECT_OP_TYPE* VECT_OP_FUNC(Col)( VECT_OP_TYPE* m, unsigned ci, unsigned rn, unsigned cn )
  3. {
  4. assert(ci<cn);
  5. return m + (ci*rn);
  6. }
  7. VECT_OP_TYPE* VECT_OP_FUNC(Row)( VECT_OP_TYPE* m, unsigned ri, unsigned rn, unsigned cn )
  8. {
  9. assert(ri<rn);
  10. return m + ri;
  11. }
  12. VECT_OP_TYPE* VECT_OP_FUNC(ElePtr)( VECT_OP_TYPE* m, unsigned ri, unsigned ci, unsigned rn, unsigned cn )
  13. {
  14. assert(ri<rn && ci<cn);
  15. return m + (ci*rn) + ri;
  16. }
  17. VECT_OP_TYPE VECT_OP_FUNC(Ele)( VECT_OP_TYPE* m, unsigned ri, unsigned ci, unsigned rn, unsigned cn )
  18. { return *VECT_OP_FUNC(ElePtr)(m,ri,ci,rn,cn); }
  19. void VECT_OP_FUNC(Set)( VECT_OP_TYPE* m, unsigned ri, unsigned ci, unsigned rn, unsigned cn, VECT_OP_TYPE v )
  20. { *(VECT_OP_FUNC(ElePtr)(m,ri,ci,rn,cn)) = v; }
  21. const VECT_OP_TYPE* VECT_OP_FUNC(CCol)( const VECT_OP_TYPE* m, unsigned ci, unsigned rn, unsigned cn )
  22. {
  23. assert(ci<cn);
  24. return m + (ci*rn);
  25. }
  26. const VECT_OP_TYPE* VECT_OP_FUNC(CRow)( const VECT_OP_TYPE* m, unsigned ri, unsigned rn, unsigned cn )
  27. {
  28. assert(ri<rn);
  29. return m + ri;
  30. }
  31. const VECT_OP_TYPE* VECT_OP_FUNC(CElePtr)( const VECT_OP_TYPE* m, unsigned ri, unsigned ci, unsigned rn, unsigned cn )
  32. {
  33. assert(ri<rn && ci<cn);
  34. return m + (ci*rn) + ri;
  35. }
  36. VECT_OP_TYPE VECT_OP_FUNC(CEle)( const VECT_OP_TYPE* m, unsigned ri, unsigned ci, unsigned rn, unsigned cn )
  37. { return *VECT_OP_FUNC(CElePtr)(m,ri,ci,rn,cn); }
  38. VECT_OP_TYPE* VECT_OP_FUNC(Fill)( VECT_OP_TYPE* dbp, unsigned dn, VECT_OP_TYPE value )
  39. {
  40. const VECT_OP_TYPE* dep = dbp + dn;
  41. VECT_OP_TYPE* dp = dbp;
  42. if( value == 0 )
  43. memset(dbp,0,(dep-dbp)*sizeof(VECT_OP_TYPE));
  44. else
  45. {
  46. while( dbp < dep )
  47. *dbp++ = value;
  48. }
  49. return dp;
  50. }
  51. VECT_OP_TYPE* VECT_OP_FUNC(Zero)( VECT_OP_TYPE* dbp, unsigned dn )
  52. {
  53. memset( dbp, 0, sizeof(VECT_OP_TYPE)*dn);
  54. return dbp;
  55. }
  56. VECT_OP_TYPE* VECT_OP_FUNC(Move)( VECT_OP_TYPE* bp, unsigned bn, const VECT_OP_TYPE* sp )
  57. {
  58. memmove(bp,sp,sizeof(VECT_OP_TYPE)*bn);
  59. return bp;
  60. }
  61. VECT_OP_TYPE* VECT_OP_FUNC(Copy)( VECT_OP_TYPE* bp, unsigned bn, const VECT_OP_TYPE* sp )
  62. {
  63. memcpy(bp,sp,sizeof(VECT_OP_TYPE)*bn);
  64. return bp;
  65. }
  66. VECT_OP_TYPE* VECT_OP_FUNC(CopyN)( VECT_OP_TYPE* bp, unsigned bn, unsigned d_stride, const VECT_OP_TYPE* sp, unsigned s_stride )
  67. {
  68. VECT_OP_TYPE* dbp = bp;
  69. const VECT_OP_TYPE* ep = bp + (bn*d_stride);
  70. for(; bp < ep; bp += d_stride, sp += s_stride )
  71. *bp = *sp;
  72. return dbp;
  73. }
  74. VECT_OP_TYPE* VECT_OP_FUNC(CopyU)( VECT_OP_TYPE* bp, unsigned bn, const unsigned* sp )
  75. {
  76. VECT_OP_TYPE* dbp = bp;
  77. const VECT_OP_TYPE* ep = bp + bn;
  78. VECT_OP_TYPE* dp = bp;
  79. while( dp < ep )
  80. *dp++ = (VECT_OP_TYPE)*sp++;
  81. return dbp;
  82. }
  83. VECT_OP_TYPE* VECT_OP_FUNC(CopyI)( VECT_OP_TYPE* dbp, unsigned dn, const int* sp )
  84. {
  85. const VECT_OP_TYPE* dep = dbp + dn;
  86. VECT_OP_TYPE* dp = dbp;
  87. while( dp < dep )
  88. *dp++ = (VECT_OP_TYPE)*sp++;
  89. return dbp;
  90. }
  91. VECT_OP_TYPE* VECT_OP_FUNC(CopyF)( VECT_OP_TYPE* dbp, unsigned dn, const float* sp )
  92. {
  93. const VECT_OP_TYPE* dep = dbp + dn;
  94. VECT_OP_TYPE* dp = dbp;
  95. while( dp < dep )
  96. *dp++ = (VECT_OP_TYPE)*sp++;
  97. return dbp;
  98. }
  99. VECT_OP_TYPE* VECT_OP_FUNC(CopyD)( VECT_OP_TYPE* dbp, unsigned dn, const double* sp )
  100. {
  101. const VECT_OP_TYPE* dep = dbp + dn;
  102. VECT_OP_TYPE* dp = dbp;
  103. while( dp < dep )
  104. *dp++ = (VECT_OP_TYPE)*sp++;
  105. return dbp;
  106. }
  107. VECT_OP_TYPE* VECT_OP_FUNC(CopyS)( VECT_OP_TYPE* dbp, unsigned dn, const cmSample_t* sp )
  108. {
  109. const VECT_OP_TYPE* dep = dbp + dn;
  110. VECT_OP_TYPE* dp = dbp;
  111. while( dp < dep )
  112. *dp++ = (VECT_OP_TYPE)*sp++;
  113. return dbp;
  114. }
  115. VECT_OP_TYPE* VECT_OP_FUNC(CopyR)( VECT_OP_TYPE* dbp, unsigned dn, const cmReal_t* sp )
  116. {
  117. const VECT_OP_TYPE* dep = dbp + dn;
  118. VECT_OP_TYPE* dp = dbp;
  119. while( dp < dep )
  120. *dp++ = (VECT_OP_TYPE)*sp++;
  121. return dbp;
  122. }
  123. VECT_OP_TYPE* VECT_OP_FUNC(CopyStride)( VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE* sp, unsigned srcStride )
  124. {
  125. const VECT_OP_TYPE* dep = dbp + dn;
  126. VECT_OP_TYPE* dp = dbp;
  127. for(; dp < dep; sp += srcStride )
  128. *dp++ = *sp;
  129. return dbp;
  130. }
  131. VECT_OP_TYPE* VECT_OP_FUNC(Shrink)( VECT_OP_TYPE* s, unsigned sn, const VECT_OP_TYPE* t, unsigned tn )
  132. {
  133. assert( s <= t && t <= (s+sn) );
  134. assert( s <= (t+tn) && (t+tn) <= (s+sn));
  135. //VECT_OP_FUNC(Move)(s,sn - ((t - s) + tn),t+tn);
  136. VECT_OP_FUNC(Move)((VECT_OP_TYPE*)t,(sn - ((t+tn)-s)) + 1,t+tn);
  137. return s;
  138. }
  139. VECT_OP_TYPE* VECT_OP_FUNC(Expand)( VECT_OP_TYPE* s, unsigned sn, const VECT_OP_TYPE* t, unsigned tn )
  140. {
  141. assert( s <= t && t <= s+sn );
  142. unsigned i = t - s;
  143. s = cmMemResizeP(VECT_OP_TYPE,s,sn+tn);
  144. t = s + i;
  145. assert( t + tn + sn - i == s + sn + tn );
  146. VECT_OP_FUNC(Move)(((VECT_OP_TYPE*)t)+tn,sn-i,t);
  147. return s;
  148. }
  149. VECT_OP_TYPE* VECT_OP_FUNC(Replace)(VECT_OP_TYPE* s, unsigned* sn, const VECT_OP_TYPE* t, unsigned tn, const VECT_OP_TYPE* u, unsigned un )
  150. {
  151. // if s is empty and t[tn] is empty
  152. if( s == NULL && tn == 0 )
  153. {
  154. if( un == 0 )
  155. return s;
  156. s = cmMemAllocZ(VECT_OP_TYPE,un);
  157. VECT_OP_FUNC(Copy)(s,un,u);
  158. if( sn != NULL )
  159. *sn = un;
  160. return s;
  161. }
  162. assert( s!=NULL && t != NULL );
  163. assert( (u!=NULL && un>0) || (u==NULL && un==0) );
  164. if( (tn==0 && un==0) || (t==NULL && u==NULL))
  165. return s;
  166. // if the area to replace is greater than the area to insert ...
  167. if( tn > un )
  168. {
  169. VECT_OP_FUNC(Shrink)(s,*sn,t+un,tn-un); // ... then shrink the buffer
  170. *sn -= tn-un;
  171. }
  172. else
  173. // if the area to insert is greater than the area to replace ...
  174. if( un > tn )
  175. {
  176. unsigned offs = t - s;
  177. s = VECT_OP_FUNC(Expand)(s,*sn,t+tn,un-tn); // ... then expand the buffer
  178. t = s + offs;
  179. *sn += un-tn;
  180. }
  181. assert(t+un <= s+(*sn));
  182. if( u!=NULL )
  183. VECT_OP_FUNC(Copy)((VECT_OP_TYPE*)t,un,u);
  184. return s;
  185. }
  186. VECT_OP_TYPE* VECT_OP_FUNC(Rotate)( VECT_OP_TYPE* dbp, unsigned dn, int shiftCnt )
  187. {
  188. VECT_OP_TYPE* dep = dbp + dn;
  189. int i = 0;
  190. unsigned k = 0;
  191. int n = dep - dbp;
  192. VECT_OP_TYPE t1 = dbp[i];
  193. for(k=0; k<n; ++k)
  194. {
  195. int j;
  196. j = (i + shiftCnt) % n;
  197. if( j<0 )
  198. j += n;
  199. VECT_OP_TYPE t2 = dbp[j];
  200. dbp[j] = t1;
  201. t1 = t2;
  202. i = j;
  203. }
  204. return dbp;
  205. }
  206. VECT_OP_TYPE* VECT_OP_FUNC(RotateM)( VECT_OP_TYPE* dbp, unsigned drn, unsigned dcn, const VECT_OP_TYPE* sbp, int rShiftCnt, int cShiftCnt )
  207. {
  208. int j;
  209. while( rShiftCnt < 0 )
  210. rShiftCnt += drn;
  211. while( cShiftCnt < 0 )
  212. cShiftCnt += dcn;
  213. int m = rShiftCnt % drn;
  214. int n = cShiftCnt % dcn;
  215. for(j=0; j<dcn; ++j,++n)
  216. {
  217. if(n==dcn)
  218. n = 0;
  219. // cnt from dst position to end of column
  220. unsigned cn = drn - m;
  221. // copy from top of src col to bottom of dst column
  222. VECT_OP_FUNC(Copy)(dbp + (n*drn) + m, cn, sbp );
  223. sbp+=cn;
  224. if( cn < drn )
  225. {
  226. // copy from bottom of src col to top of dst column
  227. VECT_OP_FUNC(Copy)(dbp + (n*drn), drn-cn, sbp );
  228. sbp += drn-cn;
  229. }
  230. }
  231. return dbp;
  232. }
  233. VECT_OP_TYPE* VECT_OP_FUNC(Shift)( VECT_OP_TYPE* dbp, unsigned dn, int shiftCnt, VECT_OP_TYPE fillValue )
  234. {
  235. VECT_OP_TYPE* dep = dbp + dn;
  236. VECT_OP_TYPE* rp = dbp;
  237. unsigned n = dep - dbp;
  238. if( shiftCnt == 0 )
  239. return dbp;
  240. if( abs(shiftCnt) >= n )
  241. return VECT_OP_FUNC(Fill)(dbp,dn,fillValue);
  242. if( shiftCnt > 0 )
  243. {
  244. const VECT_OP_TYPE* sbp = dep - (shiftCnt+1);
  245. const VECT_OP_TYPE* sep = dbp;
  246. VECT_OP_TYPE* dp = dbp + (n-1);
  247. while( sbp >= sep )
  248. *dp-- = *sbp--;
  249. while(dbp <= dp )
  250. *dbp++ = fillValue;
  251. }
  252. else
  253. {
  254. const VECT_OP_TYPE* sbp = dbp + abs(shiftCnt);
  255. while( sbp < dep )
  256. *dbp++ = *sbp++;
  257. while(dbp<dep)
  258. *dbp++ = fillValue;
  259. }
  260. return rp;
  261. }
  262. VECT_OP_TYPE* VECT_OP_FUNC(Flip)( VECT_OP_TYPE* dbp, unsigned dn)
  263. {
  264. VECT_OP_TYPE* p0 = dbp;
  265. VECT_OP_TYPE* p1 = dbp + dn - 1;
  266. while( p0 < p1 )
  267. {
  268. VECT_OP_TYPE t = *p0;
  269. *p0++ = *p1;
  270. *p1-- = t;
  271. }
  272. return dbp;
  273. }
  274. VECT_OP_TYPE* VECT_OP_FUNC(SubVS)( VECT_OP_TYPE* bp, unsigned n, VECT_OP_TYPE v )
  275. {
  276. const VECT_OP_TYPE* ep = bp + n;
  277. VECT_OP_TYPE* dp = bp;
  278. while( dp < ep )
  279. *dp++ -= v;
  280. return bp;
  281. }
  282. VECT_OP_TYPE* VECT_OP_FUNC(SubVV)( VECT_OP_TYPE* bp, unsigned n, const VECT_OP_TYPE* v )
  283. {
  284. const VECT_OP_TYPE* ep = bp + n;
  285. VECT_OP_TYPE* dp = bp;
  286. while( dp < ep )
  287. *dp++ -= *v++;
  288. return bp;
  289. }
  290. VECT_OP_TYPE* VECT_OP_FUNC(SubVVS)( VECT_OP_TYPE* bp, unsigned n, const VECT_OP_TYPE* v, VECT_OP_TYPE s )
  291. {
  292. const VECT_OP_TYPE* ep = bp + n;
  293. VECT_OP_TYPE* dp = bp;
  294. while( dp < ep )
  295. *dp++ = *v++ - s;
  296. return bp;
  297. }
  298. VECT_OP_TYPE* VECT_OP_FUNC(SubVVNN)(VECT_OP_TYPE* dp, unsigned dn, unsigned dnn, const VECT_OP_TYPE* v, unsigned n )
  299. {
  300. const VECT_OP_TYPE* ep = dp + (dn*dnn);
  301. VECT_OP_TYPE* dbp = dp;
  302. for(; dp < ep; dp+=dnn, v+=n )
  303. *dp -= *v;
  304. return dbp;
  305. }
  306. VECT_OP_TYPE* VECT_OP_FUNC(SubVVV)( VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE* sb0p, const VECT_OP_TYPE* sb1p )
  307. {
  308. const VECT_OP_TYPE* dep = dbp + dn;
  309. VECT_OP_TYPE* dp = dbp;
  310. while( dbp < dep )
  311. *dbp++ = *sb0p++ - *sb1p++;
  312. return dp;
  313. }
  314. VECT_OP_TYPE* VECT_OP_FUNC(SubVSV)( VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE s0, const VECT_OP_TYPE* sb1p )
  315. {
  316. const VECT_OP_TYPE* dep = dbp + dn;
  317. VECT_OP_TYPE* dp = dbp;
  318. while( dbp < dep )
  319. *dbp++ = s0 - *sb1p++;
  320. return dp;
  321. }
  322. VECT_OP_TYPE* VECT_OP_FUNC(AddVS)( VECT_OP_TYPE* bp, unsigned n, VECT_OP_TYPE v )
  323. {
  324. const VECT_OP_TYPE* ep = bp + n;
  325. VECT_OP_TYPE* dp = bp;
  326. while( dp < ep )
  327. *dp++ += v;
  328. return bp;
  329. }
  330. VECT_OP_TYPE* VECT_OP_FUNC(AddVV)( VECT_OP_TYPE* bp, unsigned bn, const VECT_OP_TYPE* v )
  331. {
  332. const VECT_OP_TYPE* ep = bp + bn;
  333. VECT_OP_TYPE* dp = bp;
  334. while( dp < ep )
  335. *dp++ += *v++;
  336. return bp;
  337. }
  338. VECT_OP_TYPE* VECT_OP_FUNC(AddVVS)( VECT_OP_TYPE* bp, unsigned bn, const VECT_OP_TYPE* v, VECT_OP_TYPE s )
  339. {
  340. const VECT_OP_TYPE* ep = bp + bn;
  341. VECT_OP_TYPE* dp = bp;
  342. while( dp < ep )
  343. *dp++ = *v++ + s;
  344. return bp;
  345. }
  346. VECT_OP_TYPE* VECT_OP_FUNC(AddVVNN)(VECT_OP_TYPE* dp, unsigned dn, unsigned dnn, const VECT_OP_TYPE* v, unsigned n )
  347. {
  348. const VECT_OP_TYPE* ep = dp + (dn*dnn);
  349. VECT_OP_TYPE* dbp = dp;
  350. for(; dp < ep; v+=n, dp+=dnn )
  351. *dp += *v;
  352. return dbp;
  353. }
  354. VECT_OP_TYPE* VECT_OP_FUNC(AddVVV)( VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE* sb0p, const VECT_OP_TYPE* sb1p )
  355. {
  356. const VECT_OP_TYPE* dep = dbp + dn;
  357. VECT_OP_TYPE* dp = dbp;
  358. while( dbp < dep )
  359. *dbp++ = *sb0p++ + *sb1p++;
  360. return dp;
  361. }
  362. VECT_OP_TYPE* VECT_OP_FUNC(MultVVV)( VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE* sb0p, const VECT_OP_TYPE* sb1p )
  363. {
  364. const VECT_OP_TYPE* dep = dbp + dn;
  365. VECT_OP_TYPE* dp = dbp;
  366. while( dbp < dep )
  367. *dbp++ = *sb0p++ * *sb1p++;
  368. return dp;
  369. }
  370. VECT_OP_TYPE* VECT_OP_FUNC(MultVV)( VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE* sbp )
  371. {
  372. const VECT_OP_TYPE* dep = dbp + dn;
  373. VECT_OP_TYPE* dp = dbp;
  374. while( dbp < dep )
  375. *dbp++ *= *sbp++;
  376. return dp;
  377. }
  378. VECT_OP_TYPE* VECT_OP_FUNC(MultVVNN)(VECT_OP_TYPE* dp, unsigned dn, unsigned dnn, const VECT_OP_TYPE* v, unsigned n )
  379. {
  380. const VECT_OP_TYPE* ep = dp + (dn*dnn);
  381. VECT_OP_TYPE* dbp = dp;
  382. for(; dp < ep; v+=n, dp+=dnn )
  383. *dp *= *v;
  384. return dbp;
  385. }
  386. VECT_OP_TYPE* VECT_OP_FUNC(MultVS)( VECT_OP_TYPE* dbp, unsigned dn, VECT_OP_TYPE s )
  387. {
  388. const VECT_OP_TYPE* dep = dbp + dn;
  389. VECT_OP_TYPE* dp = dbp;
  390. while( dbp < dep )
  391. *dbp++ *= s;
  392. return dp;
  393. }
  394. VECT_OP_TYPE* VECT_OP_FUNC(MultVVS)( VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE* sbp, VECT_OP_TYPE s )
  395. {
  396. const VECT_OP_TYPE* dep = dbp + dn;
  397. VECT_OP_TYPE* dp = dbp;
  398. while( dbp < dep )
  399. *dbp++ = *sbp++ * s;
  400. return dp;
  401. }
  402. VECT_OP_TYPE* VECT_OP_FUNC(MultVaVS)( VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE* sbp, VECT_OP_TYPE s )
  403. {
  404. const VECT_OP_TYPE* dep = dbp + dn;
  405. VECT_OP_TYPE* dp = dbp;
  406. while( dbp < dep )
  407. *dbp++ += *sbp++ * s;
  408. return dp;
  409. }
  410. VECT_OP_TYPE* VECT_OP_FUNC(MultSumVVS)(VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE* sbp, VECT_OP_TYPE s )
  411. {
  412. const VECT_OP_TYPE* dep = dbp + dn;
  413. VECT_OP_TYPE* dp = dbp;
  414. while( dbp < dep )
  415. *dbp++ += *sbp++ * s;
  416. return dp;
  417. }
  418. VECT_OP_TYPE* VECT_OP_FUNC(DivVVS)( VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE* sb0p, VECT_OP_TYPE s1 )
  419. {
  420. const VECT_OP_TYPE* dep = dbp + dn;
  421. VECT_OP_TYPE* dp = dbp;
  422. while( dbp < dep )
  423. *dbp++ = *sb0p++ / s1;
  424. return dp;
  425. }
  426. VECT_OP_TYPE* VECT_OP_FUNC(DivVV)( VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE* sb0p )
  427. {
  428. const VECT_OP_TYPE* dep = dbp + dn;
  429. VECT_OP_TYPE* dp = dbp;
  430. while( dbp < dep )
  431. *dbp++ /= *sb0p++;
  432. return dp;
  433. }
  434. VECT_OP_TYPE* VECT_OP_FUNC(DivVVV)( VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE* sb0p, const VECT_OP_TYPE* sb1p )
  435. {
  436. const VECT_OP_TYPE* dep = dbp + dn;
  437. VECT_OP_TYPE* dp = dbp;
  438. while( dbp < dep )
  439. *dbp++ = *sb0p++ / *sb1p++;
  440. return dp;
  441. }
  442. VECT_OP_TYPE* VECT_OP_FUNC(DivVVNN)(VECT_OP_TYPE* dp, unsigned dn, unsigned dnn, const VECT_OP_TYPE* v, unsigned n )
  443. {
  444. const VECT_OP_TYPE* ep = dp + (dn*dnn);
  445. VECT_OP_TYPE* dbp = dp;
  446. for(; dp < ep; v+=n, dp+=dnn )
  447. *dp /= *v;
  448. return dbp;
  449. }
  450. VECT_OP_TYPE* VECT_OP_FUNC(DivVS)( VECT_OP_TYPE* dbp, unsigned dn, VECT_OP_TYPE s )
  451. {
  452. const VECT_OP_TYPE* dep = dbp + dn;
  453. VECT_OP_TYPE* dp = dbp;
  454. while( dbp < dep )
  455. *dbp++ /= s;
  456. return dp;
  457. }
  458. VECT_OP_TYPE* VECT_OP_FUNC(DivVSV)( VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE s0, const VECT_OP_TYPE* sb1p )
  459. {
  460. const VECT_OP_TYPE* dep = dbp + dn;
  461. VECT_OP_TYPE* dp = dbp;
  462. while( dbp < dep )
  463. *dbp++ = s0 / *sb1p++;
  464. return dp;
  465. }
  466. VECT_OP_TYPE* VECT_OP_FUNC(DivVVZ)( VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE* sb0p )
  467. {
  468. const VECT_OP_TYPE* dep = dbp + dn;
  469. VECT_OP_TYPE* dp = dbp;
  470. for(; dbp < dep; ++sb0p )
  471. if( *sb0p == 0 )
  472. *dbp++ = 0;
  473. else
  474. *dbp++ /= *sb0p;
  475. return dp;
  476. }
  477. VECT_OP_TYPE* VECT_OP_FUNC(DivVVVZ)( VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE* sb0p, const VECT_OP_TYPE* sb1p )
  478. {
  479. const VECT_OP_TYPE* dep = dbp + dn;
  480. VECT_OP_TYPE* dp = dbp;
  481. for(; dbp < dep; ++sb0p,++sb1p )
  482. if( *sb1p == 0 )
  483. *dbp++ = 0;
  484. else
  485. *dbp++ = *sb0p / *sb1p;
  486. return dp;
  487. }
  488. VECT_OP_TYPE* VECT_OP_FUNC(DivMS)( VECT_OP_TYPE* dp, unsigned drn, unsigned dcn, const VECT_OP_TYPE* sp )
  489. {
  490. unsigned i;
  491. for(i=0; i<dcn; ++i)
  492. VECT_OP_FUNC(DivVS)( dp + i*drn, drn, sp[i] );
  493. return dp;
  494. }
  495. VECT_OP_TYPE VECT_OP_FUNC(Sum)( const VECT_OP_TYPE* bp, unsigned n )
  496. {
  497. const VECT_OP_TYPE* ep = bp + n;
  498. VECT_OP_TYPE s = 0;
  499. while( bp < ep )
  500. s += *bp++;
  501. return s;
  502. }
  503. VECT_OP_TYPE VECT_OP_FUNC(SumN)( const VECT_OP_TYPE* bp, unsigned n, unsigned stride )
  504. {
  505. const VECT_OP_TYPE* ep = bp + (n*stride);
  506. VECT_OP_TYPE s = 0;
  507. for(; bp < ep; bp += stride )
  508. s += *bp;
  509. return s;
  510. }
  511. VECT_OP_TYPE* VECT_OP_FUNC(SumM)(const VECT_OP_TYPE* sp, unsigned srn, unsigned scn, VECT_OP_TYPE* dp )
  512. {
  513. unsigned i;
  514. for(i=0; i<scn; ++i)
  515. dp[i] = VECT_OP_FUNC(Sum)(sp + (i*srn), srn );
  516. return dp;
  517. }
  518. VECT_OP_TYPE* VECT_OP_FUNC(SumMN)(const VECT_OP_TYPE* sp, unsigned srn, unsigned scn, VECT_OP_TYPE* dp )
  519. {
  520. unsigned i;
  521. for(i=0; i<srn; ++i)
  522. dp[i] = VECT_OP_FUNC(SumN)(sp + i, scn, srn );
  523. return dp;
  524. }
  525. VECT_OP_TYPE* VECT_OP_FUNC(Abs)( VECT_OP_TYPE* dbp, unsigned dn )
  526. {
  527. unsigned i;
  528. for(i=0; i<dn; ++i)
  529. if( dbp[i]<0 )
  530. dbp[i] = -dbp[i];
  531. return dbp;
  532. }
  533. // mi is a target value - it holds the number of elements in ap[an] which must be be less than the median value.
  534. // If the initial array contains an even number of values then the median value is formed by averaging the two center values.
  535. // In this case *evenFlPtr is set and used to indicate that the center-upper value must be found during undwinding.
  536. VECT_OP_TYPE VECT_OP_FUNC(MedianSearch)( unsigned mi, const VECT_OP_TYPE* ap, unsigned an, bool* evenFlPtr )
  537. {
  538. VECT_OP_TYPE x = ap[0]; // pick a random value as a potential median value
  539. VECT_OP_TYPE a1[ an ]; // values below x
  540. VECT_OP_TYPE a3[ an ]; // values above x
  541. unsigned a1n = 0;
  542. unsigned a2n = 0; // values equal to x
  543. unsigned a3n = 0;
  544. const VECT_OP_TYPE* abp = ap;
  545. const VECT_OP_TYPE* aep = abp + an;
  546. for(; abp < aep; ++abp )
  547. {
  548. if( *abp < x )
  549. a1[a1n++] = *abp;
  550. else
  551. {
  552. if( *abp > x )
  553. a3[a3n++] = *abp;
  554. else
  555. ++a2n;
  556. }
  557. }
  558. //printf("%i : %i %i %i\n",mi,a1n,a2n,a3n);
  559. // there are more values below x (mi remains the target split point)
  560. if( a1n > mi )
  561. {
  562. x = VECT_OP_FUNC(MedianSearch)(mi,a1,a1n,evenFlPtr);
  563. }
  564. else
  565. {
  566. // the target was located
  567. if( a1n+a2n >= mi )
  568. {
  569. // if a1n alone matches mi then the max value in a1[] holds the median value otherwise x is the median
  570. if(a1n>=1 && a1n==mi)
  571. {
  572. VECT_OP_TYPE mv = VECT_OP_FUNC(Max)(a1,a1n,1);
  573. x = *evenFlPtr ? (mv+x)/2 : mv;
  574. *evenFlPtr = false;
  575. }
  576. // if the evenFl is set then the closest value above the median (x) must be located
  577. if( *evenFlPtr )
  578. {
  579. // if the next greater value is in a2[]
  580. if( a2n > 1 && (a1n+a2n) > mi )
  581. *evenFlPtr = false;
  582. else
  583. // if the next greater value is in a3[]
  584. if( a3n > 1 )
  585. {
  586. x = (x + VECT_OP_FUNC(Min)(a3,a3n,1))/2;
  587. *evenFlPtr = false;
  588. }
  589. }
  590. // no need for unwind processing - all the possibilities at this level have been exhausted
  591. return x;
  592. }
  593. else
  594. {
  595. // There are more values above x - the median must therefore be in a3[].
  596. // Reset mi cmcounting for the fact that we know that there are
  597. // a1n+a2n values below the lowest value in a3.
  598. x = VECT_OP_FUNC(MedianSearch)(mi - (a1n+a2n), a3, a3n, evenFlPtr );
  599. }
  600. }
  601. if( *evenFlPtr )
  602. {
  603. // find the first value greater than x
  604. while( ap < aep && *ap <= x )
  605. ++ap;
  606. if( ap < aep )
  607. {
  608. VECT_OP_TYPE v = *ap++;
  609. // find the nearest value greater than x
  610. for(; ap < aep; ++ap )
  611. if( *ap > x && ((*ap - x) < (v-x)))
  612. v = *ap;
  613. x = (v + x)/2;
  614. *evenFlPtr = false;
  615. }
  616. }
  617. return x;
  618. }
  619. VECT_OP_TYPE VECT_OP_FUNC(Median)( const VECT_OP_TYPE* bp, unsigned n )
  620. {
  621. bool evenFl = cmIsEvenU(n);
  622. unsigned medIdx = evenFl ? n/2 : (n+1)/2;
  623. return VECT_OP_FUNC(MedianSearch)( medIdx, bp, n, &evenFl );
  624. }
  625. unsigned VECT_OP_FUNC(MinIndex)( const VECT_OP_TYPE* bp, unsigned n, unsigned stride )
  626. {
  627. const VECT_OP_TYPE* ep = bp + (n*stride);
  628. if( bp >= ep )
  629. return cmInvalidIdx;
  630. const VECT_OP_TYPE* p = bp;
  631. const VECT_OP_TYPE* mp = bp;
  632. bp+=stride;
  633. for(; bp < ep; bp+=stride )
  634. if( *bp < *mp )
  635. mp = bp;
  636. return (mp - p)/stride;
  637. }
  638. unsigned VECT_OP_FUNC(MaxIndex)( const VECT_OP_TYPE* bp, unsigned n, unsigned stride )
  639. {
  640. const VECT_OP_TYPE* ep = bp + (n*stride);
  641. if( bp >= ep )
  642. return cmInvalidIdx;
  643. const VECT_OP_TYPE* p = bp;
  644. const VECT_OP_TYPE* mp = bp;
  645. bp+=stride;
  646. for(; bp < ep; bp+=stride )
  647. if( *bp > *mp )
  648. mp = bp;
  649. return (mp - p)/stride;
  650. }
  651. VECT_OP_TYPE VECT_OP_FUNC(Min)( const VECT_OP_TYPE* bp, unsigned n, unsigned stride )
  652. {
  653. unsigned i;
  654. if((i = VECT_OP_FUNC(MinIndex)(bp,n,stride)) == cmInvalidIdx )
  655. {
  656. assert(0);
  657. return 0;
  658. }
  659. return bp[i*stride];
  660. }
  661. VECT_OP_TYPE VECT_OP_FUNC(Max)( const VECT_OP_TYPE* bp, unsigned n, unsigned stride )
  662. {
  663. unsigned i;
  664. if((i = VECT_OP_FUNC(MaxIndex)(bp,n,stride)) == cmInvalidIdx )
  665. {
  666. assert(0);
  667. return 0;
  668. }
  669. return bp[i*stride];
  670. }
  671. VECT_OP_TYPE* VECT_OP_FUNC(MinVV)( VECT_OP_TYPE* dp, unsigned dn, const VECT_OP_TYPE* sp )
  672. {
  673. unsigned i;
  674. for(i=0; i<dn; ++i)
  675. if( sp[i] < dp[i] )
  676. dp[i] = sp[i];
  677. return dp;
  678. }
  679. VECT_OP_TYPE* VECT_OP_FUNC(MaxVV)( VECT_OP_TYPE* dp, unsigned dn, const VECT_OP_TYPE* sp )
  680. {
  681. unsigned i;
  682. for(i=0; i<dn; ++i)
  683. if( sp[i] > dp[i] )
  684. dp[i] = sp[i];
  685. return dp;
  686. }
  687. unsigned* VECT_OP_FUNC(MinIndexM)( unsigned* dp, const VECT_OP_TYPE* sp, unsigned srn, unsigned scn )
  688. {
  689. unsigned i = 0;
  690. for(i=0; i<scn; ++i)
  691. dp[i] = VECT_OP_FUNC(MinIndex)(sp + (i*srn), srn, 1 );
  692. return dp;
  693. }
  694. unsigned* VECT_OP_FUNC(MaxIndexM)( unsigned* dp, const VECT_OP_TYPE* sp, unsigned srn, unsigned scn )
  695. {
  696. unsigned i = 0;
  697. for(i=0; i<scn; ++i)
  698. dp[i] = VECT_OP_FUNC(MaxIndex)(sp + (i*srn), srn, 1 );
  699. return dp;
  700. }
  701. VECT_OP_TYPE VECT_OP_FUNC(Mode)( const VECT_OP_TYPE* sp, unsigned sn )
  702. {
  703. unsigned n[sn];
  704. VECT_OP_TYPE v[sn];
  705. unsigned i,j,k = 0;
  706. unsigned n0 = 0; // idx of most freq occurring ele
  707. unsigned n1 = -1; // idx of 2nd most freq occurring ele
  708. for(i=0; i<sn; ++i)
  709. {
  710. // find sp[i] in v[]
  711. for(j=0; j<k; ++j)
  712. if( sp[i] == v[j] )
  713. {
  714. ++n[j];
  715. break;
  716. }
  717. // sp[i] was not found in v[]
  718. if( k == j )
  719. {
  720. v[j] = sp[i];
  721. n[j] = 1;
  722. ++k;
  723. }
  724. // n[j] holds frq of sp[i]
  725. // do nothing if j is already most freq
  726. if( j != n0 )
  727. {
  728. // if j is new most freq
  729. if( n[j] > n[n0] )
  730. {
  731. n1 = n0;
  732. n0 = j;
  733. }
  734. else
  735. // if j is 2nd most freq
  736. if( (n1==-1) || (n[j] > n[n1]) )
  737. n1 = j;
  738. }
  739. // if diff between two most freq is greater than remaining ele's
  740. if( (n1!=-1) && (n[n0]-n[n1]) >= (sn-i) )
  741. break;
  742. }
  743. // if there are no ele's with same count
  744. if( n[n0] > n[n1] )
  745. return v[n0];
  746. // break tie between ele's with same count be returning min value
  747. // (this is the same as Matlab tie break criteria)
  748. j = 0;
  749. for(i=1; i<k; ++i)
  750. if( (n[i] > n[j]) || (n[i] == n[j] && v[i] < v[j]) )
  751. j=i;
  752. return v[j];
  753. }
  754. unsigned VECT_OP_FUNC(Find)( const VECT_OP_TYPE* sp, unsigned sn, VECT_OP_TYPE key )
  755. {
  756. const VECT_OP_TYPE* sbp = sp;
  757. const VECT_OP_TYPE* ep = sp + sn;
  758. while( sp<ep )
  759. if( *sp++ == key )
  760. break;
  761. if( sp==ep )
  762. return cmInvalidIdx;
  763. return (sp-1) - sbp;
  764. }
  765. unsigned VECT_OP_FUNC(Count)( const VECT_OP_TYPE* sp, unsigned sn, VECT_OP_TYPE key )
  766. {
  767. unsigned cnt = 0;
  768. const VECT_OP_TYPE* ep = sp + sn;
  769. while( sp<ep )
  770. if( *sp++ == key )
  771. ++cnt;
  772. return cnt;
  773. }
  774. VECT_OP_TYPE* VECT_OP_FUNC(ReplaceLte)( VECT_OP_TYPE* dp, unsigned dn, const VECT_OP_TYPE* sp, VECT_OP_TYPE lteKeyVal, VECT_OP_TYPE replaceVal )
  775. {
  776. VECT_OP_TYPE* rp = dp;
  777. const VECT_OP_TYPE* ep = dp + dn;
  778. for(; dp < ep; ++sp )
  779. *dp++ = *sp <= lteKeyVal ? replaceVal : *sp;
  780. return rp;
  781. }
  782. VECT_OP_TYPE* VECT_OP_FUNC(Diag)( VECT_OP_TYPE* dbp, unsigned n, const VECT_OP_TYPE* sbp )
  783. {
  784. unsigned i,j;
  785. for(i=0,j=0; i<n && j<n; ++i,++j)
  786. dbp[ (i*n) + j ] = sbp[i];
  787. return dbp;
  788. }
  789. VECT_OP_TYPE* VECT_OP_FUNC(DiagZ)(VECT_OP_TYPE* dbp, unsigned n, const VECT_OP_TYPE* sbp )
  790. {
  791. VECT_OP_FUNC(Fill)(dbp,n*n,0);
  792. return VECT_OP_FUNC(Diag)(dbp,n,sbp);
  793. }
  794. VECT_OP_TYPE* VECT_OP_FUNC(Identity)( VECT_OP_TYPE* dbp, unsigned rn, unsigned cn )
  795. {
  796. unsigned i,j;
  797. for(i=0,j=0; i<cn && j<rn; ++i,++j)
  798. dbp[ (i*rn) + j ] = 1;
  799. return dbp;
  800. }
  801. VECT_OP_TYPE* VECT_OP_FUNC(IdentityZ)( VECT_OP_TYPE* dbp, unsigned rn, unsigned cn )
  802. {
  803. VECT_OP_FUNC(Fill)(dbp,rn*cn,0);
  804. return VECT_OP_FUNC(Identity)(dbp,rn,cn);
  805. }
  806. VECT_OP_TYPE* VECT_OP_FUNC(Transpose)( VECT_OP_TYPE* dbp, const VECT_OP_TYPE* sp, unsigned srn, unsigned scn )
  807. {
  808. VECT_OP_TYPE* dp = dbp;
  809. const VECT_OP_TYPE* dep = dbp + (srn*scn);
  810. while( dbp < dep )
  811. {
  812. const VECT_OP_TYPE* sbp = sp++;
  813. const VECT_OP_TYPE* sep = sbp + (srn*scn);
  814. for(; sbp < sep; sbp+=srn )
  815. *dbp++ = *sbp;
  816. }
  817. return dp;
  818. }
  819. VECT_OP_TYPE VECT_OP_FUNC(Seq)( VECT_OP_TYPE* dbp, unsigned dn, VECT_OP_TYPE beg, VECT_OP_TYPE incr )
  820. {
  821. const VECT_OP_TYPE* dep = dbp + dn;
  822. unsigned i = 0;
  823. for(; dbp<dep; ++i)
  824. *dbp++ = beg + (incr*i);
  825. return beg + (incr*i);
  826. }
  827. void VECT_OP_FUNC(FnThresh)( const VECT_OP_TYPE* xV, unsigned xN, unsigned wndN, VECT_OP_TYPE* yV, unsigned yStride, VECT_OP_TYPE (*fnPtr)(const VECT_OP_TYPE*, unsigned) )
  828. {
  829. int i0 = cmIsOddU(wndN) ? (wndN-1)/2 : wndN/2;
  830. int i1 = cmIsOddU(wndN) ? (wndN-1)/2 : wndN/2 - 1;
  831. int i,j;
  832. i0 = -i0;
  833. if( fnPtr == NULL )
  834. fnPtr = &(VECT_OP_FUNC(Median));
  835. for(i=0; i<xN; ++i,++i0,++i1)
  836. {
  837. j = (i*yStride);
  838. if( i0 < 0 )
  839. if( i1 >= xN )
  840. yV[j] = (*fnPtr)(xV,xN);
  841. else
  842. yV[j] = (*fnPtr)(xV,i1+1);
  843. else if( i1 >= xN )
  844. yV[j] = (*fnPtr)(xV+i0,xN-i0);
  845. else
  846. yV[j] = (*fnPtr)(xV+i0,wndN);
  847. }
  848. }
  849. void VECT_OP_FUNC(MedianFilt)( const VECT_OP_TYPE* xV, unsigned xN, unsigned wndN, VECT_OP_TYPE* yV, unsigned yStride )
  850. {
  851. int i0 = cmIsOddU(wndN) ? (wndN-1)/2 : wndN/2;
  852. int i1 = cmIsOddU(wndN) ? (wndN-1)/2 : wndN/2 - 1;
  853. int i,j;
  854. VECT_OP_TYPE tV[ wndN ];
  855. i0 = -i0;
  856. VECT_OP_FUNC(Fill)(tV,wndN,0);
  857. for(i=0; i<xN; ++i,++i0,++i1)
  858. {
  859. j = (i*yStride);
  860. // note that the position of the zero padding in tV[]
  861. // does not matter because the median calcluation does
  862. // not make any assumptions about the order of the argument
  863. // vector.
  864. if( i0 < 0 )
  865. {
  866. VECT_OP_FUNC(Copy)(tV,wndN+i0,xV);
  867. VECT_OP_FUNC(Fill)(tV+wndN+i0,labs(i0),0);
  868. //VECT_OP_FUNC(Print)(NULL,1,wndN,tV,-1,-1);
  869. yV[j] = VECT_OP_FUNC(Median)(tV,wndN);
  870. continue;
  871. }
  872. if( i1 >= xN )
  873. {
  874. VECT_OP_FUNC(Copy)(tV,wndN-(i1-xN+1),xV+i0);
  875. VECT_OP_FUNC(Fill)(tV+wndN-(i1-xN+1),i1-xN+1,0);
  876. //VECT_OP_FUNC(Print)(NULL,1,wndN,tV,-1,-1);
  877. yV[j] = VECT_OP_FUNC(Median)(tV,wndN);
  878. continue;
  879. }
  880. //VECT_OP_FUNC(Print)(NULL,1,wndN,xV+i0,-1,-1);
  881. yV[j] = VECT_OP_FUNC(Median)(xV+i0,wndN);
  882. }
  883. }
  884. unsigned* VECT_OP_FUNC(LevEditDistAllocMtx)(unsigned maxN)
  885. {
  886. maxN += 1;
  887. unsigned* m = cmMemAllocZ(unsigned,maxN*maxN);
  888. unsigned* p = m;
  889. unsigned i;
  890. // initialize the comparison matrix with the default costs in the
  891. // first row and column
  892. // (Note that this matrix is not oriented in column major order like most 'cm' matrices.)
  893. for(i=0; i<maxN; ++i)
  894. {
  895. p[i] = i; // 0th row
  896. p[ i * maxN ] = i; // 0th col
  897. }
  898. return m;
  899. }
  900. double VECT_OP_FUNC(LevEditDist)(unsigned mtxMaxN, unsigned* m, const VECT_OP_TYPE* s0, int n0, const VECT_OP_TYPE* s1, int n1, unsigned maxN )
  901. {
  902. mtxMaxN += 1;
  903. assert( n0 < mtxMaxN && n1 < mtxMaxN );
  904. int v = 0;
  905. unsigned i;
  906. // Note that m[maxN,maxN] is not oriented in column major order like most 'cm' matrices.
  907. for(i=1; i<n0+1; ++i)
  908. {
  909. unsigned ii = i * mtxMaxN; // current row
  910. unsigned i_1 = ii - mtxMaxN; // previous row
  911. unsigned j;
  912. for( j=1; j<n1+1; ++j)
  913. {
  914. int cost = s0[i-1] == s1[j-1] ? 0 : 1;
  915. //m[i][j] = min( m[i-1][j] + 1, min( m[i][j-1] + 1, m[i-1][j-1] + cost ) );
  916. m[ ii + j ] = v = cmMin( m[ i_1 + j] + 1, cmMin( m[ ii + j - 1] + 1, m[ i_1 + j - 1 ] + cost ) );
  917. }
  918. }
  919. return (double) v / maxN;
  920. }
  921. double VECT_OP_FUNC(LevEditDistWithCostThresh)( int mtxMaxN, unsigned* m, const VECT_OP_TYPE* s0, int n0, const VECT_OP_TYPE* s1, int n1, double maxCost, unsigned maxN )
  922. {
  923. mtxMaxN += 1;
  924. int v = 0;
  925. maxCost = cmMin(1.0,cmMax(0.0,maxCost));
  926. int iMaxCost = ceil( maxCost * maxN );
  927. assert( iMaxCost > 0 && maxCost > 0 );
  928. // If the two strings are different lengths and the min possible distance is
  929. // greater than the threshold then return the threshold as the cost.
  930. // (Note: For strings of different length the min possible distance is the
  931. // difference in length between the two strings).
  932. if( abs(n0-n1) > iMaxCost )
  933. return maxCost;
  934. int i;
  935. // for each row in the matrix ...
  936. for(i=1; i<n0+1; ++i)
  937. {
  938. int ii = i * mtxMaxN; // current row
  939. int i_1 = ii - mtxMaxN; // previous row
  940. // Limit the row to (2*iMaxCost)+1 diagnal strip.
  941. // This strip is based on the idea that the best case can be precomputed for
  942. // all matrix elements in advance - where the best case for position i,j is:
  943. // abs(i-j). This can be justified based on the idea that the least possible
  944. // distance between two strings of length i and j is abs(i-1). The minimum least
  945. // possible distance is therefore found on the matrix diagnal and grows as the
  946. // distance from the diagnal increases.
  947. int ji = cmMax( 1, i - iMaxCost );
  948. int jn = cmMin(iMaxCost + i, n1) + 1;
  949. int j;
  950. // fill in (max cost + 1) as the value in the column before the starting column
  951. // (it will be referred to during the first computation in this row)
  952. if( ji >= 2 )
  953. m[ ii + (ji-1) ] = iMaxCost + 1;
  954. // for each column in the diagnal stripe - beginning with the leftmost column.
  955. for( j=ji; j<jn; ++j)
  956. {
  957. int cost = s0[i-1] == s1[j-1] ? 0 : 1;
  958. m[ ii + j ] = v = cmMin( m[ i_1 + j] + 1, cmMin( m[ ii + j - 1] + 1, m[ i_1 + j - 1 ] + cost ) );
  959. }
  960. // fill in (max cost + 1) in the column following the last column
  961. // (it will be referred to during computation of the following row)
  962. if( j < n1+1 )
  963. m[ii + j] = iMaxCost + 1;
  964. }
  965. assert( v >= 0 );
  966. return cmMin( maxCost , (double) v / maxN);
  967. }
  968. #endif