Bookstore

 Computer Math and Games in Pascal (preview) Lazarus Handbook

Author Topic: algorithms: which method or operator must i use??  (Read 9133 times)

BobDog

• Jr. Member
• Posts: 89
Re: algorithms: which method or operator must i use??
« Reply #15 on: May 19, 2021, 05:06:14 pm »

Another example, this time using case to guide the tangents round.
I assume that random (without randomize) always yields the same set of points.
On here the circles resemble the outline of the jeepers creepers 2 fellow's head.
Tested on 3.2.0 and 3.0.4, giving the same set of random numbers.
This is lazarus independent, I use a Dev Pascal ide.
Code: Pascal  [Select][+][-]
1.
2. program tangents;
3.
4. uses
5.   Math,
6.   graph;
7.
8. type
9.   v2 = object
10.   public
11.     x: double;
12.     y: double;
13.   end;
14.
15. type
16.   circles = object
17.     ctr: v2;
18.     r: integer;
19.     col: integer;
20.     idx:integer; // only to show numbers
21.   end;
22.
23. type
24.   linesegment = object
25.     s, f: v2;
26.   end;
27.
28. const
29.   same = 0;
30.   cross = 1;
31.   near = 2;
32.   far = 3;
33.
34. function dot(v1: v2; v3: v2): double;
35. var
36.   d1, d2, v1x, v1y, v2x, v2y: double;
37. begin
38.   d1 := Sqrt(v1.x * v1.x + v1.y * v1.y);
39.   d2 := Sqrt(v3.x * v3.x + v3.y * v3.y);
40.   v1x := v1.x / d1;
41.   v1y := v1.y / d1;
42.   v2x := v3.x / d2;
43.   v2y := v3.y / d2;
44.   exit(v1x * v2x + v1y * v2y);
45. end;
46.
47. function drawline(x: double; y: double; angle: double; length: double): v2;
48. var
49.   z: v2;
50. begin
51.   angle := angle * 0.0174532925199433;  // deg to rad
52.   z.x := x + length * Cos(angle);
53.   z.y := y - length * Sin(angle);
54.   exit(z);
55. end;
56.
57. function isleft(L: linesegment; p: v2): integer;
58. begin
59.   exit(-Sign((L.s.x - L.f.x) * (p.y - L.f.y) - (p.x - L.f.x) * (L.s.y - L.f.y)));
60. end;
61.
62. function segmentintersections(L1: linesegment; L2: linesegment): integer;
63. begin
64.   if isleft(L2, L1.s) = isleft(L2, L1.f) then
65.     exit(0);
66.   if isleft(L1, L2.s) = isleft(L1, L2.f) then
67.     exit(0);
68.   exit(1);
69. end;
70.
71. function distance(a: v2; b: v2): double;
72. begin
73.   exit(sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y)));
74. end;
75.
76. function tangent(C: circles; D: circles; flag: integer; position: integer; ref: v2): linesegment;
77. var
78.   z: array[0 .. 8] of linesegment;
79.   n, idx, a, b, i: integer;
80.   dst, dd: double;
81.   X, Y, Y2, v, v3: v2;
82.   CL, VL: linesegment;
83. label
84.   lbl;
85. begin
86.   n := 0;
87.   idx := 0;
88.   i := 0;
89.   dd := 0.0;
90.   dst := 0.0;
91.   if position = near then
92.     dst := 5000.0
93.   else
94.     dst := -5000.0;
95.
96.   for a := 0 to 360 do
97.   begin
98.     for b := 0 to 360 do
99.     begin
100.       v := drawline(C.ctr.x, C.ctr.y, a, C.r);
101.       v3 := drawline(D.ctr.x, D.ctr.y, b, D.r);
102.       X.x := (v.x - v3.x);
103.       X.y := (v.y - v3.y);
104.       Y.x := (c.ctr.x - v.x);
105.       Y.y := (c.ctr.y - v.y);
106.       Y2.x := (D.ctr.x - v3.x);
107.       Y2.y := (D.ctr.y - v3.y);
108.       if (Abs(dot(Y, X)) < 0.01) and (Abs(dot(Y2, X)) < 0.01) then
109.       begin
110.         CL.s := C.ctr;
111.         CL.f := D.ctr;
112.         VL.s := v;
113.         VL.f := v3;
114.         I := segmentintersections(CL, VL);
115.         if I = flag then
116.         begin
117.           dd := distance(v, ref);
118.           if (position = far) and (dst < dd) then
119.           begin
120.             dst := dd;
121.             idx := n;
122.           end;
123.
124.           if (position = near) and (dst >= dd) then
125.           begin
126.             dst := dd;
127.             idx := n;
128.           end;
129.
130.           z[n].s := v;
131.           z[n].f := v3;
132.           n += 1;
133.           goto lbl;
134.         end;//i
135.       end; //abs
136.     end; // a
137.     lbl: ;
138.   end; // b
139.   exit(z[idx]);
140. end;//function
141.
142. procedure setcases(var cflag:integer;val1:integer;var nflag:integer;val2:integer);
143. begin
144. cflag:=val1;
145. nflag:=val2;
146. end;
147.
148. var
149.   s: linesegment;
150.   gd, gm: smallint;
151.   c: array of circles;
152.   z, rad, Count, n, nflag, cflag: integer;
153.   ref, d: v2; //reference point as screen centre.
154.    id:string;
155. begin
156.   {==========  set up graph =========}
157.   gd := D8bit;
158.   gm := m800x600;
159.   InitGraph(gd, gm, ' ');
160.   if GraphResult <> grok then
161.     halt;
162.   SetTextStyle(SansSerifFont, HorizDir, 1);
163.   setbkcolor(15);
164.   cleardevice;
165.   {===================================}
166.
167.   ref.x := 400;
168.   ref.y := 300;
169.   Count := 0;
170.   cflag := cross;
171.   nflag := far;
172.   // set up some random circles around the centre
173.   for z := 0 to 360 do
174.     if z mod 20 = 0 then
175.     begin
176.       setlength(c, Count + 1);
177.       rad := 170 + random(120);
178.       d := drawline(ref.x, ref.y, z, rad);
179.       c[Count].ctr.x := d.x;
180.       c[Count].ctr.y := d.y;
181.       c[Count].r := 10+random(10);
182.       c[Count].col := 1+random(12);
183.       Count += 1;
184.     end;//for z
185.
186.   for n := 0 to high(c) - 1 do // draw the circles
187.   begin
188.     setfillstyle(solidfill, c[n].col);
189.     setcolor(c[n].col);
190.     fillellipse(trunc(c[n].ctr.x), trunc(c[n].ctr.y), c[n].r, c[n].r);
191.      str((n),id);
192.     OutTextXY(trunc(c[n].ctr.x), trunc(c[n].ctr.y)+20,id);
193.
194.   end;
195.
196.   for n := 0 to high(c) - 1 - 1 do    // tangents all the way round
197.   begin
198.
199.       case(n) of
200.       0: setcases(cflag,cross,nflag,near);
201.       1: setcases(cflag,cross,nflag,far);
202.       2:setcases(cflag,cross,nflag,near);
203.       3: setcases(cflag,same,nflag,far);
204.       4:setcases(cflag,cross,nflag,far);
205.       5:setcases(cflag,cross,nflag,near);
206.       6:setcases(cflag,same,nflag,far);
207.       7:setcases(cflag,cross,nflag,far);
208.       8:setcases(cflag,cross,nflag,near);
209.       9:setcases(cflag,cross,nflag,far);
210.       10:setcases(cflag,cross,nflag,near);
211.       11:setcases(cflag,same,nflag,far);
212.       15:setcases(cflag,cross,nflag,far);
213.       16:setcases(cflag,cross,nflag,near);
214.       end;
215.
216.     s := tangent(c[n], c[n + 1], cflag, nflag, ref);
217.     setcolor(green);
218.     line(trunc(s.s.x), trunc(s.s.y), trunc(s.f.x), trunc(s.f.y));
219.   end;
220.   s := tangent(c[n + 1], c[0], cross, far, ref);   //last pair
221.   line(trunc(s.s.x), trunc(s.s.y), trunc(s.f.x), trunc(s.f.y));
222.   writeln('DONE . . .  please press return');
224.   closegraph;
225.
226. end.

