Add 3D

[yakir12]

An ellipse could exist in 3D too. Would be cool to add a method for that.

[feanor12]

I don't know a method for that, but maybe this can be done by rotating the 3D data into a 2D plane. This would make the algorithms reusable.

[KronosTheLate]

That is indeed a general strategy: Use the SVD decomposition (very close to performing a PrincipalComponentDecomposition, or PCA) to determine a normal to the plane in which you have your circle, project the points onto that plane, and use a 2D algorithm. There are however also algorithms that do this for 3D directly.

I am writing here because I am doing a master project, and a fellow student that sits right next to me is working on implementing all kinds of circle fitting algorithms. Some are 2D, some are 3D, but all are intended to work on 3D data. I was just wondering: is this package open to PRs implementing new methods? And before I potentially get to that, I was also wondering: The current API and the examples in the readme seem less-than-ideal to me. If I start implementing such algorithms, would this package be open to reworking the API, revamping the readme, and making a proper documentation page? All in all, the "look and feel" of the package might come out with significant differences.

Further discussion about details is probably best to have in specific pull-requests, e.g. one for API, one for documentation, and one for each new algorithm. But before I start such efforts, I just wanted to check: Is there appetite for such changes?

[feanor12]

Contributions are always welcome!

Regarding the API the initial one was just a collection of functions. I tried to get them into a more streamlined version by looking at other fitting packages. Which resulted in methods like fit or coef. If a new API is required we can for sure adapt and maybe extend the old one to keep it compatible with older code.

More documentation would also be great. As there were not that many methods fully implemented/optimized I just made a Readme, but this can be extended to a documenter.jl page for sure.

[KronosTheLate]

That sounds very good to my ears, and I am finding your reply quite encouraging! I am thinking that using StatsAPI.jl, defining methods for fit, and setting the algorithm as a type (similar to sol = solve(prob, NewtonRaphson()) from SciML's NonlinearSolve.jl) would be a good way to go. I am not 100% sure, but aligning the API with some of the major Julia packages can hardly be a bad idea.

I will see if I get around to it. In the worst case, I expect to at least just provide the functions in their current form to this package. In the best case, I will make PRs for an API rework, new algorithms, and a Documenter.jl page ^_^

[feanor12]

Sounds good to me. :) I already added a setup for the Documenter.jl page yesterday. https://feanor12.github.io/CircleFit.jl/dev/

[KronosTheLate]

Very nice! I can see that I am hard-pressed for time, and will be unlikely to make progress on this in the near future. But it is not forgotten!

[KronosTheLate]

My project is now concluded. I have some other personal business in the upcoming time, and I am not at all certain I will find time for this any time soon. But I just spent an hour translating the method based on conformal geometric algebra (CGA), so that at least a minimal working implementation will come from this. The following code is coarsely translated from python, and no effort has been put into aligning the UI with this package.

using LinearAlgebra

# Simplified helper function to extract parameters from the multivector
function extract_mv_parameters2(mv::Vector{<:Real})
    cvx, cvy, cvz = mv[1:3]   # e2^e3, e3^e1, e1^e2
    cgx, cgy, cgz = mv[4:6]   # e0^e1, e0^e2, e0^e3
    cmx, cmy, cmz = mv[7:9]   # e1^einf, e2^einf, e3^einf
    cgw = -mv[10]             # e0^einf

    num = cvx^2 + cvy^2 + cvz^2 - cgw^2 + 2 * (cgx * cmx + cgy * cmy + cgz * cmz)
    denom = cgx^2 + cgy^2 + cgz^2
    r = sqrt(num / denom)

    n = -[cgx, cgy, cgz]
    nn = n / norm(n)

    B = [cgw -cvz cvy; cvz cgw -cvx; -cvy cvx cgw]
    c = (B * n) / sum(n .^ 2)

    return c, nn, r
end

