qrdec.c

00001 /*Copyright (C) 2008-2009  Timothy B. Terriberry (tterribe@xiph.org)
00002   You can redistribute this library and/or modify it under the terms of the
00003    GNU Lesser General Public License as published by the Free Software
00004    Foundation; either version 2.1 of the License, or (at your option) any later
00005    version.*/
00006 #include <stdlib.h>
00007 #include <limits.h>
00008 #include <string.h>
00009 #include <time.h>
00010 #include "qrcode.h"
00011 #include "qrdec.h"
00012 #include "bch15_5.h"
00013 #include "rs.h"
00014 #include "isaac.h"
00015 #include "util.h"
00016 #include "binarize.h"
00017 #include "image.h"
00018 #include "error.h"
00019 #include "svg.h"
00020 
00021 typedef int qr_line[3];
00022 
00023 typedef struct qr_finder_cluster qr_finder_cluster;
00024 typedef struct qr_finder_edge_pt  qr_finder_edge_pt;
00025 typedef struct qr_finder_center   qr_finder_center;
00026 
00027 typedef struct qr_aff qr_aff;
00028 typedef struct qr_hom qr_hom;
00029 
00030 typedef struct qr_finder qr_finder;
00031 
00032 typedef struct qr_hom_cell      qr_hom_cell;
00033 typedef struct qr_sampling_grid qr_sampling_grid;
00034 typedef struct qr_pack_buf      qr_pack_buf;
00035 
00036 /*The number of bits in an int.
00037   Note the cast to (int): this prevents this value from "promoting" whole
00038    expressions to an (unsigned) size_t.*/
00039 #define QR_INT_BITS    ((int)sizeof(int)*CHAR_BIT)
00040 #define QR_INT_LOGBITS (QR_ILOG(QR_INT_BITS))
00041 
00042 /*A 14 bit resolution for a homography ensures that the ideal module size for a
00043    version 40 code differs from that of a version 39 code by at least 2.*/
00044 #define QR_HOM_BITS (14)
00045 
00046 /*The number of bits of sub-module precision to use when searching for
00047    alignment patterns.
00048   Two bits allows an alignment pattern to be found even if the modules have
00049    been eroded by up to 50% (due to blurring, etc.).
00050   This must be at least one, since it affects the dynamic range of the
00051    transforms, and we sample at half-module resolution to compute a bounding
00052    quadrilateral for the code.*/
00053 #define QR_ALIGN_SUBPREC (2)
00054 
00055 
00056 /* collection of finder lines */
00057 typedef struct qr_finder_lines {
00058     qr_finder_line *lines;
00059     int nlines, clines;
00060 } qr_finder_lines;
00061 
00062 
00063 struct qr_reader {
00064     /*The GF(256) representation used in Reed-Solomon decoding.*/
00065     rs_gf256  gf;
00066     /*The random number generator used by RANSAC.*/
00067     isaac_ctx isaac;
00068     /* current finder state, horizontal and vertical lines */
00069     qr_finder_lines finder_lines[2];
00070 };
00071 
00072 
00073 /*Initializes a client reader handle.*/
00074 static void qr_reader_init (qr_reader *reader)
00075 {
00076     /*time_t now;
00077       now=time(NULL);
00078       isaac_init(&_reader->isaac,&now,sizeof(now));*/
00079     isaac_init(&reader->isaac, NULL, 0);
00080     rs_gf256_init(&reader->gf, QR_PPOLY);
00081 }
00082 
00083 /*Allocates a client reader handle.*/
00084 qr_reader *_zbar_qr_create (void)
00085 {
00086     qr_reader *reader = (qr_reader*)calloc(1, sizeof(*reader));
00087     qr_reader_init(reader);
00088     return(reader);
00089 }
00090 
00091 /*Frees a client reader handle.*/
00092 void _zbar_qr_destroy (qr_reader *reader)
00093 {
00094     zprintf(1, "max finder lines = %dx%d\n",
00095             reader->finder_lines[0].clines,
00096             reader->finder_lines[1].clines);
00097     if(reader->finder_lines[0].lines)
00098         free(reader->finder_lines[0].lines);
00099     if(reader->finder_lines[1].lines)
00100         free(reader->finder_lines[1].lines);
00101     free(reader);
00102 }
00103 
00104 /* reset finder state between scans */
00105 void _zbar_qr_reset (qr_reader *reader)
00106 {
00107     reader->finder_lines[0].nlines = 0;
00108     reader->finder_lines[1].nlines = 0;
00109 }
00110 
00111 
00112 /*A cluster of lines crossing a finder pattern (all in the same direction).*/
00113 struct qr_finder_cluster{
00114   /*Pointers to the lines crossing the pattern.*/
00115   qr_finder_line **lines;
00116   /*The number of lines in the cluster.*/
00117   int              nlines;
00118 };
00119 
00120 
00121 /*A point on the edge of a finder pattern.
00122   These are obtained from the endpoints of the lines crossing this particular
00123    pattern.*/
00124 struct qr_finder_edge_pt{
00125   /*The location of the edge point.*/
00126   qr_point pos;
00127   /*A label classifying which edge this belongs to:
00128     0: negative u edge (left)
00129     1: positive u edge (right)
00130     2: negative v edge (top)
00131     3: positive v edge (bottom)*/
00132   int      edge;
00133   /*The (signed) perpendicular distance of the edge point from a line parallel
00134      to the edge passing through the finder center, in (u,v) coordinates.
00135     This is also re-used by RANSAC to store inlier flags.*/
00136   int      extent;
00137 };
00138 
00139 
00140 /*The center of a finder pattern obtained from the crossing of one or more
00141    clusters of horizontal finder lines with one or more clusters of vertical
00142    finder lines.*/
00143 struct qr_finder_center{
00144   /*The estimated location of the finder center.*/
00145   qr_point           pos;
00146   /*The list of edge points from the crossing lines.*/
00147   qr_finder_edge_pt *edge_pts;
00148   /*The number of edge points from the crossing lines.*/
00149   int                nedge_pts;
00150 };
00151 
00152 
00153 static int qr_finder_vline_cmp(const void *_a,const void *_b){
00154   const qr_finder_line *a;
00155   const qr_finder_line *b;
00156   a=(const qr_finder_line *)_a;
00157   b=(const qr_finder_line *)_b;
00158   return ((a->pos[0]>b->pos[0])-(a->pos[0]<b->pos[0])<<1)+
00159    (a->pos[1]>b->pos[1])-(a->pos[1]<b->pos[1]);
00160 }
00161 
00162 /*Clusters adjacent lines into groups that are large enough to be crossing a
00163    finder pattern (relative to their length).
00164   _clusters:  The buffer in which to store the clusters found.
00165   _neighbors: The buffer used to store the lists of lines in each cluster.
00166   _lines:     The list of lines to cluster.
00167               Horizontal lines must be sorted in ascending order by Y
00168                coordinate, with ties broken by X coordinate.
00169               Vertical lines must be sorted in ascending order by X coordinate,
00170                with ties broken by Y coordinate.
00171   _nlines:    The number of lines in the set of lines to cluster.
00172   _v:         0 for horizontal lines, or 1 for vertical lines.
00173   Return: The number of clusters.*/
00174 static int qr_finder_cluster_lines(qr_finder_cluster *_clusters,
00175  qr_finder_line **_neighbors,qr_finder_line *_lines,int _nlines,int _v){
00176   unsigned char   *mark;
00177   qr_finder_line **neighbors;
00178   int              nneighbors;
00179   int              nclusters;
00180   int              i;
00181   /*TODO: Kalman filters!*/
00182   mark=(unsigned char *)calloc(_nlines,sizeof(*mark));
00183   neighbors=_neighbors;
00184   nclusters=0;
00185   for(i=0;i<_nlines-1;i++)if(!mark[i]){
00186     int len;
00187     int j;
00188     nneighbors=1;
00189     neighbors[0]=_lines+i;
00190     len=_lines[i].len;
00191     for(j=i+1;j<_nlines;j++)if(!mark[j]){
00192       const qr_finder_line *a;
00193       const qr_finder_line *b;
00194       int                   thresh;
00195       a=neighbors[nneighbors-1];
00196       b=_lines+j;
00197       /*The clustering threshold is proportional to the size of the lines,
00198          since minor noise in large areas can interrupt patterns more easily
00199          at high resolutions.*/
00200       thresh=a->len+7>>2;
00201       if(abs(a->pos[1-_v]-b->pos[1-_v])>thresh)break;
00202       if(abs(a->pos[_v]-b->pos[_v])>thresh)continue;
00203       if(abs(a->pos[_v]+a->len-b->pos[_v]-b->len)>thresh)continue;
00204       if(a->boffs>0&&b->boffs>0&&
00205        abs(a->pos[_v]-a->boffs-b->pos[_v]+b->boffs)>thresh){
00206         continue;
00207       }
00208       if(a->eoffs>0&&b->eoffs>0&&
00209        abs(a->pos[_v]+a->len+a->eoffs-b->pos[_v]-b->len-b->eoffs)>thresh){
00210         continue;
00211       }
00212       neighbors[nneighbors++]=_lines+j;
00213       len+=b->len;
00214     }
00215     /*We require at least three lines to form a cluster, which eliminates a
00216        large number of false positives, saving considerable decoding time.
00217       This should still be sufficient for 1-pixel codes with no noise.*/
00218     if(nneighbors<3)continue;
00219     /*The expected number of lines crossing a finder pattern is equal to their
00220        average length.
00221       We accept the cluster if size is at least 1/3 their average length (this
00222        is a very small threshold, but was needed for some test images).*/
00223     len=((len<<1)+nneighbors)/(nneighbors<<1);
00224     if(nneighbors*(5<<QR_FINDER_SUBPREC)>=len){
00225       _clusters[nclusters].lines=neighbors;
00226       _clusters[nclusters].nlines=nneighbors;
00227       for(j=0;j<nneighbors;j++)mark[neighbors[j]-_lines]=1;
00228       neighbors+=nneighbors;
00229       nclusters++;
00230     }
00231   }
00232   free(mark);
00233   return nclusters;
00234 }
00235 
00236 /*Adds the coordinates of the edge points from the lines contained in the
00237    given list of clusters to the list of edge points for a finder center.
00238   Only the edge point position is initialized.
00239   The edge label and extent are set by qr_finder_edge_pts_aff_classify()
00240    or qr_finder_edge_pts_hom_classify().
00241   _edge_pts:   The buffer in which to store the edge points.
00242   _nedge_pts:  The current number of edge points in the buffer.
00243   _neighbors:  The list of lines in the cluster.
00244   _nneighbors: The number of lines in the list of lines in the cluster.
00245   _v:          0 for horizontal lines and 1 for vertical lines.
00246   Return: The new total number of edge points.*/
00247 static int qr_finder_edge_pts_fill(qr_finder_edge_pt *_edge_pts,int _nedge_pts,
00248  qr_finder_cluster **_neighbors,int _nneighbors,int _v){
00249   int i;
00250   for(i=0;i<_nneighbors;i++){
00251     qr_finder_cluster *c;
00252     int                j;
00253     c=_neighbors[i];
00254     for(j=0;j<c->nlines;j++){
00255       qr_finder_line *l;
00256       l=c->lines[j];
00257       if(l->boffs>0){
00258         _edge_pts[_nedge_pts].pos[0]=l->pos[0];
00259         _edge_pts[_nedge_pts].pos[1]=l->pos[1];
00260         _edge_pts[_nedge_pts].pos[_v]-=l->boffs;
00261         _nedge_pts++;
00262       }
00263       if(l->eoffs>0){
00264         _edge_pts[_nedge_pts].pos[0]=l->pos[0];
00265         _edge_pts[_nedge_pts].pos[1]=l->pos[1];
00266         _edge_pts[_nedge_pts].pos[_v]+=l->len+l->eoffs;
00267         _nedge_pts++;
00268       }
00269     }
00270   }
00271   return _nedge_pts;
00272 }
00273 
00274 static int qr_finder_center_cmp(const void *_a,const void *_b){
00275   const qr_finder_center *a;
00276   const qr_finder_center *b;
00277   a=(const qr_finder_center *)_a;
00278   b=(const qr_finder_center *)_b;
00279   return ((b->nedge_pts>a->nedge_pts)-(b->nedge_pts<a->nedge_pts)<<2)+
00280    ((a->pos[1]>b->pos[1])-(a->pos[1]<b->pos[1])<<1)+
00281    (a->pos[0]>b->pos[0])-(a->pos[0]<b->pos[0]);
00282 }
00283 
00284 /*Determine if a horizontal line crosses a vertical line.
00285   _hline: The horizontal line.
00286   _vline: The vertical line.
00287   Return: A non-zero value if the lines cross, or zero if they do not.*/
00288 static int qr_finder_lines_are_crossing(const qr_finder_line *_hline,
00289  const qr_finder_line *_vline){
00290   return
00291    _hline->pos[0]<=_vline->pos[0]&&_vline->pos[0]<_hline->pos[0]+_hline->len&&
00292    _vline->pos[1]<=_hline->pos[1]&&_hline->pos[1]<_vline->pos[1]+_vline->len;
00293 }
00294 
00295 /*Finds horizontal clusters that cross corresponding vertical clusters,
00296    presumably corresponding to a finder center.
00297   _center:     The buffer in which to store putative finder centers.
00298   _edge_pts:   The buffer to use for the edge point lists for each finder
00299                 center.
00300   _hclusters:  The clusters of horizontal lines crossing finder patterns.
00301   _nhclusters: The number of horizontal line clusters.
00302   _vclusters:  The clusters of vertical lines crossing finder patterns.
00303   _nvclusters: The number of vertical line clusters.
00304   Return: The number of putative finder centers.*/
00305 static int qr_finder_find_crossings(qr_finder_center *_centers,
00306  qr_finder_edge_pt *_edge_pts,qr_finder_cluster *_hclusters,int _nhclusters,
00307  qr_finder_cluster *_vclusters,int _nvclusters){
00308   qr_finder_cluster **hneighbors;
00309   qr_finder_cluster **vneighbors;
00310   unsigned char      *hmark;
00311   unsigned char      *vmark;
00312   int                 ncenters;
00313   int                 i;
00314   int                 j;
00315   hneighbors=(qr_finder_cluster **)malloc(_nhclusters*sizeof(*hneighbors));
00316   vneighbors=(qr_finder_cluster **)malloc(_nvclusters*sizeof(*vneighbors));
00317   hmark=(unsigned char *)calloc(_nhclusters,sizeof(*hmark));
00318   vmark=(unsigned char *)calloc(_nvclusters,sizeof(*vmark));
00319   ncenters=0;
00320   /*TODO: This may need some re-working.
00321     We should be finding groups of clusters such that _all_ horizontal lines in
00322      _all_ horizontal clusters in the group cross _all_ vertical lines in _all_
00323      vertical clusters in the group.
00324     This is equivalent to finding the maximum bipartite clique in the
00325      connectivity graph, which requires linear progamming to solve efficiently.
00326     In principle, that is easy to do, but a realistic implementation without
00327      floating point is a lot of work (and computationally expensive).
00328     Right now we are relying on a sufficient border around the finder patterns
00329      to prevent false positives.*/
00330   for(i=0;i<_nhclusters;i++)if(!hmark[i]){
00331     qr_finder_line *a;
00332     qr_finder_line *b;
00333     int             nvneighbors;
00334     int             nedge_pts;
00335     int             y;
00336     a=_hclusters[i].lines[_hclusters[i].nlines>>1];
00337     y=nvneighbors=0;
00338     for(j=0;j<_nvclusters;j++)if(!vmark[j]){
00339       b=_vclusters[j].lines[_vclusters[j].nlines>>1];
00340       if(qr_finder_lines_are_crossing(a,b)){
00341         vmark[j]=1;
00342         y+=(b->pos[1]<<1)+b->len;
00343         if(b->boffs>0&&b->eoffs>0)y+=b->eoffs-b->boffs;
00344         vneighbors[nvneighbors++]=_vclusters+j;
00345       }
00346     }
00347     if(nvneighbors>0){
00348       qr_finder_center *c;
00349       int               nhneighbors;
00350       int               x;
00351       x=(a->pos[0]<<1)+a->len;
00352       if(a->boffs>0&&a->eoffs>0)x+=a->eoffs-a->boffs;
00353       hneighbors[0]=_hclusters+i;
00354       nhneighbors=1;
00355       j=nvneighbors>>1;
00356       b=vneighbors[j]->lines[vneighbors[j]->nlines>>1];
00357       for(j=i+1;j<_nhclusters;j++)if(!hmark[j]){
00358         a=_hclusters[j].lines[_hclusters[j].nlines>>1];
00359         if(qr_finder_lines_are_crossing(a,b)){
00360           hmark[j]=1;
00361           x+=(a->pos[0]<<1)+a->len;
00362           if(a->boffs>0&&a->eoffs>0)x+=a->eoffs-a->boffs;
00363           hneighbors[nhneighbors++]=_hclusters+j;
00364         }
00365       }
00366       c=_centers+ncenters++;
00367       c->pos[0]=(x+nhneighbors)/(nhneighbors<<1);
00368       c->pos[1]=(y+nvneighbors)/(nvneighbors<<1);
00369       c->edge_pts=_edge_pts;
00370       nedge_pts=qr_finder_edge_pts_fill(_edge_pts,0,
00371        hneighbors,nhneighbors,0);
00372       nedge_pts=qr_finder_edge_pts_fill(_edge_pts,nedge_pts,
00373        vneighbors,nvneighbors,1);
00374       c->nedge_pts=nedge_pts;
00375       _edge_pts+=nedge_pts;
00376     }
00377   }
00378   free(vmark);
00379   free(hmark);
00380   free(vneighbors);
00381   free(hneighbors);
00382   /*Sort the centers by decreasing numbers of edge points.*/
00383   qsort(_centers,ncenters,sizeof(*_centers),qr_finder_center_cmp);
00384   return ncenters;
00385 }
00386 
00387 /*Locates a set of putative finder centers in the image.
00388   First we search for horizontal and vertical lines that have
00389    (dark:light:dark:light:dark) runs with size ratios of roughly (1:1:3:1:1).
00390   Then we cluster them into groups such that each subsequent pair of endpoints
00391    is close to the line before it in the cluster.
00392   This will locate many line clusters that don't cross a finder pattern, but
00393    qr_finder_find_crossings() will filter most of them out.
00394   Where horizontal and vertical clusters cross, a prospective finder center is
00395    returned.
00396   _centers:  Returns a pointer to a freshly-allocated list of finder centers.
00397              This must be freed by the caller.
00398   _edge_pts: Returns a pointer to a freshly-allocated list of edge points
00399               around those centers.
00400              This must be freed by the caller.
00401   _img:      The binary image to search.
00402   _width:    The width of the image.
00403   _height:   The height of the image.
00404   Return: The number of putative finder centers located.*/
00405 static int qr_finder_centers_locate(qr_finder_center **_centers,
00406  qr_finder_edge_pt **_edge_pts, qr_reader *reader,
00407  int _width,int _height){
00408   qr_finder_line     *hlines = reader->finder_lines[0].lines;
00409   int                 nhlines = reader->finder_lines[0].nlines;
00410   qr_finder_line     *vlines = reader->finder_lines[1].lines;
00411   int                 nvlines = reader->finder_lines[1].nlines;
00412 
00413   qr_finder_line    **hneighbors;
00414   qr_finder_cluster  *hclusters;
00415   int                 nhclusters;
00416   qr_finder_line    **vneighbors;
00417   qr_finder_cluster  *vclusters;
00418   int                 nvclusters;
00419   int                 ncenters;
00420 
00421   /*Cluster the detected lines.*/
00422   hneighbors=(qr_finder_line **)malloc(nhlines*sizeof(*hneighbors));
00423   /*We require more than one line per cluster, so there are at most nhlines/2.*/
00424   hclusters=(qr_finder_cluster *)malloc((nhlines>>1)*sizeof(*hclusters));
00425   nhclusters=qr_finder_cluster_lines(hclusters,hneighbors,hlines,nhlines,0);
00426   /*We need vertical lines to be sorted by X coordinate, with ties broken by Y
00427      coordinate, for clustering purposes.
00428     We scan the image in the opposite order for cache efficiency, so sort the
00429      lines we found here.*/
00430   qsort(vlines,nvlines,sizeof(*vlines),qr_finder_vline_cmp);
00431   vneighbors=(qr_finder_line **)malloc(nvlines*sizeof(*vneighbors));
00432   /*We require more than one line per cluster, so there are at most nvlines/2.*/
00433   vclusters=(qr_finder_cluster *)malloc((nvlines>>1)*sizeof(*vclusters));
00434   nvclusters=qr_finder_cluster_lines(vclusters,vneighbors,vlines,nvlines,1);
00435   /*Find line crossings among the clusters.*/
00436   if(nhclusters>=3&&nvclusters>=3){
00437     qr_finder_edge_pt  *edge_pts;
00438     qr_finder_center   *centers;
00439     int                 nedge_pts;
00440     int                 i;
00441     nedge_pts=0;
00442     for(i=0;i<nhclusters;i++)nedge_pts+=hclusters[i].nlines;
00443     for(i=0;i<nvclusters;i++)nedge_pts+=vclusters[i].nlines;
00444     nedge_pts<<=1;
00445     edge_pts=(qr_finder_edge_pt *)malloc(nedge_pts*sizeof(*edge_pts));
00446     centers=(qr_finder_center *)malloc(
00447      QR_MINI(nhclusters,nvclusters)*sizeof(*centers));
00448     ncenters=qr_finder_find_crossings(centers,edge_pts,
00449      hclusters,nhclusters,vclusters,nvclusters);
00450     *_centers=centers;
00451     *_edge_pts=edge_pts;
00452   }
00453   else ncenters=0;
00454   free(vclusters);
00455   free(vneighbors);
00456   free(hclusters);
00457   free(hneighbors);
00458   return ncenters;
00459 }
00460 
00461 
00462 
00463 static void qr_point_translate(qr_point _point,int _dx,int _dy){
00464   _point[0]+=_dx;
00465   _point[1]+=_dy;
00466 }
00467 
00468 static unsigned qr_point_distance2(const qr_point _p1,const qr_point _p2){
00469   return (_p1[0]-_p2[0])*(_p1[0]-_p2[0])+(_p1[1]-_p2[1])*(_p1[1]-_p2[1]);
00470 }
00471 
00472 /*Returns the cross product of the three points, which is positive if they are
00473    in CCW order (in a right-handed coordinate system), and 0 if they're
00474    colinear.*/
00475 static int qr_point_ccw(const qr_point _p0,
00476  const qr_point _p1,const qr_point _p2){
00477   return (_p1[0]-_p0[0])*(_p2[1]-_p0[1])-(_p1[1]-_p0[1])*(_p2[0]-_p0[0]);
00478 }
00479 
00480 
00481 
00482 /*Evaluates a line equation at a point.
00483   _line: The line to evaluate.
00484   _x:    The X coordinate of the point.
00485   _y:    The y coordinate of the point.
00486   Return: The value of the line equation _line[0]*_x+_line[1]*_y+_line[2].*/
00487 static int qr_line_eval(qr_line _line,int _x,int _y){
00488   return _line[0]*_x+_line[1]*_y+_line[2];
00489 }
00490 
00491 /*Computes a line passing through the given point using the specified second
00492    order statistics.
00493   Given a line defined by the equation
00494     A*x+B*y+C = 0 ,
00495    the least squares fit to n points (x_i,y_i) must satisfy the two equations
00496     A^2 + (Syy - Sxx)/Sxy*A*B - B^2 = 0 ,
00497     C = -(xbar*A+ybar*B) ,
00498    where
00499     xbar = sum(x_i)/n ,
00500     ybar = sum(y_i)/n ,
00501     Sxx = sum((x_i-xbar)**2) ,
00502     Sxy = sum((x_i-xbar)*(y_i-ybar)) ,
00503     Syy = sum((y_i-ybar)**2) .
00504   The quadratic can be solved for the ratio (A/B) or (B/A):
00505     A/B = (Syy + sqrt((Sxx-Syy)**2 + 4*Sxy**2) - Sxx)/(-2*Sxy) ,
00506     B/A = (Sxx + sqrt((Sxx-Syy)**2 + 4*Sxy**2) - Syy)/(-2*Sxy) .
00507   We pick the one that leads to the larger ratio to avoid destructive
00508    cancellation (and e.g., 0/0 for horizontal or vertical lines).
00509   The above solutions correspond to the actual minimum.
00510   The other solution of the quadratic corresponds to a saddle point of the
00511    least squares objective function.
00512   _l:   Returns the fitted line values A, B, and C.
00513   _x0:  The X coordinate of the point the line is supposed to pass through.
00514   _y0:  The Y coordinate of the point the line is supposed to pass through.
00515   _sxx: The sum Sxx.
00516   _sxy: The sum Sxy.
00517   _syy: The sum Syy.
00518   _res: The maximum number of bits occupied by the product of any two of
00519          _l[0] or _l[1].
00520         Smaller numbers give less angular resolution, but allow more overhead
00521          room for computations.*/
00522 static void qr_line_fit(qr_line _l,int _x0,int _y0,
00523  int _sxx,int _sxy,int _syy,int _res){
00524   int dshift;
00525   int dround;
00526   int u;
00527   int v;
00528   int w;
00529   u=abs(_sxx-_syy);
00530   v=-_sxy<<1;
00531   w=qr_ihypot(u,v);
00532   /*Computations in later stages can easily overflow with moderate sizes, so we
00533      compute a shift factor to scale things down into a managable range.
00534     We ensure that the product of any two of _l[0] and _l[1] fits within _res
00535      bits, which allows computation of line intersections without overflow.*/
00536   dshift=QR_MAXI(0,QR_MAXI(qr_ilog(u),qr_ilog(abs(v)))+1-(_res+1>>1));
00537   dround=(1<<dshift)>>1;
00538   if(_sxx>_syy){
00539     _l[0]=v+dround>>dshift;
00540     _l[1]=u+w+dround>>dshift;
00541   }
00542   else{
00543     _l[0]=u+w+dround>>dshift;
00544     _l[1]=v+dround>>dshift;
00545   }
00546   _l[2]=-(_x0*_l[0]+_y0*_l[1]);
00547 }
00548 
00549 /*Perform a least-squares line fit to a list of points.
00550   At least two points are required.*/
00551 static void qr_line_fit_points(qr_line _l,qr_point *_p,int _np,int _res){
00552   int sx;
00553   int sy;
00554   int xmin;
00555   int xmax;
00556   int ymin;
00557   int ymax;
00558   int xbar;
00559   int ybar;
00560   int dx;
00561   int dy;
00562   int sxx;
00563   int sxy;
00564   int syy;
00565   int sshift;
00566   int sround;
00567   int i;
00568   sx=sy=0;
00569   ymax=xmax=INT_MIN;
00570   ymin=xmin=INT_MAX;
00571   for(i=0;i<_np;i++){
00572     sx+=_p[i][0];
00573     xmin=QR_MINI(xmin,_p[i][0]);
00574     xmax=QR_MAXI(xmax,_p[i][0]);
00575     sy+=_p[i][1];
00576     ymin=QR_MINI(ymin,_p[i][1]);
00577     ymax=QR_MAXI(ymax,_p[i][1]);
00578   }
00579   xbar=(sx+(_np>>1))/_np;
00580   ybar=(sy+(_np>>1))/_np;
00581   sshift=QR_MAXI(0,qr_ilog(_np*QR_MAXI(QR_MAXI(xmax-xbar,xbar-xmin),
00582    QR_MAXI(ymax-ybar,ybar-ymin)))-(QR_INT_BITS-1>>1));
00583   sround=(1<<sshift)>>1;
00584   sxx=sxy=syy=0;
00585   for(i=0;i<_np;i++){
00586     dx=_p[i][0]-xbar+sround>>sshift;
00587     dy=_p[i][1]-ybar+sround>>sshift;
00588     sxx+=dx*dx;
00589     sxy+=dx*dy;
00590     syy+=dy*dy;
00591   }
00592   qr_line_fit(_l,xbar,ybar,sxx,sxy,syy,_res);
00593 }
00594 
00595 static void qr_line_orient(qr_line _l,int _x,int _y){
00596   if(qr_line_eval(_l,_x,_y)<0){
00597     _l[0]=-_l[0];
00598     _l[1]=-_l[1];
00599     _l[2]=-_l[2];
00600   }
00601 }
00602 
00603 static int qr_line_isect(qr_point _p,const qr_line _l0,const qr_line _l1){
00604   int d;
00605   int x;
00606   int y;
00607   d=_l0[0]*_l1[1]-_l0[1]*_l1[0];
00608   if(d==0)return -1;
00609   x=_l0[1]*_l1[2]-_l1[1]*_l0[2];
00610   y=_l1[0]*_l0[2]-_l0[0]*_l1[2];
00611   if(d<0){
00612     x=-x;
00613     y=-y;
00614     d=-d;
00615   }
00616   _p[0]=QR_DIVROUND(x,d);
00617   _p[1]=QR_DIVROUND(y,d);
00618   return 0;
00619 }
00620 
00621 
00622 
00623 /*An affine homography.
00624   This maps from the image (at subpel resolution) to a square domain with
00625    power-of-two sides (of res bits) and back.*/
00626 struct qr_aff{
00627   int fwd[2][2];
00628   int inv[2][2];
00629   int x0;
00630   int y0;
00631   int res;
00632 };
00633 
00634 
00635 static void qr_aff_init(qr_aff *_aff,
00636  const qr_point _p0,const qr_point _p1,const qr_point _p2,int _res){
00637   int det;
00638   int dx1;
00639   int dy1;
00640   int dx2;
00641   int dy2;
00642   /*det is ensured to be positive by our caller.*/
00643   det=qr_point_ccw(_p0,_p1,_p2);
00644   dx1=_p1[0]-_p0[0];
00645   dx2=_p2[0]-_p0[0];
00646   dy1=_p1[1]-_p0[1];
00647   dy2=_p2[1]-_p0[1];
00648   _aff->fwd[0][0]=dx1;
00649   _aff->fwd[0][1]=dx2;
00650   _aff->fwd[1][0]=dy1;
00651   _aff->fwd[1][1]=dy2;
00652   _aff->inv[0][0]=QR_DIVROUND(dy2<<_res,det);
00653   _aff->inv[0][1]=QR_DIVROUND(-dx2<<_res,det);
00654   _aff->inv[1][0]=QR_DIVROUND(-dy1<<_res,det);
00655   _aff->inv[1][1]=QR_DIVROUND(dx1<<_res,det);
00656   _aff->x0=_p0[0];
00657   _aff->y0=_p0[1];
00658   _aff->res=_res;
00659 }
00660 
00661 /*Map from the image (at subpel resolution) into the square domain.*/
00662 static void qr_aff_unproject(qr_point _q,const qr_aff *_aff,
00663  int _x,int _y){
00664   _q[0]=_aff->inv[0][0]*(_x-_aff->x0)+_aff->inv[0][1]*(_y-_aff->y0);
00665   _q[1]=_aff->inv[1][0]*(_x-_aff->x0)+_aff->inv[1][1]*(_y-_aff->y0);
00666 }
00667 
00668 /*Map from the square domain into the image (at subpel resolution).*/
00669 static void qr_aff_project(qr_point _p,const qr_aff *_aff,
00670  int _u,int _v){
00671   _p[0]=(_aff->fwd[0][0]*_u+_aff->fwd[0][1]*_v+(1<<_aff->res-1)>>_aff->res)
00672    +_aff->x0;
00673   _p[1]=(_aff->fwd[1][0]*_u+_aff->fwd[1][1]*_v+(1<<_aff->res-1)>>_aff->res)
00674    +_aff->y0;
00675 }
00676 
00677 
00678 
00679 /*A full homography.
00680   Like the affine homography, this maps from the image (at subpel resolution)
00681    to a square domain with power-of-two sides (of res bits) and back.*/
00682 struct qr_hom{
00683   int fwd[3][2];
00684   int inv[3][2];
00685   int fwd22;
00686   int inv22;
00687   int x0;
00688   int y0;
00689   int res;
00690 };
00691 
00692 
00693 static void qr_hom_init(qr_hom *_hom,int _x0,int _y0,
00694  int _x1,int _y1,int _x2,int _y2,int _x3,int _y3,int _res){
00695   int dx10;
00696   int dx20;
00697   int dx30;
00698   int dx31;
00699   int dx32;
00700   int dy10;
00701   int dy20;
00702   int dy30;
00703   int dy31;
00704   int dy32;
00705   int a20;
00706   int a21;
00707   int a22;
00708   int b0;
00709   int b1;
00710   int b2;
00711   int s1;
00712   int s2;
00713   int r1;
00714   int r2;
00715   dx10=_x1-_x0;
00716   dx20=_x2-_x0;
00717   dx30=_x3-_x0;
00718   dx31=_x3-_x1;
00719   dx32=_x3-_x2;
00720   dy10=_y1-_y0;
00721   dy20=_y2-_y0;
00722   dy30=_y3-_y0;
00723   dy31=_y3-_y1;
00724   dy32=_y3-_y2;
00725   a20=dx32*dy10-dx10*dy32;
00726   a21=dx20*dy31-dx31*dy20;
00727   a22=dx32*dy31-dx31*dy32;
00728   /*Figure out if we need to downscale anything.*/
00729   b0=qr_ilog(QR_MAXI(abs(dx10),abs(dy10)))+qr_ilog(abs(a20+a22));
00730   b1=qr_ilog(QR_MAXI(abs(dx20),abs(dy20)))+qr_ilog(abs(a21+a22));
00731   b2=qr_ilog(QR_MAXI(QR_MAXI(abs(a20),abs(a21)),abs(a22)));
00732   s1=QR_MAXI(0,_res+QR_MAXI(QR_MAXI(b0,b1),b2)-(QR_INT_BITS-2));
00733   r1=(1<<s1)>>1;
00734   /*Compute the final coefficients of the forward transform.
00735     The 32x32->64 bit multiplies are really needed for accuracy with large
00736      versions.*/
00737   _hom->fwd[0][0]=QR_FIXMUL(dx10,a20+a22,r1,s1);
00738   _hom->fwd[0][1]=QR_FIXMUL(dx20,a21+a22,r1,s1);
00739   _hom->x0=_x0;
00740   _hom->fwd[1][0]=QR_FIXMUL(dy10,a20+a22,r1,s1);
00741   _hom->fwd[1][1]=QR_FIXMUL(dy20,a21+a22,r1,s1);
00742   _hom->y0=_y0;
00743   _hom->fwd[2][0]=a20+r1>>s1;
00744   _hom->fwd[2][1]=a21+r1>>s1;
00745   _hom->fwd22=s1>_res?a22+(r1>>_res)>>s1-_res:a22<<_res-s1;
00746   /*Now compute the inverse transform.*/
00747   b0=qr_ilog(QR_MAXI(QR_MAXI(abs(dx10),abs(dx20)),abs(dx30)))+
00748    qr_ilog(QR_MAXI(abs(_hom->fwd[0][0]),abs(_hom->fwd[1][0])));
00749   b1=qr_ilog(QR_MAXI(QR_MAXI(abs(dy10),abs(dy20)),abs(dy30)))+
00750    qr_ilog(QR_MAXI(abs(_hom->fwd[0][1]),abs(_hom->fwd[1][1])));
00751   b2=qr_ilog(abs(a22))-s1;
00752   s2=QR_MAXI(0,QR_MAXI(b0,b1)+b2-(QR_INT_BITS-3));
00753   r2=(1<<s2)>>1;
00754   s1+=s2;
00755   r1<<=s2;
00756   /*The 32x32->64 bit multiplies are really needed for accuracy with large
00757      versions.*/
00758   _hom->inv[0][0]=QR_FIXMUL(_hom->fwd[1][1],a22,r1,s1);
00759   _hom->inv[0][1]=QR_FIXMUL(-_hom->fwd[0][1],a22,r1,s1);
00760   _hom->inv[1][0]=QR_FIXMUL(-_hom->fwd[1][0],a22,r1,s1);
00761   _hom->inv[1][1]=QR_FIXMUL(_hom->fwd[0][0],a22,r1,s1);
00762   _hom->inv[2][0]=QR_FIXMUL(_hom->fwd[1][0],_hom->fwd[2][1],
00763    -QR_EXTMUL(_hom->fwd[1][1],_hom->fwd[2][0],r2),s2);
00764   _hom->inv[2][1]=QR_FIXMUL(_hom->fwd[0][1],_hom->fwd[2][0],
00765    -QR_EXTMUL(_hom->fwd[0][0],_hom->fwd[2][1],r2),s2);
00766   _hom->inv22=QR_FIXMUL(_hom->fwd[0][0],_hom->fwd[1][1],
00767    -QR_EXTMUL(_hom->fwd[0][1],_hom->fwd[1][0],r2),s2);
00768   _hom->res=_res;
00769 }
00770 
00771 
00772 /*Map from the image (at subpel resolution) into the square domain.
00773   Returns a negative value if the point went to infinity.*/
00774 static int qr_hom_unproject(qr_point _q,const qr_hom *_hom,int _x,int _y){
00775   int x;
00776   int y;
00777   int w;
00778   _x-=_hom->x0;
00779   _y-=_hom->y0;
00780   x=_hom->inv[0][0]*_x+_hom->inv[0][1]*_y;
00781   y=_hom->inv[1][0]*_x+_hom->inv[1][1]*_y;
00782   w=_hom->inv[2][0]*_x+_hom->inv[2][1]*_y
00783    +_hom->inv22+(1<<_hom->res-1)>>_hom->res;
00784   if(w==0){
00785     _q[0]=x<0?INT_MIN:INT_MAX;
00786     _q[1]=y<0?INT_MIN:INT_MAX;
00787     return -1;
00788   }
00789   else{
00790     if(w<0){
00791       x=-x;
00792       y=-y;
00793       w=-w;
00794     }
00795     _q[0]=QR_DIVROUND(x,w);
00796     _q[1]=QR_DIVROUND(y,w);
00797   }
00798   return 0;
00799 }
00800 
00801 /*Finish a partial projection, converting from homogeneous coordinates to the
00802    normal 2-D representation.
00803   In loops, we can avoid many multiplies by computing the homogeneous _x, _y,
00804    and _w incrementally, but we cannot avoid the divisions, done here.*/
00805 static void qr_hom_fproject(qr_point _p,const qr_hom *_hom,
00806  int _x,int _y,int _w){
00807   if(_w==0){
00808     _p[0]=_x<0?INT_MIN:INT_MAX;
00809     _p[1]=_y<0?INT_MIN:INT_MAX;
00810   }
00811   else{
00812     if(_w<0){
00813       _x=-_x;
00814       _y=-_y;
00815       _w=-_w;
00816     }
00817     _p[0]=QR_DIVROUND(_x,_w)+_hom->x0;
00818     _p[1]=QR_DIVROUND(_y,_w)+_hom->y0;
00819   }
00820 }
00821 
00822 #if defined(QR_DEBUG)
00823 /*Map from the square domain into the image (at subpel resolution).
00824   Currently only used directly by debug code.*/
00825 static void qr_hom_project(qr_point _p,const qr_hom *_hom,
00826  int _u,int _v){
00827   qr_hom_fproject(_p,_hom,
00828    _hom->fwd[0][0]*_u+_hom->fwd[0][1]*_v,
00829    _hom->fwd[1][0]*_u+_hom->fwd[1][1]*_v,
00830    _hom->fwd[2][0]*_u+_hom->fwd[2][1]*_v+_hom->fwd22);
00831 }
00832 #endif
00833 
00834 
00835 
00836 /*All the information we've collected about a finder pattern in the current
00837    configuration.*/
00838 struct qr_finder{
00839   /*The module size along each axis (in the square domain).*/
00840   int                size[2];
00841   /*The version estimated from the module size along each axis.*/
00842   int                eversion[2];
00843   /*The list of classified edge points for each edge.*/
00844   qr_finder_edge_pt *edge_pts[4];
00845   /*The number of edge points classified as belonging to each edge.*/
00846   int                nedge_pts[4];
00847   /*The number of inliers found after running RANSAC on each edge.*/
00848   int                ninliers[4];
00849   /*The center of the finder pattern (in the square domain).*/
00850   qr_point           o;
00851   /*The finder center information from the original image.*/
00852   qr_finder_center  *c;
00853 };
00854 
00855 
00856 static int qr_cmp_edge_pt(const void *_a,const void *_b){
00857   const qr_finder_edge_pt *a;
00858   const qr_finder_edge_pt *b;
00859   a=(const qr_finder_edge_pt *)_a;
00860   b=(const qr_finder_edge_pt *)_b;
00861   return ((a->edge>b->edge)-(a->edge<b->edge)<<1)+
00862    (a->extent>b->extent)-(a->extent<b->extent);
00863 }
00864 
00865 /*Computes the index of the edge each edge point belongs to, and its (signed)
00866    distance along the corresponding axis from the center of the finder pattern
00867    (in the square domain).
00868   The resulting list of edge points is sorted by edge index, with ties broken
00869    by extent.*/
00870 static void qr_finder_edge_pts_aff_classify(qr_finder *_f,const qr_aff *_aff){
00871   qr_finder_center *c;
00872   int               i;
00873   int               e;
00874   c=_f->c;
00875   for(e=0;e<4;e++)_f->nedge_pts[e]=0;
00876   for(i=0;i<c->nedge_pts;i++){
00877     qr_point q;
00878     int      d;
00879     qr_aff_unproject(q,_aff,c->edge_pts[i].pos[0],c->edge_pts[i].pos[1]);
00880     qr_point_translate(q,-_f->o[0],-_f->o[1]);
00881     d=abs(q[1])>abs(q[0]);
00882     e=d<<1|(q[d]>=0);
00883     _f->nedge_pts[e]++;
00884     c->edge_pts[i].edge=e;
00885     c->edge_pts[i].extent=q[d];
00886   }
00887   qsort(c->edge_pts,c->nedge_pts,sizeof(*c->edge_pts),qr_cmp_edge_pt);
00888   _f->edge_pts[0]=c->edge_pts;
00889   for(e=1;e<4;e++)_f->edge_pts[e]=_f->edge_pts[e-1]+_f->nedge_pts[e-1];
00890 }
00891 
00892 /*Computes the index of the edge each edge point belongs to, and its (signed)
00893    distance along the corresponding axis from the center of the finder pattern
00894    (in the square domain).
00895   The resulting list of edge points is sorted by edge index, with ties broken
00896    by extent.*/
00897 static void qr_finder_edge_pts_hom_classify(qr_finder *_f,const qr_hom *_hom){
00898   qr_finder_center *c;
00899   int               i;
00900   int               e;
00901   c=_f->c;
00902   for(e=0;e<4;e++)_f->nedge_pts[e]=0;
00903   for(i=0;i<c->nedge_pts;i++){
00904     qr_point q;
00905     int      d;
00906     if(qr_hom_unproject(q,_hom,
00907      c->edge_pts[i].pos[0],c->edge_pts[i].pos[1])>=0){
00908       qr_point_translate(q,-_f->o[0],-_f->o[1]);
00909       d=abs(q[1])>abs(q[0]);
00910       e=d<<1|(q[d]>=0);
00911       _f->nedge_pts[e]++;
00912       c->edge_pts[i].edge=e;
00913       c->edge_pts[i].extent=q[d];
00914     }
00915     else{
00916       c->edge_pts[i].edge=4;
00917       c->edge_pts[i].extent=q[0];
00918     }
00919   }
00920   qsort(c->edge_pts,c->nedge_pts,sizeof(*c->edge_pts),qr_cmp_edge_pt);
00921   _f->edge_pts[0]=c->edge_pts;
00922   for(e=1;e<4;e++)_f->edge_pts[e]=_f->edge_pts[e-1]+_f->nedge_pts[e-1];
00923 }
00924 
00925 /*TODO: Perhaps these thresholds should be on the module size instead?
00926   Unfortunately, I'd need real-world images of codes with larger versions to
00927    see if these thresholds are still effective, but such versions aren't used
00928    often.*/
00929 
00930 /*The amount that the estimated version numbers are allowed to differ from the
00931    real version number and still be considered valid.*/
00932 #define QR_SMALL_VERSION_SLACK (1)
00933 /*Since cell phone cameras can have severe radial distortion, the estimated
00934    version for larger versions can be off by larger amounts.*/
00935 #define QR_LARGE_VERSION_SLACK (3)
00936 
00937 /*Estimates the size of a module after classifying the edge points.
00938   _width:  The distance between UL and UR in the square domain.
00939   _height: The distance between UL and DL in the square domain.*/
00940 static int qr_finder_estimate_module_size_and_version(qr_finder *_f,
00941  int _width,int _height){
00942   qr_point offs;
00943   int      sums[4];
00944   int      nsums[4];
00945   int      usize;
00946   int      nusize;
00947   int      vsize;
00948   int      nvsize;
00949   int      uversion;
00950   int      vversion;
00951   int      e;
00952   offs[0]=offs[1]=0;
00953   for(e=0;e<4;e++)if(_f->nedge_pts[e]>0){
00954     qr_finder_edge_pt *edge_pts;
00955     int                sum;
00956     int                mean;
00957     int                n;
00958     int                i;
00959     /*Average the samples for this edge, dropping the top and bottom 25%.*/
00960     edge_pts=_f->edge_pts[e];
00961     n=_f->nedge_pts[e];
00962     sum=0;
00963     for(i=(n>>2);i<n-(n>>2);i++)sum+=edge_pts[i].extent;
00964     n=n-((n>>2)<<1);
00965     mean=QR_DIVROUND(sum,n);
00966     offs[e>>1]+=mean;
00967     sums[e]=sum;
00968     nsums[e]=n;
00969   }
00970   else nsums[e]=sums[e]=0;
00971   /*If we have samples on both sides of an axis, refine our idea of where the
00972      unprojected finder center is located.*/
00973   if(_f->nedge_pts[0]>0&&_f->nedge_pts[1]>0){
00974     _f->o[0]-=offs[0]>>1;
00975     sums[0]-=offs[0]*nsums[0]>>1;
00976     sums[1]-=offs[0]*nsums[1]>>1;
00977   }
00978   if(_f->nedge_pts[2]>0&&_f->nedge_pts[3]>0){
00979     _f->o[1]-=offs[1]>>1;
00980     sums[2]-=offs[1]*nsums[2]>>1;
00981     sums[3]-=offs[1]*nsums[3]>>1;
00982   }
00983   /*We must have _some_ samples along each axis... if we don't, our transform
00984      must be pretty severely distorting the original square (e.g., with
00985      coordinates so large as to cause overflow).*/
00986   nusize=nsums[0]+nsums[1];
00987   if(nusize<=0)return -1;
00988   /*The module size is 1/3 the average edge extent.*/
00989   nusize*=3;
00990   usize=sums[1]-sums[0];
00991   usize=((usize<<1)+nusize)/(nusize<<1);
00992   if(usize<=0)return -1;
00993   /*Now estimate the version directly from the module size and the distance
00994      between the finder patterns.
00995     This is done independently using the extents along each axis.
00996     If either falls significantly outside the valid range (1 to 40), reject the
00997      configuration.*/
00998   uversion=(_width-8*usize)/(usize<<2);
00999   if(uversion<1||uversion>40+QR_LARGE_VERSION_SLACK)return -1;
01000   /*Now do the same for the other axis.*/
01001   nvsize=nsums[2]+nsums[3];
01002   if(nvsize<=0)return -1;
01003   nvsize*=3;
01004   vsize=sums[3]-sums[2];
01005   vsize=((vsize<<1)+nvsize)/(nvsize<<1);
01006   if(vsize<=0)return -1;
01007   vversion=(_height-8*vsize)/(vsize<<2);
01008   if(vversion<1||vversion>40+QR_LARGE_VERSION_SLACK)return -1;
01009   /*If the estimated version using extents along one axis is significantly
01010      different than the estimated version along the other axis, then the axes
01011      have significantly different scalings (relative to the grid).
01012     This can happen, e.g., when we have multiple adjacent QR codes, and we've
01013      picked two finder patterns from one and the third finder pattern from
01014      another, e.g.:
01015       X---DL UL---X
01016       |....   |....
01017       X....  UR....
01018     Such a configuration might even pass any other geometric checks if we
01019      didn't reject it here.*/
01020   if(abs(uversion-vversion)>QR_LARGE_VERSION_SLACK)return -1;
01021   _f->size[0]=usize;
01022   _f->size[1]=vsize;
01023   /*We intentionally do not compute an average version from the sizes along
01024      both axes.
01025     In the presence of projective distortion, one of them will be much more
01026      accurate than the other.*/
01027   _f->eversion[0]=uversion;
01028   _f->eversion[1]=vversion;
01029   return 0;
01030 }
01031 
01032 /*Eliminate outliers from the classified edge points with RANSAC.*/
01033 static void qr_finder_ransac(qr_finder *_f,const qr_aff *_hom,
01034  isaac_ctx *_isaac,int _e){
01035   qr_finder_edge_pt *edge_pts;
01036   int                best_ninliers;
01037   int                n;
01038   edge_pts=_f->edge_pts[_e];
01039   n=_f->nedge_pts[_e];
01040   best_ninliers=0;
01041   if(n>1){
01042     int max_iters;
01043     int i;
01044     int j;
01045     /*17 iterations is enough to guarantee an outlier-free sample with more
01046        than 99% probability given as many as 50% outliers.*/
01047     max_iters=17;
01048     for(i=0;i<max_iters;i++){
01049       qr_point  q0;
01050       qr_point  q1;
01051       int       ninliers;
01052       int       thresh;
01053       int       p0i;
01054       int       p1i;
01055       int      *p0;
01056       int      *p1;
01057       int       j;
01058       /*Pick two random points on this edge.*/
01059       p0i=isaac_next_uint(_isaac,n);
01060       p1i=isaac_next_uint(_isaac,n-1);
01061       if(p1i>=p0i)p1i++;
01062       p0=edge_pts[p0i].pos;
01063       p1=edge_pts[p1i].pos;
01064       /*If the corresponding line is not within 45 degrees of the proper
01065          orientation in the square domain, reject it outright.
01066         This can happen, e.g., when highly skewed orientations cause points to
01067          be misclassified into the wrong edge.
01068         The irony is that using such points might produce a line which _does_
01069          pass the corresponding validity checks.*/
01070       qr_aff_unproject(q0,_hom,p0[0],p0[1]);
01071       qr_aff_unproject(q1,_hom,p1[0],p1[1]);
01072       qr_point_translate(q0,-_f->o[0],-_f->o[1]);
01073       qr_point_translate(q1,-_f->o[0],-_f->o[1]);
01074       if(abs(q0[_e>>1]-q1[_e>>1])>abs(q0[1-(_e>>1)]-q1[1-(_e>>1)]))continue;
01075       /*Identify the other edge points which are inliers.
01076         The squared distance should be distributed as a \Chi^2 distribution
01077          with one degree of freedom, which means for a 95% confidence the
01078          point should lie within a factor 3.8414588 ~= 4 times the expected
01079          variance of the point locations.
01080         We grossly approximate the standard deviation as 1 pixel in one
01081          direction, and 0.5 pixels in the other (because we average two
01082          coordinates).*/
01083       thresh=qr_isqrt(qr_point_distance2(p0,p1)<<2*QR_FINDER_SUBPREC+1);
01084       ninliers=0;
01085       for(j=0;j<n;j++){
01086         if(abs(qr_point_ccw(p0,p1,edge_pts[j].pos))<=thresh){
01087           edge_pts[j].extent|=1;
01088           ninliers++;
01089         }
01090         else edge_pts[j].extent&=~1;
01091       }
01092       if(ninliers>best_ninliers){
01093         for(j=0;j<n;j++)edge_pts[j].extent<<=1;
01094         best_ninliers=ninliers;
01095         /*The actual number of iterations required is
01096             log(1-\alpha)/log(1-r*r),
01097            where \alpha is the required probability of taking a sample with
01098             no outliers (e.g., 0.99) and r is the estimated ratio of inliers
01099             (e.g. ninliers/n).
01100           This is just a rough (but conservative) approximation, but it
01101            should be good enough to stop the iteration early when we find
01102            a good set of inliers.*/
01103         if(ninliers>n>>1)max_iters=(67*n-63*ninliers-1)/(n<<1);
01104       }
01105     }
01106     /*Now collect all the inliers at the beginning of the list.*/
01107     for(i=j=0;j<best_ninliers;i++)if(edge_pts[i].extent&2){
01108       if(j<i){
01109         qr_finder_edge_pt tmp;
01110         *&tmp=*(edge_pts+i);
01111         *(edge_pts+j)=*(edge_pts+i);
01112         *(edge_pts+i)=*&tmp;
01113       }
01114       j++;
01115     }
01116   }
01117   _f->ninliers[_e]=best_ninliers;
01118 }
01119 
01120 /*Perform a least-squares line fit to an edge of a finder pattern using the
01121    inliers found by RANSAC.*/
01122 static int qr_line_fit_finder_edge(qr_line _l,
01123  const qr_finder *_f,int _e,int _res){
01124   qr_finder_edge_pt *edge_pts;
01125   qr_point          *pts;
01126   int                npts;
01127   int                i;
01128   npts=_f->ninliers[_e];
01129   if(npts<2)return -1;
01130   /*We could write a custom version of qr_line_fit_points that accesses
01131      edge_pts directly, but this saves on code size and doesn't measurably slow
01132      things down.*/
01133   pts=(qr_point *)malloc(npts*sizeof(*pts));
01134   edge_pts=_f->edge_pts[_e];
01135   for(i=0;i<npts;i++){
01136     pts[i][0]=edge_pts[i].pos[0];
01137     pts[i][1]=edge_pts[i].pos[1];
01138   }
01139   qr_line_fit_points(_l,pts,npts,_res);
01140   /*Make sure the center of the finder pattern lies in the positive halfspace
01141      of the line.*/
01142   qr_line_orient(_l,_f->c->pos[0],_f->c->pos[1]);
01143   free(pts);
01144   return 0;
01145 }
01146 
01147 /*Perform a least-squares line fit to a pair of common finder edges using the
01148    inliers found by RANSAC.
01149   Unlike a normal edge fit, we guarantee that this one succeeds by creating at
01150    least one point on each edge using the estimated module size if it has no
01151    inliers.*/
01152 static void qr_line_fit_finder_pair(qr_line _l,const qr_aff *_aff,
01153  const qr_finder *_f0,const qr_finder *_f1,int _e){
01154   qr_point          *pts;
01155   int                npts;
01156   qr_finder_edge_pt *edge_pts;
01157   qr_point           q;
01158   int                n0;
01159   int                n1;
01160   int                i;
01161   n0=_f0->ninliers[_e];
01162   n1=_f1->ninliers[_e];
01163   /*We could write a custom version of qr_line_fit_points that accesses
01164      edge_pts directly, but this saves on code size and doesn't measurably slow
01165      things down.*/
01166   npts=QR_MAXI(n0,1)+QR_MAXI(n1,1);
01167   pts=(qr_point *)malloc(npts*sizeof(*pts));
01168   if(n0>0){
01169     edge_pts=_f0->edge_pts[_e];
01170     for(i=0;i<n0;i++){
01171       pts[i][0]=edge_pts[i].pos[0];
01172       pts[i][1]=edge_pts[i].pos[1];
01173     }
01174   }
01175   else{
01176     q[0]=_f0->o[0];
01177     q[1]=_f0->o[1];
01178     q[_e>>1]+=_f0->size[_e>>1]*(2*(_e&1)-1);
01179     qr_aff_project(pts[0],_aff,q[0],q[1]);
01180     n0++;
01181   }
01182   if(n1>0){
01183     edge_pts=_f1->edge_pts[_e];
01184     for(i=0;i<n1;i++){
01185       pts[n0+i][0]=edge_pts[i].pos[0];
01186       pts[n0+i][1]=edge_pts[i].pos[1];
01187     }
01188   }
01189   else{
01190     q[0]=_f1->o[0];
01191     q[1]=_f1->o[1];
01192     q[_e>>1]+=_f1->size[_e>>1]*(2*(_e&1)-1);
01193     qr_aff_project(pts[n0],_aff,q[0],q[1]);
01194     n1++;
01195   }
01196   qr_line_fit_points(_l,pts,npts,_aff->res);
01197   /*Make sure at least one finder center lies in the positive halfspace.*/
01198   qr_line_orient(_l,_f0->c->pos[0],_f0->c->pos[1]);
01199   free(pts);
01200 }
01201 
01202 static int qr_finder_quick_crossing_check(const unsigned char *_img,
01203  int _width,int _height,int _x0,int _y0,int _x1,int _y1,int _v){
01204   /*The points must be inside the image, and have a !_v:_v:!_v pattern.
01205     We don't scan the whole line initially, but quickly reject if the endpoints
01206      aren't !_v, or the midpoint isn't _v.
01207     If either end point is out of the image, or we don't encounter a _v pixel,
01208      we return a negative value, indicating the region should be considered
01209      empty.
01210     Otherwise, we return a positive value to indicate it is non-empty.*/
01211   if(_x0<0||_x0>=_width||_y0<0||_y0>=_height||
01212    _x1<0||_x1>=_width||_y1<0||_y1>=_height){
01213     return -1;
01214   }
01215   if(!_img[_y0*_width+_x0]!=_v||!_img[_y1*_width+_x1]!=_v)return 1;
01216   if(!_img[(_y0+_y1>>1)*_width+(_x0+_x1>>1)]==_v)return -1;
01217   return 0;
01218 }
01219 
01220 /*Locate the midpoint of a _v segment along a !_v:_v:!_v line from (_x0,_y0) to
01221    (_x1,_y1).
01222   All coordinates, which are NOT in subpel resolution, must lie inside the
01223    image, and the endpoints are already assumed to have the value !_v.
01224   The returned value is in subpel resolution.*/
01225 static int qr_finder_locate_crossing(const unsigned char *_img,
01226  int _width,int _height,int _x0,int _y0,int _x1,int _y1,int _v,qr_point _p){
01227   qr_point x0;
01228   qr_point x1;
01229   qr_point dx;
01230   int      step[2];
01231   int      steep;
01232   int      err;
01233   int      derr;
01234   /*Use Bresenham's algorithm to trace along the line and find the exact
01235      transitions from !_v to _v and back.*/
01236   x0[0]=_x0;
01237   x0[1]=_y0;
01238   x1[0]=_x1;
01239   x1[1]=_y1;
01240   dx[0]=abs(_x1-_x0);
01241   dx[1]=abs(_y1-_y0);
01242   steep=dx[1]>dx[0];
01243   err=0;
01244   derr=dx[1-steep];
01245   step[0]=((_x0<_x1)<<1)-1;
01246   step[1]=((_y0<_y1)<<1)-1;
01247   /*Find the first crossing from !_v to _v.*/
01248   for(;;){
01249     /*If we make it all the way to the other side, there's no crossing.*/
01250     if(x0[steep]==x1[steep])return -1;
01251     x0[steep]+=step[steep];
01252     err+=derr;
01253     if(err<<1>dx[steep]){
01254       x0[1-steep]+=step[1-steep];
01255       err-=dx[steep];
01256     }
01257     if(!_img[x0[1]*_width+x0[0]]!=_v)break;
01258   }
01259   /*Find the last crossing from _v to !_v.*/
01260   err=0;
01261   for(;;){
01262     if(x0[steep]==x1[steep])break;
01263     x1[steep]-=step[steep];
01264     err+=derr;
01265     if(err<<1>dx[steep]){
01266       x1[1-steep]-=step[1-steep];
01267       err-=dx[steep];
01268     }
01269     if(!_img[x1[1]*_width+x1[0]]!=_v)break;
01270   }
01271   /*Return the midpoint of the _v segment.*/
01272   _p[0]=(x0[0]+x1[0]+1<<QR_FINDER_SUBPREC)>>1;
01273   _p[1]=(x0[1]+x1[1]+1<<QR_FINDER_SUBPREC)>>1;
01274   return 0;
01275 }
01276 
01277 static int qr_aff_line_step(const qr_aff *_aff,qr_line _l,
01278  int _v,int _du,int *_dv){
01279   int shift;
01280   int round;
01281   int dv;
01282   int n;
01283   int d;
01284   n=_aff->fwd[0][_v]*_l[0]+_aff->fwd[1][_v]*_l[1];
01285   d=_aff->fwd[0][1-_v]*_l[0]+_aff->fwd[1][1-_v]*_l[1];
01286   if(d<0){
01287     n=-n;
01288     d=-d;
01289   }
01290   shift=QR_MAXI(0,qr_ilog(_du)+qr_ilog(abs(n))+3-QR_INT_BITS);
01291   round=(1<<shift)>>1;
01292   n=n+round>>shift;
01293   d=d+round>>shift;
01294   /*The line should not be outside 45 degrees of horizontal/vertical.
01295     TODO: We impose this restriction to help ensure the loop below terminates,
01296      but it should not technically be required.
01297     It also, however, ensures we avoid division by zero.*/
01298   if(abs(n)>=d)return -1;
01299   n=-_du*n;
01300   dv=QR_DIVROUND(n,d);
01301   if(abs(dv)>=_du)return -1;
01302   *_dv=dv;
01303   return 0;
01304 }
01305 
01306 /*Computes the Hamming distance between two bit patterns (the number of bits
01307    that differ).
01308   May stop counting after _maxdiff differences.*/
01309 static int qr_hamming_dist(unsigned _y1,unsigned _y2,int _maxdiff){
01310   unsigned y;
01311   int      ret;
01312   y=_y1^_y2;
01313   for(ret=0;ret<_maxdiff&&y;ret++)y&=y-1;
01314   return ret;
01315 }
01316 
01317 /*Retrieve a bit (guaranteed to be 0 or 1) from the image, given coordinates in
01318    subpel resolution which have not been bounds checked.*/
01319 static int qr_img_get_bit(const unsigned char *_img,int _width,int _height,
01320  int _x,int _y){
01321   _x>>=QR_FINDER_SUBPREC;
01322   _y>>=QR_FINDER_SUBPREC;
01323   return _img[QR_CLAMPI(0,_y,_height-1)*_width+QR_CLAMPI(0,_x,_width-1)]!=0;
01324 }
01325 
01326 #if defined(QR_DEBUG)
01327 #include "image.h"
01328 
01329 static void qr_finder_dump_aff_undistorted(qr_finder *_ul,qr_finder *_ur,
01330  qr_finder *_dl,qr_aff *_aff,const unsigned char *_img,int _width,int _height){
01331   unsigned char *gimg;
01332   FILE          *fout;
01333   int            lpsz;
01334   int            pixel_size;
01335   int            dim;
01336   int            min;
01337   int            max;
01338   int            u;
01339   int            y;
01340   int            i;
01341   int            j;
01342   lpsz=qr_ilog(_ur->size[0]+_ur->size[1]+_dl->size[0]+_dl->size[1])-6;
01343   pixel_size=1<<lpsz;
01344   dim=(1<<_aff->res-lpsz)+128;
01345   gimg=(unsigned char *)malloc(dim*dim*sizeof(*gimg));
01346   for(i=0;i<dim;i++)for(j=0;j<dim;j++){
01347     qr_point p;
01348     qr_aff_project(p,_aff,(j-64)<<lpsz,(i-64)<<lpsz);
01349     gimg[i*dim+j]=_img[
01350      QR_CLAMPI(0,p[1]>>QR_FINDER_SUBPREC,_height-1)*_width+
01351      QR_CLAMPI(0,p[0]>>QR_FINDER_SUBPREC,_width-1)];
01352   }
01353   {
01354     min=(_ur->o[0]-7*_ur->size[0]>>lpsz)+64;
01355     if(min<0)min=0;
01356     max=(_ur->o[0]+7*_ur->size[0]>>lpsz)+64;
01357     if(max>dim)max=dim;
01358     for(y=-7;y<=7;y++){
01359       i=(_ur->o[1]+y*_ur->size[1]>>lpsz)+64;
01360       if(i<0||i>=dim)continue;
01361       for(j=min;j<max;j++)gimg[i*dim+j]=0x7F;
01362     }
01363     min=(_ur->o[1]-7*_ur->size[1]>>lpsz)+64;
01364     if(min<0)min=0;
01365     max=(_ur->o[1]+7*_ur->size[1]>>lpsz)+64;
01366     if(max>dim)max=dim;
01367     for(u=-7;u<=7;u++){
01368       j=(_ur->o[0]+u*_ur->size[0]>>lpsz)+64;
01369       if(j<0||j>=dim)continue;
01370       for(i=min;i<max;i++)gimg[i*dim+j]=0x7F;
01371     }
01372   }
01373   {
01374     min=(_dl->o[0]-7*_dl->size[0]>>lpsz)+64;
01375     if(min<0)min=0;
01376     max=(_dl->o[0]+7*_dl->size[0]>>lpsz)+64;
01377     if(max>dim)max=dim;
01378     for(y=-7;y<=7;y++){
01379       i=(_dl->o[1]+y*_dl->size[1]>>lpsz)+64;
01380       if(i<0||i>=dim)continue;
01381       for(j=min;j<max;j++)gimg[i*dim+j]=0x7F;
01382     }
01383     min=(_dl->o[1]-7*_dl->size[1]>>lpsz)+64;
01384     if(min<0)min=0;
01385     max=(_dl->o[1]+7*_dl->size[1]>>lpsz)+64;
01386     if(max>dim)max=dim;
01387     for(u=-7;u<=7;u++){
01388       j=(_dl->o[0]+u*_dl->size[0]>>lpsz)+64;
01389       if(j<0||j>=dim)continue;
01390       for(i=min;i<max;i++)gimg[i*dim+j]=0x7F;
01391     }
01392   }
01393   fout=fopen("undistorted_aff.png","wb");
01394   image_write_png(gimg,dim,dim,fout);
01395   fclose(fout);
01396   free(gimg);
01397 }
01398 
01399 static void qr_finder_dump_hom_undistorted(qr_finder *_ul,qr_finder *_ur,
01400  qr_finder *_dl,qr_hom *_hom,const unsigned char *_img,int _width,int _height){
01401   unsigned char *gimg;
01402   FILE          *fout;
01403   int            lpsz;
01404   int            pixel_size;
01405   int            dim;
01406   int            min;
01407   int            max;
01408   int            u;
01409   int            v;
01410   int            i;
01411   int            j;
01412   lpsz=qr_ilog(_ur->size[0]+_ur->size[1]+_dl->size[0]+_dl->size[1])-6;
01413   pixel_size=1<<lpsz;
01414   dim=(1<<_hom->res-lpsz)+256;
01415   gimg=(unsigned char *)malloc(dim*dim*sizeof(*gimg));
01416   for(i=0;i<dim;i++)for(j=0;j<dim;j++){
01417     qr_point p;
01418     qr_hom_project(p,_hom,(j-128)<<lpsz,(i-128)<<lpsz);
01419     gimg[i*dim+j]=_img[
01420      QR_CLAMPI(0,p[1]>>QR_FINDER_SUBPREC,_height-1)*_width+
01421      QR_CLAMPI(0,p[0]>>QR_FINDER_SUBPREC,_width-1)];
01422   }
01423   {
01424     min=(_ur->o[0]-7*_ur->size[0]>>lpsz)+128;
01425     if(min<0)min=0;
01426     max=(_ur->o[0]+7*_ur->size[0]>>lpsz)+128;
01427     if(max>dim)max=dim;
01428     for(v=-7;v<=7;v++){
01429       i=(_ur->o[1]+v*_ur->size[1]>>lpsz)+128;
01430       if(i<0||i>=dim)continue;
01431       for(j=min;j<max;j++)gimg[i*dim+j]=0x7F;
01432     }
01433     min=(_ur->o[1]-7*_ur->size[1]>>lpsz)+128;
01434     if(min<0)min=0;
01435     max=(_ur->o[1]+7*_ur->size[1]>>lpsz)+128;
01436     if(max>dim)max=dim;
01437     for(u=-7;u<=7;u++){
01438       j=(_ur->o[0]+u*_ur->size[0]>>lpsz)+128;
01439       if(j<0||j>=dim)continue;
01440       for(i=min;i<max;i++)gimg[i*dim+j]=0x7F;
01441     }
01442   }
01443   {
01444     min=(_dl->o[0]-7*_dl->size[0]>>lpsz)+128;
01445     if(min<0)min=0;
01446     max=(_dl->o[0]+7*_dl->size[0]>>lpsz)+128;
01447     if(max>dim)max=dim;
01448     for(v=-7;v<=7;v++){
01449       i=(_dl->o[1]+v*_dl->size[1]>>lpsz)+128;
01450       if(i<0||i>=dim)continue;
01451       for(j=min;j<max;j++)gimg[i*dim+j]=0x7F;
01452     }
01453     min=(_dl->o[1]-7*_dl->size[1]>>lpsz)+128;
01454     if(min<0)min=0;
01455     max=(_dl->o[1]+7*_dl->size[1]>>lpsz)+128;
01456     if(max>dim)max=dim;
01457     for(u=-7;u<=7;u++){
01458       j=(_dl->o[0]+u*_dl->size[0]>>lpsz)+128;
01459       if(j<0||j>=dim)continue;
01460       for(i=min;i<max;i++)gimg[i*dim+j]=0x7F;
01461     }
01462   }
01463   fout=fopen("undistorted_hom.png","wb");
01464   image_write_png(gimg,dim,dim,fout);
01465   fclose(fout);
01466   free(gimg);
01467 }
01468 #endif
01469 
01470 
01471 
01472 /*A homography from one region of the grid back to the image.
01473   Unlike a qr_hom, this does not include an inverse transform and maps directly
01474    from the grid points, not a square with power-of-two sides.*/
01475 struct qr_hom_cell{
01476   int fwd[3][3];
01477   int x0;
01478   int y0;
01479   int u0;
01480   int v0;
01481 };
01482 
01483 
01484 static void qr_hom_cell_init(qr_hom_cell *_cell,int _u0,int _v0,
01485  int _u1,int _v1,int _u2,int _v2,int _u3,int _v3,int _x0,int _y0,
01486  int _x1,int _y1,int _x2,int _y2,int _x3,int _y3){
01487   int du10;
01488   int du20;
01489   int du30;
01490   int du31;
01491   int du32;
01492   int dv10;
01493   int dv20;
01494   int dv30;
01495   int dv31;
01496   int dv32;
01497   int dx10;
01498   int dx20;
01499   int dx30;
01500   int dx31;
01501   int dx32;
01502   int dy10;
01503   int dy20;
01504   int dy30;
01505   int dy31;
01506   int dy32;
01507   int a00;
01508   int a01;
01509   int a02;
01510   int a10;
01511   int a11;
01512   int a12;
01513   int a20;
01514   int a21;
01515   int a22;
01516   int i00;
01517   int i01;
01518   int i10;
01519   int i11;
01520   int i20;
01521   int i21;
01522   int i22;
01523   int b0;
01524   int b1;
01525   int b2;
01526   int shift;
01527   int round;
01528   int x;
01529   int y;
01530   int w;
01531   /*First, correct for the arrangement of the source points.
01532     We take advantage of the fact that we know the source points have a very
01533      limited dynamic range (so there will never be overflow) and a small amount
01534      of projective distortion.*/
01535   du10=_u1-_u0;
01536   du20=_u2-_u0;
01537   du30=_u3-_u0;
01538   du31=_u3-_u1;
01539   du32=_u3-_u2;
01540   dv10=_v1-_v0;
01541   dv20=_v2-_v0;
01542   dv30=_v3-_v0;
01543   dv31=_v3-_v1;
01544   dv32=_v3-_v2;
01545   /*Compute the coefficients of the forward transform from the unit square to
01546      the source point configuration.*/
01547   a20=du32*dv10-du10*dv32;
01548   a21=du20*dv31-du31*dv20;
01549   if(a20||a21)a22=du32*dv31-du31*dv32;
01550   /*If the source grid points aren't in a non-affine arrangement, there's no
01551      reason to scale everything by du32*dv31-du31*dv32.
01552     Not doing so allows a much larger dynamic range, and is the only way we can
01553      initialize a base cell that covers the whole grid.*/
01554   else a22=1;
01555   a00=du10*(a20+a22);
01556   a01=du20*(a21+a22);
01557   a10=dv10*(a20+a22);
01558   a11=dv20*(a21+a22);
01559   /*Now compute the inverse transform.*/
01560   i00=a11*a22;
01561   i01=-a01*a22;
01562   i10=-a10*a22;
01563   i11=a00*a22;
01564   i20=a10*a21-a11*a20;
01565   i21=a01*a20-a00*a21;
01566   i22=a00*a11-a01*a10;
01567   /*Invert the coefficients.
01568     Since i22 is the largest, we divide it by all the others.
01569     The quotient is often exact (e.g., when the source points contain no
01570      projective distortion), and is never zero.
01571     Hence we can use zero to signal "infinity" when the divisor is zero.*/
01572   if(i00)i00=QR_FLIPSIGNI(QR_DIVROUND(i22,abs(i00)),i00);
01573   if(i01)i01=QR_FLIPSIGNI(QR_DIVROUND(i22,abs(i01)),i01);
01574   if(i10)i10=QR_FLIPSIGNI(QR_DIVROUND(i22,abs(i10)),i10);
01575   if(i11)i11=QR_FLIPSIGNI(QR_DIVROUND(i22,abs(i11)),i11);
01576   if(i20)i20=QR_FLIPSIGNI(QR_DIVROUND(i22,abs(i20)),i20);
01577   if(i21)i21=QR_FLIPSIGNI(QR_DIVROUND(i22,abs(i21)),i21);
01578   /*Now compute the map from the unit square into the image.*/
01579   dx10=_x1-_x0;
01580   dx20=_x2-_x0;
01581   dx30=_x3-_x0;
01582   dx31=_x3-_x1;
01583   dx32=_x3-_x2;
01584   dy10=_y1-_y0;
01585   dy20=_y2-_y0;
01586   dy30=_y3-_y0;
01587   dy31=_y3-_y1;
01588   dy32=_y3-_y2;
01589   a20=dx32*dy10-dx10*dy32;
01590   a21=dx20*dy31-dx31*dy20;
01591   a22=dx32*dy31-dx31*dy32;
01592   /*Figure out if we need to downscale anything.*/
01593   b0=qr_ilog(QR_MAXI(abs(dx10),abs(dy10)))+qr_ilog(abs(a20+a22));
01594   b1=qr_ilog(QR_MAXI(abs(dx20),abs(dy20)))+qr_ilog(abs(a21+a22));
01595   b2=qr_ilog(QR_MAXI(QR_MAXI(abs(a20),abs(a21)),abs(a22)));
01596   shift=QR_MAXI(0,QR_MAXI(QR_MAXI(b0,b1),b2)-(QR_INT_BITS-3-QR_ALIGN_SUBPREC));
01597   round=(1<<shift)>>1;
01598   /*Compute the final coefficients of the forward transform.*/
01599   a00=QR_FIXMUL(dx10,a20+a22,round,shift);
01600   a01=QR_FIXMUL(dx20,a21+a22,round,shift);
01601   a10=QR_FIXMUL(dy10,a20+a22,round,shift);
01602   a11=QR_FIXMUL(dy20,a21+a22,round,shift);
01603   /*And compose the two transforms.
01604     Since we inverted the coefficients above, we divide by them here instead
01605      of multiplying.
01606     This lets us take advantage of the full dynamic range.
01607     Note a zero divisor is really "infinity", and thus the quotient should also
01608      be zero.*/
01609   _cell->fwd[0][0]=(i00?QR_DIVROUND(a00,i00):0)+(i10?QR_DIVROUND(a01,i10):0);
01610   _cell->fwd[0][1]=(i01?QR_DIVROUND(a00,i01):0)+(i11?QR_DIVROUND(a01,i11):0);
01611   _cell->fwd[1][0]=(i00?QR_DIVROUND(a10,i00):0)+(i10?QR_DIVROUND(a11,i10):0);
01612   _cell->fwd[1][1]=(i01?QR_DIVROUND(a10,i01):0)+(i11?QR_DIVROUND(a11,i11):0);
01613   _cell->fwd[2][0]=(i00?QR_DIVROUND(a20,i00):0)+(i10?QR_DIVROUND(a21,i10):0)
01614    +(i20?QR_DIVROUND(a22,i20):0)+round>>shift;
01615   _cell->fwd[2][1]=(i01?QR_DIVROUND(a20,i01):0)+(i11?QR_DIVROUND(a21,i11):0)
01616    +(i21?QR_DIVROUND(a22,i21):0)+round>>shift;
01617   _cell->fwd[2][2]=a22+round>>shift;
01618   /*Mathematically, a02 and a12 are exactly zero.
01619     However, that concentrates all of the rounding error in the (_u3,_v3)
01620      corner; we compute offsets which distribute it over the whole range.*/
01621   x=_cell->fwd[0][0]*du10+_cell->fwd[0][1]*dv10;
01622   y=_cell->fwd[1][0]*du10+_cell->fwd[1][1]*dv10;
01623   w=_cell->fwd[2][0]*du10+_cell->fwd[2][1]*dv10+_cell->fwd[2][2];
01624   a02=dx10*w-x;
01625   a12=dy10*w-y;
01626   x=_cell->fwd[0][0]*du20+_cell->fwd[0][1]*dv20;
01627   y=_cell->fwd[1][0]*du20+_cell->fwd[1][1]*dv20;
01628   w=_cell->fwd[2][0]*du20+_cell->fwd[2][1]*dv20+_cell->fwd[2][2];
01629   a02+=dx20*w-x;
01630   a12+=dy20*w-y;
01631   x=_cell->fwd[0][0]*du30+_cell->fwd[0][1]*dv30;
01632   y=_cell->fwd[1][0]*du30+_cell->fwd[1][1]*dv30;
01633   w=_cell->fwd[2][0]*du30+_cell->fwd[2][1]*dv30+_cell->fwd[2][2];
01634   a02+=dx30*w-x;
01635   a12+=dy30*w-y;
01636   _cell->fwd[0][2]=a02+2>>2;
01637   _cell->fwd[1][2]=a12+2>>2;
01638   _cell->x0=_x0;
01639   _cell->y0=_y0;
01640   _cell->u0=_u0;
01641   _cell->v0=_v0;
01642 }
01643 
01644 /*Finish a partial projection, converting from homogeneous coordinates to the
01645    normal 2-D representation.
01646   In loops, we can avoid many multiplies by computing the homogeneous _x, _y,
01647    and _w incrementally, but we cannot avoid the divisions, done here.*/
01648 static void qr_hom_cell_fproject(qr_point _p,const qr_hom_cell *_cell,
01649  int _x,int _y,int _w){
01650   if(_w==0){
01651     _p[0]=_x<0?INT_MIN:INT_MAX;
01652     _p[1]=_y<0?INT_MIN:INT_MAX;
01653   }
01654   else{
01655     if(_w<0){
01656       _x=-_x;
01657       _y=-_y;
01658       _w=-_w;
01659     }
01660     _p[0]=QR_DIVROUND(_x,_w)+_cell->x0;
01661     _p[1]=QR_DIVROUND(_y,_w)+_cell->y0;
01662   }
01663 }
01664 
01665 static void qr_hom_cell_project(qr_point _p,const qr_hom_cell *_cell,
01666  int _u,int _v,int _res){
01667   _u-=_cell->u0<<_res;
01668   _v-=_cell->v0<<_res;
01669   qr_hom_cell_fproject(_p,_cell,
01670    _cell->fwd[0][0]*_u+_cell->fwd[0][1]*_v+(_cell->fwd[0][2]<<_res),
01671    _cell->fwd[1][0]*_u+_cell->fwd[1][1]*_v+(_cell->fwd[1][2]<<_res),
01672    _cell->fwd[2][0]*_u+_cell->fwd[2][1]*_v+(_cell->fwd[2][2]<<_res));
01673 }
01674 
01675 
01676 
01677 /*Retrieves the bits corresponding to the alignment pattern template centered
01678    at the given location in the original image (at subpel precision).*/
01679 static unsigned qr_alignment_pattern_fetch(qr_point _p[5][5],int _x0,int _y0,
01680  const unsigned char *_img,int _width,int _height){
01681   unsigned v;
01682   int      i;
01683   int      j;
01684   int      k;
01685   int      dx;
01686   int      dy;
01687   dx=_x0-_p[2][2][0];
01688   dy=_y0-_p[2][2][1];
01689   v=0;
01690   for(k=i=0;i<5;i++)for(j=0;j<5;j++,k++){
01691     v|=qr_img_get_bit(_img,_width,_height,_p[i][j][0]+dx,_p[i][j][1]+dy)<<k;
01692   }
01693   return v;
01694 }
01695 
01696 /*Searches for an alignment pattern near the given location.*/
01697 static int qr_alignment_pattern_search(qr_point _p,const qr_hom_cell *_cell,
01698  int _u,int _v,int _r,const unsigned char *_img,int _width,int _height){
01699   qr_point c[4];
01700   int      nc[4];
01701   qr_point p[5][5];
01702   qr_point pc;
01703   unsigned best_match;
01704   int      best_dist;
01705   int      bestx;
01706   int      besty;
01707   unsigned match;
01708   int      dist;
01709   int      u;
01710   int      v;
01711   int      x0;
01712   int      y0;
01713   int      w0;
01714   int      x;
01715   int      y;
01716   int      w;
01717   int      dxdu;
01718   int      dydu;
01719   int      dwdu;
01720   int      dxdv;
01721   int      dydv;
01722   int      dwdv;
01723   int      dx;
01724   int      dy;
01725   int      i;
01726   int      j;
01727   /*Build up a basic template using _cell to control shape and scale.
01728     We project the points in the template back to the image just once, since if
01729      the alignment pattern has moved, we don't really know why.
01730     If it's because of radial distortion, or the code wasn't flat, or something
01731      else, there's no reason to expect that a re-projection around each
01732      subsequent search point would be any closer to the actual shape than our
01733      first projection.
01734     Therefore we simply slide this template around, as is.*/
01735   u=(_u-2)-_cell->u0;
01736   v=(_v-2)-_cell->v0;
01737   x0=_cell->fwd[0][0]*u+_cell->fwd[0][1]*v+_cell->fwd[0][2];
01738   y0=_cell->fwd[1][0]*u+_cell->fwd[1][1]*v+_cell->fwd[1][2];
01739   w0=_cell->fwd[2][0]*u+_cell->fwd[2][1]*v+_cell->fwd[2][2];
01740   dxdu=_cell->fwd[0][0];
01741   dydu=_cell->fwd[1][0];
01742   dwdu=_cell->fwd[2][0];
01743   dxdv=_cell->fwd[0][1];
01744   dydv=_cell->fwd[1][1];
01745   dwdv=_cell->fwd[2][1];
01746   for(i=0;i<5;i++){
01747     x=x0;
01748     y=y0;
01749     w=w0;
01750     for(j=0;j<5;j++){
01751       qr_hom_cell_fproject(p[i][j],_cell,x,y,w);
01752       x+=dxdu;
01753       y+=dydu;
01754       w+=dwdu;
01755     }
01756     x0+=dxdv;
01757     y0+=dydv;
01758     w0+=dwdv;
01759   }
01760   bestx=p[2][2][0];
01761   besty=p[2][2][1];
01762   best_match=qr_alignment_pattern_fetch(p,bestx,besty,_img,_width,_height);
01763   best_dist=qr_hamming_dist(best_match,0x1F8D63F,25);
01764   if(best_dist>0){
01765     u=_u-_cell->u0;
01766     v=_v-_cell->v0;
01767     x=_cell->fwd[0][0]*u+_cell->fwd[0][1]*v+_cell->fwd[0][2]<<QR_ALIGN_SUBPREC;
01768     y=_cell->fwd[1][0]*u+_cell->fwd[1][1]*v+_cell->fwd[1][2]<<QR_ALIGN_SUBPREC;
01769     w=_cell->fwd[2][0]*u+_cell->fwd[2][1]*v+_cell->fwd[2][2]<<QR_ALIGN_SUBPREC;
01770     /*Search an area at most _r modules around the target location, in
01771        concentric squares..*/
01772     for(i=1;i<_r<<QR_ALIGN_SUBPREC;i++){
01773       int side_len;
01774       side_len=(i<<1)-1;
01775       x-=dxdu+dxdv;
01776       y-=dydu+dydv;
01777       w-=dwdu+dwdv;
01778       for(j=0;j<4*side_len;j++){
01779         int      dir;
01780         qr_hom_cell_fproject(pc,_cell,x,y,w);
01781         match=qr_alignment_pattern_fetch(p,pc[0],pc[1],_img,_width,_height);
01782         dist=qr_hamming_dist(match,0x1F8D63F,best_dist+1);
01783         if(dist<best_dist){
01784           best_match=match;
01785           best_dist=dist;
01786           bestx=pc[0];
01787           besty=pc[1];
01788         }
01789         if(j<2*side_len){
01790           dir=j>=side_len;
01791           x+=_cell->fwd[0][dir];
01792           y+=_cell->fwd[1][dir];
01793           w+=_cell->fwd[2][dir];
01794         }
01795         else{
01796           dir=j>=3*side_len;
01797           x-=_cell->fwd[0][dir];
01798           y-=_cell->fwd[1][dir];
01799           w-=_cell->fwd[2][dir];
01800         }
01801         if(!best_dist)break;
01802       }
01803       if(!best_dist)break;
01804     }
01805   }
01806   /*If the best result we got was sufficiently bad, reject the match.
01807     If we're wrong and we include it, we can grossly distort the nearby
01808      region, whereas using the initial starting point should at least be
01809      consistent with the geometry we already have.*/
01810   if(best_dist>6){
01811     _p[0]=p[2][2][0];
01812     _p[1]=p[2][2][1];
01813     return -1;
01814   }
01815   /*Now try to get a more accurate location of the pattern center.*/
01816   dx=bestx-p[2][2][0];
01817   dy=besty-p[2][2][1];
01818   memset(nc,0,sizeof(nc));
01819   memset(c,0,sizeof(c));
01820   /*We consider 8 lines across the finder pattern in turn.
01821     If we actually found a symmetric pattern along that line, search for its
01822      exact center in the image.
01823     There are plenty more lines we could use if these don't work, but if we've
01824      found anything remotely close to an alignment pattern, we should be able
01825      to use most of these.*/
01826   for(i=0;i<8;i++){
01827     static const unsigned MASK_TESTS[8][2]={
01828       {0x1040041,0x1000001},{0x0041040,0x0001000},
01829       {0x0110110,0x0100010},{0x0011100,0x0001000},
01830       {0x0420084,0x0400004},{0x0021080,0x0001000},
01831       {0x0006C00,0x0004400},{0x0003800,0x0001000},
01832     };
01833     static const unsigned char MASK_COORDS[8][2]={
01834       {0,0},{1,1},{4,0},{3,1},{2,0},{2,1},{0,2},{1,2}
01835     };
01836     if((best_match&MASK_TESTS[i][0])==MASK_TESTS[i][1]){
01837       int x0;
01838       int y0;
01839       int x1;
01840       int y1;
01841       x0=p[MASK_COORDS[i][1]][MASK_COORDS[i][0]][0]+dx>>QR_FINDER_SUBPREC;
01842       if(x0<0||x0>=_width)continue;
01843       y0=p[MASK_COORDS[i][1]][MASK_COORDS[i][0]][1]+dy>>QR_FINDER_SUBPREC;
01844       if(y0<0||y0>=_height)continue;
01845       x1=p[4-MASK_COORDS[i][1]][4-MASK_COORDS[i][0]][0]+dx>>QR_FINDER_SUBPREC;
01846       if(x1<0||x1>=_width)continue;
01847       y1=p[4-MASK_COORDS[i][1]][4-MASK_COORDS[i][0]][1]+dy>>QR_FINDER_SUBPREC;
01848       if(y1<0||y1>=_height)continue;
01849       if(!qr_finder_locate_crossing(_img,_width,_height,x0,y0,x1,y1,i&1,pc)){
01850         int w;
01851         int cx;
01852         int cy;
01853         cx=pc[0]-bestx;
01854         cy=pc[1]-besty;
01855         if(i&1){
01856           /*Weight crossings around the center dot more highly, as they are
01857              generally more reliable.*/
01858           w=3;
01859           cx+=cx<<1;
01860           cy+=cy<<1;
01861         }
01862         else w=1;
01863         nc[i>>1]+=w;
01864         c[i>>1][0]+=cx;
01865         c[i>>1][1]+=cy;
01866       }
01867     }
01868   }
01869   /*Sum offsets from lines in orthogonal directions.*/
01870   for(i=0;i<2;i++){
01871     int a;
01872     int b;
01873     a=nc[i<<1];
01874     b=nc[i<<1|1];
01875     if(a&&b){
01876       int w;
01877       w=QR_MAXI(a,b);
01878       c[i<<1][0]=QR_DIVROUND(w*(b*c[i<<1][0]+a*c[i<<1|1][0]),a*b);
01879       c[i<<1][1]=QR_DIVROUND(w*(b*c[i<<1][1]+a*c[i<<1|1][1]),a*b);
01880       nc[i<<1]=w<<1;
01881     }
01882     else{
01883       c[i<<1][0]+=c[i<<1|1][0];
01884       c[i<<1][1]+=c[i<<1|1][1];
01885       nc[i<<1]+=b;
01886     }
01887   }
01888   /*Average offsets from pairs of orthogonal lines.*/
01889   c[0][0]+=c[2][0];
01890   c[0][1]+=c[2][1];
01891   nc[0]+=nc[2];
01892   /*If we actually found any such lines, apply the adjustment.*/
01893   if(nc[0]){
01894     dx=QR_DIVROUND(c[0][0],nc[0]);
01895     dy=QR_DIVROUND(c[0][1],nc[0]);
01896     /*But only if it doesn't make things worse.*/
01897     match=qr_alignment_pattern_fetch(p,bestx+dx,besty+dy,_img,_width,_height);
01898     dist=qr_hamming_dist(match,0x1F8D63F,best_dist+1);
01899     if(dist<=best_dist){
01900       bestx+=dx;
01901       besty+=dy;
01902     }
01903   }
01904   _p[0]=bestx;
01905   _p[1]=besty;
01906   return 0;
01907 }
01908 
01909 static int qr_hom_fit(qr_hom *_hom,qr_finder *_ul,qr_finder *_ur,
01910  qr_finder *_dl,qr_point _p[4],const qr_aff *_aff,isaac_ctx *_isaac,
01911  const unsigned char *_img,int _width,int _height){
01912   qr_point *b;
01913   int       nb;
01914   int       cb;
01915   qr_point *r;
01916   int       nr;
01917   int       cr;
01918   qr_line   l[4];
01919   qr_point  q;
01920   qr_point  p;
01921   int       ox;
01922   int       oy;
01923   int       ru;
01924   int       rv;
01925   int       dru;
01926   int       drv;
01927   int       bu;
01928   int       bv;
01929   int       dbu;
01930   int       dbv;
01931   int       rx;
01932   int       ry;
01933   int       drxi;
01934   int       dryi;
01935   int       drxj;
01936   int       dryj;
01937   int       rdone;
01938   int       nrempty;
01939   int       rlastfit;
01940   int       bx;
01941   int       by;
01942   int       dbxi;
01943   int       dbyi;
01944   int       dbxj;
01945   int       dbyj;
01946   int       bdone;
01947   int       nbempty;
01948   int       blastfit;
01949   int       shift;
01950   int       round;
01951   int       version4;
01952   int       brx;
01953   int       bry;
01954   int       i;
01955   /*We attempt to correct large-scale perspective distortion by fitting lines
01956      to the edge of the code area.
01957     We could also look for an alignment pattern now, but that wouldn't work for
01958      version 1 codes, which have no alignment pattern.
01959     Even if the code is supposed to have one, there's go guarantee we'd find it
01960      intact.*/
01961   /*Fitting lines is easy for the edges on which we have two finder patterns.
01962     After the fit, UL is guaranteed to be on the proper side, but if either of
01963      the other two finder patterns aren't, something is wrong.*/
01964   qr_finder_ransac(_ul,_aff,_isaac,0);
01965   qr_finder_ransac(_dl,_aff,_isaac,0);
01966   qr_line_fit_finder_pair(l[0],_aff,_ul,_dl,0);
01967   if(qr_line_eval(l[0],_dl->c->pos[0],_dl->c->pos[1])<0||
01968    qr_line_eval(l[0],_ur->c->pos[0],_ur->c->pos[1])<0){
01969     return -1;
01970   }
01971   qr_finder_ransac(_ul,_aff,_isaac,2);
01972   qr_finder_ransac(_ur,_aff,_isaac,2);
01973   qr_line_fit_finder_pair(l[2],_aff,_ul,_ur,2);
01974   if(qr_line_eval(l[2],_dl->c->pos[0],_dl->c->pos[1])<0||
01975    qr_line_eval(l[2],_ur->c->pos[0],_ur->c->pos[1])<0){
01976     return -1;
01977   }
01978   /*The edges which only have one finder pattern are more difficult.
01979     We start by fitting a line to the edge of the one finder pattern we do
01980      have.
01981     This can fail due to an insufficient number of sample points, and even if
01982      it succeeds can be fairly inaccurate, because all of the points are
01983      clustered in one corner of the QR code.
01984     If it fails, we just use an axis-aligned line in the affine coordinate
01985      system.
01986     Then we walk along the edge of the entire code, looking for
01987      light:dark:light patterns perpendicular to the edge.
01988     Wherever we find one, we take the center of the dark portion as an
01989      additional sample point.
01990     At the end, we re-fit the line using all such sample points found.*/
01991   drv=_ur->size[1]>>1;
01992   qr_finder_ransac(_ur,_aff,_isaac,1);
01993   if(qr_line_fit_finder_edge(l[1],_ur,1,_aff->res)>=0){
01994     if(qr_line_eval(l[1],_ul->c->pos[0],_ul->c->pos[1])<0||
01995      qr_line_eval(l[1],_dl->c->pos[0],_dl->c->pos[1])<0){
01996       return -1;
01997     }
01998     /*Figure out the change in ru for a given change in rv when stepping along
01999        the fitted line.*/
02000     if(qr_aff_line_step(_aff,l[1],1,drv,&dru)<0)return -1;
02001   }
02002   else dru=0;
02003   ru=_ur->o[0]+3*_ur->size[0]-2*dru;
02004   rv=_ur->o[1]-2*drv;
02005   dbu=_dl->size[0]>>1;
02006   qr_finder_ransac(_dl,_aff,_isaac,3);
02007   if(qr_line_fit_finder_edge(l[3],_dl,3,_aff->res)>=0){
02008     if(qr_line_eval(l[3],_ul->c->pos[0],_ul->c->pos[1])<0||
02009      qr_line_eval(l[3],_ur->c->pos[0],_ur->c->pos[1])<0){
02010       return -1;
02011     }
02012     /*Figure out the change in bv for a given change in bu when stepping along
02013        the fitted line.*/
02014     if(qr_aff_line_step(_aff,l[3],0,dbu,&dbv)<0)return -1;
02015   }
02016   else dbv=0;
02017   bu=_dl->o[0]-2*dbu;
02018   bv=_dl->o[1]+3*_dl->size[1]-2*dbv;
02019   /*Set up the initial point lists.*/
02020   nr=rlastfit=_ur->ninliers[1];
02021   cr=nr+(_dl->o[1]-rv+drv-1)/drv;
02022   r=(qr_point *)malloc(cr*sizeof(*r));
02023   for(i=0;i<_ur->ninliers[1];i++){
02024     memcpy(r[i],_ur->edge_pts[1][i].pos,sizeof(r[i]));
02025   }
02026   nb=blastfit=_dl->ninliers[3];
02027   cb=nb+(_ur->o[0]-bu+dbu-1)/dbu;
02028   b=(qr_point *)malloc(cb*sizeof(*b));
02029   for(i=0;i<_dl->ninliers[3];i++){
02030     memcpy(b[i],_dl->edge_pts[3][i].pos,sizeof(b[i]));
02031   }
02032   /*Set up the step parameters for the affine projection.*/
02033   ox=(_aff->x0<<_aff->res)+(1<<_aff->res-1);
02034   oy=(_aff->y0<<_aff->res)+(1<<_aff->res-1);
02035   rx=_aff->fwd[0][0]*ru+_aff->fwd[0][1]*rv+ox;
02036   ry=_aff->fwd[1][0]*ru+_aff->fwd[1][1]*rv+oy;
02037   drxi=_aff->fwd[0][0]*dru+_aff->fwd[0][1]*drv;
02038   dryi=_aff->fwd[1][0]*dru+_aff->fwd[1][1]*drv;
02039   drxj=_aff->fwd[0][0]*_ur->size[0];
02040   dryj=_aff->fwd[1][0]*_ur->size[0];
02041   bx=_aff->fwd[0][0]*bu+_aff->fwd[0][1]*bv+ox;
02042   by=_aff->fwd[1][0]*bu+_aff->fwd[1][1]*bv+oy;
02043   dbxi=_aff->fwd[0][0]*dbu+_aff->fwd[0][1]*dbv;
02044   dbyi=_aff->fwd[1][0]*dbu+_aff->fwd[1][1]*dbv;
02045   dbxj=_aff->fwd[0][1]*_dl->size[1];
02046   dbyj=_aff->fwd[1][1]*_dl->size[1];
02047   /*Now step along the lines, looking for new sample points.*/
02048   nrempty=nbempty=0;
02049   for(;;){
02050     int ret;
02051     int x0;
02052     int y0;
02053     int x1;
02054     int y1;
02055     /*If we take too many steps without encountering a non-zero pixel, assume
02056        we have wandered off the edge and stop looking before we hit the other
02057        side of the quiet region.
02058       Otherwise, stop when the lines cross (if they do so inside the affine
02059        region) or come close to crossing (outside the affine region).
02060       TODO: We don't have any way of detecting when we've wandered into the
02061        code interior; we could stop if the outside sample ever shows up dark,
02062        but this could happen because of noise in the quiet region, too.*/
02063     rdone=rv>=QR_MINI(bv,_dl->o[1]+bv>>1)||nrempty>14;
02064     bdone=bu>=QR_MINI(ru,_ur->o[0]+ru>>1)||nbempty>14;
02065     if(!rdone&&(bdone||rv<bu)){
02066       x0=rx+drxj>>_aff->res+QR_FINDER_SUBPREC;
02067       y0=ry+dryj>>_aff->res+QR_FINDER_SUBPREC;
02068       x1=rx-drxj>>_aff->res+QR_FINDER_SUBPREC;
02069       y1=ry-dryj>>_aff->res+QR_FINDER_SUBPREC;
02070       if(nr>=cr){
02071         cr=cr<<1|1;
02072         r=(qr_point *)realloc(r,cr*sizeof(*r));
02073       }
02074       ret=qr_finder_quick_crossing_check(_img,_width,_height,x0,y0,x1,y1,1);
02075       if(!ret){
02076         ret=qr_finder_locate_crossing(_img,_width,_height,x0,y0,x1,y1,1,r[nr]);
02077       }
02078       if(ret>=0){
02079         if(!ret){
02080           qr_aff_unproject(q,_aff,r[nr][0],r[nr][1]);
02081           /*Move the current point halfway towards the crossing.
02082             We don't move the whole way to give us some robustness to noise.*/
02083           ru=ru+q[0]>>1;
02084           /*But ensure that rv monotonically increases.*/
02085           if(q[1]+drv>rv)rv=rv+q[1]>>1;
02086           rx=_aff->fwd[0][0]*ru+_aff->fwd[0][1]*rv+ox;
02087           ry=_aff->fwd[1][0]*ru+_aff->fwd[1][1]*rv+oy;
02088           nr++;
02089           /*Re-fit the line to update the step direction periodically.*/
02090           if(nr>QR_MAXI(1,rlastfit+(rlastfit>>2))){
02091             qr_line_fit_points(l[1],r,nr,_aff->res);
02092             if(qr_aff_line_step(_aff,l[1],1,drv,&dru)>=0){
02093               drxi=_aff->fwd[0][0]*dru+_aff->fwd[0][1]*drv;
02094               dryi=_aff->fwd[1][0]*dru+_aff->fwd[1][1]*drv;
02095             }
02096             rlastfit=nr;
02097           }
02098         }
02099         else nrempty=0;
02100       }
02101       else nrempty++;
02102       ru+=dru;
02103       /*Our final defense: if we overflow, stop.*/
02104       if(rv+drv>rv)rv+=drv;
02105       else nrempty=INT_MAX;
02106       rx+=drxi;
02107       ry+=dryi;
02108     }
02109     else if(!bdone){
02110       x0=bx+dbxj>>_aff->res+QR_FINDER_SUBPREC;
02111       y0=by+dbyj>>_aff->res+QR_FINDER_SUBPREC;
02112       x1=bx-dbxj>>_aff->res+QR_FINDER_SUBPREC;
02113       y1=by-dbyj>>_aff->res+QR_FINDER_SUBPREC;
02114       if(nb>=cb){
02115         cb=cb<<1|1;
02116         b=(qr_point *)realloc(b,cb*sizeof(*b));
02117       }
02118       ret=qr_finder_quick_crossing_check(_img,_width,_height,x0,y0,x1,y1,1);
02119       if(!ret){
02120         ret=qr_finder_locate_crossing(_img,_width,_height,x0,y0,x1,y1,1,b[nb]);
02121       }
02122       if(ret>=0){
02123         if(!ret){
02124           qr_aff_unproject(q,_aff,b[nb][0],b[nb][1]);
02125           /*Move the current point halfway towards the crossing.
02126             We don't move the whole way to give us some robustness to noise.*/
02127           /*But ensure that bu monotonically increases.*/
02128           if(q[0]+dbu>bu)bu=bu+q[0]>>1;
02129           bv=bv+q[1]>>1;
02130           bx=_aff->fwd[0][0]*bu+_aff->fwd[0][1]*bv+ox;
02131           by=_aff->fwd[1][0]*bu+_aff->fwd[1][1]*bv+oy;
02132           nb++;
02133           /*Re-fit the line to update the step direction periodically.*/
02134           if(nb>QR_MAXI(1,blastfit+(blastfit>>2))){
02135             qr_line_fit_points(l[3],b,nb,_aff->res);
02136             if(qr_aff_line_step(_aff,l[3],0,dbu,&dbv)>=0){
02137               dbxi=_aff->fwd[0][0]*dbu+_aff->fwd[0][1]*dbv;
02138               dbyi=_aff->fwd[1][0]*dbu+_aff->fwd[1][1]*dbv;
02139             }
02140             blastfit=nb;
02141           }
02142         }
02143         nbempty=0;
02144       }
02145       else nbempty++;
02146       /*Our final defense: if we overflow, stop.*/
02147       if(bu+dbu>bu)bu+=dbu;
02148       else nbempty=INT_MAX;
02149       bv+=dbv;
02150       bx+=dbxi;
02151       by+=dbyi;
02152     }
02153     else break;
02154   }
02155   /*Fit the new lines.
02156     If we _still_ don't have enough sample points, then just use an
02157      axis-aligned line from the affine coordinate system (e.g., one parallel
02158      to the opposite edge in the image).*/
02159   if(nr>1)qr_line_fit_points(l[1],r,nr,_aff->res);
02160   else{
02161     qr_aff_project(p,_aff,_ur->o[0]+3*_ur->size[0],_ur->o[1]);
02162     shift=QR_MAXI(0,
02163      qr_ilog(QR_MAXI(abs(_aff->fwd[0][1]),abs(_aff->fwd[1][1])))
02164      -(_aff->res+1>>1));
02165     round=(1<<shift)>>1;
02166     l[1][0]=_aff->fwd[1][1]+round>>shift;
02167     l[1][1]=-_aff->fwd[0][1]+round>>shift;
02168     l[1][2]=-(l[1][0]*p[0]+l[1][1]*p[1]);
02169   }
02170   free(r);
02171   if(nb>1)qr_line_fit_points(l[3],b,nb,_aff->res);
02172   else{
02173     qr_aff_project(p,_aff,_dl->o[0],_dl->o[1]+3*_dl->size[1]);
02174     shift=QR_MAXI(0,
02175      qr_ilog(QR_MAXI(abs(_aff->fwd[0][1]),abs(_aff->fwd[1][1])))
02176      -(_aff->res+1>>1));
02177     round=(1<<shift)>>1;
02178     l[3][0]=_aff->fwd[1][0]+round>>shift;
02179     l[3][1]=-_aff->fwd[0][0]+round>>shift;
02180     l[3][2]=-(l[1][0]*p[0]+l[1][1]*p[1]);
02181   }
02182   free(b);
02183   for(i=0;i<4;i++){
02184     if(qr_line_isect(_p[i],l[i&1],l[2+(i>>1)])<0)return -1;
02185     /*It's plausible for points to be somewhat outside the image, but too far
02186        and too much of the pattern will be gone for it to be decodable.*/
02187     if(_p[i][0]<-_width<<QR_FINDER_SUBPREC||
02188      _p[i][0]>=_width<<QR_FINDER_SUBPREC+1||
02189      _p[i][1]<-_height<<QR_FINDER_SUBPREC||
02190      _p[i][1]>=_height<<QR_FINDER_SUBPREC+1){
02191       return -1;
02192     }
02193   }
02194   /*By default, use the edge intersection point for the bottom-right corner.*/
02195   brx=_p[3][0];
02196   bry=_p[3][1];
02197   /*However, if our average version estimate is greater than 1, NOW we try to
02198      search for an alignment pattern.
02199     We get a much better success rate by doing this after our initial attempt
02200      to promote the transform to a homography than before.
02201     You might also think it would be more reliable to use the interior finder
02202      pattern edges, since the outer ones may be obscured or damaged, and it
02203      would save us a reprojection below, since they would form a nice square
02204      with the location of the alignment pattern, but this turns out to be a bad
02205      idea.
02206     Non-linear distortion is usually maximal on the outside edge, and thus
02207      estimating the grid position from points on the interior means we might
02208      get mis-aligned by the time we reach the edge.*/
02209   version4=_ul->eversion[0]+_ul->eversion[1]+_ur->eversion[0]+_dl->eversion[1];
02210   if(version4>4){
02211     qr_hom_cell cell;
02212     qr_point    p3;
02213     int         dim;
02214     dim=17+version4;
02215     qr_hom_cell_init(&cell,0,0,dim-1,0,0,dim-1,dim-1,dim-1,
02216      _p[0][0],_p[0][1],_p[1][0],_p[1][1],
02217      _p[2][0],_p[2][1],_p[3][0],_p[3][1]);
02218     if(qr_alignment_pattern_search(p3,&cell,dim-7,dim-7,4,
02219      _img,_width,_height)>=0){
02220       int c21;
02221       int dx21;
02222       int dy21;
02223       int mask;
02224       int w;
02225       /*There's no real need to update the bounding box corner, and in fact we
02226          actively perform worse if we do.
02227         Clearly it was good enough for us to find this alignment pattern, so
02228          it should be good enough to use for grid initialization.
02229         The point of doing the search was to get more accurate version
02230          estimates and a better chance of decoding the version and format info.
02231         This is particularly important for small versions that have no encoded
02232          version info, since any mismatch in version renders the code
02233          undecodable.*/
02234       /*We do, however, need four points in a square to initialize our
02235          homography, so project the point from the alignment center to the
02236          corner of the code area.*/
02237       c21=_p[2][0]*_p[1][1]-_p[2][1]*_p[1][0];
02238       dx21=_p[2][0]-_p[1][0];
02239       dy21=_p[2][1]-_p[1][1];
02240       w=(dim-7)*c21
02241        +(dim-13)*(_p[0][0]*dy21-_p[0][1]*dx21)+6*(p3[0]*dy21-p3[1]*dx21);
02242       mask=QR_SIGNMASK(w);
02243       w=abs(w);
02244       brx=(int)QR_DIVROUND(QR_EXTMUL((dim-7)*_p[0][0],p3[0]*dy21,
02245        QR_EXTMUL((dim-13)*p3[0],c21-_p[0][1]*dx21,
02246        QR_EXTMUL(6*_p[0][0],c21-p3[1]*dx21,0)))+mask^mask,w);
02247       bry=(int)QR_DIVROUND(QR_EXTMUL((dim-7)*_p[0][1],-p3[1]*dx21,
02248        QR_EXTMUL((dim-13)*p3[1],c21+_p[0][0]*dy21,
02249        QR_EXTMUL(6*_p[0][1],c21+p3[0]*dy21,0)))+mask^mask,w);
02250     }
02251   }
02252   /*Now we have four points that map to a square: initialize the projection.*/
02253   qr_hom_init(_hom,_p[0][0],_p[0][1],_p[1][0],_p[1][1],
02254    _p[2][0],_p[2][1],brx,bry,QR_HOM_BITS);
02255   return 0;
02256 }
02257 
02258 
02259 
02260 /*The BCH(18,6,3) codes are only used for version information, which must lie
02261    between 7 and 40 (inclusive).*/
02262 static const unsigned BCH18_6_CODES[34]={
02263                                                           0x07C94,
02264   0x085BC,0x09A99,0x0A4D3,0x0BBF6,0x0C762,0x0D847,0x0E60D,0x0F928,
02265   0x10B78,0x1145D,0x12A17,0x13532,0x149A6,0x15683,0x168C9,0x177EC,
02266   0x18EC4,0x191E1,0x1AFAB,0x1B08E,0x1CC1A,0x1D33F,0x1ED75,0x1F250,
02267   0x209D5,0x216F0,0x228BA,0x2379F,0x24B0B,0x2542E,0x26A64,0x27541,
02268   0x28C69
02269 };
02270 
02271 /*Corrects a BCH(18,6,3) code word.
02272   _y: Contains the code word to be checked on input, and the corrected value on
02273        output.
02274   Return: The number of errors.
02275           If more than 3 errors are detected, returns a negative value and
02276            performs no correction.*/
02277 static int bch18_6_correct(unsigned *_y){
02278   unsigned x;
02279   unsigned y;
02280   int      nerrs;
02281   y=*_y;
02282   /*Check the easy case first: see if the data bits were uncorrupted.*/
02283   x=y>>12;
02284   if(x>=7&&x<=40){
02285     nerrs=qr_hamming_dist(y,BCH18_6_CODES[x-7],4);
02286     if(nerrs<4){
02287       *_y=BCH18_6_CODES[x-7];
02288       return nerrs;
02289     }
02290   }
02291   /*Exhaustive search is faster than field operations in GF(19).*/
02292   for(x=0;x<34;x++)if(x+7!=y>>12){
02293     nerrs=qr_hamming_dist(y,BCH18_6_CODES[x],4);
02294     if(nerrs<4){
02295       *_y=BCH18_6_CODES[x];
02296       return nerrs;
02297     }
02298   }
02299   return -1;
02300 }
02301 
02302 #if 0
02303 static unsigned bch18_6_encode(unsigned _x){
02304   return (-(_x&1)&0x01F25)^(-(_x>>1&1)&0x0216F)^(-(_x>>2&1)&0x042DE)^
02305    (-(_x>>3&1)&0x085BC)^(-(_x>>4&1)&0x10B78)^(-(_x>>5&1)&0x209D5);
02306 }
02307 #endif
02308 
02309 /*Reads the version bits near a finder module and decodes the version number.*/
02310 static int qr_finder_version_decode(qr_finder *_f,const qr_hom *_hom,
02311  const unsigned char *_img,int _width,int _height,int _dir){
02312   qr_point q;
02313   unsigned v;
02314   int      x0;
02315   int      y0;
02316   int      w0;
02317   int      dxi;
02318   int      dyi;
02319   int      dwi;
02320   int      dxj;
02321   int      dyj;
02322   int      dwj;
02323   int      ret;
02324   int      i;
02325   int      j;
02326   int      k;
02327   v=0;
02328   q[_dir]=_f->o[_dir]-7*_f->size[_dir];
02329   q[1-_dir]=_f->o[1-_dir]-3*_f->size[1-_dir];
02330   x0=_hom->fwd[0][0]*q[0]+_hom->fwd[0][1]*q[1];
02331   y0=_hom->fwd[1][0]*q[0]+_hom->fwd[1][1]*q[1];
02332   w0=_hom->fwd[2][0]*q[0]+_hom->fwd[2][1]*q[1]+_hom->fwd22;
02333   dxi=_hom->fwd[0][1-_dir]*_f->size[1-_dir];
02334   dyi=_hom->fwd[1][1-_dir]*_f->size[1-_dir];
02335   dwi=_hom->fwd[2][1-_dir]*_f->size[1-_dir];
02336   dxj=_hom->fwd[0][_dir]*_f->size[_dir];
02337   dyj=_hom->fwd[1][_dir]*_f->size[_dir];
02338   dwj=_hom->fwd[2][_dir]*_f->size[_dir];
02339   for(k=i=0;i<6;i++){
02340     int x;
02341     int y;
02342     int w;
02343     x=x0;
02344     y=y0;
02345     w=w0;
02346     for(j=0;j<3;j++,k++){
02347       qr_point p;
02348       qr_hom_fproject(p,_hom,x,y,w);
02349       v|=qr_img_get_bit(_img,_width,_height,p[0],p[1])<<k;
02350       x+=dxj;
02351       y+=dyj;
02352       w+=dwj;
02353     }
02354     x0+=dxi;
02355     y0+=dyi;
02356     w0+=dwi;
02357   }
02358   ret=bch18_6_correct(&v);
02359   /*TODO: I'd certainly hope the order the version bits is accessed in is
02360      well-defined, but I seem to have images for two different codes with the
02361      same version using two different orders?
02362     Maybe the other one is a model 1 code?
02363     Even if I parse the version here, I can't decode the rest of the code.
02364     If this is really needed, we should just re-order the bits.*/
02365 #if 0
02366   if(ret<0){
02367     /*17 16 15 14 13 12 11 10  9  8  7  6  5  4  3  2  1  0
02368        0  3  6  9 12 15  1  4  7 10 13 16  2  5  8 11 14 17
02369       17 13  9  5  1 -3 10  6  2 -2 -6-10  3 -1 -5 -9-13-17*/
02370     v=0;
02371     for(k=i=0;i<3;i++){
02372       p[_dir]=_f->o[_dir]+_f->size[_dir]*(-5-i);
02373       for(j=0;j<6;j++,k++){
02374         qr_point q;
02375         p[1-_dir]=_f->o[1-_dir]+_f->size[1-_dir]*(2-j);
02376         qr_hom_project(q,_hom,p[0],p[1]);
02377         v|=qr_img_get_bit(_img,_width,_height,q[0],q[1])<<k;
02378       }
02379     }
02380     ret=bch18_6_correct(&v);
02381   }
02382 #endif
02383   return ret>=0?(int)(v>>12):ret;
02384 }
02385 
02386 /*Reads the format info bits near the finder modules and decodes them.*/
02387 static int qr_finder_fmt_info_decode(qr_finder *_ul,qr_finder *_ur,
02388  qr_finder *_dl,const qr_hom *_hom,
02389  const unsigned char *_img,int _width,int _height){
02390   qr_point p;
02391   unsigned lo[2];
02392   unsigned hi[2];
02393   int      u;
02394   int      v;
02395   int      x;
02396   int      y;
02397   int      w;
02398   int      dx;
02399   int      dy;
02400   int      dw;
02401   int      fmt_info[4];
02402   int      count[4];
02403   int      nerrs[4];
02404   int      nfmt_info;
02405   int      besti;
02406   int      imax;
02407   int      di;
02408   int      i;
02409   int      k;
02410   /*Read the bits around the UL corner.*/
02411   lo[0]=0;
02412   u=_ul->o[0]+5*_ul->size[0];
02413   v=_ul->o[1]-3*_ul->size[1];
02414   x=_hom->fwd[0][0]*u+_hom->fwd[0][1]*v;
02415   y=_hom->fwd[1][0]*u+_hom->fwd[1][1]*v;
02416   w=_hom->fwd[2][0]*u+_hom->fwd[2][1]*v+_hom->fwd22;
02417   dx=_hom->fwd[0][1]*_ul->size[1];
02418   dy=_hom->fwd[1][1]*_ul->size[1];
02419   dw=_hom->fwd[2][1]*_ul->size[1];
02420   for(k=i=0;;i++){
02421     /*Skip the timing pattern row.*/
02422     if(i!=6){
02423       qr_hom_fproject(p,_hom,x,y,w);
02424       lo[0]|=qr_img_get_bit(_img,_width,_height,p[0],p[1])<<k++;
02425       /*Don't advance q in the last iteration... we'll start the next loop from
02426          the current position.*/
02427       if(i>=8)break;
02428     }
02429     x+=dx;
02430     y+=dy;
02431     w+=dw;
02432   }
02433   hi[0]=0;
02434   dx=-_hom->fwd[0][0]*_ul->size[0];
02435   dy=-_hom->fwd[1][0]*_ul->size[0];
02436   dw=-_hom->fwd[2][0]*_ul->size[0];
02437   while(i-->0){
02438     x+=dx;
02439     y+=dy;
02440     w+=dw;
02441     /*Skip the timing pattern column.*/
02442     if(i!=6){
02443       qr_hom_fproject(p,_hom,x,y,w);
02444       hi[0]|=qr_img_get_bit(_img,_width,_height,p[0],p[1])<<k++;
02445     }
02446   }
02447   /*Read the bits next to the UR corner.*/
02448   lo[1]=0;
02449   u=_ur->o[0]+3*_ur->size[0];
02450   v=_ur->o[1]+5*_ur->size[1];
02451   x=_hom->fwd[0][0]*u+_hom->fwd[0][1]*v;
02452   y=_hom->fwd[1][0]*u+_hom->fwd[1][1]*v;
02453   w=_hom->fwd[2][0]*u+_hom->fwd[2][1]*v+_hom->fwd22;
02454   dx=-_hom->fwd[0][0]*_ur->size[0];
02455   dy=-_hom->fwd[1][0]*_ur->size[0];
02456   dw=-_hom->fwd[2][0]*_ur->size[0];
02457   for(k=0;k<8;k++){
02458     qr_hom_fproject(p,_hom,x,y,w);
02459     lo[1]|=qr_img_get_bit(_img,_width,_height,p[0],p[1])<<k;
02460     x+=dx;
02461     y+=dy;
02462     w+=dw;
02463   }
02464   /*Read the bits next to the DL corner.*/
02465   hi[1]=0;
02466   u=_dl->o[0]+5*_dl->size[0];
02467   v=_dl->o[1]-3*_dl->size[1];
02468   x=_hom->fwd[0][0]*u+_hom->fwd[0][1]*v;
02469   y=_hom->fwd[1][0]*u+_hom->fwd[1][1]*v;
02470   w=_hom->fwd[2][0]*u+_hom->fwd[2][1]*v+_hom->fwd22;
02471   dx=_hom->fwd[0][1]*_dl->size[1];
02472   dy=_hom->fwd[1][1]*_dl->size[1];
02473   dw=_hom->fwd[2][1]*_dl->size[1];
02474   for(k=8;k<15;k++){
02475     qr_hom_fproject(p,_hom,x,y,w);
02476     hi[1]|=qr_img_get_bit(_img,_width,_height,p[0],p[1])<<k;
02477     x+=dx;
02478     y+=dy;
02479     w+=dw;
02480   }
02481   /*For the 8th bit we have 3 samples... use the majority value.
02482     TODO: The DL bit appears to be wrong as much as right? Guess it's not
02483      really a third copy after all, but doesn't appear to be used for data.
02484   i=((lo[0]>>7&1)+(lo[1]>>7&1)+(hi[1]>>7&1)>>1)<<7;
02485   lo[0]=lo[0]&~0x80|i;
02486   lo[1]=lo[1]&~0x80|i;
02487   hi[1]&=~0x80;*/
02488   /*For the remaining bits we have two samples... try them in all
02489      combinations and pick the most popular valid code, breaking ties using
02490      the number of bit errors.*/
02491   imax=2<<(hi[0]!=hi[1]);
02492   di=1+(lo[0]==lo[1]);
02493   nfmt_info=0;
02494   for(i=0;i<imax;i+=di){
02495     unsigned v;
02496     int      ret;
02497     int      j;
02498     v=(lo[i&1]|hi[i>>1])^0x5412;
02499     ret=bch15_5_correct(&v);
02500     v>>=10;
02501     if(ret<0)ret=4;
02502     for(j=0;;j++){
02503       if(j>=nfmt_info){
02504         fmt_info[j]=v;
02505         count[j]=1;
02506         nerrs[j]=ret;
02507         nfmt_info++;
02508         break;
02509       }
02510       if(fmt_info[j]==(int)v){
02511         count[j]++;
02512         if(ret<nerrs[j])nerrs[j]=ret;
02513         break;
02514       }
02515     }
02516   }
02517   besti=0;
02518   for(i=1;i<nfmt_info;i++){
02519     if(nerrs[besti]>3&&nerrs[i]<=3||
02520      count[i]>count[besti]||count[i]==count[besti]&&nerrs[i]<nerrs[besti]){
02521       besti=i;
02522     }
02523   }
02524   return nerrs[besti]<4?fmt_info[besti]:-1;
02525 }
02526 
02527 
02528 
02529 /*The grid used to sample the image bits.
02530   The grid is divided into separate cells bounded by finder patterns and/or
02531    alignment patterns, and a separate map back to the original image is
02532    constructed for each cell.
02533   All of these structural elements, as well as the timing patterns, version
02534    info, and format info, are marked in fpmask so they can easily be skipped
02535    during decode.*/
02536 struct qr_sampling_grid{
02537   qr_hom_cell    *cells[6];
02538   unsigned       *fpmask;
02539   int             cell_limits[6];
02540   int             ncells;
02541 };
02542 
02543 
02544 /*Mark a given region as belonging to the function pattern.*/
02545 static void qr_sampling_grid_fp_mask_rect(qr_sampling_grid *_grid,int _dim,
02546  int _u,int _v,int _w,int _h){
02547   int i;
02548   int j;
02549   int stride;
02550   stride=_dim+QR_INT_BITS-1>>QR_INT_LOGBITS;
02551   /*Note that we store bits column-wise, since that's how they're read out of
02552      the grid.*/
02553   for(j=_u;j<_u+_w;j++)for(i=_v;i<_v+_h;i++){
02554     _grid->fpmask[j*stride+(i>>QR_INT_LOGBITS)]|=1<<(i&QR_INT_BITS-1);
02555   }
02556 }
02557 
02558 /*Determine if a given grid location is inside the function pattern.*/
02559 static int qr_sampling_grid_is_in_fp(const qr_sampling_grid *_grid,int _dim,
02560  int _u,int _v){
02561   return _grid->fpmask[_u*(_dim+QR_INT_BITS-1>>QR_INT_LOGBITS)
02562    +(_v>>QR_INT_LOGBITS)]>>(_v&QR_INT_BITS-1)&1;
02563 }
02564 
02565 /*The spacing between alignment patterns after the second for versions >= 7.
02566   We could compact this more, but the code to access it would eliminate the
02567    gains.*/
02568 static const unsigned char QR_ALIGNMENT_SPACING[34]={
02569   16,18,20,22,24,26,28,
02570   20,22,24,24,26,28,28,
02571   22,24,24,26,26,28,28,
02572   24,24,26,26,26,28,28,
02573   24,26,26,26,28,28
02574 };
02575 
02576 static inline void qr_svg_points(const char *cls,
02577                                  qr_point *p,
02578                                  int n)
02579 {
02580     svg_path_start(cls, 1, 0, 0);
02581     int i;
02582     for(i = 0; i < n; i++, p++)
02583         svg_path_moveto(SVG_ABS, p[0][0], p[0][1]);
02584     svg_path_end();
02585 }
02586 
02587 /*Initialize the sampling grid for each region of the code.
02588   _version:  The (decoded) version number.
02589   _ul_pos:   The location of the UL finder pattern.
02590   _ur_pos:   The location of the UR finder pattern.
02591   _dl_pos:   The location of the DL finder pattern.
02592   _p:        On input, contains estimated positions of the four corner modules.
02593              On output, contains a bounding quadrilateral for the code.
02594   _img:      The binary input image.
02595   _width:    The width of the input image.
02596   _height:   The height of the input image.
02597   Return: 0 on success, or a negative value on error.*/
02598 static void qr_sampling_grid_init(qr_sampling_grid *_grid,int _version,
02599  const qr_point _ul_pos,const qr_point _ur_pos,const qr_point _dl_pos,
02600  qr_point _p[4],const unsigned char *_img,int _width,int _height){
02601   qr_hom_cell          base_cell;
02602   int                  align_pos[7];
02603   int                  dim;
02604   int                  nalign;
02605   int                  i;
02606   dim=17+(_version<<2);
02607   nalign=(_version/7)+2;
02608   /*Create a base cell to bootstrap the alignment pattern search.*/
02609   qr_hom_cell_init(&base_cell,0,0,dim-1,0,0,dim-1,dim-1,dim-1,
02610    _p[0][0],_p[0][1],_p[1][0],_p[1][1],_p[2][0],_p[2][1],_p[3][0],_p[3][1]);
02611   /*Allocate the array of cells.*/
02612   _grid->ncells=nalign-1;
02613   _grid->cells[0]=(qr_hom_cell *)malloc(
02614    (nalign-1)*(nalign-1)*sizeof(*_grid->cells[0]));
02615   for(i=1;i<_grid->ncells;i++)_grid->cells[i]=_grid->cells[i-1]+_grid->ncells;
02616   /*Initialize the function pattern mask.*/
02617   _grid->fpmask=(unsigned *)calloc(dim,
02618    (dim+QR_INT_BITS-1>>QR_INT_LOGBITS)*sizeof(*_grid->fpmask));
02619   /*Mask out the finder patterns (and separators and format info bits).*/
02620   qr_sampling_grid_fp_mask_rect(_grid,dim,0,0,9,9);
02621   qr_sampling_grid_fp_mask_rect(_grid,dim,0,dim-8,9,8);
02622   qr_sampling_grid_fp_mask_rect(_grid,dim,dim-8,0,8,9);
02623   /*Mask out the version number bits.*/
02624   if(_version>6){
02625     qr_sampling_grid_fp_mask_rect(_grid,dim,0,dim-11,6,3);
02626     qr_sampling_grid_fp_mask_rect(_grid,dim,dim-11,0,3,6);
02627   }
02628   /*Mask out the timing patterns.*/
02629   qr_sampling_grid_fp_mask_rect(_grid,dim,9,6,dim-17,1);
02630   qr_sampling_grid_fp_mask_rect(_grid,dim,6,9,1,dim-17);
02631   /*If we have no alignment patterns (e.g., this is a version 1 code), just use
02632      the base cell and hope it's good enough.*/
02633   if(_version<2)memcpy(_grid->cells[0],&base_cell,sizeof(base_cell));
02634   else{
02635     qr_point *q;
02636     qr_point *p;
02637     int       j;
02638     int       k;
02639     q=(qr_point *)malloc(nalign*nalign*sizeof(*q));
02640     p=(qr_point *)malloc(nalign*nalign*sizeof(*p));
02641     /*Initialize the alignment pattern position list.*/
02642     align_pos[0]=6;
02643     align_pos[nalign-1]=dim-7;
02644     if(_version>6){
02645       int d;
02646       d=QR_ALIGNMENT_SPACING[_version-7];
02647       for(i=nalign-1;i-->1;)align_pos[i]=align_pos[i+1]-d;
02648     }
02649     /*Three of the corners use a finder pattern instead of a separate
02650        alignment pattern.*/
02651     q[0][0]=3;
02652     q[0][1]=3;
02653     p[0][0]=_ul_pos[0];
02654     p[0][1]=_ul_pos[1];
02655     q[nalign-1][0]=dim-4;
02656     q[nalign-1][1]=3;
02657     p[nalign-1][0]=_ur_pos[0];
02658     p[nalign-1][1]=_ur_pos[1];
02659     q[(nalign-1)*nalign][0]=3;
02660     q[(nalign-1)*nalign][1]=dim-4;
02661     p[(nalign-1)*nalign][0]=_dl_pos[0];
02662     p[(nalign-1)*nalign][1]=_dl_pos[1];
02663     /*Scan for alignment patterns using a diagonal sweep.*/
02664     for(k=1;k<2*nalign-1;k++){
02665       int jmin;
02666       int jmax;
02667       jmax=QR_MINI(k,nalign-1)-(k==nalign-1);
02668       jmin=QR_MAXI(0,k-(nalign-1))+(k==nalign-1);
02669       for(j=jmin;j<=jmax;j++){
02670         qr_hom_cell *cell;
02671         int          u;
02672         int          v;
02673         int          k;
02674         i=jmax-(j-jmin);
02675         k=i*nalign+j;
02676         u=align_pos[j];
02677         v=align_pos[i];
02678         q[k][0]=u;
02679         q[k][1]=v;
02680         /*Mask out the alignment pattern.*/
02681         qr_sampling_grid_fp_mask_rect(_grid,dim,u-2,v-2,5,5);
02682         /*Pick a cell to use to govern the alignment pattern search.*/
02683         if(i>1&&j>1){
02684           qr_point p0;
02685           qr_point p1;
02686           qr_point p2;
02687           /*Each predictor is basically a straight-line extrapolation from two
02688              neighboring alignment patterns (except possibly near the opposing
02689              finder patterns).*/
02690           qr_hom_cell_project(p0,_grid->cells[i-2]+j-1,u,v,0);
02691           qr_hom_cell_project(p1,_grid->cells[i-2]+j-2,u,v,0);
02692           qr_hom_cell_project(p2,_grid->cells[i-1]+j-2,u,v,0);
02693           /*Take the median of the predictions as the search center.*/
02694           QR_SORT2I(p0[0],p1[0]);
02695           QR_SORT2I(p0[1],p1[1]);
02696           QR_SORT2I(p1[0],p2[0]);
02697           QR_SORT2I(p1[1],p2[1]);
02698           QR_SORT2I(p0[0],p1[0]);
02699           QR_SORT2I(p0[1],p1[1]);
02700           /*We need a cell that has the target point at a known (u,v) location.
02701             Since our cells don't have inverses, just construct one from our
02702              neighboring points.*/
02703           cell=_grid->cells[i-1]+j-1;
02704           qr_hom_cell_init(cell,
02705            q[k-nalign-1][0],q[k-nalign-1][1],q[k-nalign][0],q[k-nalign][1],
02706            q[k-1][0],q[k-1][1],q[k][0],q[k][1],
02707            p[k-nalign-1][0],p[k-nalign-1][1],p[k-nalign][0],p[k-nalign][1],
02708            p[k-1][0],p[k-1][1],p1[0],p1[1]);
02709         }
02710         else if(i>1&&j>0)cell=_grid->cells[i-2]+j-1;
02711         else if(i>0&&j>1)cell=_grid->cells[i-1]+j-2;
02712         else cell=&base_cell;
02713         /*Use a very small search radius.
02714           A large displacement here usually means a false positive (e.g., when
02715            the real alignment pattern is damaged or missing), which can
02716            severely distort the projection.*/
02717         qr_alignment_pattern_search(p[k],cell,u,v,2,_img,_width,_height);
02718         if(i>0&&j>0){
02719           qr_hom_cell_init(_grid->cells[i-1]+j-1,
02720            q[k-nalign-1][0],q[k-nalign-1][1],q[k-nalign][0],q[k-nalign][1],
02721            q[k-1][0],q[k-1][1],q[k][0],q[k][1],
02722            p[k-nalign-1][0],p[k-nalign-1][1],p[k-nalign][0],p[k-nalign][1],
02723            p[k-1][0],p[k-1][1],p[k][0],p[k][1]);
02724         }
02725       }
02726     }
02727     qr_svg_points("align", p, nalign * nalign);
02728     free(q);
02729     free(p);
02730   }
02731   /*Set the limits over which each cell is used.*/
02732   memcpy(_grid->cell_limits,align_pos+1,
02733    (_grid->ncells-1)*sizeof(*_grid->cell_limits));
02734   _grid->cell_limits[_grid->ncells-1]=dim;
02735   /*Produce a bounding square for the code (to mark finder centers with).
02736     Because of non-linear distortion, this might not actually bound the code,
02737      but it should be good enough.
02738     I don't think it's worth computing a convex hull or anything silly like
02739      that.*/
02740   qr_hom_cell_project(_p[0],_grid->cells[0]+0,-1,-1,1);
02741   qr_hom_cell_project(_p[1],_grid->cells[0]+_grid->ncells-1,(dim<<1)-1,-1,1);
02742   qr_hom_cell_project(_p[2],_grid->cells[_grid->ncells-1]+0,-1,(dim<<1)-1,1);
02743   qr_hom_cell_project(_p[3],_grid->cells[_grid->ncells-1]+_grid->ncells-1,
02744    (dim<<1)-1,(dim<<1)-1,1);
02745   /*Clamp the points somewhere near the image (this is really just in case a
02746      corner is near the plane at infinity).*/
02747   for(i=0;i<4;i++){
02748     _p[i][0]=QR_CLAMPI(-_width<<QR_FINDER_SUBPREC,_p[i][0],
02749      _width<<QR_FINDER_SUBPREC+1);
02750     _p[i][1]=QR_CLAMPI(-_height<<QR_FINDER_SUBPREC,_p[i][1],
02751      _height<<QR_FINDER_SUBPREC+1);
02752   }
02753   /*TODO: Make fine adjustments using the timing patterns.
02754     Possible strategy: scan the timing pattern at QR_ALIGN_SUBPREC (or finer)
02755      resolution, use dynamic programming to match midpoints between
02756      transitions to the ideal grid locations.*/
02757 }
02758 
02759 static void qr_sampling_grid_clear(qr_sampling_grid *_grid){
02760   free(_grid->fpmask);
02761   free(_grid->cells[0]);
02762 }
02763 
02764 
02765 
02766 #if defined(QR_DEBUG)
02767 static void qr_sampling_grid_dump(qr_sampling_grid *_grid,int _version,
02768  const unsigned char *_img,int _width,int _height){
02769   unsigned char *gimg;
02770   FILE          *fout;
02771   int            dim;
02772   int            u;
02773   int            v;
02774   int            x;
02775   int            y;
02776   int            w;
02777   int            i;
02778   int            j;
02779   int            r;
02780   int            s;
02781   dim=17+(_version<<2)+8<<QR_ALIGN_SUBPREC;
02782   gimg=(unsigned char *)malloc(dim*dim*sizeof(*gimg));
02783   for(i=0;i<dim;i++)for(j=0;j<dim;j++){
02784     qr_hom_cell *cell;
02785     if(i>=(4<<QR_ALIGN_SUBPREC)&&i<=dim-(5<<QR_ALIGN_SUBPREC)&&
02786      j>=(4<<QR_ALIGN_SUBPREC)&&j<=dim-(5<<QR_ALIGN_SUBPREC)&&
02787      ((!(i&(1<<QR_ALIGN_SUBPREC)-1))^(!(j&(1<<QR_ALIGN_SUBPREC)-1)))){
02788       gimg[i*dim+j]=0x7F;
02789     }
02790     else{
02791       qr_point p;
02792       u=(j>>QR_ALIGN_SUBPREC)-4;
02793       v=(i>>QR_ALIGN_SUBPREC)-4;
02794       for(r=0;r<_grid->ncells-1;r++)if(u<_grid->cell_limits[r])break;
02795       for(s=0;s<_grid->ncells-1;s++)if(v<_grid->cell_limits[s])break;
02796       cell=_grid->cells[s]+r;
02797       u=j-(cell->u0+4<<QR_ALIGN_SUBPREC);
02798       v=i-(cell->v0+4<<QR_ALIGN_SUBPREC);
02799       x=cell->fwd[0][0]*u+cell->fwd[0][1]*v+(cell->fwd[0][2]<<QR_ALIGN_SUBPREC);
02800       y=cell->fwd[1][0]*u+cell->fwd[1][1]*v+(cell->fwd[1][2]<<QR_ALIGN_SUBPREC);
02801       w=cell->fwd[2][0]*u+cell->fwd[2][1]*v+(cell->fwd[2][2]<<QR_ALIGN_SUBPREC);
02802       qr_hom_cell_fproject(p,cell,x,y,w);
02803       gimg[i*dim+j]=_img[
02804        QR_CLAMPI(0,p[1]>>QR_FINDER_SUBPREC,_height-1)*_width+
02805        QR_CLAMPI(0,p[0]>>QR_FINDER_SUBPREC,_width-1)];
02806     }
02807   }
02808   for(v=0;v<17+(_version<<2);v++)for(u=0;u<17+(_version<<2);u++){
02809     if(qr_sampling_grid_is_in_fp(_grid,17+(_version<<2),u,v)){
02810       j=u+4<<QR_ALIGN_SUBPREC;
02811       i=v+4<<QR_ALIGN_SUBPREC;
02812       gimg[(i-1)*dim+j-1]=0x7F;
02813       gimg[(i-1)*dim+j]=0x7F;
02814       gimg[(i-1)*dim+j+1]=0x7F;
02815       gimg[i*dim+j-1]=0x7F;
02816       gimg[i*dim+j+1]=0x7F;
02817       gimg[(i+1)*dim+j-1]=0x7F;
02818       gimg[(i+1)*dim+j]=0x7F;
02819       gimg[(i+1)*dim+j+1]=0x7F;
02820     }
02821   }
02822   fout=fopen("grid.png","wb");
02823   image_write_png(gimg,dim,dim,fout);
02824   fclose(fout);
02825   free(gimg);
02826 }
02827 #endif
02828 
02829 /*Generate the data mask corresponding to the given mask pattern.*/
02830 static void qr_data_mask_fill(unsigned *_mask,int _dim,int _pattern){
02831   int stride;
02832   int i;
02833   int j;
02834   stride=_dim+QR_INT_BITS-1>>QR_INT_LOGBITS;
02835   /*Note that we store bits column-wise, since that's how they're read out of
02836      the grid.*/
02837   switch(_pattern){
02838     /*10101010 i+j+1&1
02839       01010101
02840       10101010
02841       01010101*/
02842     case 0:{
02843       int m;
02844       m=0x55;
02845       for(j=0;j<_dim;j++){
02846         memset(_mask+j*stride,m,stride*sizeof(*_mask));
02847         m^=0xFF;
02848       }
02849     }break;
02850     /*11111111 i+1&1
02851       00000000
02852       11111111
02853       00000000*/
02854     case 1:memset(_mask,0x55,_dim*stride*sizeof(*_mask));break;
02855     /*10010010 (j+1)%3&1
02856       10010010
02857       10010010
02858       10010010*/
02859     case 2:{
02860       unsigned m;
02861       m=0xFF;
02862       for(j=0;j<_dim;j++){
02863         memset(_mask+j*stride,m&0xFF,stride*sizeof(*_mask));
02864         m=m<<8|m>>16;
02865       }
02866     }break;
02867     /*10010010 (i+j+1)%3&1
02868       00100100
02869       01001001
02870       10010010*/
02871     case 3:{
02872       unsigned mi;
02873       unsigned mj;
02874       mj=0;
02875       for(i=0;i<(QR_INT_BITS+2)/3;i++)mj|=1<<3*i;
02876       for(j=0;j<_dim;j++){
02877         mi=mj;
02878         for(i=0;i<stride;i++){
02879           _mask[j*stride+i]=mi;
02880           mi=mi>>QR_INT_BITS%3|mi<<3-QR_INT_BITS%3;
02881         }
02882         mj=mj>>1|mj<<2;
02883       }
02884     }break;
02885     /*11100011 (i>>1)+(j/3)+1&1
02886       11100011
02887       00011100
02888       00011100*/
02889     case 4:{
02890       unsigned m;
02891       m=7;
02892       for(j=0;j<_dim;j++){
02893         memset(_mask+j*stride,(0xCC^-(m&1))&0xFF,stride*sizeof(*_mask));
02894         m=m>>1|m<<5;
02895       }
02896     }break;
02897     /*11111111 !((i*j)%6)
02898       10000010
02899       10010010
02900       10101010*/
02901     case 5:{
02902       for(j=0;j<_dim;j++){
02903         unsigned m;
02904         m=0;
02905         for(i=0;i<6;i++)m|=!((i*j)%6)<<i;
02906         for(i=6;i<QR_INT_BITS;i<<=1)m|=m<<i;
02907         for(i=0;i<stride;i++){
02908           _mask[j*stride+i]=m;
02909           m=m>>QR_INT_BITS%6|m<<6-QR_INT_BITS%6;
02910         }
02911       }
02912     }break;
02913     /*11111111 (i*j)%3+i*j+1&1
02914       11100011
02915       11011011
02916       10101010*/
02917     case 6:{
02918       for(j=0;j<_dim;j++){
02919         unsigned m;
02920         m=0;
02921         for(i=0;i<6;i++)m|=((i*j)%3+i*j+1&1)<<i;
02922         for(i=6;i<QR_INT_BITS;i<<=1)m|=m<<i;
02923         for(i=0;i<stride;i++){
02924           _mask[j*stride+i]=m;
02925           m=m>>QR_INT_BITS%6|m<<6-QR_INT_BITS%6;
02926         }
02927       }
02928     }break;
02929     /*10101010 (i*j)%3+i+j+1&1
02930       00011100
02931       10001110
02932       01010101*/
02933     default:{
02934       for(j=0;j<_dim;j++){
02935         unsigned m;
02936         m=0;
02937         for(i=0;i<6;i++)m|=((i*j)%3+i+j+1&1)<<i;
02938         for(i=6;i<QR_INT_BITS;i<<=1)m|=m<<i;
02939         for(i=0;i<stride;i++){
02940           _mask[j*stride+i]=m;
02941           m=m>>QR_INT_BITS%6|m<<6-QR_INT_BITS%6;
02942         }
02943       }
02944     }break;
02945   }
02946 }
02947 
02948 static void qr_sampling_grid_sample(const qr_sampling_grid *_grid,
02949  unsigned *_data_bits,int _dim,int _fmt_info,
02950  const unsigned char *_img,int _width,int _height){
02951   int stride;
02952   int u0;
02953   int u1;
02954   int j;
02955   /*We initialize the buffer with the data mask and XOR bits into it as we read
02956      them out of the image instead of unmasking in a separate step.*/
02957   qr_data_mask_fill(_data_bits,_dim,_fmt_info&7);
02958   stride=_dim+QR_INT_BITS-1>>QR_INT_LOGBITS;
02959   u0=0;
02960   svg_path_start("sampling-grid", 1, 0, 0);
02961   /*We read data cell-by-cell to avoid having to constantly change which
02962      projection we're using as we read each bit.
02963     This (and the position-dependent data mask) is the reason we buffer the
02964      bits we read instead of converting them directly to codewords here.
02965     Note that bits are stored column-wise, since that's how we'll scan them.*/
02966   for(j=0;j<_grid->ncells;j++){
02967     int i;
02968     int v0;
02969     int v1;
02970     u1=_grid->cell_limits[j];
02971     v0=0;
02972     for(i=0;i<_grid->ncells;i++){
02973       qr_hom_cell *cell;
02974       int          x0;
02975       int          y0;
02976       int          w0;
02977       int          u;
02978       int          du;
02979       int          dv;
02980       v1=_grid->cell_limits[i];
02981       cell=_grid->cells[i]+j;
02982       du=u0-cell->u0;
02983       dv=v0-cell->v0;
02984       x0=cell->fwd[0][0]*du+cell->fwd[0][1]*dv+cell->fwd[0][2];
02985       y0=cell->fwd[1][0]*du+cell->fwd[1][1]*dv+cell->fwd[1][2];
02986       w0=cell->fwd[2][0]*du+cell->fwd[2][1]*dv+cell->fwd[2][2];
02987       for(u=u0;u<u1;u++){
02988         int x;
02989         int y;
02990         int w;
02991         int v;
02992         x=x0;
02993         y=y0;
02994         w=w0;
02995         for(v=v0;v<v1;v++){
02996           /*Skip doing all the divisions and bounds checks if the bit is in the
02997              function pattern.*/
02998           if(!qr_sampling_grid_is_in_fp(_grid,_dim,u,v)){
02999             qr_point p;
03000             qr_hom_cell_fproject(p,cell,x,y,w);
03001             _data_bits[u*stride+(v>>QR_INT_LOGBITS)]^=
03002              qr_img_get_bit(_img,_width,_height,p[0],p[1])<<(v&QR_INT_BITS-1);
03003             svg_path_moveto(SVG_ABS, p[0], p[1]);
03004           }
03005           x+=cell->fwd[0][1];
03006           y+=cell->fwd[1][1];
03007           w+=cell->fwd[2][1];
03008         }
03009         x0+=cell->fwd[0][0];
03010         y0+=cell->fwd[1][0];
03011         w0+=cell->fwd[2][0];
03012       }
03013       v0=v1;
03014     }
03015     u0=u1;
03016   }
03017   svg_path_end();
03018 }
03019 
03020 /*Arranges the sample bits read by qr_sampling_grid_sample() into bytes and
03021    groups those bytes into Reed-Solomon blocks.
03022   The individual block pointers are destroyed by this routine.*/
03023 static void qr_samples_unpack(unsigned char **_blocks,int _nblocks,
03024  int _nshort_data,int _nshort_blocks,const unsigned *_data_bits,
03025  const unsigned *_fp_mask,int _dim){
03026   unsigned bits;
03027   int      biti;
03028   int      stride;
03029   int      blocki;
03030   int      blockj;
03031   int      i;
03032   int      j;
03033   stride=_dim+QR_INT_BITS-1>>QR_INT_LOGBITS;
03034   /*If _all_ the blocks are short, don't skip anything (see below).*/
03035   if(_nshort_blocks>=_nblocks)_nshort_blocks=0;
03036   /*Scan columns in pairs from right to left.*/
03037   bits=0;
03038   for(blocki=blockj=biti=0,j=_dim-1;j>0;j-=2){
03039     unsigned data1;
03040     unsigned data2;
03041     unsigned fp_mask1;
03042     unsigned fp_mask2;
03043     int      nbits;
03044     int      l;
03045     /*Scan up a pair of columns.*/
03046     nbits=(_dim-1&QR_INT_BITS-1)+1;
03047     l=j*stride;
03048     for(i=stride;i-->0;){
03049       data1=_data_bits[l+i];
03050       fp_mask1=_fp_mask[l+i];
03051       data2=_data_bits[l+i-stride];
03052       fp_mask2=_fp_mask[l+i-stride];
03053       while(nbits-->0){
03054         /*Pull a bit from the right column.*/
03055         if(!(fp_mask1>>nbits&1)){
03056           bits=bits<<1|data1>>nbits&1;
03057           biti++;
03058         }
03059         /*Pull a bit from the left column.*/
03060         if(!(fp_mask2>>nbits&1)){
03061           bits=bits<<1|data2>>nbits&1;
03062           biti++;
03063         }
03064         /*If we finished a byte, drop it in a block.*/
03065         if(biti>=8){
03066           biti-=8;
03067           *_blocks[blocki++]++=(unsigned char)(bits>>biti);
03068           /*For whatever reason, the long blocks are at the _end_ of the list,
03069              instead of the beginning.
03070             Even worse, the extra bytes they get come at the end of the data
03071              bytes, before the parity bytes.
03072             Hence the logic here: when we've filled up the data portion of the
03073              short blocks, skip directly to the long blocks for the next byte.
03074             It's also the reason we increment _blocks[blocki] on each store,
03075              instead of just indexing with blockj (after this iteration the
03076              number of bytes in each block differs).*/
03077           if(blocki>=_nblocks)blocki=++blockj==_nshort_data?_nshort_blocks:0;
03078         }
03079       }
03080       nbits=QR_INT_BITS;
03081     }
03082     j-=2;
03083     /*Skip the column with the vertical timing pattern.*/
03084     if(j==6)j--;
03085     /*Scan down a pair of columns.*/
03086     l=j*stride;
03087     for(i=0;i<stride;i++){
03088       data1=_data_bits[l+i];
03089       fp_mask1=_fp_mask[l+i];
03090       data2=_data_bits[l+i-stride];
03091       fp_mask2=_fp_mask[l+i-stride];
03092       nbits=QR_MINI(_dim-(i<<QR_INT_LOGBITS),QR_INT_BITS);
03093       while(nbits-->0){
03094         /*Pull a bit from the right column.*/
03095         if(!(fp_mask1&1)){
03096           bits=bits<<1|data1&1;
03097           biti++;
03098         }
03099         data1>>=1;
03100         fp_mask1>>=1;
03101         /*Pull a bit from the left column.*/
03102         if(!(fp_mask2&1)){
03103           bits=bits<<1|data2&1;
03104           biti++;
03105         }
03106         data2>>=1;
03107         fp_mask2>>=1;
03108         /*If we finished a byte, drop it in a block.*/
03109         if(biti>=8){
03110           biti-=8;
03111           *_blocks[blocki++]++=(unsigned char)(bits>>biti);
03112           /*See comments on the "up" loop for the reason behind this mess.*/
03113           if(blocki>=_nblocks)blocki=++blockj==_nshort_data?_nshort_blocks:0;
03114         }
03115       }
03116     }
03117   }
03118 }
03119 
03120 
03121 /*Bit reading code blatantly stolen^W^Wadapted from libogg/libtheora (because
03122    I've already debugged it and I know it works).
03123   Portions (C) Xiph.Org Foundation 1994-2008, BSD-style license.*/
03124 struct qr_pack_buf{
03125   const unsigned char *buf;
03126   int                  endbyte;
03127   int                  endbit;
03128   int                  storage;
03129 };
03130 
03131 
03132 static void qr_pack_buf_init(qr_pack_buf *_b,
03133  const unsigned char *_data,int _ndata){
03134   _b->buf=_data;
03135   _b->storage=_ndata;
03136   _b->endbyte=_b->endbit=0;
03137 }
03138 
03139 /*Assumes 0<=_bits<=16.*/
03140 static int qr_pack_buf_read(qr_pack_buf *_b,int _bits){
03141   const unsigned char *p;
03142   unsigned             ret;
03143   int                  m;
03144   int                  d;
03145   m=16-_bits;
03146   _bits+=_b->endbit;
03147   d=_b->storage-_b->endbyte;
03148   if(d<=2){
03149     /*Not the main path.*/
03150     if(d*8<_bits){
03151       _b->endbyte+=_bits>>3;
03152       _b->endbit=_bits&7;
03153       return -1;
03154     }
03155     /*Special case to avoid reading p[0] below, which might be past the end of
03156        the buffer; also skips some useless accounting.*/
03157     else if(!_bits)return 0;
03158   }
03159   p=_b->buf+_b->endbyte;
03160   ret=p[0]<<8+_b->endbit;
03161   if(_bits>8){
03162     ret|=p[1]<<_b->endbit;
03163     if(_bits>16)ret|=p[2]>>8-_b->endbit;
03164   }
03165   _b->endbyte+=_bits>>3;
03166   _b->endbit=_bits&7;
03167   return (ret&0xFFFF)>>m;
03168 }
03169 
03170 static int qr_pack_buf_avail(const qr_pack_buf *_b){
03171   return (_b->storage-_b->endbyte<<3)-_b->endbit;
03172 }
03173 
03174 
03175 /*The characters available in QR_MODE_ALNUM.*/
03176 static const unsigned char QR_ALNUM_TABLE[45]={
03177   '0','1','2','3','4','5','6','7','8','9',
03178   'A','B','C','D','E','F','G','H','I','J',
03179   'K','L','M','N','O','P','Q','R','S','T',
03180   'U','V','W','X','Y','Z',' ','$','%','*',
03181   '+','-','.','/',':'
03182 };
03183 
03184 static int qr_code_data_parse(qr_code_data *_qrdata,int _version,
03185  const unsigned char *_data,int _ndata){
03186   qr_pack_buf qpb;
03187   int         centries;
03188   int         len_bits_idx;
03189   /*Entries are stored directly in the struct during parsing.
03190     Caller cleans up any allocated data on failure.*/
03191   _qrdata->entries=NULL;
03192   _qrdata->nentries=0;
03193   _qrdata->sa_size=0;
03194   centries=0;
03195   /*The versions are divided into 3 ranges that each use a different number of
03196      bits for length fields.*/
03197   len_bits_idx=(_version>9)+(_version>26);
03198   qr_pack_buf_init(&qpb,_data,_ndata);
03199   /*While we have enough bits to read a mode...*/
03200   while(qr_pack_buf_avail(&qpb)>=4){
03201     qr_code_data_entry *entry;
03202     int                 mode;
03203     mode=qr_pack_buf_read(&qpb,4);
03204     /*Mode 0 is a terminator.*/
03205     if(!mode)break;
03206     if(_qrdata->nentries>=centries){
03207       centries=centries<<1|1;
03208       _qrdata->entries=(qr_code_data_entry *)realloc(_qrdata->entries,
03209        centries*sizeof(*_qrdata->entries));
03210     }
03211     entry=_qrdata->entries+_qrdata->nentries++;
03212     /*Set the mode to an invalid value until we allocate a buffer for it.
03213       This ensures we don't try to free it on clean-up until then.*/
03214     entry->mode=-1;
03215     switch(mode){
03216       /*The number of bits used to encode the character count for each version
03217          range and each data mode.*/
03218       static const unsigned char LEN_BITS[3][4]={
03219         {10, 9, 8, 8},
03220         {12,11,16,10},
03221         {14,13,16,12}
03222       };
03223       case QR_MODE_NUM:{
03224         unsigned char *buf;
03225         unsigned       bits;
03226         int            len;
03227         int            count;
03228         int            rem;
03229         len=qr_pack_buf_read(&qpb,LEN_BITS[len_bits_idx][0]);
03230         if(len<0)return -1;
03231         /*Check to see if there are enough bits left now, so we don't have to
03232            in the decode loop.*/
03233         count=len/3;
03234         rem=len%3;
03235         if(qr_pack_buf_avail(&qpb)<10*count+7*(rem>>1&1)+4*(rem&1))return -1;
03236         entry->mode=mode;
03237         entry->payload.data.buf=buf=(unsigned char *)malloc(len*sizeof(*buf));
03238         entry->payload.data.len=len;
03239         /*Read groups of 3 digits encoded in 10 bits.*/
03240         while(count-->0){
03241           bits=qr_pack_buf_read(&qpb,10);
03242           if(bits>=1000)return -1;
03243           *buf++=(unsigned char)('0'+bits/100);
03244           bits%=100;
03245           *buf++=(unsigned char)('0'+bits/10);
03246           *buf++=(unsigned char)('0'+bits%10);
03247         }
03248         /*Read the last two digits encoded in 7 bits.*/
03249         if(rem>1){
03250           bits=qr_pack_buf_read(&qpb,7);
03251           if(bits>=100)return -1;
03252           *buf++=(unsigned char)('0'+bits/10);
03253           *buf++=(unsigned char)('0'+bits%10);
03254         }
03255         /*Or the last one digit encoded in 4 bits.*/
03256         else if(rem){
03257           bits=qr_pack_buf_read(&qpb,4);
03258           if(bits>=10)return -1;
03259           *buf++=(unsigned char)('0'+bits);
03260         }
03261       }break;
03262       case QR_MODE_ALNUM:{
03263         unsigned char *buf;
03264         unsigned       bits;
03265         int            len;
03266         int            count;
03267         int            rem;
03268         len=qr_pack_buf_read(&qpb,LEN_BITS[len_bits_idx][1]);
03269         if(len<0)return -1;
03270         /*Check to see if there are enough bits left now, so we don't have to
03271            in the decode loop.*/
03272         count=len>>1;
03273         rem=len&1;
03274         if(qr_pack_buf_avail(&qpb)<11*count+6*rem)return -1;
03275         entry->mode=mode;
03276         entry->payload.data.buf=buf=(unsigned char *)malloc(len*sizeof(*buf));
03277         entry->payload.data.len=len;
03278         /*Read groups of two characters encoded in 11 bits.*/
03279         while(count-->0){
03280           bits=qr_pack_buf_read(&qpb,11);
03281           if(bits>=2025)return -1;
03282           *buf++=QR_ALNUM_TABLE[bits/45];
03283           *buf++=QR_ALNUM_TABLE[bits%45];
03284           len-=2;
03285         }
03286         /*Read the last character encoded in 6 bits.*/
03287         if(rem){
03288           bits=qr_pack_buf_read(&qpb,6);
03289           if(bits>=45)return -1;
03290           *buf++=QR_ALNUM_TABLE[bits];
03291         }
03292       }break;
03293       /*Structured-append header.*/
03294       case QR_MODE_STRUCT:{
03295         int bits;
03296         entry->mode=mode;
03297         bits=qr_pack_buf_read(&qpb,16);
03298         if(bits<0)return -1;
03299         /*We also save a copy of the data in _qrdata for easy reference when
03300            grouping structured-append codes.
03301           If for some reason the code has multiple S-A headers, last one wins
03302            (TODO: should we return an error instead?).*/
03303         _qrdata->sa_index=entry->payload.sa.sa_index=
03304          (unsigned char)(bits>>12&0xF);
03305         _qrdata->sa_size=entry->payload.sa.sa_size=
03306          (unsigned char)((bits>>8&0xF)+1);
03307         _qrdata->sa_parity=entry->payload.sa.sa_parity=
03308          (unsigned char)(bits&0xFF);
03309       }break;
03310       case QR_MODE_BYTE:{
03311         unsigned char *buf;
03312         int            len;
03313         len=qr_pack_buf_read(&qpb,LEN_BITS[len_bits_idx][2]);
03314         if(len<0)return -1;
03315         /*Check to see if there are enough bits left now, so we don't have to
03316            in the decode loop.*/
03317         if(qr_pack_buf_avail(&qpb)<len<<3)return -1;
03318         entry->mode=mode;
03319         entry->payload.data.buf=buf=(unsigned char *)malloc(len*sizeof(*buf));
03320         entry->payload.data.len=len;
03321         while(len-->0)*buf++=(unsigned char)qr_pack_buf_read(&qpb,8);
03322       }break;
03323       /*FNC1 first position marker.*/
03324       case QR_MODE_FNC1_1ST:entry->mode=mode;break;
03325       /*Extended Channel Interpretation data.*/
03326       case QR_MODE_ECI:{
03327         unsigned val;
03328         int      bits;
03329         /*ECI uses a variable-width encoding similar to UTF-8*/
03330         bits=qr_pack_buf_read(&qpb,8);
03331         if(bits<0)return -1;
03332         /*One byte:*/
03333         if(!(bits&0x80))val=bits;
03334         /*Two bytes:*/
03335         else if(!(bits&0x40)){
03336           val=bits&0x3F<<8;
03337           bits=qr_pack_buf_read(&qpb,8);
03338           if(bits<0)return -1;
03339           val|=bits;
03340         }
03341         /*Three bytes:*/
03342         else if(!(bits&0x20)){
03343           val=bits&0x1F<<16;
03344           bits=qr_pack_buf_read(&qpb,16);
03345           if(bits<0)return -1;
03346           val|=bits;
03347           /*Valid ECI values are 0...999999.*/
03348           if(val>=1000000)return -1;
03349         }
03350         /*Invalid lead byte.*/
03351         else return -1;
03352         entry->mode=mode;
03353         entry->payload.eci=val;
03354       }break;
03355       case QR_MODE_KANJI:{
03356         unsigned char *buf;
03357         unsigned       bits;
03358         int            len;
03359         len=qr_pack_buf_read(&qpb,LEN_BITS[len_bits_idx][3]);
03360         if(len<0)return -1;
03361         /*Check to see if there are enough bits left now, so we don't have to
03362            in the decode loop.*/
03363         if(qr_pack_buf_avail(&qpb)<13*len)return -1;
03364         entry->mode=mode;
03365         entry->payload.data.buf=buf=(unsigned char *)malloc(2*len*sizeof(*buf));
03366         entry->payload.data.len=2*len;
03367         /*Decode 2-byte SJIS characters encoded in 13 bits.*/
03368         while(len-->0){
03369           bits=qr_pack_buf_read(&qpb,13);
03370           bits=(bits/0xC0<<8|bits%0xC0)+0x8140;
03371           if(bits>=0xA000)bits+=0x4000;
03372           /*Are values 0xXX7F, 0xXXFD...0xXXFF always invalid?
03373             Should we reject them here?*/
03374           *buf++=(unsigned char)(bits>>8);
03375           *buf++=(unsigned char)(bits&0xFF);
03376         }
03377       }break;
03378       /*FNC1 second position marker.*/
03379       case QR_MODE_FNC1_2ND:entry->mode=mode;break;
03380       /*Unknown mode number:*/
03381       default:{
03382         /*Unfortunately, because we have to understand the format of a mode to
03383            know how many bits it occupies, we can't skip unknown modes.
03384           Therefore we have to fail.*/
03385         return -1;
03386       }break;
03387     }
03388   }
03389   /*TODO: If there was a S-A header, we should compute the parity of this
03390      code; how are non-data modes handled (ECI, FNC1)?*/
03391   _qrdata->self_parity=0;
03392   /*Success.*/
03393   _qrdata->entries=(qr_code_data_entry *)realloc(_qrdata->entries,
03394    _qrdata->nentries*sizeof(*_qrdata->entries));
03395   return 0;
03396 }
03397 
03398 static void qr_code_data_clear(qr_code_data *_qrdata){
03399   int i;
03400   for(i=0;i<_qrdata->nentries;i++){
03401     if(QR_MODE_HAS_DATA(_qrdata->entries[i].mode)){
03402       free(_qrdata->entries[i].payload.data.buf);
03403     }
03404   }
03405   free(_qrdata->entries);
03406 }
03407 
03408 
03409 void qr_code_data_list_init(qr_code_data_list *_qrlist){
03410   _qrlist->qrdata=NULL;
03411   _qrlist->nqrdata=_qrlist->cqrdata=0;
03412 }
03413 
03414 void qr_code_data_list_clear(qr_code_data_list *_qrlist){
03415   int i;
03416   for(i=0;i<_qrlist->nqrdata;i++)qr_code_data_clear(_qrlist->qrdata+i);
03417   free(_qrlist->qrdata);
03418   qr_code_data_list_init(_qrlist);
03419 }
03420 
03421 static void qr_code_data_list_add(qr_code_data_list *_qrlist,
03422  qr_code_data *_qrdata){
03423   if(_qrlist->nqrdata>=_qrlist->cqrdata){
03424     _qrlist->cqrdata=_qrlist->cqrdata<<1|1;
03425     _qrlist->qrdata=(qr_code_data *)realloc(_qrlist->qrdata,
03426      _qrlist->cqrdata*sizeof(*_qrlist->qrdata));
03427   }
03428   memcpy(_qrlist->qrdata+_qrlist->nqrdata++,_qrdata,sizeof(*_qrdata));
03429 }
03430 
03431 #if 0
03432 static const unsigned short QR_NCODEWORDS[40]={
03433     26,  44,  70, 100, 134, 172, 196, 242, 292, 346,
03434    404, 466, 532, 581, 655, 733, 815, 901, 991,1085,
03435   1156,1258,1364,1474,1588,1706,1828,1921,2051,2185,
03436   2323,2465,2611,2761,2876,3034,3196,3362,3532,3706
03437 };
03438 #endif
03439 
03440 /*The total number of codewords in a QR code.*/
03441 static int qr_code_ncodewords(unsigned _version){
03442   unsigned nalign;
03443   /*This is 24-27 instructions on ARM in thumb mode, or a 26-32 byte savings
03444      over just using a table (not counting the instructions that would be
03445      needed to do the table lookup).*/
03446   if(_version==1)return 26;
03447   nalign=(_version/7)+2;
03448   return (_version<<4)*(_version+8)
03449    -(5*nalign)*(5*nalign-2)+36*(_version<7)+83>>3;
03450 }
03451 
03452 #if 0
03453 /*The number of parity bytes per Reed-Solomon block for each version and error
03454    correction level.*/
03455 static const unsigned char QR_RS_NPAR[40][4]={
03456   { 7,10,13,17},{10,16,22,28},{15,26,18,22},{20,18,26,16},
03457   {26,24,18,22},{18,16,24,28},{20,18,18,26},{24,22,22,26},
03458   {30,22,20,24},{18,26,24,28},{20,30,28,24},{24,22,26,28},
03459   {26,22,24,22},{30,24,20,24},{22,24,30,24},{24,28,24,30},
03460   {28,28,28,28},{30,26,28,28},{28,26,26,26},{28,26,30,28},
03461   {28,26,28,30},{28,28,30,24},{30,28,30,30},{30,28,30,30},
03462   {26,28,30,30},{28,28,28,30},{30,28,30,30},{30,28,30,30},
03463   {30,28,30,30},{30,28,30,30},{30,28,30,30},{30,28,30,30},
03464   {30,28,30,30},{30,28,30,30},{30,28,30,30},{30,28,30,30},
03465   {30,28,30,30},{30,28,30,30},{30,28,30,30},{30,28,30,30}
03466 };
03467 #endif
03468 
03469 /*Bulk data for the number of parity bytes per Reed-Solomon block.*/
03470 static const unsigned char QR_RS_NPAR_VALS[71]={
03471   /*[ 0]*/ 7,10,13,17,
03472   /*[ 4]*/10,16,22, 28,26,26, 26,22, 24,22,22, 26,24,18,22,
03473   /*[19]*/15,26,18, 22,24, 30,24,20,24,
03474   /*[28]*/18,16,24, 28, 28, 28,28,30,24,
03475   /*[37]*/20,18, 18,26, 24,28,24, 30,26,28, 28, 26,28,30, 30,22,20,24,
03476   /*[55]*/20,18,26,16,
03477   /*[59]*/20,30,28, 24,22,26, 28,26, 30,28,30,30
03478 };
03479 
03480 /*An offset into QR_RS_NPAR_DATA for each version that gives the number of
03481    parity bytes per Reed-Solomon block for each error correction level.*/
03482 static const unsigned char QR_RS_NPAR_OFFS[40]={
03483    0, 4,19,55,15,28,37,12,51,39,
03484   59,62,10,24,22,41,31,44, 7,65,
03485   47,33,67,67,48,32,67,67,67,67,
03486   67,67,67,67,67,67,67,67,67,67
03487 };
03488 
03489 /*The number of Reed-Solomon blocks for each version and error correction
03490    level.*/
03491 static const unsigned char QR_RS_NBLOCKS[40][4]={
03492   { 1, 1, 1, 1},{ 1, 1, 1, 1},{ 1, 1, 2, 2},{ 1, 2, 2, 4},
03493   { 1, 2, 4, 4},{ 2, 4, 4, 4},{ 2, 4, 6, 5},{ 2, 4, 6, 6},
03494   { 2, 5, 8, 8},{ 4, 5, 8, 8},{ 4, 5, 8,11},{ 4, 8,10,11},
03495   { 4, 9,12,16},{ 4, 9,16,16},{ 6,10,12,18},{ 6,10,17,16},
03496   { 6,11,16,19},{ 6,13,18,21},{ 7,14,21,25},{ 8,16,20,25},
03497   { 8,17,23,25},{ 9,17,23,34},{ 9,18,25,30},{10,20,27,32},
03498   {12,21,29,35},{12,23,34,37},{12,25,34,40},{13,26,35,42},
03499   {14,28,38,45},{15,29,40,48},{16,31,43,51},{17,33,45,54},
03500   {18,35,48,57},{19,37,51,60},{19,38,53,63},{20,40,56,66},
03501   {21,43,59,70},{22,45,62,74},{24,47,65,77},{25,49,68,81}
03502 };
03503 
03504 /*Attempts to fully decode a QR code.
03505   _qrdata:   Returns the parsed code data.
03506   _gf:       Used for Reed-Solomon error correction.
03507   _ul_pos:   The location of the UL finder pattern.
03508   _ur_pos:   The location of the UR finder pattern.
03509   _dl_pos:   The location of the DL finder pattern.
03510   _version:  The (decoded) version number.
03511   _fmt_info: The decoded format info.
03512   _img:      The binary input image.
03513   _width:    The width of the input image.
03514   _height:   The height of the input image.
03515   Return: 0 on success, or a negative value on error.*/
03516 static int qr_code_decode(qr_code_data *_qrdata,const rs_gf256 *_gf,
03517  const qr_point _ul_pos,const qr_point _ur_pos,const qr_point _dl_pos,
03518  int _version,int _fmt_info,
03519  const unsigned char *_img,int _width,int _height){
03520   qr_sampling_grid   grid;
03521   unsigned          *data_bits;
03522   unsigned char    **blocks;
03523   unsigned char     *block_data;
03524   int                nblocks;
03525   int                nshort_blocks;
03526   int                ncodewords;
03527   int                block_sz;
03528   int                ecc_level;
03529   int                ndata;
03530   int                npar;
03531   int                dim;
03532   int                ret;
03533   int                i;
03534   /*Read the bits out of the image.*/
03535   qr_sampling_grid_init(&grid,_version,_ul_pos,_ur_pos,_dl_pos,_qrdata->bbox,
03536    _img,_width,_height);
03537 #if defined(QR_DEBUG)
03538   qr_sampling_grid_dump(&grid,_version,_img,_width,_height);
03539 #endif
03540   dim=17+(_version<<2);
03541   data_bits=(unsigned *)malloc(
03542    dim*(dim+QR_INT_BITS-1>>QR_INT_LOGBITS)*sizeof(*data_bits));
03543   qr_sampling_grid_sample(&grid,data_bits,dim,_fmt_info,_img,_width,_height);
03544   /*Group those bits into Reed-Solomon codewords.*/
03545   ecc_level=(_fmt_info>>3)^1;
03546   nblocks=QR_RS_NBLOCKS[_version-1][ecc_level];
03547   npar=*(QR_RS_NPAR_VALS+QR_RS_NPAR_OFFS[_version-1]+ecc_level);
03548   ncodewords=qr_code_ncodewords(_version);
03549   block_sz=ncodewords/nblocks;
03550   nshort_blocks=nblocks-(ncodewords%nblocks);
03551   blocks=(unsigned char **)malloc(nblocks*sizeof(*blocks));
03552   block_data=(unsigned char *)malloc(ncodewords*sizeof(*block_data));
03553   blocks[0]=block_data;
03554   for(i=1;i<nblocks;i++)blocks[i]=blocks[i-1]+block_sz+(i>nshort_blocks);
03555   qr_samples_unpack(blocks,nblocks,block_sz-npar,nshort_blocks,
03556    data_bits,grid.fpmask,dim);
03557   qr_sampling_grid_clear(&grid);
03558   free(blocks);
03559   free(data_bits);
03560   /*Perform the error correction.*/
03561   ndata=0;
03562   ncodewords=0;
03563   ret=0;
03564   for(i=0;i<nblocks;i++){
03565     int block_szi;
03566     int ndatai;
03567     block_szi=block_sz+(i>=nshort_blocks);
03568     if(rs_correct(_gf,QR_M0,block_data+ncodewords,block_szi,npar,NULL,0)<0){
03569       ret=-1;
03570       break;
03571     }
03572     ndatai=block_szi-npar;
03573     memmove(block_data+ndata,block_data+ncodewords,ndatai*sizeof(*block_data));
03574     ncodewords+=block_szi;
03575     ndata+=ndatai;
03576   }
03577   /*Parse the corrected bitstream.*/
03578   if(ret>=0){
03579     ret=qr_code_data_parse(_qrdata,_version,block_data,ndata);
03580     /*We could return any partially decoded data, but then we'd have to have
03581        API support for that; a mode ignoring ECC errors might also be useful.*/
03582     if(ret<0)qr_code_data_clear(_qrdata);
03583     _qrdata->version=_version;
03584     _qrdata->ecc_level=ecc_level;
03585   }
03586   free(block_data);
03587   return ret;
03588 }
03589 
03590 /*Searches for an arrangement of these three finder centers that yields a valid
03591    configuration.
03592   _c: On input, the three finder centers to consider in any order.
03593   Return: The detected version number, or a negative value on error.*/
03594 static int qr_reader_try_configuration(qr_reader *_reader,
03595  qr_code_data *_qrdata,const unsigned char *_img,int _width,int _height,
03596  qr_finder_center *_c[3]){
03597   int      ci[7];
03598   unsigned maxd;
03599   int      ccw;
03600   int      i0;
03601   int      i;
03602   /*Sort the points in counter-clockwise order.*/
03603   ccw=qr_point_ccw(_c[0]->pos,_c[1]->pos,_c[2]->pos);
03604   /*Colinear points can't be the corners of a quadrilateral.*/
03605   if(!ccw)return -1;
03606   /*Include a few extra copies of the cyclical list to avoid mods.*/
03607   ci[6]=ci[3]=ci[0]=0;
03608   ci[4]=ci[1]=1+(ccw<0);
03609   ci[5]=ci[2]=2-(ccw<0);
03610   /*Assume the points farthest from each other are the opposite corners, and
03611      find the top-left point.*/
03612   maxd=qr_point_distance2(_c[1]->pos,_c[2]->pos);
03613   i0=0;
03614   for(i=1;i<3;i++){
03615     unsigned d;
03616     d=qr_point_distance2(_c[ci[i+1]]->pos,_c[ci[i+2]]->pos);
03617     if(d>maxd){
03618       i0=i;
03619       maxd=d;
03620     }
03621   }
03622   /*However, try all three possible orderings, just to be sure (a severely
03623      skewed projection could move opposite corners closer than adjacent).*/
03624   for(i=i0;i<i0+3;i++){
03625     qr_aff    aff;
03626     qr_hom    hom;
03627     qr_finder ul;
03628     qr_finder ur;
03629     qr_finder dl;
03630     int       res;
03631     int       ur_version;
03632     int       dl_version;
03633     int       fmt_info;
03634     ul.c=_c[ci[i]];
03635     ur.c=_c[ci[i+1]];
03636     dl.c=_c[ci[i+2]];
03637     /*Estimate the module size and version number from the two opposite corners.
03638       The module size is not constant in the image, so we compute an affine
03639        projection from the three points we have to a square domain, and
03640        estimate it there.
03641       Although it should be the same along both axes, we keep separate
03642        estimates to account for any remaining projective distortion.*/
03643     res=QR_INT_BITS-2-QR_FINDER_SUBPREC-qr_ilog(QR_MAXI(_width,_height)-1);
03644     qr_aff_init(&aff,ul.c->pos,ur.c->pos,dl.c->pos,res);
03645     qr_aff_unproject(ur.o,&aff,ur.c->pos[0],ur.c->pos[1]);
03646     qr_finder_edge_pts_aff_classify(&ur,&aff);
03647     if(qr_finder_estimate_module_size_and_version(&ur,1<<res,1<<res)<0)continue;
03648     qr_aff_unproject(dl.o,&aff,dl.c->pos[0],dl.c->pos[1]);
03649     qr_finder_edge_pts_aff_classify(&dl,&aff);
03650     if(qr_finder_estimate_module_size_and_version(&dl,1<<res,1<<res)<0)continue;
03651     /*If the estimated versions are significantly different, reject the
03652        configuration.*/
03653     if(abs(ur.eversion[1]-dl.eversion[0])>QR_LARGE_VERSION_SLACK)continue;
03654     qr_aff_unproject(ul.o,&aff,ul.c->pos[0],ul.c->pos[1]);
03655     qr_finder_edge_pts_aff_classify(&ul,&aff);
03656     if(qr_finder_estimate_module_size_and_version(&ul,1<<res,1<<res)<0||
03657      abs(ul.eversion[1]-ur.eversion[1])>QR_LARGE_VERSION_SLACK||
03658      abs(ul.eversion[0]-dl.eversion[0])>QR_LARGE_VERSION_SLACK){
03659       continue;
03660     }
03661 #if defined(QR_DEBUG)
03662     qr_finder_dump_aff_undistorted(&ul,&ur,&dl,&aff,_img,_width,_height);
03663 #endif
03664     /*If we made it this far, upgrade the affine homography to a full
03665        homography.*/
03666     if(qr_hom_fit(&hom,&ul,&ur,&dl,_qrdata->bbox,&aff,
03667      &_reader->isaac,_img,_width,_height)<0){
03668       continue;
03669     }
03670     qr_hom_unproject(ul.o,&hom,ul.c->pos[0],ul.c->pos[1]);
03671     qr_hom_unproject(ur.o,&hom,ur.c->pos[0],ur.c->pos[1]);
03672     qr_hom_unproject(dl.o,&hom,dl.c->pos[0],dl.c->pos[1]);
03673     qr_finder_edge_pts_hom_classify(&ur,&hom);
03674     if(qr_finder_estimate_module_size_and_version(&ur,
03675      ur.o[0]-ul.o[0],ur.o[0]-ul.o[0])<0){
03676       continue;
03677     }
03678     qr_finder_edge_pts_hom_classify(&dl,&hom);
03679     if(qr_finder_estimate_module_size_and_version(&dl,
03680      dl.o[1]-ul.o[1],dl.o[1]-ul.o[1])<0){
03681       continue;
03682     }
03683 #if defined(QR_DEBUG)
03684     qr_finder_dump_hom_undistorted(&ul,&ur,&dl,&hom,_img,_width,_height);
03685 #endif
03686     /*If we have a small version (less than 7), there's no encoded version
03687        information.
03688       If the estimated version on the two corners matches and is sufficiently
03689        small, we assume this is the case.*/
03690     if(ur.eversion[1]==dl.eversion[0]&&ur.eversion[1]<7){
03691       /*We used to do a whole bunch of extra geometric checks for small
03692          versions, because with just an affine correction, it was fairly easy
03693          to estimate two consistent module sizes given a random configuration.
03694         However, now that we're estimating a full homography, these appear to
03695          be unnecessary.*/
03696 #if 0
03697       static const signed char LINE_TESTS[12][6]={
03698         /*DL left, UL > 0, UR > 0*/
03699         {2,0,0, 1,1, 1},
03700         /*DL right, UL > 0, UR < 0*/
03701         {2,1,0, 1,1,-1},
03702         /*UR top, UL > 0, DL > 0*/
03703         {1,2,0, 1,2, 1},
03704         /*UR bottom, UL > 0, DL < 0*/
03705         {1,3,0, 1,2,-1},
03706         /*UR left, DL < 0, UL < 0*/
03707         {1,0,2,-1,0,-1},
03708         /*UR right, DL > 0, UL > 0*/
03709         {1,1,2, 1,0, 1},
03710         /*DL top, UR < 0, UL < 0*/
03711         {2,2,1,-1,0,-1},
03712         /*DL bottom, UR > 0, UL > 0*/
03713         {2,3,1, 1,0, 1},
03714         /*UL left, DL > 0, UR > 0*/
03715         {0,0,2, 1,1, 1},
03716         /*UL right, DL > 0, UR < 0*/
03717         {0,1,2, 1,1,-1},
03718         /*UL top, UR > 0, DL > 0*/
03719         {0,2,1, 1,2, 1},
03720         /*UL bottom, UR > 0, DL < 0*/
03721         {0,3,1, 1,2,-1}
03722       };
03723       qr_finder *f[3];
03724       int        j;
03725       /*Start by decoding the format information.
03726         This is cheap, but unlikely to reject invalid configurations.
03727         56.25% of all bitstrings are valid, and we mix and match several pieces
03728          until we find a valid combination, so our real chances of finding a
03729          valid codeword in random bits are even higher.*/
03730       fmt_info=qr_finder_fmt_info_decode(&ul,&ur,&dl,&aff,_img,_width,_height);
03731       if(fmt_info<0)continue;
03732       /*Now we fit lines to the edges of each finder pattern and check to make
03733          sure the centers of the other finder patterns lie on the proper side.*/
03734       f[0]=&ul;
03735       f[1]=&ur;
03736       f[2]=&dl;
03737       for(j=0;j<12;j++){
03738         const signed char *t;
03739         qr_line            l0;
03740         int               *p;
03741         t=LINE_TESTS[j];
03742         qr_finder_ransac(f[t[0]],&aff,&_reader->isaac,t[1]);
03743         /*We may not have enough points to fit a line accurately here.
03744           If not, we just skip the test.*/
03745         if(qr_line_fit_finder_edge(l0,f[t[0]],t[1],res)<0)continue;
03746         p=f[t[2]]->c->pos;
03747         if(qr_line_eval(l0,p[0],p[1])*t[3]<0)break;
03748         p=f[t[4]]->c->pos;
03749         if(qr_line_eval(l0,p[0],p[1])*t[5]<0)break;
03750       }
03751       if(j<12)continue;
03752       /*All tests passed.*/
03753 #endif
03754       ur_version=ur.eversion[1];
03755     }
03756     else{
03757       /*If the estimated versions are significantly different, reject the
03758          configuration.*/
03759       if(abs(ur.eversion[1]-dl.eversion[0])>QR_LARGE_VERSION_SLACK)continue;
03760       /*Otherwise we try to read the actual version data from the image.
03761         If the real version is not sufficiently close to our estimated version,
03762          then we assume there was an unrecoverable decoding error (so many bit
03763          errors we were within 3 errors of another valid code), and throw that
03764          value away.
03765         If no decoded version could be sufficiently close, we don't even try.*/
03766       if(ur.eversion[1]>=7-QR_LARGE_VERSION_SLACK){
03767         ur_version=qr_finder_version_decode(&ur,&hom,_img,_width,_height,0);
03768         if(abs(ur_version-ur.eversion[1])>QR_LARGE_VERSION_SLACK)ur_version=-1;
03769       }
03770       else ur_version=-1;
03771       if(dl.eversion[0]>=7-QR_LARGE_VERSION_SLACK){
03772         dl_version=qr_finder_version_decode(&dl,&hom,_img,_width,_height,1);
03773         if(abs(dl_version-dl.eversion[0])>QR_LARGE_VERSION_SLACK)dl_version=-1;
03774       }
03775       else dl_version=-1;
03776       /*If we got at least one valid version, or we got two and they match,
03777          then we found a valid configuration.*/
03778       if(ur_version>=0){
03779         if(dl_version>=0&&dl_version!=ur_version)continue;
03780       }
03781       else if(dl_version<0)continue;
03782       else ur_version=dl_version;
03783     }
03784     qr_finder_edge_pts_hom_classify(&ul,&hom);
03785     if(qr_finder_estimate_module_size_and_version(&ul,
03786      ur.o[0]-dl.o[0],dl.o[1]-ul.o[1])<0||
03787      abs(ul.eversion[1]-ur.eversion[1])>QR_SMALL_VERSION_SLACK||
03788      abs(ul.eversion[0]-dl.eversion[0])>QR_SMALL_VERSION_SLACK){
03789       continue;
03790     }
03791     fmt_info=qr_finder_fmt_info_decode(&ul,&ur,&dl,&hom,_img,_width,_height);
03792     if(fmt_info<0)continue;
03793     if(qr_code_decode(_qrdata,&_reader->gf,ul.c->pos,ur.c->pos,dl.c->pos,
03794      ur_version,fmt_info,_img,_width,_height)<0){
03795       /*TODO: Maybe somebody flipped the code?
03796         We'd still get a valid version, and probably valid (but incorrect)
03797          format info.
03798         After we've come this far, it should be a simple matter to check.*/
03799       continue;
03800     }
03801     return ur_version;
03802   }
03803   return -1;
03804 }
03805 
03806 void qr_reader_match_centers(qr_reader *_reader,qr_code_data_list *_qrlist,
03807  qr_finder_center *_centers,int _ncenters,
03808  const unsigned char *_img,int _width,int _height){
03809   /*The number of centers should be small, so an O(n^3) exhaustive search of
03810      which ones go together should be reasonable.*/
03811   unsigned char *mark;
03812   int            i;
03813   int            j;
03814   int            k;
03815   mark=(unsigned char *)calloc(_ncenters,sizeof(*mark));
03816   for(i=0;i<_ncenters;i++){
03817     /*TODO: We might be able to accelerate this step significantly by
03818        considering the remaining finder centers in a more intelligent order,
03819        based on the first finder center we just chose.*/
03820     for(j=i+1;!mark[i]&&j<_ncenters;j++){
03821       for(k=j+1;!mark[j]&&k<_ncenters;k++)if(!mark[k]){
03822         qr_finder_center *c[3];
03823         qr_code_data      qrdata;
03824         int               version;
03825         c[0]=_centers+i;
03826         c[1]=_centers+j;
03827         c[2]=_centers+k;
03828         version=qr_reader_try_configuration(_reader,&qrdata,
03829          _img,_width,_height,c);
03830         if(version>=0){
03831           int ninside;
03832           int l;
03833           /*Add the data to the list.*/
03834           qr_code_data_list_add(_qrlist,&qrdata);
03835           /*Convert the bounding box we're returning to the user to normal
03836              image coordinates.*/
03837           for(l=0;l<4;l++){
03838             _qrlist->qrdata[_qrlist->nqrdata-1].bbox[l][0]>>=QR_FINDER_SUBPREC;
03839             _qrlist->qrdata[_qrlist->nqrdata-1].bbox[l][1]>>=QR_FINDER_SUBPREC;
03840           }
03841           /*Mark these centers as used.*/
03842           mark[i]=mark[j]=mark[k]=1;
03843           /*Find any other finder centers located inside this code.*/
03844           for(l=ninside=0;l<_ncenters;l++)if(!mark[l]){
03845             if(qr_point_ccw(qrdata.bbox[0],qrdata.bbox[1],_centers[l].pos)>=0&&
03846              qr_point_ccw(qrdata.bbox[1],qrdata.bbox[3],_centers[l].pos)>=0&&
03847              qr_point_ccw(qrdata.bbox[3],qrdata.bbox[2],_centers[l].pos)>=0&&
03848              qr_point_ccw(qrdata.bbox[2],qrdata.bbox[0],_centers[l].pos)>=0){
03849               mark[l]=2;
03850               ninside++;
03851             }
03852           }
03853           if(ninside>=3){
03854             /*We might have a "Double QR": a code inside a code.
03855               Copy the relevant centers to a new array and do a search confined
03856                to that subset.*/
03857             qr_finder_center *inside;
03858             inside=(qr_finder_center *)malloc(ninside*sizeof(*inside));
03859             for(l=ninside=0;l<_ncenters;l++){
03860               if(mark[l]==2)*&inside[ninside++]=*&_centers[l];
03861             }
03862             qr_reader_match_centers(_reader,_qrlist,inside,ninside,
03863              _img,_width,_height);
03864             free(inside);
03865           }
03866           /*Mark _all_ such centers used: codes cannot partially overlap.*/
03867           for(l=0;l<_ncenters;l++)if(mark[l]==2)mark[l]=1;
03868         }
03869       }
03870     }
03871   }
03872   free(mark);
03873 }
03874 
03875 int _zbar_qr_found_line (qr_reader *reader,
03876                          int dir,
03877                          const qr_finder_line *line)
03878 {
03879     /* minimally intrusive brute force version */
03880     qr_finder_lines *lines = &reader->finder_lines[dir];
03881 
03882     if(lines->nlines >= lines->clines) {
03883         lines->clines *= 2;
03884         lines->lines = realloc(lines->lines,
03885                                ++lines->clines * sizeof(*lines->lines));
03886     }
03887 
03888     memcpy(lines->lines + lines->nlines++, line, sizeof(*line));
03889 
03890     return(0);
03891 }
03892 
03893 static inline void qr_svg_centers (const qr_finder_center *centers,
03894                                    int ncenters)
03895 {
03896     int i, j;
03897     svg_path_start("centers", 1, 0, 0);
03898     for(i = 0; i < ncenters; i++)
03899         svg_path_moveto(SVG_ABS, centers[i].pos[0], centers[i].pos[1]);
03900     svg_path_end();
03901 
03902     svg_path_start("edge-pts", 1, 0, 0);
03903     for(i = 0; i < ncenters; i++) {
03904         const qr_finder_center *cen = centers + i;
03905         for(j = 0; j < cen->nedge_pts; j++)
03906             svg_path_moveto(SVG_ABS,
03907                             cen->edge_pts[j].pos[0], cen->edge_pts[j].pos[1]);
03908     }
03909     svg_path_end();
03910 }
03911 
03912 int _zbar_qr_decode (qr_reader *reader,
03913                      zbar_image_scanner_t *iscn,
03914                      zbar_image_t *img)
03915 {
03916     int nqrdata = 0;
03917     qr_finder_edge_pt *edge_pts = NULL;
03918     qr_finder_center *centers = NULL;
03919 
03920     if(reader->finder_lines[0].nlines < 9 ||
03921        reader->finder_lines[1].nlines < 9)
03922         return(0);
03923 
03924     svg_group_start("finder", 0, 1. / (1 << QR_FINDER_SUBPREC), 0, 0, 0);
03925 
03926     int ncenters = qr_finder_centers_locate(&centers, &edge_pts, reader, 0, 0);
03927 
03928     zprintf(14, "%dx%d finders, %d centers:\n",
03929             reader->finder_lines[0].nlines,
03930             reader->finder_lines[1].nlines,
03931             ncenters);
03932     qr_svg_centers(centers, ncenters);
03933 
03934     if(ncenters >= 3) {
03935         void *bin = qr_binarize(img->data, img->width, img->height);
03936 
03937         qr_code_data_list qrlist;
03938         qr_code_data_list_init(&qrlist);
03939 
03940         qr_reader_match_centers(reader, &qrlist, centers, ncenters,
03941                                 bin, img->width, img->height);
03942 
03943         if(qrlist.nqrdata > 0)
03944             nqrdata = qr_code_data_list_extract_text(&qrlist, iscn, img);
03945 
03946         qr_code_data_list_clear(&qrlist);
03947         free(bin);
03948     }
03949     svg_group_end();
03950 
03951     if(centers)
03952         free(centers);
03953     if(edge_pts)
03954         free(edge_pts);
03955     return(nqrdata);
03956 }
03957