marcov

• Hero Member
• Posts: 9458
• FPC developer.
Re: algorithms: which method or operator must i use??
« Reply #16 on: May 19, 2021, 05:07:49 pm »
Call randomize on the first line of your main program.

BobDog

• Jr. Member
• Posts: 89
Re: algorithms: which method or operator must i use??
« Reply #17 on: May 19, 2021, 05:32:19 pm »

Thanks marcov.
But I want the random set to suit the tangents in case() of, so here I don't want to use randomize.

y.ivanov

• Full Member
• Posts: 197
Re: algorithms: which method or operator must i use??
« Reply #18 on: May 19, 2021, 05:53:42 pm »
*snip*
If I understand correctly, _2CExtTangents calculates the "external similitude center".
This will be undefined if the circles have the same radius. In your code the "eps" prevents a divide by zero. Perhaps this is the reason for the required +1.

This is still puzzling however, if I change +1 to -1 it fails, as it does using random radii.

To overcome the lack of an external homothetic center, I suggest using a simple right triangle logic instead of eps adjustments or like.

Sample code:
Code: Pascal  [Select][+][-]
1. procedure N_2CExtTangents(c0: TPoint; r0: Double; c1: TPoint; r1: Double;
2.   var t1, t2, t3, t4: TPoint; aMemo: TMemo = nil);
3. var
4.   dist, dx, dy: Double;
5.   yleg, xleg: LongInt;
6. begin
7.   // has external homothetic center?
8.   if not SameValue(r0, r1) then
9.     _2CExtTangents(c0, r0, c1, r1, t1, t2, t3, t4, aMemo)
10.   else
11.   begin
12.     yleg := 0; // assume x leg zero
13.     xleg := 0; // assume y leg zero
14.     dx := (c1.X - c0.X); // centers, delta x
15.     dy := (c1.Y - c0.Y); // centers, delta y
16.     if dy = 0 then // slope of pi/2, 3*pi/2, xleg zero
17.       yleg := Round(r0) // sin(pi/2) * r0 = r0
18.     else if dx = 0 then // slope of 0, pi, yleg zero
19.       xleg := Round(r0) // cos(0) * r0 = r0
20.     else // both legs nonzero
21.     begin
22.       dist := sqrt(dx * dx + dy * dy); // dist. between c0, c1
23.       yleg := Round(r0 * dx / dist); // r0 * sin(x)
24.       xleg := Round(r0 * dy / dist); // r0 * cos(x)
25.     end;
26.     // points of tangents due to symmetry
27.     t1 := Point(c0.X + xleg, c0.Y - yleg);
28.     t2 := Point(c0.X - xleg, c0.Y + yleg);
29.     t3 := Point(c1.X + xleg, c1.Y - yleg);
30.     t4 := Point(c1.X - xleg, c1.Y + yleg);
31.
32.     if Assigned(aMemo) then
33.       {show your debug info};
34.   end;
35. end;
36.

