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°1Code : 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
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
