[2] | 1 | ;------------------------------------------------------------ |
---|
| 2 | ;------------------------------------------------------------ |
---|
| 3 | ;------------------------------------------------------------ |
---|
| 4 | ;+ |
---|
| 5 | ; |
---|
[150] | 6 | ; @file_comments |
---|
| 7 | ; Construct the triangulation array. |
---|
[2] | 8 | ; |
---|
[226] | 9 | ; The idea is: construct a list of triangle which link points between them. |
---|
[150] | 10 | ; This is automatically done by the function TRIANGULATE |
---|
| 11 | ; Here: |
---|
[226] | 12 | ; we consider the fact that points are disposed on a grid (regular or not, |
---|
| 13 | ; but not unstructured, that is to say that points are written following a |
---|
| 14 | ; rectangular matrix). A easy way to do triangles between all points is then: |
---|
[2] | 15 | ; |
---|
[163] | 16 | ; for each point (i,j) of the matrix -except those of the last line and of |
---|
[150] | 17 | ; the last column- we call rectangle (i,j) the rectangle made of the four |
---|
| 18 | ; points (i,j), (i+1,j), (i,j+1), (i+1,j+1). To trace all triangle, we just |
---|
| 19 | ; have to trace the 2 triangles contained in rectangles (i,j) |
---|
[2] | 20 | ; |
---|
[150] | 21 | ; We notice that each rectangle (i,j) have 2 diagonals (it is true... Make a |
---|
| 22 | ; drawing to make sure!!), so there are two possible choice for each rectangle |
---|
| 23 | ; we want to cut in 2 triangles... |
---|
[226] | 24 | ; |
---|
[150] | 25 | ; It is thanks to this choice that we will be able to trace coast with right |
---|
| 26 | ; angles. At each angle of coast remarkable by the existence of an unique land |
---|
| 27 | ; point or of an unique ocean point on one of the four summit of a rectangle (i,j), |
---|
| 28 | ; we have to cut the rectangle following the diagonal passing by this point. |
---|
[226] | 29 | ; |
---|
[150] | 30 | ; @categories |
---|
[157] | 31 | ; Graphics |
---|
[2] | 32 | ; |
---|
[163] | 33 | ; @param MASKENTREE {in}{optional}{type=2d array} |
---|
[226] | 34 | ; It is a 2d array which will serve to mask the field we will trace after with CONTOUR, |
---|
[2] | 35 | ; ...TRIANGULATION=triangule(mask) |
---|
[150] | 36 | ; If this argument is not specified, the function use tmask |
---|
[2] | 37 | ; |
---|
[150] | 38 | ; @keyword BASIC |
---|
| 39 | ; Specify that the mask is on a basic grid (use the triangulation for vertical cuts and hovmoellers) |
---|
[29] | 40 | ; |
---|
[150] | 41 | ; @keyword KEEP_CONT |
---|
| 42 | ; To keep the triangulation even on the continents |
---|
[2] | 43 | ; |
---|
[163] | 44 | ; @keyword COINMONTE {type=array} |
---|
[226] | 45 | ; To obtain the array of "ascending land corner" to be treated with |
---|
| 46 | ; completecointerre.pro in the variable array instead of make it pass by the global |
---|
[150] | 47 | ; variable twin_corners_up. |
---|
[29] | 48 | ; |
---|
[163] | 49 | ; @keyword COINDESCEND {type=array} |
---|
| 50 | ; See COINMONTE |
---|
[2] | 51 | ; |
---|
[150] | 52 | ; @returns |
---|
| 53 | ; res: tableau 2d (3,nbre de triangles). |
---|
[226] | 54 | ; Each line of res represent indexes of points constituting summits of a triangle. |
---|
[150] | 55 | ; See how we trace triangles in definetri.pro |
---|
[2] | 56 | ; |
---|
[150] | 57 | ; @uses |
---|
| 58 | ; common.pro |
---|
| 59 | ; different.pro |
---|
| 60 | ; definetri.pro |
---|
[2] | 61 | ; |
---|
[150] | 62 | ; @restrictions |
---|
[226] | 63 | ; Datas whose we want to do the contour must be disposed in a matrix. |
---|
| 64 | ; On the other hand, in the matrix, the points's arrangement can not be |
---|
[150] | 65 | ; irregular. If it is, use TRIANGULE. |
---|
[2] | 66 | ; |
---|
[150] | 67 | ; @history |
---|
[157] | 68 | ; Sebastien Masson (smasson\@lodyc.jussieu.fr) |
---|
[150] | 69 | ; 26/4/1999 |
---|
[2] | 70 | ; |
---|
[150] | 71 | ; @version |
---|
| 72 | ; $Id$ |
---|
[2] | 73 | ; |
---|
[150] | 74 | ; @todo |
---|
[226] | 75 | ; seb L.267->268 je ne pense pas que ce soit ce que tu voulais dire mais |
---|
[150] | 76 | ; c'est la traduction de ce qu'il y avait écrit. Correction si besoin. |
---|
[2] | 77 | ;- |
---|
| 78 | ;------------------------------------------------------------ |
---|
| 79 | ;------------------------------------------------------------ |
---|
| 80 | ;------------------------------------------------------------ |
---|
[29] | 81 | FUNCTION triangule_c, maskentree, COINMONTE = coinmonte, COINDESCEND = coindescend, BASIC = basic, KEEP_CONT = keep_cont |
---|
[114] | 82 | ; |
---|
[150] | 83 | compile_opt idl2, strictarrsubs |
---|
[114] | 84 | ; |
---|
[150] | 85 | tempsun = systime(1) ; For key_performance |
---|
[29] | 86 | ;--------------------------------------------------------- |
---|
| 87 | @cm_4mesh |
---|
[150] | 88 | IF NOT keyword_set(key_forgetold) THEN BEGIN |
---|
[29] | 89 | @updatenew |
---|
[150] | 90 | ENDIF |
---|
[2] | 91 | ;------------------------------------------------------------ |
---|
[150] | 92 | ; Is the mask given or do we have to take tmask? |
---|
[2] | 93 | ;------------------------------------------------------------ |
---|
| 94 | ; |
---|
[150] | 95 | msk = maskentree |
---|
| 96 | taille = size(msk) |
---|
| 97 | nx = taille[1] |
---|
| 98 | ny = taille[2] |
---|
[2] | 99 | ; |
---|
[150] | 100 | IF n_elements(keep_cont) EQ 0 THEN keep_cont = 1-key_irregular |
---|
[2] | 101 | ;------------------------------------------------------------ |
---|
[150] | 102 | if keyword_set(key_periodic)*(nx EQ jpi) $ |
---|
[226] | 103 | AND NOT keyword_set(basic) then BEGIN |
---|
[150] | 104 | msk = [msk, msk[0, *]] |
---|
| 105 | nx = nx+1 |
---|
| 106 | ENDIF |
---|
[2] | 107 | ;------------------------------------------------------------ |
---|
[226] | 108 | ; We will find the list of rectangles (i,j)(located by their left |
---|
| 109 | ; bottom corner) we have to cut following a descendant diagonal. |
---|
[150] | 110 | ; We will call this list : pts_downward |
---|
[226] | 111 | ; |
---|
[150] | 112 | pts_downward = 0 |
---|
[2] | 113 | |
---|
[150] | 114 | ; We construct the test which allow to find this triangle : |
---|
[2] | 115 | ; |
---|
| 116 | ; |
---|
| 117 | ; shift(msk, 0, -1)------------shift(msk, -1, -1) |
---|
| 118 | ; | | |
---|
| 119 | ; | | |
---|
| 120 | ; | | |
---|
| 121 | ; | | |
---|
| 122 | ; msk---------------------shift(msk, -1, 0) |
---|
| 123 | ; |
---|
[150] | 124 | sum1 = msk+shift(msk, -1, 0)+shift(msk, -1, -1) ;points which surround the left top point. |
---|
| 125 | sum2 = msk+shift(msk, 0, -1)+shift(msk, -1, -1) ;points which surround the right bottom point. |
---|
[2] | 126 | |
---|
| 127 | |
---|
[150] | 128 | tempdeux = systime(1) ; For key_performance =2 |
---|
| 129 | ; The left top land point surrounded by ocean points |
---|
| 130 | liste = where( (4-sum1)*(1-shift(msk, 0, -1)) EQ 1 ) |
---|
| 131 | if liste[0] NE -1 THEN pts_downward = [pts_downward,liste ] |
---|
| 132 | ; The left top ocean point surrounded by land points |
---|
| 133 | liste = where( (1-sum1)*shift(msk, 0, -1) EQ 1) |
---|
| 134 | if liste[0] NE -1 THEN pts_downward = [pts_downward,liste ] |
---|
| 135 | ; The right bottom land point surrounded by ocean points |
---|
| 136 | liste = where( (4-sum2)*(1-shift(msk, -1, 0)) EQ 1) |
---|
| 137 | if liste[0] NE -1 THEN pts_downward = [pts_downward,liste ] |
---|
| 138 | ; The right bottom ocean point surrounded by land points |
---|
| 139 | liste = where( (1-sum2)*shift(msk, -1, 0) EQ 1) |
---|
| 140 | if liste[0] NE -1 THEN pts_downward = [pts_downward,liste ] |
---|
| 141 | undefine, liste |
---|
[2] | 142 | ; |
---|
[150] | 143 | IF testvar(var = key_performance) EQ 2 THEN $ |
---|
| 144 | print, 'temps triangule: trouver pts_downward', systime(1)-tempdeux |
---|
[2] | 145 | ; |
---|
[150] | 146 | if (NOT keyword_set(basic)) OR keyword_set(coinmonte) OR keyword_set(coindescend) then begin |
---|
| 147 | tempdeux = systime(1) ; For key_performance =2 |
---|
| 148 | ;2 land points in ascendant diagonal with 2 ocean points in descendant diagonal. |
---|
| 149 | coinmont = where( (1-msk)*(1-shift(msk, -1, -1)) $ |
---|
| 150 | *(shift(msk, 0, -1)*shift(msk, -1, 0) EQ 1) ) |
---|
| 151 | if coinmont[0] NE -1 THEN pts_downward = [pts_downward, coinmont] |
---|
[2] | 152 | ; |
---|
[150] | 153 | IF testvar(var = key_performance) EQ 2 THEN $ |
---|
| 154 | print, 'temps triangule: trouver coinmont', systime(1)-tempdeux |
---|
| 155 | tempdeux = systime(1) ; pour key_performance =2 |
---|
[2] | 156 | ; |
---|
[150] | 157 | coindesc = where( ((1-shift(msk, 0, -1))*(1-shift(msk, -1, 0)) $ |
---|
| 158 | *msk*shift(msk, -1, -1) EQ 1) ) |
---|
[2] | 159 | ; |
---|
[150] | 160 | ;2 land points in descendant diagonal with 2 ocean points in ascendant diagonal. |
---|
| 161 | IF testvar(var = key_performance) EQ 2 THEN $ |
---|
| 162 | print, 'temps triangule: trouver coindesc', systime(1)-tempdeux |
---|
[2] | 163 | ; |
---|
[150] | 164 | ENDIF |
---|
[2] | 165 | ; |
---|
[226] | 166 | if n_elements(pts_downward) EQ 1 then BEGIN |
---|
[150] | 167 | tempdeux = systime(1) ; For key_performance =2 |
---|
[2] | 168 | ; |
---|
[150] | 169 | triang = definetri(nx, ny) |
---|
[2] | 170 | ; |
---|
[150] | 171 | IF testvar(var = key_performance) EQ 2 THEN $ |
---|
| 172 | print, 'temps triangule: definetri', systime(1)-tempdeux |
---|
| 173 | coinmont = -1 |
---|
| 174 | coindesc = -1 |
---|
[226] | 175 | ENDIF ELSE BEGIN |
---|
[150] | 176 | tempdeux = systime(1) ; For key_performance =2 |
---|
| 177 | pts_downward = pts_downward[1:n_elements(pts_downward)-1] |
---|
| 178 | pts_downward = pts_downward[uniq(pts_downward, sort(pts_downward))] |
---|
[226] | 179 | ; None rectangle can have an element of the last column or of the |
---|
[150] | 180 | ; last line as left bottom corner. |
---|
| 181 | ; so we have to remove these points if they has been selected in pts_downward. |
---|
[226] | 182 | derniere_colonne = (lindgen(ny)+1)*nx-1 |
---|
| 183 | derniere_ligne = lindgen(nx)+(ny-1)*nx |
---|
[150] | 184 | pts_downward =different(pts_downward,derniere_colonne ) |
---|
| 185 | pts_downward =different(pts_downward,derniere_ligne ) |
---|
| 186 | if (NOT keyword_set(basic)) OR keyword_set(coinmonte) OR keyword_set(coindescend) then begin |
---|
| 187 | if coinmont[0] NE -1 then begin |
---|
[2] | 188 | coinmont =different(coinmont,derniere_colonne ) |
---|
| 189 | coinmont =different(coinmont,derniere_ligne ) |
---|
[150] | 190 | endif |
---|
| 191 | if coindesc[0] NE -1 then begin |
---|
[2] | 192 | coindesc =different(coindesc,derniere_colonne ) |
---|
| 193 | coindesc =different(coindesc,derniere_ligne ) |
---|
[150] | 194 | endif |
---|
[226] | 195 | ENDIF ELSE BEGIN |
---|
[150] | 196 | coinmont = -1 |
---|
| 197 | coindesc = -1 |
---|
[226] | 198 | ENDELSE |
---|
[150] | 199 | IF testvar(var = key_performance) EQ 2 THEN $ |
---|
| 200 | print, 'temps triangule: menage ds pts_downward coinmont et coindesc', systime(1)-tempdeux |
---|
[2] | 201 | ; |
---|
[150] | 202 | tempdeux = systime(1) ; For key_performance =2 |
---|
| 203 | if pts_downward[0] EQ -1 then triang = definetri(nx, ny) $ |
---|
| 204 | ELSE triang = definetri(nx, ny, pts_downward) |
---|
| 205 | IF testvar(var = key_performance) EQ 2 THEN $ |
---|
| 206 | print, 'temps triangule: definetri', systime(1)-tempdeux |
---|
[226] | 207 | ENDELSE |
---|
[2] | 208 | ;------------------------------------------------------------ |
---|
[150] | 209 | ; We delete land points which only contain land points. |
---|
[2] | 210 | ;------------------------------------------------------------ |
---|
| 211 | ; |
---|
[226] | 212 | ; |
---|
[150] | 213 | if (NOT keyword_set(basic)) AND (NOT keyword_set(keep_cont)) then begin |
---|
| 214 | tempdeux = systime(1) ; For key_performance =2 |
---|
| 215 | ; We delete rectangles which are entirely in the land. |
---|
| 216 | recdsterre = where((1-msk)*(1-shift(msk, -1, 0))*(1-shift(msk, 0, -1))*(1-shift(msk, -1, -1)) EQ 1) |
---|
| 217 | IF testvar(var = key_performance) EQ 2 THEN $ |
---|
| 218 | print, 'temps triangule: tous les recdsterre', systime(1)-tempdeux |
---|
[2] | 219 | |
---|
[150] | 220 | ; We do an other sort : |
---|
| 221 | ; We have to do not remove rectangles which only have one common summit. |
---|
[2] | 222 | ; t1 = systime(1) |
---|
[150] | 223 | indice = intarr(nx, ny) |
---|
| 224 | trimask = intarr(nx, ny) |
---|
| 225 | trimask[0:nx-2, 0:ny-2] = 1 |
---|
[226] | 226 | IF recdsterre[0] NE -1 then BEGIN |
---|
[150] | 227 | tempdeux = systime(1) ; For key_performance =2 |
---|
| 228 | indice[recdsterre] = 1 |
---|
| 229 | if NOT keyword_set(basic) then begin |
---|
| 230 | vire1 = 0 |
---|
| 231 | vire2 = 0 |
---|
| 232 | while (vire1[0] NE -1 OR vire2[0] NE -1) ne 0 do begin |
---|
[226] | 233 | ; Delete rectangles we have to remove from recsterre (in fact those we have |
---|
| 234 | ; to keep although they are entirely in the land). |
---|
[150] | 235 | vire1 = where( (indice*shift(indice, -1, -1) $ |
---|
| 236 | *(1-shift(indice, 0, -1))*(1-shift(indice, -1, 0))*trimask) EQ 1) |
---|
[226] | 237 | if vire1[0] NE -1 THEN BEGIN |
---|
[150] | 238 | indice[vire1] = 0 |
---|
[2] | 239 | ; indice[vire1+nx+1] = 0 |
---|
[150] | 240 | endif |
---|
[226] | 241 | |
---|
[150] | 242 | vire2 = where( ((1-indice)*(1-shift(indice, -1, -1)) $ |
---|
| 243 | *shift(indice, 0, -1)*shift(indice, -1, 0)*trimask) EQ 1) |
---|
[226] | 244 | if vire2[0] NE -1 THEN BEGIN |
---|
[150] | 245 | indice[vire2+1] = 0 |
---|
[2] | 246 | ; indice[vire2+nx] = 0 |
---|
[150] | 247 | endif |
---|
| 248 | endwhile |
---|
| 249 | IF testvar(var = key_performance) EQ 2 THEN $ |
---|
| 250 | print, 'temps triangule: trier les recdsterre', systime(1)-tempdeux |
---|
| 251 | endif |
---|
| 252 | indice[*, ny-1] = 1 ; The last column and the last line |
---|
| 253 | indice[nx-1, *] = 1 ; can not define any rectangle. |
---|
[2] | 254 | ; |
---|
[150] | 255 | tempdeux = systime(1) ; For key_performance =2 |
---|
| 256 | recgarde = where(indice EQ 0) |
---|
| 257 | ; We recuperate numbers of triangles we will keep. |
---|
| 258 | trigarde = 2*[recgarde-recgarde/nx] |
---|
| 259 | trigarde = transpose(temporary(trigarde)) |
---|
| 260 | trigarde = [trigarde, trigarde+1] |
---|
[226] | 261 | ; |
---|
[150] | 262 | triang = triang[*, temporary(trigarde[*])] |
---|
| 263 | IF testvar(var = key_performance) EQ 2 THEN $ |
---|
[29] | 264 | print, 'temps triangule: virer les triangle de la liste', systime(1)-tempdeux |
---|
[150] | 265 | endif |
---|
| 266 | endif |
---|
[226] | 267 | ; print, 'temps tri triangles', systime(1)-t1 |
---|
[2] | 268 | ;------------------------------------------------------------ |
---|
[226] | 269 | ; When key_periodic equal 1, triang is a list of indexes's array which |
---|
[150] | 270 | ; have a surplus column. |
---|
[226] | 271 | ; We have to put it back to the initial matrix by putting indexes of |
---|
[150] | 272 | ; the last column equal to these of the last column... |
---|
[2] | 273 | ;------------------------------------------------------------ |
---|
[150] | 274 | tempdeux = systime(1) ; For key_performance =2 |
---|
| 275 | if keyword_set(key_periodic)*(nx-1 EQ jpi) $ |
---|
[226] | 276 | AND NOT keyword_set(basic) then BEGIN |
---|
[150] | 277 | indicey = triang/nx |
---|
| 278 | indicex = triang-indicey*nx |
---|
| 279 | nx = nx-1 |
---|
| 280 | liste = where(indicex EQ nx) |
---|
| 281 | if liste[0] NE -1 then indicex[liste] = 0 |
---|
| 282 | triang = indicex+nx*indicey |
---|
| 283 | nx = nx+1 |
---|
| 284 | if coinmont[0] NE -1 then begin |
---|
| 285 | indicey = coinmont/nx |
---|
| 286 | indicex = coinmont-indicey*nx |
---|
| 287 | nx = nx-1 |
---|
| 288 | liste = where(indicex EQ nx) |
---|
| 289 | if liste[0] NE -1 THEN indicex[liste] = 0 |
---|
| 290 | coinmont = indicex+nx*indicey |
---|
| 291 | nx = nx+1 |
---|
| 292 | endif |
---|
| 293 | if coindesc[0] NE -1 then begin |
---|
| 294 | indicey = coindesc/nx |
---|
| 295 | indicex = coindesc-indicey*nx |
---|
| 296 | nx = nx-1 |
---|
| 297 | liste = where(indicex EQ nx) |
---|
| 298 | if liste[0] NE -1 THEN indicex[liste] = 0 |
---|
| 299 | coindesc = indicex+nx*indicey |
---|
| 300 | nx = nx+1 |
---|
| 301 | endif |
---|
| 302 | endif |
---|
| 303 | IF testvar(var = key_performance) EQ 2 THEN $ |
---|
| 304 | print, 'temps triangule: finitions', systime(1)-tempdeux |
---|
[2] | 305 | |
---|
| 306 | ;------------------------------------------------------------ |
---|
[150] | 307 | if keyword_set(coinmonte) THEN coinmonte = coinmont ELSE twin_corners_up = coinmont |
---|
| 308 | if keyword_set(coindescend) THEN coindescend = coindesc ELSE twin_corners_dn = coindesc |
---|
[2] | 309 | ;------------------------------------------------------------ |
---|
[150] | 310 | IF NOT keyword_set(key_forgetold) THEN BEGIN |
---|
[29] | 311 | @updateold |
---|
[226] | 312 | ENDIF |
---|
[2] | 313 | |
---|
[226] | 314 | IF keyword_set(key_performance) THEN print, 'temps triangule', systime(1)-tempsun |
---|
[2] | 315 | |
---|
[150] | 316 | return, triang |
---|
[2] | 317 | |
---|
[226] | 318 | END |
---|