vw_smooth.F File Reference

#include "implicit_f.inc"

Go to the source code of this file.

Functions/Subroutines
subroutine	vw_smooth (npt, ntarget, x, y)

Function/Subroutine Documentation

vw_smooth()

subroutine vw_smooth	(	integer, intent(inout)	npt,
		integer, intent(in)	ntarget,
		intent(inout)	x,
		intent(inout)	y )

Definition at line 29 of file vw_smooth.F.

C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
C-----------------------------------------------
      INTEGER                     ,INTENT(INOUT)  :: NPT      ! number of input function points
      INTEGER                     ,INTENT(IN)     :: NTARGET  ! target length of abscissa vector
      my_real ,DIMENSION(NPT)     ,INTENT(INOUT)  :: x,y
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER :: I,I0,I1,I2,IDX,IMIN,IPREV,INEXT,ITER,NITER,IZERO
      my_real :: amin,aa,s,t
      INTEGER ,DIMENSION(NPT) :: PREV,NEXT
      my_real ,DIMENSION(NPT) :: area
c-----------------------------------------------
c     smooth input curve using Visvalingam-Whyatt algorithm 
c     reduce number of points from NPT to NTARGET
c     first and last points are never eliminated, as well as (0,0)
c=======================================================================
      niter = npt - ntarget
      IF (niter < 1) RETURN
      
      prev(1)   = 0
      next(1)   = 2
      area(1)   = infinity / ten
      prev(npt) = npt - 1
      next(npt) = 10000000
      area(npt) = infinity / ten 
c
      DO i = 2,npt-1
        i1 = i-1
        i2 = i+1
        prev(i) = i1
        next(i) = i2
        IF (x(i) == zero .and. y(i) == zero) THEN
          area(i) = infinity / ten
        ELSE
          area(i) = x(i1)*y(i)  + x(i)*y(i2) + x(i2)*y(i1)
     .            - x(i1)*y(i2) - x(i)*y(i1) - x(i2)*y(i)
          area(i)  = abs(area(i))
        END IF
      END DO
c----------------------------------------------------      
c     eliminate points with min area until NPT = target number of points
c----------------------------------------------------      
      imin = 1
      DO iter = 1,niter                           
        amin = infinity
        idx  = 1
        DO i=1,npt
          IF (area(i) < amin) THEN
            amin = area(i)
            imin = i
          ENDIF
        END DO      
c     
        iprev = prev(imin)
        inext = next(imin)      
        next(iprev) = inext
        prev(inext) = iprev     
        area(imin)  = infinity
c
        ! recalculate PREV POINT AREA
        i0 = iprev
        i1 = prev(iprev)
        i2 = inext
        IF (i1 > 0) THEN
          aa = x(i0)*y(i0) + x(i0)*y(i2) + x(i2)*y(i1) 
     .       - x(i0)*y(i2) - x(i0)*y(i1) - x(i2)*y(i0)
          area(i0) = abs(aa)
        END IF
c
        ! recalculate NEXT POINT AREA
        i0 = inext                                                                 
        i1 = iprev                                                                 
        i2 = next(inext)                                               
        IF (i2 < npt) THEN                                                   
          aa = x(i0)*y(i0) + x(i0)*y(i2) + x(i2)*y(i1) 
     .       - x(i0)*y(i2) - x(i0)*y(i1) - x(i2)*y(i0)
          area(i0) = abs(aa)
        END IF
      END DO   ! ITER
c------------   
      idx = 0
      DO i=1,npt
        IF (area(i) < infinity) THEN
          idx = idx + 1
          x(idx) = x(i)
          y(idx) = y(i)
        END IF
      END DO
      npt = 0
      izero = 0
      DO i = 1,idx
        s = x(i)
        t = y(i)
        IF (s <  zero .and. t >  zero .or.
     .      s >  zero .and. t <  zero .or.
     .      s == zero .and. t /= zero) THEN
          IF (izero == 1) CONTINUE
          s = zero
          t = zero
          izero = 1
        ELSE IF (s < zero .and. t < zero .and.
     .           x(i+1) > zero .and. y(i+1) > zero) THEN ! .or. 
!     .           S > ZERO .and. T > ZERO .and.
!     .           X(I-1) < ZERO .and. Y(I-1) < ZERO .and. I > 1) THEN
          IF (izero == 1) CONTINUE
          s = zero
          t = zero
          izero = 1
        END IF
        npt = npt + 1
        x(npt) = s
        y(npt) = t
      END DO
c-----------
      RETURN