Principal Component Analysis

[IdeasVacuum]

Has anyone used PCA to determine X,Y,Z direction vectors (axes of a Cartesian coordinate system) for a mesh model (points set) of arbitrary location/orientation in space?

[PB]

Aren't these the English forums?

[DarkDragon]

Sure, I did it once, I will look for it in my holidays which will begin next week.

[IdeasVacuum]

Hi Dark Dragon, that will be great thank you. Too late to save my hair though, it's nearly 100% gray already.........

[idle]

If you only need to find the axis of orientation maybe you could use Moments of Inertia

[IdeasVacuum]

Hello idle

The mesh models are complex, irregular shapes unfortunately.

Edit: Ah, the Moment of inertia tensor?

[idle]

I don't know if that would be a problem, I've never tried to use moments in 3D to obtain an objects orientation.

[IdeasVacuum]

Thus far, I have not got very far! What I have discovered is that if a I use a simple data set (cube), then I get a good value for the centroid, but for a complex data set, the centroid is never correct. The datum axes produced are parallel to the global axes, when they should not be parallel to any global plane as the data set is at a compound angle (always).

idle tried a different approach producing second moments (Many thanks idle). This didn't produce a local coordinate system either but I have learnt a lot about PB macros in the process.

[IdeasVacuum]

....tried a different approach, move the data set (vertex points) to global datum, rotate the points by n degrees, test the volume of their bounding box (BB in global coordinates), repeat rotates until smallest volume is found - a kind of 'virtual caliper' approach. The results are still not good enough unfortunately, even though rotating about all three axis and by ever smaller amounts.

[DarkDragon]

Hello,

I had to recode it, because I couldn't find the sources anymore. It can calculate the eigenvalues and with them it will give you the eigenvectors, which are the directional vectors of the PCA. The list of eigenvectors is already sorted decending by corresponding eigenvalue, so you just need to cut eigenvectors from the end to get a featurevector.

Code: Select all

; #############################################
; #                                           #
; # Author:   Daniel Brall                    #
; # Website:  http://www.bradan.eu/           #
; #                                           #
; # Thanks to PureFan, who helped me a lot!   #
; #                                           #
; #############################################

#PRECISION = 0.00000005 ; 0.00..005 is more exact in floating point than 0.000...001

