← Back to branch summary

~ubuntu-branches/ubuntu/karmic/scilab/karmic

« back to all changes in this revision

Viewing changes to routines/metanet/hullcvex.f

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2002-03-21 16:57:43 UTC
  • Revision ID: james.westby@ubuntu.com-20020321165743-e9mv12c1tb1plztg
Tags: upstream-2.6
Import upstream version 2.6

Show diffs side-by-side

added added

removed removed

Lines of Context:
routines/metanet/hullcvex.f
 
1
c     from algorithm 523, acm.
 
2
      subroutine hullcvex(n,nn,xx,nhull,iwork,in,ih,il)
 
3
      implicit doubleprecision (a-h,o-z)
 
4
      dimension xx(2,n),in(n),ih(n),iwork(nn),il(nn)
 
5
      ncount=0
 
6
      nhull=0
 
7
      do 1,i=1,nn
 
8
         iwork(i)=0
 
9
 1    continue
 
10
c ncount is total number of points passed to split
 
11
      n1=n+1
 
12
      do 2 i=1,n
 
13
         j=n1-i
 
14
         in(j)=i
 
15
 2    continue
 
16
c array in contains indices 1-n in reverse order
 
17
      do 10 m=4,n
 
18
         call convex(n,xx,m,in,iwork,iwork(n+1),ih,nhull,il)
 
19
         ik=il(1)
 
20
         do 6 i=1,nhull
 
21
            j=ih(ik)
 
22
            iwork(i)=j
 
23
            ik=il(ik)
 
24
 6       continue
 
25
 10   continue
 
26
      end
 
27
      subroutine husplit(n,x,m,in,ii,jj,s,iabv,na,maxa,ibel,
 
28
     1  nb,maxb)
 
29
c this subroutine takes the m points of array x whose subscripts are in
 
30
c  array in
 
31
c and partitions them by the line joining the two points in array x whose  
 
32
c subscripts are ii and jj. the subscripts of the points above the line
 
33
c are put into array iabv, and the subscripts of the points below are put
 
34
c into array ibel. na and nb are,respectively, the number of points above 
 
35
c the line and the number below. maxa and maxb are the subscripts for array
 
36
c x of the point furthest above the line and the point furthest below, 
 
37
c respectively. if either subset is null the corresponding 
 
38
c subscript (maxa or maxb) is set to zero
 
39
c input
 
40
c n    integer           total number of data points
 
41
c x    real array (2,n)  (x,y) co-ordinates of the data
 
42
c m    integer           number of points in input subset
 
43
c in   integer array (m) subscripts for array x of the 
 
44
c                        points in the input subset
 
45
c ii   integer           subscript for array x of one point
 
46
c                        on the partitioning line
 
47
c jj   integer           subscript for array x of another
 
48
c                        point on the partitioning line
 
49
c s    integer           switch to determine output. refer
 
50
c                        to comments below
 
51
c output
 
52
c iabv integer array (m) subscripts for array x of the
 
53
c                        points above the partitioning line
 
54
c na   integer           number of elements in iabv
 
55
c maxa integer           subscript for array x of point
 
56
c                        furthest above the line. set to
 
57
c                        zero if na is zero
 
58
c ibel integer array (m) subscripts for array x of the
 
59
c                        points below the partitioning line
 
60
c nb   integer           number of elements in ibel
 
61
c maxb integer           subscript for array x of point
 
62
c                        furthest below the line. set to
 
63
c                        zero if nb is zero
 
64
      implicit doubleprecision (a-h,o-z)
 
65
      dimension x(2,n)
 
66
      dimension in(m),iabv(m),ibel(m)
 
67
      integer s
 
68
c if s = 2 dont save ibel,nb,maxb.
 
69
c if s =-2 dont save iabv,na,maxa.
 
70
c otherwise save everything
 
71
c if s is positive the array being partitioned is above the initial 
 
72
c partitioning line. if it is negative, then the set of points is below.
 
73
      logical t
 
74
      t=.false.
 
75
c check to see if the line is vertical
 
76
      if(x(1,jj).ne.x(1,ii))goto 1
 
77
      xt=x(1,ii)
 
78
      dir=sign(1.D0,x(2,jj)-x(2,ii))*sign(1.D0,dble(s))
 
79
      t=.true.
 
80
      goto 2
 
81
1     a=(x(2,jj)-x(2,ii))/(x(1,jj)-x(1,ii))
 
82
      b=x(2,ii)-a*x(1,ii)
 
83
2     up=0.
 
84
      na=0
 
85
      maxa=0
 
86
      down=0.
 
87
      nb=0
 
88
      maxb=0
 
89
      do 6 i=1,m
 
90
        is=in(i)
 
91
        if(t)goto 3
 
92
        z=x(2,is)-a*x(1,is)-b
 
93
        goto 4
 
94
3       z=dir*(x(1,is)-xt)
 
95
4       if(z.le.0.)goto 5
 
96
c the point is above the line
 
97
        if(s.eq.-2)goto 6
 
