[7541] | 1 | |
---|
| 2 | MODULE polygones |
---|
| 3 | ! |
---|
| 4 | ! This module is there to handle polygones. It is an essential tool in interpolation methods. |
---|
| 5 | ! |
---|
| 6 | USE defprec |
---|
| 7 | |
---|
| 8 | IMPLICIT NONE |
---|
| 9 | |
---|
| 10 | PRIVATE |
---|
| 11 | PUBLIC :: polygones_pointinside, polygones_extend, polygones_intersection, polygones_cleanup, & |
---|
| 12 | & polygones_area, polygones_crossing, polygones_convexhull |
---|
| 13 | ! |
---|
| 14 | REAL(r_std),PARAMETER :: zero=0.0 |
---|
| 15 | ! |
---|
| 16 | CONTAINS |
---|
| 17 | ! |
---|
| 18 | !======================================================================================================= |
---|
| 19 | ! |
---|
| 20 | SUBROUTINE polygones_pointinside(nvert_in, poly, point_x, point_y, inside) |
---|
| 21 | ! |
---|
| 22 | ! Routine to determine of point (point_x, point_y) is inside of polygone poly |
---|
| 23 | ! Mathode based on docume http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html |
---|
| 24 | ! |
---|
| 25 | ! Arguments |
---|
| 26 | INTEGER(i_std), INTENT(in) :: nvert_in |
---|
| 27 | REAL(r_std), DIMENSION(:,:), INTENT(in) :: poly |
---|
| 28 | REAL(r_std), INTENT(in) :: point_x, point_y |
---|
| 29 | LOGICAL, INTENT(out) :: inside |
---|
| 30 | ! |
---|
| 31 | ! Local |
---|
| 32 | INTEGER(i_std) :: i, j, nvert |
---|
| 33 | REAL(r_std) :: xonline |
---|
| 34 | ! |
---|
| 35 | ! |
---|
| 36 | inside = .FALSE. |
---|
| 37 | ! |
---|
| 38 | nvert=SIZE(poly,DIM=1) |
---|
| 39 | IF ( nvert < nvert_in) THEN |
---|
| 40 | WRITE(*,*) "Input polygone is smaller than the size proposed in the interface." |
---|
| 41 | STOP "polygones_pointinside" |
---|
| 42 | ENDIF |
---|
| 43 | IF (SIZE(poly,DIM=2) .NE. 2) THEN |
---|
| 44 | WRITE(*,*) "This cannot be a polygone :", SIZE(poly) |
---|
| 45 | STOP "polygones_pointinside" |
---|
| 46 | ENDIF |
---|
| 47 | ! |
---|
| 48 | j = nvert_in |
---|
| 49 | DO i=1,nvert_in |
---|
| 50 | IF ( (poly(i,2) > point_y) .NEQV. (poly(j,2) > point_y) ) THEN |
---|
| 51 | xonline = (poly(j,1)-poly(i,1))*(point_y-poly(i,2))/(poly(j,2)-poly(i,2))+poly(i,1) |
---|
| 52 | IF ( point_x < xonline ) THEN |
---|
| 53 | inside = (.NOT. inside) |
---|
| 54 | ENDIF |
---|
| 55 | ENDIF |
---|
| 56 | j = i |
---|
| 57 | ENDDO |
---|
| 58 | ! |
---|
| 59 | END SUBROUTINE polygones_pointinside |
---|
| 60 | ! |
---|
| 61 | !======================================================================================================= |
---|
| 62 | ! |
---|
| 63 | SUBROUTINE polygones_lineintersect(la, lb, intersection, point_x, point_y) |
---|
| 64 | ! |
---|
| 65 | ! Simple alogorith tho determine the intersections of 2 lines. |
---|
| 66 | ! |
---|
| 67 | ! ARGUMENTS |
---|
| 68 | REAL(r_std), DIMENSION(2,2), INTENT(in) :: la, lb |
---|
| 69 | LOGICAL, INTENT(out) :: intersection |
---|
| 70 | REAL(r_std), INTENT(out) :: point_x, point_y |
---|
| 71 | ! |
---|
| 72 | ! Local |
---|
| 73 | REAL(r_std) :: den, ua, ub |
---|
| 74 | REAL(r_std) :: x1, x2, x3, x4, y1, y2, y3, y4 |
---|
| 75 | ! |
---|
| 76 | intersection = .FALSE. |
---|
| 77 | ! |
---|
| 78 | x1 = la(1,1) |
---|
| 79 | y1 = la(1,2) |
---|
| 80 | x2 = la(2,1) |
---|
| 81 | y2 = la(2,2) |
---|
| 82 | |
---|
| 83 | x3 = lb(1,1) |
---|
| 84 | y3 = lb(1,2) |
---|
| 85 | x4 = lb(2,1) |
---|
| 86 | y4 = lb(2,2) |
---|
| 87 | |
---|
| 88 | den = (y4-y3)*(x2-x1)-(x4-x3)*(y2-y1) |
---|
| 89 | |
---|
| 90 | IF ( ABS(den) > EPSILON(den) ) THEN |
---|
| 91 | ua = ((x4-x3)*(y1-y3)-(y4-y3)*(x1-x3))/den |
---|
| 92 | ub = ((x2-x1)*(y1-y3)-(y2-y1)*(x1-x3))/den |
---|
| 93 | IF ( ua >= 0 .AND. ua <= 1 .AND. ub >= 0 .AND. ub <= 1) THEN |
---|
| 94 | intersection = .TRUE. |
---|
| 95 | point_x = x1 + ua*(x2-x1) |
---|
| 96 | point_y = y1 + ua*(y2-y1) |
---|
| 97 | ENDIF |
---|
| 98 | ENDIF |
---|
| 99 | |
---|
| 100 | END SUBROUTINE polygones_lineintersect |
---|
| 101 | ! |
---|
| 102 | !======================================================================================================= |
---|
| 103 | ! |
---|
| 104 | SUBROUTINE polygones_extend(nvert_in, poly_in, nbdots, nvert_out, poly_out) |
---|
| 105 | ! |
---|
| 106 | ! This function simply add nbdots onto each side of the polygone. |
---|
| 107 | ! |
---|
| 108 | ! Arguments |
---|
| 109 | INTEGER(i_std), INTENT(in) :: nvert_in |
---|
| 110 | REAL(r_std), DIMENSION(:,:), INTENT(in) :: poly_in |
---|
| 111 | INTEGER(i_std), INTENT(in) :: nbdots |
---|
| 112 | INTEGER(i_std), INTENT(out) :: nvert_out |
---|
| 113 | REAL(r_std), DIMENSION(:,:), INTENT(out) :: poly_out |
---|
| 114 | ! |
---|
| 115 | ! Local |
---|
| 116 | ! |
---|
| 117 | INTEGER(i_std) :: nvert, nvert_tmp, i, j, id, ipos |
---|
| 118 | REAL(r_std) :: xs, xe, ys, ye |
---|
| 119 | ! |
---|
| 120 | nvert=SIZE(poly_in,DIM=1) |
---|
| 121 | IF ( nvert < nvert_in ) THEN |
---|
| 122 | WRITE(*,*) "Input polygone is smaller than the size proposed in the interface." |
---|
| 123 | STOP "polygones_extend" |
---|
| 124 | ENDIF |
---|
| 125 | IF (SIZE(poly_in,DIM=2) .NE. 2) THEN |
---|
| 126 | WRITE(*,*) "Poly_in cannot be a polygone :", SIZE(poly_in) |
---|
| 127 | STOP "polygones_extend" |
---|
| 128 | ENDIF |
---|
| 129 | nvert_tmp=SIZE(poly_out,DIM=1) |
---|
| 130 | IF (SIZE(poly_out,DIM=2) .NE. 2) THEN |
---|
| 131 | WRITE(*,*) "Poly_out cannot be a polygone :", SIZE(poly_out) |
---|
| 132 | STOP "polygones_extend" |
---|
| 133 | ENDIF |
---|
| 134 | ! |
---|
| 135 | IF ( nvert_tmp < nvert_in*nbdots) THEN |
---|
| 136 | WRITE(*,*) "Poly_out is too small for extending it :", SIZE(poly_out), " we would need :", nvert_in*nbdots |
---|
| 137 | STOP "polygones_extend" |
---|
| 138 | ENDIF |
---|
| 139 | nvert_out = nvert_in*nbdots |
---|
| 140 | ! |
---|
| 141 | j = nvert_in |
---|
| 142 | ipos = 0 |
---|
| 143 | DO i=1,nvert_in |
---|
| 144 | xs = poly_in(j,1) |
---|
| 145 | xe = poly_in(i,1) |
---|
| 146 | ys = poly_in(j,2) |
---|
| 147 | ye = poly_in(i,2) |
---|
| 148 | IF (xs == xe ) THEN |
---|
| 149 | DO id=1,nbdots |
---|
| 150 | ipos = ipos + 1 |
---|
| 151 | poly_out(ipos,1) = xs |
---|
| 152 | poly_out(ipos,2) = ys + (id-1)*(ye-ys)/nbdots |
---|
| 153 | ENDDO |
---|
| 154 | ELSE IF (ys == ye ) THEN |
---|
| 155 | DO id=1,nbdots |
---|
| 156 | ipos = ipos + 1 |
---|
| 157 | poly_out(ipos,1) = xs + (id-1)*(xe-xs)/nbdots |
---|
| 158 | poly_out(ipos,2) = ys |
---|
| 159 | ENDDO |
---|
| 160 | ELSE |
---|
| 161 | DO id=1,nbdots |
---|
| 162 | ipos = ipos + 1 |
---|
| 163 | poly_out(ipos,2) = ys + (id-1)*(ye-ys)/nbdots |
---|
| 164 | poly_out(ipos,1) = (xs-xe)*(poly_out(ipos,2)-ye)/(ys-ye)+xe |
---|
| 165 | ENDDO |
---|
| 166 | ENDIF |
---|
| 167 | j = i |
---|
| 168 | ENDDO |
---|
| 169 | ! |
---|
| 170 | END SUBROUTINE polygones_extend |
---|
| 171 | ! |
---|
| 172 | !======================================================================================================= |
---|
| 173 | ! |
---|
| 174 | SUBROUTINE polygones_intersection(nvert_a, poly_a, nvert_b, poly_b, nbpts, poly_out) |
---|
| 175 | ! |
---|
| 176 | ! Find the polygon resulting from the intersection of poly_a and poly_b |
---|
| 177 | ! |
---|
| 178 | ! Arguments |
---|
| 179 | INTEGER(i_std), INTENT(in) :: nvert_a |
---|
| 180 | REAL(r_std), DIMENSION(:,:), INTENT(in) :: poly_a |
---|
| 181 | INTEGER(i_std), INTENT(in) :: nvert_b |
---|
| 182 | REAL(r_std), DIMENSION(:,:), INTENT(in) :: poly_b |
---|
| 183 | INTEGER(i_std), INTENT(out) :: nbpts |
---|
| 184 | REAL(r_std), DIMENSION(:,:), INTENT(out) :: poly_out |
---|
| 185 | ! |
---|
| 186 | ! Local |
---|
| 187 | INTEGER(i_std) :: nvert_tmpa, nvert_tmpb, nvert_out, i, j |
---|
| 188 | LOGICAL :: inside |
---|
| 189 | INTEGER(i_std) :: in_a, in_b |
---|
| 190 | INTEGER(i_std), ALLOCATABLE, DIMENSION(:) :: ptin_a, ptin_b |
---|
| 191 | ! |
---|
| 192 | nvert_tmpa=SIZE(poly_a,DIM=1) |
---|
| 193 | IF ( nvert_tmpa < nvert_a) THEN |
---|
| 194 | WRITE(*,*) "Input polygone (poly_a) is smaller than the size proposed in the interface." |
---|
| 195 | STOP "polygones_intersection" |
---|
| 196 | ENDIF |
---|
| 197 | IF (SIZE(poly_a,DIM=2) .NE. 2) THEN |
---|
| 198 | WRITE(*,*) "Poly_a cannot be a polygone :", SIZE(poly_a) |
---|
| 199 | STOP "polygones_intersection" |
---|
| 200 | ENDIF |
---|
| 201 | ALLOCATE(ptin_a(nvert_a)) |
---|
| 202 | nvert_tmpb=SIZE(poly_b,DIM=1) |
---|
| 203 | IF ( nvert_tmpb < nvert_a) THEN |
---|
| 204 | WRITE(*,*) "Input polygone (poly_b) is smaller than the size proposed in the interface." |
---|
| 205 | STOP "polygones_intersection" |
---|
| 206 | ENDIF |
---|
| 207 | IF (SIZE(poly_b,DIM=2) .NE. 2) THEN |
---|
| 208 | WRITE(*,*) "Poly_b cannot be a polygone :", SIZE(poly_b) |
---|
| 209 | STOP "polygones_intersection" |
---|
| 210 | ENDIF |
---|
| 211 | ALLOCATE(ptin_b(nvert_b)) |
---|
| 212 | ! |
---|
| 213 | ! |
---|
| 214 | in_a = 0 |
---|
| 215 | DO i=1,nvert_a |
---|
| 216 | CALL polygones_pointinside(nvert_b, poly_b, poly_a(i,1), poly_a(i,2), inside) |
---|
| 217 | IF ( inside ) THEN |
---|
| 218 | in_a = in_a + 1 |
---|
| 219 | ptin_a(in_a) = i |
---|
| 220 | ENDIF |
---|
| 221 | ENDDO |
---|
| 222 | ! |
---|
| 223 | in_b = 0 |
---|
| 224 | DO i=1,nvert_b |
---|
| 225 | CALL polygones_pointinside(nvert_a, poly_a, poly_b(i,1), poly_b(i,2), inside) |
---|
| 226 | IF ( inside ) THEN |
---|
| 227 | in_b = in_b + 1 |
---|
| 228 | ptin_b(in_b) = i |
---|
| 229 | ENDIF |
---|
| 230 | ENDDO |
---|
| 231 | ! |
---|
| 232 | nbpts = in_a + in_b |
---|
| 233 | ! |
---|
| 234 | nvert_out=SIZE(poly_out,DIM=1) |
---|
| 235 | IF (SIZE(poly_out,DIM=2) .NE. 2) THEN |
---|
| 236 | WRITE(*,*) "Poly_out cannot be a polygone :", SIZE(poly_out) |
---|
| 237 | STOP "polygones_intersection" |
---|
| 238 | ENDIF |
---|
| 239 | IF ( nbpts <= nvert_out ) THEN |
---|
| 240 | i = 0 |
---|
| 241 | DO j=1,in_a |
---|
| 242 | i = i+1 |
---|
| 243 | poly_out(i,:) = poly_a(ptin_a(j),:) |
---|
| 244 | ENDDO |
---|
| 245 | DO j=1,in_b |
---|
| 246 | i = i+1 |
---|
| 247 | poly_out(i,:) = poly_b(ptin_b(j),:) |
---|
| 248 | ENDDO |
---|
| 249 | ELSE |
---|
| 250 | WRITE(*,*) "The intersection polygone has", nbpts, "points and that is more than available in poly_out", SIZE(poly_out) |
---|
| 251 | STOP "polygones_intersection" |
---|
| 252 | ENDIF |
---|
| 253 | ! |
---|
| 254 | DEALLOCATE(ptin_a, ptin_b) |
---|
| 255 | ! |
---|
| 256 | END SUBROUTINE polygones_intersection |
---|
| 257 | ! |
---|
| 258 | !======================================================================================================= |
---|
| 259 | ! |
---|
| 260 | SUBROUTINE polygones_cleanup(nvert_in, poly_in, nvert_out, poly_out) |
---|
| 261 | ! |
---|
| 262 | ! Clean_up the polygone by deleting points which are redundant and ordering based |
---|
| 263 | ! on the proximity of the points. Beware this does not give a convex hull to the surface |
---|
| 264 | ! inside the polygone. |
---|
| 265 | ! |
---|
| 266 | ! Arguments |
---|
| 267 | INTEGER(i_std), INTENT(in) :: nvert_in |
---|
| 268 | REAL(r_std), DIMENSION(:,:), INTENT(in) :: poly_in |
---|
| 269 | INTEGER(i_std), INTENT(out) :: nvert_out |
---|
| 270 | REAL(r_std), DIMENSION(:,:), INTENT(out) :: poly_out |
---|
| 271 | ! |
---|
| 272 | ! Local |
---|
| 273 | ! |
---|
| 274 | INTEGER(i_std) :: nvert, nvert_tmpin, nvert_tmp |
---|
| 275 | INTEGER(i_std) :: i, j |
---|
| 276 | INTEGER(i_std), DIMENSION(1) :: ismin |
---|
| 277 | REAL(r_std), ALLOCATABLE, DIMENSION(:) :: dist |
---|
| 278 | REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: poly_tmp |
---|
| 279 | ! |
---|
| 280 | nvert_tmpin=SIZE(poly_in,DIM=1) |
---|
| 281 | IF ( nvert_tmpin < nvert_in ) THEN |
---|
| 282 | WRITE(*,*) "Input polygone is smaller than the size proposed in the interface." |
---|
| 283 | STOP "polygones_cleanup" |
---|
| 284 | ENDIF |
---|
| 285 | IF (SIZE(poly_in,DIM=2) .NE. 2) THEN |
---|
| 286 | WRITE(*,*) "Poly_in cannot be a polygone :", SIZE(poly_in) |
---|
| 287 | STOP "polygones_cleanup" |
---|
| 288 | ENDIF |
---|
| 289 | nvert_tmp=SIZE(poly_out,DIM=1) |
---|
| 290 | IF (SIZE(poly_out,DIM=2) .NE. 2) THEN |
---|
| 291 | WRITE(*,*) "Poly_out cannot be a polygone :", SIZE(poly_out) |
---|
| 292 | STOP "polygones_cleanup" |
---|
| 293 | ENDIF |
---|
| 294 | ! |
---|
| 295 | ALLOCATE(poly_tmp(nvert_in,2)) |
---|
| 296 | ALLOCATE(dist(nvert_in)) |
---|
| 297 | ! |
---|
| 298 | nvert = nvert_in-1 |
---|
| 299 | poly_tmp(1:nvert,:) = poly_in(2:nvert_in,:) |
---|
| 300 | poly_out(1,:) = poly_in(1,:) |
---|
| 301 | nvert_out = 1 |
---|
| 302 | ! |
---|
| 303 | DO WHILE ( nvert > 0 ) |
---|
| 304 | ! |
---|
| 305 | DO j=1,nvert |
---|
| 306 | dist(j) = SQRT((poly_out(nvert_out,1)-poly_tmp(j,1))**2+(poly_out(nvert_out,2)-poly_tmp(j,2))**2) |
---|
| 307 | ENDDO |
---|
| 308 | ! |
---|
| 309 | ismin = MINLOC(dist(1:nvert)) |
---|
| 310 | ! |
---|
| 311 | IF ( dist(ismin(1)) > EPSILON(dist(1)) ) THEN |
---|
| 312 | ! |
---|
| 313 | ! The smallest distance between 2 vertices is larger than zero : |
---|
| 314 | ! Add the vertex to poly_out and delete it from poly_tmp |
---|
| 315 | ! |
---|
| 316 | nvert_out = nvert_out + 1 |
---|
| 317 | IF ( nvert_out > nvert_tmp) THEN |
---|
| 318 | WRITE(*,*) "Output polygone too small" |
---|
| 319 | STOP "polygones_cleanup" |
---|
| 320 | ENDIF |
---|
| 321 | poly_out(nvert_out,:) = poly_tmp(ismin(1),:) |
---|
| 322 | poly_tmp(ismin(1):nvert-1,:) = poly_tmp(ismin(1)+1:nvert,:) |
---|
| 323 | ELSE |
---|
| 324 | ! |
---|
| 325 | ! Else the vertex already exists in poly_out and thus we just need |
---|
| 326 | ! to delete it from poly_tmp |
---|
| 327 | ! |
---|
| 328 | poly_tmp(ismin(1):nvert-1,:) = poly_tmp(ismin(1)+1:nvert,:) |
---|
| 329 | ENDIF |
---|
| 330 | ! |
---|
| 331 | ! One less vertices exists in poly_tmp |
---|
| 332 | ! |
---|
| 333 | nvert = nvert-1 |
---|
| 334 | ! |
---|
| 335 | ENDDO |
---|
| 336 | |
---|
| 337 | DEALLOCATE(poly_tmp, dist) |
---|
| 338 | |
---|
| 339 | END SUBROUTINE polygones_cleanup |
---|
| 340 | ! |
---|
| 341 | !======================================================================================================= |
---|
| 342 | ! |
---|
| 343 | SUBROUTINE polygones_area(nvert_in, poly_in, dx, dy, area) |
---|
| 344 | ! |
---|
| 345 | ! Find the polygon resulting from the intersection of poly_a and poly_b |
---|
| 346 | ! |
---|
| 347 | ! Arguments |
---|
| 348 | INTEGER(i_std), INTENT(in) :: nvert_in |
---|
| 349 | REAL(r_std), DIMENSION(:,:), INTENT(in) :: poly_in |
---|
| 350 | REAL(r_std), INTENT(in) :: dx, dy |
---|
| 351 | REAL(r_std), INTENT(out) :: area |
---|
| 352 | ! |
---|
| 353 | ! Local |
---|
| 354 | ! |
---|
| 355 | INTEGER(i_std) :: nvert, i, j |
---|
| 356 | ! |
---|
| 357 | nvert = SIZE(poly_in,DIM=1) |
---|
| 358 | IF ( nvert < nvert_in) THEN |
---|
| 359 | WRITE(*,*) "Input polygone is smaller than the size proposed in the interface." |
---|
| 360 | STOP "polygones_area" |
---|
| 361 | ENDIF |
---|
| 362 | IF (SIZE(poly_in,DIM=2) .NE. 2) THEN |
---|
| 363 | WRITE(*,*) "Poly_in cannot be a polygone :", SIZE(poly_in) |
---|
| 364 | STOP "polygones_area" |
---|
| 365 | ENDIF |
---|
| 366 | ! |
---|
| 367 | area = zero |
---|
| 368 | j = nvert_in |
---|
| 369 | DO i=1,nvert_in |
---|
| 370 | !! area = area + (dy*(poly_in(j,2)+poly_in(i,2))/2.0)*(dx*(poly_in(j,1)-poly_in(i,1))) |
---|
| 371 | area = area + dy*dx/2.0*(poly_in(j,2)+poly_in(i,2))*(poly_in(j,1)-poly_in(i,1)) |
---|
| 372 | j = i |
---|
| 373 | ENDDO |
---|
| 374 | ! |
---|
| 375 | area = ABS(area) |
---|
| 376 | ! |
---|
| 377 | END SUBROUTINE polygones_area |
---|
| 378 | ! |
---|
| 379 | !======================================================================================================= |
---|
| 380 | ! |
---|
| 381 | SUBROUTINE polygones_crossing(nvert_a, poly_a, nvert_b, poly_b, nbpts, crossings) |
---|
| 382 | ! |
---|
| 383 | ! Find the polygon resulting from the intersection of poly_a and poly_b |
---|
| 384 | ! |
---|
| 385 | ! Arguments |
---|
| 386 | INTEGER(i_std), INTENT(in) :: nvert_a |
---|
| 387 | REAL(r_std), DIMENSION(:,:), INTENT(in) :: poly_a |
---|
| 388 | INTEGER(i_std), INTENT(in) :: nvert_b |
---|
| 389 | REAL(r_std), DIMENSION(:,:), INTENT(in) :: poly_b |
---|
| 390 | INTEGER(i_std), INTENT(out) :: nbpts |
---|
| 391 | REAL(r_std), DIMENSION(:,:), INTENT(out) :: crossings |
---|
| 392 | ! |
---|
| 393 | ! Local |
---|
| 394 | ! |
---|
| 395 | INTEGER(i_std) :: nvert, ia, ja, ib, jb |
---|
| 396 | REAL(r_std), DIMENSION(2,2) :: la, lb |
---|
| 397 | REAL(r_std) :: x, y |
---|
| 398 | LOGICAL :: intersect |
---|
| 399 | ! |
---|
| 400 | nvert=SIZE(poly_a,DIM=1) |
---|
| 401 | IF ( nvert < nvert_a) THEN |
---|
| 402 | WRITE(*,*) "Input polygone (poly_a) is smaller than the size proposed in the interface." |
---|
| 403 | STOP "polygones_crossing" |
---|
| 404 | ENDIF |
---|
| 405 | IF (SIZE(poly_a,DIM=2) .NE. 2) THEN |
---|
| 406 | WRITE(*,*) "Poly_a cannot be a polygone :", SIZE(poly_a) |
---|
| 407 | STOP "polygones_crossing" |
---|
| 408 | ENDIF |
---|
| 409 | nvert=SIZE(poly_b,DIM=1) |
---|
| 410 | IF ( nvert < nvert_b) THEN |
---|
| 411 | WRITE(*,*) "Input polygone (poly_b) is smaller than the size proposed in the interface." |
---|
| 412 | STOP "polygones_crossing" |
---|
| 413 | ENDIF |
---|
| 414 | IF (SIZE(poly_b,DIM=2) .NE. 2) THEN |
---|
| 415 | WRITE(*,*) "Poly_b cannot be a polygone :", SIZE(poly_b) |
---|
| 416 | STOP "polygones_crossing" |
---|
| 417 | ENDIF |
---|
| 418 | nvert=SIZE(crossings,DIM=1) |
---|
| 419 | IF (SIZE(crossings,DIM=2) .NE. 2) THEN |
---|
| 420 | WRITE(*,*) "Poly_a cannot be a polygone :", SIZE(poly_a) |
---|
| 421 | STOP "polygones_crossing" |
---|
| 422 | ENDIF |
---|
| 423 | ! |
---|
| 424 | nbpts = 0 |
---|
| 425 | ! |
---|
| 426 | ja = nvert_a |
---|
| 427 | DO ia=1,nvert_a |
---|
| 428 | ! |
---|
| 429 | ! Put one side of the polygone a into la |
---|
| 430 | ! |
---|
| 431 | la(1,1) = poly_a(ja,1) |
---|
| 432 | la(1,2) = poly_a(ja,2) |
---|
| 433 | la(2,1) = poly_a(ia,1) |
---|
| 434 | la(2,2) = poly_a(ia,2) |
---|
| 435 | ! |
---|
| 436 | jb = nvert_b |
---|
| 437 | DO ib=1,nvert_b |
---|
| 438 | ! |
---|
| 439 | ! Put one side of the polygone b into lb |
---|
| 440 | ! |
---|
| 441 | lb(1,1) = poly_b(jb,1) |
---|
| 442 | lb(1,2) = poly_b(jb,2) |
---|
| 443 | lb(2,1) = poly_b(ib,1) |
---|
| 444 | lb(2,2) = poly_b(ib,2) |
---|
| 445 | ! |
---|
| 446 | CALL polygones_lineintersect(la, lb, intersect, x, y) |
---|
| 447 | ! |
---|
| 448 | IF ( intersect ) THEN |
---|
| 449 | nbpts = nbpts + 1 |
---|
| 450 | IF ( nbpts <= nvert ) THEN |
---|
| 451 | crossings(nbpts, 1) = x |
---|
| 452 | crossings(nbpts, 2) = y |
---|
| 453 | ELSE |
---|
| 454 | WRITE(*,*) "Polygone to write the intersection points is too small", nvert, nbpts |
---|
| 455 | STOP "polygones_crossing" |
---|
| 456 | ENDIF |
---|
| 457 | ENDIF |
---|
| 458 | ! |
---|
| 459 | jb = ib |
---|
| 460 | ENDDO |
---|
| 461 | ja = ia |
---|
| 462 | ENDDO |
---|
| 463 | ! |
---|
| 464 | |
---|
| 465 | END SUBROUTINE polygones_crossing |
---|
| 466 | ! |
---|
| 467 | !======================================================================================================= |
---|
| 468 | ! |
---|
| 469 | SUBROUTINE polygones_convexhull(nvert_in, poly_in, nvert_out, poly_out) |
---|
| 470 | ! |
---|
| 471 | ! Routine orders points in poly_in so that they form a convex shape. This is based on the Graham scan |
---|
| 472 | ! algorith as implemented by Alan Miller : http://jblevins.org/mirror/amiller/ |
---|
| 473 | ! |
---|
| 474 | ! The output polygone can be smaller than the input one as we are computing the convex envelope. |
---|
| 475 | ! |
---|
| 476 | ! Arguments |
---|
| 477 | INTEGER(i_std), INTENT(in) :: nvert_in |
---|
| 478 | REAL(r_std), DIMENSION(:,:), INTENT(in) :: poly_in |
---|
| 479 | INTEGER(i_std), INTENT(out) :: nvert_out |
---|
| 480 | REAL(r_std), DIMENSION(:,:), INTENT(out) :: poly_out |
---|
| 481 | ! |
---|
| 482 | ! Local |
---|
| 483 | ! |
---|
| 484 | INTEGER(i_std) :: nvert, nvert_tmp |
---|
| 485 | INTEGER(i_std), ALLOCATABLE, DIMENSION(:):: iwk, next, vertex |
---|
| 486 | REAL(r_std) :: xmin, temp, dy, dx2, dy1, dy2 |
---|
| 487 | REAL(r_std) :: x1, x2, y1, y2 |
---|
| 488 | REAL(r_std) :: xmax, ymin, ymax |
---|
| 489 | REAL(r_std) :: dx, dx1 |
---|
| 490 | REAL(r_std) :: dmax, dmax1, dmax2 |
---|
| 491 | REAL(r_std) :: dist, dmin |
---|
| 492 | INTEGER(i_std) :: i, i1, i2, j, jp1, jp2, i2save, i3, i2next |
---|
| 493 | LOGICAL :: points_todo |
---|
| 494 | INTEGER(i_std), PARAMETER :: nextinc=20 |
---|
| 495 | ! |
---|
| 496 | nvert_tmp=SIZE(poly_in,DIM=1) |
---|
| 497 | IF ( nvert_tmp < nvert_in) THEN |
---|
| 498 | WRITE(*,*) "Input polygone is smaller than the size proposed in the interface." |
---|
| 499 | STOP "polygones_convexhull" |
---|
| 500 | ENDIF |
---|
| 501 | IF ( nvert_tmp <= 2 ) THEN |
---|
| 502 | WRITE(*,*) "The input polygone is too small to compute a convex hull." |
---|
| 503 | STOP "polygones_convexhull" |
---|
| 504 | ENDIF |
---|
| 505 | IF (SIZE(poly_in,DIM=2) .NE. 2) THEN |
---|
| 506 | WRITE(*,*) "Poly_a cannot be a polygone :", SIZE(poly_in) |
---|
| 507 | STOP "polygones_convexhull" |
---|
| 508 | ENDIF |
---|
| 509 | IF (SIZE(poly_out,DIM=2) .NE. 2) THEN |
---|
| 510 | WRITE(*,*) "Poly_b cannot be a polygone :", SIZE(poly_out) |
---|
| 511 | STOP "polygones_convexhull" |
---|
| 512 | ENDIF |
---|
| 513 | ! |
---|
| 514 | ALLOCATE(vertex(nvert_in)) |
---|
| 515 | ALLOCATE(iwk(nvert_in)) |
---|
| 516 | ALLOCATE(next(nvert_in*nextinc)) |
---|
| 517 | ! |
---|
| 518 | ! |
---|
| 519 | ! Choose the points with smallest & largest x- values as the |
---|
| 520 | ! first two vertices of the polygon. |
---|
| 521 | ! |
---|
| 522 | IF (poly_in(1,1) > poly_in(nvert_in,1)) THEN |
---|
| 523 | vertex(1) = nvert_in |
---|
| 524 | vertex(2) = 1 |
---|
| 525 | xmin = poly_in(nvert_in,1) |
---|
| 526 | xmax = poly_in(1,1) |
---|
| 527 | ELSE |
---|
| 528 | vertex(1) = 1 |
---|
| 529 | vertex(2) = nvert_in |
---|
| 530 | xmin = poly_in(1,1) |
---|
| 531 | xmax = poly_in(nvert_in,1) |
---|
| 532 | END IF |
---|
| 533 | |
---|
| 534 | DO i = 2, nvert_in-1 |
---|
| 535 | temp = poly_in(i,1) |
---|
| 536 | IF (temp < xmin) THEN |
---|
| 537 | vertex(1) = i |
---|
| 538 | xmin = temp |
---|
| 539 | ELSE IF (temp > xmax) THEN |
---|
| 540 | vertex(2) = i |
---|
| 541 | xmax = temp |
---|
| 542 | END IF |
---|
| 543 | END DO |
---|
| 544 | ! |
---|
| 545 | ! Special case, xmax = xmin. |
---|
| 546 | ! |
---|
| 547 | IF (xmax == xmin) THEN |
---|
| 548 | IF (poly_in(1,2) > poly_in(nvert_in,2)) THEN |
---|
| 549 | vertex(1) = nvert_in |
---|
| 550 | vertex(2) = 1 |
---|
| 551 | ymin = poly_in(nvert_in,2) |
---|
| 552 | ymax = poly_in(1,2) |
---|
| 553 | ELSE |
---|
| 554 | vertex(1) = 1 |
---|
| 555 | vertex(2) = nvert_in |
---|
| 556 | ymin = poly_in(1,2) |
---|
| 557 | ymax = poly_in(nvert_in,2) |
---|
| 558 | END IF |
---|
| 559 | |
---|
| 560 | DO i = 2, nvert_in-1 |
---|
| 561 | temp = poly_in(i,2) |
---|
| 562 | IF (temp < ymin) THEN |
---|
| 563 | vertex(1) = i |
---|
| 564 | ymin = temp |
---|
| 565 | ELSE IF (temp > ymax) THEN |
---|
| 566 | vertex(2) = i |
---|
| 567 | ymax = temp |
---|
| 568 | END IF |
---|
| 569 | END DO |
---|
| 570 | |
---|
| 571 | nvert = 2 |
---|
| 572 | IF (ymax == ymin) nvert = 1 |
---|
| 573 | RETURN |
---|
| 574 | END IF |
---|
| 575 | ! |
---|
| 576 | ! Set up two initial lists of points; those points above & those below the |
---|
| 577 | ! line joining the first two vertices. next(i) will hold the pointer to the |
---|
| 578 | ! point furthest from the line joining vertex(i) to vertex(i+1) on the left |
---|
| 579 | ! hand side. |
---|
| 580 | ! |
---|
| 581 | i1 = vertex(1) |
---|
| 582 | i2 = vertex(2) |
---|
| 583 | iwk(i1) = -1 |
---|
| 584 | iwk(i2) = -1 |
---|
| 585 | dx = xmax - xmin |
---|
| 586 | y1 = poly_in(i1,2) |
---|
| 587 | dy = poly_in(i2,2) - y1 |
---|
| 588 | dmax = zero |
---|
| 589 | dmin = zero |
---|
| 590 | next(1) = -1 |
---|
| 591 | next(2) = -1 |
---|
| 592 | ! |
---|
| 593 | DO i = 1, nvert_in |
---|
| 594 | IF (i == vertex(1) .OR. i == vertex(2)) CYCLE |
---|
| 595 | dist = (poly_in(i,2) - y1)*dx - (poly_in(i,1) - xmin)*dy |
---|
| 596 | IF (dist > zero) THEN |
---|
| 597 | iwk(i1) = i |
---|
| 598 | i1 = i |
---|
| 599 | IF (dist > dmax) THEN |
---|
| 600 | next(1) = i |
---|
| 601 | dmax = dist |
---|
| 602 | END IF |
---|
| 603 | ELSE IF (dist < zero) THEN |
---|
| 604 | iwk(i2) = i |
---|
| 605 | i2 = i |
---|
| 606 | IF (dist < dmin) THEN |
---|
| 607 | next(2) = i |
---|
| 608 | dmin = dist |
---|
| 609 | END IF |
---|
| 610 | END IF |
---|
| 611 | END DO |
---|
| 612 | ! |
---|
| 613 | ! Ends of lists are indicated by pointers to -ve positions. |
---|
| 614 | ! |
---|
| 615 | iwk(i1) = -1 |
---|
| 616 | iwk(i2) = -1 |
---|
| 617 | nvert = 2 |
---|
| 618 | ! |
---|
| 619 | j = 1 |
---|
| 620 | ! |
---|
| 621 | ! Start of main process. |
---|
| 622 | ! |
---|
| 623 | ! Introduce new vertex between vertices j & j+1, if one has been found. |
---|
| 624 | ! Otherwise increase j. Exit if no more vertices. |
---|
| 625 | ! |
---|
| 626 | ! |
---|
| 627 | points_todo = .TRUE. |
---|
| 628 | DO WHILE ( points_todo ) |
---|
| 629 | ! |
---|
| 630 | DO WHILE ( points_todo .AND. next(j) < 0 ) |
---|
| 631 | IF (j == nvert) points_todo = .FALSE. |
---|
| 632 | j = j + 1 |
---|
| 633 | ENDDO |
---|
| 634 | ! |
---|
| 635 | IF ( points_todo ) THEN |
---|
| 636 | ! |
---|
| 637 | jp1 = j + 1 |
---|
| 638 | IF ( jp1 >= nvert_in*nextinc) THEN |
---|
| 639 | STOP "polygones_convexhull : please increase nextinc" |
---|
| 640 | ENDIF |
---|
| 641 | ! |
---|
| 642 | DO i = nvert, jp1, -1 |
---|
| 643 | vertex(i+1) = vertex(i) |
---|
| 644 | next(i+1) = next(i) |
---|
| 645 | END DO |
---|
| 646 | jp2 = jp1 + 1 |
---|
| 647 | nvert = nvert + 1 |
---|
| 648 | IF (jp2 > nvert) jp2 = 1 |
---|
| 649 | i1 = vertex(j) |
---|
| 650 | i2 = next(j) |
---|
| 651 | i3 = vertex(jp2) |
---|
| 652 | vertex(jp1) = i2 |
---|
| 653 | ! |
---|
| 654 | ! Process the list of points associated with vertex j. New list at vertex j |
---|
| 655 | ! consists of those points to the left of the line joining it to the new |
---|
| 656 | ! vertex (j+1). Similarly for the list at the new vertex. |
---|
| 657 | ! Points on or to the right of these lines are dropped. |
---|
| 658 | ! |
---|
| 659 | x1 = poly_in(i1,1) |
---|
| 660 | x2 = poly_in(i2,1) |
---|
| 661 | y1 = poly_in(i1,2) |
---|
| 662 | y2 = poly_in(i2,2) |
---|
| 663 | ! |
---|
| 664 | dx1 = x2 - x1 |
---|
| 665 | dx2 = poly_in(i3,1) - x2 |
---|
| 666 | dy1 = y2 - y1 |
---|
| 667 | dy2 = poly_in(i3,2) - y2 |
---|
| 668 | dmax1 = zero |
---|
| 669 | dmax2 = zero |
---|
| 670 | next(j) = -1 |
---|
| 671 | next(jp1) = -1 |
---|
| 672 | i2save = i2 |
---|
| 673 | i2next = iwk(i2) |
---|
| 674 | i = iwk(i1) |
---|
| 675 | iwk(i1) = -1 |
---|
| 676 | iwk(i2) = -1 |
---|
| 677 | ! |
---|
| 678 | DO WHILE ( i > 0 ) |
---|
| 679 | IF (i /= i2save) THEN |
---|
| 680 | dist = (poly_in(i,2) - y1)*dx1 - (poly_in(i,1) - x1)*dy1 |
---|
| 681 | IF (dist > zero) THEN |
---|
| 682 | iwk(i1) = i |
---|
| 683 | i1 = i |
---|
| 684 | IF (dist > DMAX1) THEN |
---|
| 685 | next(j) = i |
---|
| 686 | dmax1 = dist |
---|
| 687 | END IF |
---|
| 688 | ELSE |
---|
| 689 | dist = (poly_in(i,2) - y2)*dx2 - (poly_in(i,1) - x2)*dy2 |
---|
| 690 | IF (dist > zero) THEN |
---|
| 691 | iwk(i2) = i |
---|
| 692 | i2 = i |
---|
| 693 | IF (dist > dmax2) THEN |
---|
| 694 | next(jp1) = i |
---|
| 695 | dmax2 = dist |
---|
| 696 | END IF |
---|
| 697 | END IF |
---|
| 698 | END IF |
---|
| 699 | i = iwk(i) |
---|
| 700 | ELSE |
---|
| 701 | i = i2next |
---|
| 702 | END IF |
---|
| 703 | ! |
---|
| 704 | ! Get next point from old list at vertex j. |
---|
| 705 | ! |
---|
| 706 | ENDDO |
---|
| 707 | ! |
---|
| 708 | ! End lists with -ve values. |
---|
| 709 | ! |
---|
| 710 | iwk(i1) = -1 |
---|
| 711 | iwk(i2) = -1 |
---|
| 712 | ! |
---|
| 713 | ENDIF |
---|
| 714 | ENDDO |
---|
| 715 | ! |
---|
| 716 | ! Copy the polygone in he right order to the output variable |
---|
| 717 | ! |
---|
| 718 | nvert_out = nvert |
---|
| 719 | nvert_tmp=SIZE(poly_out,DIM=1) |
---|
| 720 | IF ( nvert_tmp < nvert_out) THEN |
---|
| 721 | WRITE(*,*) "Ouptput polygone is smaller than the size proposed in the interface." |
---|
| 722 | STOP "polygones_convexhull" |
---|
| 723 | ENDIF |
---|
| 724 | ! |
---|
| 725 | DO i=1,nvert |
---|
| 726 | poly_out(i,:) = poly_in(vertex(i),:) |
---|
| 727 | ENDDO |
---|
| 728 | ! |
---|
| 729 | DEALLOCATE(vertex, iwk, next) |
---|
| 730 | ! |
---|
| 731 | END SUBROUTINE polygones_convexhull |
---|
| 732 | |
---|
| 733 | END MODULE polygones |
---|