VTwin

• Hero Member
• Posts: 1037
• Former Turbo Pascal 3 user
Re: algorithms: which method or operator must i use??
« Reply #19 on: May 19, 2021, 08:16:56 pm »
To overcome the lack of an external homothetic center, I suggest using a simple right triangle logic instead of eps adjustments or like.

Nice idea.

With a slight modification:

Code: Pascal  [Select][+][-]
1.     t1 := Point(c0.X - xleg, c0.Y + yleg);
2.     t2 := Point(c0.X + xleg, c0.Y - yleg);
3.     t3 := Point(c1.X - xleg, c1.Y + yleg);
4.     t4 := Point(c1.X + xleg, c1.Y - yleg);

this works except for the "crossed" pulleys in the "Preset" configuration. All "+1" are removed.
« Last Edit: May 19, 2021, 08:31:27 pm by VTwin »
“Talk is cheap. Show me the code.” -Linus Torvalds

Free Pascal Compiler 3.2.0
macOS 10.13.6: Lazarus 2.0.12 (64 bit Cocoa)
Ubuntu 18.04.3: Lazarus 2.0.12 (64 bit on VBox)
Windows 7 Pro SP1: Lazarus 2.0.12 (64 bit on VBox)

y.ivanov

• Full Member
• Posts: 197
Re: algorithms: which method or operator must i use??
« Reply #20 on: May 21, 2021, 01:52:38 am »
*snip*
With a slight modification:

