1/*************************************************************************
2 * Copyright (c) 2011 AT&T Intellectual Property
3 * All rights reserved. This program and the accompanying materials
4 * are made available under the terms of the Eclipse Public License v1.0
5 * which accompanies this distribution, and is available at
6 * http://www.eclipse.org/legal/epl-v10.html
7 *
8 * Contributors: See CVS logs. Details at http://www.graphviz.org/
9 *************************************************************************/
10
11#include <math.h>
12#include <stdlib.h>
13#include "types.h"
14#include "globals.h"
15#include "general.h"
16#include "SparseMatrix.h"
17#include "edge_bundling.h"
18#include "ink.h"
19
20double ink_count;
21
22static point_t addPoint (point_t a, point_t b)
23{
24 a.x += b.x;
25 a.y += b.y;
26 return a;
27}
28
29static point_t subPoint (point_t a, point_t b)
30{
31 a.x -= b.x;
32 a.y -= b.y;
33 return a;
34}
35
36static point_t scalePoint (point_t a, double d)
37{
38 a.x *= d;
39 a.y *= d;
40 return a;
41}
42
43static double dotPoint(point_t a, point_t b){
44 return a.x*b.x + a.y*b.y;
45}
46
47
48point_t Origin;
49
50/* sumLengths:
51 */
52static double sumLengths_avoid_bad_angle(point_t* points, int npoints, point_t end, point_t meeting, real angle_param)
53{
54 /* avoid sharp turns, we want cos_theta to be as close to -1 as possible */
55 int i;
56 double len0, len, sum = 0;
57 double diff_x, diff_y, diff_x0, diff_y0;
58 double cos_theta, cos_max = -10;
59
60 diff_x0 = end.x-meeting.x;
61 diff_y0 = end.y-meeting.y;
62 len0 = sum = sqrt(diff_x0*diff_x0+diff_y0*diff_y0);
63
64 // distance form each of 'points' till 'meeting'
65 for (i=0; i<npoints; i++) {
66 diff_x = points[i].x-meeting.x;
67 diff_y = points[i].y-meeting.y;
68 len = sqrt(diff_x*diff_x+diff_y*diff_y);
69 sum += len;
70 cos_theta = (diff_x0*diff_x + diff_y0*diff_y)/MAX((len*len0),0.00001);
71 cos_max = MAX(cos_max, cos_theta);
72 }
73
74 // distance of single line from 'meeting' to 'end'
75 return sum*(cos_max + angle_param);/* straight line gives angle_param - 1, turning angle of 180 degree gives angle_param + 1 */
76}
77static double sumLengths(point_t* points, int npoints, point_t end, point_t meeting)
78{
79 int i;
80 double sum = 0;
81 double diff_x, diff_y;
82
83 // distance form each of 'points' till 'meeting'
84 for (i=0; i<npoints; i++) {
85 diff_x = points[i].x-meeting.x;
86 diff_y = points[i].y-meeting.y;
87 sum += sqrt(diff_x*diff_x+diff_y*diff_y);
88 }
89 // distance of single line from 'meeting' to 'end'
90 diff_x = end.x-meeting.x;
91 diff_y = end.y-meeting.y;
92 sum += sqrt(diff_x*diff_x+diff_y*diff_y);
93 return sum;
94}
95
96/* bestInk:
97 */
98static double bestInk(point_t* points, int npoints, point_t begin, point_t end, double prec, point_t *meet, real angle_param)
99{
100 point_t first, second, third, fourth, diff, meeting;
101 double value1, value2, value3, value4;
102
103 first = begin;
104 fourth = end;
105
106 do {
107 diff = subPoint(fourth,first);
108 second = addPoint(first,scalePoint(diff,1.0/3.0));
109 third = addPoint(first,scalePoint(diff,2.0/3.0));
110
111 if (angle_param < 1){
112 value1 = sumLengths(points, npoints, end, first);
113 value2 = sumLengths(points, npoints, end, second);
114 value3 = sumLengths(points, npoints, end, third);
115 value4 = sumLengths(points, npoints, end, fourth);
116 } else {
117 value1 = sumLengths_avoid_bad_angle(points, npoints, end, first, angle_param);
118 value2 = sumLengths_avoid_bad_angle(points, npoints, end, second, angle_param);
119 value3 = sumLengths_avoid_bad_angle(points, npoints, end, third, angle_param);
120 value4 = sumLengths_avoid_bad_angle(points, npoints, end, fourth, angle_param);
121 }
122
123 if (value1<value2) {
124 if (value1<value3) {
125 if (value1<value4) {
126 // first is smallest
127 fourth = second;
128 }
129 else {
130 // fourth is smallest
131 first = third;
132 }
133 }
134 else {
135 if (value3<value4) {
136 // third is smallest
137 first = second;
138 }
139 else {
140 // fourth is smallest
141 first = third;
142 }
143 }
144 }
145 else {
146 if (value2<value3) {
147 if (value2<value4) {
148 // second is smallest
149 fourth = third;
150 }
151 else {
152 // fourth is smallest
153 first = third;
154 }
155 }
156 else {
157 if (value3<value4) {
158 // third is smallest
159 first = second;
160 }
161 else {
162 // fourth is smallest
163 first = third;
164 }
165 }
166 }
167 } while (fabs(value1-value4)/(MIN(value1,value4)+1e-10)>prec && dotPoint(diff, diff) > 1.e-20);
168
169 meeting = scalePoint(addPoint(first,fourth),0.5);
170 *meet = meeting;
171
172 return sumLengths(points, npoints, end, meeting);
173
174}
175
176static double project_to_line(point_t pt, point_t left, point_t right, real angle){
177 /* pt
178 ^ ^
179 . \ \
180 . \ \
181 d . a\ \
182 . \ \
183 . \ \
184 . c \ alpha \ b
185 .<------left:0 ----------------------------> right:1. Find the projection of pt on the left--right line. If the turning angle is small,
186 | |
187 |<-------f---------
188 we should get a negative number. Let a := left->pt, b := left->right, then we are calculating:
189 c = |a| cos(a,b)/|b| b
190 d = a - c
191 f = -ctan(alpha)*|d|/|b| b
192 and the project of alpha degree on the left->right line is
193 c-f = |a| cos(a,b)/|b| b - -ctan(alpha)*|d|/|b| b
194 = (|a| a.b/(|a||b|) + ctan(alpha)|a-c|)/|b| b
195 = (a.b/|b| + ctan(alpha)|a-c|)/|b| b
196 the dimentionless projection is:
197 a.b/|b|^2 + ctan(alpha)|a-c|/|b|
198 = a.b/|b|^2 + ctan(alpha)|d|/|b|
199 */
200
201
202 point_t b, a;
203 real bnorm, dnorm;
204 real alpha, ccord;
205
206 if (angle <=0 || angle >= M_PI) return 2;/* return outside of the interval which should be handled as a sign of infeasible turning angle */
207 alpha = angle;
208
209 assert(alpha > 0 && alpha < M_PI);
210 b = subPoint(right, left);
211 a = subPoint(pt, left);
212 bnorm = MAX(1.e-10, dotPoint(b, b));
213 ccord = dotPoint(b, a)/bnorm;
214 dnorm = dotPoint(a,a)/bnorm - ccord*ccord;
215 if (alpha == M_PI/2){
216 return ccord;
217 } else {
218 // assert(dnorm >= MIN(-1.e-5, -1.e-5*bnorm));
219 return ccord + sqrt(MAX(0, dnorm))/tan(alpha);
220 }
221}
222
223
224
225
226
227
228
229
230
231/* ink:
232 * Compute minimal ink used the input edges are bundled.
233 * Assumes tails all occur on one side and heads on the other.
234 */
235double ink(pedge* edges, int numEdges, int *pick, double *ink0, point_t *meet1, point_t *meet2, real angle_param, real angle)
236{
237 int i;
238 point_t begin, end, mid, diff;
239 pedge e;
240 real *x;
241 point_t* sources = N_NEW(numEdges, point_t);
242 point_t* targets = N_NEW(numEdges, point_t);
243 double inkUsed;
244 //double eps = 1.0e-2;
245 double eps = 1.0e-2;
246 double cend = 0, cbegin = 0;
247 double wgt = 0;
248
249 // fprintf(stderr,"in ink code ========\n");
250 ink_count += numEdges;
251
252 *ink0 = 0;
253
254 /* canonicalize so that edges 1,2,3 and 3,2,1 gives the same optimal ink */
255 if (pick) vector_sort_int(numEdges, pick, TRUE);
256
257 begin = end = Origin;
258 for (i = 0; i < numEdges; i++) {
259 if (pick) {
260 e = edges[pick[i]];
261 } else {
262 e = edges[i];
263 }
264 x = e->x;
265 sources[i].x = x[0];
266 sources[i].y = x[1];
267 targets[i].x = x[e->dim*e->npoints - e->dim];
268 targets[i].y = x[e->dim*e->npoints - e->dim + 1];
269 (*ink0) += sqrt((sources[i].x - targets[i].x)*(sources[i].x - targets[i].x) + (sources[i].y - targets[i].y)*(sources[i].y - targets[i].y));
270 begin = addPoint (begin, scalePoint(sources[i], e->wgt));
271 end = addPoint (end, scalePoint(targets[i], e->wgt));
272 wgt += e->wgt;
273 //fprintf(stderr,"source={%f,%f}, target = {%f,%f}\n",sources[i].x, sources[i].y,
274 //targets[i].x, targets[i].y);
275 }
276
277 begin = scalePoint (begin, 1.0/wgt);
278 end = scalePoint (end, 1.0/wgt);
279
280
281 if (numEdges == 1){
282 *meet1 = begin;
283 *meet2 = end;
284 //fprintf(stderr,"ink used = %f\n",*ink0);
285 free (sources);
286 free (targets);
287 return *ink0;
288 }
289
290 /* shift the begin and end point to avoid sharp turns */
291 for (i = 0; i < numEdges; i++) {
292 if (pick) {
293 e = edges[pick[i]];
294 } else {
295 e = edges[i];
296 }
297 x = e->x;
298 sources[i].x = x[0];
299 sources[i].y = x[1];
300 targets[i].x = x[e->dim*e->npoints - e->dim];
301 targets[i].y = x[e->dim*e->npoints - e->dim + 1];
302 /* begin(1) ----------- mid(0) */
303 if (i == 0){
304 cbegin = project_to_line(sources[i], begin, end, angle);
305 cend = project_to_line(targets[i], end, begin, angle);
306 } else {
307 cbegin = MAX(cbegin, project_to_line(sources[i], begin, end, angle));
308 cend = MAX(cend, project_to_line(targets[i], end, begin, angle));
309 }
310 }
311
312 if (angle > 0 && angle < M_PI){
313 if (cbegin + cend > 1 || cbegin > 1 || cend > 1){
314 /* no point can be found that satisfies the angular constraints, so we give up and set ink to a large value */
315 if (Verbose && 0) fprintf(stderr,"no point satisfying any angle constraints can be found. cbeg=%f cend=%f\n",cbegin,cend);
316 inkUsed = 1000*(*ink0);
317 *meet1 = *meet2 = mid;
318 free (sources);
319 free (targets);
320 return inkUsed;
321 }
322 /* make sure the turning angle is no more than alpha degree */
323 cbegin = MAX(0, cbegin);/* make sure the new adjusted point is with in [begin,end] internal */
324 diff = subPoint(end, begin);
325 begin = addPoint(begin, scalePoint(diff, cbegin));
326
327 cend = MAX(0, cend);/* make sure the new adjusted point is with in [end,begin] internal */
328 end = subPoint(end, scalePoint(diff, cend));
329 }
330 mid = scalePoint (addPoint(begin,end),0.5);
331
332 inkUsed = (bestInk (sources, numEdges, begin, mid, eps, meet1, angle_param)
333 + bestInk (targets, numEdges, end, mid, eps, meet2, angle_param));
334 //fprintf(stderr,"beg={%f,%f}, meet1={%f,%f}, meet2={%f,%f}, mid={%f,%f}, end={%f,%f}\n",begin.x, begin.y, meet1->x, meet1->y, meet2->x, meet2->y,
335 //mid.x, mid.y, end.x, end.y);
336
337 //fprintf(stderr,"ink used = %f\n",inkUsed);
338 // fprintf(stderr,"{cb,ce}={%f, %f} end={%f,%f}, meet={%f,%f}, mid={%f,%f}\n",cbegin, cend, end.x, end.y, meet2->x, meet2->y, mid.x, mid.y);
339 free (sources);
340 free (targets);
341 return inkUsed;
342}
343
344
345double ink1(pedge e){
346
347
348 real *x, xx, yy;
349
350 real ink0 = 0;
351
352 x = e->x;
353 xx = x[0] - x[e->dim*e->npoints - e->dim];
354 yy = x[1] - x[e->dim*e->npoints - e->dim + 1];
355 ink0 += sqrt(xx*xx + yy*yy);
356 return ink0;
357}
358