source: trunk/SRC/ToBeReviewed/TRIANGULATION/triangule_c.pro @ 325

;+
;
; @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.pro
; different.pro
; definetri.pro
;
; @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