[2] | 1 | ;------------------------------------------------------------ |
---|
| 2 | ;------------------------------------------------------------ |
---|
| 3 | ;------------------------------------------------------------ |
---|
| 4 | ;+ |
---|
| 5 | ; NAME:triangule_c |
---|
| 6 | ; |
---|
| 7 | ; PURPOSE:construit le tableau de triangulation. |
---|
| 8 | ; |
---|
| 9 | ; L''idee est de |
---|
| 10 | ; construire une liste de triangles qui relient les points entre |
---|
| 11 | ; eux. Ceci est fait automatiquement avec la fonction TRIANGULATE. |
---|
| 12 | ; ICI: |
---|
| 13 | ; on tient compte du fait que les points sont disposes sur une grille |
---|
| 14 | ; (reguliere ou pas, mais pas destructuree, cad que les points sont |
---|
| 15 | ; ecrits suivant une matrice rectangulaire). Un moyen tres simple de |
---|
| 16 | ; faire des triangles entre tous les points est alors: |
---|
| 17 | ; |
---|
| 18 | ; pour chaque point (i,j) de la matrice -sauf ceux de la derniere |
---|
| 19 | ; ligne et de la derniere colonne- on on appelle le rectangle |
---|
| 20 | ; (i,j) le rectangle forme par les 4 points (i,j), (i+1,j), |
---|
| 21 | ; (i,j+1), (i+1,j+1). Pour tracer tous les triangles, il suffit de |
---|
| 22 | ; tracer les 2 triangles contenus ds les rectangles (i,j) |
---|
| 23 | ; |
---|
| 24 | ; au passage on remarque que chaque rectangle (i,j) possede 2 diagonales (si |
---|
| 25 | ; si faites un dessin c''est vrai), il y a donc 2 choix possibles pour |
---|
| 26 | ; chaque rectangles qd on veut le couper en 2 triangles... |
---|
| 27 | ; |
---|
| 28 | ; C''est grace a ce choix que l''on va pouvoir tracer les cotes avec |
---|
| 29 | ; des angles droits. A chaque angle de cote remarquable par |
---|
| 30 | ; l''existance d''un unique point terre ou d''un unique point mer sur |
---|
| 31 | ; les 4 cotes d''un rectangle (i,j), il faut couper le rectangle |
---|
| 32 | ; suivant la diagonale qui qui passe par le point singulier. |
---|
| 33 | ; |
---|
| 34 | ; CATEGORY:pour faire de beaux graphiques masques |
---|
| 35 | ; |
---|
| 36 | ; CALLING SEQUENCE:res=triangule([mask]) |
---|
| 37 | ; |
---|
| 38 | ; INPUTS:optionnel:mask c''est le tableau 2d qui sevira a masquer le |
---|
| 39 | ; champ que l''on tracera apres avec CONTOUR, |
---|
| 40 | ; ...TRIANGULATION=triangule(mask) |
---|
| 41 | ; si cet argument n''est pas specifie, la function utilise tmask. |
---|
| 42 | ; |
---|
| 43 | ; KEYWORD PARAMETERS: |
---|
[29] | 44 | ; |
---|
| 45 | ; /BASIC: specifie que le masque est sur une grille basice |
---|
[2] | 46 | ; (utiliser pour la triangulation ds les coupes verticales et |
---|
| 47 | ; des hovmoellers) |
---|
| 48 | ; |
---|
[29] | 49 | ; /KEEP_CONT: to keep the triangulation even on the continents |
---|
| 50 | ; |
---|
[2] | 51 | ; COINMONTE=tableau, pour obtenir le tableau de "coins de terre |
---|
| 52 | ; montant" a traiter avec completecointerre.pro ds la variable |
---|
| 53 | ; tableau plutot que de la faire passer par la variable globale |
---|
[29] | 54 | ; twin_corners_up. |
---|
[2] | 55 | ; |
---|
| 56 | ; COINDESCEND=tableau cf COINMONTE |
---|
| 57 | ; |
---|
| 58 | ; OUTPUTS: |
---|
| 59 | ; res: tableau 2d (3,nbre de triangles). |
---|
| 60 | ; chaque ligne de res represente les indices des points |
---|
| 61 | ; constituants les sommets d''un triangle. |
---|
| 62 | ; cf. comment on trace les triangles ds dessinetri.pro |
---|
| 63 | ; |
---|
| 64 | ; COMMON BLOCKS: |
---|
| 65 | ; common.pro different.pro definetri.pro |
---|
| 66 | ; |
---|
| 67 | ; SIDE EFFECTS: |
---|
| 68 | ; |
---|
| 69 | ; RESTRICTIONS:les donnees dont un veut ensuite faire le contour |
---|
| 70 | ; doivent etre disposees dans une matrice. Par contre dans la matrice, |
---|
| 71 | ; la disposition des points peut ne pas etre irreguliere. Si les |
---|
| 72 | ; donnees sont disposees completement de facon irreguliere, utiliser |
---|
| 73 | ; TRIANGULE. |
---|
| 74 | ; |
---|
| 75 | ; EXAMPLE: |
---|
| 76 | ; |
---|
| 77 | ; MODIFICATION HISTORY: Sebastien Masson (smasson@lodyc.jussieu.fr) |
---|
| 78 | ; 26/4/1999 |
---|
| 79 | ;- |
---|
| 80 | ;------------------------------------------------------------ |
---|
| 81 | ;------------------------------------------------------------ |
---|
| 82 | ;------------------------------------------------------------ |
---|
[29] | 83 | FUNCTION triangule_c, maskentree, COINMONTE = coinmonte, COINDESCEND = coindescend, BASIC = basic, KEEP_CONT = keep_cont |
---|
[114] | 84 | ; |
---|
| 85 | compile_opt idl2, strictarrsubs |
---|
| 86 | ; |
---|
[2] | 87 | tempsun = systime(1) ; pour key_performance |
---|
[29] | 88 | ;--------------------------------------------------------- |
---|
| 89 | @cm_4mesh |
---|
| 90 | IF NOT keyword_set(key_forgetold) THEN BEGIN |
---|
| 91 | @updatenew |
---|
| 92 | ENDIF |
---|
[2] | 93 | ;------------------------------------------------------------ |
---|
| 94 | ; le masque est donne ou il faut prendre tmask? |
---|
| 95 | ;------------------------------------------------------------ |
---|
| 96 | ; |
---|
| 97 | msk = maskentree |
---|
| 98 | taille = size(msk) |
---|
| 99 | nx = taille[1] |
---|
| 100 | ny = taille[2] |
---|
| 101 | ; |
---|
[29] | 102 | IF n_elements(keep_cont) EQ 0 THEN keep_cont = 1-key_irregular |
---|
[2] | 103 | ;------------------------------------------------------------ |
---|
[29] | 104 | if keyword_set(key_periodic)*(nx EQ jpi) $ |
---|
| 105 | AND NOT keyword_set(basic) then BEGIN |
---|
[2] | 106 | msk = [msk, msk[0, *]] |
---|
| 107 | nx = nx+1 |
---|
| 108 | ENDIF |
---|
| 109 | ;------------------------------------------------------------ |
---|
| 110 | ; on va trouver la liste des rectangles (i,j) (reperes par leur coin |
---|
| 111 | ; en bas a gauche) qu''il faut couper suivant une diagonale descendante |
---|
| 112 | ; on appellera cette liste : pts_downward |
---|
| 113 | ; |
---|
| 114 | pts_downward = 0 |
---|
| 115 | |
---|
| 116 | ; on construit le test qui permet de trouver un tel triangle: |
---|
| 117 | ; |
---|
| 118 | ; |
---|
| 119 | ; shift(msk, 0, -1)------------shift(msk, -1, -1) |
---|
| 120 | ; | | |
---|
| 121 | ; | | |
---|
| 122 | ; | | |
---|
| 123 | ; | | |
---|
| 124 | ; msk---------------------shift(msk, -1, 0) |
---|
| 125 | ; |
---|
| 126 | sum1 = msk+shift(msk, -1, 0)+shift(msk, -1, -1) ;pts qui entourrent le pt en haut a gauche |
---|
| 127 | sum2 = msk+shift(msk, 0, -1)+shift(msk, -1, -1) ;pts qui entourrent le pt en bas a droite |
---|
| 128 | |
---|
| 129 | |
---|
| 130 | tempdeux = systime(1) ; pour key_performance =2 |
---|
| 131 | ; pt terre en haut a gauche entoure de pts mer |
---|
| 132 | liste = where( (4-sum1)*(1-shift(msk, 0, -1)) EQ 1 ) |
---|
| 133 | if liste[0] NE -1 THEN pts_downward = [pts_downward,liste ] |
---|
| 134 | ; pt mer en haut a gauche entoure de pts terre |
---|
| 135 | liste = where( (1-sum1)*shift(msk, 0, -1) EQ 1) |
---|
| 136 | if liste[0] NE -1 THEN pts_downward = [pts_downward,liste ] |
---|
| 137 | ; pt terre en bas a droite entoure de pts mer |
---|
| 138 | liste = where( (4-sum2)*(1-shift(msk, -1, 0)) EQ 1) |
---|
| 139 | if liste[0] NE -1 THEN pts_downward = [pts_downward,liste ] |
---|
| 140 | ; pt mer en bas a droite entoure de pts terre |
---|
| 141 | liste = where( (1-sum2)*shift(msk, -1, 0) EQ 1) |
---|
| 142 | if liste[0] NE -1 THEN pts_downward = [pts_downward,liste ] |
---|
| 143 | undefine, liste |
---|
| 144 | ; |
---|
| 145 | IF testvar(var = key_performance) EQ 2 THEN $ |
---|
| 146 | print, 'temps triangule: trouver pts_downward', systime(1)-tempdeux |
---|
| 147 | ; |
---|
[29] | 148 | if (NOT keyword_set(basic)) OR keyword_set(coinmonte) OR keyword_set(coindescend) then begin |
---|
[2] | 149 | tempdeux = systime(1) ; pour key_performance =2 |
---|
| 150 | ;2 points terre en diagonale montante avec 2 points mer sur la diagonale descendante |
---|
| 151 | coinmont = where( (1-msk)*(1-shift(msk, -1, -1)) $ |
---|
| 152 | *(shift(msk, 0, -1)*shift(msk, -1, 0) EQ 1) ) |
---|
| 153 | if coinmont[0] NE -1 THEN pts_downward = [pts_downward, coinmont] |
---|
| 154 | ; |
---|
| 155 | IF testvar(var = key_performance) EQ 2 THEN $ |
---|
| 156 | print, 'temps triangule: trouver coinmont', systime(1)-tempdeux |
---|
| 157 | tempdeux = systime(1) ; pour key_performance =2 |
---|
| 158 | ; |
---|
| 159 | ;2 points terre en diagonale descendante avec 2 points mer sur la diagonale montante |
---|
| 160 | coindesc = where( ((1-shift(msk, 0, -1))*(1-shift(msk, -1, 0)) $ |
---|
| 161 | *msk*shift(msk, -1, -1) EQ 1) ) |
---|
| 162 | ; |
---|
| 163 | IF testvar(var = key_performance) EQ 2 THEN $ |
---|
| 164 | print, 'temps triangule: trouver coindesc', systime(1)-tempdeux |
---|
| 165 | ; |
---|
[29] | 166 | ENDIF |
---|
[2] | 167 | ; |
---|
| 168 | if n_elements(pts_downward) EQ 1 then BEGIN |
---|
| 169 | tempdeux = systime(1) ; pour key_performance =2 |
---|
| 170 | ; |
---|
| 171 | triang = definetri(nx, ny) |
---|
| 172 | ; |
---|
| 173 | IF testvar(var = key_performance) EQ 2 THEN $ |
---|
| 174 | print, 'temps triangule: definetri', systime(1)-tempdeux |
---|
| 175 | coinmont = -1 |
---|
| 176 | coindesc = -1 |
---|
| 177 | ENDIF ELSE BEGIN |
---|
| 178 | tempdeux = systime(1) ; pour key_performance =2 |
---|
| 179 | pts_downward = pts_downward[1:n_elements(pts_downward)-1] |
---|
[114] | 180 | pts_downward = pts_downward[uniq(pts_downward, sort(pts_downward))] |
---|
[2] | 181 | ; aucun rectangle ne peut avoir comme coin en bas a gauche un element |
---|
| 182 | ; de la derniere colonne ou de la derniere ligne. |
---|
| 183 | ; il faut donc enlever ces points si ils ont ete selectionnes dans |
---|
| 184 | ; pts_downward. |
---|
| 185 | derniere_colonne = (lindgen(ny)+1)*nx-1 |
---|
| 186 | derniere_ligne = lindgen(nx)+(ny-1)*nx |
---|
| 187 | pts_downward =different(pts_downward,derniere_colonne ) |
---|
| 188 | pts_downward =different(pts_downward,derniere_ligne ) |
---|
[29] | 189 | if (NOT keyword_set(basic)) OR keyword_set(coinmonte) OR keyword_set(coindescend) then begin |
---|
[2] | 190 | if coinmont[0] NE -1 then begin |
---|
| 191 | coinmont =different(coinmont,derniere_colonne ) |
---|
| 192 | coinmont =different(coinmont,derniere_ligne ) |
---|
| 193 | endif |
---|
| 194 | if coindesc[0] NE -1 then begin |
---|
| 195 | coindesc =different(coindesc,derniere_colonne ) |
---|
| 196 | coindesc =different(coindesc,derniere_ligne ) |
---|
| 197 | endif |
---|
| 198 | ENDIF ELSE BEGIN |
---|
| 199 | coinmont = -1 |
---|
| 200 | coindesc = -1 |
---|
| 201 | ENDELSE |
---|
| 202 | IF testvar(var = key_performance) EQ 2 THEN $ |
---|
| 203 | print, 'temps triangule: menage ds pts_downward coinmont et coindesc', systime(1)-tempdeux |
---|
| 204 | ; |
---|
| 205 | tempdeux = systime(1) ; pour key_performance =2 |
---|
| 206 | if pts_downward[0] EQ -1 then triang = definetri(nx, ny) $ |
---|
| 207 | ELSE triang = definetri(nx, ny, pts_downward) |
---|
| 208 | IF testvar(var = key_performance) EQ 2 THEN $ |
---|
| 209 | print, 'temps triangule: definetri', systime(1)-tempdeux |
---|
| 210 | ENDELSE |
---|
| 211 | ;------------------------------------------------------------ |
---|
| 212 | ; on vire les triangles qui ne contiennent que des points terre |
---|
| 213 | ;------------------------------------------------------------ |
---|
| 214 | ; |
---|
| 215 | ; tres bonne idee qui ne marche pas encore a 200% avec IDL 5.2 |
---|
| 216 | ; ca devrait aller mieux dans les prochaines versions d''IDL... |
---|
| 217 | ; |
---|
[29] | 218 | if (NOT keyword_set(basic)) AND (NOT keyword_set(keep_cont)) then begin |
---|
| 219 | tempdeux = systime(1) ; pour key_performance =2 |
---|
[2] | 220 | ; on enleve les rectangles qui sont entierement dans la terre |
---|
[29] | 221 | recdsterre = where((1-msk)*(1-shift(msk, -1, 0))*(1-shift(msk, 0, -1))*(1-shift(msk, -1, -1)) EQ 1) |
---|
| 222 | IF testvar(var = key_performance) EQ 2 THEN $ |
---|
| 223 | print, 'temps triangule: tous les recdsterre', systime(1)-tempdeux |
---|
[2] | 224 | |
---|
| 225 | ; en attendant une version qui marche parfaitement, on est contraint |
---|
| 226 | ; de faire un nouveau tri: |
---|
| 227 | ; il ne faut pas enlever les rectangles qui n''ont qu''un sommet en |
---|
| 228 | ; commun. |
---|
| 229 | ; t1 = systime(1) |
---|
[29] | 230 | indice = intarr(nx, ny) |
---|
| 231 | trimask = intarr(nx, ny) |
---|
| 232 | trimask[0:nx-2, 0:ny-2] = 1 |
---|
| 233 | IF recdsterre[0] NE -1 then BEGIN |
---|
| 234 | tempdeux = systime(1) ; pour key_performance =2 |
---|
| 235 | indice[recdsterre] = 1 |
---|
| 236 | ; if NOT keyword_set(basic) then begin |
---|
[2] | 237 | vire1 = 0 |
---|
| 238 | vire2 = 0 |
---|
| 239 | while (vire1[0] NE -1 OR vire2[0] NE -1) ne 0 do begin |
---|
| 240 | ; vire sont les rectangles qu''il faut retirer de recsterre (en fait |
---|
| 241 | ; qu''il faut garder bien qu''ils soient entirement dans la terre) |
---|
| 242 | vire1 = where( (indice*shift(indice, -1, -1) $ |
---|
| 243 | *(1-shift(indice, 0, -1))*(1-shift(indice, -1, 0))*trimask) EQ 1) |
---|
| 244 | if vire1[0] NE -1 THEN BEGIN |
---|
| 245 | indice[vire1] = 0 |
---|
| 246 | ; indice[vire1+nx+1] = 0 |
---|
| 247 | endif |
---|
| 248 | |
---|
| 249 | vire2 = where( ((1-indice)*(1-shift(indice, -1, -1)) $ |
---|
| 250 | *shift(indice, 0, -1)*shift(indice, -1, 0)*trimask) EQ 1) |
---|
| 251 | if vire2[0] NE -1 THEN BEGIN |
---|
| 252 | indice[vire2+1] = 0 |
---|
| 253 | ; indice[vire2+nx] = 0 |
---|
| 254 | endif |
---|
| 255 | endwhile |
---|
| 256 | IF testvar(var = key_performance) EQ 2 THEN $ |
---|
| 257 | print, 'temps triangule: trier les recdsterre', systime(1)-tempdeux |
---|
[29] | 258 | ; endif |
---|
| 259 | indice[*, ny-1] = 1 ; la deriere colonne te la derniere ligne |
---|
| 260 | indice[nx-1, *] = 1 ; ne peuvent definir de rectangle. |
---|
[2] | 261 | ; |
---|
[29] | 262 | tempdeux = systime(1) ; pour key_performance =2 |
---|
| 263 | recgarde = where(indice EQ 0) |
---|
[2] | 264 | ; on recupere les numeros des triangles que l'' on va garder |
---|
[29] | 265 | trigarde = 2*[recgarde-recgarde/nx] |
---|
| 266 | trigarde = transpose(temporary(trigarde)) |
---|
| 267 | trigarde = [trigarde, trigarde+1] |
---|
[2] | 268 | ; |
---|
[29] | 269 | triang = triang[*, temporary(trigarde[*])] |
---|
| 270 | IF testvar(var = key_performance) EQ 2 THEN $ |
---|
| 271 | print, 'temps triangule: virer les triangle de la liste', systime(1)-tempdeux |
---|
| 272 | endif |
---|
[2] | 273 | endif |
---|
| 274 | ; print, 'temps tri triangles', systime(1)-t1 |
---|
| 275 | ;------------------------------------------------------------ |
---|
[29] | 276 | ; quand key_periodic eq 1, triang est une liste d''indice d'un |
---|
[2] | 277 | ; tableau qui a une colonne de trop. |
---|
| 278 | ; il faut ramener ca a la matrice initiale en mettant les indivces de |
---|
| 279 | ; la derniere colonne egaux a ceux de la derniere colonne... |
---|
| 280 | ;------------------------------------------------------------ |
---|
| 281 | tempdeux = systime(1) ; pour key_performance =2 |
---|
[29] | 282 | if keyword_set(key_periodic)*(nx-1 EQ jpi) $ |
---|
| 283 | AND NOT keyword_set(basic) then BEGIN |
---|
[2] | 284 | indicey = triang/nx |
---|
| 285 | indicex = triang-indicey*nx |
---|
| 286 | nx = nx-1 |
---|
| 287 | liste = where(indicex EQ nx) |
---|
| 288 | if liste[0] NE -1 then indicex[liste] = 0 |
---|
| 289 | triang = indicex+nx*indicey |
---|
| 290 | nx = nx+1 |
---|
| 291 | if coinmont[0] NE -1 then begin |
---|
| 292 | indicey = coinmont/nx |
---|
| 293 | indicex = coinmont-indicey*nx |
---|
| 294 | nx = nx-1 |
---|
| 295 | liste = where(indicex EQ nx) |
---|
| 296 | if liste[0] NE -1 THEN indicex[liste] = 0 |
---|
| 297 | coinmont = indicex+nx*indicey |
---|
| 298 | nx = nx+1 |
---|
| 299 | endif |
---|
| 300 | if coindesc[0] NE -1 then begin |
---|
| 301 | indicey = coindesc/nx |
---|
| 302 | indicex = coindesc-indicey*nx |
---|
| 303 | nx = nx-1 |
---|
| 304 | liste = where(indicex EQ nx) |
---|
| 305 | if liste[0] NE -1 THEN indicex[liste] = 0 |
---|
| 306 | coindesc = indicex+nx*indicey |
---|
| 307 | nx = nx+1 |
---|
| 308 | endif |
---|
| 309 | endif |
---|
| 310 | IF testvar(var = key_performance) EQ 2 THEN $ |
---|
| 311 | print, 'temps triangule: finitions', systime(1)-tempdeux |
---|
| 312 | |
---|
| 313 | ;------------------------------------------------------------ |
---|
[29] | 314 | if keyword_set(coinmonte) THEN coinmonte = coinmont ELSE twin_corners_up = coinmont |
---|
| 315 | if keyword_set(coindescend) THEN coindescend = coindesc ELSE twin_corners_dn = coindesc |
---|
[2] | 316 | ;------------------------------------------------------------ |
---|
[29] | 317 | IF NOT keyword_set(key_forgetold) THEN BEGIN |
---|
| 318 | @updateold |
---|
| 319 | ENDIF |
---|
[2] | 320 | |
---|
| 321 | IF keyword_set(key_performance) THEN print, 'temps triangule', systime(1)-tempsun |
---|
| 322 | |
---|
| 323 | return, triang |
---|
| 324 | |
---|
| 325 | END |
---|