98
        na=na+1
 
99
        iabv(na)=is
 
100
        if(z.lt.up)goto 6
 
101
        up=z
 
102
        maxa=na
 
103
        goto 6
 
104
5       if(s.eq.2)goto 6
 
105
        if(z.ge.0.)goto 6
 
106
c the point is below the line
 
107
        nb=nb+1
 
108
        ibel(nb)=is
 
109
        if(z.gt.down)goto 6
 
110
        down=z
 
111
        maxb=nb
 
112
6     continue
 
113
      return
 
114
      end
 
115
      subroutine convex(n,x,m,in,ia,ib,ih,nh,il)
 
116
c this subroutine determines which of the m points of array x whose subscripts
 
117
c are in array in are vertices of the minimum area convex polygon containing
 
118
c the m points. the subscripts of the vertices are placed in array ih in the
 
119
c order they are found. nh is the number of elements in array ih and  array il.
 
120
c array il is a linked list giving the order of the elements of array ih
 
121
c in a counter clockwise direction. this algorithm corresponds to a
 
122
c preorder traversal of a certain binary tree. each vertex of the binary tree
 
123
c represents a subset of the m points. at each step the subset of points 
 
124
c corresponding to the current vertex of the tree is partitioned by a line
 
125
c joining two vertices of the convex polygon. the left son vertex in the binary
 
126
c tree represents the subset of points above the partitioning line and the 
 
127
c right son vertex, the subset below the line. the leaves of the tree represent
 
128
c either null subsets or subsets inside a triangle whose vertices coincide 
 
129
c with vertices of the convex polygon.
 
130
c formal parameters
 
131
c input
 
132
c n  integer           total number of data points
 
133
c x  real array (2,n)  (x,y) co-ordinates of the data
 
134
c m  integer           number of points in the input subset
 
135
c in integer array (m) subscripts for array x of the points
 
136
c                      in the input subset
 
137
c work area
 
138
c ia integer array (m) subscripts for array x of left son
 
139
c                      subsets. see comments after dimension
 
140
c                      statements
 
141
c ib integer array (m) subscripts for array x of right son
 
142
c                      subsets
 
143
c output
 
144
c ih integer array (m) subscripts for array x of the
 
145
c                      vertices of the convex hull
 
146
c nh integer           number of elements in array ih and
 
147
c                      array il. same as number of vertices
 
148
c                      of the convex polygon
 
149
c il integer array (m) a linked list giving in order in a
 
150
c                      counter-clockwise direction the
 
151
c                      elements of array ih
 
152
      implicit doubleprecision (a-h,o-z)
 
153
      dimension x(2,n)
 
154
      dimension in(m),ia(m),ib(m),ih(m),il(m)
 
155
c the upper end of array ia is used to store temporarily the sizes of the 
 
156
c subsets which correspond to right son vertices, while traversing down 
 
157
c the left sons when on the left half of the tree, and to store the sizes 
 
158
c of the left sons while traversing the right sons(down the right half)
 
159
      logical maxe,mine
 
160
      if(m.eq.1)goto 22
 
161
      il(1)=2
 
162
      il(2)=1
 
163
      kn=in(1)
 
164
      kx=in(2)
 
165
      if(m.eq.2)goto 21
 
166
      mp1=m+1
 
167
      min=1
 
168
      mx=1
 
169
      kx=in(1)
 
170
      maxe=.false.
 
171
      mine=.false.
 
172
c find two vertices of the convex hull for the initial partition
 
173
      do 6 i=2,m
 
174
        j=in(i)
 
175
        if(x(1,j)-x(1,kx))3,1,2
 
176
1       maxe=.true.
 
177
        goto 3
 
178
2       maxe=.false.
 
179
        mx=i
 
180
        kx=j
 
181
3       if(x(1,j)-x(1,kn))5,4,6
 
182
4       mine=.true.
 
183
        goto 6
 
184
5       mine=.false.
 
185
        min=i
 
186
        kn=j
 
187
6     continue
 
188
c if the max and min are equal, all m points lie on a vertical line
 
189
      if(kx.eq.kn)goto 18
 
190
c if maxe (or mine) has the value true there are several
 
191
c maxima (or minima) with equal first coordinates
 
192
      if(maxe.or.mine)goto 23
 
193
7     ih(1)=kx
 
194
      ih(2)=kn
 
195
      nh=3
 
196
      inh=1
 
197
      nib=1
 
198
      ma=m
 
199
      in(mx)=in(m)
 
200
      in(m)=kx
 
201
      mm=m-2
 
202
      if(min.eq.m)min=mx
 
203
      in(min)=in(m-1)
 
204
      in(m-1)=kn
 
205
c begin by partitioning the root of the tree
 
206
      call husplit(n,x,mm,in,ih(1),ih(2),0,ia,mb,mxa,ib,ia(ma),
 
207
     1  mxbb)
 
208
c first traverse the left half of the tree start with the left son
 
209
8     nib=nib+ia(ma)
 
210
      ma=ma-1
 
