pdollar0.d

   1 /*
   2  * The $P+ Point-Cloud Recognizer (JavaScript version)
   3  *
   4  *  Radu-Daniel Vatavu, Ph.D.
   5  *  University Stefan cel Mare of Suceava
   6  *  Suceava 720229, Romania
   7  *  vatavu@eed.usv.ro
   8  *
   9  * The academic publication for the $P+ recognizer, and what should be
  10  * used to cite it, is:
  11  *
  12  *     Vatavu, R.-D. (2017). Improving gesture recognition accuracy on
  13  *     touch screens for users with low vision. Proceedings of the ACM
  14  *     Conference on Human Factors in Computing Systems (CHI '17). Denver,
  15  *     Colorado (May 6-11, 2017). New York: ACM Press, pp. 4667-4679.
  16  *     https://dl.acm.org/citation.cfm?id=3025941
  17  *
  18  * This software is distributed under the "New BSD License" agreement:
  19  *
  20  * Copyright (C) 2017-2018, Radu-Daniel Vatavu and Jacob O. Wobbrock. All
  21  * rights reserved. Last updated July 14, 2018.
  22  *
  23  * Redistribution and use in source and binary forms, with or without
  24  * modification, are permitted provided that the following conditions are met:
  25  *    * Redistributions of source code must retain the above copyright
  26  *      notice, this list of conditions and the following disclaimer.
  27  *    * Redistributions in binary form must reproduce the above copyright
  28  *      notice, this list of conditions and the following disclaimer in the
  29  *      documentation and/or other materials provided with the distribution.
  30  *    * Neither the name of the University Stefan cel Mare of Suceava, nor the
  31  *      names of its contributors may be used to endorse or promote products
  32  *      derived from this software without specific prior written permission.
  33  *
  34  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
  35  * IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
  36  * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
  37  * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL Radu-Daniel Vatavu OR Lisa Anthony
  38  * OR Jacob O. Wobbrock BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
  39  * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
  40  * OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
  41  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
  42  * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  43  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  44  * SUCH DAMAGE.
  45  */
  46 module pdollar0;
  47
  48 // use $P instead of $P+?
  49 // doesn't work right yet (something is wrong with the weights)
  50 // maybe this is due to the different resampler?
  51 // or due to the bug in angle calculation?
  52 version = use_older_dollar_p;
  53
  54 import iv.alice;
  55 import iv.cmdcon;
  56 import iv.vfs;
  57
  58
  59 // ////////////////////////////////////////////////////////////////////////// //
  60 struct DPResult {
  61   string name;
  62   float score;
  63   @property bool valid () const nothrow @safe @nogc { pragma(inline, true); import std.math : isNaN; return !score.isNaN; }
  64 }
  65
  66
  67 // ////////////////////////////////////////////////////////////////////////// //
  68 struct DPPoint {
  69   float x=0.0f, y=0.0f;
  70   uint id = 0; // stroke ID to which this point belongs (1,2,...)
  71   version(use_older_dollar_p) {} else {
  72     float angle=0.0f;
  73   }
  74 }
  75
  76
  77 final class DPPointCloud {
  78 public:
  79   enum NumPoints = 32;
  80
  81 public:
  82   string name;
  83   DPPoint[NumPoints] points = DPPoint(0.0f, 0.0f);
  84
  85 private:
  86   this () {}
  87
  88 public:
  89   this (string aname, const(DPPoint)[] pts...) nothrow @nogc {
  90     name = aname;
  91     resample(points, pts);
  92     scale(points);
  93     translateTo(points, DPPoint(0.0f, 0.0f));
  94     version(use_older_dollar_p) {} else computeNormalizedTurningAngles(points);
  95   }
  96
  97   void save (VFile fl) const {
  98     string nn = name;
  99     if (nn.length > 65535) nn = nn[0..65535]; // fuck you
 100     if (nn.length > 254) {
 101       fl.writeNum!ubyte(255);
 102       fl.writeNum!ushort(cast(ushort)nn.length);
 103     } else {
 104       fl.writeNum!ubyte(cast(ubyte)nn.length);
 105     }
 106     fl.rawWriteExact(nn);
 107     fl.writeNum!ubyte(NumPoints);
 108     foreach (const ref DPPoint pt; points[]) {
 109       fl.writeNum!float(cast(float)pt.x);
 110       fl.writeNum!float(cast(float)pt.y);
 111       fl.writeNum!uint(pt.id);
 112     }
 113   }
 114
 115   void load (VFile fl) {
 116     char[] nn;
 117     ubyte len = fl.readNum!ubyte;
 118     if (len == 255) {
 119       ushort xlen = fl.readNum!ushort;
 120       nn.length = xlen;
 121     } else {
 122       nn.length = len;
 123     }
 124     fl.rawReadExact(nn);
 125     name = cast(string)nn; // it is safe to cast here
 126     if (fl.readNum!ubyte != NumPoints) throw new Exception("invalid number of points in cloud");
 127     foreach (ref DPPoint pt; points[]) {
 128       pt.x = cast(float)fl.readNum!float;
 129       pt.y = cast(float)fl.readNum!float;
 130       pt.id = fl.readNum!uint;
 131     }
 132     // perform last $P+ step
 133     version(use_older_dollar_p) {} else computeNormalizedTurningAngles(points);
 134   }
 135
 136   static DPPointCloud loadNew (VFile fl) {
 137     auto res = new DPPointCloud();
 138     res.load(fl);
 139     return res;
 140   }
 141 }
 142
 143
 144 // ////////////////////////////////////////////////////////////////////////// //
 145 final class DPGestureList {
 146 private:
 147   DPPointCloud[] gestures;
 148
 149 public:
 150   this () nothrow @nogc {}
 151
 152   void clear () {
 153     delete gestures;
 154   }
 155
 156   string[] knownGestureNames () {
 157     string[] res;
 158     mainloop: foreach (DPPointCloud gst; gestures) {
 159       foreach (string s; res) if (s == gst.name) continue mainloop;
 160       res ~= gst.name;
 161     }
 162     return res;
 163   }
 164
 165   void appendGesture (string aname, const(DPPoint)[] pts...) {
 166     appendGesture(new DPPointCloud(aname, pts));
 167   }
 168
 169   // you can add as many gestures with same name as you want to
 170   void appendGesture (DPPointCloud gg) {
 171     import core.memory : GC;
 172     if (gg is null) return;
 173     auto optr = gestures.ptr;
 174     gestures ~= gg;
 175     if (gestures.ptr !is optr) {
 176       optr = gestures.ptr;
 177       if (optr is GC.addrOf(optr)) GC.setAttr(optr, GC.BlkAttr.NO_INTERIOR);
 178     }
 179   }
 180
 181   // remove *ALL* gestures with the given name
 182   void removeGesture (const(char)[] name) {
 183     usize idx = 0;
 184     while (idx < gestures.length) {
 185       if (gestures[idx].name == name) {
 186         foreach (immutable c; idx+1..gestures.length) gestures[c-1] = gestures[c];
 187         gestures[$-1] = null;
 188         gestures.length -= 1;
 189         gestures.assumeSafeAppend;
 190       } else {
 191         ++idx;
 192       }
 193     }
 194   }
 195
 196   // if you have more than one gesture with the same name, keep increasing `idx` to get more, until you'll get `null`
 197   DPPointCloud findGesture (const(char)[] name, uint idx=0) nothrow @nogc {
 198     foreach (DPPointCloud g; gestures) {
 199       if (g.name == name) {
 200         if (idx-- == 0) return g;
 201       }
 202     }
 203     return null;
 204   }
 205
 206   DPResult recognize (const(DPPoint)[] origpoints) nothrow @nogc {
 207     import std.algorithm : max;
 208     if (origpoints.length < 2) return DPResult.init;
 209     DPPoint[DPPointCloud.NumPoints] points;
 210     resample(points, origpoints);
 211     scale(points);
 212     translateTo(points, DPPoint(0.0f, 0.0f));
 213     version(use_older_dollar_p) {} else computeNormalizedTurningAngles(points);
 214     float b = float.infinity;
 215     int u = -1;
 216     foreach (immutable idx, DPPointCloud gst; gestures) {
 217       version(use_older_dollar_p) {
 218         immutable float d = greedyCloudMatch(points, gst);
 219       } else {
 220         immutable float d0 = cloudDistance(points, gst.points);
 221         immutable float d1 = cloudDistance(gst.points, points);
 222         immutable float d = (d0 < d1 ? d0 : d1);
 223       }
 224       if (d < b) {
 225         b = d; // best (least) distance
 226         u = cast(int)idx; // point-cloud
 227       }
 228     }
 229     if (u == -1) return DPResult();
 230     { import core.stdc.stdio; printf("b=%f (%f) (%f)\n", b, (b-2.0)/-2.0, (b > 1.0f ? 1.0f/b : 1.0f)); }
 231     version(use_older_dollar_p) {
 232       return (u != -1 ? DPResult(gestures[u].name, max((b-2.0f)/-2.0f, 0.0f)) : DPResult.init);
 233     } else {
 234       //return (u != -1 ? DPResult(gestures[u].name, max((b-2.0f)/-2.0f, 0.0f)) : DPResult.init);
 235       return (u != -1 ? DPResult(gestures[u].name, (b > 1.0f ? 1.0f/b : 1.0f)) : DPResult.init);
 236     }
 237   }
 238
 239 public:
 240   enum Signature = "DOLP8LB0";
 241
 242   void save (VFile fl) const {
 243     fl.rawWriteExact(Signature);
 244     fl.writeNum!uint(cast(uint)gestures.length);
 245     foreach (const DPPointCloud gst; gestures) gst.save(fl);
 246   }
 247
 248   void load (VFile fl) {
 249     delete gestures;
 250     char[Signature.length] sign;
 251     fl.rawReadExact(sign);
 252     if (sign[0..$-1] != Signature[0..$-1]) throw new Exception("invalid $P library signature");
 253     if (sign[$-1] != Signature[$-1]) throw new Exception("invalid $P library version");
 254     uint count = fl.readNum!uint;
 255     if (count > uint.max/16) throw new Exception("too many gestures in library");
 256     gestures.reserve(count);
 257     foreach (immutable idx; 0..count) appendGesture(DPPointCloud.loadNew(fl));
 258   }
 259 }
 260
 261
 262 // ////////////////////////////////////////////////////////////////////////// //
 263 private:
 264
 265 version(use_older_dollar_p) {
 266   float greedyCloudMatch(size_t n) (in ref DPPoint[n] points, in DPPointCloud P) nothrow @nogc {
 267     import std.algorithm : min;
 268     import std.math : floor, pow;
 269     enum e = cast(float)0.50;
 270     int step = cast(int)floor(pow(points.length, 1-e));
 271     assert(step > 0);
 272     float minv = float.infinity;
 273     for (int i = 0; i < points.length; i += step) {
 274       auto d1 = cloudDistanceOld(points, P.points, i);
 275       auto d2 = cloudDistanceOld(P.points, points, i);
 276       minv = min(minv, d1, d2);
 277     }
 278     return minv;
 279   }
 280
 281   float cloudDistanceOld(size_t n) (in ref DPPoint[n] pts1, in ref DPPoint[n] pts2, int start) nothrow @nogc {
 282     import std.math : sqrt;
 283     bool[n] matched = false;
 284     float sum = 0;
 285     int i = start;
 286     do {
 287       int index = -1;
 288       float minv = float.infinity;
 289       foreach (immutable j, bool mtv; matched[]) {
 290         if (!mtv) {
 291           auto d = distanceSqr(pts1[i], pts2[j]);
 292           if (d < minv) { minv = d; index = cast(int)j; }
 293         }
 294       }
 295       assert(index >= 0);
 296       matched[index] = true;
 297       float weight = cast(float)1-((i-start+cast(int)pts1.length)%cast(int)pts1.length)/cast(float)pts1.length;
 298       sum += weight*sqrt(minv);
 299       i = (i+1)%cast(int)pts1.length;
 300     } while (i != start);
 301     return sum;
 302   }
 303 } else {
 304   // $P+
 305   float cloudDistance(size_t n) (in ref DPPoint[n] pts1, in ref DPPoint[n] pts2) nothrow @nogc {
 306     import std.math : sqrt;
 307     bool[n] matched = false;
 308     size_t mtcount = 0;
 309     float sum = 0.0f;
 310
 311     foreach (immutable i; 0..pts1.length) {
 312       int index = -1;
 313       float minv = float.infinity;
 314       foreach (immutable j; 0..pts2.length) {
 315         immutable float d = distanceWithAngleSqr(pts1[i], pts2[j]);
 316         if (d < minv) {
 317           minv = d;
 318           index = cast(int)j;
 319         }
 320       }
 321       if (!matched[index]) {
 322         ++mtcount;
 323         matched[index] = true;
 324       }
 325       sum += cast(float)sqrt(minv);
 326     }
 327
 328     if (mtcount != matched.length) {
 329       foreach (immutable j, bool mtv; matched[]) {
 330         if (!mtv) {
 331           float minv = float.infinity;
 332           foreach (immutable i; 0..pts1.length) {
 333             immutable float d = distanceWithAngleSqr(pts1[i], pts2[j]);
 334             if (d < minv) minv = d;
 335           }
 336           sum += cast(float)sqrt(minv);
 337         }
 338       }
 339     }
 340
 341     return sum;
 342   }
 343 }
 344
 345
 346 void resample(size_t n) (ref DPPoint[n] newpoints, in DPPoint[] points) nothrow @nogc {
 347   import std.algorithm : max, min;
 348   import std.math : isNaN;
 349   assert(n > 0);
 350   assert(points.length > 1);
 351   immutable float I = pathLength(points)/(n-1); // interval length
 352   float D = 0.0f;
 353   uint nppos = 0;
 354   newpoints[nppos++] = points[0];
 355   foreach (immutable idx; 1..points.length) {
 356     if (points[idx].id == points[idx-1].id) {
 357       auto d = distance(points[idx-1], points[idx]);
 358       if (D+d >= I) {
 359         DPPoint firstPoint = points[idx-1];
 360         while (D+d >= I) {
 361           // add interpolated point
 362           float t = min(max((I-D)/d, 0.0f), 1.0f);
 363           if (isNaN(t)) t = 0.5f;
 364           newpoints[nppos++] = DPPoint((cast(float)1-t)*firstPoint.x+t*points[idx].x,
 365                                        (cast(float)1-t)*firstPoint.y+t*points[idx].y, points[idx].id);
 366           // update partial length
 367           d = D+d-I;
 368           D = 0.0f;
 369           firstPoint = newpoints[nppos-1];
 370         }
 371         D = d;
 372       } else {
 373         D += d;
 374       }
 375     }
 376   }
 377   if (nppos == n-1) {
 378     // sometimes we fall a rounding-error short of adding the last point, so add it if so
 379     newpoints[nppos++] = DPPoint(points[$-1].x, points[$-1].y, points[$-1].id);
 380   }
 381   assert(nppos == n);
 382 }
 383
 384
 385 void scale(size_t n) (ref DPPoint[n] points) nothrow @nogc {
 386   import std.algorithm : max, min;
 387   float minX = +float.infinity;
 388   float minY = +float.infinity;
 389   float maxX = -float.infinity;
 390   float maxY = -float.infinity;
 391   foreach (ref DPPoint p; points[]) {
 392     minX = min(minX, p.x);
 393     minY = min(minY, p.y);
 394     maxX = max(maxX, p.x);
 395     maxY = max(maxY, p.y);
 396   }
 397   immutable float size = max(maxX-minX, maxY-minY);
 398   foreach (ref DPPoint p; points[]) {
 399     p.x = (p.x-minX)/size;
 400     p.y = (p.y-minY)/size;
 401   }
 402 }
 403
 404
 405 // translates points' centroid
 406 void translateTo(size_t n) (ref DPPoint[n] points, DPPoint pt) nothrow @nogc {
 407   auto c = centroid(points);
 408   foreach (ref DPPoint p; points[]) {
 409     p.x = p.x+pt.x-c.x;
 410     p.y = p.y+pt.y-c.y;
 411   }
 412 }
 413
 414
 415 DPPoint centroid(size_t n) (in ref DPPoint[n] points) nothrow @nogc {
 416   float x = 0.0f, y = 0.0f;
 417   foreach (const ref DPPoint p; points[]) {
 418     x += p.x;
 419     y += p.y;
 420   }
 421   immutable float pl = cast(float)1/cast(float)points.length;
 422   return DPPoint(x*pl, y*pl);
 423 }
 424
 425
 426 // modifies `points` angles
 427 void computeNormalizedTurningAngles(size_t n) (ref DPPoint[n] points) nothrow @nogc {
 428   import std.math : acos, PI;
 429   points[0].angle = 0.0f;
 430   foreach (immutable idx; 1..n-1) {
 431     immutable float dx = (points[idx+1].x-points[idx].x)*(points[idx].x-points[idx-1].x);
 432     immutable float dy = (points[idx+1].y-points[idx].y)*(points[idx].y-points[idx-1].y);
 433     immutable float dn = distance(points[idx+1], points[idx])*distance(points[idx], points[idx-1]);
 434     //float cosangle = Math.max(-1.0, Math.min(1.0, (dx + dy) / dn)); // ensure [-1,+1]
 435     float cosangle = (dx+dy)/dn;
 436     // ensure [-1,+1]
 437     if (cosangle < -1.0f) cosangle = -1.0f; else if (cosangle > 1.0f) cosangle = 1.0f;
 438     immutable float angle = acos(cosangle)/cast(float)PI; // normalized angle
 439     points[idx].angle = angle;
 440   }
 441   // last point
 442   points[n-1].angle = 0.0f;
 443 }
 444
 445
 446 // length traversed by a point path
 447 float pathLength (const(DPPoint)[] points) nothrow @nogc {
 448   float d = 0;
 449   foreach (immutable idx; 1..points.length) {
 450     if (points[idx].id == points[idx-1].id) d += distance(points[idx-1], points[idx]);
 451   }
 452   return d;
 453 }
 454
 455
 456 // Euclidean distance between two points
 457 float distance() (in auto ref DPPoint p1, in auto ref DPPoint p2) nothrow @nogc {
 458   import std.math : sqrt;
 459   immutable float dx = p2.x-p1.x;
 460   immutable float dy = p2.y-p1.y;
 461   return cast(float)sqrt(dx*dx+dy*dy);
 462 }
 463
 464
 465 version(use_older_dollar_p) {
 466   // Euclidean distance between two points
 467   float distanceSqr() (in auto ref DPPoint p1, in auto ref DPPoint p2) nothrow @nogc {
 468     immutable float dx = p2.x-p1.x;
 469     immutable float dy = p2.y-p1.y;
 470     return cast(float)(dx*dx+dy*dy);
 471   }
 472 } else {
 473   float distanceWithAngleSqr() (in auto ref DPPoint p1, in auto ref DPPoint p2) nothrow @nogc {
 474     import std.math : sqrt;
 475     immutable float dx = p2.x-p1.x;
 476     immutable float dy = p2.y-p1.y;
 477     immutable float da = p2.angle-p1.angle;
 478     //return cast(float)sqrt(dx*dx+dy*dy+da*da);
 479     return cast(float)(dx*dx+dy*dy+da*da);
 480   }
 481 }