chiark / gitweb /
cca777f8cff22ded6570402e06cfbb3fe2fe5b59
[moebius2.git] / energy.c
1 /*
2  * We try to find an optimal triangle grid
3  */
4
5 #include "common.h"
6 #include "minimise.h"
7 #include "mgraph.h"
8
9 double vertex_areas[N], vertex_mean_edge_lengths[N], edge_lengths[N][V6];
10
11 static double best_energy= DBL_MAX;
12
13 static void addcost(double *energy, double tweight, double tcost, int pr);
14 #define COST(weight, compute) addcost(&energy, (weight), (compute), printing)
15
16 void energy_init(void) {
17 }
18
19 /*---------- main energy computation and subroutines ----------*/
20
21 double compute_energy(const struct Vertices *vs) {
22   static int bests_unprinted;
23   
24   double energy;
25   int printing;
26
27   compute_edge_lengths(vs->a);
28   compute_vertex_areas(vs->a);
29   energy= 0;
30
31   printing= printing_check(pr_cost,0);
32
33   if (printing) printf("%15lld c>e |", evaluations);
34
35   if (XBITS==3) {
36     COST(  3e2,   line_bending_cost(vs->a));
37     COST(  1e3,   edge_length_variation_cost(vs->a));
38     COST( 0.4e3,  rim_proximity_cost(vs->a));
39     COST(  1e6,   edge_angle_cost(vs->a));
40 //    COST(  1e1,   small_triangles_cost(vs->a));
41     COST(  1e12,   noncircular_rim_cost(vs->a));
42     stop_epsilon= 1e-6;
43   } else if (XBITS==4) {
44     COST(  3e2,   line_bending_cost(vs->a));
45     COST(  3e3,   edge_length_variation_cost(vs->a));
46     COST( 9.0e1,  rim_proximity_cost(vs->a)); // 5e1 is too much
47                                                  // 2.5e1 is too little
48     // 0.2e1 grows compared to previous ?
49     // 0.6e0 shrinks compared to previous ?
50     COST(  1e12,   edge_angle_cost(vs->a));
51     COST(  1e12,   noncircular_rim_cost(vs->a));
52     stop_epsilon= 1e-5;
53   } else {
54     abort();
55   }
56
57   if (printing) printf("| total %# e |", energy);
58
59   if (energy < best_energy) {
60     FILE *best_f;
61     int r;
62
63     if (printing) {
64       printf(" BEST");
65       if (bests_unprinted) printf(" [%4d]",bests_unprinted);
66       bests_unprinted= 0;
67     } else {
68       bests_unprinted++;
69     }
70
71     best_f= fopen(best_file_tmp,"wb");  if (!best_f) diee("fopen new out");
72     r= fwrite(vs->a,sizeof(vs->a),1,best_f); if (r!=1) diee("fwrite");
73     if (fclose(best_f)) diee("fclose new best");
74     if (rename(best_file_tmp,best_file)) diee("rename install new best");
75
76     best_energy= energy;
77   }
78   if (printing) {
79     putchar('\n');
80     flushoutput();
81   }
82
83   evaluations++;
84   return energy;
85 }
86
87 static void addcost(double *energy, double tweight, double tcost, int pr) {
88   double tenergy= tweight * tcost;
89   if (pr) printf(" %# e x %g > %# e* |", tcost, tweight, tenergy);
90   *energy += tenergy;
91 }
92
93 /*---------- Precomputations ----------*/
94
95 void compute_edge_lengths(const Vertices vertices) {
96   int v1,e,v2;
97
98   FOR_EDGE(v1,e,v2)
99     edge_lengths[v1][e]= hypotD(vertices[v1],vertices[v2]);
100 }
101
102 void compute_vertex_areas(const Vertices vertices) {
103   int v0,v1,v2, e1,e2;
104 //  int k;
105
106   FOR_VERTEX(v0) {
107     double total= 0.0, edges_total=0;
108     int count= 0;
109
110     FOR_VEDGE(v0,e1,v1) {
111       e2= (e1+1) % V6;
112       v2= EDGE_END2(v0,e2);
113       if (v2<0) continue;
114
115       edges_total += edge_lengths[v0][e1];
116
117 //      double e1v[D3], e2v[D3], av[D3];
118 //      K {
119 //      e1v[k]= vertices[v1][k] - vertices[v0][k];
120 //      e2v[k]= vertices[v2][k] - vertices[v0][k];
121 //      }
122 //      xprod(av, e1v, e2v);
123 //      total += magnD(av);
124
125       count++;
126     }
127     vertex_areas[v0]= total / count;
128     vertex_mean_edge_lengths[v0]= edges_total / count;
129   }
130 }
131
132 /*---------- Edgewise vertex displacement ----------*/
133
134   /*
135    * Definition:
136    *
137    *    At each vertex Q, in each direction e:
138    *
139    *                                  e
140    *                           Q ----->----- R
141    *                      _,-'\__/
142    *                  _,-'       delta
143    *               P '
144    *
145    *                      r
146    *       cost    = delta          (we use r=3)
147    *           Q,e
148    *
149    *
150    * Calculation:
151    *
152    *      Let vector A = PQ
153    *                 B = QR
154    *
155    *                   -1   A . B
156    *      delta =  tan     -------
157    *                      | A x B |
158    *
159    *      which is always in the range 0..pi because the denominator
160    *      is nonnegative.  We add epsilon to |AxB| to avoid division
161    *      by zero.
162    *
163    *                     r
164    *      cost    = delta
165    *          Q,e
166    */
167
168 double line_bending_cost(const Vertices vertices) {
169   static const double axb_epsilon= 1e-6;
170   static const double exponent_r= 3;
171
172   int pi,e,qi,ri, k;
173   double  a[D3], b[D3], axb[D3];
174   double total_cost= 0;
175
176   FOR_EDGE(qi,e,ri) {
177     pi= EDGE_END2(qi,(e+3)%V6); if (pi<0) continue;
178
179     K a[k]= -vertices[pi][k] + vertices[qi][k];
180     K b[k]= -vertices[qi][k] + vertices[ri][k];
181
182     xprod(axb,a,b);
183
184     double delta= atan2(magnD(axb) + axb_epsilon, dotprod(a,b));
185     double cost= pow(delta,exponent_r);
186
187     if (!e && !(qi & ~XMASK))
188       cost *= 10;
189
190     total_cost += cost;
191   }
192   return total_cost;
193 }
194
195 /*---------- edge length variation ----------*/
196
197   /*
198    * Definition:
199    *
200    *    See the diagram above.
201    *                                r
202    *       cost    = ( |PQ| - |QR| )
203    *           Q,e
204    */
205
206 double edge_length_variation_cost(const Vertices vertices) {
207   double diff, cost= 0, exponent_r= 2;
208   int q, e,r, eback;
209
210   FOR_EDGE(q,e,r) {
211     eback= edge_reverse(q,e);
212     diff= edge_lengths[q][e] - edge_lengths[q][eback];
213     cost += pow(diff,exponent_r);
214   }
215   return cost;
216 }
217
218 /*---------- rim proximity cost ----------*/
219
220 static void find_nearest_oncircle(double oncircle[D3], const double p[D3]) {
221   /* By symmetry, nearest point on circle is the one with
222    * the same angle subtended at the z axis. */
223   oncircle[0]= p[0];
224   oncircle[1]= p[1];
225   oncircle[2]= 0;
226   double mult= 1.0/ magnD(oncircle);
227   oncircle[0] *= mult;
228   oncircle[1] *= mult;
229 }
230
231 double rim_proximity_cost(const Vertices vertices) {
232   double oncircle[3], cost=0;
233   int v;
234
235   FOR_VERTEX(v) {
236     int y= v >> YSHIFT;
237     int nominal_edge_distance= y <= Y/2 ? y : Y-1-y;
238     if (nominal_edge_distance==0) continue;
239
240     find_nearest_oncircle(oncircle, vertices[v]);
241
242     cost +=
243       vertex_mean_edge_lengths[v] *
244       (nominal_edge_distance*nominal_edge_distance) /
245       (hypotD2(vertices[v], oncircle) + 1e-6);
246   }
247   return cost;
248 }
249
250 /*---------- noncircular rim cost ----------*/
251
252 double noncircular_rim_cost(const Vertices vertices) {
253   int vy,vx,v;
254   double cost= 0.0;
255   double oncircle[3];
256
257   FOR_RIM_VERTEX(vy,vx,v) {
258     find_nearest_oncircle(oncircle, vertices[v]);
259
260     double d2= hypotD2(vertices[v], oncircle);
261     cost += d2*d2;
262   }
263   return cost;
264 }
265
266 /*---------- triangle bad normals cost ----------*/
267
268   /*
269    *
270    *                Q `-_
271    *              / |    `-_
272    *             /  |       `-.
273    *            /   |           S
274    *           /    |      _,-'
275    *          /     |  _,-'
276    *         /    , P '
277    *        /  ,-'
278    *       /,-'
279    *      /'
280    *     R
281    *
282    *  Let delta =  angle between two triangles' normals
283    *
284    *  Giving energy contribution:
285    *
286    *                                   2
287    *    E             =  F   .  delta
288    *     vd, edge PQ      vd
289    */
290
291 double edge_angle_cost(const Vertices vertices) {
292   double sq[D3], pq[D3], qr[D3], sqp[D3], pqr[D3], rs[D3];
293   double x[D3], y[D3];
294   int pi,e,qi,ri,si, k;
295 //  double our_epsilon=1e-6;
296   double total_cost= 0;
297   
298   FOR_EDGE(pi,e,qi) {
299 //    if (!(RIM_VERTEX_P(pi) || RIM_VERTEX_P(qi))) continue;
300
301     si= EDGE_END2(pi,(e+V6-1)%V6);  if (si<0) continue;
302     ri= EDGE_END2(pi,(e   +1)%V6);  if (ri<0) continue;
303
304     K {
305       sq[k]= vertices[si][k] - vertices[qi][k];
306       pq[k]= vertices[pi][k] - vertices[qi][k];
307       qr[k]= vertices[qi][k] - vertices[ri][k];
308     }
309     xprod(sqp, pq,sq);
310     xprod(pqr, pq,qr);
311
312     double dot= dotprod(sqp,pqr);
313     xprod(x,sqp,pqr);
314
315     K rs[k]= -vertices[ri][k] + vertices[si][k];
316     xprod(y, pq,rs);
317
318     double delta=
319       pow(atan2(magnD(x), dot), 2.0) * pow(magnD2(pq), 2.0) /
320       (pow(magnD(y), 0.3) + 1e-6);
321     double cost= pow(delta, 2.0);
322
323 //double cost= pow(magnD(spqxpqr), 3);
324 //assert(dot>=-1 && dot <=1);
325 //double cost= 1-dot;
326     total_cost += cost;
327   }
328   
329   return total_cost;
330 }
331
332 /*---------- small triangles cost ----------*/
333
334   /*
335    *
336    *                Q `-_
337    *              / |    `-_
338    *             /  |       `-.
339    *            /   |           S
340    *           /    |      _,-'
341    *          /     |  _,-'
342    *         /    , P '
343    *        /  ,-'
344    *       /,-'
345    *      /'
346    *     R
347    *
348    *  Let delta =  angle between two triangles' normals
349    *
350    *  Giving energy contribution:
351    *
352    *                                   2
353    *    E             =  F   .  delta
354    *     vd, edge PQ      vd
355    */
356
357 double small_triangles_cost(const Vertices vertices) {
358   double pq[D3], ps[D3];
359   double x[D3];
360   int pi,e,qi,si, k;
361 //  double our_epsilon=1e-6;
362   double total_cost= 0;
363   
364   FOR_EDGE(pi,e,qi) {
365 //    if (!(RIM_VERTEX_P(pi) || RIM_VERTEX_P(qi))) continue;
366
367     si= EDGE_END2(pi,(e+V6-1)%V6);  if (si<0) continue;
368
369     K {
370       pq[k]= vertices[qi][k] - vertices[pi][k];
371       ps[k]= vertices[si][k] - vertices[pi][k];
372     }
373     xprod(x, pq,ps);
374
375     double cost= 1/(magnD2(x) + 0.01);
376
377 //double cost= pow(magnD(spqxpqr), 3);
378 //assert(dot>=-1 && dot <=1);
379 //double cost= 1-dot;
380     total_cost += cost;
381   }
382   
383   return total_cost;
384 }