;------------------------------------------------------------ ;------------------------------------------------------------ ;------------------------------------------------------------ ;+ ; 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