Peter Popov wrote:

> >I have posted another  image to p.b.i that shows a "height field"
> >above a plane with 6 points. This the lowest point in this height
> >field shows where the "best" centre point is for a circle that comes
> >close to the given points.
>
> Very illustrative, thanks!

You're welcome. Interesting problem btw.

> >The "height field" can be considered as an "error field".
> >The height above "ground" gives the least "error" generated if
> >placing the circle below.
>
> So basically you find the standard error of the distances between the
> points and the circle? Nifty :)

;)

First I  thought about doing linear regression analysis on the points
to get a vector for a tangent line of the circle. Then one could
construct
a normal vector to a tangent line of the circle. This vector could then
be used as a direction vector for a line to place the circle centre on.

But then I realized that it could be difficult to place the line and
thereafter to place the circle centre on it.

The linear regression gave birth to the idea of using the sum of
squares of  the errors as a measure of how suitable a certain
point would be for the circle centre.


> >The code generating this image also estimates which radius to use
> >for the circle. And the code can also be used with more or less
> >points.
>
> Interesting, how do you do that, by sampling?

Mean distance from the point in question to the other points.

(See code below.)

The problem with this code is that one has to know which region
to do the search. This could be frustrating if many of the points
lies close to a straight line. That would cause the region to
expand very much because the centre would approach infinity.

Therefore I'm working on macros to implement the Nelder-
Mead simplex method to locate the minima of the "error field".

Tor Olav
--
mailto:tor### [at] hotmailcom
http://www.crosswinds.net/~tok/tokrays.html

#version 3.1;

global_settings {
  #max_trace_level 1
  ambient_light 2
}

#include "colors.inc"


#macro SumSqdErrors(Point, PointArray, Radius)

  #local NrOfPoints = dimension_size(PointArray, 1);
  #local ErrorSum = 0;

  #local PtCnt = 0;
  #while(PtCnt < NrOfPoints)
    #local ErrorSum = ErrorSum +
           pow(vlength(Point - PointArray[PtCnt]) - Radius, 2);
    #local PtCnt = PtCnt + 1;
  #end // while

  ErrorSum

#end // macro SumSqdErrors


#macro MeanDist(Point, PointArray)

  #local NrOfPoints = dimension_size(PointArray, 1);
  #local DistSum = 0;

  #local PtCnt = 0;
  #while(PtCnt < NrOfPoints)
    #local ThisDist = vlength(Point - PointArray[PtCnt]);
    #local DistSum = DistSum + ThisDist;
    #local PtCnt = PtCnt + 1;
  #end // while

  (DistSum/NrOfPoints)

#end // macro MeanDist


#macro StoreErrors(PointArray, ErrorArray, RadiusArray, xMin, xMax, zMin,
zMax)

  #local xNr = dimension_size(ErrorArray, 1);
  #local zNr = dimension_size(ErrorArray, 2);

  #local dx = (xMax - xMin)/xNr;
  #local dz = (zMax - zMin)/zNr;

  #local xPos = xMin;
  #local xCnt = 0;
  #while(xCnt < xNr)
    #local zPos = zMin;
    #local zCnt = 0;
    #while(zCnt < zNr)
      #local Point = <xPos, 0, zPos>;
      #local BestRadius = MeanDist(Point, PointArray);
      #local MinimumError = SumSqdErrors(Point, PointArray, BestRadius);
      #declare RadiusArray[xCnt][zCnt] = BestRadius;
      #declare ErrorArray[xCnt][zCnt] = Point + y*MinimumError;
      #local zCnt = zCnt + 1;
      #local zPos = zPos + dz;
    #end // while
    #local xCnt = xCnt + 1;
    #local xPos = xPos + dx;
  #end // while

#end // macro StoreErrors


#macro LocateBest(ErrorArray, RadiusArray, xBest, zBest)

  #local xNr = dimension_size(ErrorArray, 1);
  #local zNr = dimension_size(ErrorArray, 2);
  #declare xBest = 0;
  #declare zBest = 0;
  #local LeastError = ErrorArray[0][0].y;

  #local xCnt = 0;
  #while(xCnt < xNr)
    #local zCnt = 0;
    #while(zCnt < zNr)
      #local CurrentError = ErrorArray[xCnt][zCnt].y;
      #if(CurrentError < LeastError)
        #local LeastError = CurrentError;
        #declare xBest = xCnt;
        #declare zBest = zCnt;
      #end // if
      #local zCnt = zCnt + 1;
    #end // while
    #local xCnt = xCnt + 1;
  #end // while

#end // macro LocateBest


#macro PlotErrors(ErrorArray, SphereRadius)

  #local xNr = dimension_size(ErrorArray, 1);
  #local zNr = dimension_size(ErrorArray, 2);

  union {
    #local xCnt = 0;
    #while(xCnt < xNr)
      #local zCnt = 0;
      #while(zCnt < zNr)
        sphere { ErrorArray[xCnt][zCnt], SphereRadius }
        #local zCnt = zCnt + 1;
      #end // while
      #local xCnt = xCnt + 1;
    #end // while
  }

#end // macro PlotErrors


#macro PtPlot(PointArray, SphereRadius)

  #local NrOfPoints = dimension_size(PointArray, 1);

  union {
    #local PtCnt = 0;
    #while(PtCnt < NrOfPoints)
      sphere { PointArray[PtCnt], SphereRadius }
      #local PtCnt = PtCnt + 1;
    #end // while
  }

#end // macro PtPlot


#declare SplinePts =
array[6] {
  <1.0, 0, 1.9>,
  <2.0, 0, 1.1>,
  <1.7, 0, 1.5>,
  <2.9, 0, 1.7>,
  <3.1, 0, 3.0>,
  <3.1, 0, 3.4>
}


#declare xCount = 80;
#declare zCount = 80;

#declare ErrArray = array[xCount][zCount]
#declare RadArray = array[xCount][zCount]

#declare xMin = -1;
#declare xMax =  5;
#declare zMin = -1;
#declare zMax =  5;

StoreErrors(SplinePts, ErrArray, RadArray, xMin, xMax, zMin, zMax)

#declare xBest = 0;
#declare zBest = 0;

LocateBest(ErrArray, RadArray, xBest, zBest)

#declare BestPoint = ErrArray[xBest][zBest];


object {
  PtPlot(SplinePts, 0.06)
  pigment { color Red*3 }
}

sphere {
  BestPoint, 0.06
  pigment { color Green }
}

torus {
  RadArray[xBest][zBest], 0.02
  translate (x + z)*BestPoint
  pigment { color White*2 }
}

object {
  PlotErrors(ErrArray, 0.03*1.2)
  pigment { color Yellow }
  scale <1, 1/3.5, 1>
}


camera {
//  location <-3, 0.7, 7>
//  location <4.5, 0.8, 9>
  location <-1.5, 4.5*1.3, 9>*0.85
//  location BestPoint + 7*y
  look_at BestPoint
}

background { color Blue/2 }

light_source { <10,  100, 10>, 1 }
light_source { <10, -100, 10>, 1 }

Post a reply to this message