libcm is a C development framework with an emphasis on audio signal processing applications.
Vous ne pouvez pas sélectionner plus de 25 sujets Les noms de sujets doivent commencer par une lettre ou un nombre, peuvent contenir des tirets ('-') et peuvent comporter jusqu'à 35 caractères.

cmVectOpsRICode.h 29KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240
  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. bool VECT_OP_FUNC(IsEqual)( const VECT_OP_TYPE* s0p, const VECT_OP_TYPE* s1p, unsigned sn )
  702. {
  703. const VECT_OP_TYPE* ep = s0p + sn;
  704. for(; s0p < ep; ++s0p,++s1p )
  705. if( *s0p != *s1p )
  706. return false;
  707. return true;
  708. }
  709. bool VECT_OP_FUNC(IsClose)( const VECT_OP_TYPE* s0p, const VECT_OP_TYPE* s1p, unsigned sn, double pct )
  710. {
  711. const VECT_OP_TYPE* ep = s0p + sn;
  712. for(; s0p < ep; ++s0p,++s1p )
  713. {
  714. double d = abs(*s1p - *s0p);
  715. double s = cmMin(*s1p,*s0p);
  716. if( d*100.0/s > pct )
  717. return false;
  718. }
  719. return true;
  720. }
  721. VECT_OP_TYPE VECT_OP_FUNC(Mode)( const VECT_OP_TYPE* sp, unsigned sn )
  722. {
  723. unsigned n[sn];
  724. VECT_OP_TYPE v[sn];
  725. unsigned i,j,k = 0;
  726. unsigned n0 = 0; // idx of most freq occurring ele
  727. unsigned n1 = -1; // idx of 2nd most freq occurring ele
  728. for(i=0; i<sn; ++i)
  729. {
  730. // find sp[i] in v[]
  731. for(j=0; j<k; ++j)
  732. if( sp[i] == v[j] )
  733. {
  734. ++n[j];
  735. break;
  736. }
  737. // sp[i] was not found in v[]
  738. if( k == j )
  739. {
  740. v[j] = sp[i];
  741. n[j] = 1;
  742. ++k;
  743. }
  744. // n[j] holds frq of sp[i]
  745. // do nothing if j is already most freq
  746. if( j != n0 )
  747. {
  748. // if j is new most freq
  749. if( n[j] > n[n0] )
  750. {
  751. n1 = n0;
  752. n0 = j;
  753. }
  754. else
  755. // if j is 2nd most freq
  756. if( (n1==-1) || (n[j] > n[n1]) )
  757. n1 = j;
  758. }
  759. // if diff between two most freq is greater than remaining ele's
  760. if( (n1!=-1) && (n[n0]-n[n1]) >= (sn-i) )
  761. break;
  762. }
  763. // if there are no ele's with same count
  764. if( n[n0] > n[n1] )
  765. return v[n0];
  766. // break tie between ele's with same count be returning min value
  767. // (this is the same as Matlab tie break criteria)
  768. j = 0;
  769. for(i=1; i<k; ++i)
  770. if( (n[i] > n[j]) || (n[i] == n[j] && v[i] < v[j]) )
  771. j=i;
  772. return v[j];
  773. }
  774. unsigned VECT_OP_FUNC(Find)( const VECT_OP_TYPE* sp, unsigned sn, VECT_OP_TYPE key )
  775. {
  776. const VECT_OP_TYPE* sbp = sp;
  777. const VECT_OP_TYPE* ep = sp + sn;
  778. while( sp<ep )
  779. if( *sp++ == key )
  780. break;
  781. if( sp==ep )
  782. return cmInvalidIdx;
  783. return (sp-1) - sbp;
  784. }
  785. unsigned VECT_OP_FUNC(Count)( const VECT_OP_TYPE* sp, unsigned sn, VECT_OP_TYPE key )
  786. {
  787. unsigned cnt = 0;
  788. const VECT_OP_TYPE* ep = sp + sn;
  789. while( sp<ep )
  790. if( *sp++ == key )
  791. ++cnt;
  792. return cnt;
  793. }
  794. 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 )
  795. {
  796. VECT_OP_TYPE* rp = dp;
  797. const VECT_OP_TYPE* ep = dp + dn;
  798. for(; dp < ep; ++sp )
  799. *dp++ = *sp <= lteKeyVal ? replaceVal : *sp;
  800. return rp;
  801. }
  802. VECT_OP_TYPE* VECT_OP_FUNC(Diag)( VECT_OP_TYPE* dbp, unsigned n, const VECT_OP_TYPE* sbp )
  803. {
  804. unsigned i,j;
  805. for(i=0,j=0; i<n && j<n; ++i,++j)
  806. dbp[ (i*n) + j ] = sbp[i];
  807. return dbp;
  808. }
  809. VECT_OP_TYPE* VECT_OP_FUNC(DiagZ)(VECT_OP_TYPE* dbp, unsigned n, const VECT_OP_TYPE* sbp )
  810. {
  811. VECT_OP_FUNC(Fill)(dbp,n*n,0);
  812. return VECT_OP_FUNC(Diag)(dbp,n,sbp);
  813. }
  814. VECT_OP_TYPE* VECT_OP_FUNC(Identity)( VECT_OP_TYPE* dbp, unsigned rn, unsigned cn )
  815. {
  816. unsigned i,j;
  817. for(i=0,j=0; i<cn && j<rn; ++i,++j)
  818. dbp[ (i*rn) + j ] = 1;
  819. return dbp;
  820. }
  821. VECT_OP_TYPE* VECT_OP_FUNC(IdentityZ)( VECT_OP_TYPE* dbp, unsigned rn, unsigned cn )
  822. {
  823. VECT_OP_FUNC(Fill)(dbp,rn*cn,0);
  824. return VECT_OP_FUNC(Identity)(dbp,rn,cn);
  825. }
  826. VECT_OP_TYPE* VECT_OP_FUNC(Transpose)( VECT_OP_TYPE* dbp, const VECT_OP_TYPE* sp, unsigned srn, unsigned scn )
  827. {
  828. VECT_OP_TYPE* dp = dbp;
  829. const VECT_OP_TYPE* dep = dbp + (srn*scn);
  830. while( dbp < dep )
  831. {
  832. const VECT_OP_TYPE* sbp = sp++;
  833. const VECT_OP_TYPE* sep = sbp + (srn*scn);
  834. for(; sbp < sep; sbp+=srn )
  835. *dbp++ = *sbp;
  836. }
  837. return dp;
  838. }
  839. VECT_OP_TYPE VECT_OP_FUNC(Seq)( VECT_OP_TYPE* dbp, unsigned dn, VECT_OP_TYPE beg, VECT_OP_TYPE incr )
  840. {
  841. const VECT_OP_TYPE* dep = dbp + dn;
  842. unsigned i = 0;
  843. for(; dbp<dep; ++i)
  844. *dbp++ = beg + (incr*i);
  845. return beg + (incr*i);
  846. }
  847. 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) )
  848. {
  849. int i0 = cmIsOddU(wndN) ? (wndN-1)/2 : wndN/2;
  850. int i1 = cmIsOddU(wndN) ? (wndN-1)/2 : wndN/2 - 1;
  851. int i,j;
  852. i0 = -i0;
  853. if( fnPtr == NULL )
  854. fnPtr = &(VECT_OP_FUNC(Median));
  855. for(i=0; i<xN; ++i,++i0,++i1)
  856. {
  857. j = (i*yStride);
  858. if( i0 < 0 )
  859. if( i1 >= xN )
  860. yV[j] = (*fnPtr)(xV,xN);
  861. else
  862. yV[j] = (*fnPtr)(xV,i1+1);
  863. else if( i1 >= xN )
  864. yV[j] = (*fnPtr)(xV+i0,xN-i0);
  865. else
  866. yV[j] = (*fnPtr)(xV+i0,wndN);
  867. }
  868. }
  869. void VECT_OP_FUNC(MedianFilt)( const VECT_OP_TYPE* xV, unsigned xN, unsigned wndN, VECT_OP_TYPE* yV, unsigned yStride )
  870. {
  871. int i0 = cmIsOddU(wndN) ? (wndN-1)/2 : wndN/2;
  872. int i1 = cmIsOddU(wndN) ? (wndN-1)/2 : wndN/2 - 1;
  873. int i,j;
  874. VECT_OP_TYPE tV[ wndN ];
  875. i0 = -i0;
  876. VECT_OP_FUNC(Fill)(tV,wndN,0);
  877. for(i=0; i<xN; ++i,++i0,++i1)
  878. {
  879. j = (i*yStride);
  880. // note that the position of the zero padding in tV[]
  881. // does not matter because the median calcluation does
  882. // not make any assumptions about the order of the argument
  883. // vector.
  884. if( i0 < 0 )
  885. {
  886. VECT_OP_FUNC(Copy)(tV,wndN+i0,xV);
  887. VECT_OP_FUNC(Fill)(tV+wndN+i0,labs(i0),0);
  888. //VECT_OP_FUNC(Print)(NULL,1,wndN,tV,-1,-1);
  889. yV[j] = VECT_OP_FUNC(Median)(tV,wndN);
  890. continue;
  891. }
  892. if( i1 >= xN )
  893. {
  894. VECT_OP_FUNC(Copy)(tV,wndN-(i1-xN+1),xV+i0);
  895. VECT_OP_FUNC(Fill)(tV+wndN-(i1-xN+1),i1-xN+1,0);
  896. //VECT_OP_FUNC(Print)(NULL,1,wndN,tV,-1,-1);
  897. yV[j] = VECT_OP_FUNC(Median)(tV,wndN);
  898. continue;
  899. }
  900. //VECT_OP_FUNC(Print)(NULL,1,wndN,xV+i0,-1,-1);
  901. yV[j] = VECT_OP_FUNC(Median)(xV+i0,wndN);
  902. }
  903. }
  904. unsigned* VECT_OP_FUNC(LevEditDistAllocMtx)(unsigned maxN)
  905. {
  906. maxN += 1;
  907. unsigned* m = cmMemAllocZ(unsigned,maxN*maxN);
  908. unsigned* p = m;
  909. unsigned i;
  910. // initialize the comparison matrix with the default costs in the
  911. // first row and column
  912. // (Note that this matrix is not oriented in column major order like most 'cm' matrices.)
  913. for(i=0; i<maxN; ++i)
  914. {
  915. p[i] = i; // 0th row
  916. p[ i * maxN ] = i; // 0th col
  917. }
  918. return m;
  919. }
  920. 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 )
  921. {
  922. mtxMaxN += 1;
  923. assert( n0 < mtxMaxN && n1 < mtxMaxN );
  924. int v = 0;
  925. unsigned i;
  926. // Note that m[maxN,maxN] is not oriented in column major order like most 'cm' matrices.
  927. for(i=1; i<n0+1; ++i)
  928. {
  929. unsigned ii = i * mtxMaxN; // current row
  930. unsigned i_1 = ii - mtxMaxN; // previous row
  931. unsigned j;
  932. for( j=1; j<n1+1; ++j)
  933. {
  934. int cost = s0[i-1] == s1[j-1] ? 0 : 1;
  935. //m[i][j] = min( m[i-1][j] + 1, min( m[i][j-1] + 1, m[i-1][j-1] + cost ) );
  936. m[ ii + j ] = v = cmMin( m[ i_1 + j] + 1, cmMin( m[ ii + j - 1] + 1, m[ i_1 + j - 1 ] + cost ) );
  937. }
  938. }
  939. return (double) v / maxN;
  940. }
  941. 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 )
  942. {
  943. mtxMaxN += 1;
  944. int v = 0;
  945. maxCost = cmMin(1.0,cmMax(0.0,maxCost));
  946. int iMaxCost = ceil( maxCost * maxN );
  947. assert( iMaxCost > 0 && maxCost > 0 );
  948. // If the two strings are different lengths and the min possible distance is
  949. // greater than the threshold then return the threshold as the cost.
  950. // (Note: For strings of different length the min possible distance is the
  951. // difference in length between the two strings).
  952. if( abs(n0-n1) > iMaxCost )
  953. return maxCost;
  954. int i;
  955. // for each row in the matrix ...
  956. for(i=1; i<n0+1; ++i)
  957. {
  958. int ii = i * mtxMaxN; // current row
  959. int i_1 = ii - mtxMaxN; // previous row
  960. // Limit the row to (2*iMaxCost)+1 diagnal strip.
  961. // This strip is based on the idea that the best case can be precomputed for
  962. // all matrix elements in advance - where the best case for position i,j is:
  963. // abs(i-j). This can be justified based on the idea that the least possible
  964. // distance between two strings of length i and j is abs(i-1). The minimum least
  965. // possible distance is therefore found on the matrix diagnal and grows as the
  966. // distance from the diagnal increases.
  967. int ji = cmMax( 1, i - iMaxCost );
  968. int jn = cmMin(iMaxCost + i, n1) + 1;
  969. int j;
  970. // fill in (max cost + 1) as the value in the column before the starting column
  971. // (it will be referred to during the first computation in this row)
  972. if( ji >= 2 )
  973. m[ ii + (ji-1) ] = iMaxCost + 1;
  974. // for each column in the diagnal stripe - beginning with the leftmost column.
  975. for( j=ji; j<jn; ++j)
  976. {
  977. int cost = s0[i-1] == s1[j-1] ? 0 : 1;
  978. m[ ii + j ] = v = cmMin( m[ i_1 + j] + 1, cmMin( m[ ii + j - 1] + 1, m[ i_1 + j - 1 ] + cost ) );
  979. }
  980. // fill in (max cost + 1) in the column following the last column
  981. // (it will be referred to during computation of the following row)
  982. if( j < n1+1 )
  983. m[ii + j] = iMaxCost + 1;
  984. }
  985. assert( v >= 0 );
  986. return cmMin( maxCost , (double) v / maxN);
  987. }
  988. #endif