Tracer des surfaces de niveau ?

Sujets variés concernant le développement en PureBasic
Répondre
Dräc
Messages : 526
Inscription : dim. 29/août/2004 0:45

Tracer des surfaces de niveau ?

Message par Dräc »

Je cherche à tracer les surfaces de niveau d’un plan discret (données contenues dans un tableau de valeur).
Pour ne pas réinventer la roue, quelqu’un connaît-il un algorithme rapide (fonctions 3D excluses)?

Exemple avec un tableau 6x6:
111111
122221
123321
123321
122221
111111

on devrait voir apparaître des disques concentriques.
Haut
erix14
Messages : 480
Inscription : sam. 27/mars/2004 16:44
Contact :
Contacter erix14

Message par erix14 »

J'avais fait ça il y a un an, si ça peut aider... :D

;/ Programme Erix14
 ;/ PureBasic 3.94
 #Longueur = 200
#Largeur = 200
#Window = 0
Global hWnd
Structure UnPoint
     Epaisseur.f
     NormalX.f
     NormalY.f
     NormalZ.f
EndStructure
Dim Buffet3D. UnPoint ( #Longueur , #Largeur )
Procedure RectangleArrondi3D(PositionX,PositionY,longueur,largeur,rayon,Epaisseur)
     Zone = rayon-Epaisseur
      For j = 0 To rayon : JJ = j + PositionY; coin supérieur gauche
           For i = 0 To rayon : II = i + PositionX
               x = i - rayon : y = j - rayon : Distance.f = Sqr (x*x + y*y)
                If Distance <= rayon And Distance >= Zone
                    D.l = Distance-Zone
                    Buffet3D(II,JJ)\NormalX = (x*D)/(Epaisseur*Distance)
                    Buffet3D(II,JJ)\NormalY = (y*D)/(Epaisseur*Distance)
                    Buffet3D(II,JJ)\NormalZ = Sqr (Epaisseur*Epaisseur-D*D)/Epaisseur
                    Buffet3D(II,JJ)\Epaisseur = Buffet3D(II,JJ)\NormalZ*Epaisseur
                ElseIf Distance < Zone
                    Buffet3D(II,JJ)\Epaisseur = Epaisseur
                    Buffet3D(II,JJ)\NormalX = 0
                    Buffet3D(II,JJ)\NormalY = 0
                    Buffet3D(II,JJ)\NormalZ = 1
                Endif
           Next
      Next
      For j = rayon To rayon*2 : JJ = PositionY + j + largeur - 2*rayon; coin inférieur droit
           For i = rayon To rayon*2 : II = PositionX + i + longueur - 2*rayon
               x = i - rayon : y = j - rayon : Distance.f = Sqr (x*x + y*y)
                If Distance <= rayon And Distance >= Zone
                    D.l = Distance-Zone
                    Buffet3D(II,JJ)\NormalX = (x*D)/(Epaisseur*Distance)
                    Buffet3D(II,JJ)\NormalY = (y*D)/(Epaisseur*Distance)
                    Buffet3D(II,JJ)\NormalZ = Sqr (Epaisseur*Epaisseur-D*D)/Epaisseur
                    Buffet3D(II,JJ)\Epaisseur = Buffet3D(II,JJ)\NormalZ*Epaisseur
                ElseIf Distance < Zone
                    Buffet3D(II,JJ)\Epaisseur = Epaisseur
                    Buffet3D(II,JJ)\NormalX = 0
                    Buffet3D(II,JJ)\NormalY = 0
                    Buffet3D(II,JJ)\NormalZ = 1
                Endif
           Next
      Next
      For j = 0 To rayon : JJ = PositionY + j; coin supérieur droit
           For i = rayon To rayon*2 : II = PositionX + i + longueur - 2*rayon
               x = i - rayon : y = j - rayon : Distance.f = Sqr (x*x + y*y)
                If Distance <= rayon And Distance >= Zone
                    D.l = Distance-Zone
                    Buffet3D(II,JJ)\NormalX = (x*D)/(Epaisseur*Distance)
                    Buffet3D(II,JJ)\NormalY = (y*D)/(Epaisseur*Distance)
                    Buffet3D(II,JJ)\NormalZ = Sqr (Epaisseur*Epaisseur-D*D)/Epaisseur
                    Buffet3D(II,JJ)\Epaisseur = Buffet3D(II,JJ)\NormalZ*Epaisseur
                ElseIf Distance < Zone
                    Buffet3D(II,JJ)\Epaisseur = Epaisseur
                    Buffet3D(II,JJ)\NormalX = 0
                    Buffet3D(II,JJ)\NormalY = 0
                    Buffet3D(II,JJ)\NormalZ = 1
                Endif
           Next
      Next
      For j = rayon To rayon*2 : JJ = PositionY + j + largeur - 2*rayon; coin inférieur gauche
           For i = 0 To rayon : II = PositionX + i
               x = i - rayon : y = j - rayon : Distance.f = Sqr (x*x + y*y)
                If Distance <= rayon And Distance >= Zone
                    D.l = Distance-Zone
                    Buffet3D(II,JJ)\NormalX = (x*D)/(Epaisseur*Distance)
                    Buffet3D(II,JJ)\NormalY = (y*D)/(Epaisseur*Distance)
                    Buffet3D(II,JJ)\NormalZ = Sqr (Epaisseur*Epaisseur-D*D)/Epaisseur
                    Buffet3D(II,JJ)\Epaisseur = Buffet3D(II,JJ)\NormalZ*Epaisseur
                ElseIf Distance < Zone
                    Buffet3D(II,JJ)\Epaisseur = Epaisseur
                    Buffet3D(II,JJ)\NormalX = 0
                    Buffet3D(II,JJ)\NormalY = 0
                    Buffet3D(II,JJ)\NormalZ = 1
                Endif
           Next
      Next
     RR = PositionY + rayon
      For y=rayon To largeur-rayon
          JJ = PositionY + y
           For x=0 To Epaisseur
               II = PositionX + x
               Buffet3D(II,JJ)\Epaisseur = Buffet3D(II,RR)\Epaisseur; côté gauche
               Buffet3D(II,JJ)\NormalX = Buffet3D(II,RR)\NormalX
               Buffet3D(II,JJ)\NormalY = Buffet3D(II,RR)\NormalY
               Buffet3D(II,JJ)\NormalZ = Buffet3D(II,RR)\NormalZ
               x1 = II + longueur - Epaisseur
               Buffet3D(x1,JJ)\Epaisseur = Buffet3D(x1,RR)\Epaisseur; côté droit
               Buffet3D(x1,JJ)\NormalX = Buffet3D(x1,RR)\NormalX
               Buffet3D(x1,JJ)\NormalY = Buffet3D(x1,RR)\NormalY
               Buffet3D(x1,JJ)\NormalZ = Buffet3D(x1,RR)\NormalZ
           Next
           For x=Epaisseur To longueur-Epaisseur; centre
               II = PositionX + x
               Buffet3D(II,JJ)\Epaisseur = Epaisseur
               Buffet3D(II,JJ)\NormalX = 0
               Buffet3D(II,JJ)\NormalY = 0
               Buffet3D(II,JJ)\NormalZ = 1
           Next
      Next
     RR = PositionX + rayon
      For x=rayon To longueur-rayon
          II = PositionX + x
           For y=0 To Epaisseur
               JJ = PositionY + y
               Buffet3D(II,JJ)\Epaisseur = Buffet3D(RR,JJ)\Epaisseur; côté supérieur
               Buffet3D(II,JJ)\NormalX = Buffet3D(RR,JJ)\NormalX
               Buffet3D(II,JJ)\NormalY = Buffet3D(RR,JJ)\NormalY
               Buffet3D(II,JJ)\NormalZ = Buffet3D(RR,JJ)\NormalZ
               y1 = JJ + largeur - Epaisseur
               Buffet3D(II,y1)\Epaisseur = Buffet3D(RR,y1)\Epaisseur; côté inférieur
               Buffet3D(II,y1)\NormalX = Buffet3D(RR,y1)\NormalX
               Buffet3D(II,y1)\NormalY = Buffet3D(RR,y1)\NormalY
               Buffet3D(II,y1)\NormalZ = Buffet3D(RR,y1)\NormalZ
           Next
      Next
      For x=rayon To longueur-rayon
          II = PositionX + x
           For y=Epaisseur To rayon
               JJ = PositionY + y
               Buffet3D(II,JJ)\Epaisseur = Epaisseur; centre supérieur
               Buffet3D(II,JJ)\NormalX = 0
               Buffet3D(II,JJ)\NormalY = 0
               Buffet3D(II,JJ)\NormalZ = 1
           Next
           For y=largeur-rayon To largeur-Epaisseur
               JJ = PositionY + y
               Buffet3D(II,JJ)\Epaisseur = Epaisseur; centre inférieur
               Buffet3D(II,JJ)\NormalX = 0
               Buffet3D(II,JJ)\NormalY = 0
               Buffet3D(II,JJ)\NormalZ = 1
           Next
      Next
EndProcedure
Procedure Rendu3D(longueur,largeur,couleur)
     LR = Red(couleur) : LG = Green(couleur) : LB = Blue (couleur); Couleur du rendu
     LX.f = 30 : LY.f = 30 : LZ.f = 100; Position de la lampe ponctuel
     PR = 128 : PG = 128 : PB = 128; Lumière de la lampe ponctuel
     AR = 32 : AG = 32 : AB = 32; Lumière d'ambiance
      For y=0 To largeur
           For x=0 To longueur
                If Buffet3D(x,y)\Epaisseur = 0 : Continue : Endif
               Distance.f = Sqr( Pow (x-LX,2)+ Pow (y-LY,2)+ Pow (Buffet3D(x,y)\Epaisseur-LZ,2))
               DirX.f = (x-LX)/Distance
               DirY.f = (y-LY)/Distance
               DirZ.f = (Buffet3D(x,y)\Epaisseur-LZ)/Distance
               K.f = -(Buffet3D(x,y)\NormalX*DirX + Buffet3D(x,y)\NormalY*DirY + Buffet3D(x,y)\NormalZ*DirZ)
               r = LR + AR + K*PR : If r > 255 : r = 255 : Endif : If r < 0 : r = 0 : Endif
               g = LG + AG + K*PG : If g > 255 : g = 255 : Endif : If g < 0 : g = 0 : Endif
               b = LB + AB + K*PB : If b > 255 : b = 255 : Endif : If b < 0 : b = 0 : Endif
                Plot (x,y,RGB(r,g,b))
           Next
      Next
EndProcedure
;- Debut du programme
hWnd = OpenWindow ( #Window , 0, 0, #Longueur , #Largeur , #PB_Window_SystemMenu | #PB_Window_Invisible | #PB_Window_ScreenCentered , " RectangleArrondi3D ")
CreateGadgetList (WindowID())
;{/ Image
 CreateImage (0, #Longueur , #Largeur )
StartDrawing (ImageOutput())
RectangleArrondi3D(0,0, #Longueur-1 , #Largeur-1 ,100,100)
Rendu3D( #Longueur-1 , #Largeur-1 , RGB (128,128,0))
StopDrawing ();}
ImageGadget (0, 0, 0, 0, 0, ImageID())
HideWindow ( #Window ,0)
Repeat
      Select WaitWindowEvent ()
           Case #PB_Event_CloseWindow
               quit = #True
      EndSelect
Until quit
Haut
Avatar de l’utilisateur
djes
Messages : 4252
Inscription : ven. 11/févr./2005 17:34
Localisation : Arras, France

Message par djes »

Ouahow, y'a des choses bien là!

Erix14> Est-ce que ton code peut servir à tracer des courbes d'isovaleurs?
Haut
Dräc
Messages : 526
Inscription : dim. 29/août/2004 0:45

Message par Dräc »

Oui, je me souviens :)
Cependant il est très spécifique. Mais je le garde en tete, surtout pour le traitement des dégradés.

Ma demande porte sur le tracer des isocontours que l’on colorise par la suite.
J’aurai souhaité savoir si quelqu’un avec des tuyaux sur le sujet ?

De mon coté, j’ai trouvé cet algorithme sur le net:
http://astronomy.swin.edu.au/~pbourke/p ... on/conrec/
Haut
erix14
Messages : 480
Inscription : sam. 27/mars/2004 16:44
Contact :
Contacter erix14

Message par erix14 »

djes> Rendu3D(...) transformer un tableau de normals et de hauteurs de points en une image, c'est à l'utilisateur d'alimenter ce tableau (qui s'appelle Buffet3D) Dans mon programme c'est la Procedure RectangleArrondi3D qui alimentent ce tableau. On peut très bien imaginer une procédure pour tracer des courbes, mais Rendu3D ne gère pas les ombres... :cry:
C'est une procédure à améliorer... :D
Haut
Avatar de l’utilisateur
djes
Messages : 4252
Inscription : ven. 11/févr./2005 17:34
Localisation : Arras, France

Message par djes »

Non mais c'est génial! Ca fait des lustres que je cherche un bon algo pour ça! Vous savez que ça se vend cher ce genre de trucs?! :P

Merci!
Haut
Dräc
Messages : 526
Inscription : dim. 29/août/2004 0:45

Message par Dräc »

djes a écrit :Non mais c'est génial! Ca fait des lustres que je cherche un bon algo pour ça! Vous savez que ça se vend cher ce genre de trucs?! :P

Merci!
Tu en veux encore du gratuit : Voici toute une librairie de fonctions de tracé
http://www.astro.caltech.edu/~tjp/pgplot/

avec la présentation des fonctions (dont celles d’isocontour)
http://www.sns.ias.edu/Main/linux/pgplo ... tines.html

Il faut surtout remercier ceux qui partagent ;)
Haut
Avatar de l’utilisateur
djes
Messages : 4252
Inscription : ven. 11/févr./2005 17:34
Localisation : Arras, France

Message par djes »

Tu as raison, aussi je vais faire un petit don moi aussi :)
Haut
Dräc
Messages : 526
Inscription : dim. 29/août/2004 0:45

Message par Dräc »

Voici deux traductions en PureBasic du code CONREC (http://astronomy.swin.edu.au/~pbourke/p ... on/conrec/).

Le code commence part renseigner les tableaux de données, i.e. d(i,j), x(i), y(j) et z(k) où
d(i,j) = f( x(i), y(j) ), f la fonction étudiée (ici f(x,y) = 1/( (x^2+(y-0.842)*(y+0.842))^2 + ( x*(y-0.842) + x*(y-0.842))^2, exemple 1 de l’article)

z(k) contient les isocontours à afficher


Je qualifierais de littérale la première version du code. C’est à dire que c’est la plus proche du code original. Cette version permet donc de faire facilement le lien entre l’article et le code.

La deuxième version cherche à rendre indépendante la routine CONREC vis à vis des noms de tableau (et éviter ainsi de duplique les données en mémoire).
Elle se débrouille donc à lire directement les données grâce aux adresses des tableaux.

Pour cela, je propose une méthode valable pour les tableaux de 1 et 2 dimensions (extensible à n dimension) qui n’est rien d’autre que ce que le compilo écrit lui-même.

Cette méthode utilise une Structure « array » et une procédure « Array( d.array, i, j) »

Il est claire que cette dernière version cherche à compenser humblement un manque du langage :)

Version n°1

Code : Tout sélectionner

#Longueur = 200 
#Largeur = 200 
#Window = 0

Procedure.f min(x.f, y.f)
  If x < y
    ProcedureReturn x
  Else
    ProcedureReturn y
  EndIf
EndProcedure

Procedure.f max(x.f, y.f)
  If x > y
    ProcedureReturn x
  Else
    ProcedureReturn y
  EndIf
EndProcedure

Procedure.f float(x.l)
  ProcedureReturn x
EndProcedure

xmin.f = -1.5
xmax.f = 1.5
nx = 100

ymin.f = -1.5
ymax.f = 1.5
ny = 100

nc = 20

Dim d.f(nx, ny)
Dim x.f(nx)
Dim y.f(ny)
Dim z.f(nc)
Dim colorz.f(nc)

For i = 0 To nx
 x(i) = xmin + (xmax - xmin) * i / float(nx)
Next
For j = 0 To ny
 y(j) = ymax - (ymax - ymin) * j / float(ny)
Next

zmax.f = -float(999999999)
zmin.f = float(999999999)
For i = 0 To nx
  For j = 0 To ny
    d(i,j)= 1./( Pow(x(i)*x(i)+(y(j)-0.842)*(y(j)+0.842),2) + Pow( x(i)*(y(j)-0.842) + x(i)*(y(j)-0.842),2))
    zmin = min(zmin, d(i,j))
    zmax = max(zmax, d(i,j))
  Next
Next

zmax = 2.25
zmin = 0

For k = 0 To nc
 z(k) = (zmax - zmin) * k / float(nc+1)
 colorz(k) = Int(255 * k / float(nc+1))
Next

Global xmax, xmin, ymax, ymin, zmax, zmin

Dim h.f(4)
Dim sh(4)
Dim xh.f(4)
Dim yh.f(4)

Dim im(3)
Dim jm(3)
Dim castab(2,2,2)

; Initialisation des tableaux
Restore Donnees
For i = 0 To 3
  Read im(i)
Next
For i = 0 To 3
  Read jm(i)
Next

For i = 0 To 2
  For j = 0 To 2
    For k = 0 To 2
       Read castab(i, j, k)
    Next
  Next
Next

Procedure.f xsect(p1, p2)
  ProcedureReturn (h(p2)*xh(p1)-h(p1)*xh(p2))/(h(p2)-h(p1))
EndProcedure
Procedure.f ysect(p1, p2)
  ProcedureReturn (h(p2)*yh(p1)-h(p1)*yh(p2))/(h(p2)-h(p1))
EndProcedure

Procedure vecout(x1.f, y1.f, x2.f, y2.f, k)
  x1 = #Longueur * (x1-xmin)/(xmax-xmin)
  y1 = #Largeur * (y1-ymin)/(ymax-ymin)
  x2 = #Longueur * (x2-xmin)/(xmax-xmin)
  y2 = #Largeur * (y2-ymin)/(ymax-ymin)
  
  LineXY(x1, y1, x2, y2, colorz(k))
  ;Debug Str(x1)+"/"+Str(y1)+"/"+Str(x2)+"/"+Str(y2)+"/"+Str(z)
EndProcedure

;=============================================================================
;
;     CONREC is a contouring subroutine for rectangularily spaced data.
;
;     It emits calls to a line drawing subroutine supplied by the user
;     which draws a contour map corresponding to real*4data on a randomly
;     spaced rectangular grid. The coordinates emitted are in the same
;     units given in the x() and y() arrays.
;
;     Any number of contour levels may be specified but they must be
;     in order of increasing value.
;
;     This version is an adaptation of the original Paul Bourke subroutine:
;     http://astronomy.swin.edu.au/~pbourke/projection/conrec/ 
;
;=============================================================================
Procedure conrec(ilb, iub, jlb, jub, nc)
  For j= jub-1 To jlb Step -1
    For i=ilb To iub-1
      temp1.f = min(d(i, j), d(i, j+1))
      temp2.f = min(d(i+1, j), d(i+1, j+1))
      dmin.f = min(temp1,temp2)
      temp1 = max(d(i,j),d(i,j+1))
      temp2 = max(d(i+1,j),d(i+1,j+1))
      dmax.f = max(temp1,temp2)

      If dmax>=z(0) And dmin<=z(nc) 
        For k=0 To nc
        
	        If z(k)>=dmin And z(k)<=dmax
            For m=4 To 0 Step -1
	           If m
            	h(m) = d(i+im(m-1), j+jm(m-1))-z(k)
            	xh(m) = x(i+im(m-1))
              yh(m) = y(j+jm(m-1))
	           Else
           		h(0) = 0.25*(h(1)+h(2)+h(3)+h(4))
	            xh(0)=0.5*(x(i)+x(i+1))
              yh(0)=0.5*(y(j)+y(j+1)) 
             EndIf
             
	           If h(m)>0
               sh(m) = 1
             ElseIf h(m) < 0
               sh(m) = -1
             Else
               sh(m) = 0
             EndIf
           Next
	    ;=================================================================
	    ;
	    ; Note: at this stage the relative heights of the corners and the
	    ; centre are in the h array, and the corresponding coordinates are
	    ; in the xh and yh arrays. The centre of the box is indexed by 0
	    ; and the 4 corners by 1 to 4 as shown below.
	    ; Each triangle is then indexed by the parameter m, and the 3
	    ; vertices of each triangle are indexed by parameters m1,m2,and
	    ; m3.
	    ; It is assumed that the centre of the box is always vertex 2
	    ; though this isimportant only when all 3 vertices lie exactly on
	    ; the same contour level, in which case only the side of the box
	    ; is drawn.
	    ;
	    ;
	    ;      vertex 4 +-------------------+ vertex 3
	    ;               | \               / |
	    ;               |   \    m-3    /   |
	    ;               |     \       /     |
	    ;               |       \   /       |
	    ;               |  m=2    X   m=2   |       the centre is vertex 0
	    ;               |       /   \       |
	    ;               |     /       \     |
	    ;               |   /    m=1    \   |
	    ;               | /               \ |
	    ;      vertex 1 +-------------------+ vertex 2
	    ;
	    ;
	    ;
	    ;               Scan each triangle in the box
	    ;
	    ;=================================================================             
	    For m=1 To 4
	      m1 = m
	      m2 = 0
	      If m<>4
          m3 = m+1
	      Else
          m3 = 1
        EndIf
	      
        case_value = castab(sh(m1)+1, sh(m2)+1, sh(m3)+1);
 	      If case_value<>0
    Select case_value
		  ;===========================================================
		  ;     Case 1 - Line between vertices 1 and 2
		  ;===========================================================
		Case 1
		  x1.f=xh(m1);
		  y1.f=yh(m1);
		  x2.f=xh(m2);
		  y2.f=yh(m2);
		  ;===========================================================
		  ;     Case 2 - Line between vertices 2 and 3
		  ;===========================================================
		Case 2
		  x1.f=xh(m2);
		  y1.f=yh(m2);
		  x2.f=xh(m3);
		  y2.f=yh(m3);
		  ;===========================================================
		  ;     Case 3 - Line between vertices 3 and 1
		  ;===========================================================
		Case 3
		  x1.f=xh(m3);
		  y1.f=yh(m3);
		  x2.f=xh(m1);
		  y2.f=yh(m1);
		  ;===========================================================
		  ;     Case 4 - Line between vertex 1 and side 2-3
		  ;===========================================================
		Case 4
		  x1.f=xh(m1);
		  y1.f=yh(m1);
		  x2.f=xsect(m2,m3);
		  y2.f=ysect(m2,m3);
		  ;===========================================================
		  ;     Case 5 - Line between vertex 2 and side 3-1
		  ;===========================================================
		Case 5
		  x1.f=xh(m2);
		  y1.f=yh(m2);
		  x2.f=xsect(m3,m1);
		  y2.f=ysect(m3,m1);
		  ;===========================================================
		  ;     Case 6 - Line between vertex 3 and side 1-2
		  ;===========================================================
		Case 6
		  x1.f=xh(m3);
		  y1.f=yh(m3);
		  x2.f=xsect(m1,m2);
		  y2.f=ysect(m1,m2);
		  ;===========================================================
		  ;     Case 7 - Line between sides 1-2 and 2-3
		  ;===========================================================
		Case 7
		  x1.f=xsect(m1,m2);
		  y1.f=ysect(m1,m2);
		  x2.f=xsect(m2,m3);
		  y2.f=ysect(m2,m3);
		  ;===========================================================
		  ;     Case 8 - Line between sides 2-3 and 3-1
		  ;===========================================================
		Case 8
		  x1.f=xsect(m2,m3);
		  y1.f=ysect(m2,m3);
		  x2.f=xsect(m3,m1);
		  y2.f=ysect(m3,m1);
		  ;===========================================================
		  ;     Case 9 - Line between sides 3-1 and 1-2
		  ;===========================================================
		Case 9
		  x1.f=xsect(m3,m1);
		  y1.f=ysect(m3,m1);
		  x2.f=xsect(m1,m2);
		  y2.f=ysect(m1,m2);
		EndSelect       
		;=============================================================
		; Put your graphic processing code here 
		;=============================================================
    vecout(x1,y1,x2,y2, k)
          EndIf
        
            Next              
          EndIf
        Next
      EndIf
      
    Next
  Next
    
EndProcedure

DataSection
Donnees:
Data.l 0, 1, 1, 0; im()
Data.l 0, 0, 1, 1; jm()
Data.l 0,0,8, 0,2,5, 7,6,9; case_value
Data.l 0,3,4, 1,3,1, 4,3,0
Data.l 9,6,7, 5,2,0, 8,0,0
EndDataSection

;- Debut du programme 

hWnd = OpenWindow ( #Window , 0, 0, #Longueur , #Largeur , #PB_Window_SystemMenu | #PB_Window_Invisible | #PB_Window_ScreenCentered , " Isocontour") 
CreateGadgetList (WindowID()) 
CreateImage (0, #Longueur , #Largeur )
 
StartDrawing (ImageOutput()) 
conrec(0, nx, 0, ny, nc)
StopDrawing ()

ImageGadget (0, 0, 0, 0, 0, ImageID()) 
HideWindow ( #Window ,0) 
Repeat 
      Select WaitWindowEvent () 
           Case #PB_Event_CloseWindow 
               Quit = #True 
      EndSelect 
Until Quit
Version n°2

Code : Tout sélectionner

#Longueur = 200 
#Largeur = 200 
#Window = 0

Procedure.f min(x.f, y.f)
  If x < y
    ProcedureReturn x
  Else
    ProcedureReturn y
  EndIf
EndProcedure

Procedure.f max(x.f, y.f)
  If x > y
    ProcedureReturn x
  Else
    ProcedureReturn y
  EndIf
EndProcedure

Procedure.f float(x.l)
  ProcedureReturn x
EndProcedure

Structure array
  Adrss.l
  NbColum.l
  ElementSize.l
EndStructure

Structure conrec
 d.array; d               ! matrix of data to contour
 ilb.l; ilb,iub,jlb,jub ! index bounds of data matrix
 iub.l
 jlb.l
 jub.l
 x.array; x               ! data matrix column coordinates
 y.array; y               ! data matrix row coordinates
 nc.l; nc              ! number of contour levels
 z.array
EndStructure

xmin.f = -1.5
xmax.f = 1.5
nx = 100

ymin.f = -1.5
ymax.f = 1.5
ny = 100

nc = 20

Dim d.f(nx, ny)
Dim x.f(nx)
Dim y.f(ny)
Dim z.f(nc)
Dim colorz.f(nc)

For i = 0 To nx
 x(i) = xmin + (xmax - xmin) * i / float(nx)
Next
For j = 0 To ny
 y(j) = ymax - (ymax - ymin) * j / float(ny)
Next

zmax.f = -float(999999999)
zmin.f = float(999999999)
For i = 0 To nx
  For j = 0 To ny
    d(i,j)= 1./( Pow(x(i)*x(i)+(y(j)-0.842)*(y(j)+0.842),2) + Pow( x(i)*(y(j)-0.842) + x(i)*(y(j)-0.842),2))
    zmin = min(zmin, d(i,j))
    zmax = max(zmax, d(i,j))
  Next
Next

zmax = 2.25
zmin = 0

For k = 0 To nc
 z(k) = (zmax - zmin) * k / float(nc+1)
 colorz(k) = Int(255 * k / float(nc+1))
Next

Global xmax, xmin, ymax, ymin, zmax, zmin

Dim h.f(4)
Dim sh(4)
Dim xh.f(4)
Dim yh.f(4)

Dim im(3)
Dim jm(3)
Dim castab(2,2,2)

; Initialisation des tableaux
Restore Donnees
For i = 0 To 3
  Read im(i)
Next
For i = 0 To 3
  Read jm(i)
Next

For i = 0 To 2
  For j = 0 To 2
    For k = 0 To 2
       Read castab(i, j, k)
    Next
  Next
Next

Procedure.f xsect(p1, p2)
  ProcedureReturn (h(p2)*xh(p1)-h(p1)*xh(p2))/(h(p2)-h(p1))
EndProcedure
Procedure.f ysect(p1, p2)
  ProcedureReturn (h(p2)*yh(p1)-h(p1)*yh(p2))/(h(p2)-h(p1))
EndProcedure

Procedure vecout(x1.f, y1.f, x2.f, y2.f, k)
  x1 = #Longueur * (x1-xmin)/(xmax-xmin)
  y1 = #Largeur * (y1-ymin)/(ymax-ymin)
  x2 = #Longueur * (x2-xmin)/(xmax-xmin)
  y2 = #Largeur * (y2-ymin)/(ymax-ymin)
  
  LineXY(x1, y1, x2, y2, colorz(k))
  ;Debug Str(x1)+"/"+Str(y1)+"/"+Str(x2)+"/"+Str(y2)+"/"+Str(z)
EndProcedure

Procedure.l Array(*array.array, i, j)
 ProcedureReturn *array\Adrss + ((*array\NbColum+1) * i + j ) * *array\ElementSize
EndProcedure

Procedure.l Initconrec(*d, ilb, iub, jlb, jub, *x, *y, nc, *z)
 *map.conrec = AllocateMemory(SizeOf(conrec))
 
 *map\d\Adrss = *d
 *map\d\NbColum = jub
 *map\d\ElementSize = SizeOf(FLOAT)
 
 *map\ilb = ilb
 *map\iub = iub
 *map\jlb = jlb
 *map\jub = jub
 
 *map\x\Adrss = *x
 *map\x\NbColum = 0
 *map\x\ElementSize = SizeOf(FLOAT)
 
 *map\y\Adrss = *y
 *map\y\NbColum = 0
 *map\y\ElementSize = SizeOf(FLOAT)
 
 *map\nc = nc
 
 *map\z\Adrss = *z
 *map\z\NbColum = 0
 *map\z\ElementSize = SizeOf(FLOAT)
 ProcedureReturn *map
EndProcedure
;=============================================================================
;
;     CONREC is a contouring subroutine for rectangularily spaced data.
;
;     It emits calls to a line drawing subroutine supplied by the user
;     which draws a contour map corresponding to real*4data on a randomly
;     spaced rectangular grid. The coordinates emitted are in the same
;     units given in the x() and y() arrays.
;
;     Any number of contour levels may be specified but they must be
;     in order of increasing value.
;
;     This version is an adaptation of the original Paul Bourke subroutine:
;     http://astronomy.swin.edu.au/~pbourke/projection/conrec/ 
;
;=============================================================================
Procedure conrec(*map.conrec)
  *z0.FLOAT= Array(*map\z, 0, 0)
  *znc.FLOAT =Array(*map\z, *map\nc, 0)
      
  For j= *map\jub-1 To *map\jlb Step -1
    For i=*map\ilb To *map\iub-1
    
      *d1.FLOAT = Array(*map\d,i, j)
      *d2.FLOAT = Array(*map\d,i, j+1)
      *d3.FLOAT = Array(*map\d,i+1, j)
      *d4.FLOAT = Array(*map\d,i+1, j+1)
      
      temp1.f = min(*d1\f, *d2\f)
      temp2.f = min(*d3\f, *d4\f)
      dmin.f = min(temp1,temp2)
      temp1 = max(*d1\f, *d2\f)
      temp2 = max(*d3\f, *d4\f)
      dmax.f = max (temp1,temp2)

      If dmax>=*z0\f And dmin<=*znc\f 
        For k=0 To *map\nc
        *z.FLOAT= Array(*map\z, k, 0)
	        If *z\f>=dmin And *z\f<=dmax
            For m=4 To 0 Step -1
	           If m
             *d.FLOAT= Array(*map\d, i+im(m-1), j+jm(m-1))
             *x.FLOAT= Array(*map\x, i+im(m-1), 0)
             *y.FLOAT= Array(*map\y, j+jm(m-1), 0)
            	h(m) = *d\f-*z\f
            	xh(m) = *x\f
              yh(m) = *y\f
	           Else
             *x1.FLOAT= Array(*map\x, i, 0)
             *y1.FLOAT= Array(*map\y, j, 0)
             *x2.FLOAT = Array(*map\x, i+1, 0)
             *y2.FLOAT = Array(*map\y, j+1, 0)
           		h(0) = 0.25*(h(1)+h(2)+h(3)+h(4))
	            xh(0)=0.5*(*x1\f+*x2\f)
              yh(0)=0.5*(*y1\f+*y2\f) 
             EndIf
             
	           If h(m)>0
               sh(m) = 1
             ElseIf h(m) < 0
               sh(m) = -1
             Else
               sh(m) = 0
             EndIf
           Next
	    ;=================================================================
	    ;
	    ; Note: at this stage the relative heights of the corners and the
	    ; centre are in the h array, and the corresponding coordinates are
	    ; in the xh and yh arrays. The centre of the box is indexed by 0
	    ; and the 4 corners by 1 to 4 as shown below.
	    ; Each triangle is then indexed by the parameter m, and the 3
	    ; vertices of each triangle are indexed by parameters m1,m2,and
	    ; m3.
	    ; It is assumed that the centre of the box is always vertex 2
	    ; though this isimportant only when all 3 vertices lie exactly on
	    ; the same contour level, in which case only the side of the box
	    ; is drawn.
	    ;
	    ;
	    ;      vertex 4 +-------------------+ vertex 3
	    ;               | \               / |
	    ;               |   \    m-3    /   |
	    ;               |     \       /     |
	    ;               |       \   /       |
	    ;               |  m=2    X   m=2   |       the centre is vertex 0
	    ;               |       /   \       |
	    ;               |     /       \     |
	    ;               |   /    m=1    \   |
	    ;               | /               \ |
	    ;      vertex 1 +-------------------+ vertex 2
	    ;
	    ;
	    ;
	    ;               Scan each triangle in the box
	    ;
	    ;=================================================================             
	    For m=1 To 4
	      m1 = m
	      m2 = 0
	      If m<>4
          m3 = m+1
	      Else
          m3 = 1
        EndIf
	      
        case_value = castab(sh(m1)+1, sh(m2)+1, sh(m3)+1);
 	      If case_value<>0
    Select case_value
		  ;===========================================================
		  ;     Case 1 - Line between vertices 1 and 2
		  ;===========================================================
		Case 1
		  x1.f=xh(m1);
		  y1.f=yh(m1);
		  x2.f=xh(m2);
		  y2.f=yh(m2);
		  ;===========================================================
		  ;     Case 2 - Line between vertices 2 and 3
		  ;===========================================================
		Case 2
		  x1.f=xh(m2);
		  y1.f=yh(m2);
		  x2.f=xh(m3);
		  y2.f=yh(m3);
		  ;===========================================================
		  ;     Case 3 - Line between vertices 3 and 1
		  ;===========================================================
		Case 3
		  x1.f=xh(m3);
		  y1.f=yh(m3);
		  x2.f=xh(m1);
		  y2.f=yh(m1);
		  ;===========================================================
		  ;     Case 4 - Line between vertex 1 and side 2-3
		  ;===========================================================
		Case 4
		  x1.f=xh(m1);
		  y1.f=yh(m1);
		  x2.f=xsect(m2,m3);
		  y2.f=ysect(m2,m3);
		  ;===========================================================
		  ;     Case 5 - Line between vertex 2 and side 3-1
		  ;===========================================================
		Case 5
		  x1.f=xh(m2);
		  y1.f=yh(m2);
		  x2.f=xsect(m3,m1);
		  y2.f=ysect(m3,m1);
		  ;===========================================================
		  ;     Case 6 - Line between vertex 3 and side 1-2
		  ;===========================================================
		Case 6
		  x1.f=xh(m3);
		  y1.f=yh(m3);
		  x2.f=xsect(m1,m2);
		  y2.f=ysect(m1,m2);
		  ;===========================================================
		  ;     Case 7 - Line between sides 1-2 and 2-3
		  ;===========================================================
		Case 7
		  x1.f=xsect(m1,m2);
		  y1.f=ysect(m1,m2);
		  x2.f=xsect(m2,m3);
		  y2.f=ysect(m2,m3);
		  ;===========================================================
		  ;     Case 8 - Line between sides 2-3 and 3-1
		  ;===========================================================
		Case 8
		  x1.f=xsect(m2,m3);
		  y1.f=ysect(m2,m3);
		  x2.f=xsect(m3,m1);
		  y2.f=ysect(m3,m1);
		  ;===========================================================
		  ;     Case 9 - Line between sides 3-1 and 1-2
		  ;===========================================================
		Case 9
		  x1.f=xsect(m3,m1);
		  y1.f=ysect(m3,m1);
		  x2.f=xsect(m1,m2);
		  y2.f=ysect(m1,m2);
		EndSelect       
		;=============================================================
		; Put your graphic processing code here 
		;=============================================================
    vecout(x1,y1,x2,y2, k)
          EndIf
        
            Next              
          EndIf
        Next
      EndIf
      
    Next
  Next
    
EndProcedure

DataSection
Donnees:
Data.l 0, 1, 1, 0; im()
Data.l 0, 0, 1, 1; jm()
Data.l 0,0,8, 0,2,5, 7,6,9
Data.l 0,3,4, 1,3,1, 4,3,0
Data.l 9,6,7, 5,2,0, 8,0,0
EndDataSection

;- Debut du programme 

*map.conrec = Initconrec(@d(), 0, nx, 0, ny, @x(), @y(), nc,@z())

hWnd = OpenWindow ( #Window , 0, 0, #Longueur , #Largeur , #PB_Window_SystemMenu | #PB_Window_Invisible | #PB_Window_ScreenCentered , " Isocontour") 
CreateGadgetList (WindowID()) 
CreateImage (0, #Longueur , #Largeur )
 
StartDrawing (ImageOutput()) 
conrec(*map)
StopDrawing ()

ImageGadget (0, 0, 0, 0, 0, ImageID()) 
HideWindow ( #Window ,0) 
Repeat 
      Select WaitWindowEvent () 
           Case #PB_Event_CloseWindow 
               Quit = #True 
      EndSelect 
Until Quit
Dernière modification par Dräc le sam. 05/nov./2005 13:05, modifié 1 fois.
Haut
Avatar de l’utilisateur
djes
Messages : 4252
Inscription : ven. 11/févr./2005 17:34
Localisation : Arras, France

Message par djes »

Tu te surpasses! Ca pourra servir à un ami qui fait de la géophysique. :)
Par contre, j'ai un array index out of bound sur la première version quand le débogueur est activé.

Comme promis, j'ai donné quelques codes 3d dans la partie 3d du forum :)
Haut
Dräc
Messages : 526
Inscription : dim. 29/août/2004 0:45

Message par Dräc »

Exact, corrigé!
Haut
Répondre

Revenir à « Discussion générale »