binarize.c 22 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757
  1. /*Copyright (C) 2008-2009 Timothy B. Terriberry (tterribe@xiph.org)
  2. You can redistribute this library and/or modify it under the terms of the
  3. GNU Lesser General Public License as published by the Free Software
  4. Foundation; either version 2.1 of the License, or (at your option) any later
  5. version.*/
  6. //#include <stdlib.h>
  7. //#include <math.h>
  8. //#include <string.h>
  9. #include "util.h"
  10. #include "image.h"
  11. #include "binarize.h"
  12. #if 0
  13. /*Binarization based on~\cite{GPP06}.
  14. @ARTICLE{GPP06,
  15. author="Basilios Gatos and Ioannis E. Pratikakis and Stavros J. Perantonis",
  16. title="Adaptive Degraded Document Image Binarization",
  17. journal="Pattern Recognition",
  18. volume=39,
  19. number=3,
  20. pages="317-327",
  21. month=Mar,
  22. year=2006
  23. }*/
  24. #if 0
  25. /*Applies a 5x5 Wiener filter to the image, in-place, emphasizing differences
  26. where the local variance is small, and de-emphasizing them where it is
  27. large.*/
  28. void qr_wiener_filter(unsigned char *_img,int _width,int _height){
  29. unsigned *m_buf[8];
  30. unsigned *sn2_buf[8];
  31. unsigned char g;
  32. int x;
  33. int y;
  34. if(_width<=0||_height<=0)return;
  35. m_buf[0]=(unsigned *)malloc((_width+4<<3)*sizeof(*m_buf));
  36. sn2_buf[0]=(unsigned *)malloc((_width+4<<3)*sizeof(*sn2_buf));
  37. for(y=1;y<8;y++){
  38. m_buf[y]=m_buf[y-1]+_width+4;
  39. sn2_buf[y]=sn2_buf[y-1]+_width+4;
  40. }
  41. for(y=-4;y<_height;y++){
  42. unsigned *pm;
  43. unsigned *psn2;
  44. int i;
  45. int j;
  46. pm=m_buf[y+2&7];
  47. psn2=sn2_buf[y+2&7];
  48. for(x=-4;x<_width;x++){
  49. unsigned m;
  50. unsigned m2;
  51. m=m2=0;
  52. if(y>=0&&y<_height-4&&x>=0&&x<_width-4)for(i=0;i<5;i++)for(j=0;j<5;j++){
  53. g=_img[(y+i)*_width+x+j];
  54. m+=g;
  55. m2+=g*g;
  56. }
  57. else for(i=0;i<5;i++)for(j=0;j<5;j++){
  58. g=_img[QR_CLAMPI(0,y+i,_height-1)*_width+QR_CLAMPI(0,x+j,_width-1)];
  59. m+=g;
  60. m2+=g*g;
  61. }
  62. pm[x+4]=m;
  63. psn2[x+4]=(m2*25-m*m);
  64. }
  65. pm=m_buf[y&7];
  66. if(y>=0)for(x=0;x<_width;x++){
  67. int sn2;
  68. sn2=sn2_buf[y&7][x+2];
  69. if(sn2){
  70. int vn3;
  71. int m;
  72. /*Gatos et al. give the expression
  73. mu+(s2-v2)*(g-mu)/s2 ,
  74. which we reduce to
  75. mu+(s2-v2)*g/s2-(s2-v2)*mu/s2 ,
  76. g-(v2/s2)*g+(v2/s2)*mu ,
  77. g+(mu-g)*(v2/s2) .
  78. However, s2 is much noisier than v2, and dividing by it often gives
  79. extremely large adjustments, causing speckle near edges.
  80. Therefore we limit the ratio (v2/s2) to lie between 0 and 1.*/
  81. vn3=0;
  82. for(i=-2;i<3;i++){
  83. psn2=sn2_buf[y+i&7];
  84. for(j=0;j<5;j++)vn3+=psn2[x+j];
  85. }
  86. m=m_buf[y&7][x+2];
  87. vn3=vn3+1023>>10;
  88. sn2=25*sn2+1023>>10;
  89. if(vn3<sn2){
  90. int a;
  91. g=_img[y*_width+x];
  92. a=(m-25*g)*vn3;
  93. sn2*=25;
  94. _img[y*_width+x]=QR_CLAMP255(g+QR_DIVROUND(a,sn2));
  95. }
  96. else _img[y*_width+x]=(unsigned char)(((m<<1)+25)/50);
  97. }
  98. }
  99. }
  100. free(sn2_buf[0]);
  101. free(m_buf[0]);
  102. }
  103. #else
  104. /*Applies a 3x3 Wiener filter to the image, in-place, emphasizing differences
  105. where the local variance is small, and de-emphasizing them where it is
  106. large.*/
  107. void qr_wiener_filter(unsigned char *_img,int _width,int _height){
  108. unsigned *m_buf[4];
  109. unsigned *sn2_buf[4];
  110. unsigned char g;
  111. int x;
  112. int y;
  113. if(_width<=0||_height<=0)return;
  114. m_buf[0]=(unsigned *)malloc((_width+2<<2)*sizeof(*m_buf));
  115. sn2_buf[0]=(unsigned *)malloc((_width+2<<2)*sizeof(*sn2_buf));
  116. for(y=1;y<4;y++){
  117. m_buf[y]=m_buf[y-1]+_width+2;
  118. sn2_buf[y]=sn2_buf[y-1]+_width+2;
  119. }
  120. for(y=-2;y<_height;y++){
  121. unsigned *pm;
  122. unsigned *psn2;
  123. int i;
  124. int j;
  125. pm=m_buf[y+1&3];
  126. psn2=sn2_buf[y+1&3];
  127. for(x=-2;x<_width;x++){
  128. unsigned m;
  129. unsigned m2;
  130. m=m2=0;
  131. if(y>=0&&y<_height-2&&x>=0&&x<_width-2)for(i=0;i<3;i++)for(j=0;j<3;j++){
  132. g=_img[(y+i)*_width+x+j];
  133. m+=g;
  134. m2+=g*g;
  135. }
  136. else for(i=0;i<3;i++)for(j=0;j<3;j++){
  137. g=_img[QR_CLAMPI(0,y+i,_height-1)*_width+QR_CLAMPI(0,x+j,_width-1)];
  138. m+=g;
  139. m2+=g*g;
  140. }
  141. pm[x+2]=m;
  142. psn2[x+2]=(m2*9-m*m);
  143. }
  144. pm=m_buf[y&3];
  145. if(y>=0)for(x=0;x<_width;x++){
  146. int sn2;
  147. sn2=sn2_buf[y&3][x+1];
  148. if(sn2){
  149. int m;
  150. int vn3;
  151. /*Gatos et al. give the expression
  152. mu+(s2-v2)*(g-mu)/s2 ,
  153. which we reduce to
  154. mu+(s2-v2)*g/s2-(s2-v2)*mu/s2 ,
  155. g-(v2/s2)*g+(v2/s2)*mu ,
  156. g+(mu-g)*(v2/s2) .
  157. However, s2 is much noisier than v2, and dividing by it often gives
  158. extremely large adjustments, causing speckle near edges.
  159. Therefore we limit the ratio (v2/s2) to lie between 0 and 1.*/
  160. vn3=0;
  161. for(i=-1;i<2;i++){
  162. psn2=sn2_buf[y+i&3];
  163. for(j=0;j<3;j++)vn3+=psn2[x+j];
  164. }
  165. m=m_buf[y&3][x+1];
  166. vn3=vn3+31>>5;
  167. sn2=9*sn2+31>>5;
  168. if(vn3<sn2){
  169. int a;
  170. g=_img[y*_width+x];
  171. a=m-9*g;
  172. sn2*=9;
  173. _img[y*_width+x]=QR_CLAMP255(g+QR_DIVROUND(a,sn2));
  174. }
  175. else _img[y*_width+x]=(unsigned char)(((m<<1)+9)/18);
  176. }
  177. }
  178. }
  179. free(sn2_buf[0]);
  180. free(m_buf[0]);
  181. }
  182. #endif
  183. /*Computes a (conservative) foreground mask using the adaptive binarization
  184. threshold given in~\cite{SP00}, but knocking the threshold parameter down to
  185. k=0.2.
  186. Note on dynamic range: we assume _width*_height<=0x1000000 (24 bits).
  187. Returns the average background value.
  188. @ARTICLE{SP00,
  189. author="Jaakko J. Sauvola and Matti Pietik\"{a}inen",
  190. title="Adaptive Document Image Binarization",
  191. volume=33,
  192. number=2,
  193. pages="225--236",
  194. month=Feb,
  195. year=2000
  196. }*/
  197. static void qr_sauvola_mask(unsigned char *_mask,unsigned *_b,int *_nb,
  198. const unsigned char *_img,int _width,int _height){
  199. unsigned b;
  200. int nb;
  201. b=0;
  202. nb=0;
  203. if(_width>0&&_height>0){
  204. unsigned *col_sums;
  205. unsigned *col2_sums;
  206. int logwindw;
  207. int logwindh;
  208. int windw;
  209. int windh;
  210. int y0offs;
  211. int y1offs;
  212. unsigned g;
  213. unsigned g2;
  214. int x;
  215. int y;
  216. /*We keep the window size fairly large to ensure it doesn't fit completely
  217. inside the center of a finder pattern of a version 1 QR code at full
  218. resolution.*/
  219. for(logwindw=4;logwindw<8&&(1<<logwindw)<(_width+7>>3);logwindw++);
  220. for(logwindh=4;logwindh<8&&(1<<logwindh)<(_height+7>>3);logwindh++);
  221. windw=1<<logwindw;
  222. windh=1<<logwindh;
  223. col_sums=(unsigned *)malloc(_width*sizeof(*col_sums));
  224. col2_sums=(unsigned *)malloc(_width*sizeof(*col2_sums));
  225. /*Initialize sums down each column.*/
  226. for(x=0;x<_width;x++){
  227. g=_img[x];
  228. g2=g*g;
  229. col_sums[x]=(g<<logwindh-1)+g;
  230. col2_sums[x]=(g2<<logwindh-1)+g2;
  231. }
  232. for(y=1;y<(windh>>1);y++){
  233. y1offs=QR_MINI(y,_height-1)*_width;
  234. for(x=0;x<_width;x++){
  235. g=_img[y1offs+x];
  236. col_sums[x]+=g;
  237. col2_sums[x]+=g*g;
  238. }
  239. }
  240. for(y=0;y<_height;y++){
  241. unsigned m;
  242. unsigned m2;
  243. int x0;
  244. int x1;
  245. /*Initialize the sums over the window.*/
  246. m=(col_sums[0]<<logwindw-1)+col_sums[0];
  247. m2=(col2_sums[0]<<logwindw-1)+col2_sums[0];
  248. for(x=1;x<(windw>>1);x++){
  249. x1=QR_MINI(x,_width-1);
  250. m+=col_sums[x1];
  251. m2+=col2_sums[x1];
  252. }
  253. for(x=0;x<_width;x++){
  254. int d;
  255. /*Perform the test against the threshold T = (m/n)*(1+k*(s/R-1)),
  256. where n=windw*windh, s=sqrt((m2-(m*m)/n)/n), and R=128.
  257. We don't actually compute the threshold directly, as that would
  258. require a square root.
  259. Instead we perform the equivalent test:
  260. (m/n)*(m/n)*(m2/n-(m/n)*(m/n))/16 > (((1/k)*g-((1-k)/k)*(m/n))*32)**2
  261. R is split up across each side of the inequality to maximize the
  262. dynamic range available for the right hand side, which requires
  263. 31 bits in the worst case.*/
  264. /*(m/n)*(1+(1/5)*(sqrt((m2-m*m/n)/n)/128-1)) > g
  265. m*(1+(1/5)*(sqrt((m2-m*m/n)/n)/128-1)) > g*n
  266. m*sqrt((m2-m*m/n)/n) > 5*g*n-4*m<<7
  267. m*m*(m2*n-m*m) > (5*g*n-4*m<<7)**2*n*n || 5*g*n-4*m < 0 */
  268. g=_img[y*_width+x];
  269. d=(5*g<<logwindw+logwindh)-4*m;
  270. if(d>=0){
  271. unsigned mm;
  272. unsigned mms2;
  273. unsigned d2;
  274. mm=(m>>logwindw)*(m>>logwindh);
  275. mms2=(m2-mm>>logwindw+logwindh)*(mm>>logwindw+logwindh)+15>>4;
  276. d2=d>>logwindw+logwindh-5;
  277. d2*=d2;
  278. if(d2>=mms2){
  279. /*Update the background average.*/
  280. b+=g;
  281. nb++;
  282. _mask[y*_width+x]=0;
  283. }
  284. else _mask[y*_width+x]=0xFF;
  285. }
  286. else _mask[y*_width+x]=0xFF;
  287. /*Update the window sums.*/
  288. if(x+1<_width){
  289. x0=QR_MAXI(0,x-(windw>>1));
  290. x1=QR_MINI(x+(windw>>1),_width-1);
  291. m+=col_sums[x1]-col_sums[x0];
  292. m2+=col2_sums[x1]-col2_sums[x0];
  293. }
  294. }
  295. /*Update the column sums.*/
  296. if(y+1<_height){
  297. y0offs=QR_MAXI(0,y-(windh>>1))*_width;
  298. y1offs=QR_MINI(y+(windh>>1),_height-1)*_width;
  299. for(x=0;x<_width;x++){
  300. g=_img[y0offs+x];
  301. col_sums[x]-=g;
  302. col2_sums[x]-=g*g;
  303. g=_img[y1offs+x];
  304. col_sums[x]+=g;
  305. col2_sums[x]+=g*g;
  306. }
  307. }
  308. }
  309. free(col2_sums);
  310. free(col_sums);
  311. }
  312. *_b=b;
  313. *_nb=nb;
  314. }
  315. /*Interpolates a background image given the source and a conservative
  316. foreground mask.
  317. If the current window contains no foreground pixels, the average background
  318. value over the whole image is used.
  319. Note on dynamic range: we assume _width*_height<=0x8000000 (23 bits).
  320. Returns the average difference between the foreground and the interpolated
  321. background.*/
  322. static void qr_interpolate_background(unsigned char *_dst,
  323. int *_delta,int *_ndelta,const unsigned char *_img,const unsigned char *_mask,
  324. int _width,int _height,unsigned _b,int _nb){
  325. int delta;
  326. int ndelta;
  327. delta=ndelta=0;
  328. if(_width>0&&_height>0){
  329. unsigned *col_sums;
  330. unsigned *ncol_sums;
  331. int logwindw;
  332. int logwindh;
  333. int windw;
  334. int windh;
  335. int y0offs;
  336. int y1offs;
  337. unsigned b;
  338. unsigned g;
  339. int x;
  340. int y;
  341. b=_nb>0?((_b<<1)+_nb)/(_nb<<1):0xFF;
  342. for(logwindw=4;logwindw<8&&(1<<logwindw)<(_width+15>>4);logwindw++);
  343. for(logwindh=4;logwindh<8&&(1<<logwindh)<(_height+15>>4);logwindh++);
  344. windw=1<<logwindw;
  345. windh=1<<logwindh;
  346. col_sums=(unsigned *)malloc(_width*sizeof(*col_sums));
  347. ncol_sums=(unsigned *)malloc(_width*sizeof(*ncol_sums));
  348. /*Initialize sums down each column.*/
  349. for(x=0;x<_width;x++){
  350. if(!_mask[x]){
  351. g=_img[x];
  352. col_sums[x]=(g<<logwindh-1)+g;
  353. ncol_sums[x]=(1<<logwindh-1)+1;
  354. }
  355. else col_sums[x]=ncol_sums[x]=0;
  356. }
  357. for(y=1;y<(windh>>1);y++){
  358. y1offs=QR_MINI(y,_height-1)*_width;
  359. for(x=0;x<_width;x++)if(!_mask[y1offs+x]){
  360. col_sums[x]+=_img[y1offs+x];
  361. ncol_sums[x]++;
  362. }
  363. }
  364. for(y=0;y<_height;y++){
  365. unsigned n;
  366. unsigned m;
  367. int x0;
  368. int x1;
  369. /*Initialize the sums over the window.*/
  370. m=(col_sums[0]<<logwindw-1)+col_sums[0];
  371. n=(ncol_sums[0]<<logwindw-1)+ncol_sums[0];
  372. for(x=1;x<(windw>>1);x++){
  373. x1=QR_MINI(x,_width-1);
  374. m+=col_sums[x1];
  375. n+=ncol_sums[x1];
  376. }
  377. for(x=0;x<_width;x++){
  378. if(!_mask[y*_width+x])g=_img[y*_width+x];
  379. else{
  380. g=n>0?((m<<1)+n)/(n<<1):b;
  381. delta+=(int)g-_img[y*_width+x];
  382. ndelta++;
  383. }
  384. _dst[y*_width+x]=(unsigned char)g;
  385. /*Update the window sums.*/
  386. if(x+1<_width){
  387. x0=QR_MAXI(0,x-(windw>>1));
  388. x1=QR_MINI(x+(windw>>1),_width-1);
  389. m+=col_sums[x1]-col_sums[x0];
  390. n+=ncol_sums[x1]-ncol_sums[x0];
  391. }
  392. }
  393. /*Update the column sums.*/
  394. if(y+1<_height){
  395. y0offs=QR_MAXI(0,y-(windh>>1))*_width;
  396. y1offs=QR_MINI(y+(windh>>1),_height-1)*_width;
  397. for(x=0;x<_width;x++){
  398. if(!_mask[y0offs+x]){
  399. col_sums[x]-=_img[y0offs+x];
  400. ncol_sums[x]--;
  401. }
  402. if(!_mask[y1offs+x]){
  403. col_sums[x]+=_img[y1offs+x];
  404. ncol_sums[x]++;
  405. }
  406. }
  407. }
  408. }
  409. free(ncol_sums);
  410. free(col_sums);
  411. }
  412. *_delta=delta;
  413. *_ndelta=ndelta;
  414. }
  415. /*Parameters of the logistic sigmoid function that defines the threshold based
  416. on the background intensity.
  417. They should all be between 0 and 1.*/
  418. #define QR_GATOS_Q (0.7)
  419. #define QR_GATOS_P1 (0.5)
  420. #define QR_GATOS_P2 (0.8)
  421. /*Compute the final binarization mask according to Gatos et al.'s
  422. method~\cite{GPP06}.*/
  423. static void qr_gatos_mask(unsigned char *_mask,const unsigned char *_img,
  424. const unsigned char *_background,int _width,int _height,
  425. unsigned _b,int _nb,int _delta,int _ndelta){
  426. unsigned thresh[256];
  427. unsigned g;
  428. double delta;
  429. double b;
  430. int x;
  431. int y;
  432. /*Construct a lookup table for the thresholds.
  433. This bit uses floating point, but doesn't need to do much calculation, so
  434. emulation should be fine.*/
  435. b=_nb>0?(_b+0.5)/_nb:0xFF;
  436. delta=_ndelta>0?(_delta+0.5)/_ndelta:0xFF;
  437. for(g=0;g<256;g++){
  438. double d;
  439. d=QR_GATOS_Q*delta*(QR_GATOS_P2+(1-QR_GATOS_P2)/
  440. (1+exp(2*(1+QR_GATOS_P1)/(1-QR_GATOS_P1)-4*g/(b*(1-QR_GATOS_P1)))));
  441. if(d<1)d=1;
  442. else if(d>0xFF)d=0xFF;
  443. thresh[g]=(unsigned)floor(d);
  444. }
  445. /*Apply the adaptive threshold.*/
  446. for(y=0;y<_height;y++)for(x=0;x<_width;x++){
  447. g=_background[y*_width+x];
  448. /*_background[y*_width+x]=thresh[g];*/
  449. _mask[y*_width+x]=(unsigned char)(-(g-_img[y*_width+x]>thresh[g])&0xFF);
  450. }
  451. /*{
  452. FILE *fout;
  453. fout=fopen("thresh.png","wb");
  454. image_write_png(_background,_width,_height,fout);
  455. fclose(fout);
  456. }*/
  457. }
  458. /*Binarizes a grayscale image.*/
  459. void qr_binarize(unsigned char *_img,int _width,int _height){
  460. unsigned char *mask;
  461. unsigned char *background;
  462. unsigned b;
  463. int nb;
  464. int delta;
  465. int ndelta;
  466. /*qr_wiener_filter(_img,_width,_height);
  467. {
  468. FILE *fout;
  469. fout=fopen("wiener.png","wb");
  470. image_write_png(_img,_width,_height,fout);
  471. fclose(fout);
  472. }*/
  473. mask=(unsigned char *)malloc(_width*_height*sizeof(*mask));
  474. qr_sauvola_mask(mask,&b,&nb,_img,_width,_height);
  475. /*{
  476. FILE *fout;
  477. fout=fopen("foreground.png","wb");
  478. image_write_png(mask,_width,_height,fout);
  479. fclose(fout);
  480. }*/
  481. background=(unsigned char *)malloc(_width*_height*sizeof(*mask));
  482. qr_interpolate_background(background,&delta,&ndelta,
  483. _img,mask,_width,_height,b,nb);
  484. /*{
  485. FILE *fout;
  486. fout=fopen("background.png","wb");
  487. image_write_png(background,_width,_height,fout);
  488. fclose(fout);
  489. }*/
  490. qr_gatos_mask(_img,_img,background,_width,_height,b,nb,delta,ndelta);
  491. free(background);
  492. free(mask);
  493. }
  494. #else
  495. /*The above algorithms are computationally expensive, and do not work as well
  496. as the simple algorithm below.
  497. Sauvola by itself does an excellent job of classifying regions outside the
  498. QR code as background, which greatly reduces the chance of false alarms.
  499. However, it also tends to over-shrink isolated black dots inside the code,
  500. making them easy to miss with even slight mis-alignment.
  501. Since the Gatos method uses Sauvola as input to its background interpolation
  502. method, it cannot possibly mark any pixels as foreground which Sauvola
  503. classified as background, and thus suffers from the same problem.
  504. The following simple adaptive threshold method does not have this problem,
  505. though it produces essentially random noise outside the QR code region.
  506. QR codes are structured well enough that this does not seem to lead to any
  507. actual false alarms in practice, and it allows many more codes to be
  508. detected and decoded successfully than the Sauvola or Gatos binarization
  509. methods.*/
  510. /*A simplified adaptive thresholder.
  511. This compares the current pixel value to the mean value of a (large) window
  512. surrounding it.*/
  513. unsigned char *qr_binarize(const unsigned char *_img,int _width,int _height){
  514. unsigned char *mask = NULL;
  515. if(_width>0&&_height>0){
  516. unsigned *col_sums;
  517. int logwindw;
  518. int logwindh;
  519. int windw;
  520. int windh;
  521. int y0offs;
  522. int y1offs;
  523. unsigned g;
  524. int x;
  525. int y;
  526. mask=(unsigned char *)malloc(_width*_height*sizeof(*mask));
  527. /*We keep the window size fairly large to ensure it doesn't fit completely
  528. inside the center of a finder pattern of a version 1 QR code at full
  529. resolution.*/
  530. for(logwindw=4;logwindw<8&&(1<<logwindw)<(_width+7>>3);logwindw++);
  531. for(logwindh=4;logwindh<8&&(1<<logwindh)<(_height+7>>3);logwindh++);
  532. windw=1<<logwindw;
  533. windh=1<<logwindh;
  534. col_sums=(unsigned *)malloc(_width*sizeof(*col_sums));
  535. /*Initialize sums down each column.*/
  536. for(x=0;x<_width;x++){
  537. g=_img[x];
  538. col_sums[x]=(g<<logwindh-1)+g;
  539. }
  540. for(y=1;y<(windh>>1);y++){
  541. y1offs=QR_MINI(y,_height-1)*_width;
  542. for(x=0;x<_width;x++){
  543. g=_img[y1offs+x];
  544. col_sums[x]+=g;
  545. }
  546. }
  547. for(y=0;y<_height;y++){
  548. unsigned m;
  549. int x0;
  550. int x1;
  551. /*Initialize the sum over the window.*/
  552. m=(col_sums[0]<<logwindw-1)+col_sums[0];
  553. for(x=1;x<(windw>>1);x++){
  554. x1=QR_MINI(x,_width-1);
  555. m+=col_sums[x1];
  556. }
  557. for(x=0;x<_width;x++){
  558. /*Perform the test against the threshold T = (m/n)-D,
  559. where n=windw*windh and D=3.*/
  560. g=_img[y*_width+x];
  561. mask[y*_width+x]=-(g+3<<logwindw+logwindh<m)&0xFF;
  562. /*Update the window sum.*/
  563. if(x+1<_width){
  564. x0=QR_MAXI(0,x-(windw>>1));
  565. x1=QR_MINI(x+(windw>>1),_width-1);
  566. m+=col_sums[x1]-col_sums[x0];
  567. }
  568. }
  569. /*Update the column sums.*/
  570. if(y+1<_height){
  571. y0offs=QR_MAXI(0,y-(windh>>1))*_width;
  572. y1offs=QR_MINI(y+(windh>>1),_height-1)*_width;
  573. for(x=0;x<_width;x++){
  574. col_sums[x]-=_img[y0offs+x];
  575. col_sums[x]+=_img[y1offs+x];
  576. }
  577. }
  578. }
  579. free(col_sums);
  580. }
  581. #if defined(QR_DEBUG)
  582. {
  583. FILE *fout;
  584. fout=fopen("binary.png","wb");
  585. image_write_png(_img,_width,_height,fout);
  586. fclose(fout);
  587. }
  588. #endif
  589. return(mask);
  590. }
  591. unsigned char *qr_binarize2(unsigned char *_img, int _width, int _height) {
  592. unsigned char *mask = NULL;
  593. if (_width > 0 && _height > 0) {
  594. unsigned *col_sums;
  595. int logwindw;
  596. int logwindh;
  597. int windw;
  598. int windh;
  599. int y0offs;
  600. int y1offs;
  601. unsigned g;
  602. int x;
  603. int y;
  604. int r_line = 0;
  605. int w_line = 0;
  606. int line_cnt = 0;
  607. int line_max = 0;
  608. int line_offset = 0;
  609. /*We keep the window size fairly large to ensure it doesn't fit completely
  610. inside the center of a finder pattern of a version 1 QR code at full
  611. resolution.*/
  612. for (logwindw = 4; logwindw < 8 && (1 << logwindw) < (_width + 7 >> 3); logwindw++);
  613. for (logwindh = 4; logwindh < 8 && (1 << logwindh) < (_height + 7 >> 3); logwindh++);
  614. windw = 1 << logwindw;
  615. windh = 1 << logwindh;
  616. col_sums = (unsigned *)malloc(_width * sizeof(*col_sums));
  617. mask = (unsigned char *)malloc(((windh >> 1) + 2)*_width * sizeof(*mask));
  618. line_max = (windh >> 1) + 2;
  619. /*Initialize sums down each column.*/
  620. for (x = 0; x < _width; x++) {
  621. g = _img[x];
  622. col_sums[x] = (g << logwindh - 1) + g;
  623. }
  624. for (y = 1; y < (windh >> 1); y++) {
  625. y1offs = QR_MINI(y, _height - 1)*_width;
  626. for (x = 0; x < _width; x++) {
  627. g = _img[y1offs + x];
  628. col_sums[x] += g;
  629. }
  630. }
  631. for (y = 0; y < _height; y++) {
  632. unsigned m;
  633. int x0;
  634. int x1;
  635. /*Initialize the sum over the window.*/
  636. m = (col_sums[0] << logwindw - 1) + col_sums[0];
  637. for (x = 1; x < (windw >> 1); x++) {
  638. x1 = QR_MINI(x, _width - 1);
  639. m += col_sums[x1];
  640. }
  641. for (x = 0; x < _width; x++) {
  642. /*Perform the test against the threshold T = (m/n)-D,
  643. where n=windw*windh and D=3.*/
  644. g = _img[y*_width + x];
  645. mask[w_line*_width + x] = -(g + 3 << logwindw + logwindh < m) & 0xFF;
  646. /*Update the window sum.*/
  647. if (x + 1 < _width) {
  648. x0 = QR_MAXI(0, x - (windw >> 1));
  649. x1 = QR_MINI(x + (windw >> 1), _width - 1);
  650. m += col_sums[x1] - col_sums[x0];
  651. }
  652. }
  653. w_line++;
  654. line_cnt++;
  655. if (w_line >= line_max)
  656. {
  657. w_line = 0;
  658. }
  659. /*Update the column sums.*/
  660. if (y + 1 < _height) {
  661. y0offs = QR_MAXI(0, y - (windh >> 1))*_width;
  662. y1offs = QR_MINI(y + (windh >> 1), _height - 1)*_width;
  663. for (x = 0; x < _width; x++) {
  664. col_sums[x] -= _img[y0offs + x];
  665. col_sums[x] += _img[y1offs + x];
  666. }
  667. if (y0offs)
  668. {
  669. memcpy(&_img[line_offset * _width], &mask[r_line *_width], _width);
  670. r_line++;
  671. line_offset++;
  672. line_cnt--;
  673. if (r_line >= line_max)
  674. {
  675. r_line = 0;
  676. }
  677. }
  678. }
  679. else
  680. {
  681. for (int i = 0; i < line_cnt; i++)
  682. {
  683. memcpy(&_img[line_offset * _width], &mask[r_line *_width], _width);
  684. r_line++;
  685. line_offset++;
  686. if (r_line >= line_max)
  687. {
  688. r_line = 0;
  689. }
  690. }
  691. }
  692. }
  693. free(col_sums);
  694. free(mask);
  695. mask = _img;
  696. }
  697. #if defined(QR_DEBUG)
  698. {
  699. FILE *fout;
  700. fout = fopen("binary.png", "wb");
  701. image_write_png(_img, _width, _height, fout);
  702. fclose(fout);
  703. }
  704. #endif
  705. return(mask);
  706. }
  707. #endif
  708. #if defined(TEST_BINARIZE)
  709. #include <stdio.h>
  710. #include "image.c"
  711. int main(int _argc,char **_argv){
  712. unsigned char *img;
  713. int width;
  714. int height;
  715. int x;
  716. int y;
  717. if(_argc<2){
  718. fprintf(stderr,"usage: %s <image>.png\n",_argv[0]);
  719. return EXIT_FAILURE;
  720. }
  721. /*width=1182;
  722. height=1181;
  723. img=(unsigned char *)malloc(width*height*sizeof(*img));
  724. for(y=0;y<height;y++)for(x=0;x<width;x++){
  725. img[y*width+x]=(unsigned char)(-((x&1)^(y&1))&0xFF);
  726. }*/
  727. {
  728. FILE *fin;
  729. fin=fopen(_argv[1],"rb");
  730. image_read_png(&img,&width,&height,fin);
  731. fclose(fin);
  732. }
  733. qr_binarize(img,width,height);
  734. /*{
  735. FILE *fout;
  736. fout=fopen("binary.png","wb");
  737. image_write_png(img,width,height,fout);
  738. fclose(fout);
  739. }*/
  740. free(img);
  741. return EXIT_SUCCESS;
  742. }
  743. #endif