source: trunk/TRIANGULATION/triangule_c.pro @ 2

;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
;+
; NAME:triangule_c
;
; PURPOSE:construit le tableau de triangulation.
;
; L''idee est de
; construire une liste de triangles qui relient les points entre
; eux. Ceci est fait automatiquement avec la fonction TRIANGULATE.
;  ICI:
; on tient compte du fait que les points sont disposes sur une grille
; (reguliere ou pas, mais pas destructuree, cad que les points sont
; ecrits suivant une matrice rectangulaire). Un moyen tres simple de
; faire des triangles entre tous les points est alors:
;
;     pour chaque point (i,j) de la matrice -sauf ceux de la derniere
;     ligne et de la derniere colonne- on on appelle le rectangle
;     (i,j) le rectangle forme par les 4 points (i,j), (i+1,j),
;     (i,j+1), (i+1,j+1). Pour tracer tous les triangles, il suffit de
;     tracer les 2 triangles contenus ds les rectangles (i,j)
;
; au passage on remarque que chaque rectangle (i,j) possede 2 diagonales (si
; si faites un dessin c''est vrai), il y a donc 2 choix possibles pour
; chaque rectangles qd on veut le couper en 2 triangles...
; 
; C''est grace a ce choix que l''on va pouvoir tracer les cotes avec
; des angles droits. A chaque angle de cote remarquable par
; l''existance d''un unique point terre ou d''un unique point mer sur
; les 4 cotes d''un rectangle (i,j), il faut couper le rectangle
; suivant la diagonale qui qui passe par le point singulier.
; 
; CATEGORY:pour faire de beaux graphiques masques
;
; CALLING SEQUENCE:res=triangule([mask])
;
; INPUTS:optionnel:mask c''est le tableau 2d qui sevira a masquer le
; champ que l''on tracera apres avec CONTOUR,
; ...TRIANGULATION=triangule(mask)
; si cet argument n''est pas specifie, la function utilise tmask.
;
; KEYWORD PARAMETERS:
;       /REGULIER: specifie que le masque est sur une grille reguliere
;       (utiliser pour la triangulation ds les coupes verticales et
;       des hovmoellers)
;
;       COINMONTE=tableau, pour obtenir le tableau de "coins de terre
;       montant" a traiter avec completecointerre.pro ds la variable
;       tableau plutot que de la faire passer par la variable globale
;       cointerremont.
;
;       COINDESCEND=tableau cf COINMONTE
;
; OUTPUTS:
;       res: tableau 2d (3,nbre de triangles).
;    chaque ligne de res represente les indices des points
;    constituants les sommets d''un triangle.
;    cf. comment on trace les triangles ds dessinetri.pro
;
; COMMON BLOCKS:
;       common.pro different.pro definetri.pro
;
; SIDE EFFECTS:
;
; RESTRICTIONS:les donnees dont un veut ensuite faire le contour
; doivent etre disposees dans une matrice. Par contre dans la matrice,
; la disposition des points peut ne pas etre irreguliere. Si les
; donnees sont disposees completement de facon irreguliere, utiliser
; TRIANGULE.
;
; EXAMPLE:
;
; MODIFICATION HISTORY: Sebastien Masson (smasson@lodyc.jussieu.fr)
;                       26/4/1999
;-
;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
FUNCTION triangule_c, maskentree, COINMONTE = coinmonte, COINDESCEND = coindescend $
                  , REGULIER = regulier, PERIODIQUE = periodique
   tempsun = systime(1)         ; pour key_performance
@common
   ;;
;------------------------------------------------------------
; le masque est donne ou il faut prendre tmask?
;------------------------------------------------------------
;
   msk = maskentree
   taille = size(msk)
   nx = taille[1]
   ny = taille[2]
;
;------------------------------------------------------------
   if n_elements(periodique) EQ 0 then periodique = keyword_set(key_periodique)
   if keyword_set(key_periodique) and keyword_set(periodique) $
    AND NOT keyword_set(regulier) then BEGIN 
      msk = [msk, msk[0, *]]
      nx = nx+1
   ENDIF
;------------------------------------------------------------
; on va trouver la liste des rectangles (i,j) (reperes par leur coin
; en bas a gauche) qu''il faut couper suivant une diagonale descendante
; on appellera cette liste : pts_downward
; 
   pts_downward = 0

; on construit le test qui permet de trouver un tel triangle:
;
;
;       shift(msk,  0, -1)------------shift(msk, -1, -1)
;              |                             |
;              |                             |
;              |                             |
;              |                             |
;             msk---------------------shift(msk, -1,  0)
;
   sum1 = msk+shift(msk, -1, 0)+shift(msk, -1, -1) ;pts qui entourrent le pt en haut a gauche
   sum2 = msk+shift(msk, 0, -1)+shift(msk, -1, -1) ;pts qui entourrent le pt en bas a droite


   tempdeux = systime(1)        ; pour key_performance =2
