source: trunk/SRC/ToBeReviewed/TRIANGULATION/tracemask.pro @ 373

;+
;
; @file_comments
; Draw contours of a mask
;
; @categories
; Utilities
;
; @param MASKENTREE {in}{required}
; 2d array specifying the mask
;
; @param XIN {in}{required}
; 2d array specifying longitude coordinates.
;
; @param YIN {in}{required}
; 2d array specifying latitude coordinates.
;
; @keyword COAST_COLOR {default=0}
; The color of the coastline.
; Default is black (0).
;
; @keyword COAST_THICK {default=1}
; The thick of the trait to trace continents
;
; @keyword OVERPLOT
; To do a plot over an other one.
;
; @keyword _EXTRA
; Used to pass keywords
;
; @uses
; <pro>common</pro>
;
; @history
; Sebastien Masson (smasson\@lodyc.jussieu.fr)
;
; @version
; $Id$
;
;-
PRO tracemask, maskentree, xin, yin $
             , COAST_COLOR=coast_color, COAST_THICK=coast_thick $
             , OVERPLOT=overplot, _EXTRA=ex
;
  compile_opt idl2, strictarrsubs
;
   if keyword_set(overplot) then return
;---------------------------------------------------------
@cm_general
  IF NOT keyword_set(key_forgetold) THEN BEGIN
@updatekwd
  ENDIF
;---------------------------------------------------------
   tempsun = systime(1)         ; For key_performance
; We avoid edging problems:
   tempdeux = systime(1)        ; For key_performance =2
   tailleentree = size(maskentree)
   nx = tailleentree[1]+1
   ny = tailleentree[2]+1
; we check the input axis
  IF n_elements(xin) EQ 0 THEN xentree = findgen(nx-1) ELSE xentree = xin
  IF (size(xentree))[0] EQ 1 THEN xentree = xentree#replicate(1,ny-1)
  IF n_elements(yin) EQ 0 THEN yentree = findgen(ny-1) ELSE yentree = yin
  IF (size(yentree))[0] EQ 1 THEN yentree = replicate(1,nx-1)#yentree
; We enlarge the mask by 1 column to the left an d1 line to the bottom
   mask = intarr(tailleentree[1]+1, tailleentree[2]+1)
   mask[1:tailleentree[1], 1:tailleentree[2]] = maskentree
; The 2 first columns are identical.
   mask[0, 1:tailleentree[2]] = maskentree[0, *]
; The 2 first lines are identical.
   mask[1:tailleentree[1], 0] = maskentree[*, 0]
; We calculate the position following x of points which will serve to trace the mask. They are situated between each points of the mask, exept for the last column we can not calculate and so we put at max (!x.range).
   xrange = !x.range[sort(!x.range)] ; if REVERSE_X is used
   xentree = .5*(xentree+shift(xentree, -1, 0))
   IF not keyword_set(overplot) THEN xentree[nx-2, *] = xrange[1] $
   ELSE xentree[nx-2, *] = xentree[nx-3, *]
; we sill
   xentree = xrange[0] > xentree < xrange[1]
; we enlarge the array
   xf = fltarr(nx, ny)
   xf[1:nx-1, 1:ny-1] = xentree
   IF not keyword_set(overplot) THEN xf[0, *] = xrange[0] $
   ELSE xf[0, *] = xf[1, *]
   xf[1:nx-1, 0] = xentree[*, 0]
;
   yinverse = yentree[0, 0] GT yentree[0, ny-2]
   yrange = !y.range[sort(!y.range)]
   yentree = .5*(yentree+shift(yentree, 0, -1))
   IF not keyword_set(overplot) THEN BEGIN
      if yinverse then yentree[*, ny-2] = yrange[0] ELSE yentree[*, ny-2] = yrange[1]
   ENDIF ELSE yentree[*, ny-2] = yentree[*, ny-3]
   yentree = yrange[0] > yentree < yrange[1]
   yf = fltarr(nx, ny)
   yf[1:nx-1, 1:ny-1] = yentree
   yf[0, 1:ny-1] = yentree[0, *]
   IF not keyword_set(overplot) THEN BEGIN
      if yinverse then yf[*, 0] = yrange[1] ELSE yf[*, 0] = yrange[0]
   ENDIF ELSE yentree[*, 0] = yentree[*, 1]
;
   IF testvar(var = key_performance) EQ 2 THEN $
    print, 'temps tracemask: determination du mask et de ses coordonnees', systime(1)-tempdeux
;
; We trace vertical segments:
;
   tempdeux = systime(1)        ; For key_performance =2
   liste = where(mask+shift(mask, -1, 0) EQ 1)
   IF liste[0] NE -1 THEN BEGIN
; We recuperate lx and ly which are indexes in a 2d array of points given by list
      ly = liste/nx & lx = temporary(liste)-nx*ly
      indice = where(ly NE 0) ; We do not take points concerning
; the first line because in this case, the point j-1 is not defined
      if indice[0] NE -1 then begin
         lx = lx[indice] & ly = ly[temporary(indice)]
         IF testvar(var = key_performance) EQ 2 THEN $
          print, 'temps tracemask: liste traits verticaux', systime(1)-tempdeux
         tempdeux = systime(1)  ; For key_performance =2
; loop on concerned points and drawing of the segment.
; comments: we use plots instead of plot because plots is faster.
         for pt = 0L, n_elements(lx)-1 do BEGIN
            i = lx[pt] & j = ly[pt]
            plots, [xf[i, j-1], xf[i, j]], [yf[i, j-1], yf[i, j]] $
              , color = coast_color, thick = coast_thick, _extra = ex
            if pt LT 5 then begin
            endif
         endfor
         IF testvar(var = key_performance) EQ 2 THEN $
          print, 'temps tracemask: trace traits verticaux', systime(1)-tempdeux
      endif
   ENDIF
;
; We trace horizontal segments:
;
   tempdeux = systime(1)        ; For key_performance =2
   liste = where(mask+shift(mask, 0, -1) EQ 1)
   IF liste[0] NE -1 THEN BEGIN
      ly = liste/nx & lx = temporary(liste)-nx*ly
      indice = where(lx NE 0)   ; We do not take point sof the first column.
      if indice[0] EQ -1 then return
      lx = lx[indice] & ly = ly[temporary(indice)]
      IF testvar(var = key_performance) EQ 2 THEN $
       print, 'temps tracemask: liste traits horizontaux', systime(1)-tempdeux
      tempdeux = systime(1)     ; For key_performance =2
      for pt = 0L, n_elements(lx)-1 do BEGIN
         i = lx[pt] & j = ly[pt]
         plots, [xf[i-1, j], xf[i, j]], [yf[i-1, j], yf[i, j]] $
           , color = coast_color, thick = coast_thick, _extra = ex
      endfor
      IF testvar(var = key_performance) EQ 2 THEN $
       print, 'temps tracemask: trace traits horizontaux', systime(1)-tempdeux
   endif
;
;
   if keyword_set(key_performance) THEN print, 'temps tracemask', systime(1)-tempsun

   return
end