;- Polynomialfunctions
;{
Structure Monom
  coeff.d
  exponent.i
EndStructure

Structure Polynomial
  List Monoms.Monom()
EndStructure

Procedure TidyUpPolynomial(*Poly.Polynomial)
  ForEach *Poly\Monoms()
    If *Poly\Monoms()\coeff = 0.0
      DeleteElement(*Poly\Monoms())
      ResetList(*Poly\Monoms())
    EndIf
  Next
EndProcedure

Procedure CopyPolynomial(*Dest.Polynomial, *Src.Polynomial)
  ClearList(*Dest\Monoms())

  ForEach *Src\Monoms()
    AddElement(*Dest\Monoms())
    *Dest\Monoms() = *Src\Monoms()
  Next

  TidyUpPolynomial(*Dest)
EndProcedure

Procedure SetPolynomialCoeff(*Poly.Polynomial, Exponent.i, Coeff.d)
  Protected Found.i

  Found = #False
  ForEach *Poly\Monoms()
    If *Poly\Monoms()\exponent = Exponent
      *Poly\Monoms()\coeff = Coeff
      Found                = #True
      Break
    EndIf
  Next

  If Not Found
    AddElement(*Poly\Monoms())
    *Poly\Monoms()\coeff    = Coeff
    *Poly\Monoms()\exponent = Exponent
  EndIf

  TidyUpPolynomial(*Poly)
EndProcedure

Procedure ClearPolynomial(*Poly.Polynomial)
  ClearStructure(*Poly, Polynomial)
  InitializeStructure(*Poly, Polynomial)
EndProcedure

Procedure AddPolynomials(*Dest.Polynomial, *Src.Polynomial)
  Protected Found.i

  ForEach *Src\Monoms()
    Found = #False

    ForEach *Dest\Monoms()
      If *Dest\Monoms()\exponent = *Src\Monoms()\exponent
        Found = #True
        *Dest\Monoms()\coeff + *Src\Monoms()\coeff
        Break
      EndIf
    Next

    If Not Found
      AddElement(*Dest\Monoms())
      *Dest\Monoms()\coeff    = *Src\Monoms()\coeff
      *Dest\Monoms()\exponent = *Src\Monoms()\exponent
    EndIf
  Next

  TidyUpPolynomial(*Dest)
EndProcedure

Procedure SubPolynomials(*Dest.Polynomial, *Src.Polynomial)
  Protected Found.i

  ForEach *Src\Monoms()
    Found = #False

    ForEach *Dest\Monoms()
      If *Dest\Monoms()\exponent = *Src\Monoms()\exponent
        Found = #True
        *Dest\Monoms()\coeff - *Src\Monoms()\coeff
        Break
      EndIf
    Next

    If Not Found
      AddElement(*Dest\Monoms())
      *Dest\Monoms()\coeff    = -1 * *Src\Monoms()\coeff
      *Dest\Monoms()\exponent = *Src\Monoms()\exponent
    EndIf
  Next

  TidyUpPolynomial(*Dest)
EndProcedure

Procedure MulPolynomials(*Dest.Polynomial, *Src.Polynomial)
  Protected ResultPolynomial.Polynomial
  Protected Found.i

  ForEach *Src\Monoms()
    ForEach *Dest\Monoms()
      Found = #False

      ForEach ResultPolynomial\Monoms()
        If ResultPolynomial\Monoms()\exponent = *Dest\Monoms()\exponent + *Src\Monoms()\exponent
          Found = #True
          ResultPolynomial\Monoms()\coeff  + *Dest\Monoms()\coeff * *Src\Monoms()\coeff
          Break
        EndIf
      Next

      If Not Found
        AddElement(ResultPolynomial\Monoms())
        ResultPolynomial\Monoms()\exponent = *Dest\Monoms()\exponent + *Src\Monoms()\exponent
        ResultPolynomial\Monoms()\coeff    = *Dest\Monoms()\coeff * *Src\Monoms()\coeff
      EndIf
    Next
  Next

  CopyPolynomial(*Dest, @ResultPolynomial)

  TidyUpPolynomial(*Dest)
EndProcedure

Procedure MulPolynomial(*Dest.Polynomial, Factor.d)
  ForEach *Dest\Monoms()
    *Dest\Monoms()\coeff * Factor
  Next
  TidyUpPolynomial(*Dest)
EndProcedure

Procedure DerivePolynomial(*Dest.Polynomial, *Src.Polynomial)
  ClearList(*Dest\Monoms())
  
  ForEach *Src\Monoms()
    If *Src\Monoms()\exponent <> 0
      AddElement(*Dest\Monoms())
      *Dest\Monoms()\coeff = *Src\Monoms()\exponent * *Src\Monoms()\coeff
      *Dest\Monoms()\exponent = *Src\Monoms()\exponent - 1
    EndIf
  Next

  TidyUpPolynomial(*Dest)
EndProcedure

Procedure.i PolynomialDegree(*Polynomial.Polynomial)
  Protected Result.i = 0

  ForEach *Polynomial\Monoms()
    If *Polynomial\Monoms()\exponent > Result
      Result = *Polynomial\Monoms()\exponent
    EndIf
  Next

  ProcedureReturn Result
EndProcedure

Procedure.d PolynomialValueAt(*Polynomial.Polynomial, x.d, Degree.i = -1)
  Protected Result.d

  Result = 0.0

  If Degree = -1
    ForEach *Polynomial\Monoms()
      Result + *Polynomial\Monoms()\coeff * Pow(x, *Polynomial\Monoms()\exponent)
    Next
  Else
    ForEach *Polynomial\Monoms()
      If *Polynomial\Monoms()\exponent = Degree
        Result + *Polynomial\Monoms()\coeff * Pow(x, *Polynomial\Monoms()\exponent)
      EndIf
    Next
  EndIf

  ProcedureReturn Result
EndProcedure

Procedure.d BisectionNull(*Polynomial.Polynomial, IntervalStart.d, IntervalEnd.d)
  Protected Middle.d
  Protected CompareValue.d
  Protected Value1.d
  Protected Value2.d
  
  Repeat
    Middle = IntervalStart + (IntervalEnd - IntervalStart) * 0.5
    
    Value1 = PolynomialValueAt(*Polynomial, IntervalStart)
    Value2 = PolynomialValueAt(*Polynomial, IntervalEnd)
    
    CompareValue = PolynomialValueAt(*Polynomial, Middle)
    
    If Value1 < Value2
      If CompareValue < 0.0
        IntervalStart = Middle
      ElseIf CompareValue > 0.0
        IntervalEnd = Middle
      Else
        Break
      EndIf
    Else
      If CompareValue < 0.0
        IntervalEnd   = Middle
      ElseIf CompareValue > 0.0
        IntervalStart = Middle
      Else
        Break
      EndIf
    EndIf
  Until Abs(IntervalEnd - IntervalStart) <= 0.00000000000005 Or CompareValue = 0.0
  
  ProcedureReturn Middle
EndProcedure

Procedure FindNulls(*Polynomial.Polynomial, List Points.d())
  Protected NewList ExtremePoints.d()
  Protected Derivation.Polynomial
  Protected IntervalStart.d, IntervalEnd.d
  Protected V1.d
  Protected V2.d
  Protected Null1.d, Null2.d
  Protected Degree.i
  Protected DegreeCoeff.d
  
  ClearList(Points())
  
  Degree      = PolynomialDegree(*Polynomial)
  DegreeCoeff = PolynomialValueAt(*Polynomial, 1.0, Degree)

  If DegreeCoeff <> 0.0
    MulPolynomial(*Polynomial, 1.0 / DegreeCoeff)
  EndIf
  
  V1 = PolynomialValueAt(*Polynomial, 1.0, Degree - 1)
  V2 = PolynomialValueAt(*Polynomial, 1.0, Degree - 2)
  
  Degree = PolynomialDegree(*Polynomial)
  If Degree > 2
    DegreeCoeff = 1.0 / Degree
    IntervalStart = DegreeCoeff * (-1.0 * V1 - (Degree - 1) * Sqr(V1 * V1 - (2 * Degree * V2) / (Degree - 1)))
    IntervalEnd   = DegreeCoeff * (-1.0 * V1 + (Degree - 1) * Sqr(V1 * V1 - (2 * Degree * V2) / (Degree - 1)))
    
    DerivePolynomial(@Derivation, *Polynomial)
    
    FindNulls(@Derivation, ExtremePoints())
    
    ; we also need the last interval point as extreme point
    AddElement(ExtremePoints())
    ExtremePoints() = IntervalEnd + 0.5
    
    SortStructuredList(ExtremePoints(), #PB_Sort_Ascending, 0, #PB_Sort_Double)
    
    V1 = -Infinity()
    ; remove all duplicates
    ForEach ExtremePoints()
      If ExtremePoints() <= V1
        DeleteElement(ExtremePoints())
      Else
        V1 = ExtremePoints()
      EndIf
    Next
    
    V2 = IntervalStart
    ForEach ExtremePoints()
      V1 = V2
      V2 = ExtremePoints()
      If (PolynomialValueAt(*Polynomial, V1) <= 0.0) XOr (PolynomialValueAt(*Polynomial, V2) < 0.0)
        ; now search for null inside this area
        AddElement(Points())
        Points() = BisectionNull(*Polynomial, V1, V2)
      EndIf
    Next
  ElseIf Degree = 2
    ; parable
    V1 * 0.5

    Null1 = -1.0 * V1 - Sqr(V1 * V1 - V2)
    Null2 = -1.0 * V1 + Sqr(V1 * V1 - V2)

    If Not IsNAN(Null1) And Not IsInfinity(Null1)
      AddElement(Points())
      Points() = Null1
    EndIf

    If Abs(Null1 - Null2) > 0.0 And Not IsNAN(Null2) And Not IsInfinity(Null2)
      AddElement(Points())
      Points() = Null2
    EndIf
  ElseIf Degree = 1
    Null1 = -1.0 * V2 / V1

    If Not IsNAN(Null1) And Not IsInfinity(Null1)
      AddElement(Points())
      Points() = Null1
    EndIf
  Else
    ; now everything is either null or something else.. what to do now?
  EndIf
  
  SortStructuredList(Points(), #PB_Sort_Ascending, 0, #PB_Sort_Double)
  
  V1 = -Infinity()
  ; remove all duplicates
  ForEach Points()
    If Points() <= V1
      DeleteElement(Points())
    Else
      V1 = Points()
    EndIf
  Next
EndProcedure

Procedure PolynomialFindNulls(*Polynomial.Polynomial, List Nulls.d())
  Protected PolynomialCopy.Polynomial

  ClearList(Nulls())

  CopyPolynomial(@PolynomialCopy, *Polynomial)

  FindNulls(@PolynomialCopy, Nulls())
EndProcedure
;}

;- Transformation for the linear equation system
;{
Procedure.i ClampI(Value.i, Min.i, Max.i)
  If Value < Min : Value = Min : EndIf
  If Value > Max : Value = Max : EndIf
  
  ProcedureReturn Value
EndProcedure

Procedure RowEchelonForm(Array Matrix.d(2), Start.i = 0)
  Protected x.i
  Protected y.i
  Protected z.i
  Protected Value.d
  Protected Row.i
  Protected StartX.i = ClampI(Start, 0, ArraySize(Matrix(), 1))
  Protected StartY.i = ClampI(Start, 0, ArraySize(Matrix(), 2))

  If ArraySize(Matrix(), 2) < ArraySize(Matrix(), 1)
    ProcedureReturn 0
  EndIf

  If StartX >= ArraySize(Matrix(), 1) Or StartY >= ArraySize(Matrix(), 2)
    If Abs(Matrix(StartX, StartY)) > #PRECISION
      For x = 0 To ArraySize(Matrix(), 1)
        Matrix(x, StartY) / Matrix(StartX, StartY)
      Next x
    EndIf
    ProcedureReturn 1
  EndIf

  For y = Start To ArraySize(Matrix(), 1)
    If Abs(Matrix(StartX, y)) > #PRECISION
      Row = y
      Break
    EndIf
  Next y

  If Row > 0
    ; swap this row with the first one
    For x = 0 To ArraySize(Matrix(), 1)
      Value             = Matrix(x, StartY)
      Matrix(x, StartY) = Matrix(x, Row)
      Matrix(x, Row)    = Value
    Next x
  EndIf

  If Abs(Matrix(StartX, StartY)) > #PRECISION
    For y = Start + 1 To ArraySize(Matrix(), 2)
      If Abs(Matrix(StartX, y)) > #PRECISION

        Value = Matrix(StartX, y) / Matrix(StartX, StartY)

        For z = 0 To ArraySize(Matrix(), 1)
          Matrix(z, y) - Value * Matrix(z, StartY)
        Next z

      EndIf
    Next y
  EndIf

  If RowEchelonForm(Matrix(), StartY + 1)
    If Abs(Matrix(StartX + 1, StartY + 1)) > #PRECISION
      For y = Start To 0 Step - 1
        Value = Matrix(StartX + 1, y) / Matrix(StartY + 1, StartY + 1)

        For z = 0 To ArraySize(Matrix(), 1)
          Matrix(z, y) - Value * Matrix(z, StartY + 1)
        Next z
      Next y
    EndIf

    Value = Matrix(StartX, StartY)
    If Abs(Value) > #PRECISION
      Value = 1.0 / Value
      For z = 0 To ArraySize(Matrix(), 1)
        Matrix(z, StartY) * Value
      Next z
    EndIf

    ProcedureReturn 1
  EndIf
EndProcedure
;}

;- Statistical methods including eigen-values and eigen-vectors
;{
Structure Vector
  Array v.d(1)
EndStructure

Procedure.d Mean(*Buffer.Double, BufferDimensions.i, BufferSize.i, UsedDimension.i = 0)
  Protected Result.d = 0.0
  Protected Count.i
  Protected StepSize.i = BufferDimensions * SizeOf(Double)
  Protected *BufferEnd = *Buffer + BufferSize
  Protected *Cursor.Double

  Count = (*BufferEnd - *Buffer) / (BufferDimensions * SizeOf(Double))

  If UsedDimension >= BufferDimensions Or Count <= 0
    ProcedureReturn 0.0
  EndIf

  *Cursor = *Buffer + UsedDimension * SizeOf(Double)
  While *Cursor < *BufferEnd
    Result + *Cursor\d

    *Cursor + StepSize
  Wend

  ProcedureReturn Result / Count
EndProcedure

Procedure.d Variance(*Buffer.Double, BufferDimensions.i, BufferSize.i, Mean.d, UsedDimension.i = 0)
  Protected Result.d = 0.0
  Protected Value.d
  Protected Count.i
  Protected StepSize.i = BufferDimensions * SizeOf(Double)
  Protected *BufferEnd = *Buffer + BufferSize
  Protected *Cursor.Double

  Count = (*BufferEnd - *Buffer) / (BufferDimensions * SizeOf(Double))

  If UsedDimension >= BufferDimensions Or Count <= 0
    ProcedureReturn 0.0
  EndIf

  If Count <= 15
    ; our sample is too small
    ; at a minimum amount we will use n instead of n - 1
    Count + 1
  EndIf

  *Cursor = *Buffer + UsedDimension * SizeOf(Double)
  While *Cursor < *BufferEnd
    Value = (*Cursor\d - Mean)
    Result + Value * Value

    *Cursor + StepSize
  Wend

  ProcedureReturn Result / (Count - 1)
EndProcedure

Procedure.d StdDev(Variance.d)
  ProcedureReturn Sqr(Variance)
EndProcedure

Procedure.d CoVariance(*Buffer.Double, BufferDimensions.i, BufferSize.i, MeanX.d, MeanY.d, UsedDimensionX.i = 0, UsedDimensionY.i = 1)
  Protected Result.d = 0.0
  Protected Count.i
  Protected StepSize.i = BufferDimensions * SizeOf(Double)
  Protected *BufferEnd = *Buffer + BufferSize
  Protected *CursorX.Double
  Protected *CursorY.Double

  Count = (*BufferEnd - *Buffer) / (BufferDimensions * SizeOf(Double))

  If UsedDimensionX >= BufferDimensions Or UsedDimensionY >= BufferDimensions Or Count <= 0
    ProcedureReturn 0.0
  EndIf

  If Count <= 15
    ; our sample is too small
    ; at a minimum amount we will use n instead of n - 1
    Count + 1
  EndIf

  *CursorX = *Buffer + UsedDimensionX * SizeOf(Double)
  *CursorY = *Buffer + UsedDimensionY * SizeOf(Double)
  While *CursorX < *BufferEnd And *CursorY < *BufferEnd
    Result + (*CursorX\d - MeanX) * (*CursorY\d - MeanY)

    *CursorX + StepSize
    *CursorY + StepSize
  Wend

  ProcedureReturn Result / (Count - 1)
EndProcedure

Procedure CoVarianceMatrix(*Buffer.Double, BufferDimensions.i, BufferSize.i, Array ResultMatrix.d(2))
  Protected Dim MeanValues.d(BufferDimensions - 1)
  Protected X.i, Y.i

  If ArraySize(ResultMatrix(), 1) <> ArraySize(MeanValues()) Or ArraySize(ResultMatrix(), 2) <> ArraySize(MeanValues())
    ProcedureReturn 0
  EndIf

  For X = 0 To ArraySize(MeanValues())
    MeanValues(X) = Mean(*Buffer, BufferDimensions, BufferSize, X)
  Next X

  For X = 0 To ArraySize(MeanValues())
    For Y = 0 To ArraySize(MeanValues())
      ResultMatrix(X, Y) = CoVariance(*Buffer, BufferDimensions, BufferSize, MeanValues(X), MeanValues(Y), X, Y)
    Next Y
  Next X

  ProcedureReturn 1
EndProcedure

Procedure PolyEigenDet(*Polynomial.Polynomial, Array Matrix.Polynomial(2))
  Protected i.i
  Protected j.i
  Protected k.i
  Protected Sign.d = 1.0
  Protected Polynomial2.Polynomial
  
  ; it should always be quadratical and *Polynom should be initialized
  If ArraySize(Matrix(), 2) <> ArraySize(Matrix(), 1) Or *Polynomial = #Null
    ProcedureReturn 0
  EndIf
  
  ClearPolynomial(Polynomial2)

  ; trivials
  If ArraySize(Matrix(), 1) = 0
    ProcedureReturn CopyPolynomial(*Polynomial, @Matrix(0, 0))
  ElseIf ArraySize(Matrix(), 1) = 1
    ; as this is just a copy we can directly act on the matrix
    MulPolynomials(@Matrix(0, 0), @Matrix(1, 1))
    MulPolynomials(@Matrix(0, 1), @Matrix(1, 0))
    SubPolynomials(@Matrix(0, 0), @Matrix(0, 1))
    
    ProcedureReturn CopyPolynomial(*Polynomial, @Matrix(0, 0))
  EndIf

  ; get the submatrix
  Protected Dim SubMatrix.Polynomial(ArraySize(Matrix(), 1) - 1, ArraySize(Matrix(), 1) - 1)
  
  For i = 0 To ArraySize(Matrix(), 1)
    
    ; clear all polynoms
    For j = 0 To ArraySize(SubMatrix(), 1)
      For k = 0 To ArraySize(SubMatrix(), 2)
        ClearPolynomial(@SubMatrix(j, k))
      Next k
    Next j
    
    ; copy the matrix
    For j = 0 To ArraySize(Matrix(), 1) - 1
      For k = 0 To i - 1
        CopyPolynomial(@SubMatrix(j, k), @Matrix(1 + j, k))
      Next
      For k = i To ArraySize(Matrix(), 2) - 1
        CopyPolynomial(@SubMatrix(j, k), @Matrix(1 + j, k + 1))
      Next
    Next
    
    PolyEigenDet(@Polynomial2, @SubMatrix())
    
    MulPolynomials(@Polynomial2, @Matrix(0, i))
    MulPolynomial(@Polynomial2, Sign)
    
    AddPolynomials(*Polynomial, @Polynomial2)
    Sign * -1.0 ; swap the sign
  Next
EndProcedure

Procedure GetEigenValues(List Values.d(), Array Matrix.d(2))
  Protected Polynomial.Polynomial
  Protected x.i
  Protected y.i
  
  ClearPolynomial(@Polynomial)
  
  ; it should always be quadratical
  If ArraySize(Matrix(), 2) <> ArraySize(Matrix(), 1)
    ProcedureReturn 0
  EndIf
  
  Protected Dim PolyMatrix.Polynomial(ArraySize(Matrix(), 1), ArraySize(Matrix(), 2))
  
  For x = 0 To ArraySize(PolyMatrix(), 1)
    For y = 0 To ArraySize(PolyMatrix(), 2)
      ClearPolynomial(@PolyMatrix(x, y))
      SetPolynomialCoeff(@PolyMatrix(x, y), 0, Matrix(x, y))
      
      If x = y
        SetPolynomialCoeff(@PolyMatrix(x, y), 1, -1.0)
      EndIf
    Next y
  Next x
  
  PolyEigenDet(@Polynomial, PolyMatrix())

  ; now calculate the nulls
  PolynomialFindNulls(@Polynomial, Values())
  
  ProcedureReturn 1
EndProcedure

Procedure.d GetEigenVectors(List Vectors.Vector(), List Values.d(), Array Matrix.d(2))
  Protected x.i
  Protected y.i
  Protected z.i
  Protected DiagSize.i
  Protected Found.i
  Protected Value.d
  Protected Dim ValueMatrix.d(ArraySize(Matrix(), 1), ArraySize(Matrix(), 2))
  
  ClearList(Vectors())
  
  SortStructuredList(Values(), #PB_Sort_Descending, 0, #PB_Sort_Double)
  
  ForEach Values()
    ; generate the matrix for this value
    For x = 0 To ArraySize(Matrix(), 1)
      For y = 0 To ArraySize(Matrix(), 2)
        ValueMatrix(x, y) = Matrix(x, y)
        If x = y
          ValueMatrix(x, y) - Values()
        EndIf
      Next y
    Next x
    
    If RowEchelonForm(ValueMatrix())
      DiagSize = ArraySize(ValueMatrix(), 1)
      If DiagSize > ArraySize(ValueMatrix(), 2)
        DiagSize = ArraySize(ValueMatrix(), 2)
      EndIf
      
      For x = 0 To DiagSize
        If Abs(ValueMatrix(x, x)) <= 0.00000005
          ; now move all lines down
          For y = ArraySize(ValueMatrix(), 2) To x + 1 Step -1
            For z = 0 To ArraySize(ValueMatrix(), 1)
              ValueMatrix(z, y) = ValueMatrix(z, y - 1)
            Next z
          Next y
          
          ; and set the value of this point to -1.0
          ValueMatrix(x, x) = -1.0
          
          ; check if this vector is already inside the list of vectors
          Found = #False
          ForEach Vectors()
            Found = #True
            For y = 0 To x
              If Abs(Vectors()\v(y) - ValueMatrix(x, y)) > 0.00000005
                Found = #False
                Break
              EndIf
            Next y
            
            If Found
              Break
            EndIf
          Next
          
          If Not Found
            ; copy this column into our list of vectors
            LastElement(Vectors())
            AddElement(Vectors())
            
            Dim Vectors()\v.d(ArraySize(ValueMatrix(), 2))
            
            For y = 0 To x
              Vectors()\v(y) = ValueMatrix(x, y)
            Next y
          EndIf
        EndIf
      Next x
    EndIf
  Next
  
  ; normalize each vector
  ForEach Vectors()
    Value = 0.0
    For y = 0 To ArraySize(Matrix(), 2)
      Value + Vectors()\v(y) * Vectors()\v(y)
    Next y
    If Value <> 0.0
      Value = 1.0 / Sqr(Value)
      For y = 0 To ArraySize(Matrix(), 2)
        Vectors()\v(y) * Value
      Next y
    EndIf
  Next
  
  ProcedureReturn 1
EndProcedure
;}

;- Example:
;{
Define k.i
Define i.i
Define Text.s

Dim TestData.d(149, 2)

; values taken from the iris dataset included with the R programming language:
TestData(  0, 0) =   5.1 : TestData(  0, 1) =   3.5 : TestData(  0, 2) =   1.4
TestData(  1, 0) =   4.9 : TestData(  1, 1) =   3.0 : TestData(  1, 2) =   1.4
TestData(  2, 0) =   4.7 : TestData(  2, 1) =   3.2 : TestData(  2, 2) =   1.3
TestData(  3, 0) =   4.6 : TestData(  3, 1) =   3.1 : TestData(  3, 2) =   1.5
TestData(  4, 0) =   5.0 : TestData(  4, 1) =   3.6 : TestData(  4, 2) =   1.4
TestData(  5, 0) =   5.4 : TestData(  5, 1) =   3.9 : TestData(  5, 2) =   1.7
TestData(  6, 0) =   4.6 : TestData(  6, 1) =   3.4 : TestData(  6, 2) =   1.4
TestData(  7, 0) =   5.0 : TestData(  7, 1) =   3.4 : TestData(  7, 2) =   1.5
TestData(  8, 0) =   4.4 : TestData(  8, 1) =   2.9 : TestData(  8, 2) =   1.4
TestData(  9, 0) =   4.9 : TestData(  9, 1) =   3.1 : TestData(  9, 2) =   1.5
TestData( 10, 0) =   5.4 : TestData( 10, 1) =   3.7 : TestData( 10, 2) =   1.5
TestData( 11, 0) =   4.8 : TestData( 11, 1) =   3.4 : TestData( 11, 2) =   1.6
TestData( 12, 0) =   4.8 : TestData( 12, 1) =   3.0 : TestData( 12, 2) =   1.4
TestData( 13, 0) =   4.3 : TestData( 13, 1) =   3.0 : TestData( 13, 2) =   1.1
TestData( 14, 0) =   5.8 : TestData( 14, 1) =   4.0 : TestData( 14, 2) =   1.2
TestData( 15, 0) =   5.7 : TestData( 15, 1) =   4.4 : TestData( 15, 2) =   1.5
TestData( 16, 0) =   5.4 : TestData( 16, 1) =   3.9 : TestData( 16, 2) =   1.3
TestData( 17, 0) =   5.1 : TestData( 17, 1) =   3.5 : TestData( 17, 2) =   1.4
TestData( 18, 0) =   5.7 : TestData( 18, 1) =   3.8 : TestData( 18, 2) =   1.7
TestData( 19, 0) =   5.1 : TestData( 19, 1) =   3.8 : TestData( 19, 2) =   1.5
TestData( 20, 0) =   5.4 : TestData( 20, 1) =   3.4 : TestData( 20, 2) =   1.7
TestData( 21, 0) =   5.1 : TestData( 21, 1) =   3.7 : TestData( 21, 2) =   1.5
TestData( 22, 0) =   4.6 : TestData( 22, 1) =   3.6 : TestData( 22, 2) =   1.0
TestData( 23, 0) =   5.1 : TestData( 23, 1) =   3.3 : TestData( 23, 2) =   1.7
TestData( 24, 0) =   4.8 : TestData( 24, 1) =   3.4 : TestData( 24, 2) =   1.9
TestData( 25, 0) =   5.0 : TestData( 25, 1) =   3.0 : TestData( 25, 2) =   1.6
TestData( 26, 0) =   5.0 : TestData( 26, 1) =   3.4 : TestData( 26, 2) =   1.6
TestData( 27, 0) =   5.2 : TestData( 27, 1) =   3.5 : TestData( 27, 2) =   1.5
TestData( 28, 0) =   5.2 : TestData( 28, 1) =   3.4 : TestData( 28, 2) =   1.4
TestData( 29, 0) =   4.7 : TestData( 29, 1) =   3.2 : TestData( 29, 2) =   1.6
TestData( 30, 0) =   4.8 : TestData( 30, 1) =   3.1 : TestData( 30, 2) =   1.6
TestData( 31, 0) =   5.4 : TestData( 31, 1) =   3.4 : TestData( 31, 2) =   1.5
TestData( 32, 0) =   5.2 : TestData( 32, 1) =   4.1 : TestData( 32, 2) =   1.5
TestData( 33, 0) =   5.5 : TestData( 33, 1) =   4.2 : TestData( 33, 2) =   1.4
TestData( 34, 0) =   4.9 : TestData( 34, 1) =   3.1 : TestData( 34, 2) =   1.5
TestData( 35, 0) =   5.0 : TestData( 35, 1) =   3.2 : TestData( 35, 2) =   1.2
TestData( 36, 0) =   5.5 : TestData( 36, 1) =   3.5 : TestData( 36, 2) =   1.3
TestData( 37, 0) =   4.9 : TestData( 37, 1) =   3.6 : TestData( 37, 2) =   1.4
TestData( 38, 0) =   4.4 : TestData( 38, 1) =   3.0 : TestData( 38, 2) =   1.3
TestData( 39, 0) =   5.1 : TestData( 39, 1) =   3.4 : TestData( 39, 2) =   1.5
TestData( 40, 0) =   5.0 : TestData( 40, 1) =   3.5 : TestData( 40, 2) =   1.3
TestData( 41, 0) =   4.5 : TestData( 41, 1) =   2.3 : TestData( 41, 2) =   1.3
TestData( 42, 0) =   4.4 : TestData( 42, 1) =   3.2 : TestData( 42, 2) =   1.3
TestData( 43, 0) =   5.0 : TestData( 43, 1) =   3.5 : TestData( 43, 2) =   1.6
TestData( 44, 0) =   5.1 : TestData( 44, 1) =   3.8 : TestData( 44, 2) =   1.9
TestData( 45, 0) =   4.8 : TestData( 45, 1) =   3.0 : TestData( 45, 2) =   1.4
TestData( 46, 0) =   5.1 : TestData( 46, 1) =   3.8 : TestData( 46, 2) =   1.6
TestData( 47, 0) =   4.6 : TestData( 47, 1) =   3.2 : TestData( 47, 2) =   1.4
TestData( 48, 0) =   5.3 : TestData( 48, 1) =   3.7 : TestData( 48, 2) =   1.5
TestData( 49, 0) =   5.0 : TestData( 49, 1) =   3.3 : TestData( 49, 2) =   1.4
TestData( 50, 0) =   7.0 : TestData( 50, 1) =   3.2 : TestData( 50, 2) =   4.7
TestData( 51, 0) =   6.4 : TestData( 51, 1) =   3.2 : TestData( 51, 2) =   4.5
TestData( 52, 0) =   6.9 : TestData( 52, 1) =   3.1 : TestData( 52, 2) =   4.9
TestData( 53, 0) =   5.5 : TestData( 53, 1) =   2.3 : TestData( 53, 2) =   4.0
TestData( 54, 0) =   6.5 : TestData( 54, 1) =   2.8 : TestData( 54, 2) =   4.6
TestData( 55, 0) =   5.7 : TestData( 55, 1) =   2.8 : TestData( 55, 2) =   4.5
TestData( 56, 0) =   6.3 : TestData( 56, 1) =   3.3 : TestData( 56, 2) =   4.7
TestData( 57, 0) =   4.9 : TestData( 57, 1) =   2.4 : TestData( 57, 2) =   3.3
TestData( 58, 0) =   6.6 : TestData( 58, 1) =   2.9 : TestData( 58, 2) =   4.6
TestData( 59, 0) =   5.2 : TestData( 59, 1) =   2.7 : TestData( 59, 2) =   3.9
TestData( 60, 0) =   5.0 : TestData( 60, 1) =   2.0 : TestData( 60, 2) =   3.5
TestData( 61, 0) =   5.9 : TestData( 61, 1) =   3.0 : TestData( 61, 2) =   4.2
TestData( 62, 0) =   6.0 : TestData( 62, 1) =   2.2 : TestData( 62, 2) =   4.0
TestData( 63, 0) =   6.1 : TestData( 63, 1) =   2.9 : TestData( 63, 2) =   4.7
TestData( 64, 0) =   5.6 : TestData( 64, 1) =   2.9 : TestData( 64, 2) =   3.6
TestData( 65, 0) =   6.7 : TestData( 65, 1) =   3.1 : TestData( 65, 2) =   4.4
TestData( 66, 0) =   5.6 : TestData( 66, 1) =   3.0 : TestData( 66, 2) =   4.5
TestData( 67, 0) =   5.8 : TestData( 67, 1) =   2.7 : TestData( 67, 2) =   4.1
TestData( 68, 0) =   6.2 : TestData( 68, 1) =   2.2 : TestData( 68, 2) =   4.5
TestData( 69, 0) =   5.6 : TestData( 69, 1) =   2.5 : TestData( 69, 2) =   3.9
TestData( 70, 0) =   5.9 : TestData( 70, 1) =   3.2 : TestData( 70, 2) =   4.8
TestData( 71, 0) =   6.1 : TestData( 71, 1) =   2.8 : TestData( 71, 2) =   4.0
TestData( 72, 0) =   6.3 : TestData( 72, 1) =   2.5 : TestData( 72, 2) =   4.9
TestData( 73, 0) =   6.1 : TestData( 73, 1) =   2.8 : TestData( 73, 2) =   4.7
TestData( 74, 0) =   6.4 : TestData( 74, 1) =   2.9 : TestData( 74, 2) =   4.3
TestData( 75, 0) =   6.6 : TestData( 75, 1) =   3.0 : TestData( 75, 2) =   4.4
TestData( 76, 0) =   6.8 : TestData( 76, 1) =   2.8 : TestData( 76, 2) =   4.8
TestData( 77, 0) =   6.7 : TestData( 77, 1) =   3.0 : TestData( 77, 2) =   5.0
TestData( 78, 0) =   6.0 : TestData( 78, 1) =   2.9 : TestData( 78, 2) =   4.5
TestData( 79, 0) =   5.7 : TestData( 79, 1) =   2.6 : TestData( 79, 2) =   3.5
TestData( 80, 0) =   5.5 : TestData( 80, 1) =   2.4 : TestData( 80, 2) =   3.8
TestData( 81, 0) =   5.5 : TestData( 81, 1) =   2.4 : TestData( 81, 2) =   3.7
TestData( 82, 0) =   5.8 : TestData( 82, 1) =   2.7 : TestData( 82, 2) =   3.9
TestData( 83, 0) =   6.0 : TestData( 83, 1) =   2.7 : TestData( 83, 2) =   5.1
TestData( 84, 0) =   5.4 : TestData( 84, 1) =   3.0 : TestData( 84, 2) =   4.5
TestData( 85, 0) =   6.0 : TestData( 85, 1) =   3.4 : TestData( 85, 2) =   4.5
TestData( 86, 0) =   6.7 : TestData( 86, 1) =   3.1 : TestData( 86, 2) =   4.7
TestData( 87, 0) =   6.3 : TestData( 87, 1) =   2.3 : TestData( 87, 2) =   4.4
TestData( 88, 0) =   5.6 : TestData( 88, 1) =   3.0 : TestData( 88, 2) =   4.1
TestData( 89, 0) =   5.5 : TestData( 89, 1) =   2.5 : TestData( 89, 2) =   4.0
TestData( 90, 0) =   5.5 : TestData( 90, 1) =   2.6 : TestData( 90, 2) =   4.4
TestData( 91, 0) =   6.1 : TestData( 91, 1) =   3.0 : TestData( 91, 2) =   4.6
TestData( 92, 0) =   5.8 : TestData( 92, 1) =   2.6 : TestData( 92, 2) =   4.0
TestData( 93, 0) =   5.0 : TestData( 93, 1) =   2.3 : TestData( 93, 2) =   3.3
TestData( 94, 0) =   5.6 : TestData( 94, 1) =   2.7 : TestData( 94, 2) =   4.2
TestData( 95, 0) =   5.7 : TestData( 95, 1) =   3.0 : TestData( 95, 2) =   4.2
TestData( 96, 0) =   5.7 : TestData( 96, 1) =   2.9 : TestData( 96, 2) =   4.2
TestData( 97, 0) =   6.2 : TestData( 97, 1) =   2.9 : TestData( 97, 2) =   4.3
TestData( 98, 0) =   5.1 : TestData( 98, 1) =   2.5 : TestData( 98, 2) =   3.0
TestData( 99, 0) =   5.7 : TestData( 99, 1) =   2.8 : TestData( 99, 2) =   4.1
TestData(100, 0) =   6.3 : TestData(100, 1) =   3.3 : TestData(100, 2) =   6.0
TestData(101, 0) =   5.8 : TestData(101, 1) =   2.7 : TestData(101, 2) =   5.1
TestData(102, 0) =   7.1 : TestData(102, 1) =   3.0 : TestData(102, 2) =   5.9
TestData(103, 0) =   6.3 : TestData(103, 1) =   2.9 : TestData(103, 2) =   5.6
TestData(104, 0) =   6.5 : TestData(104, 1) =   3.0 : TestData(104, 2) =   5.8
TestData(105, 0) =   7.6 : TestData(105, 1) =   3.0 : TestData(105, 2) =   6.6
TestData(106, 0) =   4.9 : TestData(106, 1) =   2.5 : TestData(106, 2) =   4.5
TestData(107, 0) =   7.3 : TestData(107, 1) =   2.9 : TestData(107, 2) =   6.3
TestData(108, 0) =   6.7 : TestData(108, 1) =   2.5 : TestData(108, 2) =   5.8
TestData(109, 0) =   7.2 : TestData(109, 1) =   3.6 : TestData(109, 2) =   6.1
TestData(110, 0) =   6.5 : TestData(110, 1) =   3.2 : TestData(110, 2) =   5.1
TestData(111, 0) =   6.4 : TestData(111, 1) =   2.7 : TestData(111, 2) =   5.3
TestData(112, 0) =   6.8 : TestData(112, 1) =   3.0 : TestData(112, 2) =   5.5
TestData(113, 0) =   5.7 : TestData(113, 1) =   2.5 : TestData(113, 2) =   5.0
TestData(114, 0) =   5.8 : TestData(114, 1) =   2.8 : TestData(114, 2) =   5.1
TestData(115, 0) =   6.4 : TestData(115, 1) =   3.2 : TestData(115, 2) =   5.3
TestData(116, 0) =   6.5 : TestData(116, 1) =   3.0 : TestData(116, 2) =   5.5
TestData(117, 0) =   7.7 : TestData(117, 1) =   3.8 : TestData(117, 2) =   6.7
TestData(118, 0) =   7.7 : TestData(118, 1) =   2.6 : TestData(118, 2) =   6.9
TestData(119, 0) =   6.0 : TestData(119, 1) =   2.2 : TestData(119, 2) =   5.0
TestData(120, 0) =   6.9 : TestData(120, 1) =   3.2 : TestData(120, 2) =   5.7
TestData(121, 0) =   5.6 : TestData(121, 1) =   2.8 : TestData(121, 2) =   4.9
TestData(122, 0) =   7.7 : TestData(122, 1) =   2.8 : TestData(122, 2) =   6.7
TestData(123, 0) =   6.3 : TestData(123, 1) =   2.7 : TestData(123, 2) =   4.9
TestData(124, 0) =   6.7 : TestData(124, 1) =   3.3 : TestData(124, 2) =   5.7
TestData(125, 0) =   7.2 : TestData(125, 1) =   3.2 : TestData(125, 2) =   6.0
TestData(126, 0) =   6.2 : TestData(126, 1) =   2.8 : TestData(126, 2) =   4.8
TestData(127, 0) =   6.1 : TestData(127, 1) =   3.0 : TestData(127, 2) =   4.9
TestData(128, 0) =   6.4 : TestData(128, 1) =   2.8 : TestData(128, 2) =   5.6
TestData(129, 0) =   7.2 : TestData(129, 1) =   3.0 : TestData(129, 2) =   5.8
TestData(130, 0) =   7.4 : TestData(130, 1) =   2.8 : TestData(130, 2) =   6.1
TestData(131, 0) =   7.9 : TestData(131, 1) =   3.8 : TestData(131, 2) =   6.4
TestData(132, 0) =   6.4 : TestData(132, 1) =   2.8 : TestData(132, 2) =   5.6
TestData(133, 0) =   6.3 : TestData(133, 1) =   2.8 : TestData(133, 2) =   5.1
TestData(134, 0) =   6.1 : TestData(134, 1) =   2.6 : TestData(134, 2) =   5.6
TestData(135, 0) =   7.7 : TestData(135, 1) =   3.0 : TestData(135, 2) =   6.1
TestData(136, 0) =   6.3 : TestData(136, 1) =   3.4 : TestData(136, 2) =   5.6
TestData(137, 0) =   6.4 : TestData(137, 1) =   3.1 : TestData(137, 2) =   5.5
TestData(138, 0) =   6.0 : TestData(138, 1) =   3.0 : TestData(138, 2) =   4.8
TestData(139, 0) =   6.9 : TestData(139, 1) =   3.1 : TestData(139, 2) =   5.4
TestData(140, 0) =   6.7 : TestData(140, 1) =   3.1 : TestData(140, 2) =   5.6
TestData(141, 0) =   6.9 : TestData(141, 1) =   3.1 : TestData(141, 2) =   5.1
TestData(142, 0) =   5.8 : TestData(142, 1) =   2.7 : TestData(142, 2) =   5.1
TestData(143, 0) =   6.8 : TestData(143, 1) =   3.2 : TestData(143, 2) =   5.9
TestData(144, 0) =   6.7 : TestData(144, 1) =   3.3 : TestData(144, 2) =   5.7
TestData(145, 0) =   6.7 : TestData(145, 1) =   3.0 : TestData(145, 2) =   5.2
TestData(146, 0) =   6.3 : TestData(146, 1) =   2.5 : TestData(146, 2) =   5.0
TestData(147, 0) =   6.5 : TestData(147, 1) =   3.0 : TestData(147, 2) =   5.2
TestData(148, 0) =   6.2 : TestData(148, 1) =   3.4 : TestData(148, 2) =   5.4
TestData(149, 0) =   5.9 : TestData(149, 1) =   3.0 : TestData(149, 2) =   5.1

Dim MeanValue.d(ArraySize(TestData(), 2))
Dim VarValue.d(ArraySize(TestData(), 2))
Dim StdDevValue.d(ArraySize(TestData(), 2))
Dim CovMatrix.d(ArraySize(TestData(), 2), ArraySize(TestData(), 2))

Debug "Data:"
For k = 0 To ArraySize(TestData(), 1)
  Text = ""
  For i = 0 To ArraySize(TestData(), 2)
    Text + RSet(StrD(TestData(k, i), 2), 10, " ")
  Next i
  Debug Text
Next k
Debug ""

Debug "Means:"
For k = 0 To ArraySize(MeanValue())
  MeanValue(k) = Mean(@TestData(), ArraySize(TestData(), 2) + 1, (ArraySize(TestData(), 2) + 1) * (ArraySize(TestData(), 1) + 1) * SizeOf(Double), k)
  Debug MeanValue(k)
Next k
Debug ""

Debug "Variances:"
For k = 0 To ArraySize(VarValue())
  VarValue(k) = Variance(@TestData(), ArraySize(TestData(), 2) + 1, (ArraySize(TestData(), 2) + 1) * (ArraySize(TestData(), 1) + 1) * SizeOf(Double), MeanValue(k), k)
  Debug VarValue(k)
Next k
Debug ""

Debug "Std.Dev.:"
For k = 0 To ArraySize(StdDevValue())
  StdDevValue(k) = StdDev(VarValue(k))
  Debug StdDevValue(k)
Next k
Debug ""

; subtract the mean of all components
For k = 0 To ArraySize(TestData(), 1)
  For i = 0 To ArraySize(TestData(), 2)
    TestData(k, i) - MeanValue(i)
  Next i
Next k

; calculate the covariance matrix
If CoVarianceMatrix(@TestData(), ArraySize(TestData(), 2) + 1, (ArraySize(TestData(), 2) + 1) * (ArraySize(TestData(), 1) + 1) * SizeOf(Double), CovMatrix())
  Debug "CoVarianceMatrix:"
  For k = 0 To ArraySize(CovMatrix(), 1)
    Text.s = ""
    For i = 0 To ArraySize(CovMatrix(), 1)
      Text + RSet(StrD(CovMatrix(i, k)), 20, " ")
    Next i
    Debug Text
  Next k
  Debug ""
EndIf

NewList Values.d()

If GetEigenValues(Values(), CovMatrix())
  Debug "Eigen-Values:"
  ForEach Values()
    Debug Values()
  Next
  Debug ""
EndIf

NewList Vectors.Vector()

If GetEigenVectors(Vectors(), Values(), CovMatrix())
  Debug "Eigen-Vectors (Sorted by value, so you can directly cut out the feature vector):"
  For k = 0 To ArraySize(CovMatrix(), 2)
    Text.s = ""
    ForEach Vectors()
      Text + RSet(StrD(Vectors()\v(k), 4), 10, " ")
    Next
    Debug Text
  Next k
  Debug ""
EndIf

; > eigen(cov(iris[1:3]))
; $values
; [1] 3.69111979 0.24137727 0.05945372
; 
; $vectors
;             [,1]       [,2]       [,3]
; [1,]  0.38983343  0.6392233 -0.6628903
; [2,] -0.09100801  0.7430587  0.6630093
; [3,]  0.91637735 -0.1981349  0.3478435

;}

[IdeasVacuum]

Wow, thank you DarkDragon. I can see it will take me a while to get my head around it, there is no resemblance to my efforts (which is a good thing!).

[J@ckWhiteIII]

Hi,
I must say that I do not really understand what your code does, DarkDragon. I came to this topic as I was searching for "derivation". Now I've read Polynomial in your code several times, like AddPolynomial etc.
But this topic is about vectors.
So here my questions:
- is this code able to derive polynomials?
- if yes, how are they set?
- what does the code do to the vectors?

I'd be glad to hear from you.
Thank you in advance.

[DarkDragon]

Hi,
I must say that I do not really understand what your code does, DarkDragon. I came to this topic as I was searching for "derivation". Now I've read Polynomial in your code several times, like AddPolynomial etc.
But this topic is about vectors.
So here my questions:
- is this code able to derive polynomials?
- if yes, how are they set?
- what does the code do to the vectors?

I'd be glad to hear from you.
Thank you in advance.

1. Yes, it is able to derive polynomials in real number space through the procedure DerivePolynomial.
2. The polynomials are simple lists of monomes, each consisting of a coefficient and an exponent.
3. Uhm explaining that in detail would take too much time for me now, but the main part of the code is determining the eigenvalues and eigenvectors. The eigenvectors are basically "the orientation" of the data stored inside the vectors. At the end you'll get a transformation matrix and if you apply it to the data stored inside the vectors you'll get a centered and aligned point cloud back, which can be compared to other point clouds without rotational variance or translational variance.

[Psychophanta]

Nice!

To calculate centroid of a planar irregular poligon y tried to translate this one, but never worked:
http://pier.guillen.com.mx/algorithms/0 ... troide.htm
Is there an algorithm to do that, because i am not able to find any functional one.

[IdeasVacuum]

Depends on your purpose - if it doesn't have to be highly accurate, just calculate the bounding box and from that the centroid. The PCA code finds the centroid (STL files for example) but unfortunately for anything more complex/greater mesh density than the test data, the code hangs or crashes. I tried feeding it filtered values (e.g. every tenth triangle) but never got it to work reliably. At the time my project had run out of man hours so I had to abandon the PCA code, but I would very much like to get it working.