Code: Pascal  [Select][+][-]
1.     t1 := Point(c0.X - xleg, c0.Y + yleg);
2.     t2 := Point(c0.X + xleg, c0.Y - yleg);
3.     t3 := Point(c1.X - xleg, c1.Y + yleg);
4.     t4 := Point(c1.X + xleg, c1.Y - yleg);

this works except for the "crossed" pulleys in the "Preset" configuration. All "+1" are removed.

The correction is just for the "preset", the mistake into my suggestion was the sign omission. Furthermore, the original _2CExtTangents() swaps two circles when r0<r1. That disturbs the LTR order of the tangents. The updated snippet is:
Code: Pascal  [Select][+][-]
1. procedure N_2CExtTangents(c0: TPoint; r0: Double; c1: TPoint; r1: Double;
2.   var t1, t2, t3, t4: TPoint; aMemo: TMemo = nil);
3. var
4.   dist, dx, dy: Double;
5.   yleg, xleg: LongInt;
6. begin
7.   // has external homothetic center?
8.   if not SameValue(r0, r1) then
9.   begin
10.     if r0 < r1 then // keep the ltr order!
11.       _2CExtTangents(c0, r0, c1, r1, t2, t1, t4, t3, aMemo) else
12.       _2CExtTangents(c0, r0, c1, r1, t1, t2, t3, t4, aMemo)
13.   end
14.   else
15.   begin
16.     yleg := 0; // assume x leg zero
17.     xleg := 0; // assume y leg zero
18.     dx := (c1.X - c0.X); // centers, delta x
19.     dy := (c1.Y - c0.Y); // centers, delta y
20.     if dy = 0 then // slope of pi/2, 3*pi/2, xleg zero
21.       yleg := Round(Sign(dx) * r0) // sin(pi/2) * r0 = r0
22.     else if dx = 0 then // slope of 0, pi, yleg zero
23.       xleg := Round(Sign(dy) * r0) // cos(0) * r0 = r0
24.     else // both legs nonzero
25.     begin
26.       dist := sqrt(dx * dx + dy * dy); // dist. between c0, c1
27.       yleg := Round(r0 * dx / dist); // r0 * sin(x)
28.       xleg := Round(r0 * dy / dist); // r0 * cos(x)
29.     end;
30.     // points of tng due to symmetry
31.     t1 := Point(c0.X + xleg, c0.Y - yleg);
32.     t2 := Point(c0.X - xleg, c0.Y + yleg);
33.     t3 := Point(c1.X + xleg, c1.Y - yleg);
34.     t4 := Point(c1.X - xleg, c1.Y + yleg);
35.
36.     if Assigned(aMemo) then
37.       {show debug info};
38.
39.   end;
40. end;
41.

I'm also attaching a replacement for the TfrmMain.cmdPulleysClick() method. It works according to the idea into my first reply on the topic. With the updated N_2CExtTangents() it behaves fairly well with both same and different radii.
It does not give the best solution (with the most curves) but is a straightforward one with a view just over 2 pulleys at a time and no backtracking.

