parametric, non-rigid registration

[Loane Le Gall]

Hi everyone, 

I'm trying to make a non-rigid registration to make point-to-point correspondences between the reference and my all meshes. I succeeded but it's really not perfect. I'm building the ssm of a scapula and I'd like the details, especially at the glenoid, to be precise enough.
From what I've understood from the tutorials, this happens at the level of the kernels, but I tested several combinations and chose the one that seemed best, but it's not great.

I'm attaching my code and a picture of what I can get as a result (scapula in pink is the reuslt of non-rigid registration).

Thank you in advance!

(ps: sorry for my english but i'm a french student.. :) )

 

Here my code :

val zeroMean = Field(EuclideanSpace3D, (_: Point[_3D]) => EuclideanVector.zeros[_3D])
    val kernelCoarse = GaussianKernel3D(86*10)
    val kernelFine = GaussianKernel3D(30*3)
    val kernelfinal = kernelCoarse+kernelFine
    val kernel = DiagonalKernel3D(kernelfinal, outputDim = 3)
    val gp = GaussianProcess(zeroMean, kernel)
    val interpolator = TriangleMeshInterpolator3D[EuclideanVector[_3D]]()
    val lowRankGP = LowRankGaussianProcess.approximateGPCholesky(
      reference,
      gp,
      relativeTolerance = 0.005,
      interpolator = interpolator
    )
    val modelGroup = ui.createGroup("model")
    val gpView = ui.addTransformation(modelGroup, lowRankGP, "gp")

    def doRegistration(
                        lowRankGP: LowRankGaussianProcess[_3D, EuclideanVector[_3D]],
                        referenceMesh: TriangleMesh[_3D],
                        targetICP: TriangleMesh[_3D],
                        initialCoefficients: DenseVector[Double],
                        regParameters: RegistrationParameters): DenseVector[Double] = {
      val transformationSpace = GaussianProcessTransformationSpace(lowRankGP)
      val fixedImage = referenceMesh.operations.toDistanceImage
      val movingImage = targetICP.operations.toDistanceImage
      val sampler = FixedPointsUniformMeshSampler3D(referenceMesh, numberOfPoints = regParameters.numberOfSampledPoints)
      val metric = MeanSquaresMetric(
        fixedImage,
        movingImage,
        transformationSpace,
        sampler)
      val optimizer = LBFGSOptimizer(maxNumberOfIterations = regParameters.numberOfIterations)
      val regularizer = L2Regularizer(transformationSpace)

      val registration = Registration(
        metric,
        regularizer,
        regularizationWeight = regParameters.regularizationWeight,
        optimizer
      )
      val registrationIterator = registration.iterator(initialCoefficients)
      val visualizingRegistrationIterator = for ((it, itnum) <- registrationIterator.zipWithIndex) yield {
        //println(s"object value in iteration $itnum is ${it.value}")
        gpView.coefficients = it.parameters
        it
      }
      val registrationResult = visualizingRegistrationIterator.toSeq.last
      registrationResult.parameters
    }
    val initialCoefficients = DenseVector.zeros[Double](lowRankGP.rank)
    val allPoints = reference.pointSet.numberOfPoints
    case class RegistrationParameters(regularizationWeight: Double, numberOfIterations: Int, numberOfSampledPoints: Int)
    val registrationParameters = Seq(
      RegistrationParameters(regularizationWeight = 1e-1, numberOfIterations = 50, numberOfSampledPoints = 1000),
      RegistrationParameters(regularizationWeight = 1e-2, numberOfIterations = 50, numberOfSampledPoints = 2000),
      RegistrationParameters(regularizationWeight = 1e-4, numberOfIterations = 100, numberOfSampledPoints = 4000),
      RegistrationParameters(regularizationWeight = 1e-6, numberOfIterations = 100, numberOfSampledPoints = 10000)
    )

[Loane Le Gall]

can anyone help me please :) 

I've tried to change, to modify everything I thought was the problem, but it doesn't change anything, or even worse!

thank you in advance

[Dennis Madsen]

Hi Loane, 

Finally found time to look at your problem. 
To me, it mainly seems to be a problem with kernel design. Beside the sigma value, you might also need to scale the individual kernels (see below code). 
If you want to describe all the intricate details of the target you will need a GPMM with hundreds or maybe even thousands of basis functions (model.rank). 
Also, I find that in the example the NearestNeighbourInterpolator worked better for me. 

For the kernel (sigma) values to choose, have a look at the guide I wrote here: https://dennismadsen.me/posts/how-to-shape-model-part4/ 