211
9     if(mxa.eq.0)goto 11
 
212
      il(nh)=il(inh)
 
213
      il(inh)=nh
 
214
      ih(nh)=ia(mxa)
 
215
      ia(mxa)=ia(mb)
 
216
      mb=mb-1
 
217
      nh=nh+1
 
218
      if(mb.eq.0)goto 10
 
219
      ilinh=il(inh)
 
220
      call husplit(n,x,mb,ia,ih(inh),ih(ilinh),1,ia,mbb,mxa,
 
221
     1  ib(nib),ia(ma),mxb)
 
222
      mb=mbb
 
223
      goto 8
 
224
c then the right son
 
225
10    inh=il(inh)
 
226
11    inh=il(inh)
 
227
      ma=ma+1
 
228
      nib=nib-ia(ma)
 
229
      if(ma.ge.m)goto 12
 
230
      if(ia(ma).eq.0)goto 11
 
231
      ilinh=il(inh)
 
232
c on the left side of the tree, the right son of a right son must represent 
 
233
c a subset of points which is inside a triangle with vertices which are also 
 
234
c vertices of the convex polygon and hence the subset may be neglected.
 
235
      call husplit(n,x,ia(ma),ib(nib),ih(inh),ih(ilinh),2,ia,
 
236
     1  mb,mxa,ib(nib),mbb,mxb)
 
237
      ia(ma)=mbb
 
238
      goto 9
 
239
c now traverse the right half of the tree
 
240
12    mxb=mxbb
 
241
      ma=m
 
242
      mb=ia(ma)
 
243
      nia=1
 
244
      ia(ma)=0
 
245
c start with the right son
 
246
13    nia=nia+ia(ma)
 
247
      ma=ma-1
 
248
14    if(mxb.eq.0)goto 16
 
249
      il(nh)=il(inh)
 
250
      il(inh)=nh
 
251
      ih(nh)=ib(mxb)
 
252
      ib(mxb)=ib(mb)
 
253
      mb=mb-1
 
254
      nh=nh+1
 
255
      if(mb.eq.0)goto 15
 
256
      ilinh=il(inh)
 
257
      call husplit(n,x,mb,ib(nib),ih(inh),ih(ilinh),-1,ia(nia),
 
258
     1  ia(ma),mxa,ib(nib),mbb,mxb)
 
259
      mb=mbb
 
260
      goto 13
 
261
c then the left son
 
262
15    inh=il(inh)
 
263
16    inh=il(inh)
 
264
      ma=ma+1
 
265
      nia=nia-ia(ma)
 
266
      if(ma.eq.mp1)goto 17
 
267
      if(ia(ma).eq.0)goto 16
 
268
      ilinh=il(inh)
 
269
c on the right side of the tree, the left son of a left son must represent 
 
270
c a subset of points which is inside a triangle with vertices which are 
 
271
c also vertices of the convex polygon and hence the subset may be neglected.
 
272
      call husplit(n,x,ia(ma),ia(nia),ih(inh),ih(ilinh),-2,
 
273
     1  ia(nia),mbb,mxa,ib(nib),mb,mxb)
 
274
      goto 14
 
275
17    nh=nh-1
 
276
      return
 
277
c all the special cases are handled down here
 
278
c if all the points lie on a vertical line
 
279
18    kx=in(1)
 
280
      kn=in(1)
 
281
      do 20 i=1,m
 
282
        j=in(i)
 
283
        if(x(2,j).le.x(2,kx))goto 19
 
284
        mx=i
 
285
        kx=j
 
286
19      if(x(2,j).ge.x(2,kn))goto 20
 
287
        min=i
 
288
        kn=j
 
289
20    continue
 
290
      if(kx.eq.kn)goto 22
 
291
c if there are only two points
 
292
21    ih(1)=kx
 
293
      ih(2)=kn
 
294
      nh=3
 
295
      if((x(1,kn).eq.x(1,kx)).and.(x(2,kn).eq.x(2,kx)))nh=2
 
296
      goto 17
 
297
c if there is only one point
 
298
22    nh=2
 
299
      ih(1)=in(1)
 
300
      il(1)=1
 
301
      goto 17
 
302
c multiple extremes are handled here if there are several points 
 
303
c with the (same) largest first coordinate
 
304
23    if(.not.maxe)goto 25
 
305
      do 24 i=1,m
 
306
        j=in(i)
 
307
        if(x(1,j).ne.x(1,kx))goto 24
 
308
        if(x(2,j).le.x(2,kx))goto 24
 
309
        mx=i
 
310
        kx=j
 
311
24    continue
 
312
c if there are several points with the (same) smallest first coordinate
 
313
25    if(.not.mine)goto 7
 
314
      do 26 i=1,m
 
315
        j=in(i)
 
316
        if(x(1,j).ne.x(1,kn))goto 26
 
317
        if(x(2,j).ge.x(2,kn))goto 26
 
318
        min=i
 
319
        kn=j
 
320
26    continue
 
321
      goto 7
 
322
      end