; pt terre en haut a gauche entoure de pts mer
   liste = where( (4-sum1)*(1-shift(msk, 0, -1)) EQ 1 )
   if liste[0] NE -1 THEN pts_downward = [pts_downward,liste ]
; pt mer en haut a gauche entoure de pts terre
   liste = where( (1-sum1)*shift(msk, 0, -1) EQ 1)
   if liste[0] NE -1 THEN pts_downward = [pts_downward,liste ]
; pt terre en bas a droite entoure de pts mer
   liste = where( (4-sum2)*(1-shift(msk, -1,  0)) EQ 1)
   if liste[0] NE -1 THEN pts_downward = [pts_downward,liste ]
; pt mer en bas a droite entoure de pts terre
   liste = where( (1-sum2)*shift(msk, -1,  0) EQ 1)
   if liste[0] NE -1 THEN pts_downward = [pts_downward,liste ]
   undefine, liste
;
   IF testvar(var = key_performance) EQ 2 THEN $
    print, 'temps triangule: trouver pts_downward', systime(1)-tempdeux
;
   if NOT keyword_set(regulier) then begin
      tempdeux = systime(1)     ; pour key_performance =2
;2 points terre en diagonale montante avec 2 points mer sur la diagonale descendante
      coinmont = where( (1-msk)*(1-shift(msk, -1, -1)) $
                        *(shift(msk, 0, -1)*shift(msk, -1,  0) EQ 1) )
      if coinmont[0] NE -1 THEN pts_downward = [pts_downward, coinmont]
;
      IF testvar(var = key_performance) EQ 2 THEN $
       print, 'temps triangule: trouver coinmont', systime(1)-tempdeux
      tempdeux = systime(1)     ; pour key_performance =2
;
;2 points terre en diagonale descendante avec 2 points mer sur la diagonale montante
      coindesc = where( ((1-shift(msk,  0, -1))*(1-shift(msk, -1, 0)) $
                         *msk*shift(msk, -1, -1) EQ 1) )
;
      IF testvar(var = key_performance) EQ 2 THEN $
       print, 'temps triangule: trouver coindesc', systime(1)-tempdeux
;
   endif
;
   if n_elements(pts_downward) EQ 1 then BEGIN 
      tempdeux = systime(1)     ; pour key_performance =2
;
      triang = definetri(nx, ny)
;
      IF testvar(var = key_performance) EQ 2 THEN $
       print, 'temps triangule: definetri', systime(1)-tempdeux
      coinmont = -1
      coindesc = -1
   ENDIF ELSE BEGIN 
      tempdeux = systime(1)     ; pour key_performance =2
      pts_downward = pts_downward[1:n_elements(pts_downward)-1]
      pts_downward = pts_downward(uniq(pts_downward, sort(pts_downward)))
; aucun rectangle ne peut avoir comme coin en bas a gauche un element
; de la derniere colonne ou de la derniere ligne.
; il faut donc enlever ces points si ils ont ete selectionnes dans
; pts_downward.
      derniere_colonne = (lindgen(ny)+1)*nx-1 
      derniere_ligne = lindgen(nx)+(ny-1)*nx 
      pts_downward =different(pts_downward,derniere_colonne )
      pts_downward =different(pts_downward,derniere_ligne )
      if NOT keyword_set(regulier) then begin
         if coinmont[0] NE -1 then begin
            coinmont =different(coinmont,derniere_colonne )
            coinmont =different(coinmont,derniere_ligne )
         endif
         if coindesc[0] NE -1 then begin
            coindesc =different(coindesc,derniere_colonne )
            coindesc =different(coindesc,derniere_ligne )
         endif
      ENDIF ELSE BEGIN 
         coinmont = -1
         coindesc = -1
      ENDELSE 
      IF testvar(var = key_performance) EQ 2 THEN $
       print, 'temps triangule: menage ds pts_downward coinmont et coindesc', systime(1)-tempdeux
;
      tempdeux = systime(1)     ; pour key_performance =2
      if  pts_downward[0] EQ -1 then triang = definetri(nx, ny) $
      ELSE triang = definetri(nx, ny, pts_downward)
      IF testvar(var = key_performance) EQ 2 THEN $
       print, 'temps triangule: definetri', systime(1)-tempdeux
   ENDELSE 
;------------------------------------------------------------
; on vire les triangles qui ne contiennent que des points terre
;------------------------------------------------------------
;
;  tres bonne idee qui ne marche pas encore a 200% avec IDL 5.2
;  ca devrait aller mieux dans les prochaines versions d''IDL...
;
   tempdeux = systime(1)        ; pour key_performance =2
