Re: iso curve extraction

Francois,

First, I am sorry for giving credit to another person on the board
mistakenly. Your help is making me feel optimistic about what I am
trying to do.

Thanks for the description. For orientation, I am looking at local
unit normal pointing from f<0 (inside) to f>0 (outside) for each
segment of the iso curve. This part is somewhat tricky. If my
orientation calculation fails, all of my subsequent calculation will
suffer.

The degenerate case where two curves pass thru one element makes life
difficult. I believe however, that this problem should be overcome by
trying a finer mesh. Is't this corret? (thats why you are computing

Yes it's correct : a fine meshing usually avoids that situation

value at the center right?). The three points in the element cas - One
at the corner and two others is basically same as 4-points and two
segments right?

Right


One possibility is that after I generate a list of segments using your
subroutine, I somehow "clean the list" to have unique set of points
that are connected correctly. I have not idea, however, how to do it.
for example, there could be "disconnected" curves - In my case the
number of disconnected regions is not likely to be more than 2-3.

If anyone has any suggestions about handling orientations and
disconneced curves, please let me know.

You can do the job easily in several successive steps
- building up all segments in calling TITQ4N on all the meshes. The
result could be :
. a number of segments NSEG
. a number of nodes NN equal to 2*NSEG (two nodes by segment)
. the coordinates of the nodes (array XY(2,NN))
- two other integer arrays are necessary to manage segments :
. ISEG(2,NSEG) : the two node indexes defining each segment
. IUSE(NSEG) indicating whether of not each segment has been used to
build up a curve
IUSE must be initialized to 0
- eliminating double nodes
- building up curves

To realize the two last steps, you could use routines similar to the
following ones :

SUBROUTINE eliminate_double(nn,xy,nseg,iseg)

! after building up all isocurve segments, one normally gets :
! - a number of nodes nn equals to 2*nseg
! - the coordinates of the nodes
! To compute the connection between segments, one starts by the
elimination
! of double nodes

INTEGER,INTENT(in) :: nn ! number of nodes (= 2*nseg)
REAL ,INTENT(in) :: xy(2,nn) ! coordinates
INTEGER,INTENT(in) :: nseg ! number of segments built up
by calling TITQ4N
INTEGER,INTENT(out) :: iseg(2,nseg) ! the indexes of the nodes
composing each segment

INTEGER :: i,j,k
REAL :: epsilon=1.e-5

DO i=1,nseg
iseg(1,i)=2*i-1
iseg(2,i)=2*i
ENDDO

DO i=nn,1,-1
DO j=1,i-1
IF((xy(1,i)-xy(1,j))**2+(xy(2,i)-xy(2,j))**2 < epsilon) THEN ! i
and j are identical nodes
DO k=1,nseg ! replacing the node i by the node j
IF(iseg(1,k) == i) iseg(1,k)=j
IF(iseg(2,k) == i) iseg(2,k)=j
ENDDO
EXIT
ENDIF
ENDDO
ENDDO

END SUBROUTINE

SUBROUTINE extract_curve(nseg,iseg,iuse,nc,ic)

INTEGER,INTENT(in) :: nseg ! number of segments
INTEGER,INTENT(in) :: iseg(2,nseg) ! node indexes defining each
segment
INTEGER,INTENT(inout) :: iuse(nseg) ! use index (must be
initialized to 0 by the caller routine)
INTEGER,INTENT(out) :: nc ! number of nodes composing
the iso curve
INTEGER,INTENT(out) :: ic(*) ! node indexes (maximum
size : nseg+1)

INTEGER :: k,j

nc=0

DO k=1,nseg
IF(iuse(k) == 1) CYCLE
nc=2
ic(1)=iseg(1,k)
ic(2)=iseg(2,k)
iuse(k)=1 ! the segment k is declared "used"
DO j=k+1,nseg
IF(iuse(j) == 1) CYCLE
IF(iseg(1,j) == ic(1)) THEN
! the segment j is connected to the beginning of the curve by
its first node
nc=nc+1
ic(2:nc)=ic(1:nc-1) ! shift to allow the insertion
of a value
ic(1)=iseg(2,j)
iuse(j)=1
ELSE IF(iseg(2,j) == ic(1)) THEN
! the segment j is connected to the beginning of the curve by
its second node
nc=nc+1
ic(2:nc)=ic(1:nc-1) ! shift to allow the insertion
of a value
ic(1)=iseg(1,j)
iuse(j)=1
ELSE IF(iseg(1,j) == ic(nc)) THEN
! the segment j is connected to the end of the curve by its
first node
nc=nc+1
ic(nc)=iseg(2,j)
iuse(j)=1
ELSE IF(iseg(2,j) == ic(nc)) THEN
! the segment j is connected to the end of the curve by its
second node
nc=nc+1
ic(nc)=iseg(1,j)
iuse(j)=1
ENDIF
ENDDO
EXIT
ENDDO
END SUBROUTINE

Take care : this is a non tested example of programming !

The routine extract_curve must be called in a loop. One exits the loop
only if NC=0 which means that all segments have been used.

Each call to curve_extract builds up an iso curve. This iso curve is
closed if the first and last node indexes are equal.

Nodes are listed successively but the algorithm does not determine the
orientation (rotation in the direct or reverse sens depending on the
definition of the first selected segment).


.

Relevant Pages

  • Re: Calculus XOR Probability
    ... When you use these broken line segment approximations, ... the segments have endpoints on the curve, ... the curve between those endpoints is parallel to the segment. ... tangent lines at any corner point. ...
    (sci.math)
  • Re: Discrete-time finite-state Markov model
    ... The 2 most common algorithms of fitting polylines to a curve is by ... angle deviation or by chordal deviation. ... To approximate a curve by angle deviation each new line segment is ...
    (comp.dsp)
  • Re: iso curve extraction
    ... segment of the iso curve. ... To compute the connection between segments, ... its first node ...
    (comp.lang.fortran)
  • Re: iso curve extraction
    ... segment of the iso curve. ... To compute the connection between segments, ...
    (comp.lang.fortran)
  • Re: iso curve extraction
    ... CCC The curve has a point in the segment I,I+1 ... I will take some time to see what this routine does. ... dy - increments in x and y directions for the grid ...
    (comp.lang.fortran)
Loading