Code: Pascal  [Select][+][-]
1. procedure TfrmMain.cmdPulleysClick(Sender: TObject);
2. var
3.   H, W, nCount, k, n: integer;
4.   c, r, g, b: TColor;
5.
6.   cA, cB: TPoint; // centre
7.   rA, rB: Integer; // radii
8.   I, PI, J, M: Integer;
9.   LSide, ASide: Boolean;
10.   tng: array of array of TPoint;
11.
12.   function OnSegment(p, q, r: TPoint): Boolean;
13.   begin
14.     Result := (q.x <= Max(p.x, r.x)) and (q.x >= Min(p.x, r.x)) and
15.       (q.y <= Max(p.y, r.y)) and (q.y >= Min(p.y, r.y));
16.   end;
17.
18.   function Orientation(p, q, r: TPoint): Integer; // sideOfPoint()?
19.   var
20.     V: LongInt;
21.   begin
22.     V := (q.y - p.y) * (r.x - q.x) - (q.x - p.x) * (r.y - q.y);
23.     if V = 0 then
24.       Result := 0 // on a line
25.     else if V > 0 then
26.       Result := 1 // clock wise
27.     else
28.       Result := -1; // counterclock wise
29.   end;
30.
31.   function Intersects(p1, q1, p2, q2: TPoint): Boolean;
32.   var
33.     o1, o2, o3, o4: Integer;
34.   begin
35.     o1 := Orientation(p1, q1, p2);
36.     o2 := Orientation(p1, q1, q2);
37.     o3 := Orientation(p2, q2, p1);
38.     o4 := Orientation(p2, q2, q1);
39.
40.     if (o1 <> o2) and (o3 <> o4) then
41.       Result := True
42.     else if
43.       ((o1 = 0) and onSegment(p1, p2, q1)) or
44.       ((o2 = 0) and onSegment(p1, q2, q1)) or
45.       ((o3 = 0) and onSegment(p2, p1, q2)) or
46.       ((o4 = 0) and onSegment(p2, q1, q2))
47.     then
48.       Result := True
49.     else
50.       Result := False;
51.   end;
52.
53.   procedure LoadCircle(I: Integer; out C: TPoint; out R: Integer);
54.   var IM: Integer;
55.   begin
56.     IM := I mod nPoints;
57.     C := Point(Round(points[IM].X), Round(points[IM].Y));
58.     R := Round(yarichaps[IM]);
59.   end;
60.
61.   procedure ShowPulley(ACenter: TPoint; ARadii: Integer; ALabel: String; WithColor: TColor);
62.   begin
63.     imgInput.Canvas.Pen.Color := WithColor;
64.     imgInput.Canvas.Pen.Width := 2;
65.     imgInput.Canvas.Ellipse(ACenter.X - ARadii, ACenter.Y - ARadii,
66.       ACenter.X + ARadii, ACenter.Y + ARadii);
67.     imgInput.Canvas.TextOut(ACenter.X - 5, ACenter.Y - 5, ALabel);
68.   end;
69.
70.   procedure SelTangent(I: Integer; ATangent: Integer);
71.   var
72.     J: Integer;
73.   begin
74.     for J := 0 to 3 do
75.       if J <> ATangent then
76.         tng[I, J] := Point(-1, -1)
77.       else if tng[I, J].X <> -1 then
78.       begin
79.         pout[I * 2] := tng[I, J];
80.         pout[I * 2 + 1] := tng[I, J + 4];
81.       end;
82.   end;
83.
84. begin
85.   cmdResetClick(Sender);
86.   lastAction := cmdPulleysClick;
87.   Screen.Cursor := crHourGlass;
88.   STOP := False;
89.   memoDebug.Clear;
90.   H := 500; W := 500;
91.   imgInput.Picture.Bitmap.Height := H;
92.   imgInput.Picture.Bitmap.Width := W;
93.   imgInput.Canvas.Brush.Color := clWhite;
94.   imgInput.Canvas.FillRect(imgInput.ClientRect);
95.   imgInput.Canvas.Pen.Mode := pmCopy;
96.   imgInput.Canvas.Brush.Style := bsClear;
97.
98.   //==================================
99.   SetLength(tng, nPoints, 8);
100.
101.   // Calculate tangents
102.   PI := Pred(nPoints);
103.   LoadCircle(PI, cA, rA);
104.   for I := 0 to Pred(nPoints) do
105.   begin
106.     LoadCircle(I, cB, rB);
107.     N_2CExtTangents(cA, rA, cB, rB, tng[I, 0], tng[I, 1], tng[I, 4], tng[I, 5]);
108.     _2CIntTangents(cA, rA, cB, rB, tng[I, 2], tng[I, 3], tng[I, 6], tng[I, 7]);
109.     ShowPulley(cB, rB, I.ToString, clLtGray);
110.     cA := cB;  rA := rB;
111.     PI := I;
112.   end;
113.
114.   // Remove intersecting tangents
115.   for I := Pred(nPoints) downto 0 do
116.   begin
117.     PI := Pred(I); if PI < 0 then PI := Pred(nPoints);
118.     for J := 0 to 3 do
119.     begin
120.       M := 0;
121.       while (M < 4) and not
122.         Intersects(tng[I, J], tng[I, J + 4], tng[PI, M], tng[PI, M + 4])
123.       do
124.         Inc(M);
125.       if M < 4 then
126.         tng[I, J] := Point(-1, -1); // mark deleted
127.     end;
128.   end;
129.
130.   // Select a tangent for drawing
131.   ASide := False;
132.   PI := Pred(nPoints);
133.   for I := 0 to Pred(nPoints) do
134.   begin
135.     J := 0; M := 0;
136.     J := J + IfThen(tng[I, 0].X <> -1, 1); // left ext
137.     J := J + IfThen(tng[I, 2].X <> -1, 1); // left int
138.     M := M + IfThen(tng[I, 1].X <> -1, 1); // right ext
139.     M := M + IfThen(tng[I, 3].X <> -1, 1); // right int
140.
141.     if J = M then
142.     begin
143.       LSide := ASide;
144.       ASide := not ASide;
145.     end
146.     else
147.       LSide := J > M;
148.
149.     if LSide then  // left side
150.     begin
151.       if (tng[PI, 3].X <> -1) then
152.         SelTangent(PI, 3) else
153.         SelTangent(PI, 0);
154.       tng[I, 1] := Point(-1, -1);
155.       tng[I, 3] := Point(-1, -1);
156.     end
157.     else // right side
158.     begin
159.       if (tng[PI, 2].X <> -1) then
160.         SelTangent(PI, 2) else
161.         SelTangent(PI, 1);
162.       tng[I, 0] := Point(-1, -1);
163.       tng[I, 2] := Point(-1, -1);
164.     end;
165.     PI := I;
166.   end;
167.   nCount := nPoints;
168.   //=============================
169.
170.   //////////////////
171.   imgInput.Canvas.Pen.Color := clBlue;
172.   imgInput.Canvas.Pen.Width := 2;
173.   //if False then
174.   for k:=0 to nCount-1 do
175.   begin
176.     r := random(255); g := random(255); b := random(255);
177.     if (r < 32) then r := 32; if (r > 240) then r := 240;
178.     if (g < 32) then g := 32; if (g > 240) then g := 240;
179.     if (b < 32) then b := 32; if (b > 240) then b := 240;
180.     c := RGB(r, g, b);
181.     imgInput.Canvas.Pen.Color := c;
182.     n := 2*k + 0;
183.     imgInput.Canvas.Ellipse(pout[n].X-2, pout[n].Y-2, pout[n].X+2, pout[n].Y+2);
184.     imgInput.Canvas.MoveTo(pout[n].X, pout[n].Y);
185.     n := 2*k + 1;
186.     imgInput.Canvas.LineTo(pout[n].X, pout[n].Y);
187.     imgInput.Canvas.Ellipse(pout[n].X-2, pout[n].Y-2, pout[n].X+2, pout[n].Y+2);
188.   end;
189.   if chkFileOutput.Checked then
190.     imgInput.Picture.Bitmap.SaveToFile(ExtractFilePath(Application.ExeName) + 'zzzPulleys-' + IntToHex(RandSeed, 8) + '.bmp');
191.   Screen.Cursor := crDefault;
192.   STOP := True;
193. end;
194.
« Last Edit: May 21, 2021, 03:16:46 pm by y.ivanov »