; on enleve les rectangles qui sont entierement dans la terre
   recdsterre = where((1-msk)*(1-shift(msk, -1, 0))*(1-shift(msk, 0, -1))*(1-shift(msk, -1, -1)) EQ 1)
   IF testvar(var = key_performance) EQ 2 THEN $
    print, 'temps triangule: tous les recdsterre', systime(1)-tempdeux

; en attendant une version qui marche parfaitement, on est contraint
; de faire un nouveau tri:
; il ne faut pas enlever les rectangles qui n''ont qu''un sommet en
; commun.
; t1 = systime(1)
   indice = intarr(nx, ny)
   trimask = intarr(nx, ny)
   trimask[0:nx-2, 0:ny-2] = 1
   IF recdsterre[0] NE -1 then BEGIN 
      tempdeux = systime(1)     ; pour key_performance =2
      indice[recdsterre] = 1
      if NOT keyword_set(regulier) then begin
         vire1 = 0
         vire2 = 0
         while (vire1[0] NE -1 OR vire2[0] NE -1) ne 0 do begin
; vire sont les rectangles qu''il faut retirer de recsterre (en fait
; qu''il faut garder bien qu''ils soient entirement dans la terre)  
            vire1 = where( (indice*shift(indice, -1, -1) $
                            *(1-shift(indice, 0, -1))*(1-shift(indice, -1, 0))*trimask) EQ 1)
            if vire1[0] NE -1 THEN BEGIN 
               indice[vire1] = 0
;               indice[vire1+nx+1] = 0
            endif
            
            vire2 = where( ((1-indice)*(1-shift(indice, -1, -1)) $
                            *shift(indice, 0, -1)*shift(indice, -1, 0)*trimask) EQ 1)
            if vire2[0] NE -1 THEN BEGIN 
               indice[vire2+1] = 0
;               indice[vire2+nx] = 0
            endif
         endwhile
         IF testvar(var = key_performance) EQ 2 THEN $
          print, 'temps triangule: trier les recdsterre', systime(1)-tempdeux
      endif
      indice[*, ny-1] = 1       ; la deriere colonne te la derniere ligne
      indice[nx-1, *] = 1       ; ne peuvent definirde rectangle.
;
      tempdeux = systime(1)     ; pour key_performance =2
      recgarde = where(indice EQ 0)
; on recupere les numeros des triangles que l'' on va garder
      trigarde = 2*[recgarde-recgarde/nx]
      trigarde = [trigarde, trigarde+1]
      trigarde = trigarde[sort(trigarde)]
; 
      triang = triang[*, trigarde]
      IF testvar(var = key_performance) EQ 2 THEN $
       print, 'temps triangule: virer les triangle de la liste', systime(1)-tempdeux
   endif
; print, 'temps tri triangles', systime(1)-t1 
;------------------------------------------------------------
; quand key_periodique eq 1, triang est une liste d''indice d'un
; tableau qui a une colonne de trop.
; il faut ramener ca a la matrice initiale en mettant les indivces de
; la derniere colonne egaux a ceux de la derniere colonne...
;------------------------------------------------------------
   tempdeux = systime(1)        ; pour key_performance =2
   if keyword_set(key_periodique) and keyword_set(periodique) $
    AND NOT keyword_set(regulier) then BEGIN 
      indicey = triang/nx
      indicex = triang-indicey*nx
      nx = nx-1
      liste = where(indicex EQ nx)
      if liste[0] NE -1 then indicex[liste] = 0
      triang = indicex+nx*indicey
      nx = nx+1
      if coinmont[0] NE -1 then begin
         indicey = coinmont/nx
         indicex = coinmont-indicey*nx
         nx = nx-1
         liste = where(indicex EQ nx)
         if liste[0] NE -1 THEN indicex[liste] = 0
         coinmont = indicex+nx*indicey
         nx = nx+1
      endif
      if coindesc[0] NE -1 then begin
         indicey = coindesc/nx
         indicex = coindesc-indicey*nx
         nx = nx-1
         liste = where(indicex EQ nx)
         if liste[0] NE -1 THEN indicex[liste] = 0
         coindesc = indicex+nx*indicey
         nx = nx+1
      endif
   endif
   IF testvar(var = key_performance) EQ 2 THEN $
    print, 'temps triangule: finitions', systime(1)-tempdeux

;------------------------------------------------------------
   if arg_present(coinmonte) THEN coinmonte = coinmont ELSE cointerremont = coinmont
   if arg_present(coindescend) THEN coindescend = coindesc ELSE cointerredesc = coindesc
;
; 
;------------------------------------------------------------


   IF keyword_set(key_performance) THEN print, 'temps triangule', systime(1)-tempsun 

   return, triang

END