source: trunk/SRC/ToBeReviewed/TRIANGULATION/ciseauxtri.pro @ 318

;+
;
; @file_comments
; Delete arrays which do not have to be drawn thanks to 2 tests:
;     1) Corners of the array must be in the window
;     2) Lengths of side of triangles expressed in normalized
;        coordinates must not surpass a sill length.
;
; @categories
;
; @param TRIANG
;
; @param GLAM
;
; @param GPHI
;
; @keyword _EXTRA
; Used to pass keywords
;
; @uses
; common.pro
;
; @history
; Sebastien Masson (smasson\@lodyc.jussieu.fr)
;                       20/2/99
;
; @version
; $Id$
;
;-
;
FUNCTION ciseauxtri, triang, glam, gphi, _EXTRA = ex
;
  compile_opt idl2, strictarrsubs
;
@cm_4mesh
  IF NOT keyword_set(key_forgetold) THEN BEGIN
@updatenew
  ENDIF
;---------------------------------------------------------
  IF NOT keyword_set(key_periodic) AND NOT keyword_set(key_irregular) $
    AND NOT (!map.projection LE 7 AND !map.projection NE 0) $
    AND NOT (!map.projection EQ 14 OR !map.projection EQ 15 $
             OR !map.projection EQ 18) THEN return, triang
;
   tempsun = systime(1)         ; For key_performance
;
   taille = size(glam)
   nx = taille[1]
   ny = taille[2]

   tempdeux = systime(1)        ; For key_performance =2
   z = convert_coord(glam[*],gphi[*],/data,/to_normal)
   x = z[0, *]
   y = z[1, *]
   tempvar = SIZE(TEMPORARY(z)) ; delete z
   IF testvar(var = key_performance) EQ 2 THEN $
    print, 'temps ciseauxtri: convert_coord data to_normal', systime(1)-tempdeux
;
; Beware, following the projection, some points x or y can become NaN
; (see points behind the Earth in an orthographic projection).
; In this case, we have to remove all triangle which contain one of these points.
;
   if (!map.projection LE 7 AND !map.projection NE 0) $
    OR !map.projection EQ 14 OR !map.projection EQ 15 OR !map.projection EQ 18 then begin
      tempdeux = systime(1)     ; For key_performance =2
;
      test = (x*y)[triang]
      test = finite(temporary(test), /nan)
      test = total(temporary(test), 1)
      ind = where(temporary(test) EQ 0)
;
      if ind[0] NE -1 then triang = triang[*, temporary(ind)] ELSE return,  -1
      trichanged = 1b
      IF testvar(var = key_performance) EQ 2 THEN $
       print, 'temps ciseauxtri: recherche points a NAN', systime(1)-tempdeux
   endif
;
   seuil = 5 < (min([nx, ny])-2)
;
; Now we delete triangles whose one side has a size superior to 1/sill from the domain reserved for the drawing.
;

   if keyword_set(key_periodic)  then begin
      tempdeux = systime(1)     ; For key_performance =2
;
      xtri = x[triang]
      xtri = xtri-shift(xtri, 1, 0)
      testx = abs(temporary(xtri)) GT ((!p.position[2]-!p.position[0])/seuil)
      testx = total(temporary(testx), 1)
;
      ytri = y[triang]
      ytri = ytri-shift(ytri, 1, 0)
      testy = abs(temporary(ytri)) GT ((!p.position[3]-!p.position[1])/seuil)
      testy = total(temporary(testy), 1)
;
      test = temporary(testx)+temporary(testy)
      ind=where(temporary(test) EQ 0)
;
      IF testvar(var = key_performance) EQ 2 THEN $
       print, 'temps ciseauxtri: trouver les triangles trop grands', systime(1)-tempdeux
;
      tempdeux = systime(1)     ; For key_performance =2
      if ind[0] NE -1 then triang = triang[*, temporary(ind)] ELSE return,  -1
      trichanged = 1b
      IF testvar(var = key_performance) EQ 2 THEN $
       print, 'temps ciseauxtri: virer les triangles trop grands', systime(1)-tempdeux
   endif
;
; We delete all triangle which are out of the valid frame
;
    if keyword_set(key_irregular) then begin

       tempdeux = systime(1)     ; For key_performance =2
       xtri = x[triang]
       test1 = xtri GE !p.position[0]
       test2 = xtri LE !p.position[2]
       undefine, xtri
       testx = temporary(test1)*temporary(test2)
       testx = total(temporary(testx), 1)