# Main function to fit a circle in 3D space using CGA
"""
Takes in points in 3D space and outputs a fitted circle in 3D based on a method in Conformal Geometric Algebra (CGA).

    Input is points in 3D euclidean space with each row being a point, and each column is the xyz coordinate.
    The ouput is: centrum, radius, normal_vector
    where circle_cga is the CGA representation of the fitted circle, in the IPNS representation!

    This method is solved by calculating the stacked matrix: M_stacked, which is solved with SVD.

    Example usage:
    >>> centrum, radius, normal_vector, circle_cga = circlefitCGA3D(points)
"""
function circlefitCGA3D(points)
    N = size(points, 1)
    M_stacked = zeros(Float64, 5 * N, 10)

    for i in 1:N
        x4 = 0.5 * sum(points[i, :] .^ 2)

        px, py, pz = points[i, :]
        vec_cross = [0 -pz py; pz 0 -px; -py px 0]

        Mi = [
            -vec_cross  -x4 * I(3)  I(3)  zeros(3, 1);
            zeros(1, 3)  -points[i, :]'  zeros(1, 3)  1;
            zeros(1, 3)  zeros(1, 3)  points[i, :]'  -x4
        ]

        M_stacked[5*(i-1)+1:5*i, :] = Mi
    end

    V = svd(M_stacked, full=false).V
    circle_og = V[:, end]

    return extract_mv_parameters2(circle_og)
end

# Wrapper function to accept separate x, y, z vectors
circlefitCGA3D(xs, ys, zs) = circlefitCGA3D(hcat(xs, ys, zs))

Some code for testing the function

"""
Generates points on a circle in 3D and outputs the 3D euclidean coordinates of the generated points.

The function takes in a normal vector, which expresses the direction of the circle, meaning the direction of the plane the circle lies in,
    radius of the circle, offsets to displace the circle center in 3D, how many points are to be generated, and how large the random noise to be added is.
    The output is a matrix of size num_points x 3, with the format:
    [[x1, y1, z1],
    [x2, y2, z2],
    [x3, y3, z3]...]

    The noise is independent Gaussian noise on the x, y, and z-coordinate.

    Example usage:
    >>> circle_points = generate_circle_points([1.0, 1.0, 1.0], 2.0, [1.0, 1.0, 1.0], 100, 0.05)
"""
    function generate_circle_points(normal_vector::Vector{<:Real}, radius::Real, offsets::Vector{<:Real}, num_points::Integer, noise_scale::Real)
    # Generate points on a 2D circle
    theta = range(0, stop=2π*(num_points-1)/num_points, length=num_points)
    x_circle = radius * cos.(theta)
    y_circle = radius * sin.(theta)

    # Normalize the normal vector
    normal_vector /= norm(normal_vector)

    # Calculate the axis-angle representation for rotation
    # Reference: https://en.wikipedia.org/wiki/Rotation_matrix#Axis_and_angle
    angle = acos(normal_vector[3])  # Angle between z-axis and normal vector
    axis = cross([0.0, 0.0, 1.0], normal_vector)  # Axis of rotation
    axis /= norm(axis)

    # Rotation matrix
    c = cos(angle)
    s = sin(angle)
    t = 1 - c
    rotation_matrix = [
        t*axis[1]^2 + c                 t*axis[1]*axis[2] - s*axis[3]   t*axis[1]*axis[3] + s*axis[2]
        t*axis[1]*axis[2] + s*axis[3]   t*axis[2]^2 + c                 t*axis[2]*axis[3] - s*axis[1]
        t*axis[1]*axis[3] - s*axis[2]   t*axis[2]*axis[3] + s*axis[1]   t*axis[3]^2 + c
    ]

    # Rotate points to align with the specified normal vector
    rotated_points = rotation_matrix * hcat(x_circle, y_circle, zeros(num_points))'
    x_circle, y_circle, z_circle = rotated_points[1, :], rotated_points[2, :], rotated_points[3, :]

    # Set random seed for reproducibility of 'random' noise. Can be commented out.
    # Random.seed!(42)

    # Add noise to the points and the constant offset
    x_circle += randn(num_points) * noise_scale .+ offsets[1]
    y_circle += randn(num_points) * noise_scale .+ offsets[2]
    z_circle += randn(num_points) * noise_scale .+ offsets[3]

    return hcat(x_circle, y_circle, z_circle)
end

let
    # Example usage
    normal_vector = [1, 1, 1]  # Normal vector
    radius_og = 1
    num_points = 100
    offsets = [1, -1, 0.25]

    noise_scale = 0.05
    points = generate_circle_points(normal_vector, radius_og, offsets, num_points, noise_scale)

    # 'Default method for circle fitting, which is SVD on M_stacked
    # centrum, radius, normal_vector = circlefitCGA3D(points)
    centrum, normal_vector, radius = circlefitCGA3D(points)
    @show radius
    @show centrum
    @show normal_vector
end

[feanor12]

Very nice! Do you want to provide the pull request? In case, you are ok with adding your code I can also generate a pull request.

[KronosTheLate]

Feel free to add it - this is rather far back on my priority list right now, so it would have been a while before I got to it xD