Seeking Advice on Implementing a Bidirectional Mean Squares Metric

[Sybren van Rijn]

Hi everyone,

First off, I want to say how much I appreciate this community! 

I'm currently trying to implement a bidirectional mean squares metric for deterministic optimization and was wondering if anyone here has experimented with something similar. If so, I’d love to hear about your experiences or see any implementations you might be willing to share!

Here’s what I’ve tried so far. I implemented the metric as shown below, but I’m encountering two main issues:

1. The computed metric values seem way too high.
2. It produces surface intersections, which I suspect might indicate an underlying issue with how I’m computing the metric.

Here’s my current code:

import scalismo.common.*
import scalismo.geometry.{NDSpace, _3D}
import scalismo.numerics.Sampler
import scalismo.registration.{MeanSquaresMetric, RegistrationMetric}
import scalismo.transformations.TransformationSpace
import breeze.linalg.DenseVector

/**
* A bidirectional mean squares metric for image registration.
*
* @param fixedImage The fixed image in the registration process.
* @param movingImage The moving image in the registration process.
* @param transformationSpace The transformation space used for the registration.
* @param fixedSampler The sampler for the fixed image.
* @param movingSampler The sampler for the moving image.
*/

case class BidirectionalMeanSquaresMetric(
fixedImage: DifferentiableField[_3D, Float],
movingImage: DifferentiableField[_3D, Float],
transformationSpace: TransformationSpace[_3D],
fixedSampler: Sampler[_3D],
movingSampler: Sampler[_3D]
) extends RegistrationMetric[_3D] {

// Metric for the reference to target direction
private val metricReferenceToTarget = MeanSquaresMetric(fixedImage, movingImage, transformationSpace, fixedSampler)
// Metric for the target to reference direction
private val metricTargetToReference = MeanSquaresMetric(movingImage, fixedImage, transformationSpace, movingSampler)

/**
* Computes the value of the bidirectional mean squares metric.
*
* @param parameters The parameters for the transformation.
* @return The computed value of the metric.
*/
override def value(parameters: DenseVector[Double]): Double = {
val mseReferenceToTarget = metricReferenceToTarget.value(parameters)
val mseTargetToReference = metricTargetToReference.value(parameters)
(mseReferenceToTarget + mseTargetToReference) / 2.0
}

/**
* Computes the value and derivative of the bidirectional mean squares metric.
*
* @param parameters The parameters for the transformation.
* @return The computed value and derivative of the metric.
*/
override def valueAndDerivative(parameters: DenseVector[Double]): RegistrationMetric.ValueAndDerivative = {
val valueAndDerivative1: RegistrationMetric.ValueAndDerivative = metricReferenceToTarget.valueAndDerivative(parameters)
val valueAndDerivative2: RegistrationMetric.ValueAndDerivative = metricTargetToReference.valueAndDerivative(parameters)
val combinedValue = (valueAndDerivative1.value + valueAndDerivative2.value) / 2.0
val combinedDerivative = (valueAndDerivative1.derivative + valueAndDerivative2.derivative) / 2.0
RegistrationMetric.ValueAndDerivative(combinedValue, combinedDerivative)
}

/**
* Computes the derivative of the bidirectional mean squares metric.
*
* @param parameters The parameters for the transformation.
* @return The computed derivative of the metric.
*/
override def derivative(parameters: DenseVector[Double]): DenseVector[Double] = {
val valueAndDeriv = valueAndDerivative(parameters)
valueAndDeriv.derivative
}

implicit val ndSpace: NDSpace[_3D] = implicitly[NDSpace[_3D]]
}

I’d really appreciate it if anyone could point me in the right direction or share tips on what might be going wrong here.

Thanks in advance and kind regards,
Sybren

[Marcel Luethi]

Dear Sybren

I haven't worked with this part of scalismo for a long time and therefore you should take my question with the appropriate caution.

Maybe somebody more knowledgeable can chime in.

As far as I remember, computing analytic derivatives for the bidirectional metric is difficult. In one direction, where the take the distance to the target surface, computing the derivative o the full transformation is easy. The other direction would require recomputing the distance map in every gradient step (as the distances change when the model is update). This is why we have always only used the bidirectional distance in gradient free approaches, such as ICP or MCMC. Could this explain your problem?

Best regards,

Marcel

--
You received this message because you are subscribed to the Google Groups "scalismo" group.
To unsubscribe from this group and stop receiving emails from it, send an email to scalismo+u...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/scalismo/ffd4d81b-331c-4a85-8a89-b966732a31c1n%40googlegroups.com.