;
       ytri = y[triang]
       test1 = ytri GE !p.position[1]
       test2 = ytri LE !p.position[3]
       undefine, ytri
       testy = temporary(test1)*temporary(test2)
       testy = total(temporary(testy), 1)
;
       test = temporary(testx)*temporary(testy);
;
       ind=where(temporary(test) NE 0)
;
       if ind[0] NE -1 then triang = triang[*, temporary(ind)] ELSE return,  -1
       trichanged = 1b
       IF testvar(var = key_performance) EQ 2 THEN $
        print, 'temps ciseauxtri: virer les triangles hors du cadre', systime(1)-tempdeux
    ENDIF
;
;        Last sort.
;
   if keyword_set(trichanged) then BEGIN
;
; We have to check that triangles we keep do not form a triangulation
; in which 2 triangles have a common summit without have a common side.
;
; We constitute the list of rectangles we want to keep (we keep every
; rectangle containing a triangle)
;
; In definetri, we have construct triangles just so the first and the
; last summit are those of the diagonale of the rectangle defined by
; the mesh size. To find from which rectangle a triangle comes, we look
; for the min of the index following x and following y of each triangle.
; Then we go by again this indexion following x and y in an indexion
; following nx*ny/
;
      tempdeux = systime(1)     ; For key_performance =2
; y indexes of rectangles
      indytriang = (triang[0, *]/nx) < (triang[2, *]/nx)
; x indexes of rectangles
      indxtriang0 = triang[0, *]-nx*(triang[0, *]/nx)
      indxtriang2 = triang[2, *]-nx*(triang[2, *]/nx)
      indxmin = indxtriang0 < indxtriang2
; Beware in the case where the grid is periodic and where we have all points
; following x, triangles which assure the periodicity in x have indexes
; following x which are 0 and nx-1. They belong to rectangles of the column
; nx-1 and not to column 0.
      if keyword_set(key_periodic) AND nx EQ jpi then BEGIN
      indxmax = indxtriang0 > indxtriang2
      indxtriang = indxmin + (nx-1)*(indxmin EQ 0 AND indxmax EQ (nx-1))
      ENDIF ELSE indxtriang = indxmin
      listrect = nx*indytriang+indxtriang
      IF testvar(var = key_performance) EQ 2 THEN $
       print, 'temps ciseauxtri: liste des rectangles', systime(1)-tempdeux
;
; Now we have this list, we will make sure that we do not have triangles
; with only a common summit.
;
      test = bytarr(nx, ny)
      test[listrect] = 1
      dejavire = 1b-test
;
      tempdeux = systime(1)     ; For key_performance =2
      vire1 = 0
      vire2 = 0
      while (vire1[0] NE -1 OR vire2[0] NE -1) ne 0 do begin
         vire1 = where( (test*shift(test, -1, -1) $
                         *(1-shift(test, 0, -1))*(1-shift(test, -1, 0))) EQ 1)
         if vire1[0] NE -1 THEN test[vire1] = 0 ; We delete the rectangle
         vire2 = where( ((1-test)*(1-shift(test, -1, -1)) $
                         *shift(test, 0, -1)*shift(test, -1, 0)) EQ 1)
; We delete the top rectangle (same x index but equal to 1)
         if vire2[0] NE -1 THEN test[vire2+nx] = 0
      ENDWHILE
;stop
      test = test+temporary(dejavire)
;
      avirer = where(temporary(test) EQ 0)
      IF testvar(var = key_performance) EQ 2 THEN $
       print, 'temps ciseauxtri: determination des rectangles a virer', systime(1)-tempdeux
;
      if avirer[0] NE -1 then begin
      tempdeux = systime(1)     ; For key_performance =2
      indnx = n_elements(listrect)
      indny = n_elements(avirer)
      ind = listrect[*]#replicate(1l, indny)
      ind = ind EQ replicate(1, indnx)#avirer[*]
      if indny GT 1 then ind = total(ind, 2)
      ind = where(ind EQ 0)
      if ind[0] NE -1 then triang = triang[*, ind] ELSE return,  -1
      endif
      IF testvar(var = key_performance) EQ 2 THEN $
       print, 'temps ciseauxtri: derniere retouche de la triangulation', systime(1)-tempdeux
   endif
;
;
   if keyword_set(key_performance) THEN print, 'temps ciseauxtri', systime(1)-tempsun
;
   return,  triang
end