Math - reducing dots to lines (or vectors)

[Michael Vogel]

I'm working hard to solve a problem for about two weeks and tried a lot of things now...

...however, the puzzle I want to solve contains a lot of dots within an image which should be summarized to (directed) lines automatically.

As a first step I made an example which has red dots like in the code below. A "correct" result would calculate three lines, one horizontal to the first larger dot followed by an ascending (between the larger dots) and a descending line (between the second larger dot to the last one).

Later on, the source dots wouldn't be ordered in such a perfect way. Similar results should also be found even when the distances between the dots vary (xerror and yerror) or even when an additional group of dots would be there (like an additional "line" at the top of the screen).

Actually I tried to calculate the different slopes of lines using the least square method but I didn't find some usable results for now. So maybe someone has a better idea which could help solving this riddle.

Code: Select all

xerrors=	0;10
yerrors=	0;20
epsilon.f=0.011

; Define

	Structure myData
		x.f
		y.f
	EndStructure

	Structure resData
		yest.f
		resi.f
	EndStructure

	Structure TrendData
		yint.f
		slope.f
	EndStructure

; EndDefine

Procedure leastSquare(Array mInput.myData(1),Array mOutput.resData(1),*dat.TrendData,len=0)

	Protected a.i,sumX.f,sumY.f,sumXY.f,sumXX.f,slope.f,yintercept.f,sumRes.f

	If len<1
		len = ArraySize(mInput())
	EndIf

	ReDim mOutput.resData(len)

	For a=1 To len
		sumX + mInput(a)\x
		sumY + mInput(a)\y
		sumXY + mInput(a)\x * mInput(a)\y
		sumXX + mInput(a)\x * mInput(a)\x
	Next

	slope.f = ((sumX*sumY)-(len * sumXY)) / ((sumX*sumX) - (len*sumXX));
	yintercept = (sumY - slope*sumX) / len;

	For a=1 To len
		yestimate = slope * mInput(a)\x + yintercept;
		resi = mInput(a)\y - yestimate;
		SUMres + resi*resi;
		moutput(a)\yest = yestimate
		moutput(a)\resi = resi
	Next

	*dat\slope = slope
	*dat\yint = yintercept

EndProcedure

width = 800
height = 600
points = 50

Dim MyInput.myData(points)
Dim MyOutput.resData(points)
myTrend.TrendData
oldTrend.TrendData

OpenWindow(0,0,0,width,height,"test")
StartDrawing(WindowOutput(0))
For a = 1 To points
	x = a*15+Random(xerrors)
	If a<points*0.5
		y=500+Random(yerrors)
		size=1
	ElseIf a<points*0.9
		y=500-k*10+Random(yerrors)
		k+1
		size=2
	Else
		y=500-k*10+Random(yerrors)
		k-2
		size=3
	EndIf
	MyInput(a)\x = x
	MyInput(a)\y = y
	If y > 0 And y < height And lpx
		LineXY(lpx,lpy,x,y,RGB(255,0,0))
		Circle(x,y,3+2*Bool(oldsize<>size),RGB(255,0,0))
	EndIf
	oldsize=size
	lpx=x
	lpy=y
Next

l=10
;l=points

