source: trunk/TRIANGULATION/ciseauxtri.pro @ 2

Last change on this file since 2 was 2, checked in by opalod, 22 years ago

Initial revision
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 8.2 KB
Line 
1;------------------------------------------------------------
2;------------------------------------------------------------
3;------------------------------------------------------------
4;+
5; NAME:
6;
7; PURPOSE:vire les tableaux qui ne doivent pas etre dessines grace a 2
8; tests:
9;     1) les coins du tableau doivent etre ds la fenetre
10;     2) les clongeurs des cotes des triangfles exprimes en
11;     coordonnees normalisesne doivent pas depasser une certaine
12;     longueur seuil.
13;
14; CATEGORY:
15;
16; CALLING SEQUENCE:
17;
18; INPUTS:
19;
20; KEYWORD PARAMETERS:
21;
22; OUTPUTS:
23;
24; COMMON BLOCKS:
25;       common.pro
26;
27; SIDE EFFECTS:
28;
29; RESTRICTIONS:
30;
31; EXAMPLE:
32;
33; MODIFICATION HISTORY: Sebastien Masson (smasson@lodyc.jussieu.fr)
34;                       20/2/99
35;-
36;------------------------------------------------------------
37;------------------------------------------------------------
38;------------------------------------------------------------
39FUNCTION ciseauxtri, triang, glam, gphi, TOUT = tout, _EXTRA = ex
40@common
41   tempsun = systime(1)         ; pour key_performance
42   taille = size(glam)
43   nx = taille[1]
44   ny = taille[2]
45
46   tempdeux = systime(1)        ; pour key_performance =2
47   z = convert_coord(glam(*),gphi(*),/data,/to_normal)
48   x = z(0, *)
49   y = z(1, *)
50   tempvar = SIZE(TEMPORARY(z)) ; delete z
51   IF testvar(var = key_performance) EQ 2 THEN $
52    print, 'temps ciseauxtri: convert_coord data to_normal', systime(1)-tempdeux
53;
54; attention, suivant la projection certains points x ou y peuvent
55; devenir NaN (cf. points deriere la terre ds une projection
56; orthographique) il faut dans ce cas enlever tous les triangles qui
57; contiennent un de ces points
58;
59   if (!map.projection LE 7 AND !map.projection NE 0) $
60    OR !map.projection EQ 14 OR !map.projection EQ 15 OR !map.projection EQ 18 then begin
61      tempdeux = systime(1)     ; pour key_performance =2
62
63      ind = where(( finite((x*y)[triang[0, *]])$ ; points NaN ds la premiere colonne
64                    *finite((x*y)[triang[1, *]]) $ ; points NaN ds la 2eme colonne
65                    *finite((x*y)[triang[2, *]]) ) EQ 1) ; points NaN ds la 3eme colonne
66
67      if ind[0] NE -1 then triang = triang[*, ind] ELSE return,  -1
68      IF testvar(var = key_performance) EQ 2 THEN $
69       print, 'temps ciseauxtri: recherche points a NAN', systime(1)-tempdeux
70   endif
71;
72
73   seuil = 5.
74;
75; maintenant on vire les triangles dont un des cotes a une taille
76; superieure a 1/seuil du domaine reserve au dessin.
77;
78
79   if keyword_set(key_periodique)  then begin
80      tempdeux = systime(1)     ; pour key_performance =2
81      ind=where( (abs(x[triang[1,*]]-x[triang[0,*]]) LE (!p.position[2]-!p.position[0])/seuil) $
82                 AND (abs(x[triang[2,*]]-x[triang[1,*]]) LE (!p.position[2]-!p.position[0])/seuil) $
83                 AND (abs(x[triang[0,*]]-x[triang[2,*]]) LE (!p.position[2]-!p.position[0])/seuil) $
84;
85      AND (abs(y[triang[0,*]]-y[triang[1,*]]) LE (!p.position[3]-!p.position[1])/seuil) $
86       AND (abs(y[triang[1,*]]-y[triang[2,*]]) LE (!p.position[3]-!p.position[1])/seuil) $
87       AND (abs(y[triang[2,*]]-y[triang[0,*]]) LE (!p.position[3]-!p.position[1])/seuil) )
88
89      IF testvar(var = key_performance) EQ 2 THEN $
90       print, 'temps ciseauxtri: trouver les triangles trop grands', systime(1)-tempdeux
91;
92      tempdeux = systime(1)     ; pour key_performance =2
93      if ind[0] NE -1 then triang = triang[*, ind] ELSE return,  -1
94      IF testvar(var = key_performance) EQ 2 THEN $
95       print, 'temps ciseauxtri: virer les triangles trop grands', systime(1)-tempdeux
96   endif
97;
98; on supprime les triangles qui sont un peu trop hors du cadre
99; valable qd /TOUT est active
100;
101;    if keyword_set(key_irregular) then begin
102
103;       tempdeux = systime(1)     ; pour key_performance =2
104;       seuil = 0
105;       seuil = .1*max([!p.position[2]-!p.position[0],!p.position[3]-!p.position[1]])
106;       ind=where(x[triang[0,*]] GE !p.position[0]-seuil AND x[triang[0,*]] LE !p.position[2]+seuil $
107;                 AND x[triang[1,*]] GE !p.position[0]-seuil AND x[triang[1,*]] LE !p.position[2]+seuil $
108;                 AND x[triang[2,*]] GE !p.position[0]-seuil AND x[triang[2,*]] LE !p.position[2]+seuil $
109; ;
110;       AND y[triang[0,*]] GE !p.position[1]-seuil AND y[triang[0,*]] LE !p.position[3]+seuil $
111;        AND y[triang[1,*]] GE !p.position[1]-seuil AND y[triang[1,*]] LE !p.position[3]+seuil $
112;        AND y[triang[2,*]] GE !p.position[1]-seuil AND y[triang[2,*]] LE !p.position[3]+seuil )
113;       IF testvar(var = key_performance) EQ 2 THEN $
114;        print, 'temps ciseauxtri: trouver les triangles hors du cadre', systime(1)-tempdeux
115; ;
116;       tempdeux = systime(1)     ; pour key_performance =2
117;       if ind[0] NE -1 then triang = triang[*, ind] ELSE return,  -1
118;       IF testvar(var = key_performance) EQ 2 THEN $
119;        print, 'temps ciseauxtri: virer les triangles hors du cadre', systime(1)-tempdeux
120;    ENDIF
121;
122;        dernier tri
123;
124   if n_elements(ind) GT 1 then BEGIN
125;
126; il faut verifier que les triangles que l''on garde ne
127; forment pas une triangulation dans laquelle 2 triangles ont un
128; sommet en commun sans avoir de cote de commun.
129;
130; on constitue la liste des rectangles que l''on souhaite garder (on
131; garde un rectangle des qu''il y a un triangle dedans)
132;
133; dans definetri, on a construit les triangles tels que le premier et
134; le dernier sommets soient ceux de la diagonale du rectangle definit
135; par le maillage.
136; pour retouver de quel rectangle provient un triangle, on cherche le
137; min de l''indice suivant x et suivant y de chaque triangle. Apres on
138; repasse cette indissage suivant x et y en indicage suivant nx*ny
139;
140      tempdeux = systime(1)     ; pour key_performance =2
141; indices en y des rectangles
142      indytriang = (triang[0, *]/nx) < (triang[2, *]/nx)
143; indices en x des rectangles
144      indxtriang0 = triang[0, *]-nx*(triang[0, *]/nx)
145      indxtriang2 = triang[2, *]-nx*(triang[2, *]/nx)
146      indxmin = indxtriang0 < indxtriang2
147; attention dans le cas ou la grille est periodique et ou on a tous
148; les points suivant x, les triangles qui assurent la periodicite en x
149; ont des indices suivant x qui sont 0 et nx-1. Ils appatienent au
150; rectangles de la colonne nx-1 et non de la colonne 0
151      if keyword_set(key_periodique) AND nx EQ jpi then BEGIN
152      indxmax = indxtriang0 > indxtriang2
153      indxtriang = indxmin + (nx-1)*(indxmin EQ 0 AND indxmax EQ (nx-1))
154      ENDIF ELSE indxtriang = indxmin
155      listrect = nx*indytriang+indxtriang
156      IF testvar(var = key_performance) EQ 2 THEN $
157       print, 'temps ciseauxtri: liste des rectangles', systime(1)-tempdeux
158;
159; maintenant qu''on a cette liste on va s''assuter que l''on a pas de
160; triangles qui n''ont qu''on sommet en commun.
161;
162      test = bytarr(nx, ny)
163      test[listrect] = 1
164      dejavire = 1b-test
165;
166      tempdeux = systime(1)     ; pour key_performance =2
167      vire1 = 0
168      vire2 = 0
169      while (vire1[0] NE -1 OR vire2[0] NE -1) ne 0 do begin
170         vire1 = where( (test*shift(test, -1, -1) $
171                         *(1-shift(test, 0, -1))*(1-shift(test, -1, 0))) EQ 1)
172         if vire1[0] NE -1 THEN test[vire1] = 0 ; on vire le rectangle
173         vire2 = where( ((1-test)*(1-shift(test, -1, -1)) $
174                         *shift(test, 0, -1)*shift(test, -1, 0)) EQ 1)
175; on vire le rectangle du dessus (meme indice x mais egale a 1)
176         if vire2[0] NE -1 THEN test[vire2+nx] = 0
177      ENDWHILE
178;stop
179      test = test+temporary(dejavire)
180;
181      avirer = where(temporary(test) EQ 0)
182      IF testvar(var = key_performance) EQ 2 THEN $
183       print, 'temps ciseauxtri: determinationdes rectangles a virer', systime(1)-tempdeux
184;
185      if avirer[0] NE -1 then begin
186      tempdeux = systime(1)     ; pour key_performance =2
187      indnx = n_elements(listrect)
188      indny = n_elements(avirer)
189      ind = listrect[*]#replicate(1l, indny)
190      ind = ind EQ replicate(1, indnx)#avirer[*]
191      if indny GT 1 then ind = total(ind, 2)
192      ind = where(ind EQ 0)
193      if ind[0] NE -1 then triang = triang[*, ind] ELSE return,  -1
194      endif
195      IF testvar(var = key_performance) EQ 2 THEN $
196       print, 'temps ciseauxtri: derniere retouche de la triangulation', systime(1)-tempdeux
197   endif
198;   
199;
200   if keyword_set(key_performance) THEN print, 'temps ciseauxtri', systime(1)-tempsun
201   return,  triang
202end
Note: See TracBrowser for help on using the repository browser.