;+
;
; @file_comments
; Construct the triangulation array.
;
; The idea is: construct a list of triangle which link points between them.
; This is automatically done by the function TRIANGULATE
; Here:
; we consider the fact that points are disposed on a grid (regular or not,
; but not unstructured, that is to say that points are written following a
; rectangular matrix). A easy way to do triangles between all points is then:
;
; for each point (i,j) of the matrix -except those of the last line and of
; the last column- we call rectangle (i,j) the rectangle made of the four
; points (i,j), (i+1,j), (i,j+1), (i+1,j+1). To trace all triangle, we just
; have to trace the 2 triangles contained in rectangles (i,j)
;
; We notice that each rectangle (i,j) have 2 diagonals (it is true... Make a
; drawing to make sure!!), so there are two possible choice for each rectangle
; we want to cut in 2 triangles...
;
; It is thanks to this choice that we will be able to trace coast with right
; angles. At each angle of coast remarkable by the existence of an unique land
; point or of an unique ocean point on one of the four summit of a rectangle (i,j),
; we have to cut the rectangle following the diagonal passing by this point.
;
; @categories
; Graphics
;
; @param MASKENTREE {in}{optional}{type=2d array}
; It is a 2d array which will serve to mask the field we will trace after with CONTOUR,
; ...TRIANGULATION=triangule(mask)
; If this argument is not specified, the function use tmask
;
; @keyword BASIC
; Specify that the mask is on a basic grid (use the triangulation for vertical cuts and hovmoellers)
;
; @keyword KEEP_CONT
; To keep the triangulation even on the continents
;
; @keyword COINMONTE {type=array}
; To obtain the array of "ascending land corner" to be treated with
; completecointerre.pro in the variable array instead of make it pass by the global
; variable twin_corners_up.
;
; @keyword COINDESCEND {type=array}
; See COINMONTE
;
; @returns
; res: tableau 2d (3,nbre de triangles).
; Each line of res represent indexes of points constituting summits of a triangle.
; See how we trace triangles in definetri.pro
;
; @uses
; common
; different
; definetri
;
; @restrictions
; Data whose we want to do the contour must be disposed in a matrix.
; On the other hand, in the matrix, the points's arrangement can not be
; irregular. If it is, use TRIANGULE.
;
; @history
; Sebastien Masson (smasson\@lodyc.jussieu.fr)
; 26/4/1999
;
; @version
; $Id$
;
; @todo
; seb L.267->268 je ne pense pas que ce soit ce que tu voulais dire mais
; c'est la traduction de ce qu'il y avait écrit. Correction si besoin.
;
;-
FUNCTION triangule_c, maskentree $
, COINMONTE=coinmonte, COINDESCEND=coindescend $
, BASIC=basic, KEEP_CONT=keep_cont
;
compile_opt idl2, strictarrsubs
;
tempsun = systime(1) ; For key_performance
;---------------------------------------------------------
@cm_4mesh
IF NOT keyword_set(key_forgetold) THEN BEGIN
@updatenew
ENDIF
;------------------------------------------------------------
; Is the mask given or do we have to take tmask?
;------------------------------------------------------------
;
msk = maskentree
taille = size(msk)
nx = taille[1]
ny = taille[2]
;
IF n_elements(keep_cont) EQ 0 THEN keep_cont = 1-key_irregular
;------------------------------------------------------------
if keyword_set(key_periodic)*(nx EQ jpi) $
AND NOT keyword_set(basic) then BEGIN
msk = [msk, msk[0, *]]
nx = nx+1
ENDIF
;------------------------------------------------------------
; We will find the list of rectangles (i,j)(located by their left
; bottom corner) we have to cut following a descendant diagonal.
; We will call this list : pts_downward
;
pts_downward = 0
; We construct the test which allow to find this triangle :
;
;
; shift(msk, 0, -1)------------shift(msk, -1, -1)
; | |
; | |
; | |
; | |
; msk---------------------shift(msk, -1, 0)
;
sum1 = msk+shift(msk, -1, 0)+shift(msk, -1, -1) ;points which surround the left top point.
sum2 = msk+shift(msk, 0, -1)+shift(msk, -1, -1) ;points which surround the right bottom point.
tempdeux = systime(1) ; For key_performance =2
; The left top land point surrounded by ocean points
liste = where( (4-sum1)*(1-shift(msk, 0, -1)) EQ 1 )
if liste[0] NE -1 THEN pts_downward = [pts_downward,liste ]
; The left top ocean point surrounded by land points
liste = where( (1-sum1)*shift(msk, 0, -1) EQ 1)
if liste[0] NE -1 THEN pts_downward = [pts_downward,liste ]
; The right bottom land point surrounded by ocean points
liste = where( (4-sum2)*(1-shift(msk, -1, 0)) EQ 1)
if liste[0] NE -1 THEN pts_downward = [pts_downward,liste ]
; The right bottom ocean point surrounded by land points
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(basic)) OR keyword_set(coinmonte) OR keyword_set(coindescend) then begin
tempdeux = systime(1) ; For key_performance =2
;2 land points in ascendant diagonal with 2 ocean points in descendant diagonal.
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
;
coindesc = where( ((1-shift(msk, 0, -1))*(1-shift(msk, -1, 0)) $
*msk*shift(msk, -1, -1) EQ 1) )
;
;2 land points in descendant diagonal with 2 ocean points in ascendant diagonal.
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) ; For 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) ; For key_performance =2
pts_downward = pts_downward[1:n_elements(pts_downward)-1]
pts_downward = pts_downward[uniq(pts_downward, sort(pts_downward))]
; None rectangle can have an element of the last column or of the
; last line as left bottom corner.
; so we have to remove these points if they has been selected in 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(basic)) OR keyword_set(coinmonte) OR keyword_set(coindescend) 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) ; For 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
;------------------------------------------------------------
; We delete land points which only contain land points.
;------------------------------------------------------------
;
;
if (NOT keyword_set(basic)) AND (NOT keyword_set(keep_cont)) then begin
tempdeux = systime(1) ; For key_performance =2
; We delete rectangles which are entirely in the land.
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
; We do an other sort :
; We have to do not remove rectangles which only have one common summit.
; 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) ; For key_performance =2
indice[recdsterre] = 1
if NOT keyword_set(basic) then begin
vire1 = 0
vire2 = 0
while (vire1[0] NE -1 OR vire2[0] NE -1) ne 0 do begin
; Delete rectangles we have to remove from recsterre (in fact those we have
; to keep although they are entirely in the land).
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 ; The last column and the last line
indice[nx-1, *] = 1 ; can not define any rectangle.
;
tempdeux = systime(1) ; For key_performance =2
recgarde = where(indice EQ 0)
; We recuperate numbers of triangles we will keep.
trigarde = 2*[recgarde-recgarde/nx]
trigarde = transpose(temporary(trigarde))
trigarde = [trigarde, trigarde+1]
;
triang = triang[*, temporary(trigarde[*])]
IF testvar(var = key_performance) EQ 2 THEN $
print, 'temps triangule: virer les triangle de la liste', systime(1)-tempdeux
endif
endif
; print, 'temps tri triangles', systime(1)-t1
;------------------------------------------------------------
; When key_periodic equal 1, triang is a list of indexes's array which
; have a surplus column.
; We have to put it back to the initial matrix by putting indexes of
; the last column equal to these of the last column...
;------------------------------------------------------------
tempdeux = systime(1) ; For key_performance =2
if keyword_set(key_periodic)*(nx-1 EQ jpi) $
AND NOT keyword_set(basic) 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 keyword_set(coinmonte) THEN coinmonte = coinmont ELSE twin_corners_up = coinmont
if keyword_set(coindescend) THEN coindescend = coindesc ELSE twin_corners_dn = coindesc
;------------------------------------------------------------
IF NOT keyword_set(key_forgetold) THEN BEGIN
@updateold
ENDIF
IF keyword_set(key_performance) THEN print, 'temps triangule', systime(1)-tempsun
return, triang
END