Repeat
	leastSquare(MyInput(),MyOutput(),myTrend,l)
	draw=Bool(Abs(myTrend\slope-oldTrend\slope)>epsilon)
	If draw
		line+1
		Debug myTrend\slope-oldTrend\slope
		color=Random(#White)&$7f7f7f
	EndIf
	oldTrend\slope=myTrend\slope

	If draw
		lpx=0
		For a = 1 To l
			str.s = Str(MyInput(a)\y) +  " " + Str(MyOutput(a)\yest) + " " + Str(MyOutput(a)\resi)
			;Debug str
			y = MyOutput(a)\yest
			If y > 0 And y < height And lpx
				LineXY(lpx,lpy,myinput(a)\x,y,color)
				Circle(myinput(a)\x,y,1,color)
			EndIf
			lpx = myinput(a)\x
			lpy = y
		Next
	EndIf

	l+1
Until l>points Or line>10

StopDrawing()


Repeat
	ev = WaitWindowEvent()
Until ev = #PB_Event_CloseWindow

[jacdelad]

I...don't get the task. Since we are both german, may I ask you to describe it in german again, please?

[SMaag]

maybe this!

https://de.wikipedia.org/wiki/Douglas-P ... lgorithmus

https://docs.mcneel.com/rhino/7/help/de ... lyline.htm

here the implementaion at Rosetta Code

https://rosettacode.org/wiki/Ramer-Doug ... ion#QBasic

[jacdelad]

Ah thanks, now I get it!

[Michael Vogel]

It's not easy to describe my problem - no matter in which language - but maybe the next code does show everything a little bit clearer...

There will be a base image (which shows the yellow character 'O' in the example below) which has a lot of pixel errors (like a photo with bad quality). Anyhow I am sure that a lot of dots can be detected (red boxes with black border in the example). The number of these dots and also the density of them can vary greatly.

My aim is to find clusters of dots which can be "connected" by a line (like the magenta stripes in the example below).

Code: Select all

; Define

	RandomSeed(0)

	#Text="O"

	#FontName="Segoe UI"
	#FontSize=250
	#FontType=#PB_Font_Bold;|#PB_Font_Italic

	#WX=400
	#WY=600
	#FX=40
	#FY=0

	#ColBack=	$FFFFC0
	#ColBorder=	#Black
	#ColFill=		#Yellow

	#Draw=		$FF000000
	#Color=		$00FFFFFF

	Global Dim snw(#WX)

	Global Dim Scan(#WX)
	Global Dim Dir(#WX)
	Global Dim Clst(#WX)

; EndDefine

Procedure Init()

	LoadFont (0,#FontName,#FontSize,#FontType)

	OpenWindow(0,0,0,#WX,#WY,"Check Dots",#PB_Window_SystemMenu|#PB_Window_ScreenCentered)
	CanvasGadget(0,0,0,#WX,#WY)

	CreateImage(0,#WX,#WY,32,#Red)
	StartDrawing(ImageOutput(0))
	DrawingMode(#PB_2DDrawing_AllChannels)
	Box(0,0,#WX,#WY,#ColBack)
	StopDrawing()

	StartVectorDrawing(ImageVectorOutput(0))
	VectorFont(FontID(0))
	MovePathCursor(#FX,#FY)
	AddPathText(#Text)
	VectorSourceColor(#Draw|#ColBorder)
	StrokePath(8,#PB_Path_Preserve|#PB_Path_RoundCorner)
	VectorSourceColor(#Draw|#ColFill)
	FillPath()

	w.f

	For j=0 To 3
		Select j
		Case 0
			VectorSourceColor($ff000000)
			l=3
			n=250
			w=1.55
		Case 1
			VectorSourceColor($20000000|#ColFill)
			l=5
			n=2000
			w=1.5
		Case 2
			VectorSourceColor($80000000|#ColBack)
			l=5
			n=2000
		Case 3
			n=10000
			l=25
		EndSelect

		For i=0 To n
			x=Random(#wx)
			y=Random(#wy)
			MovePathCursor(x,y)
			AddPathLine(x+Random(l),y+Random(l))
			If j=3
				VectorSourceColor($30000000|Random(#White))
			EndIf
			StrokePath(w)
		Next i
	Next j

	StopVectorDrawing()

EndProcedure
Procedure ScanNow()

	Protected x,y

	StartDrawing(ImageOutput(0))
	DrawingFont(FontID(0))

	While x<#WX
		y=0
		While y<#WY
			If Point(x,y)&#Color=0;<>#ColBack
				;Debug Hex(#ColBack)+" - "+Hex(Point(x,y)&#Color)
				snw(x)=y
				;Box(x,0,1,y,$F3E1CA)
				;Box(x-1,y,1,1,#Red)
				Scan(x)=y
				y=#WY
			Else
				y+1
			EndIf
		Wend
		x+1
	Wend

	StopDrawing()

EndProcedure
Procedure CalcNow()

	Protected x,y
	Protected ymin,ymax

	ymin=#WY
	StartDrawing(ImageOutput(0))
	DrawingFont(FontID(0))

	While x<#WX
		If Scan(x)
			If Scan(x)>ymax : ymax=Scan(x) : EndIf
			If Scan(x)<ymin : ymin=Scan(x) : EndIf
			Box(x-3,Scan(x)-3,7,7,#Black)
			Box(x-1,Scan(x)-1,3,3,#Red)
		EndIf
		x+1
	Wend

	LineXY(0,ymin,#WX,ymin,#Blue)
	LineXY(0,ymax,#WX,ymax,#Blue)

	y=ymax
	While y>=ymin
		x=0
		While x<#WX
			If Scan(x)=y
				;Box(x,Scan(x),1,1,#White)
			EndIf
			x+1
		Wend
		y-1
	Wend

	StopDrawing()


EndProcedure
Procedure FindNow()

	StartVectorDrawing(ImageVectorOutput(0))
	VectorSourceColor($80000000|#Magenta)

	MovePathCursor(54,260)
	AddPathLine(95,190)
	
	MovePathCursor(105,180)
	AddPathLine(180,150)
	
	MovePathCursor(185,150)
	AddPathLine(270,160)
	
	MovePathCursor(150,220)
	AddPathLine(230,210)
	
	StrokePath(20)

	StopVectorDrawing()


EndProcedure

Procedure Main()

	Protected i

	Init()
	ScanNow()
	CalcNow()
	FindNow()

	StartDrawing(CanvasOutput(0))
	DrawImage(ImageID(0),0,0)

	StopDrawing()

	Repeat
		Select WaitWindowEvent()
		Case #PB_Event_CloseWindow,#WM_CHAR
			End
		EndSelect
	ForEver

EndProcedure

Main()

PS: more than years have gone since the Peucker thing has been used to keep the world running...

Code: Select all

Procedure Init()

	; My little World Version 1.5 – (c) 2013 by Michael Vogel

	#Dots=1054
	#Shapes=132

	#GlobeSize=400
	#GlobeWindow=480
	#GlobeRadius=#GlobeSize/2
	#GlobeCenter=#GlobeWindow/2
	#GlobePrecision=2<<28
	#GlobeSpace=#GlobeSize/1.33

	#Radiant=#PI/180
	#Rad090=#Radiant*90
	#Rad180=#Radiant*180
	#Rad360=#Radiant*360

	#Undefined=-#True

	#WhiteWorld=$F0FFFF
	#GreenWorld=$80FF80
	#BlueWorld=$F0F0FF
	#GrayWorld=$80808080

	#Gray  = $808080
	#Blue  = $FF0000
	#White = $FFFFFF
	#Black = $000000

	Protected i,j,count

	Global GlobeRotX=20
	Global GlobeRotY=-25

	Global FlagHidden=#True
	Global FlagMask=#True
	Global FlagFill=#True
	Global FlagGrid=#True

	Structure DotType
		x.i
		y.i
		x_.i
		y_.i
		hidden.i
	EndStructure

	Structure ShapeType
		len.i
		start.i
	EndStructure

	Global Dim Dots.DotType(#Dots)
	Global Dim Shapes.ShapeType(#Shapes)

	Global Dim Polygon.l(666)
	Global PolygonCount
	Global PolygonDC

	count=0
	For i=0 To #Shapes
		Shapes(i)\start=count
		Read.a Shapes(i)\len

		For j=1 To Shapes(i)\len
			Read.a Dots(count)\x
			Read.a Dots(count)\y
			Dots(count)\x=(Dots(count)\x/256.0*2-1)*#PI*#GlobePrecision
			Dots(count)\y=(Dots(count)\y/256.0-0.5)*#PI*#GlobePrecision
			count+1
		Next j
	Next i

	DataSection
		Data.a 128,255,35,254,37,249,41,243,43,242,49,239,54,237,48,241,43,240,40,235,44,229,43,225,47,226,51,227,57,224,63,220,67,218,73,217,79,216,72,211,73,212,78,213,82,212,85,213,88,212,91,210,94,208,96,205
		Data.a 98,203,98,202,101,204,105,204,110,200,111,198,112,199,118,201,124,198,119,197,113,197,107,194,102,192,96,187,100,184,106,184,112,181,113,180,107,179,101,176,96,172,91,167,89,163,85,163,89,167,90
		Data.a 169,95,167,101,163,106,158,109,157,102,155,96,153,90,150,85,152,79,150,76,146,72,150,68,156,69,154,64,152,64,148,64,144,70,141,71,139,66,136,65,140,70,138,70,135,65,129,68,126,74,121,74,121,68,125
		Data.a 66,126,61,127,58,131,53,133,48,137,51,142,48,144,43,145,42,143,38,143,34,140,39,139,45,136,46,133,45,132,39,137,34,141,29,147,27,152,29,156,33,152,36,157,33,160,33,165,30,171,29,175,28,178,27,178
		Data.a 33,183,32,179,27,182,27,185,28,187,23,192,20,198,19,203,19,208,21,202,24,208,24,214,24,219,25,224,26,230,25,235,27,241,29,247,30,254,29,123,11,38,11,42,13,44,15,47,15,47,19,44,21,41,22,42,27,42
		Data.a 31,45,34,48,36,52,39,56,40,59,40,64,40,70,41,75,44,80,46,86,48,92,49,92,47,85,48,86,51,92,53,99,56,103,61,105,66,110,69,116,72,118,72,124,70,130,70,136,72,143,75,150,78,155,77,161,77,167,77,173
		Data.a 76,179,75,185,75,192,75,197,76,203,79,198,80,192,81,186,85,182,86,177,90,175,92,169,94,162,98,157,100,151,100,145,103,138,100,133,96,130,93,129,92,125,89,119,84,114,79,112,76,113,74,113,71,114,68
		Data.a 111,68,105,65,101,64,97,60,102,58,96,58,90,61,85,66,84,69,89,70,87,70,81,73,76,74,74,76,69,79,64,81,64,83,63,82,59,78,60,81,56,86,55,86,51,85,49,82,44,79,44,76,41,72,41,73,47,71,52,69,54,68,49,63
		Data.a 46,60,42,62,37,65,36,68,34,69,29,66,32,62,29,59,27,60,32,56,31,51,30,49,31,45,30,39,29,35,29,33,29,27,29,21,28,16,27,11,30,14,33,8,34,13,36,12,38,59,145,82,141,82,139,83,135,80,135,76,131,75,127
		Data.a 77,123,78,121,82,120,86,117,90,116,95,116,99,116,103,116,108,116,112,118,116,121,120,125,120,129,118,131,121,134,124,134,128,135,132,136,136,137,140,137,144,136,149,135,153,137,157,137,161,138,165
		Data.a 139,169,140,174,143,175,147,174,149,170,150,166,152,162,152,158,153,153,156,149,156,145,155,141,155,136,156,131,158,127,160,123,162,118,163,114,161,112,157,111,156,107,155,102,154,97,152,93,151
		Data.a 89,150,85,147,83,31,228,143,227,148,226,152,223,149,224,145,220,145,218,149,216,149,214,153,211,156,208,159,208,163,208,167,209,172,209,176,213,175,217,173,221,172,223,176,224,176,226,181,231,181
		Data.a 233,179,235,174,236,170,236,166,235,162,233,157,231,153,230,148,228,144,17,123,57,125,56,127,55,128,54,127,53,127,51,126,49,126,47,125,46,124,44,123,46,123,48,125,50,125,52,124,53,125,54,124,56
		Data.a 11,113,37,112,36,112,35,110,35,112,34,113,34,115,34,117,33,118,35,116,37,114,38,10,232,139,230,134,227,131,223,132,221,128,221,132,225,135,227,140,231,140,232,140,8,205,125,208,121,210,117,210,122
		Data.a 210,127,209,132,205,129,205,125,9,249,176,251,178,252,180,253,182,252,184,251,184,251,182,251,180,250,178,6,249,185,250,188,248,191,245,192,247,189,249,186,9,198,123,200,126,202,130,202,135,199
		Data.a 132,198,127,196,124,196,120,198,123,7,122,50,120,51,121,53,121,54,123,53,123,51,123,49,10,228,71,227,73,227,78,224,79,222,79,221,77,224,76,226,74,226,70,228,70,10,87,57,86,59,86,60,88,60,89,60,90
		Data.a 60,89,58,88,57,88,55,87,56,15,213,125,215,126,216,126,213,127,212,128,214,128,214,131,214,133,213,133,212,133,212,135,212,133,211,131,212,128,212,126,16,163,147,163,149,162,152,162,154,161,157,161
		Data.a 159,161,162,159,163,158,160,158,158,159,155,158,152,159,150,161,148,161,146,162,146,4,79,202,81,205,78,203,79,202,6,74,99,71,97,69,95,69,94,71,96,74,98,7,207,139,205,138,203,137,203,136,205,137
		Data.a 207,138,208,139,9,214,101,214,104,213,106,214,108,214,108,213,107,212,104,212,102,214,101,5,7,34,5,33,4,33,3,32,1,31,5,67,34,69,36,68,37,66,37,66,34,11,229,57,228,58,228,61,228,62,228,60,228,57
		Data.a 228,55,228,53,228,51,229,53,229,55,8,216,114,217,116,217,118,216,118,215,117,214,117,214,115,216,114,4,39,57,40,58,38,57,37,56,2,18,46,19,46,7,230,65,229,67,228,67,227,68,227,67,227,65,228,64,5
		Data.a 78,101,78,100,77,99,76,100,77,102,5,230,185,232,185,232,187,230,188,230,185,4,220,79,219,81,220,83,221,81,3,244,158,244,157,244,158,4,233,135,235,133,235,135,233,136,3,7,34,6,35,4,36,3,4,35,2,35
		Data.a 0,35,2,34,51,33,51,2,7,38,6,37,3,134,71,133,71,134,69,4,138,73,138,75,136,74,138,73,4,218,124,218,126,217,127,218,124,3,215,140,216,140,215,140,2,85,200,85,200,2,85,62,84,62,3,212,139,214,139,213
		Data.a 140,3,143,73,143,75,143,74,4,211,113,210,115,211,113,211,113,3,219,132,219,132,219,132,5,213,92,213,94,213,96,212,94,213,92,5,184,113,184,116,184,118,185,117,184,114,3,134,66,134,69,134,67,2,145
		Data.a 77,145,78,3,234,131,236,133,235,133,2,74,202,74,202,2,17,99,17,101,3,206,99,205,101,204,99,3,222,80,221,80,222,80,2,215,110,216,111,2,253,152,253,152,2,163,30,161,30,2,11,50,11,50,2,143,44,144,45
		Data.a 2,233,63,232,63,2,145,73,144,73,2,245,148,245,148,13,163,74,165,75,165,71,166,70,164,68,163,65,165,63,165,61,162,62,161,65,162,68,162,71,162,74,8,71,64,68,62,65,63,65,64,65,68,66,65,68,64,70,64
		Data.a 9,50,39,48,39,46,39,46,39,45,40,43,41,45,41,47,40,49,39,3,42,33,42,35,43,33,4,58,50,59,54,58,54,58,51,8,68,62,67,59,65,58,64,60,62,61,64,61,65,61,67,61,4,49,43,51,43,51,44,49,44,5,151,127,151,129
		Data.a 151,131,150,129,150,127,4,74,64,72,65,72,66,73,65,5,170,66,171,64,171,62,169,63,169,65,7,205,49,205,51,204,53,202,54,202,54,203,52,204,50,4,71,67,69,67,69,69,71,67,5,180,61,182,61,180,62,180,64
		Data.a 180,62,3,149,40,149,42,150,41,6,149,139,149,137,148,134,148,134,148,137,149,139,4,152,143,152,143,151,145,152,147,2,47,68,47,68,2,137,44,137,44,3,138,109,137,107,137,109,69,0,238,6,238,11,238,16
		Data.a 236,20,236,24,235,28,233,32,233,36,232,40,231,44,232,48,232,53,232,54,231,57,229,61,230,65,230,70,231,74,231,75,229,77,225,80,223,83,219,85,220,84,225,85,229,84,233,89,236,95,238,101,238,107,235
		Data.a 111,232,115,230,120,228,124,227,128,226,133,226,137,225,142,226,147,226,152,225,156,225,160,223,164,221,168,222,172,223,177,224,182,225,185,221,189,221,193,222,197,220,202,221,206,221,210,221,213
		Data.a 220,217,222,221,221,225,221,230,222,234,224,238,225,242,227,246,228,245,232,242,235,244,237,249,237,253,238,48,106,31,110,29,110,28,109,27,112,26,108,24,112,23,112,21,114,19,113,15,113,13,117,13
		Data.a 114,12,110,14,109,12,104,12,108,11,112,11,108,10,104,10,104,10,108,10,103,9,99,9,95,10,95,10,91,11,87,11,84,12,80,14,80,15,76,16,80,18,79,20,84,20,87,22,88,26,91,28,91,28,89,32,90,36,91,38,94,41
		Data.a 98,40,99,36,102,34,104,30,106,31,13,69,14,68,13,62,12,68,11,74,10,80,10,79,12,79,13,73,15,68,19,68,18,69,16,69,14,14,79,38,80,36,83,33,79,29,76,26,70,25,67,25,65,23,65,28,71,28,76,31,72,35,77,38
		Data.a 80,38,10,42,25,44,24,45,23,43,22,41,22,40,23,39,25,39,26,41,26,42,25,5,66,14,63,17,63,15,60,13,64,13,5,139,14,140,18,139,17,136,16,138,14,9,56,28,53,29,48,30,44,28,45,27,44,25,49,25,53,25,56,28
		Data.a 3,143,14,144,15,141,14,5,49,20,52,20,48,22,45,21,47,20,6,170,20,176,19,172,21,168,23,167,22,168,20,2,161,14,161,13,3,57,20,57,20,56,19,5,65,20,69,20,69,22,65,22,62,20,5,57,24,57,26,58,24,57,23,57
		Data.a 24,3,143,17,143,18,142,17,5,165,24,167,25,168,27,165,26,165,24,5,57,16,55,17,53,17,54,15,56,16,3,60,18,61,19,63,20,4,60,23,62,23,61,24,60,23,3,160,13,160,13,160,13,4,42,20,44,19,45,18,43,18,2,170
		Data.a 13,170,14,4,194,14,193,13,195,12,195,14,4,198,16,199,15,201,16,200,17,3,60,17,58,17,59,16,2,172,12,172,12,4,198,15,197,16,195,15,193,15,3,50,17,47,17,49,18,3,60,20,61,22,59,21,3,167,13,167,13,167
		Data.a 13,3,226,21,224,21,225,20,4,59,29,57,30,59,30,59,29,3,49,16,47,16,49,16,2,151,14,151,14,3,226,20,228,20,227,21,2,60,17,60,18,4,227,20,229,20,229,21,228,20,3,42,20,40,20,42,19,3,233,22,232,21,233
		Data.a 22,2,192,14,192,14,2,135,16,135,16,3,227,23,227,24,227,23,3,93,10,95,11,93,10,2,53,24,52,23,3,90,28,90,29,89,28,4,71,23,73,24,71,24,71,23,2,0,26,1,27,2,111,24,111,25,2,109,24,109,24,2,112,24,112,24
	EndDataSection

EndProcedure
Procedure PolygonDot(x,y,mode=#True)

	If mode
		Polygon(PolygonCount)=x
		PolygonCount+1
		Polygon(PolygonCount)=y
		PolygonCount+1
	Else
		PolygonCount=0
	EndIf

EndProcedure
Procedure PolygonCheck(x1,y1,x2,y2)

	Protected a1,a2
	Protected d1,d2
	Protected fall

	a1=ATan2(x1-#GlobeCenter,y1-#GlobeCenter)/#Radiant
	a2=ATan2(x2-#GlobeCenter,y2-#GlobeCenter)/#Radiant

	d1=Abs(a1-a2)

	If a1<0
		a1+360
	EndIf
	If a2<0
		a2+360
	EndIf

	d2=Abs(a1-a2)
	If d2>d1
		d2=#True
	Else
		d1=d2
		d2=#False
	EndIf

	If d1>090
		If d2
			If a1>180
				a1-360
			EndIf
			If a2>180
				a2-360
			EndIf
		EndIf

		y1=Sin((a1+a2)/2*#Radiant)*#GlobeSize+#GlobeCenter
		x1=Cos((a1+a2)/2*#Radiant)*#GlobeSize+#GlobeCenter
		PolygonDot(x1,y1)

	EndIf

EndProcedure
Procedure GlobeDot(x.f,y.f,color)

	Protected rx.f,ry.f
	Protected x_.f,y_.f

	rx=GlobeRotX*#Radiant
	ry=GlobeRotY*#Radiant

	x*#Radiant
	y*#Radiant

	x_=Sin(x-rx)*Cos(y)
	y_=Sin(y)*Cos(ry)-Cos(x-rx)*Cos(y)*Sin(ry)

	If FlagHidden=0 Or Cos(x-rx)*Cos(y)*Cos(ry)+Sin(y)*Sin(ry)>=0
		Plot(x_*#GlobeRadius+#GlobeCenter,y_*#GlobeRadius+#GlobeCenter,color)
	EndIf

EndProcedure
Procedure GlobeDraw()

	Protected i,j,n,start,stop
	Protected x.f,y.f,rx.f,ry.f

	Enumeration
		#FlagInvisible
		#FlagInit
		#FlagVisible
	EndEnumeration

	rx=GlobeRotX*#Radiant
	ry=GlobeRotY*#Radiant

	For i=0 To #Dots
		x=Dots(i)\x/#GlobePrecision
		y=Dots(i)\y/#GlobePrecision
		If FlagHidden And Cos(x-rx)*Cos(y)*Cos(ry)+Sin(y)*Sin(ry)<0
			Dots(i)\hidden=#True
		Else
			Dots(i)\hidden=#False
		EndIf
		Dots(i)\x_=(Sin(x-rx)*Cos(y))*#GlobeRadius
		Dots(i)\y_=(Sin(y)*Cos(ry)-Cos(x-rx)*Cos(y)*Sin(ry))*#GlobeRadius
	Next i

	PolygonDC=StartDrawing(CanvasOutput(0))
	Circle(#GlobeCenter,#GlobeCenter,#GlobeRadius,#BlueWorld)
	DrawingMode(#PB_2DDrawing_Transparent)
	
	Protected DotFirst,DotLast,DotMode

	If FlagFill

		For i=0 To #Shapes

			n=#False
			start=Shapes(i)\start
			stop=start+Shapes(i)\len-1

			For j=start To stop
				If Dots(j)\hidden
					x.f=Sqr(Dots(j)\x_*Dots(j)\x_+Dots(j)\y_*Dots(j)\y_)
					If x
						x=#GlobeSpace/x
						Dots(j)\x_*x
						Dots(j)\y_*x
					EndIf
				Else
					n+1
				EndIf
				Dots(j)\x_+#GlobeCenter
				Dots(j)\y_+#GlobeCenter
			Next j

			If n
				DotFirst=#Undefined
				DotLast=#Undefined
				DotMode=#True

				PolygonDot(#Undefined,#Undefined,#False)

				For j=start To stop

					If Dots(j)\hidden
						If DotMode
							DotLast=j
						Else
							PolygonDot(Dots(j)\x_,Dots(j)\y_)
							DotFirst=j
							DotLast=#Undefined
							DotMode=#True
						EndIf

					Else
						If DotLast>=0
							If DotFirst>=0
								PolygonCheck(Dots(DotFirst)\x_,Dots(DotFirst)\y_,Dots(DotLast)\x_,Dots(DotLast)\y_)
							EndIf
							PolygonDot(Dots(DotLast)\x_,Dots(DotLast)\y_)
						EndIf
						PolygonDot(Dots(j)\x_,Dots(j)\y_)
						DotLast=#Undefined;
						DotFirst=#Undefined
						DotMode=#False
					EndIf
				Next j

				If Dots(start)\hidden=#False
					If DotLast>=0
						PolygonCheck(Polygon(PolygonCount-2),Polygon(PolygonCount-1),Dots(stop)\x_,Dots(stop)\y_)
						PolygonDot(Dots(stop)\x_,Dots(stop)\y_)
					EndIf

				Else;			
					If DotMode=#False
						PolygonDot(Dots(start)\x_,Dots(start)\y_)
						PolygonCheck(Polygon(0),Polygon(1),Dots(start)\x_,Dots(start)\y_)
					Else;		
						PolygonCheck(Polygon(0),Polygon(1),Polygon(PolygonCount-2),Polygon(PolygonCount-1))
					EndIf
				EndIf

				If i<63
					n=#GreenWorld
				ElseIf i<82
					n=#BlueWorld
				Else
					n=#WhiteWorld
				EndIf

				n=CreateSolidBrush_(n)
				SelectObject_(PolygonDC,n)
				Polygon_(PolygonDC,@Polygon(),PolygonCount>>1)
				DeleteObject_(n)
			EndIf

		Next i

	Else

		For i=0 To #Shapes
			start=Shapes(i)\start
			stop=start+Shapes(i)\len-1
			For j=start To stop-1
				If FlagHidden=0 Or Dots(j)\hidden+Dots(j+1)\hidden=#False
					LineXY(Dots(j)\x_+#GlobeCenter,Dots(j)\y_+#GlobeCenter,Dots(j+1)\x_+#GlobeCenter,Dots(j+1)\y_+#GlobeCenter,$FFFF0000)
				EndIf
			Next j
		Next i

	EndIf

	If FlagGrid

		DrawingMode(#PB_2DDrawing_AlphaBlend)

		#GridLarge=15
		#GridSmall=1

		i=-180-#GridLarge
		While i<180
			i+#GridLarge
			j=-90-#GridSmall
			While j<90
				j+#GridSmall
				GlobeDot(i,j,#GrayWorld)
			Wend
		Wend

		i=-90-#GridLarge
		While i<90
			i+#GridLarge
			n=Abs(i)/30+1
			j=-180-n
			While j<180
				j+n
				GlobeDot(j,i,#GrayWorld)
			Wend
		Wend

	EndIf

	If FlagMask
		DrawingMode(#PB_2DDrawing_AlphaBlend)
		DrawImage(ImageID(0),0,0)
	EndIf

	StopDrawing()

EndProcedure

Procedure Main()

	Init()

	OpenWindow(0,0,0,#GlobeCenter<<1,#GlobeCenter<<1,"World by Michael Vogel  •  use cursor keys, 'space', 'f', 'g' and 'm'...",#PB_Window_ScreenCentered)
	CanvasGadget(0,0,0,#GlobeCenter<<1,#GlobeCenter<<1)

	AddKeyboardShortcut(0,#PB_Shortcut_Escape,#PB_Shortcut_Escape)
	AddKeyboardShortcut(0,#PB_Shortcut_Left,#PB_Shortcut_Left)
	AddKeyboardShortcut(0,#PB_Shortcut_Right,#PB_Shortcut_Right)
	AddKeyboardShortcut(0,#PB_Shortcut_Up,#PB_Shortcut_Up)
	AddKeyboardShortcut(0,#PB_Shortcut_Down,#PB_Shortcut_Down)
	AddKeyboardShortcut(0,#PB_Shortcut_Space,#PB_Shortcut_Space)
	AddKeyboardShortcut(0,#PB_Shortcut_F,#PB_Shortcut_F)
	AddKeyboardShortcut(0,#PB_Shortcut_G,#PB_Shortcut_G)
	AddKeyboardShortcut(0,#PB_Shortcut_M,#PB_Shortcut_M)

	CreateImage(0,#GlobeWindow,#GlobeWindow,32)
	StartDrawing(ImageOutput(0))
	DrawingMode(#PB_2DDrawing_AllChannels)
	Box(0,0,#GlobeCenter<<1,#GlobeCenter<<1,$FFffe0e0)
	DrawingMode(#PB_2DDrawing_Gradient|#PB_2DDrawing_AllChannels)

	BackColor($40FFFFFF)
	GradientColor(0.2,$40FFFFFF)
	GradientColor(0.75,$40000000)
	FrontColor($40000000)

	CircularGradient(#GlobeSize/3+#GlobeCenter-#GlobeRadius,#GlobeSize/3+#GlobeCenter-#GlobeRadius,#GlobeSize)
	Circle(#GlobeCenter,#GlobeCenter,#GlobeRadius)

	DrawingMode(#PB_2DDrawing_AlphaBlend|#PB_2DDrawing_Outlined)
	Circle(#GlobeCenter,#GlobeCenter,#GlobeRadius,$10000000)
	Circle(#GlobeCenter,#GlobeCenter,#GlobeRadius+1,$10000000)
	StopDrawing()

	GlobeDraw()

	Repeat
		Select WaitWindowEvent()
		Case #PB_Event_CloseWindow
			End
		Case #PB_Event_Menu
			Select EventGadget()
			Case #PB_Shortcut_Escape
				End
			Case #PB_Shortcut_Left
				GlobeRotX+5
			Case #PB_Shortcut_Right
				GlobeRotX-5
			Case #PB_Shortcut_Up
				GlobeRotY+5
			Case #PB_Shortcut_Down
				GlobeRotY-5
			Case #PB_Shortcut_Space
				FlagHidden!1
			Case #PB_Shortcut_F
				FlagFill!1
			Case #PB_Shortcut_G
				FlagGrid!1
			Case #PB_Shortcut_M
				FlagMask!1
			EndSelect
			GlobeDraw()

		EndSelect

	ForEver

EndProcedure

Main()

[SMaag]

That sounds like the Conrec Algorithm!

"Conrec is an algorithm for contouring surfaces (in other words, it finds the contour lines of a two-dimensional function)"

here is my base implementation of Conrec. But until now completely untested.

https://paulbourke.net/papers/conrec/

Code: Select all

EnableExplicit

Procedure Conrec(Array z.d(2), Array x.d(1), Array y.d(1), nc, Array contour.d(1), ilb, iub, jlb, jub, Array color.l(1))


;from: 
; http://astronomy.swin.edu.au/~pbourke/projection/conrec/conrec_vb.txt
; https://paulbourke.net/papers/conrec/

    ;=============================================================================
    ;     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*4 (double) data
    ;     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.
    ;
    ;     adapted from the fortran-77 routine CONREC.F developed by Paul D. Bourke
    ;=============================================================================
    
    ; Z(#,#)          ; matrix of Data To contour
    ; ilb,iub,jlb,jub ; index bounds of Data matrix (x-lower,x-upper,y-lower,y-upper)
    ; X(#)            ; Data matrix column coordinates
    ; Y(#)            ; Data matrix row coordinates
    ; nc              ; number of contour levels
    ; contour(#)      ; contour levels in increasing order
                         
    ; Dim m1, m2, m3, case_value As Integer
    Protected.i m1, m2, m3, case_value
    
    ;Dim dmin, dmax As Double
    ;Dim x1, x2, y1, y2 As Double
    Protected.d x1, x2, y1, y2, dmin, dmax
    
    ;Dim i, j, k, m As Integer
    Protected.i i, j, k, m
    
    ; Dim h(5) As Double
    ; Dim sh(5) As Integer
    ; Dim xh(5), yh(5) As Double
    
    Dim  h.d(5)
    Dim xh.d(5)
    Dim yh.d(5)
    Dim sh.i(5)
    
    Protected Nullcode.d = 1E+37
    
    Dim im.i(4)
    Dim jm.i(4)
    
    im(0) = 0
    im(1) = 1
    im(2) = 1
    im(3) = 0
    jm(0) = 0
    jm(1) = 0
    jm(2) = 1
    jm(3) = 1
    
    Dim castab.i(3, 3, 3)
    castab(0, 0, 0) = 0
    castab(0, 0, 1) = 0
    castab(0, 0, 2) = 8 
    castab(0, 1, 0) = 0
    castab(0, 1, 1) = 2
    castab(0, 1, 2) = 5 
    castab(0, 2, 0) = 7
    castab(0, 2, 1) = 6
    castab(0, 2, 2) = 9 
    castab(1, 0, 0) = 0
    castab(1, 0, 1) = 3
    castab(1, 0, 2) = 4 
    castab(1, 1, 0) = 1
    castab(1, 1, 1) = 3
    castab(1, 1, 2) = 1 
    castab(1, 2, 0) = 4
    castab(1, 2, 1) = 3
    castab(1, 2, 2) = 0 
    castab(2, 0, 0) = 9
    castab(2, 0, 1) = 6
    castab(2, 0, 2) = 7 
    castab(2, 1, 0) = 5
    castab(2, 1, 1) = 2
    castab(2, 1, 2) = 0 
    castab(2, 2, 0) = 8
    castab(2, 2, 1) = 0
    castab(2, 2, 2) = 0 
    
    If nc <> 0 
      For j = jub - 1 To jlb Step -1
        For i = ilb To iub - 1
          Protected.d temp1, temp2
          
          ; temp1 = Min(z(i, j), z(i, j + 1))
          If z(i, j) < z(i, j+1)
            temp1 = z(i, j)
          Else
            temp1 = z(i, j+1)
          EndIf
            
          ; temp2 = Min(z(i + 1, j), z(i + 1, j + 1))
          If z(i+1, j) < z(i+1, j+1)
            temp2 =  z(i+1, j)
          Else
            temp2 = z(i+1, j+1)
          EndIf
          
          ;dmin = Min(temp1, temp2)
          If temp1 < temp2
            dmin = temp1
          Else
            dmin = temp2
          EndIf
          
          
          ;temp1 = Max(z(i, j), z(i, j + 1))
          If z(i, j) > z(i, j+1)
            temp1 = z(i, j)
          Else
            temp1 = z(i, j+1)
          EndIf
         
          ;temp2 = Max(z(i + 1, j), z(i + 1, j + 1))
          If z(i+1, j) > z(i+1, j+1)
            temp2 =  z(i+1, j)
          Else
            temp2 = z(i+1, j+1)
          EndIf
          
          ;dmax = Max(temp1, temp2)
          If temp1 > temp2
            dmax = temp1
          Else
            dmax = temp2
          EndIf
          
    ;-------------------------------------------------------------------------
           ;extra conditional added here to insure that large values are not plotted
           ;if an area should not be contoured, values above nullcode should be entered in
           ;the matrix Z
          
    ;------------------------------------------------------------------------
          If dmax >= contour(0) And dmin <= contour(nc - 1) And dmax < Nullcode
            For k = 0 To nc - 1
              If contour(k) >= dmin And contour(k) < dmax 
                For m = 4 To 0 Step -1
                  
                  If m >0
                    h(m) = z(i + im(m - 1), j + jm(m - 1)) - contour(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 m
               
                ;=================================================================
                ;
                ; 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
                        x1 = xh(m1)
                        y1 = yh(m1)
                        x2 = xh(m2)
                        y2 = yh(m2)
                         
                      Case 2
                        ; Line between vertices 2 and 3
                        x1 = xh(m2)
                        y1 = yh(m2)
                        x2 = xh(m3)
                        y2 = yh(m3)
                        
                      Case 3
                        ;Line between vertices 3 and 1 
                        x1 = xh(m3)
                        y1 = yh(m3)
                        x2 = xh(m1)
                        y2 = yh(m1)
                         
                      Case 4
                        ; Line between vertex 1 and side 2-3
                        x1 = xh(m1)
                        y1 = yh(m1)
                        x2 = (h(m3) * xh(m2) - h(m2) * xh(m3)) / (h(m3) - h(m2))
                        y2 = (h(m3) * yh(m2) - h(m2) * yh(m3)) / (h(m3) - h(m2))
                           
                      Case 5
                        ; Line between vertex 2 and side 3-1  
                        x1 = xh(m2)
                        y1 = yh(m2)
                        x2 = (h(m1) * xh(m3) - h(m3) * xh(m1)) / (h(m1) - h(m3))
                        y2 = (h(m1) * yh(m3) - h(m3) * yh(m1)) / (h(m1) - h(m3))
                         
                      Case 6
                        ; Line between vertex 3 and side 1-2 
                        x1 = xh(m3)
                        y1 = yh(m3)
                        x2 = (h(m2) * xh(m1) - h(m1) * xh(m2)) / (h(m2) - h(m1))
                        y2 = (h(m2) * yh(m1) - h(m1) * yh(m2)) / (h(m2) - h(m1))
                         
                      Case 7
                        ; Line between sides 1-2 and 2-3  
                        x1 = (h(m2) * xh(m1) - h(m1) * xh(m2)) / (h(m2) - h(m1))
                        y1 = (h(m2) * yh(m1) - h(m1) * yh(m2)) / (h(m2) - h(m1))
                        x2 = (h(m3) * xh(m2) - h(m2) * xh(m3)) / (h(m3) - h(m2))
                        y2 = (h(m3) * yh(m2) - h(m2) * yh(m3)) / (h(m3) - h(m2))
                           
                      Case 8
                        ; Line between sides 2-3 and 3-1 
                        x1 = (h(m3) * xh(m2) - h(m2) * xh(m3)) / (h(m3) - h(m2))
                        y1 = (h(m3) * yh(m2) - h(m2) * yh(m3)) / (h(m3) - h(m2))
                        x2 = (h(m1) * xh(m3) - h(m3) * xh(m1)) / (h(m1) - h(m3))
                        y2 = (h(m1) * yh(m3) - h(m3) * yh(m1)) / (h(m1) - h(m3))
                         
                      Case 9
                        ; Line between sides 3-1 and 1-2
                        x1 = (h(m1) * xh(m3) - h(m3) * xh(m1)) / (h(m1) - h(m3))
                        y1 = (h(m1) * yh(m3) - h(m3) * yh(m1)) / (h(m1) - h(m3))
                        x2 = (h(m2) * xh(m1) - h(m1) * xh(m2)) / (h(m2) - h(m1))
                        y2 = (h(m2) * yh(m1) - h(m1) * yh(m2)) / (h(m2) - h(m1))
                        
                    EndSelect
                    ;--------------------------------------------------------------
                    ;this is where the program specific drawing routine comes in.
                    ;This specific command will work well for a properly dimensioned
                    ;vb picture box or vb form (where "Form1" is the name of the form)
                       
                    ;-------------------------------------------------------------------
                        LineXY(x1, y1, x2, y2, color(k))
                       
                    ;-------------------------------------------------------------------
                  EndIf
                Next m
              EndIf
            Next k
          EndIf
         Next i      
       Next j
    EndIf
EndProcedure

[Michael Vogel]

Thanks for this one, never heard before from Conrec

Anyhow, do you have an idea how to fill the x, y and z arrays? I thought to put the interesting dots as "1" into the z matrix ad leave all other values equal "0". But what to do with x and y?

Here's the modified code - just pasted your conrec procedure into it and created some arrays. FindNow is filling the z array and calling the conrec procedure.

Code: Select all

; Define

	RandomSeed(0)

	#Text="O"

	#FontName="Segoe UI"
	#FontSize=250
	#FontType=#PB_Font_Bold;|#PB_Font_Italic

	#WX=400
	#WY=600
	#FX=40
	#FY=0

	#ColBack=	$FFFFC0
	#ColBorder=	#Black
	#ColFill=		#Yellow

	#Draw=		$FF000000
	#Color=		$00FFFFFF

	Global Dim snw(#WX)

	Global Dim Scan(#WX)
	Global Dim Dir(#WX)
	Global Dim Clst(#WX)

; EndDefine
Procedure Conrec(Array z.d(2), Array x.d(1), Array y.d(1), nc, Array contour.d(1), ilb, iub, jlb, jub, Array color.i(1))


	;from:
	; http://astronomy.swin.edu.au/~pbourke/projection/conrec/conrec_vb.txt
	; https://paulbourke.net/papers/conrec/

	;=============================================================================
	;     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*4 (double) data
	;     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.
	;
	;     adapted from the fortran-77 routine CONREC.F developed by Paul D. Bourke
	;=============================================================================

	; Z(#,#)          ; matrix of Data To contour
	; ilb,iub,jlb,jub ; index bounds of Data matrix (x-lower,x-upper,y-lower,y-upper)
	; X(#)            ; Data matrix column coordinates
	; Y(#)            ; Data matrix row coordinates
	; nc              ; number of contour levels
	; contour(#)      ; contour levels in increasing order

	; Dim m1, m2, m3, case_value As Integer
	Protected.i m1, m2, m3, case_value

	;Dim dmin, dmax As Double
	;Dim x1, x2, y1, y2 As Double
	Protected.d x1, x2, y1, y2, dmin, dmax

	;Dim i, j, k, m As Integer
	Protected.i i, j, k, m

	; Dim h(5) As Double
	; Dim sh(5) As Integer
	; Dim xh(5), yh(5) As Double

	Dim  h.d(5)
	Dim xh.d(5)
	Dim yh.d(5)
	Dim sh.i(5)

	Protected Nullcode.d = 1E+37

	Dim im.i(4)
	Dim jm.i(4)

	im(0) = 0
	im(1) = 1
	im(2) = 1
	im(3) = 0
	jm(0) = 0
	jm(1) = 0
	jm(2) = 1
	jm(3) = 1

	Dim castab.i(3, 3, 3)
	castab(0, 0, 0) = 0
	castab(0, 0, 1) = 0
	castab(0, 0, 2) = 8
	castab(0, 1, 0) = 0
	castab(0, 1, 1) = 2
	castab(0, 1, 2) = 5
	castab(0, 2, 0) = 7
	castab(0, 2, 1) = 6
	castab(0, 2, 2) = 9
	castab(1, 0, 0) = 0
	castab(1, 0, 1) = 3
	castab(1, 0, 2) = 4
	castab(1, 1, 0) = 1
	castab(1, 1, 1) = 3
	castab(1, 1, 2) = 1
	castab(1, 2, 0) = 4
	castab(1, 2, 1) = 3
	castab(1, 2, 2) = 0
	castab(2, 0, 0) = 9
	castab(2, 0, 1) = 6
	castab(2, 0, 2) = 7
	castab(2, 1, 0) = 5
	castab(2, 1, 1) = 2
	castab(2, 1, 2) = 0
	castab(2, 2, 0) = 8
	castab(2, 2, 1) = 0
	castab(2, 2, 2) = 0

	If nc <> 0
		For j = jub - 1 To jlb Step -1
			For i = ilb To iub - 1
				Protected.d temp1, temp2

				; temp1 = Min(z(i, j), z(i, j + 1))
				If z(i, j) < z(i, j+1)
					temp1 = z(i, j)
				Else
					temp1 = z(i, j+1)
				EndIf

				; temp2 = Min(z(i + 1, j), z(i + 1, j + 1))
				If z(i+1, j) < z(i+1, j+1)
					temp2 =  z(i+1, j)
				Else
					temp2 = z(i+1, j+1)
				EndIf

				;dmin = Min(temp1, temp2)
				If temp1 < temp2
					dmin = temp1
				Else
					dmin = temp2
				EndIf


				;temp1 = Max(z(i, j), z(i, j + 1))
				If z(i, j) > z(i, j+1)
					temp1 = z(i, j)
				Else
					temp1 = z(i, j+1)
				EndIf

				;temp2 = Max(z(i + 1, j), z(i + 1, j + 1))
				If z(i+1, j) > z(i+1, j+1)
					temp2 =  z(i+1, j)
				Else
					temp2 = z(i+1, j+1)
				EndIf

				;dmax = Max(temp1, temp2)
				If temp1 > temp2
					dmax = temp1
				Else
					dmax = temp2
				EndIf

				;-------------------------------------------------------------------------
				;extra conditional added here to insure that large values are not plotted
				;if an area should not be contoured, values above nullcode should be entered in
				;the matrix Z

				;------------------------------------------------------------------------
				If dmax >= contour(0) And dmin <= contour(nc - 1) And dmax < Nullcode
					For k = 0 To nc - 1
						If contour(k) >= dmin And contour(k) < dmax
							For m = 4 To 0 Step -1

								If m >0
									h(m) = z(i + im(m - 1), j + jm(m - 1)) - contour(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 m

							;=================================================================
							;
							; 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
										x1 = xh(m1)
										y1 = yh(m1)
										x2 = xh(m2)
										y2 = yh(m2)

									Case 2
										; Line between vertices 2 and 3
										x1 = xh(m2)
										y1 = yh(m2)
										x2 = xh(m3)
										y2 = yh(m3)

									Case 3
										;Line between vertices 3 and 1
										x1 = xh(m3)
										y1 = yh(m3)
										x2 = xh(m1)
										y2 = yh(m1)

									Case 4
										; Line between vertex 1 and side 2-3
										x1 = xh(m1)
										y1 = yh(m1)
										x2 = (h(m3) * xh(m2) - h(m2) * xh(m3)) / (h(m3) - h(m2))
										y2 = (h(m3) * yh(m2) - h(m2) * yh(m3)) / (h(m3) - h(m2))

									Case 5
										; Line between vertex 2 and side 3-1
										x1 = xh(m2)
										y1 = yh(m2)
										x2 = (h(m1) * xh(m3) - h(m3) * xh(m1)) / (h(m1) - h(m3))
										y2 = (h(m1) * yh(m3) - h(m3) * yh(m1)) / (h(m1) - h(m3))

									Case 6
										; Line between vertex 3 and side 1-2
										x1 = xh(m3)
										y1 = yh(m3)
										x2 = (h(m2) * xh(m1) - h(m1) * xh(m2)) / (h(m2) - h(m1))
										y2 = (h(m2) * yh(m1) - h(m1) * yh(m2)) / (h(m2) - h(m1))

									Case 7
										; Line between sides 1-2 and 2-3
										x1 = (h(m2) * xh(m1) - h(m1) * xh(m2)) / (h(m2) - h(m1))
										y1 = (h(m2) * yh(m1) - h(m1) * yh(m2)) / (h(m2) - h(m1))
										x2 = (h(m3) * xh(m2) - h(m2) * xh(m3)) / (h(m3) - h(m2))
										y2 = (h(m3) * yh(m2) - h(m2) * yh(m3)) / (h(m3) - h(m2))

									Case 8
										; Line between sides 2-3 and 3-1
										x1 = (h(m3) * xh(m2) - h(m2) * xh(m3)) / (h(m3) - h(m2))
										y1 = (h(m3) * yh(m2) - h(m2) * yh(m3)) / (h(m3) - h(m2))
										x2 = (h(m1) * xh(m3) - h(m3) * xh(m1)) / (h(m1) - h(m3))
										y2 = (h(m1) * yh(m3) - h(m3) * yh(m1)) / (h(m1) - h(m3))

									Case 9
										; Line between sides 3-1 and 1-2
										x1 = (h(m1) * xh(m3) - h(m3) * xh(m1)) / (h(m1) - h(m3))
										y1 = (h(m1) * yh(m3) - h(m3) * yh(m1)) / (h(m1) - h(m3))
										x2 = (h(m2) * xh(m1) - h(m1) * xh(m2)) / (h(m2) - h(m1))
										y2 = (h(m2) * yh(m1) - h(m1) * yh(m2)) / (h(m2) - h(m1))

									EndSelect
									;--------------------------------------------------------------
									;this is where the program specific drawing routine comes in.
									;This specific command will work well for a properly dimensioned
									;vb picture box or vb form (where "Form1" is the name of the form)

									;-------------------------------------------------------------------

									For xxx=-1 To 1
										For yyy=-1 To 1
											LineXY(x1+xxx, y1+yyy, x2+xxx, y2+yyy, #White)
										Next yyy
									Next xxx
									LineXY(x1, y1, x2, y2, color(k))
									Debug Str(x1)+"/"+Str(y1)+" - "+Str(x2)+"/"+Str(y2)

									;-------------------------------------------------------------------
								EndIf
							Next m
						EndIf
					Next k
				EndIf
			Next i
		Next j
	EndIf
EndProcedure
Procedure Init()

	LoadFont (0,#FontName,#FontSize,#FontType)


	OpenWindow(0,0,0,#WX<<10/DesktopScaledX(1024),#WY<<10/DesktopScaledY(1024),"Check Dots",#PB_Window_SystemMenu|#PB_Window_ScreenCentered)
	CanvasGadget(0,0,0,#WX,#WY)

	CreateImage(0,#WX,#WY,32,#Red)
	StartDrawing(ImageOutput(0))
	DrawingMode(#PB_2DDrawing_AllChannels)
	Box(0,0,#WX,#WY,#ColBack)
	StopDrawing()

	StartVectorDrawing(ImageVectorOutput(0))
	VectorFont(FontID(0))
	MovePathCursor(#FX,#FY)
	AddPathText(#Text)
	VectorSourceColor(#Draw|#ColBorder)
	StrokePath(8,#PB_Path_Preserve|#PB_Path_RoundCorner)
	VectorSourceColor(#Draw|#ColFill)
	FillPath()

	w.f

	For j=0 To 3
		Select j
		Case 0
			VectorSourceColor($ff000000)
			l=3
			n=250
			w=1.55
		Case 1
			VectorSourceColor($20000000|#ColFill)
			l=5
			n=2000
			w=1.5
		Case 2
			VectorSourceColor($80000000|#ColBack)
			l=5
			n=2000
		Case 3
			n=10000
			l=25
		EndSelect

		For i=0 To n
			x=Random(#wx)
			y=Random(#wy)
			MovePathCursor(x,y)
			AddPathLine(x+Random(l),y+Random(l))
			If j=3
				VectorSourceColor($30000000|Random(#White))
			EndIf
			StrokePath(w)
		Next i
	Next j

	StopVectorDrawing()

EndProcedure
Procedure ScanNow()

	Protected x,y

	StartDrawing(ImageOutput(0))
	DrawingFont(FontID(0))

	While x<#WX
		y=0
		While y<#WY
			If Point(x,y)&#Color=0;<>#ColBack
				;Debug Hex(#ColBack)+" - "+Hex(Point(x,y)&#Color)
				snw(x)=y
				;Box(x,0,1,y,$F3E1CA)
				;Box(x-1,y,1,1,#Red)
				Scan(x)=y
				y=#WY
			Else
				y+1
			EndIf
		Wend
		x+1
	Wend

	StopDrawing()

EndProcedure
Procedure CalcNow()

	Protected x,y
	Protected ymin,ymax

	ymin=#WY
	StartDrawing(ImageOutput(0))
	DrawingFont(FontID(0))

	While x<#WX
		If Scan(x)
			If Scan(x)>ymax : ymax=Scan(x) : EndIf
			If Scan(x)<ymin : ymin=Scan(x) : EndIf
			Box(x-3,Scan(x)-3,7,7,#Black)
			Box(x-1,Scan(x)-1,3,3,#Red)
		EndIf
		x+1
	Wend

	LineXY(0,ymin,#WX,ymin,#Blue)
	LineXY(0,ymax,#WX,ymax,#Blue)

	y=ymax
	While y>=ymin
		x=0
		While x<#WX
			If Scan(x)=y
				;Box(x,Scan(x),1,1,#White)
			EndIf
			x+1
		Wend
		y-1
	Wend

	StopDrawing()


EndProcedure

Global Dim z.d(#wx,#wy)
Global Dim x.d(#wx)
Global Dim y.d(#wy)
Global Dim contour.d(1)
Global Dim color(1)
contour(1)=1
color(0)=$ff000000
color(1)=$ff0000ff

Procedure FindNow()


	x=0
	z=0
	While x<#WX
		If Scan(x)
			z+1
			z(x,Scan(x))=1
			x(z)=x
			y(z)=Scan(x)
		EndIf
		x+1
	Wend

	StartDrawing(ImageOutput(0))
	Conrec(z(),x(),y(),1,contour(),0,#wx,0,#wy,color())
	StopDrawing()

	If 0
		StartVectorDrawing(ImageVectorOutput(0))
		VectorSourceColor($80000000|#Magenta)

		MovePathCursor(54,260)
		AddPathLine(95,190)

		MovePathCursor(105,180)
		AddPathLine(180,150)

		MovePathCursor(185,150)
		AddPathLine(270,160)

		MovePathCursor(150,220)
		AddPathLine(230,210)

		StrokePath(20)

		StopVectorDrawing()
	EndIf


EndProcedure
Procedure Main()

	Protected i

	Init()
	ScanNow()
	CalcNow()
	FindNow()

	StartDrawing(CanvasOutput(0))
	DrawImage(ImageID(0),0,0)

	StopDrawing()

	Repeat
		Select WaitWindowEvent()
		Case #PB_Event_CloseWindow,#WM_CHAR
			End
		EndSelect
	ForEver

EndProcedure

Main()

[Michael Vogel]

Another twist - got no good results for Conrec - now I am scanning from each point (x) in the list a certain number of dots to the right (maxradius) and calculate the slope of a line from the point x to these dots. After that I check which slopes are seen multiple times and which is the most right dot with such a slope.

Actually the line from x to such a dot (drawn in magenta) is what I am searching for (but not in all cases, see below).

Before doing the next scan, I remove all dots from the list which are nearby the found line.

In the quick and dirty written source below (a perfect example for bad programming style) the procedure FindNow is doing this job - and it looks fine at the first moment, because I have limited the scanning part to a single point (For x=107 to 107), as soon this line will be changed the program does show a lot of unwanted lines as well.

Code: Select all

; Define

	RandomSeed(0)

	#Text="O"

	#FontName="Segoe UI"
	#FontSize=250
	#FontType=#PB_Font_Bold;|#PB_Font_Italic

	#WX=400
	#WY=600
	#FX=40
	#FY=0

	#ColBack=	$FFFFC0
	#ColBorder=	#Black
	#ColFill=		#Yellow

	#Draw=		$FF000000
	#Color=		$00FFFFFF

	Global Dim snw(#WX)

	Global Dim Scan(#WX)
	Global Dim Dir(#WX)
	Global Dim Clst(#WX)

; EndDefine
Procedure Init()

	LoadFont (0,#FontName,#FontSize,#FontType)


	OpenWindow(0,0,0,#WX<<10/DesktopScaledX(1024),#WY<<10/DesktopScaledY(1024),"Check Dots",#PB_Window_SystemMenu|#PB_Window_ScreenCentered)
	CanvasGadget(0,0,0,#WX,#WY)

	CreateImage(0,#WX,#WY,32,#Red)
	StartDrawing(ImageOutput(0))
	DrawingMode(#PB_2DDrawing_AllChannels)
	Box(0,0,#WX,#WY,#ColBack)
	StopDrawing()

	StartVectorDrawing(ImageVectorOutput(0))
	VectorFont(FontID(0))
	MovePathCursor(#FX,#FY)
	AddPathText(#Text)
	VectorSourceColor(#Draw|#ColBorder)
	StrokePath(8,#PB_Path_Preserve|#PB_Path_RoundCorner)
	VectorSourceColor(#Draw|#ColFill)
	FillPath()

	w.f

	For j=0 To 3
		Select j
		Case 0
			VectorSourceColor($ff000000)
			l=3
			n=250
			w=1.55
		Case 1
			VectorSourceColor($20000000|#ColFill)
			l=5
			n=2000
			w=1.5
		Case 2
			VectorSourceColor($80000000|#ColBack)
			l=5
			n=2000
		Case 3
			n=10000
			l=25
		EndSelect

		For i=0 To n
			x=Random(#wx)
			y=Random(#wy)
			MovePathCursor(x,y)
			AddPathLine(x+Random(l),y+Random(l))
			If j=3
				VectorSourceColor($30000000|Random(#White))
			EndIf
			StrokePath(w)
		Next i
	Next j

	StopVectorDrawing()

EndProcedure
Procedure ScanNow()

	Protected x,y

	StartDrawing(ImageOutput(0))
	DrawingFont(FontID(0))

	While x<#WX
		y=0
		While y<#WY
			If Point(x,y)&#Color=0;<>#ColBack
				;Debug Hex(#ColBack)+" - "+Hex(Point(x,y)&#Color)
				snw(x)=y
				;Box(x,0,1,y,$F3E1CA)
				;Box(x-1,y,1,1,#Red)
				Scan(x)=y
				y=#WY
			Else
				y+1
			EndIf
		Wend
		x+1
	Wend

	StopDrawing()

EndProcedure
Procedure CalcNow()

	Protected x,y
	Protected ymin,ymax

	ymin=#WY
	StartDrawing(ImageOutput(0))
	DrawingFont(FontID(0))

	While x<#WX
		If Scan(x)
			If Scan(x)>ymax : ymax=Scan(x) : EndIf
			If Scan(x)<ymin : ymin=Scan(x) : EndIf
			Box(x-3,Scan(x)-3,7,7,#Black)
			Box(x-1,Scan(x)-1,3,3,#Red)
		EndIf
		x+1
	Wend

	LineXY(0,ymin,#WX,ymin,#Blue)
	LineXY(0,ymax,#WX,ymax,#Blue)

	y=ymax
	While y>=ymin
		x=0
		While x<#WX
			If Scan(x)=y
				;Box(x,Scan(x),1,1,#White)
			EndIf
			x+1
		Wend
		y-1
	Wend

	StopDrawing()


EndProcedure

Procedure FatLine(x1,y1,x2,y2,width=1,color1=#Red,color2=#Black)

	Protected xxx,yyy


	For xxx=-width To width
		For yyy=-width To width
			If Abs(xxx)=width Or Abs(yyy)=width
				LineXY(x1+xxx,y1+yyy,x2+xxx,y2+yyy,color2)
			EndIf
		Next yyy
	Next xxx

	width-1
	For xxx=-width To width
		For yyy=-width To width
			LineXY(x1+xxx,y1+yyy,x2+xxx,y2+yyy,color1)
		Next yyy
	Next xxx

EndProcedure
Procedure InRange(value,min,max)

	If min>max
		Swap min,max
	EndIf

	ProcedureReturn Bool(value>=min And value<=max)

EndProcedure
Procedure FindNow()

	If 0
		StartVectorDrawing(ImageVectorOutput(0))
		VectorSourceColor($80000000|#Magenta)

		MovePathCursor(54,260)
		AddPathLine(95,190)

		MovePathCursor(105,180)
		AddPathLine(180,150)

		MovePathCursor(185,150)
		AddPathLine(270,160)

		MovePathCursor(150,220)
		AddPathLine(230,210)

		StrokePath(20)

		StopVectorDrawing()
	EndIf


	If 1
		StartDrawing(ImageOutput(0))

		maxslope=		100
		maxdist=		10
		slotsize=		maxdist/2
		maxradius=	100
		maxdelta=		5

		maxslots=		maxslope/slotsize

		d.i

		Protected Dim slp(maxradius)
		Protected Dim idx(maxradius)
		Protected Dim sct(maxradius)
		Protected Dim slt(maxslots*2)

		;For x=0 To #WX-1;					This is the correct loop to scan all dots
		;For x=105 To 130;					This limits the routine to scan only for some dots
		For x=107 To 107;					Just a demo which scans for a single dot to show how it should work
			If Scan(x)
				index=0
				x2=x+maxradius
				If x2>=#wx
					x2=#wx-1
				EndIf
				While x2>x
					If Scan(x2)
						d=(scan(x2)-Scan(x))*100/(x2-x)
						If Abs(d)<=maxslope
							slp(index)=d;					Steigung
							idx(index)=x2;					X-Position
							sct(index)=d/slotsize+maxslots;	Sektor
							index+1
							slt(d/slotsize+maxslots)+1;		Sektorzähler erhöhen
							LineXY(x,scan(x),x2,scan(x2),Random(#White))
						EndIf
					EndIf
					x2-1
				Wend
				If index
					; Höchste Zahl der Vektoren in zwei benachbarten Sektoren suchen...
					For i=0 To index-1
						;Debug Str(i)+": "+Str(slp(i))+"="+Str(sct(i))+", "+Str(idx(i))+": "+Str(x)+"|"+Str(Scan(x))+" - "+Str(idx(i))+"|"+Str(Scan(idx(i)))
					Next i
					Debug "---"
					max=0
					best=-1
					For i=0 To maxslots*2-1
						n=slt(i)+slt(i+1)
						If n>max
							max=n
							best=i
						EndIf
					Next i

					If best>=0 And max>5;	Minimum 5 dots in range
						Debug best
						Debug slp(best)

						i=0
						found=0
						While i<index
							d=slp(i)/slotsize+maxslots
							If found=0
								If d=best Or d=best+1
									Debug sct(i)
									;Circle(idx(i),Scan(idx(i))-30,3,#Black)
									;Circle(idx(i),Scan(idx(i))-30,2,#Magenta)
									Debug "LINE"
									FatLine(x,Scan(x),idx(i),Scan(idx(i)),6,#Magenta)
									found=idx(i)
									slope.d=(scan(found)-Scan(x))/(found-x)
									Debug slope
								EndIf
							Else
								dotx=idx(i)
								doty=(dotx-x)*slope+Scan(x)
								If Abs(doty-Scan(dotx))<6;	remove dots from list when not near ideal line
									FatLine(idx(i),Scan(idx(i)),idx(i),Scan(idx(i)),3,#White)
									Scan(idx(i))=0
									
								ElseIf d=best Or d=best+1;	(disabled) - remove dots wehn in sector
									Debug Str(x)+" - "+Str(idx(i))+" - "+Str(found)
									FatLine(idx(i),Scan(idx(i)),idx(i),Scan(idx(i)),2,#Green)
									;Scan(idx(i))=0

								ElseIf found And InRange(idx(i),x,found) And InRange(Scan(idx(i)),Scan(x),Scan(found)); (disabled)
									FatLine(idx(i),Scan(idx(i)),idx(i),Scan(idx(i)),2,#Blue)
									;Scan(idx(i))=0

								EndIf
							EndIf
							i+1
						Wend
					EndIf
				EndIf
			EndIf

		Next x
		StopDrawing()

	EndIf

EndProcedure
Procedure Main()

	Protected i

	Init()
	ScanNow()
	CalcNow()
	FindNow()

	StartDrawing(CanvasOutput(0))
	DrawImage(ImageID(0),0,0)

	StopDrawing()

	Repeat
		Select WaitWindowEvent()
		Case #PB_Event_CloseWindow,#WM_CHAR
			End
		EndSelect
	ForEver

EndProcedure

Main()