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

;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
;+
; NAME:
;
; PURPOSE:vire les tableaux qui ne doivent pas etre dessines grace a 2
; tests:
;     1) les coins du tableau doivent etre ds la fenetre
;     2) les clongeurs des cotes des triangfles exprimes en
;     coordonnees normalisesne doivent pas depasser une certaine
;     longueur seuil. 
;
; CATEGORY:
;
; CALLING SEQUENCE:
;
; INPUTS:
;
; KEYWORD PARAMETERS:
;
; OUTPUTS:
;
; COMMON BLOCKS:
;       common.pro
;
; SIDE EFFECTS:
;
; RESTRICTIONS:
;
; EXAMPLE:
;
; MODIFICATION HISTORY: Sebastien Masson (smasson@lodyc.jussieu.fr)
;                       20/2/99
;-
;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
FUNCTION ciseauxtri, triang, glam, gphi, TOUT = tout, _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)         ; pour key_performance
;
   taille = size(glam)
   nx = taille[1]
   ny = taille[2]

   tempdeux = systime(1)        ; pour 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
;
; attention, suivant la projection certains points x ou y peuvent
; devenir NaN (cf. points deriere la terre ds une projection
; orthographique) il faut dans ce cas enlever tous les triangles qui
; contiennent un de ces 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)     ; pour 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)
;
; maintenant on vire les triangles dont un des cotes a une taille
; superieure a 1/seuil du domaine reserve au dessin.
;

   if keyword_set(key_periodic)  then begin
      tempdeux = systime(1)     ; pour 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)     ; pour 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
;
; on supprime les triangles qui sont un peu trop hors du cadre
; valable qd /TOUT est active
;
    if keyword_set(key_irregular) then begin

       tempdeux = systime(1)     ; pour 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
;
;        dernier tri
;
   if keyword_set(trichanged) then BEGIN
;
; il faut verifier que les triangles que l''on garde ne
; forment pas une triangulation dans laquelle 2 triangles ont un
; sommet en commun sans avoir de cote de commun.
;
; on constitue la liste des rectangles que l''on souhaite garder (on
; garde un rectangle des qu''il y a un triangle dedans)
;
; dans definetri, on a construit les triangles tels que le premier et 
; le dernier sommets soient ceux de la diagonale du rectangle definit
; par le maillage.
; pour retouver de quel rectangle provient un triangle, on cherche le
; min de l''indice suivant x et suivant y de chaque triangle. Apres on
; repasse cette indissage suivant x et y en indicage suivant nx*ny
;
      tempdeux = systime(1)     ; pour key_performance =2
; indices en y des rectangles
      indytriang = (triang[0, *]/nx) < (triang[2, *]/nx)
; indices en x des rectangles
      indxtriang0 = triang[0, *]-nx*(triang[0, *]/nx)
      indxtriang2 = triang[2, *]-nx*(triang[2, *]/nx)
      indxmin = indxtriang0 < indxtriang2
; attention dans le cas ou la grille est periodique et ou on a tous
; les points suivant x, les triangles qui assurent la periodicite en x
; ont des indices suivant x qui sont 0 et nx-1. Ils appatienent au
; rectangles de la colonne nx-1 et non de la colonne 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
;
; maintenant qu''on a cette liste on va s''assuter que l''on a pas de
; triangles qui n''ont qu''on sommet en commun.
;
      test = bytarr(nx, ny)
      test[listrect] = 1
      dejavire = 1b-test
;
      tempdeux = systime(1)     ; pour 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 ; on vire le rectangle 
         vire2 = where( ((1-test)*(1-shift(test, -1, -1)) $
                         *shift(test, 0, -1)*shift(test, -1, 0)) EQ 1)
; on vire le rectangle du dessus (meme indice x mais egale a 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: determinationdes rectangles a virer', systime(1)-tempdeux
;
      if avirer[0] NE -1 then begin
      tempdeux = systime(1)     ; pour 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