If the scapula you use is similar to this one (https://commons.wikimedia.org/wiki/File:Human_scapula_1.stl), then it could use some adjustment at least. 

To me, the registration step seems to work fine, but remember that it only optimizes the shape parameters and not the pose. For this, an initial step using the CPD algorithm from GiNGR could maybe be a solution (https://github.com/unibas-gravis/GiNGR/blob/release-1.0/examples/DemoMultiResolution.scala). 

To check the fitting code, I just used a random sample from the GPMM as the target. 

```scala

//> using scala "3.3"

//> using dep "ch.unibas.cs.gravis::scalismo-ui:1.0-RC1"

import scalismo.io.MeshIO

import scalismo.ui.api.ScalismoUI

import java.io.File

import scalismo.kernels.GaussianKernel3D

import scalismo.kernels.DiagonalKernel3D

import scalismo.statisticalmodel.GaussianProcess3D

import scalismo.geometry.EuclideanVector

import scalismo.geometry._3D

import scalismo.common.interpolation.TriangleMeshInterpolator3D

import scalismo.statisticalmodel.LowRankGaussianProcess

import scalismo.statisticalmodel.PointDistributionModel3D

import scalismo.common.interpolation.NearestNeighborInterpolator3D

import scalismo.mesh.TriangleMesh

import breeze.linalg.DenseVector

import scalismo.registration.GaussianProcessTransformationSpace

import scalismo.numerics.FixedPointsUniformMeshSampler3D

import scalismo.registration.MeanSquaresMetric

import scalismo.numerics.LBFGSOptimizer

import scalismo.registration.L2Regularizer

import scalismo.registration.Registration

import scalismo.utils.Random.implicits.randomGenerator

import scala.util.Random

case class RegistrationParameters(regularizationWeight: Double, numberOfIterations: Int, numberOfSampledPoints: Int)

@main def nonrigidFitting() = {

println("Doing stuff")

val raw = MeshIO.readMesh(new File("Human_scapula_1.stl")).get // https://commons.wikimedia.org/wiki/File:Human_scapula_1.stl

val reference = raw.operations.decimate(10000)

// val kernelCoarse = GaussianKernel3D(86*10)

// val kernelFine = GaussianKernel3D(30*3)

val kernelCoarse = GaussianKernel3D(150/2, 15) // 150mm on longest side

val kernelMid = GaussianKernel3D(150/5, 10)

val kernelFine = GaussianKernel3D(150/10, 5)

val kernelfinal = kernelCoarse+kernelMid+kernelFine

val kernel = DiagonalKernel3D(kernelfinal, outputDim = 3)

val gp = GaussianProcess3D[EuclideanVector[_3D]](kernel)

val lowRankGP = LowRankGaussianProcess.approximateGPCholesky(

reference,

gp,

relativeTolerance = 0.01,

interpolator = NearestNeighborInterpolator3D()

)

val gpmm = PointDistributionModel3D(reference, lowRankGP)

println(s"Model with ${gpmm.rank} basis functions and the reference with ${reference.pointSet.numberOfPoints} vertices")

val ui = ScalismoUI()

val modelGroup = ui.createGroup("model")

val modelView = ui.show(modelGroup, gpmm, "model")

// val gpView = ui.addTransformation(modelGroup, lowRankGP, "gp")

modelView.shapeModelTransformationView.shapeTransformationView.coefficients = it.parameters

it

}

val registrationResult = visualizingRegistrationIterator.toSeq.last

registrationResult.parameters

}

val initialCoefficients = DenseVector.zeros[Double](gpmm.rank)

val allPoints = reference.pointSet.numberOfPoints

def createRandomDenseVector(vecSize: Int, minValue: Double = 0.0, maxValue: Double = 0.5): DenseVector[Double] = {

DenseVector.fill(vecSize) {

minValue + (maxValue - minValue) * Random.nextDouble()

}

}

val instancePars = createRandomDenseVector(gpmm.rank)

val target = gpmm.instance(instancePars) // alternative gpmm.sample()

val meshGroup = ui.createGroup("target")

ui.show(meshGroup, target, "target")

ui.show(meshGroup, reference, "reference")

val registrationParameters = Seq(

RegistrationParameters(regularizationWeight = 1e-1, numberOfIterations = 50, numberOfSampledPoints = 100),

RegistrationParameters(regularizationWeight = 1e-2, numberOfIterations = 50, numberOfSampledPoints = 500),

RegistrationParameters(regularizationWeight = 1e-4, numberOfIterations = 50, numberOfSampledPoints = 100),

RegistrationParameters(regularizationWeight = 1e-6, numberOfIterations = 50, numberOfSampledPoints = 5000)

)

val finalCoefficients = registrationParameters.foldLeft(initialCoefficients) { (coefficients, params) =>

println(params)

doRegistration(lowRankGP, reference, target, coefficients, params)

}

val finalFit = gpmm.instance(finalCoefficients)

ui.show(meshGroup, finalFit, "fit")

}

```
To get more details, use a less decimated reference. 

Hope this helps. 

Best,

Dennis