+typedef struct {
+ int v, e;
+} Traverse;
+
+static void traverse_next(Traverse *t) {
+ int v2;
+
+ if (t->v<0)
+ return;
+
+ v2= EDGE_END2(t->v, t->e);
+ if (v2>=0) {
+ int e2= edge_reverse(t->v, t->e);
+ assert(EDGE_END2(v2,e2) == t->v);
+ t->e= (e2+3) % V6;
+ }
+ t->v= v2;
+}
+
+static void interpolate_lin4point(void) {
+ /* four points P Q R S, although P and S may be missing
+ * interpolate in QR finding M. */
+ int xq,yq,eqr, vp,vq,vm,vr,vs, k;
+ double pqtarg[D3], srtarg[D3];
+
+ if (inc != 2) fail("cannot do >2x lin4point interpolation\n");
+
+ for (eqr=1; eqr!=4; eqr+=5, eqr%=V6) { /* each old edge exactly once */
+ printf("eqr=%d\n",eqr);
+ for (yq=0; yq<Y; yq+=inc) {
+ printf(" yq=%2d ",yq);
+ for (xq=((yq>>1)&1); xq<X; xq+=inc) {
+ vq= yq << YSHIFT | xq;
+ Traverse trav; trav.v=vq; trav.e=eqr;
+
+ traverse_next(&trav);
+ vm= trav.v;
+ if (vm<0) continue;
+
+ traverse_next(&trav);
+ vr= trav.v;
+ assert(vr>=0);
+
+ traverse_next(&trav);
+ traverse_next(&trav);
+ vs= trav.v;
+
+ trav.v= vq; trav.e= EDGE_OPPOSITE(eqr);
+ traverse_next(&trav);
+ traverse_next(&trav);
+ vp= trav.v;
+
+ printf(" 0x%02x-##-%02x-!%02x!-%02x-##-%02x",
+ vp&0xff,vq,vm,vr,vs&0xff);
+
+ if (vp>=0)
+ K pqtarg[k]= all.a[vq][k]*2 - all.a[vp][k];
+ else
+ K pqtarg[k]= all.a[vr][k];
+
+ if (vs>=0)
+ K srtarg[k]= all.a[vr][k]*2 - all.a[vs][k];
+ else
+ K srtarg[k]= all.a[vq][k];