diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
new file mode 100644
index 0000000..d09996c
--- /dev/null
+++ b/.github/workflows/ci.yml
@@ -0,0 +1,30 @@
+name: Scala CI
+
+on:
+  push:
+    branches: [ main, release* ]
+  pull_request:
+    branches: [ main, release* ]
+
+jobs:
+  test:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v2
+      - name: Set up JDK 11
+        uses: actions/setup-java@v2
+        with:
+          distribution: 'temurin'
+          java-version: 11
+      - name: Run tests
+        run:  |
+          sudo apt-get update
+          sudo apt-get install -y libopengl0  
+          sbt -Djava.awt.headless=true +compile test
+
+  lint:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v2
+      - name: Formatting
+        run: sbt scalafmtSbtCheck scalafmtCheck test:scalafmtCheck
diff --git a/.gitignore b/.gitignore
index ea73961..0c9d499 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,3 +1,7 @@
+*.tasty
+*.h5
+*.h5.json
+
 # bloop and metals
 .bloop
 .bsp
diff --git a/.scalafmt.conf b/.scalafmt.conf
index 233eb33..d949583 100644
--- a/.scalafmt.conf
+++ b/.scalafmt.conf
@@ -1,5 +1,13 @@
-version = "3.5.3"
-maxColumn = 120
-align = true
-align.preset = more
+version=3.8.0
+project.git = true
 runner.dialect = scala3
+align.openParenCallSite = true
+align.openParenDefnSite = true
+maxColumn = 120
+continuationIndent.defnSite = 2
+assumeStandardLibraryStripMargin = true
+danglingParentheses.preset = true
+rewrite.rules = [SortImports, SortModifiers]
+docstrings.style = Asterisk
+
+onTestFailure = "To fix this, run: sbt scalafmtAll"
diff --git a/README.md b/README.md
index 10b0c9c..550f1f8 100644
--- a/README.md
+++ b/README.md
@@ -13,50 +13,54 @@ BibTex:
 }
 ```
 
-The GiNGR framework allows one to perform non-rigid registration with an iterative algorithm that makes use of Gaussian Process Regression in each iteration to find its next state. 
+The GiNGR framework allows one to perform non-rigid registration with an iterative algorithm that uses Gaussian Process Regression in each iteration to find its next state.
 
-Existing algorithms can be converted into the GiNGR framework, and compared based on 3 properties:
+Existing algorithms can be converted into the GiNGR framework and compared based on three properties:
  
- - Kernel function: how similar should the deformation of neighbouring points be - this is determined based on their correlation
- - Correspondence estimation function: how to estimate corresponding points between the moving instance (reference) and the target.
- - Observation uncertainty: what is the noise assumption of the correspondence estimations?
+ 1. **Kernel function**: how similar should the deformation of neighboring points be - this is determined based on their correlation.
+ 2. **Correspondence estimation function**: how to calculate/estimate anatomically corresponding points between the moving instance (reference) and the target (e.g., use the closest point as a naive approach).
+ 3. **Observation uncertainty**: what is the noise assumption of the correspondence estimations?
 
-This framework contains a general library to input these 3 properties. 
+This framework contains a general library to input these three properties. 
 
-The core part of the GiNGR framework is found in `gingr/api/GingrAlgorithm` with the `update` function performing one iteration of GiNGR update.
-Different pre-implemented configuration files can be found in `gingr/api/registration/config` for CPD and ICP.
+The core part of the GiNGR framework is found in `gingr/api/GingrAlgorithm`, with the `update` function performing one iteration of a GiNGR update.
+Different pre-implemented configuration files can be found under `gingr/api/registration/config` for the CPD and the ICP algorithms.
 
 ## Installation
-To use the GiNGR framework, you can make use of the maven snapshot by including the follow in your Scala 3 script:
+To use the GiNGR framework, you can make use of the maven repository version by including the following in your Scala 3 script:
 
-```
-//> using repository "https://oss.sonatype.org/content/repositories/snapshots"
-//> using lib "ch.unibas.cs.gravis::gingr:0.1.0-SNAPSHOT"
+```scala
+//> using lib "ch.unibas.cs.gravis::gingr:1.0-RC1"
 ```
 
 If using SBT, then add the following to your `build.sbt`:
 
-```
-resolvers +=
-  "Sonatype OSS Snapshots" at "https://oss.sonatype.org/content/repositories/snapshots"
-libraryDependencies += "ch.unibas.cs.gravis" %% "gingr" % "0.1.0-SNAPSHOT"
+```scala
+libraryDependencies += "ch.unibas.cs.gravis" %% "gingr" % "1.0-RC1"
 ```
 
 Or you can install GiNGR to your local repo (_.ivy2_) by running `sbt publishLocal`.
 
-To run the examples, we make use of the VSCODE IDE. For installation help, please see https://scalismo.org/docs/Setup/vscode.
+The examples can be run with [Scala-CLI](https://scala-cli.virtuslab.org/).
+To run the first tutorial, first CD into the example directory and run:
+
+```scala
+# scala-cli project.scala DemoHelper CreateFemurGPMM.scala
+```
+
+Alternatively, the examples can be run with the VSCODE IDE. For installation help, please see https://scalismo.org/docs/Setup/vscode.
 
 After installing VSCODE:
 
  - Go to the examples folder: `cd examples`
  - Setup Code IDE with `scala-cli setup-ide .`
  - Open the IDE with `code .`
- - Now run the individual examples from the "run" menu.
+ - Now run the individual examples
 
 ## General use
-To use GiNGR, one need to specify the deformation model to use in form of a GPMM model as well as the correspondence estimation function and the uncertainty update.
+To use GiNGR, one needs to specify the deformation model as a GPMM model, the correspondence estimation function, and the uncertainty update.
 ### Define the prior model
-The creation of the GPMM is separate from the registration step. For examples, look in the `examples` folder where demo scripts have been created to compute and visualize GPMMs for an Armadillo, Bunny and a Femur bone. In the UI, the deformation model can be evaluated by sampling from it.
+The creation of the GPMM is separate from the registration step. Look in the `examples` folder where demo scripts have been created to compute and visualize GPMMs for an Armadillo, Bunny and a Femur bone. The deformation model can be evaluated in the UI by sampling from it.
 ### Configure the registration algorithm
 The next step is to define the correspondence and uncertainty estimation update. 
 For this, default configurations have been implemented for CPD and ICP. 
@@ -65,7 +69,7 @@ Simple Demo applications can be found in `examples/DemoICP` and `examples/DemoCP
 The demo scripts both perform deterministic and probabilistic registration one after the other. 
 
 ### GiNGR state
-In each iteration a new GiNGR state is computed which contains the GPMM model, the current `fit`, the target as well as all the GPMM model parameters (non-rigid and global pose).
+In each iteration, a new GiNGR state is computed which contains the GPMM model, the current `fit`, the target as well as all the GPMM model parameters (non-rigid and global pose).
 
 ### Deterministic vs Probabilistic
 The probabilistic implementation is based on the ICP-Proposal repository: https://github.com/unibas-gravis/icp-proposal
@@ -81,7 +85,7 @@ First CPD is used on a very coarse mesh (100 vertices), then CPD is used on a me
 
 
 ## Implementation of existing algorithms
-In the GiNGR code base, the basic implementations of existing algorithms can also be found for comparison. The algoriths are found under `gingr/other/algorithms`
+In the GiNGR code base, the basic implementations of existing algorithms can also be found for comparison. The algorithms are found under `gingr/other/algorithms`
 ### CPD: Coherent Point Drift (only Naïve version)
 Implementation of the CPD algorithm from https://arxiv.org/pdf/0905.2635.pdf
 ### BCPD: Bayesian Coherent Point Drift (only Naïve version)
diff --git a/build.sbt b/build.sbt
index be9aeb4..fd2baac 100644
--- a/build.sbt
+++ b/build.sbt
@@ -1,15 +1,17 @@
+ThisBuild / scalaVersion := "3.3.1"
+ThisBuild / version := "1.0-RC1"
+
 lazy val root = (project in file("."))
   .settings(
-    name          := "GiNGR",
-    organization  := "ch.unibas.cs.gravis",
-    scalaVersion  := "3.1.0",
-    scalacOptions := Seq("-deprecation"),
+    name := "GiNGR",
+    homepage := Some(url("https://github.com/unibas-gravis/GiNGR")),
+    organization := "ch.unibas.cs.gravis",
     licenses += ("Apache-2.0", url("http://www.apache.org/licenses/LICENSE-2.0")),
     scmInfo := Some(
       ScmInfo(url("https://github.com/unibas-gravis/GiNGR"), "git@github.com:unibas-gravis/GiNGR.git")
     ),
     developers := List(
-      Developer("madsendennis", "madsendennis", "dennis.madsen@unibas.ch", url("https://github.com/madsendennis"))
+      Developer("madsendennis", "madsendennis", "dennis@dentexion.com", url("https://github.com/madsendennis"))
     ),
     publishMavenStyle := true,
     publishTo := Some(
@@ -18,11 +20,9 @@ lazy val root = (project in file("."))
       else
         Opts.resolver.sonatypeStaging
     ),
-    crossScalaVersions := Seq("2.13.6", "3.1.0"),
     resolvers ++= Seq(
       Resolver.jcenterRepo,
-      Resolver.sonatypeRepo("releases"),
-      Resolver.sonatypeRepo("snapshots")
+      Resolver.sonatypeRepo("releases")
     ),
     scalacOptions ++= {
       Seq(
@@ -32,41 +32,12 @@ lazy val root = (project in file("."))
         "-language:implicitConversions"
         // disabled during the migration
         // "-Xfatal-warnings"
-      ) ++
-        (CrossVersion.partialVersion(scalaVersion.value) match {
-          case Some((3, _)) =>
-            Seq(
-              "-unchecked",
-              "-source:3.0-migration"
-            )
-          case _ =>
-            Seq(
-              "-deprecation",
-              "-Xfatal-warnings",
-              "-Wunused:imports,privates,locals",
-              "-Wvalue-discard"
-            )
-        })
+      )
     },
     javacOptions ++= Seq("-source", "1.8", "-target", "1.8"),
     libraryDependencies ++= Seq(
-      "ch.unibas.cs.gravis" %% "scalismo" % "0.91.+"
-    ),
-    libraryDependencies ++= (scalaBinaryVersion.value match {
-      case "3" =>
-        Seq(
-          "org.scala-lang.modules" %% "scala-parallel-collections" % "1.0.4"
-        )
-      case "2.13" =>
-        Seq(
-          "org.scala-lang.modules" %% "scala-parallel-collections" % "0.2.0"
-        )
-      case _ => { println(scalaBinaryVersion.value); Seq() }
-    })
+      "org.scalatest" %% "scalatest" % "3.2.18" % "test",
+      "ch.unibas.cs.gravis" %% "scalismo" % "1.0-RC1",
+      "io.spray" %% "spray-json" % "1.3.6"
+    )
   )
-//  .enablePlugins(GitVersioning)
-//  .settings(
-//    git.baseVersion := "develop",
-//    git.useGitDescribe := false,
-//    useJGit
-//  )
diff --git a/examples/CreateArmadilloGPMM.scala b/examples/CreateArmadilloGPMM.scala
index ed32ef2..620093a 100644
--- a/examples/CreateArmadilloGPMM.scala
+++ b/examples/CreateArmadilloGPMM.scala
@@ -1,5 +1,3 @@
-//> using scala "3"
-//> using lib "ch.unibas.cs.gravis::scalismo-ui:0.91.2"
 import DemoHelper.DemoDatasetLoader
 import scalismo.ui.api.ScalismoUI
 
diff --git a/examples/CreateBunnyGPMM.scala b/examples/CreateBunnyGPMM.scala
index 5a8a03f..de8a3ca 100644
--- a/examples/CreateBunnyGPMM.scala
+++ b/examples/CreateBunnyGPMM.scala
@@ -1,5 +1,3 @@
-//> using scala "3"
-//> using lib "ch.unibas.cs.gravis::scalismo-ui:0.91.2"
 import DemoHelper.DemoDatasetLoader
 import scalismo.ui.api.ScalismoUI
 
diff --git a/examples/CreateFemurGPMM.scala b/examples/CreateFemurGPMM.scala
index 4d6b090..fd103fd 100644
--- a/examples/CreateFemurGPMM.scala
+++ b/examples/CreateFemurGPMM.scala
@@ -1,5 +1,3 @@
-//> using scala "3"
-//> using lib "ch.unibas.cs.gravis::scalismo-ui:0.91.2"
 import DemoHelper.DemoDatasetLoader
 import scalismo.ui.api.ScalismoUI
 
diff --git a/examples/DemoCPD.scala b/examples/DemoCPD.scala
index 5947251..db4adca 100644
--- a/examples/DemoCPD.scala
+++ b/examples/DemoCPD.scala
@@ -1,7 +1,3 @@
-//> using scala "3"
-//> using repository "https://oss.sonatype.org/content/repositories/snapshots"
-//> using lib "ch.unibas.cs.gravis::gingr:0.1.0-SNAPSHOT"
-//> using lib "ch.unibas.cs.gravis::scalismo-ui:0.91.2"
 import DemoHelper.DemoDatasetLoader
 import DemoHelper.CallBackFunctions.{visualLogger}
 import gingr.simple.GingrInterface
diff --git a/examples/DemoHelper/CallBackFunctions.scala b/examples/DemoHelper/CallBackFunctions.scala
index fc21f5b..7bd5a3f 100644
--- a/examples/DemoHelper/CallBackFunctions.scala
+++ b/examples/DemoHelper/CallBackFunctions.scala
@@ -1,6 +1,3 @@
-//> using scala "3"
-//> using lib "ch.unibas.cs.gravis::gingr:0.1.0-SNAPSHOT"
-//> using lib "ch.unibas.cs.gravis::scalismo-ui:0.91.2"
 package DemoHelper
 
 import gingr.api.GingrRegistrationState
diff --git a/examples/DemoHelper/DemoDatasetLoader.scala b/examples/DemoHelper/DemoDatasetLoader.scala
index 3df55c7..97b6a87 100644
--- a/examples/DemoHelper/DemoDatasetLoader.scala
+++ b/examples/DemoHelper/DemoDatasetLoader.scala
@@ -1,6 +1,3 @@
-//> using scala "3"
-//> using lib "ch.unibas.cs.gravis::gingr:0.1.0-SNAPSHOT"
-//> using lib "ch.unibas.cs.gravis::scalismo:0.91.2"
 package DemoHelper
 
 import scalismo.common.interpolation.TriangleMeshInterpolator3D
@@ -47,7 +44,7 @@ trait DataSetLoader:
   ): (PointDistributionModel[_3D, TriangleMesh], Option[Seq[Landmark[_3D]]]) =
     val (ref, lm) = reference(decimate)
     val decstring = if (decimate.nonEmpty) decimate.get.toString else "full"
-    val modelFile = new File(path, s"${name}_dec-${decstring}_${kernelSelect.name}_${kernelSelect.printpars}.h5")
+    val modelFile = new File(path, s"${name}_dec-${decstring}_${kernelSelect.name}_${kernelSelect.printpars}.h5.json")
     val mdec = StatisticalModelIO.readStatisticalTriangleMeshModel3D(modelFile).getOrElse {
       println(s"Model file does not exist, creating: ${modelFile}")
       val m = SimpleTriangleModels3D.create(ref, kernelSelect, relativeTolerance = relativeTolerance)
@@ -109,8 +106,7 @@ trait DataSetLoader:
 
 
 object DemoDatasetLoader:
-  scalismo.initialize()
-  val dataPath = new File("../data")
+  val dataPath = new File("data")
   case object femur extends DataSetLoader:
     override def name: String        = "femur"
     override def path: File          = new File(dataPath, name)
diff --git a/examples/DemoICP.scala b/examples/DemoICP.scala
index ec73923..51f9030 100644
--- a/examples/DemoICP.scala
+++ b/examples/DemoICP.scala
@@ -1,7 +1,3 @@
-//> using scala "3"
-//> using repository "https://oss.sonatype.org/content/repositories/snapshots"
-//> using lib "ch.unibas.cs.gravis::gingr:0.1.0-SNAPSHOT"
-//> using lib "ch.unibas.cs.gravis::scalismo-ui:0.91.2"
 import DemoHelper.DemoDatasetLoader
 import DemoHelper.CallBackFunctions.{visualLogger}
 import gingr.simple.GingrInterface
diff --git a/examples/DemoLandmarks.scala b/examples/DemoLandmarks.scala
index db154c7..553119d 100644
--- a/examples/DemoLandmarks.scala
+++ b/examples/DemoLandmarks.scala
@@ -1,6 +1,3 @@
-//> using scala "3"
-//> using lib "ch.unibas.cs.gravis::gingr:0.1.0-SNAPSHOT"
-//> using lib "ch.unibas.cs.gravis::scalismo-ui:0.91.2"
 import DemoHelper.DemoDatasetLoader
 import DemoHelper.CallBackFunctions.{visualLogger}
 import gingr.simple.GingrInterface
diff --git a/examples/DemoMultiResolution.scala b/examples/DemoMultiResolution.scala
index 82803ad..e162c1f 100644
--- a/examples/DemoMultiResolution.scala
+++ b/examples/DemoMultiResolution.scala
@@ -1,6 +1,3 @@
-//> using scala "3"
-//> using lib "ch.unibas.cs.gravis::gingr:0.1.0-SNAPSHOT"
-//> using lib "ch.unibas.cs.gravis::scalismo-ui:0.91.2"
 import DemoHelper.DemoDatasetLoader
 import DemoHelper.CallBackFunctions.{visualLogger}
 import gingr.simple.GingrInterface
diff --git a/examples/DemoPosteriorVisualizationFemur.scala b/examples/DemoPosteriorVisualizationFemur.scala
index 1e23df4..38fdc1f 100644
--- a/examples/DemoPosteriorVisualizationFemur.scala
+++ b/examples/DemoPosteriorVisualizationFemur.scala
@@ -1,6 +1,3 @@
-//> using scala "3"
-//> using lib "ch.unibas.cs.gravis::gingr:0.1.0-SNAPSHOT"
-//> using lib "ch.unibas.cs.gravis::scalismo-ui:0.91.2"
 import DemoHelper.DemoDatasetLoader
 import gingr.api.ModelFittingParameters
 import gingr.api.helper.{LogHelper, PosteriorHelper}
diff --git a/examples/VisualizeGPMMCorrelation.scala b/examples/VisualizeGPMMCorrelation.scala
index 4529871..75d776b 100644
--- a/examples/VisualizeGPMMCorrelation.scala
+++ b/examples/VisualizeGPMMCorrelation.scala
@@ -1,6 +1,3 @@
-//> using scala "3"
-//> using lib "ch.unibas.cs.gravis::gingr:0.1.0-SNAPSHOT"
-//> using lib "ch.unibas.cs.gravis::scalismo-ui:0.91.2"
 import DemoHelper.DemoDatasetLoader
 import DemoHelper.CallBackFunctions.{visualLogger}
 import gingr.api.gpmm.GPMMTriangleMesh3D
diff --git a/examples/project.scala b/examples/project.scala
new file mode 100644
index 0000000..1f6dee1
--- /dev/null
+++ b/examples/project.scala
@@ -0,0 +1,3 @@
+//> using scala "3.3"
+//> using dep "ch.unibas.cs.gravis::scalismo-ui:1.0-RC1"
+//> using dep "ch.unibas.cs.gravis::gingr:1.0-RC1"
diff --git a/project/build.properties b/project/build.properties
index 3161d21..abbbce5 100644
--- a/project/build.properties
+++ b/project/build.properties
@@ -1 +1 @@
-sbt.version=1.6.1
+sbt.version=1.9.8
diff --git a/project/plugins.sbt b/project/plugins.sbt
index 3384980..7ba74e8 100644
--- a/project/plugins.sbt
+++ b/project/plugins.sbt
@@ -1,6 +1,3 @@
-resolvers += "sonatype-releases" at "https://oss.sonatype.org/content/repositories/releases/"
-
-addSbtPlugin("com.eed3si9n"     % "sbt-buildinfo" % "0.9.0")
-addSbtPlugin("com.typesafe.sbt" % "sbt-git"       % "1.0.0")
-addSbtPlugin("com.eed3si9n"     % "sbt-assembly"  % "0.14.10")
-addSbtPlugin("org.scalameta"    % "sbt-scalafmt"  % "2.4.3")
+addSbtPlugin("com.jsuereth" % "sbt-pgp" % "1.1.1")
+addSbtPlugin("com.eed3si9n" % "sbt-buildinfo" % "0.9.0")
+addSbtPlugin("org.scalameta" % "sbt-scalafmt" % "2.5.2")
diff --git a/src/main/scala/gingr/api/CorrespondencePairs.scala b/src/main/scala/gingr/api/CorrespondencePairs.scala
index 9480321..5af464b 100644
--- a/src/main/scala/gingr/api/CorrespondencePairs.scala
+++ b/src/main/scala/gingr/api/CorrespondencePairs.scala
@@ -18,7 +18,7 @@
 package gingr.api
 
 import scalismo.common.PointId
-import scalismo.geometry.{Point, _3D}
+import scalismo.geometry.{_3D, Point}
 
 case class CorrespondencePairs(pairs: IndexedSeq[(PointId, Point[_3D])])
 
diff --git a/src/main/scala/gingr/api/GeneralRegistrationState.scala b/src/main/scala/gingr/api/GeneralRegistrationState.scala
index 3458a8e..e2b1a0e 100644
--- a/src/main/scala/gingr/api/GeneralRegistrationState.scala
+++ b/src/main/scala/gingr/api/GeneralRegistrationState.scala
@@ -20,30 +20,30 @@ package gingr.api
 import breeze.linalg.{DenseMatrix, DenseVector}
 import gingr.api.FittingStatuses.FittingStatus
 import scalismo.common.PointId
-import scalismo.geometry.{EuclideanVector, Landmark, Point, _3D}
+import scalismo.geometry.{_3D, EuclideanVector, Landmark, Point}
 import scalismo.mesh.TriangleMesh
 import scalismo.statisticalmodel.{MultivariateNormalDistribution, PointDistributionModel}
 import scalismo.transformations._
 
 case class GeneralRegistrationState(
-    override val model: PointDistributionModel[_3D, TriangleMesh],
-    override val modelParameters: ModelFittingParameters,
-    override val modelLandmarks: Option[Seq[Landmark[_3D]]] = None,
-    override val target: TriangleMesh[_3D],
-    override val targetLandmarks: Option[Seq[Landmark[_3D]]] = None,
-    override val fit: TriangleMesh[_3D],
-    override val sigma2: Double = 1.0,
-    override val globalTransformation: GlobalTranformationType = RigidTransforms,
-    override val stepLength: Double = 1.0,
-    override val generatedBy: String = "",
-    override val iteration: Int = 0,
-    override val status: FittingStatus = FittingStatuses.None
+  override val model: PointDistributionModel[_3D, TriangleMesh],
+  override val modelParameters: ModelFittingParameters,
+  override val modelLandmarks: Option[Seq[Landmark[_3D]]] = None,
+  override val target: TriangleMesh[_3D],
+  override val targetLandmarks: Option[Seq[Landmark[_3D]]] = None,
+  override val fit: TriangleMesh[_3D],
+  override val sigma2: Double = 1.0,
+  override val globalTransformation: GlobalTranformationType = RigidTransforms,
+  override val stepLength: Double = 1.0,
+  override val generatedBy: String = "",
+  override val iteration: Int = 0,
+  override val status: FittingStatus = FittingStatuses.None
 ) extends RegistrationState[GeneralRegistrationState] {
 
   lazy val landmarkCorrespondences: IndexedSeq[(PointId, Point[_3D], MultivariateNormalDistribution)] = {
     if (modelLandmarks.nonEmpty && targetLandmarks.nonEmpty) {
-      val m             = modelLandmarks.get
-      val t             = targetLandmarks.get
+      val m = modelLandmarks.get
+      val t = targetLandmarks.get
       val commonLmNames = m.map(_.id) intersect t.map(_.id)
       commonLmNames.map { name =>
         val mPoint = m.find(_.id == name).get
@@ -61,15 +61,16 @@ case class GeneralRegistrationState(
     }
   }
 
-  /** Updates the current state with the new fit.
-    *
-    * @param next
-    *   The newly calculated shape / fit.
-    */
+  /**
+   * Updates the current state with the new fit.
+   *
+   * @param next
+   *   The newly calculated shape / fit.
+   */
   override def updateFit(next: TriangleMesh[_3D]): GeneralRegistrationState = this.copy(fit = next)
-  override def updateIteration(): GeneralRegistrationState                  = this.copy(iteration = this.iteration + 1)
-  override def clearIteration(): GeneralRegistrationState                   = this.copy(iteration = 0)
-  override def updateStatus(next: FittingStatus): GeneralRegistrationState  = this.copy(status = next)
+  override def updateIteration(): GeneralRegistrationState = this.copy(iteration = this.iteration + 1)
+  override def clearIteration(): GeneralRegistrationState = this.copy(iteration = 0)
+  override def updateStatus(next: FittingStatus): GeneralRegistrationState = this.copy(status = next)
 
   override private[api] def updateTranslation(next: EuclideanVector[_3D]): GeneralRegistrationState = {
     this.copy(modelParameters = this.modelParameters.copy(pose = this.modelParameters.pose.copy(translation = next)))
@@ -80,7 +81,7 @@ case class GeneralRegistrationState(
   }
 
   override private[api] def updateRotation(next: Rotation[_3D]): GeneralRegistrationState = {
-    val angles   = RotationSpace3D.rotMatrixToEulerAngles(next.rotationMatrix)
+    val angles = RotationSpace3D.rotMatrixToEulerAngles(next.rotationMatrix)
     val newEuler = EulerRotation(EulerAngles(angles._1, angles._2, angles._3), next.center)
     this.copy(modelParameters = this.modelParameters.copy(pose = this.modelParameters.pose.copy(rotation = newEuler)))
   }
@@ -115,29 +116,29 @@ case class GeneralRegistrationState(
 
 object GeneralRegistrationState {
   def apply(
-      model: PointDistributionModel[_3D, TriangleMesh],
-      target: TriangleMesh[_3D],
-      modelTranform: Option[TranslationAfterRotation[_3D]]
+    model: PointDistributionModel[_3D, TriangleMesh],
+    target: TriangleMesh[_3D],
+    modelTranform: Option[TranslationAfterRotation[_3D]]
   ): GeneralRegistrationState = {
     apply(model, target, RigidTransforms, modelTranform)
   }
 
   def apply(
-      model: PointDistributionModel[_3D, TriangleMesh],
-      target: TriangleMesh[_3D],
-      transform: GlobalTranformationType,
-      modelTranform: Option[TranslationAfterRotation[_3D]]
+    model: PointDistributionModel[_3D, TriangleMesh],
+    target: TriangleMesh[_3D],
+    transform: GlobalTranformationType,
+    modelTranform: Option[TranslationAfterRotation[_3D]]
   ): GeneralRegistrationState = {
     apply(model, Seq(), target, Seq(), transform, modelTranform)
   }
 
   def apply(
-      model: PointDistributionModel[_3D, TriangleMesh],
-      modelLandmarks: Seq[Landmark[_3D]],
-      target: TriangleMesh[_3D],
-      targetLandmarks: Seq[Landmark[_3D]],
-      transform: GlobalTranformationType,
-      modelTranform: Option[TranslationAfterRotation[_3D]]
+    model: PointDistributionModel[_3D, TriangleMesh],
+    modelLandmarks: Seq[Landmark[_3D]],
+    target: TriangleMesh[_3D],
+    targetLandmarks: Seq[Landmark[_3D]],
+    transform: GlobalTranformationType,
+    modelTranform: Option[TranslationAfterRotation[_3D]]
   ): GeneralRegistrationState = {
     val (t, r) = if (modelTranform.isDefined) {
       val initialAngles = RotationSpace3D.rotMatrixToEulerAngles(modelTranform.get.rotation.rotationMatrix)
@@ -178,11 +179,11 @@ object GeneralRegistrationState {
   }
 
   def apply(
-      model: PointDistributionModel[_3D, TriangleMesh],
-      modelLandmarks: Seq[Landmark[_3D]],
-      target: TriangleMesh[_3D],
-      targetLandmarks: Seq[Landmark[_3D]],
-      modelTranform: Option[TranslationAfterRotation[_3D]]
+    model: PointDistributionModel[_3D, TriangleMesh],
+    modelLandmarks: Seq[Landmark[_3D]],
+    target: TriangleMesh[_3D],
+    targetLandmarks: Seq[Landmark[_3D]],
+    modelTranform: Option[TranslationAfterRotation[_3D]]
   ): GeneralRegistrationState = {
     apply(model, modelLandmarks, target, targetLandmarks, RigidTransforms, modelTranform)
   }
diff --git a/src/main/scala/gingr/api/GingrAlgorithm.scala b/src/main/scala/gingr/api/GingrAlgorithm.scala
index 864647e..cd7450e 100644
--- a/src/main/scala/gingr/api/GingrAlgorithm.scala
+++ b/src/main/scala/gingr/api/GingrAlgorithm.scala
@@ -22,7 +22,7 @@ import gingr.api.sampling.generators.{GeneratorWrapperDeterministic, GeneratorWr
 import gingr.api.sampling.loggers.BestAndCurrentSampleLogger
 import gingr.api.sampling.{AcceptAll, Evaluator, Generator}
 import scalismo.common.PointId
-import scalismo.geometry.{Point, _3D}
+import scalismo.geometry.{_3D, Point}
 import scalismo.mesh.TriangleMesh
 import scalismo.registration.LandmarkRegistration
 import scalismo.sampling.algorithms.MetropolisHastings
@@ -43,8 +43,8 @@ import scalismo.utils.{Memoize, Random}
 import scala.util.Try
 
 case class ProbabilisticSettings[State <: GingrRegistrationState[State]](
-    evaluators: Evaluator[State],
-    randomMixture: Double = 0.5
+  evaluators: Evaluator[State],
+  randomMixture: Double = 0.5
 ) {
   require(randomMixture >= 0.0 && randomMixture <= 1.0)
 }
@@ -65,26 +65,26 @@ trait GingrRegistrationState[State] {
 trait GingrAlgorithm[State <: GingrRegistrationState[State], config <: GingrConfig] {
   val getCorrespondence: State => CorrespondencePairs
   val getUncertainty: (PointId, State) => MultivariateNormalDistribution
-  private val cashedPosterior        = Memoize(computePosterior, 10)
+  private val cashedPosterior = Memoize(computePosterior, 10)
   private val retryCounterInitialize = 10
-  private var retryCounter           = 10
+  private var retryCounter = 10
 
   def name: String
 
   def initializeState(general: GeneralRegistrationState, config: config): State
 
   def estimateRigidTransform(
-      current: TriangleMesh[_3D],
-      update: IndexedSeq[(PointId, Point[_3D])]
+    current: TriangleMesh[_3D],
+    update: IndexedSeq[(PointId, Point[_3D])]
   ): TranslationAfterScalingAfterRotation[_3D] = {
     val pairs = update.map(f => (current.pointSet.point(f._1), f._2))
-    val t     = LandmarkRegistration.rigid3DLandmarkRegistration(pairs, Point(0, 0, 0))
+    val t = LandmarkRegistration.rigid3DLandmarkRegistration(pairs, Point(0, 0, 0))
     TranslationAfterScalingAfterRotation(t.translation, Scaling(1.0), t.rotation)
   }
 
   def estimateSimilarityTransform(
-      current: TriangleMesh[_3D],
-      update: IndexedSeq[(PointId, Point[_3D])]
+    current: TriangleMesh[_3D],
+    update: IndexedSeq[(PointId, Point[_3D])]
   ): TranslationAfterScalingAfterRotation[_3D] = {
     val pairs = update.map(f => (current.pointSet.point(f._1), f._2))
     LandmarkRegistration.similarity3DLandmarkRegistration(pairs, Point(0, 0, 0))
@@ -94,34 +94,35 @@ trait GingrAlgorithm[State <: GingrRegistrationState[State], config <: GingrConf
     model.instance(state.modelParameters.shape.parameters).transform(state.modelParameters.similarityTransform)
   }
 
-  /** Runs the actual registration with the provided configuration through the passed parameters.
-    *
-    * @param initialState
-    *   State from which the registration is started.
-    * @param callBackLogger
-    *   Logger to provide call back functionality to user after each iteration
-    * @param acceptRejectLogger
-    *   Logger to use for advanced file logging
-    * @param probabilisticSettings
-    *   Evaluator to be used if probabilistic registration is set
-    * @param generators
-    *   Pass in external generators to use
-    * @param rnd
-    *   Implicit random number generator.
-    * @return
-    *   Returns the best sample of the chain given the evaluator..
-    */
+  /**
+   * Runs the actual registration with the provided configuration through the passed parameters.
+   *
+   * @param initialState
+   *   State from which the registration is started.
+   * @param callBackLogger
+   *   Logger to provide call back functionality to user after each iteration
+   * @param acceptRejectLogger
+   *   Logger to use for advanced file logging
+   * @param probabilisticSettings
+   *   Evaluator to be used if probabilistic registration is set
+   * @param generators
+   *   Pass in external generators to use
+   * @param rnd
+   *   Implicit random number generator.
+   * @return
+   *   Returns the best sample of the chain given the evaluator..
+   */
   def run(
-      initialState: State,
-      acceptRejectLogger: Option[AcceptRejectLogger[State]],
-      callBackLogger: Option[ChainStateLogger[State]],
-      probabilisticSettings: Option[ProbabilisticSettings[State]],
-      generators: Option[ProposalGenerator[State] with TransitionProbability[State]] = None
+    initialState: State,
+    acceptRejectLogger: Option[AcceptRejectLogger[State]],
+    callBackLogger: Option[ChainStateLogger[State]],
+    probabilisticSettings: Option[ProbabilisticSettings[State]],
+    generators: Option[ProposalGenerator[State] with TransitionProbability[State]] = None
   )(implicit rnd: Random): State = {
-    val evaluator             = probabilisticSettings.getOrElse(ProbabilisticSettings(AcceptAll(), randomMixture = 0.0))
+    val evaluator = probabilisticSettings.getOrElse(ProbabilisticSettings(AcceptAll(), randomMixture = 0.0))
     val registrationEvaluator = EvaluatorWrapper(probabilisticSettings.nonEmpty, evaluator.evaluators)
     val registrationGenerator = generatorCombined(probabilisticSettings, generators)
-    val bestSampleLogger      = BestAndCurrentSampleLogger[State](registrationEvaluator)
+    val bestSampleLogger = BestAndCurrentSampleLogger[State](registrationEvaluator)
     val logs =
       if (callBackLogger.isDefined) ChainStateLoggerContainer(Seq(callBackLogger.get, bestSampleLogger))
       else bestSampleLogger
@@ -133,7 +134,7 @@ trait GingrAlgorithm[State <: GingrRegistrationState[State], config <: GingrConf
       else mhChain.iterator(initialState).loggedWith(logs)
 
     var converged: Boolean = false
-    var error: Boolean     = false
+    var error: Boolean = false
     // we need to query if there is a next element, otherwise due to laziness the chain is not calculated
     var currentState: Option[GeneralRegistrationState] = None
     states
@@ -174,8 +175,8 @@ trait GingrAlgorithm[State <: GingrRegistrationState[State], config <: GingrConf
   }
 
   def generatorCombined(
-      probabilisticSettings: Option[ProbabilisticSettings[State]],
-      mixing: Option[ProposalGenerator[State] with TransitionProbability[State]]
+    probabilisticSettings: Option[ProbabilisticSettings[State]],
+    mixing: Option[ProposalGenerator[State] with TransitionProbability[State]]
   )(implicit rnd: Random): ProposalGenerator[State] with TransitionProbability[State] = {
     probabilisticSettings match {
       case Some(setting) =>
@@ -207,7 +208,7 @@ trait GingrAlgorithm[State <: GingrRegistrationState[State], config <: GingrConf
       }
     } else {
       retryCounter = math.min(retryCounterInitialize, retryCounter + 1)
-      val shapeproposal        = if (!probabilistic) posterior.get.mean else posterior.get.sample()
+      val shapeproposal = if (!probabilistic) posterior.get.mean else posterior.get.sample()
       val transformedModelInit = current.general.model.transform(current.general.modelParameters.rigidTransform)
 
       val newCoefficients = Try(
@@ -240,7 +241,7 @@ trait GingrAlgorithm[State <: GingrRegistrationState[State], config <: GingrConf
             .updateRotation(newGlobalAlignment.rotation)
             .updateScaling(ScaleParameter(globalTransform.scaling.s))
             .updateShapeParameters(ShapeParameters(alpha.get))
-          val newState  = current.updateGeneral(general)
+          val newState = current.updateGeneral(general)
           val newSigma2 = updateSigma2(newState)
           newState.updateGeneral(newState.general.updateSigma2(newSigma2))
         } else {
@@ -257,8 +258,8 @@ trait GingrAlgorithm[State <: GingrRegistrationState[State], config <: GingrConf
   }
 
   def estimateRigidTransform(
-      current: TriangleMesh[_3D],
-      update: TriangleMesh[_3D]
+    current: TriangleMesh[_3D],
+    update: TriangleMesh[_3D]
   ): TranslationAfterScalingAfterRotation[_3D] = {
     val t = LandmarkRegistration.rigid3DLandmarkRegistration(
       current.pointSet.points.toSeq.zip(update.pointSet.points.toSeq),
@@ -268,8 +269,8 @@ trait GingrAlgorithm[State <: GingrRegistrationState[State], config <: GingrConf
   }
 
   def estimateSimilarityTransform(
-      current: TriangleMesh[_3D],
-      update: TriangleMesh[_3D]
+    current: TriangleMesh[_3D],
+    update: TriangleMesh[_3D]
   ): TranslationAfterScalingAfterRotation[_3D] = {
     LandmarkRegistration.similarity3DLandmarkRegistration(
       current.pointSet.points.toSeq.zip(update.pointSet.points.toSeq),
@@ -281,7 +282,7 @@ trait GingrAlgorithm[State <: GingrRegistrationState[State], config <: GingrConf
     val correspondences = getCorrespondence(current)
     val correspondencesWithUncertainty = correspondences.pairs.map { pair =>
       val (pid, point) = pair
-      val uncertainty  = getUncertainty(pid, current)
+      val uncertainty = getUncertainty(pid, current)
       (pid, point, uncertainty)
     }
     val observationsWithUncertainty = if (current.config.useLandmarkCorrespondence) {
diff --git a/src/main/scala/gingr/api/GlobalTranformationType.scala b/src/main/scala/gingr/api/GlobalTranformationType.scala
index 08ae951..9fec908 100644
--- a/src/main/scala/gingr/api/GlobalTranformationType.scala
+++ b/src/main/scala/gingr/api/GlobalTranformationType.scala
@@ -20,5 +20,5 @@ package gingr.api
 sealed trait GlobalTranformationType
 
 case object SimilarityTransforms extends GlobalTranformationType
-case object RigidTransforms      extends GlobalTranformationType
-case object NoTransforms         extends GlobalTranformationType
+case object RigidTransforms extends GlobalTranformationType
+case object NoTransforms extends GlobalTranformationType
diff --git a/src/main/scala/gingr/api/ModelFittingParameters.scala b/src/main/scala/gingr/api/ModelFittingParameters.scala
index e9a39f6..1b7aa5d 100644
--- a/src/main/scala/gingr/api/ModelFittingParameters.scala
+++ b/src/main/scala/gingr/api/ModelFittingParameters.scala
@@ -18,7 +18,7 @@
 package gingr.api
 
 import breeze.linalg.DenseVector
-import scalismo.geometry.{EuclideanVector, Point, _3D}
+import scalismo.geometry.{_3D, EuclideanVector, Point}
 import scalismo.mesh.TriangleMesh
 import scalismo.statisticalmodel.PointDistributionModel
 import scalismo.transformations._
@@ -37,7 +37,7 @@ case class EulerAngles(phi: Double, theta: Double, psi: Double) {
 }
 
 case class EulerRotation(angles: EulerAngles, center: Point[_3D]) {
-  def rotation: Rotation[_3D]         = Rotation(angles.phi, angles.theta, angles.psi, center)
+  def rotation: Rotation[_3D] = Rotation(angles.phi, angles.theta, angles.psi, center)
   def parameters: DenseVector[Double] = DenseVector.vertcat(angles.parameters, center.toBreezeVector)
 }
 
@@ -59,14 +59,14 @@ case class ModelFittingParameters(scale: ScaleParameter, pose: PoseParameters, s
 
   def rigidTransform: TranslationAfterRotation[_3D] = {
     val translation = Translation(pose.translation)
-    val rotation    = pose.rotation.rotation
+    val rotation = pose.rotation.rotation
     TranslationAfterRotation(translation, rotation)
   }
 
   def similarityTransform: TranslationAfterScalingAfterRotation[_3D] = {
     val translation = Translation(pose.translation)
-    val rotation    = pose.rotation.rotation
-    val scaling     = Scaling[_3D](scale.s)
+    val rotation = pose.rotation.rotation
+    val scaling = Scaling[_3D](scale.s)
     TranslationAfterScalingAfterRotation(translation, scaling, rotation)
   }
 
@@ -96,11 +96,11 @@ object ModelFittingParametersJson {
     }
   }
 
-  implicit val myJsonShapePars: RootJsonFormat[ShapeParameters]        = jsonFormat1(ShapeParameters.apply)
-  implicit val myJsonAngPars: RootJsonFormat[EulerAngles]              = jsonFormat3(EulerAngles.apply)
-  implicit val myJsonRotPars: RootJsonFormat[EulerRotation]            = jsonFormat2(EulerRotation.apply)
-  implicit val myJsonPosePars: RootJsonFormat[PoseParameters]          = jsonFormat2(PoseParameters.apply)
-  implicit val myJsonScalePars: RootJsonFormat[ScaleParameter]         = jsonFormat1(ScaleParameter.apply)
+  implicit val myJsonShapePars: RootJsonFormat[ShapeParameters] = jsonFormat1(ShapeParameters.apply)
+  implicit val myJsonAngPars: RootJsonFormat[EulerAngles] = jsonFormat3(EulerAngles.apply)
+  implicit val myJsonRotPars: RootJsonFormat[EulerRotation] = jsonFormat2(EulerRotation.apply)
+  implicit val myJsonPosePars: RootJsonFormat[PoseParameters] = jsonFormat2(PoseParameters.apply)
+  implicit val myJsonScalePars: RootJsonFormat[ScaleParameter] = jsonFormat1(ScaleParameter.apply)
   implicit val myJsonModelPars: RootJsonFormat[ModelFittingParameters] = jsonFormat3(ModelFittingParameters.apply)
 }
 
@@ -128,15 +128,15 @@ object ModelFittingParameters {
   }
 
   def modelInstanceShapePose(
-      model: PointDistributionModel[_3D, TriangleMesh],
-      pars: ModelFittingParameters
+    model: PointDistributionModel[_3D, TriangleMesh],
+    pars: ModelFittingParameters
   ): TriangleMesh[_3D] = {
     model.instance(pars.shape.parameters).transform(pars.rigidTransform)
   }
 
   def modelInstanceShapePoseScale(
-      model: PointDistributionModel[_3D, TriangleMesh],
-      pars: ModelFittingParameters
+    model: PointDistributionModel[_3D, TriangleMesh],
+    pars: ModelFittingParameters
   ): TriangleMesh[_3D] = {
     val instance = modelInstanceShapePose(model, pars)
     instance.transform(pars.scaleTransform)
@@ -153,7 +153,7 @@ object ModelFittingParameters {
   }
 
   def load(file: File): ModelFittingParameters = {
-    val src  = Source.fromFile(file.toString)
+    val src = Source.fromFile(file.toString)
     val data = src.mkString.parseJson.convertTo[ModelFittingParameters]
     src.close()
     data
diff --git a/src/main/scala/gingr/api/RegistrationState.scala b/src/main/scala/gingr/api/RegistrationState.scala
index 2e40f90..ac4ca62 100644
--- a/src/main/scala/gingr/api/RegistrationState.scala
+++ b/src/main/scala/gingr/api/RegistrationState.scala
@@ -18,30 +18,31 @@
 package gingr.api
 
 import gingr.api.FittingStatuses.FittingStatus
-import scalismo.geometry.{EuclideanVector, Landmark, _3D}
+import scalismo.geometry.{_3D, EuclideanVector, Landmark}
 import scalismo.mesh.TriangleMesh
 import scalismo.statisticalmodel.PointDistributionModel
 import scalismo.transformations.Rotation
 
 trait RegistrationState[T] {
   def model: PointDistributionModel[_3D, TriangleMesh] // Prior statistical mesh model
-  def modelParameters: ModelFittingParameters          // parameters of the current fitting state in the model
-  def modelLandmarks: Option[Seq[Landmark[_3D]]]       // Landmarks on the model
-  def target: TriangleMesh[_3D]                        // Target mesh
-  def targetLandmarks: Option[Seq[Landmark[_3D]]]      // Landmarks on the target
+  def modelParameters: ModelFittingParameters // parameters of the current fitting state in the model
+  def modelLandmarks: Option[Seq[Landmark[_3D]]] // Landmarks on the model
+  def target: TriangleMesh[_3D] // Target mesh
+  def targetLandmarks: Option[Seq[Landmark[_3D]]] // Landmarks on the target
   def fit: TriangleMesh[_3D] // Current fit based on model parameters, global alignment and scaling
-  def sigma2: Double         // Global uncertainty parameter
+  def sigma2: Double // Global uncertainty parameter
   def globalTransformation: GlobalTranformationType // Type of global transformation (none, rigid, similarity)
-  def stepLength: Double                            // Step length of a single registration step (0.0 to 1.0)
-  def generatedBy: String                           // Name of generator that produced the State
+  def stepLength: Double // Step length of a single registration step (0.0 to 1.0)
+  def generatedBy: String // Name of generator that produced the State
   def iteration: Int
   def status: FittingStatus // Status
 
-  /** Updates the current state with the new fit.
-    *
-    * @param next
-    *   The newly calculated shape / fit.
-    */
+  /**
+   * Updates the current state with the new fit.
+   *
+   * @param next
+   *   The newly calculated shape / fit.
+   */
   private[api] def updateFit(next: TriangleMesh[_3D]): T
 //  private[api] def updateAlignment(next: TranslationAfterRotation[_3D]): T
   private[api] def updateTranslation(next: EuclideanVector[_3D]): T
diff --git a/src/main/scala/gingr/api/gpmm/GPMMHelper.scala b/src/main/scala/gingr/api/gpmm/GPMMHelper.scala
index 620ded7..1044452 100644
--- a/src/main/scala/gingr/api/gpmm/GPMMHelper.scala
+++ b/src/main/scala/gingr/api/gpmm/GPMMHelper.scala
@@ -21,16 +21,16 @@ import breeze.linalg.sum
 import scalismo.common.DiscreteField.ScalarMeshField
 import scalismo.common._
 import scalismo.common.interpolation.NearestNeighborInterpolator
-import scalismo.geometry.{EuclideanVector, _3D}
+import scalismo.geometry.{_3D, EuclideanVector}
 import scalismo.kernels.{DiagonalKernel, GaussianKernel, MatrixValuedPDKernel}
 import scalismo.mesh.TriangleMesh
 import scalismo.statisticalmodel.{GaussianProcess, LowRankGaussianProcess, PointDistributionModel}
 
 trait GPMM[DDomain[_3D] <: DiscreteDomain[_3D]] {
   def construct(
-      reference: DDomain[_3D],
-      kernel: MatrixValuedPDKernel[_3D],
-      relativeTolerance: Double = 0.01
+    reference: DDomain[_3D],
+    kernel: MatrixValuedPDKernel[_3D],
+    relativeTolerance: Double = 0.01
   ): PointDistributionModel[_3D, DDomain]
 }
 
@@ -38,9 +38,9 @@ object GPMM {
 
   implicit object gpmm3Dtriangle extends GPMM[TriangleMesh] {
     override def construct(
-        reference: TriangleMesh[_3D],
-        kernel: MatrixValuedPDKernel[_3D],
-        relativeTolerance: Double
+      reference: TriangleMesh[_3D],
+      kernel: MatrixValuedPDKernel[_3D],
+      relativeTolerance: Double
     ): PointDistributionModel[_3D, TriangleMesh] = {
       val gp = GaussianProcess[_3D, EuclideanVector[_3D]](kernel)
       val lowRankGP = LowRankGaussianProcess.approximateGPCholesky(
@@ -55,9 +55,9 @@ object GPMM {
 
   implicit object gpmm3Dpoints extends GPMM[UnstructuredPointsDomain] {
     override def construct(
-        reference: UnstructuredPointsDomain[_3D],
-        kernel: MatrixValuedPDKernel[_3D],
-        relativeTolerance: Double
+      reference: UnstructuredPointsDomain[_3D],
+      kernel: MatrixValuedPDKernel[_3D],
+      relativeTolerance: Double
     ): PointDistributionModel[_3D, UnstructuredPointsDomain] = {
       val gp = GaussianProcess[_3D, EuclideanVector[_3D]](kernel)
       val lowRankGP = LowRankGaussianProcess.approximateGPCholesky(
@@ -94,7 +94,7 @@ case class PointSetHelper[DDomain[_3D] <: DiscreteDomain[_3D]](reference: DDomai
 case class GaussianKernelParameters(sigma: Double, scaling: Double)
 
 case class GPMMTriangleMesh3D(reference: TriangleMesh[_3D], relativeTolerance: Double)(implicit
-    gpmm: GPMM[TriangleMesh]
+  gpmm: GPMM[TriangleMesh]
 ) {
   def Gaussian(sigma: Double, scaling: Double): PointDistributionModel[_3D, TriangleMesh] = {
     val kernel = GaussianKernel[_3D](sigma) * scaling
@@ -105,21 +105,21 @@ case class GPMMTriangleMesh3D(reference: TriangleMesh[_3D], relativeTolerance: D
     gpmm.construct(reference, DiagonalKernel(kernel, 3), relativeTolerance)
   }
   def GaussianSymmetry(sigma: Double, scaling: Double): PointDistributionModel[_3D, TriangleMesh] = {
-    val kernel     = GaussianKernel[_3D](sigma) * scaling
+    val kernel = GaussianKernel[_3D](sigma) * scaling
     val symmKernel = KernelHelper.symmetrizeKernel(kernel)
     gpmm.construct(reference, symmKernel, relativeTolerance)
   }
 
   def GaussianMixture(pars: Seq[GaussianKernelParameters]): PointDistributionModel[_3D, TriangleMesh] = {
     val kernels = pars.map(p => GaussianKernel[_3D](p.sigma) * p.scaling)
-    val kernel  = kernels.tail.foldLeft(kernels.head)(_ + _)
+    val kernel = kernels.tail.foldLeft(kernels.head)(_ + _)
     gpmm.construct(reference, DiagonalKernel(kernel, 3), relativeTolerance)
   }
 
   def AutomaticGaussian(): PointDistributionModel[_3D, TriangleMesh] = {
     val refPoints = reference.pointSet
-    val psHelper  = PointSetHelper[TriangleMesh](reference)
-    val maxDist   = psHelper.maximumPointDistance(refPoints)
+    val psHelper = PointSetHelper[TriangleMesh](reference)
+    val maxDist = psHelper.maximumPointDistance(refPoints)
     GaussianMixture(
       Seq(
         GaussianKernelParameters(maxDist / 4.0, maxDist / 8.0),
@@ -129,20 +129,20 @@ case class GPMMTriangleMesh3D(reference: TriangleMesh[_3D], relativeTolerance: D
   }
 
   def InverseLaplacian(scaling: Double): PointDistributionModel[_3D, TriangleMesh] = {
-    val m      = LaplacianHelper(reference).inverseLaplacianMatrix()
+    val m = LaplacianHelper(reference).inverseLaplacianMatrix()
     val kernel = DiagonalKernel(LookupKernel(reference, m) * scaling, 3)
     gpmm.construct(reference, kernel, relativeTolerance)
   }
 
   def InverseLaplacianDot(scaling: Double, gamma: Double): PointDistributionModel[_3D, TriangleMesh] = {
-    val m      = LaplacianHelper(reference).inverseLaplacianMatrix()
+    val m = LaplacianHelper(reference).inverseLaplacianMatrix()
     val kernel = DiagonalKernel(DotProductKernel(LookupKernel(reference, m), gamma) * scaling, 3)
     gpmm.construct(reference, kernel, relativeTolerance)
   }
 
   def computeDistanceAbsMesh(
-      model: PointDistributionModel[_3D, TriangleMesh],
-      lmId: PointId
+    model: PointDistributionModel[_3D, TriangleMesh],
+    lmId: PointId
   ): ScalarMeshField[Double] = {
     val dist: Array[Double] = model.reference.pointSet.pointIds.toSeq.map { pid =>
       val cov = model.gp.cov(lmId, pid)
diff --git a/src/main/scala/gingr/api/gpmm/KernelHelper.scala b/src/main/scala/gingr/api/gpmm/KernelHelper.scala
index 14fc762..27fd822 100644
--- a/src/main/scala/gingr/api/gpmm/KernelHelper.scala
+++ b/src/main/scala/gingr/api/gpmm/KernelHelper.scala
@@ -19,21 +19,21 @@ package gingr.api.gpmm
 
 import breeze.linalg.DenseMatrix
 import scalismo.common.{Domain, EuclideanSpace}
-import scalismo.geometry.{Point, _3D}
+import scalismo.geometry.{_3D, Point}
 import scalismo.kernels.{DiagonalKernel, MatrixValuedPDKernel, PDKernel}
 import scalismo.mesh.TriangleMesh
 
 object KernelHelper {
   def symmetrizeKernel(kernel: PDKernel[_3D]): MatrixValuedPDKernel[_3D] = {
     val xmirrored = xMirroredKernel3D(kernel)
-    val k1        = DiagonalKernel(kernel, 3)
-    val k2        = DiagonalKernel(xmirrored * -1f, xmirrored, xmirrored)
+    val k1 = DiagonalKernel(kernel, 3)
+    val k2 = DiagonalKernel(xmirrored * -1f, xmirrored, xmirrored)
     k1 + k2
   }
 }
 
 case class xMirroredKernel3D(kernel: PDKernel[_3D]) extends PDKernel[_3D] {
-  override def domain: Domain[_3D]                     = kernel.domain
+  override def domain: Domain[_3D] = kernel.domain
   override def k(x: Point[_3D], y: Point[_3D]): Double = kernel(Point(x(0) * -1.0, x(1), x(2)), y)
 }
 
@@ -67,7 +67,7 @@ case class RationalQuadratic(beta: Double) extends PDKernel[_3D] {
   override def domain: EuclideanSpace[_3D] = EuclideanSpace[_3D]
 
   override def k(x: Point[_3D], y: Point[_3D]): Double = {
-    val r  = x - y
+    val r = x - y
     val nn = r.norm2
     1 - (nn / (nn + beta2))
   }
diff --git a/src/main/scala/gingr/api/gpmm/LaplacianHelper.scala b/src/main/scala/gingr/api/gpmm/LaplacianHelper.scala
index 4beae33..d50f8a5 100644
--- a/src/main/scala/gingr/api/gpmm/LaplacianHelper.scala
+++ b/src/main/scala/gingr/api/gpmm/LaplacianHelper.scala
@@ -25,7 +25,7 @@ import scalismo.mesh.TriangleMesh
 
 case class LaplacianHelper(template: TriangleMesh[_3D]) {
   private val n = template.pointSet.numberOfPoints
-  val m         = DenseMatrix.zeros[Double](n, n)
+  val m = DenseMatrix.zeros[Double](n, n)
 
   def deg(i: Int): Int = template.triangulation.adjacentPointsForPoint(PointId(i)).length // get degree of i
   def adj(i: Int, j: Int): Boolean =
diff --git a/src/main/scala/gingr/api/gpmm/MatrixHelper.scala b/src/main/scala/gingr/api/gpmm/MatrixHelper.scala
index 3abf649..da7122f 100644
--- a/src/main/scala/gingr/api/gpmm/MatrixHelper.scala
+++ b/src/main/scala/gingr/api/gpmm/MatrixHelper.scala
@@ -17,7 +17,7 @@
 
 package gingr.api.gpmm
 
-import breeze.linalg.{DenseMatrix, diag, svd}
+import breeze.linalg.{diag, svd, DenseMatrix}
 
 object MatrixHelper {
   def pinv(m: DenseMatrix[Double], precision: Double = 0.00001): DenseMatrix[Double] = {
diff --git a/src/main/scala/gingr/api/helper/CallBackFunctions.scala b/src/main/scala/gingr/api/helper/CallBackFunctions.scala
index 40c94a7..c8143b5 100644
--- a/src/main/scala/gingr/api/helper/CallBackFunctions.scala
+++ b/src/main/scala/gingr/api/helper/CallBackFunctions.scala
@@ -24,8 +24,8 @@ import scalismo.sampling.loggers.ChainStateLogger
 object CallBackFunctions {
 
   case class SimpleLogger[State <: GingrRegistrationState[State]](
-      jsonLogger: Option[JSONStateLogger[State]] = None,
-      printUpdateFrequency: Int = 100
+    jsonLogger: Option[JSONStateLogger[State]] = None,
+    printUpdateFrequency: Int = 100
   ) extends ChainStateLogger[State] {
     var counter = 0
     override def logState(sample: State): Unit = {
diff --git a/src/main/scala/gingr/api/helper/LogHelper.scala b/src/main/scala/gingr/api/helper/LogHelper.scala
index 151ceef..57f4c94 100644
--- a/src/main/scala/gingr/api/helper/LogHelper.scala
+++ b/src/main/scala/gingr/api/helper/LogHelper.scala
@@ -18,7 +18,7 @@
 package gingr.api.helper
 
 import gingr.api.ModelFittingParameters
-import gingr.api.sampling.loggers.{JSONStateLogger, jsonLogFormat}
+import gingr.api.sampling.loggers.{jsonLogFormat, JSONStateLogger}
 import scalismo.geometry._3D
 import scalismo.mesh.TriangleMesh
 import scalismo.statisticalmodel.PointDistributionModel
@@ -28,10 +28,10 @@ import scala.annotation.tailrec
 object LogHelper {
 
   def samplesFromLog(
-      log: IndexedSeq[jsonLogFormat],
-      takeEveryN: Int = 50,
-      total: Int = 100,
-      burnIn: Int = 0
+    log: IndexedSeq[jsonLogFormat],
+    takeEveryN: Int = 50,
+    total: Int = 100,
+    burnIn: Int = 0
   ): IndexedSeq[(jsonLogFormat, Int)] = {
     @tailrec
     def getLogIndex(i: Int): Int = {
@@ -39,14 +39,14 @@ object LogHelper {
       else getLogIndex(i - 1)
     }
     println("Log length: " + log.length)
-    val indexes  = (burnIn until math.min(log.length, total) by takeEveryN).map(i => getLogIndex(i))
+    val indexes = (burnIn until math.min(log.length, total) by takeEveryN).map(i => getLogIndex(i))
     val filtered = indexes.map(i => (log(i), i))
     filtered.take(math.min(total, filtered.length))
   }
 
   def logSamples2shapes(
-      model: PointDistributionModel[_3D, TriangleMesh],
-      log: IndexedSeq[jsonLogFormat]
+    model: PointDistributionModel[_3D, TriangleMesh],
+    log: IndexedSeq[jsonLogFormat]
   ): IndexedSeq[TriangleMesh[_3D]] = {
     log.map { l =>
       val modelPars = JSONStateLogger.jsonFormatToModelFittingParameters(l)
diff --git a/src/main/scala/gingr/api/helper/PosteriorHelper.scala b/src/main/scala/gingr/api/helper/PosteriorHelper.scala
index 44ce404..105287f 100644
--- a/src/main/scala/gingr/api/helper/PosteriorHelper.scala
+++ b/src/main/scala/gingr/api/helper/PosteriorHelper.scala
@@ -17,20 +17,20 @@
 
 package gingr.api.helper
 
-import breeze.linalg.{DenseMatrix, DenseVector, trace}
+import breeze.linalg.{trace, DenseMatrix, DenseVector}
 import scalismo.common.DiscreteField.ScalarMeshField
 import scalismo.common.ScalarMeshField
 import scalismo.geometry._3D
 import scalismo.mesh.TriangleMesh
 
 object PosteriorHelper {
-  private val sampleDim  = 3
-  private val zeroVec    = DenseVector.zeros[Double](sampleDim)
+  private val sampleDim = 3
+  private val zeroVec = DenseVector.zeros[Double](sampleDim)
   private val zeroMatrix = DenseMatrix.zeros[Double](sampleDim, sampleDim)
 
   def computeDistanceMapFromMeshesTotal(
-      meshes: IndexedSeq[TriangleMesh[_3D]],
-      ref: TriangleMesh[_3D]
+    meshes: IndexedSeq[TriangleMesh[_3D]],
+    ref: TriangleMesh[_3D]
   ): ScalarMeshField[Double] = {
     val numSamples = meshes.length
 
@@ -42,8 +42,8 @@ object PosteriorHelper {
     }
 
     val variance = pointIdPointSeq.map { samples =>
-      val mean         = samples.foldLeft(zeroVec)((acc, s) => acc + s) * (1.0 / numSamples)
-      val cov          = samples.foldLeft(zeroMatrix)((acc, s) => acc + outer(s - mean, s - mean)) * (1.0 / (numSamples - 1))
+      val mean = samples.foldLeft(zeroVec)((acc, s) => acc + s) * (1.0 / numSamples)
+      val cov = samples.foldLeft(zeroMatrix)((acc, s) => acc + outer(s - mean, s - mean)) * (1.0 / (numSamples - 1))
       val sum_variance = trace(cov)
 
       (mean, cov, sum_variance)
@@ -53,9 +53,9 @@ object PosteriorHelper {
   }
 
   def computeDistanceMapFromMeshesNormal(
-      meshes: IndexedSeq[TriangleMesh[_3D]],
-      ref: TriangleMesh[_3D],
-      sumNormals: Boolean = true
+    meshes: IndexedSeq[TriangleMesh[_3D]],
+    ref: TriangleMesh[_3D],
+    sumNormals: Boolean = true
   ): ScalarMeshField[Double] = {
     val numSamples = meshes.length
 
diff --git a/src/main/scala/gingr/api/helper/RegistrationComparison.scala b/src/main/scala/gingr/api/helper/RegistrationComparison.scala
index eacdcea..08fec57 100644
--- a/src/main/scala/gingr/api/helper/RegistrationComparison.scala
+++ b/src/main/scala/gingr/api/helper/RegistrationComparison.scala
@@ -35,22 +35,22 @@ object RegistrationComparison {
   }
 
   def evaluateReconstruction2GroundTruth(
-      id: String,
-      reconstruction: TriangleMesh3D,
-      groundTruth: TriangleMesh3D
+    id: String,
+    reconstruction: TriangleMesh3D,
+    groundTruth: TriangleMesh3D
   ): (Double, Double, Double) = {
     val avgDist2Surf = MeshMetrics.avgDistance(reconstruction, groundTruth)
 
     val hausdorffDistance = MeshMetrics.hausdorffDistance(reconstruction, groundTruth)
-    val md                = maxDistance(reconstruction, groundTruth)
+    val md = maxDistance(reconstruction, groundTruth)
     println(s"ID: ${id} average2surface: ${avgDist2Surf} max: ${md}, hausdorff: ${hausdorffDistance}")
     (avgDist2Surf, md, hausdorffDistance)
   }
 
   def evaluateReconstruction2GroundTruthDouble(
-      id: String,
-      reconstruction: TriangleMesh3D,
-      groundTruth: TriangleMesh3D
+    id: String,
+    reconstruction: TriangleMesh3D,
+    groundTruth: TriangleMesh3D
   ): (Double, Double) = {
     val avgDist2Surf = (MeshMetrics
       .avgDistance(reconstruction, groundTruth) + MeshMetrics.avgDistance(groundTruth, reconstruction)) / 2.0
@@ -64,7 +64,7 @@ object RegistrationComparison {
 
     val pointsOnSample = m1.pointSet.points
     val dists = for (p <- pointsOnSample) yield {
-      val pTarget   = m2.operations.closestPointOnSurface(p).point
+      val pTarget = m2.operations.closestPointOnSurface(p).point
       val pTargetId = m2.pointSet.findClosestPoint(pTarget).id
       if (m2.operations.pointIsOnBoundary(pTargetId)) None
       else Option((pTarget - p).norm)
@@ -74,14 +74,14 @@ object RegistrationComparison {
   }
 
   def evaluateReconstruction2GroundTruthBoundaryAware(
-      id: String,
-      reconstruction: TriangleMesh3D,
-      groundTruth: TriangleMesh3D
+    id: String,
+    reconstruction: TriangleMesh3D,
+    groundTruth: TriangleMesh3D
   ): Unit = {
     val (avgDist2Surf1, maxDist2Surf1) = avgDistanceBoundaryAware(reconstruction, groundTruth)
     val (avgDist2Surf2, maxDist2Surf2) = avgDistanceBoundaryAware(groundTruth, reconstruction)
-    val avgDist2Surf                   = (avgDist2Surf1 + avgDist2Surf2) / 2.0
-    val maxDist2Surf                   = math.max(maxDist2Surf1, maxDist2Surf2)
+    val avgDist2Surf = (avgDist2Surf1 + avgDist2Surf2) / 2.0
+    val maxDist2Surf = math.max(maxDist2Surf1, maxDist2Surf2)
     println(s"ID: ${id} average2surface: ${avgDist2Surf} max: ${maxDist2Surf}")
   }
 
diff --git a/src/main/scala/gingr/api/registration/SimpleRegistrator.scala b/src/main/scala/gingr/api/registration/SimpleRegistrator.scala
index a842cc7..2bc8d07 100644
--- a/src/main/scala/gingr/api/registration/SimpleRegistrator.scala
+++ b/src/main/scala/gingr/api/registration/SimpleRegistrator.scala
@@ -33,7 +33,7 @@ import gingr.api.{
   RigidTransforms
 }
 import scalismo.common.interpolation.NearestNeighborInterpolator
-import scalismo.geometry.{Landmark, _3D}
+import scalismo.geometry.{_3D, Landmark}
 import scalismo.mesh.TriangleMesh
 import scalismo.sampling.loggers.ChainStateLogger
 import scalismo.statisticalmodel.PointDistributionModel
@@ -46,27 +46,27 @@ class SimpleRegistrator[State <: GingrRegistrationState[State], Config <: GingrC
   State,
   Config
 ]](
-    algorithm: Algorithm,
-    config: Config,
-    model: PointDistributionModel[_3D, TriangleMesh],
-    target: TriangleMesh[_3D],
-    initialModelParameterTransform: Option[TranslationAfterRotation[_3D]] = None,
-    modelLandmarks: Option[Seq[Landmark[_3D]]] = None,
-    targetLandmarks: Option[Seq[Landmark[_3D]]] = None,
-    evaluationMode: EvaluationMode = ModelToTargetEvaluation,
-    evaluatorUncertainty: Double = 1.0,
-    evaluatedPoints: Option[Int] = None,
-    logFileFittingParameters: Option[File] = None
+  algorithm: Algorithm,
+  config: Config,
+  model: PointDistributionModel[_3D, TriangleMesh],
+  target: TriangleMesh[_3D],
+  initialModelParameterTransform: Option[TranslationAfterRotation[_3D]] = None,
+  modelLandmarks: Option[Seq[Landmark[_3D]]] = None,
+  targetLandmarks: Option[Seq[Landmark[_3D]]] = None,
+  evaluationMode: EvaluationMode = ModelToTargetEvaluation,
+  evaluatorUncertainty: Double = 1.0,
+  evaluatedPoints: Option[Int] = None,
+  logFileFittingParameters: Option[File] = None
 )(implicit rnd: Random) {
 
   def runDecimated(
-      modelPoints: Int,
-      targetPoints: Int,
-      generalState: Option[GeneralRegistrationState] = None,
-      globalTransformation: GlobalTranformationType = RigidTransforms,
-      probabilistic: Boolean = false,
-      randomMixture: Double = 0.5,
-      callback: Option[ChainStateLogger[State]] = None
+    modelPoints: Int,
+    targetPoints: Int,
+    generalState: Option[GeneralRegistrationState] = None,
+    globalTransformation: GlobalTranformationType = RigidTransforms,
+    probabilistic: Boolean = false,
+    randomMixture: Double = 0.5,
+    callback: Option[ChainStateLogger[State]] = None
   ): State = {
     val initState = decimateState(generalState, globalTransformation, modelPoints, targetPoints)
     run(Option(initState), globalTransformation, probabilistic, randomMixture, callback)
@@ -82,13 +82,13 @@ class SimpleRegistrator[State <: GingrRegistrationState[State], Config <: GingrC
   }
 
   private def decimateState(
-      generalState: Option[GeneralRegistrationState],
-      globalTransformation: GlobalTranformationType,
-      modelPoints: Int,
-      targetPoints: Int
+    generalState: Option[GeneralRegistrationState],
+    globalTransformation: GlobalTranformationType,
+    modelPoints: Int,
+    targetPoints: Int
   ): GeneralRegistrationState = {
-    val newRef          = model.reference.operations.decimate(modelPoints)
-    val decimatedModel  = model.newReference(newRef, NearestNeighborInterpolator())
+    val newRef = model.reference.operations.decimate(modelPoints)
+    val decimatedModel = model.newReference(newRef, NearestNeighborInterpolator())
     val decimatedTarget = target.operations.decimate(targetPoints)
     val initState =
       if (generalState.isDefined) combineStates(generalState.get)
@@ -105,10 +105,10 @@ class SimpleRegistrator[State <: GingrRegistrationState[State], Config <: GingrC
   }
 
   def createInitialState(
-      model: PointDistributionModel[_3D, TriangleMesh],
-      target: TriangleMesh[_3D],
-      globalTransformation: GlobalTranformationType,
-      modelTranform: Option[TranslationAfterRotation[_3D]] = None
+    model: PointDistributionModel[_3D, TriangleMesh],
+    target: TriangleMesh[_3D],
+    globalTransformation: GlobalTranformationType,
+    modelTranform: Option[TranslationAfterRotation[_3D]] = None
   ): State = {
     val generalState: GeneralRegistrationState =
       if (modelLandmarks.nonEmpty && targetLandmarks.nonEmpty)
@@ -126,11 +126,11 @@ class SimpleRegistrator[State <: GingrRegistrationState[State], Config <: GingrC
   }
 
   def run(
-      generalState: Option[GeneralRegistrationState] = None,
-      globalTransformation: GlobalTranformationType = RigidTransforms,
-      probabilistic: Boolean = false,
-      randomMixture: Double = 0.5,
-      callback: Option[ChainStateLogger[State]] = None
+    generalState: Option[GeneralRegistrationState] = None,
+    globalTransformation: GlobalTranformationType = RigidTransforms,
+    probabilistic: Boolean = false,
+    randomMixture: Double = 0.5,
+    callback: Option[ChainStateLogger[State]] = None
   ): State = {
     val state =
       if (generalState.isDefined) combineStates(generalState.get)
diff --git a/src/main/scala/gingr/api/registration/config/CPD.scala b/src/main/scala/gingr/api/registration/config/CPD.scala
index 360aa58..70dadeb 100644
--- a/src/main/scala/gingr/api/registration/config/CPD.scala
+++ b/src/main/scala/gingr/api/registration/config/CPD.scala
@@ -16,13 +16,12 @@
  */
 
 package gingr.api.registration.config
-
-import breeze.linalg.{Axis, DenseMatrix, DenseVector, sum, tile}
+import breeze.linalg.{sum, tile, Axis, DenseMatrix, DenseVector}
 import gingr.api.registration.utils.PointSequenceConverter
 import gingr.api.{CorrespondencePairs, GeneralRegistrationState, GingrAlgorithm, GingrConfig, GingrRegistrationState}
 import scalismo.common.PointId
 import scalismo.geometry.Point.Point3DVectorizer
-import scalismo.geometry.{Point, _3D}
+import scalismo.geometry.{_3D, Point}
 import scalismo.statisticalmodel.MultivariateNormalDistribution
 
 import scala.collection.parallel.CollectionConverters.ImmutableSeqIsParallelizable
@@ -33,16 +32,18 @@ object CPDCorrespondence {
   def estimate[T](state: CpdRegistrationState): CorrespondencePairs = {
     val refPoints = state.general.fit.pointSet.points.toSeq
     val tarPoints = state.general.target.pointSet.points.toSeq
-    val P         = state.P
-    val P1inv     = 1.0 / sum(P, Axis._1)
+    val P = state.P
+    val Psum = sum(P, Axis._1)
+    val P1inv = Psum.map(p => 1.0 / p)
 
     val deform = refPoints.zipWithIndex.par.map { case (y, i) =>
       val xscale = tarPoints.zipWithIndex.map { case (x, j) =>
-        P1inv(i) * P(i, j) * x.toBreezeVector
+        val tmp = P1inv(i) * P(i, j)
+        x.toBreezeVector.map(p => tmp * p)
       }
       vectorizer.unvectorize(sum(xscale) - y.toBreezeVector).toVector
     }
-    val td   = refPoints.zip(deform).map { case (p, d) => p + d }
+    val td = refPoints.zip(deform).map { case (p, d) => p + d }
     val corr = state.general.fit.pointSet.pointIds.toSeq.zip(td).map(t => (t._1, t._2)).toIndexedSeq
     CorrespondencePairs(corr)
   }
@@ -56,8 +57,8 @@ case class CpdRegistrationState(general: GeneralRegistrationState, config: CpdCo
     }
     val refPoints = general.fit.pointSet.points.toSeq
     val tarPoints = general.target.pointSet.points.toSeq
-    val M         = refPoints.length
-    val N         = tarPoints.length
+    val M = refPoints.length
+    val N = tarPoints.length
     // TODO: Approximate using nyström
     val P: DenseMatrix[Double] = DenseMatrix.zeros[Double](M, N)
     refPoints.zipWithIndex.par.foreach { case (y, i) =>
@@ -68,7 +69,7 @@ case class CpdRegistrationState(general: GeneralRegistrationState, config: CpdCo
     val c =
       config.w / (1 - config.w) * math.pow((2.0 * math.Pi * general.sigma2), 3.0 / 2.0) * (M.toDouble / N.toDouble)
     val denRow = DenseMatrix(sum(P, Axis._0).t)
-    val den    = tile(denRow, M, 1) + c
+    val den = tile(denRow, M, 1) + c
 
     P /:/ den
   }
@@ -102,52 +103,52 @@ object CpdRegistrationState {
 }
 
 case class CpdConfiguration(
-    override val maxIterations: Int = 100,
-    override val threshold: Double = 1e-10,
-    override val converged: (GeneralRegistrationState, GeneralRegistrationState, Double) => Boolean =
-      (last: GeneralRegistrationState, current: GeneralRegistrationState, threshold: Double) =>
-        math.abs(last.sigma2 - current.sigma2) < threshold,
-    override val useLandmarkCorrespondence: Boolean = true,
-    initialSigma: Option[Double] = None,
-    w: Double = 0.0,
-    lambda: Double = 1.0
+  override val maxIterations: Int = 100,
+  override val threshold: Double = 1e-10,
+  override val converged: (GeneralRegistrationState, GeneralRegistrationState, Double) => Boolean =
+    (last: GeneralRegistrationState, current: GeneralRegistrationState, threshold: Double) =>
+      math.abs(last.sigma2 - current.sigma2) < threshold,
+  override val useLandmarkCorrespondence: Boolean = true,
+  initialSigma: Option[Double] = None,
+  w: Double = 0.0,
+  lambda: Double = 1.0
 ) extends GingrConfig {}
 
 class CpdRegistration(
-    override val getCorrespondence: CpdRegistrationState => CorrespondencePairs = (state: CpdRegistrationState) =>
-      CPDCorrespondence.estimate(state),
-    override val getUncertainty: (PointId, CpdRegistrationState) => MultivariateNormalDistribution =
-      (id: PointId, state: CpdRegistrationState) => {
-        val P     = state.P
-        val P1inv = 1.0 / sum(P, Axis._1)
-        MultivariateNormalDistribution(
-          DenseVector.zeros[Double](3),
-          DenseMatrix.eye[Double](3) * state.general.sigma2 * state.config.lambda * P1inv(id.id)
-        )
-      }
+  override val getCorrespondence: CpdRegistrationState => CorrespondencePairs = (state: CpdRegistrationState) =>
+    CPDCorrespondence.estimate(state),
+  override val getUncertainty: (PointId, CpdRegistrationState) => MultivariateNormalDistribution =
+    (id: PointId, state: CpdRegistrationState) => {
+      val P = state.P
+      val P1inv = sum(P, Axis._1).map(p => 1.0 / p)
+      MultivariateNormalDistribution(
+        DenseVector.zeros[Double](3),
+        DenseMatrix.eye[Double](3) * state.general.sigma2 * state.config.lambda * P1inv(id.id)
+      )
+    }
 ) extends GingrAlgorithm[CpdRegistrationState, CpdConfiguration] {
   def name = "CPD"
 
   // possibility to override the update function, or just use the base class method?
   override def updateSigma2(current: CpdRegistrationState): Double = {
     val meanUpdate = current.general.fit.pointSet.points.toSeq
-    val P          = current.P
-    val X          = PointSequenceConverter.toMatrix(current.general.target.pointSet.points.toSeq)
-    val TY         = PointSequenceConverter.toMatrix(meanUpdate)
-    val P1         = sum(P, Axis._1)
-    val Pt1        = sum(P, Axis._0)
-    val Np         = sum(P1)
-
-    val xPx: Double   = Pt1.t.dot(sum(X *:* X, Axis._1))
-    val yPy: Double   = P1.t * sum(TY *:* TY, Axis._1)
+    val P = current.P
+    val X = PointSequenceConverter.toMatrix(current.general.target.pointSet.points.toSeq)
+    val TY = PointSequenceConverter.toMatrix(meanUpdate)
+    val P1 = sum(P, Axis._1)
+    val Pt1 = sum(P, Axis._0)
+    val Np = sum(P1)
+
+    val xPx: Double = Pt1.t.dot(sum(X *:* X, Axis._1))
+    val yPy: Double = P1.t * sum(TY *:* TY, Axis._1)
     val trPXY: Double = sum(TY *:* (P * X))
-    val sigma2        = (xPx - 2 * trPXY + yPy) / (Np * 3.0)
+    val sigma2 = (xPx - 2 * trPXY + yPy) / (Np * 3.0)
     sigma2
   }
 
   override def initializeState(
-      general: GeneralRegistrationState,
-      config: CpdConfiguration
+    general: GeneralRegistrationState,
+    config: CpdConfiguration
   ): CpdRegistrationState = {
     CpdRegistrationState(general, config)
   }
diff --git a/src/main/scala/gingr/api/registration/config/ICP.scala b/src/main/scala/gingr/api/registration/config/ICP.scala
index 7ba6cac..c1f3d42 100644
--- a/src/main/scala/gingr/api/registration/config/ICP.scala
+++ b/src/main/scala/gingr/api/registration/config/ICP.scala
@@ -29,9 +29,9 @@ import scalismo.statisticalmodel.MultivariateNormalDistribution
 
 sealed trait ICPCorrespondenceMethod
 
-case object TriangularClosestPoint  extends ICPCorrespondenceMethod
+case object TriangularClosestPoint extends ICPCorrespondenceMethod
 case object AlongNormalClosestPoint extends ICPCorrespondenceMethod
-case object PointcloudClosestPoint  extends ICPCorrespondenceMethod
+case object PointcloudClosestPoint extends ICPCorrespondenceMethod
 
 object ICPCorrespondence {
   def estimate[T](state: IcpRegistrationState): CorrespondencePairs = {
@@ -52,15 +52,15 @@ object ICPCorrespondence {
 }
 
 case class IcpConfiguration(
-    override val maxIterations: Int = 100,
-    override val threshold: Double = 1e-10,
-    override val converged: (GeneralRegistrationState, GeneralRegistrationState, Double) => Boolean =
-      (last: GeneralRegistrationState, current: GeneralRegistrationState, threshold: Double) => false,
-    override val useLandmarkCorrespondence: Boolean = true,
-    initialSigma: Double = 100.0,
-    endSigma: Double = 1.0,
-    reverseCorrespondenceDirection: Boolean = false,
-    correspondenceMethod: ICPCorrespondenceMethod = TriangularClosestPoint
+  override val maxIterations: Int = 100,
+  override val threshold: Double = 1e-10,
+  override val converged: (GeneralRegistrationState, GeneralRegistrationState, Double) => Boolean =
+    (last: GeneralRegistrationState, current: GeneralRegistrationState, threshold: Double) => false,
+  override val useLandmarkCorrespondence: Boolean = true,
+  initialSigma: Double = 100.0,
+  endSigma: Double = 1.0,
+  reverseCorrespondenceDirection: Boolean = false,
+  correspondenceMethod: ICPCorrespondenceMethod = TriangularClosestPoint
 ) extends GingrConfig {
   val sigmaStep: Double = (initialSigma - endSigma) / maxIterations.toDouble
 }
@@ -85,11 +85,11 @@ object IcpRegistrationState {
 }
 
 class IcpRegistration(
-    override val getCorrespondence: IcpRegistrationState => CorrespondencePairs = (state: IcpRegistrationState) =>
-      ICPCorrespondence.estimate(state),
-    override val getUncertainty: (PointId, IcpRegistrationState) => MultivariateNormalDistribution =
-      (id: PointId, state: IcpRegistrationState) =>
-        MultivariateNormalDistribution(DenseVector.zeros[Double](3), DenseMatrix.eye[Double](3) * state.general.sigma2)
+  override val getCorrespondence: IcpRegistrationState => CorrespondencePairs = (state: IcpRegistrationState) =>
+    ICPCorrespondence.estimate(state),
+  override val getUncertainty: (PointId, IcpRegistrationState) => MultivariateNormalDistribution =
+    (id: PointId, state: IcpRegistrationState) =>
+      MultivariateNormalDistribution(DenseVector.zeros[Double](3), DenseMatrix.eye[Double](3) * state.general.sigma2)
 ) extends GingrAlgorithm[IcpRegistrationState, IcpConfiguration] {
   def name = "ICP"
   // possibility to override the update function, or just use the base class method?
@@ -99,8 +99,8 @@ class IcpRegistration(
   }
 
   override def initializeState(
-      general: GeneralRegistrationState,
-      config: IcpConfiguration
+    general: GeneralRegistrationState,
+    config: IcpConfiguration
   ): IcpRegistrationState = {
     IcpRegistrationState(general, config)
   }
diff --git a/src/main/scala/gingr/api/registration/config/Template.scala b/src/main/scala/gingr/api/registration/config/Template.scala
index 3858545..f2e941e 100644
--- a/src/main/scala/gingr/api/registration/config/Template.scala
+++ b/src/main/scala/gingr/api/registration/config/Template.scala
@@ -22,11 +22,10 @@ import scalismo.common.PointId
 import scalismo.statisticalmodel.MultivariateNormalDistribution
 
 case class TemplateConfiguration(
-    override val maxIterations: Int = 1,
-    override val threshold: Double = 1e-5,
-    override val converged: (GeneralRegistrationState, GeneralRegistrationState, Double) => Boolean =
-      (_, _, _) => false,
-    override val useLandmarkCorrespondence: Boolean = true
+  override val maxIterations: Int = 1,
+  override val threshold: Double = 1e-5,
+  override val converged: (GeneralRegistrationState, GeneralRegistrationState, Double) => Boolean = (_, _, _) => false,
+  override val useLandmarkCorrespondence: Boolean = true
 ) extends GingrConfig {}
 
 case class TemplateRegistrationState(general: GeneralRegistrationState, config: TemplateConfiguration)
@@ -44,8 +43,8 @@ object TemplateRegistrationState {
 }
 
 class TemplateRegistration(
-    override val getCorrespondence: TemplateRegistrationState => CorrespondencePairs =
-      (state: TemplateRegistrationState) => CorrespondencePairs.empty()
+  override val getCorrespondence: TemplateRegistrationState => CorrespondencePairs =
+    (state: TemplateRegistrationState) => CorrespondencePairs.empty()
 ) extends GingrAlgorithm[TemplateRegistrationState, TemplateConfiguration] {
   def name = "Template"
   override val getUncertainty: (PointId, TemplateRegistrationState) => MultivariateNormalDistribution =
@@ -53,8 +52,8 @@ class TemplateRegistration(
       MultivariateNormalDistribution(DenseVector.zeros[Double](3), DenseMatrix.eye[Double](3))
 
   override def initializeState(
-      general: GeneralRegistrationState,
-      config: TemplateConfiguration
+    general: GeneralRegistrationState,
+    config: TemplateConfiguration
   ): TemplateRegistrationState = {
     TemplateRegistrationState(general, config)
   }
diff --git a/src/main/scala/gingr/api/registration/utils/ClosestPointRegistrator.scala b/src/main/scala/gingr/api/registration/utils/ClosestPointRegistrator.scala
index 36211cb..56bf9bd 100644
--- a/src/main/scala/gingr/api/registration/utils/ClosestPointRegistrator.scala
+++ b/src/main/scala/gingr/api/registration/utils/ClosestPointRegistrator.scala
@@ -18,7 +18,7 @@
 package gingr.api.registration.utils
 
 import scalismo.common.{DiscreteDomain, PointId, UnstructuredPointsDomain}
-import scalismo.geometry.{EuclideanVector, Point, _3D}
+import scalismo.geometry.{_3D, EuclideanVector, Point}
 import scalismo.mesh.TriangleMesh
 
 trait ClosestPointRegistrator[DDomain[_3D] <: DiscreteDomain[_3D]] {
@@ -27,17 +27,17 @@ trait ClosestPointRegistrator[DDomain[_3D] <: DiscreteDomain[_3D]] {
            Additionally the average closest point distance is returned
    */
   def closestPointCorrespondence(
-      template: DDomain[_3D],
-      target: DDomain[_3D]
+    template: DDomain[_3D],
+    target: DDomain[_3D]
   ): (Seq[(PointId, Point[_3D], Double)], Double)
 
   def closestPointCorrespondenceReversal(
-      template: DDomain[_3D],
-      target: DDomain[_3D]
+    template: DDomain[_3D],
+    target: DDomain[_3D]
   ): (Seq[(PointId, Point[_3D], Double)], Double) = {
     val corr = closestPointCorrespondence(target, template)
     val inverted = corr._1.map { case (id, p, w) =>
-      val templateId  = template.pointSet.findClosestPoint(p).id
+      val templateId = template.pointSet.findClosestPoint(p).id
       val targetPoint = target.pointSet.point(id)
       (templateId, targetPoint, w)
     }
@@ -73,14 +73,14 @@ object NonRigidClosestPointRegistrator {
 
   object ClosestPointTriangleMesh3D extends ClosestPointRegistrator[TriangleMesh] {
     override def closestPointCorrespondence(
-        template: TriangleMesh[_3D],
-        target: TriangleMesh[_3D]
+      template: TriangleMesh[_3D],
+      target: TriangleMesh[_3D]
     ): (Seq[(PointId, Point[_3D], Double)], Double) = {
       var distance = 0.0
       val corr = template.pointSet.pointIds.toSeq.map { id =>
-        val p                     = template.pointSet.point(id)
+        val p = template.pointSet.point(id)
         val closestPointOnSurface = target.operations.closestPointOnSurface(p)
-        val closestPoint          = target.pointSet.findClosestPoint(closestPointOnSurface.point)
+        val closestPoint = target.pointSet.findClosestPoint(closestPointOnSurface.point)
         val w =
           if (isPointOnBoundary(closestPoint.id, target)) 0.0
           else if (
@@ -97,8 +97,8 @@ object NonRigidClosestPointRegistrator {
 
   object ClosestPointAlongNormalTriangleMesh3D extends ClosestPointRegistrator[TriangleMesh] {
     override def closestPointCorrespondence(
-        template: TriangleMesh[_3D],
-        target: TriangleMesh[_3D]
+      template: TriangleMesh[_3D],
+      target: TriangleMesh[_3D]
     ): (Seq[(PointId, Point[_3D], Double)], Double) = {
       var distance = 0.0
       val corr = template.pointSet.pointIds.toSeq.map { id =>
@@ -132,14 +132,14 @@ object NonRigidClosestPointRegistrator {
 
   object ClosestPointUnstructuredPointsDomain3D extends ClosestPointRegistrator[UnstructuredPointsDomain] {
     override def closestPointCorrespondence(
-        template: UnstructuredPointsDomain[_3D],
-        target: UnstructuredPointsDomain[_3D]
+      template: UnstructuredPointsDomain[_3D],
+      target: UnstructuredPointsDomain[_3D]
     ): (Seq[(PointId, Point[_3D], Double)], Double) = {
       var distance = 0.0
       val corr = template.pointSet.pointIds.toSeq.map { id =>
-        val p            = template.pointSet.point(id)
+        val p = template.pointSet.point(id)
         val closestPoint = target.pointSet.findClosestPoint(p)
-        val w            = 1.0
+        val w = 1.0
         distance += (p - closestPoint.point).norm
         (id, closestPoint.point, w)
       }
@@ -149,8 +149,8 @@ object NonRigidClosestPointRegistrator {
 
   object ClosestPointTriangleMesh3DSimple extends ClosestPointRegistrator[TriangleMesh] {
     override def closestPointCorrespondence(
-        template: TriangleMesh[_3D],
-        target: TriangleMesh[_3D]
+      template: TriangleMesh[_3D],
+      target: TriangleMesh[_3D]
     ): (Seq[(PointId, Point[_3D], Double)], Double) = {
       ClosestPointUnstructuredPointsDomain3D.closestPointCorrespondence(
         UnstructuredPointsDomain(template.pointSet),
diff --git a/src/main/scala/gingr/api/registration/utils/DiscreteDomainConverter.scala b/src/main/scala/gingr/api/registration/utils/DiscreteDomainConverter.scala
index 7337c75..b1a9490 100644
--- a/src/main/scala/gingr/api/registration/utils/DiscreteDomainConverter.scala
+++ b/src/main/scala/gingr/api/registration/utils/DiscreteDomainConverter.scala
@@ -25,17 +25,17 @@ import scalismo.common.{
   UnstructuredPointsDomain3D,
   Vectorizer
 }
-import scalismo.geometry.{Point, _3D}
+import scalismo.geometry.{_3D, Point}
 import scalismo.mesh.{TetrahedralMesh, TetrahedralMesh3D, TriangleMesh, TriangleMesh3D}
 
 trait DiscreteDomainConverter[DDomain[_3D] <: DiscreteDomain[_3D]] {
   def denseMatrixToDomain(
-      mat: DenseMatrix[Double],
-      reference: DDomain[_3D]
+    mat: DenseMatrix[Double],
+    reference: DDomain[_3D]
   ): DDomain[_3D]
 
   def toMatrix(dc: DDomain[_3D])(implicit
-      vectorizer: Vectorizer[Point[_3D]]
+    vectorizer: Vectorizer[Point[_3D]]
   ): DenseMatrix[Double] = {
     PointSequenceConverter.toMatrix(dc.pointSet.points.toSeq)
   }
@@ -45,8 +45,8 @@ object DiscreteDomainConverter {
 
   implicit object denseMatrixToPointDomain3D extends DiscreteDomainConverter[UnstructuredPointsDomain] {
     override def denseMatrixToDomain(
-        mat: DenseMatrix[Double],
-        reference: UnstructuredPointsDomain[_3D]
+      mat: DenseMatrix[Double],
+      reference: UnstructuredPointsDomain[_3D]
     ): UnstructuredPointsDomain[_3D] = {
       UnstructuredPointsDomain3D(PointSequenceConverter.toPointSequence(mat).toIndexedSeq)
     }
@@ -63,8 +63,8 @@ object DiscreteDomainConverter {
 
   implicit object denseMatrixToTetrahedralMesh3D extends DiscreteDomainConverter[TetrahedralMesh] {
     override def denseMatrixToDomain(
-        mat: DenseMatrix[Double],
-        reference: TetrahedralMesh[_3D]
+      mat: DenseMatrix[Double],
+      reference: TetrahedralMesh[_3D]
     ): TetrahedralMesh[_3D] = {
       TetrahedralMesh3D(PointSequenceConverter.toPointSequence(mat).toIndexedSeq, reference.tetrahedralization)
     }
diff --git a/src/main/scala/gingr/api/registration/utils/GPMMHelper.scala b/src/main/scala/gingr/api/registration/utils/GPMMHelper.scala
index 0b9a6d8..81ae36b 100644
--- a/src/main/scala/gingr/api/registration/utils/GPMMHelper.scala
+++ b/src/main/scala/gingr/api/registration/utils/GPMMHelper.scala
@@ -19,7 +19,7 @@ package gingr.api.registration.utils
 
 import scalismo.common.UnstructuredPoints
 import scalismo.common.interpolation.TriangleMeshInterpolator3D
-import scalismo.geometry.{EuclideanVector, _3D}
+import scalismo.geometry.{_3D, EuclideanVector}
 import scalismo.kernels.{DiagonalKernel, GaussianKernel}
 import scalismo.mesh.TriangleMesh
 import scalismo.statisticalmodel.{GaussianProcess, LowRankGaussianProcess, PointDistributionModel}
@@ -37,18 +37,18 @@ object GPMMHelper {
   }
 
   def automaticGPMMfromTemplate(
-      template: TriangleMesh[_3D],
-      relativeTolerance: Double = 0.1
+    template: TriangleMesh[_3D],
+    relativeTolerance: Double = 0.1
   ): PointDistributionModel[_3D, TriangleMesh] = {
     println("Constructing GPMM from template")
     val maxDist = maximumPointDistance(template.pointSet)
     val minDist = minimumPointDistance(template.pointSet)
-    val sigma1  = maxDist / 4
-    val sigma2  = maxDist / 8
-    val sigma3  = minDist * 5
-    val scale1  = sigma1 / 2
-    val scale2  = sigma2 / 2
-    val scale3  = sigma3 / 2
+    val sigma1 = maxDist / 4
+    val sigma2 = maxDist / 8
+    val sigma3 = minDist * 5
+    val scale1 = sigma1 / 2
+    val scale2 = sigma2 / 2
+    val scale3 = sigma3 / 2
 
     println(s"Maximum distance: ${maxDist}, minimum distance: ${minDist}")
     val k = DiagonalKernel(GaussianKernel[_3D](sigma1) * scale1, 3) +
diff --git a/src/main/scala/gingr/api/registration/utils/PointSequenceConverter.scala b/src/main/scala/gingr/api/registration/utils/PointSequenceConverter.scala
index c3d3cd7..6210020 100644
--- a/src/main/scala/gingr/api/registration/utils/PointSequenceConverter.scala
+++ b/src/main/scala/gingr/api/registration/utils/PointSequenceConverter.scala
@@ -20,7 +20,7 @@ package gingr.api.registration.utils
 import breeze.linalg.{CSCMatrix, DenseMatrix, DenseVector}
 import scalismo.common.{PointId, Vectorizer}
 import scalismo.geometry.Point.Point3DVectorizer
-import scalismo.geometry.{Point, Point3D, _3D}
+import scalismo.geometry.{_3D, Point, Point3D}
 
 import scala.collection.parallel.CollectionConverters.ImmutableSeqIsParallelizable
 
@@ -29,7 +29,7 @@ object PointSequenceConverter {
 
   def toMatrix(points: Seq[Point[_3D]])(implicit vectorizer: Vectorizer[Point[_3D]]): DenseMatrix[Double] = {
     val dim: Int = vectorizer.dim
-    val m        = DenseMatrix.zeros[Double](points.length, dim)
+    val m = DenseMatrix.zeros[Double](points.length, dim)
     points.zipWithIndex.par.foreach { case (p, i) =>
       m(i, ::) := vectorizer.vectorize(p).t
     }
@@ -37,11 +37,11 @@ object PointSequenceConverter {
   }
 
   def toLMMatrix(points: Seq[Point[_3D]], pointIds: Seq[PointId], n: Int)(implicit
-      vectorizer: Vectorizer[Point[_3D]]
+    vectorizer: Vectorizer[Point[_3D]]
   ): CSCMatrix[Double] = {
-    val nPoints  = points.size
+    val nPoints = points.size
     val dim: Int = vectorizer.dim
-    val csc      = CSCMatrix.zeros[Double](nPoints, dim)
+    val csc = CSCMatrix.zeros[Double](nPoints, dim)
     points.zip(pointIds).zipWithIndex.foreach { case ((p, pid), i) =>
       val v = vectorizer.vectorize(p)
       for (j <- 0 until dim) {
@@ -75,7 +75,7 @@ object PointSequenceConverter {
   }
 
   def toPointSequence(
-      v: DenseVector[Double]
+    v: DenseVector[Double]
   )(implicit vectorizer: Vectorizer[Point[_3D]]): Seq[Point[_3D]] = {
     vectorTo3Dpoints(v)
   }
diff --git a/src/main/scala/gingr/api/registration/utils/SimilarityTransformParameters.scala b/src/main/scala/gingr/api/registration/utils/SimilarityTransformParameters.scala
index 7172b07..015550c 100644
--- a/src/main/scala/gingr/api/registration/utils/SimilarityTransformParameters.scala
+++ b/src/main/scala/gingr/api/registration/utils/SimilarityTransformParameters.scala
@@ -17,7 +17,7 @@
 
 package gingr.api.registration.utils
 
-import scalismo.geometry.{Point, _3D}
+import scalismo.geometry.{_3D, Point}
 import scalismo.transformations.{
   Rotation,
   Scaling,
@@ -28,9 +28,9 @@ import scalismo.transformations.{
 
 case class SimilarityTransformParameters(s: Scaling[_3D], t: Translation[_3D], R: Rotation[_3D]) {
   val simTransform: TranslationAfterScalingAfterRotation[_3D] = TranslationAfterScalingAfterRotation(t, s, R)
-  val rigidTransform: TranslationAfterRotation[_3D]           = TranslationAfterRotation(t, R)
-  val invSimTransform                                         = simTransform.inverse
-  val invRigidTransform                                       = rigidTransform.inverse
-  def transform(points: Seq[Point[_3D]]): Seq[Point[_3D]]     = points.map(simTransform)
-  def invTransform(points: Seq[Point[_3D]]): Seq[Point[_3D]]  = points.map(invSimTransform)
+  val rigidTransform: TranslationAfterRotation[_3D] = TranslationAfterRotation(t, R)
+  val invSimTransform = simTransform.inverse
+  val invRigidTransform = rigidTransform.inverse
+  def transform(points: Seq[Point[_3D]]): Seq[Point[_3D]] = points.map(simTransform)
+  def invTransform(points: Seq[Point[_3D]]): Seq[Point[_3D]] = points.map(invSimTransform)
 }
diff --git a/src/main/scala/gingr/api/sampling/Evaluator.scala b/src/main/scala/gingr/api/sampling/Evaluator.scala
index 7b4a992..da6cd2b 100644
--- a/src/main/scala/gingr/api/sampling/Evaluator.scala
+++ b/src/main/scala/gingr/api/sampling/Evaluator.scala
@@ -39,10 +39,10 @@ case class AcceptAll[State <: GingrRegistrationState[State]]() extends Evaluator
 }
 
 case class IndependentPoints[State <: GingrRegistrationState[State]](
-    state: State,
-    uncertainty: Double,
-    mode: EvaluationMode = ModelToTargetEvaluation,
-    evaluatedPoints: Option[Int] = None
+  state: State,
+  uncertainty: Double,
+  mode: EvaluationMode = ModelToTargetEvaluation,
+  evaluatedPoints: Option[Int] = None
 ) extends Evaluator[State] {
   private val likelihoodModel = breeze.stats.distributions.Gaussian(0, uncertainty)
 
diff --git a/src/main/scala/gingr/api/sampling/Generator.scala b/src/main/scala/gingr/api/sampling/Generator.scala
index b37459d..04368d0 100644
--- a/src/main/scala/gingr/api/sampling/Generator.scala
+++ b/src/main/scala/gingr/api/sampling/Generator.scala
@@ -26,10 +26,10 @@ import scalismo.utils.Random
 class Generator[State <: GingrRegistrationState[State]](implicit rnd: Random) {
 
   val defaultTranslation = 0.1
-  val defaultRotation    = 0.01
+  val defaultRotation = 0.01
 
   def RandomShape(
-      steps: Seq[Double] = Seq(1.0, 0.1, 0.01)
+    steps: Seq[Double] = Seq(1.0, 0.1, 0.01)
   ): ProposalGenerator[State] with TransitionProbability[State] = {
     val gen = steps.map { d =>
       (1.0 / steps.length.toDouble, RandomShapeUpdateProposal[State](d, generatedBy = s"RandomShape-${d}"))
@@ -38,9 +38,9 @@ class Generator[State <: GingrRegistrationState[State]](implicit rnd: Random) {
   }
 
   def RandomRotation(
-      rotYaw: Double = defaultRotation,
-      rotPitch: Double = defaultRotation,
-      rotRoll: Double = defaultRotation
+    rotYaw: Double = defaultRotation,
+    rotPitch: Double = defaultRotation,
+    rotRoll: Double = defaultRotation
   ): ProposalGenerator[State] with TransitionProbability[State] = {
     MixtureProposal(
       0.5 *: GaussianAxisRotationProposal[State](rotYaw, YawAxis, generatedBy = s"RotationYaw-${rotYaw}") +
@@ -50,9 +50,9 @@ class Generator[State <: GingrRegistrationState[State]](implicit rnd: Random) {
   }
 
   def RandomTranslation(
-      transX: Double = defaultTranslation,
-      transY: Double = defaultTranslation,
-      transZ: Double = defaultTranslation
+    transX: Double = defaultTranslation,
+    transY: Double = defaultTranslation,
+    transZ: Double = defaultTranslation
   ): ProposalGenerator[State] with TransitionProbability[State] = {
     MixtureProposal(
       0.5 *: GaussianAxisTranslationProposal[State](transX, 0, generatedBy = s"TranslationX-${transX}") +
@@ -62,12 +62,12 @@ class Generator[State <: GingrRegistrationState[State]](implicit rnd: Random) {
   }
 
   def RandomPose(
-      rotYaw: Double = defaultRotation,
-      rotPitch: Double = defaultRotation,
-      rotRoll: Double = defaultRotation,
-      transX: Double = defaultTranslation,
-      transY: Double = defaultTranslation,
-      transZ: Double = defaultTranslation
+    rotYaw: Double = defaultRotation,
+    rotPitch: Double = defaultRotation,
+    rotRoll: Double = defaultRotation,
+    transX: Double = defaultTranslation,
+    transY: Double = defaultTranslation,
+    transZ: Double = defaultTranslation
   ): ProposalGenerator[State] with TransitionProbability[State] = {
     MixtureProposal(
       0.5 *: RandomRotation(rotYaw, rotPitch, rotRoll) +
diff --git a/src/main/scala/gingr/api/sampling/evaluators/EvaluationCaching.scala b/src/main/scala/gingr/api/sampling/evaluators/EvaluationCaching.scala
index 7a82a07..9d245fb 100644
--- a/src/main/scala/gingr/api/sampling/evaluators/EvaluationCaching.scala
+++ b/src/main/scala/gingr/api/sampling/evaluators/EvaluationCaching.scala
@@ -21,8 +21,9 @@ import gingr.api.GingrRegistrationState
 import scalismo.sampling.DistributionEvaluator
 import scalismo.utils.Memoize
 
-/** This trait can be mixe in with an DistributionEvalutor, to enable caching of the values
-  */
+/**
+ * This trait can be mixe in with an DistributionEvalutor, to enable caching of the values
+ */
 trait EvaluationCaching[State <: GingrRegistrationState[State]] {
   self: DistributionEvaluator[State] =>
 
diff --git a/src/main/scala/gingr/api/sampling/evaluators/IndependentPointDistanceEvaluator.scala b/src/main/scala/gingr/api/sampling/evaluators/IndependentPointDistanceEvaluator.scala
index ed4ef73..f2cf3f0 100644
--- a/src/main/scala/gingr/api/sampling/evaluators/IndependentPointDistanceEvaluator.scala
+++ b/src/main/scala/gingr/api/sampling/evaluators/IndependentPointDistanceEvaluator.scala
@@ -20,7 +20,7 @@ package gingr.api.sampling.evaluators
 import breeze.stats.distributions.ContinuousDistr
 import gingr.api.GingrRegistrationState
 import scalismo.common.PointId
-import scalismo.geometry.{Point, _3D}
+import scalismo.geometry.{_3D, Point}
 import scalismo.mesh.TriangleMesh3D
 import scalismo.sampling.DistributionEvaluator
 
@@ -33,15 +33,15 @@ case object TargetToModelEvaluation extends EvaluationMode
 case object SymmetricEvaluation extends EvaluationMode
 
 case class IndependentPointDistanceEvaluator[State <: GingrRegistrationState[State]](
-    sample: State,
-    likelihoodModel: ContinuousDistr[Double],
-    evaluationMode: EvaluationMode,
-    numberOfPointsForComparison: Option[Int]
+  sample: State,
+  likelihoodModel: ContinuousDistr[Double],
+  evaluationMode: EvaluationMode,
+  numberOfPointsForComparison: Option[Int]
 ) extends DistributionEvaluator[State]
     with EvaluationCaching[State] {
 
   private val instance = sample.general.fit
-  private val target   = sample.general.target
+  private val target = sample.general.target
 
   private val (instanceDecimated, targetDecimated) = numberOfPointsForComparison match {
     case Some(num) => (instance.operations.decimate(num), target.operations.decimate(num))
@@ -49,7 +49,7 @@ case class IndependentPointDistanceEvaluator[State <: GingrRegistrationState[Sta
   }
 
   private val randomPointsOnTarget: IndexedSeq[Point[_3D]] = targetDecimated.pointSet.points.toIndexedSeq
-  private val randomPointIdsOnModel: IndexedSeq[PointId]   = instanceDecimated.pointSet.pointIds.toIndexedSeq
+  private val randomPointIdsOnModel: IndexedSeq[PointId] = instanceDecimated.pointSet.pointIds.toIndexedSeq
 
   def distModelToTarget(modelSample: TriangleMesh3D): Double = {
     val pointsOnSample = randomPointIdsOnModel.map(modelSample.pointSet.point)
diff --git a/src/main/scala/gingr/api/sampling/evaluators/ModelEvaluator.scala b/src/main/scala/gingr/api/sampling/evaluators/ModelEvaluator.scala
index a06302c..6114f78 100644
--- a/src/main/scala/gingr/api/sampling/evaluators/ModelEvaluator.scala
+++ b/src/main/scala/gingr/api/sampling/evaluators/ModelEvaluator.scala
@@ -17,7 +17,7 @@
 
 package gingr.api.sampling.evaluators
 
-import breeze.linalg.{DenseVector, diag}
+import breeze.linalg.{diag, DenseVector}
 import gingr.api.GingrRegistrationState
 import scalismo.sampling.DistributionEvaluator
 import scalismo.statisticalmodel.MultivariateNormalDistribution
diff --git a/src/main/scala/gingr/api/sampling/generators/GeneratorWrapperDeterministic.scala b/src/main/scala/gingr/api/sampling/generators/GeneratorWrapperDeterministic.scala
index dc542f6..a9bd8a4 100644
--- a/src/main/scala/gingr/api/sampling/generators/GeneratorWrapperDeterministic.scala
+++ b/src/main/scala/gingr/api/sampling/generators/GeneratorWrapperDeterministic.scala
@@ -21,8 +21,8 @@ import gingr.api.GingrRegistrationState
 import scalismo.utils.Random
 
 case class GeneratorWrapperDeterministic[State <: GingrRegistrationState[State]](
-    update: (State, Boolean) => State,
-    generatedBy: String = "Deterministic"
+  update: (State, Boolean) => State,
+  generatedBy: String = "Deterministic"
 )(implicit rnd: Random)
     extends GingrGeneratorWrapper[State] {
   override def gingrPropose(current: State): State = {
diff --git a/src/main/scala/gingr/api/sampling/generators/GeneratorWrapperStochastic.scala b/src/main/scala/gingr/api/sampling/generators/GeneratorWrapperStochastic.scala
index afa1f7a..c1a9f04 100644
--- a/src/main/scala/gingr/api/sampling/generators/GeneratorWrapperStochastic.scala
+++ b/src/main/scala/gingr/api/sampling/generators/GeneratorWrapperStochastic.scala
@@ -26,9 +26,9 @@ import scalismo.utils.{Memoize, Random}
 import scala.util.Try
 
 case class GeneratorWrapperStochastic[State <: GingrRegistrationState[State]](
-    update: (State, Boolean) => State,
-    cashedPosterior: Memoize[State, Try[PointDistributionModel[_3D, TriangleMesh]]],
-    generatedBy: String = "InformedProposal"
+  update: (State, Boolean) => State,
+  cashedPosterior: Memoize[State, Try[PointDistributionModel[_3D, TriangleMesh]]],
+  generatedBy: String = "InformedProposal"
 )(implicit rnd: Random)
     extends GingrGeneratorWrapper[State] {
   override def gingrPropose(current: State): State = {
@@ -51,7 +51,7 @@ case class GeneratorWrapperStochastic[State <: GingrRegistrationState[State]](
       } else from.general.fit
       try {
         val projectedTo = posterior.get.coefficients(toMesh)
-        val logpdf      = posterior.get.gp.logpdf(projectedTo)
+        val logpdf = posterior.get.gp.logpdf(projectedTo)
         logpdf
       } catch {
         case _: Throwable => {
diff --git a/src/main/scala/gingr/api/sampling/generators/RandomPoseUpdateProposal.scala b/src/main/scala/gingr/api/sampling/generators/RandomPoseUpdateProposal.scala
index 9813f84..7849db2 100644
--- a/src/main/scala/gingr/api/sampling/generators/RandomPoseUpdateProposal.scala
+++ b/src/main/scala/gingr/api/sampling/generators/RandomPoseUpdateProposal.scala
@@ -28,9 +28,9 @@ case object PitchAxis extends RotationAxis
 case object YawAxis extends RotationAxis
 
 case class GaussianAxisRotationProposal[State <: GingrRegistrationState[State]](
-    sdevRot: Double,
-    axis: RotationAxis,
-    generatedBy: String = "RotationProposal"
+  sdevRot: Double,
+  axis: RotationAxis,
+  generatedBy: String = "RotationProposal"
 ) extends GingrGeneratorWrapper[State] {
 
   val perturbationDistr = new breeze.stats.distributions.Gaussian(0, sdevRot)
@@ -59,7 +59,7 @@ case class GaussianAxisRotationProposal[State <: GingrRegistrationState[State]](
       Double.NegativeInfinity
     } else {
       val rotParamsFrom = from.general.modelParameters.pose.rotation.angles
-      val rotParamsTo   = to.general.modelParameters.pose.rotation.angles
+      val rotParamsTo = to.general.modelParameters.pose.rotation.angles
       val residual = axis match {
         case RollAxis  => rotParamsTo.phi - rotParamsFrom.phi
         case PitchAxis => rotParamsTo.theta - rotParamsFrom.theta
@@ -71,9 +71,9 @@ case class GaussianAxisRotationProposal[State <: GingrRegistrationState[State]](
 }
 
 case class GaussianAxisTranslationProposal[State <: GingrRegistrationState[State]](
-    sdevTrans: Double,
-    axis: Int,
-    generatedBy: String = "TranslationProposal"
+  sdevTrans: Double,
+  axis: Int,
+  generatedBy: String = "TranslationProposal"
 ) extends GingrGeneratorWrapper[State] {
   require(axis < 3)
   val perturbationDistr = new breeze.stats.distributions.Gaussian(0, sdevTrans)
diff --git a/src/main/scala/gingr/api/sampling/generators/RandomShapeUpdateProposal.scala b/src/main/scala/gingr/api/sampling/generators/RandomShapeUpdateProposal.scala
index d873fb9..342c880 100644
--- a/src/main/scala/gingr/api/sampling/generators/RandomShapeUpdateProposal.scala
+++ b/src/main/scala/gingr/api/sampling/generators/RandomShapeUpdateProposal.scala
@@ -21,8 +21,8 @@ import gingr.api.{GingrRegistrationState, ShapeParameters}
 import scalismo.utils.Random
 
 case class RandomShapeUpdateProposal[State <: GingrRegistrationState[State]](
-    stdev: Double,
-    generatedBy: String = "RandomShapeUpdateProposal"
+  stdev: Double,
+  generatedBy: String = "RandomShapeUpdateProposal"
 )(implicit random: Random)
     extends GingrGeneratorWrapper[State] {
 
diff --git a/src/main/scala/gingr/api/sampling/loggers/BestAndCurrentSampleLogger.scala b/src/main/scala/gingr/api/sampling/loggers/BestAndCurrentSampleLogger.scala
index ac4b7d8..dd37928 100644
--- a/src/main/scala/gingr/api/sampling/loggers/BestAndCurrentSampleLogger.scala
+++ b/src/main/scala/gingr/api/sampling/loggers/BestAndCurrentSampleLogger.scala
@@ -26,7 +26,7 @@ class BestAndCurrentSampleLogger[A](evaluator: DistributionEvaluator[A]) extends
   private case class EvaluatedSample(sample: A, value: Double)
 
   private var bestState: Option[EvaluatedSample] = None
-  private var currentState: Option[A]            = None
+  private var currentState: Option[A] = None
 
   override def logState(sample: A): Unit = {
     val value = evaluator.logValue(sample)
diff --git a/src/main/scala/gingr/api/sampling/loggers/EmptyChainStateLogger.scala b/src/main/scala/gingr/api/sampling/loggers/EmptyChainStateLogger.scala
index 0e721c0..54f207f 100644
--- a/src/main/scala/gingr/api/sampling/loggers/EmptyChainStateLogger.scala
+++ b/src/main/scala/gingr/api/sampling/loggers/EmptyChainStateLogger.scala
@@ -27,16 +27,16 @@ case class EmptyChainStateLogger[State <: GingrRegistrationState[State]]() exten
 
 case class EmptyAcceptRejectLogger[State <: GingrRegistrationState[State]]() extends AcceptRejectLogger[State] {
   override def accept(
-      current: State,
-      sample: State,
-      generator: ProposalGenerator[State],
-      evaluator: DistributionEvaluator[State]
+    current: State,
+    sample: State,
+    generator: ProposalGenerator[State],
+    evaluator: DistributionEvaluator[State]
   ): Unit = {}
 
   override def reject(
-      current: State,
-      sample: State,
-      generator: ProposalGenerator[State],
-      evaluator: DistributionEvaluator[State]
+    current: State,
+    sample: State,
+    generator: ProposalGenerator[State],
+    evaluator: DistributionEvaluator[State]
   ): Unit = {}
 }
diff --git a/src/main/scala/gingr/api/sampling/loggers/JSONStateLogger.scala b/src/main/scala/gingr/api/sampling/loggers/JSONStateLogger.scala
index 29c4bd2..ef43707 100644
--- a/src/main/scala/gingr/api/sampling/loggers/JSONStateLogger.scala
+++ b/src/main/scala/gingr/api/sampling/loggers/JSONStateLogger.scala
@@ -34,16 +34,16 @@ import scala.collection.mutable.ListBuffer
 import scala.io.Source
 
 case class jsonLogFormat(
-    index: Int,
-    name: String,
-    logvalue: Map[String, Double],
-    status: Boolean,
-    modelParameters: Seq[Double],
-    translation: Seq[Double],
-    rotation: Seq[Double],
-    rotationCenter: Seq[Double],
-    scaling: Double,
-    datetime: String
+  index: Int,
+  name: String,
+  logvalue: Map[String, Double],
+  status: Boolean,
+  modelParameters: Seq[Double],
+  translation: Seq[Double],
+  rotation: Seq[Double],
+  rotationCenter: Seq[Double],
+  scaling: Double,
+  datetime: String
 )
 
 object JsonLoggerProtocol {
@@ -51,16 +51,16 @@ object JsonLoggerProtocol {
 }
 
 case class JSONStateLogger[State <: GingrRegistrationState[State]](
-    evaluators: Evaluator[State],
-    filePath: Option[File] = None
+  evaluators: Evaluator[State],
+  filePath: Option[File] = None
 ) extends AcceptRejectLogger[State] {
   import JsonLoggerProtocol._
 
-  private var numOfRejected: Int             = 0
-  private var numOfAccepted: Int             = 0
+  private var numOfRejected: Int = 0
+  private var numOfAccepted: Int = 0
   private var generatedBy: SortedSet[String] = SortedSet()
-  private val datetimeFormat                 = new SimpleDateFormat("yyyy-MM-dd HH:mm:ss")
-  private val productEvaluator               = evaluators.productEvaluator()
+  private val datetimeFormat = new SimpleDateFormat("yyyy-MM-dd HH:mm:ss")
+  private val productEvaluator = evaluators.productEvaluator()
 
   filePath.map { f =>
     if (f.getParentFile != null) {
@@ -99,10 +99,10 @@ case class JSONStateLogger[State <: GingrRegistrationState[State]](
   }
 
   override def accept(
-      current: State,
-      sample: State,
-      generator: ProposalGenerator[State],
-      evaluator: DistributionEvaluator[State]
+    current: State,
+    sample: State,
+    generator: ProposalGenerator[State],
+    evaluator: DistributionEvaluator[State]
   ): Unit = {
     val evalValue = mapEvaluators(sample)
     log += jsonLogFormat(
@@ -122,10 +122,10 @@ case class JSONStateLogger[State <: GingrRegistrationState[State]](
 
   // The rejected state will contain the same parameters as the previous accepted state, so no need to double store all the information
   override def reject(
-      current: State,
-      sample: State,
-      generator: ProposalGenerator[State],
-      evaluator: DistributionEvaluator[State]
+    current: State,
+    sample: State,
+    generator: ProposalGenerator[State],
+    evaluator: DistributionEvaluator[State]
   ): Unit = {
     val evalValue = mapEvaluators(sample)
     log += jsonLogFormat(
@@ -208,7 +208,7 @@ object JSONStateLogger {
 
   def loadLog(filePath: File): IndexedSeq[jsonLogFormat] = {
     println(s"Loading JSON log file: ${filePath.toString}")
-    val src  = Source.fromFile(filePath.toString)
+    val src = Source.fromFile(filePath.toString)
     val data = src.mkString.parseJson.convertTo[IndexedSeq[jsonLogFormat]]
     src.close()
     data
diff --git a/src/main/scala/gingr/other/algorithms/BCPDRegistration.scala b/src/main/scala/gingr/other/algorithms/BCPDRegistration.scala
index 30ee6a7..3fe9eb3 100644
--- a/src/main/scala/gingr/other/algorithms/BCPDRegistration.scala
+++ b/src/main/scala/gingr/other/algorithms/BCPDRegistration.scala
@@ -19,21 +19,21 @@ package gingr.other.algorithms
 
 import gingr.other.algorithms.cpd.BCPD
 import scalismo.common.{DiscreteDomain, DiscreteField, DomainWarp, Vectorizer}
-import scalismo.geometry.{Point, _3D}
+import scalismo.geometry.{_3D, Point}
 import scalismo.kernels.PDKernel
 
 class BCPDRegistration[DDomain[_3D] <: DiscreteDomain[_3D]](
-    template: DDomain[_3D],
-    target: DDomain[_3D],
-    w: Double = 0,         // Outlier, [0,1]
-    lambda: Double = 2.0,  // Noise scaling, R+
-    gamma: Double = 1.0,   // Initial noise scaling, R+
-    k: Double = 1.0,       // Dirichlet distribution parameter
-    kernel: PDKernel[_3D], // Positive semi-def kernel
-    max_iterations: Int = 100
+  template: DDomain[_3D],
+  target: DDomain[_3D],
+  w: Double = 0, // Outlier, [0,1]
+  lambda: Double = 2.0, // Noise scaling, R+
+  gamma: Double = 1.0, // Initial noise scaling, R+
+  k: Double = 1.0, // Dirichlet distribution parameter
+  kernel: PDKernel[_3D], // Positive semi-def kernel
+  max_iterations: Int = 100
 )(implicit
-    warper: DomainWarp[_3D, DDomain],
-    vectorizer: Vectorizer[Point[_3D]]
+  warper: DomainWarp[_3D, DDomain],
+  vectorizer: Vectorizer[Point[_3D]]
 ) {
   val cpd = new BCPD(template.pointSet.points.toSeq, target.pointSet.points.toSeq, w, lambda, gamma, k, kernel)
 
diff --git a/src/main/scala/gingr/other/algorithms/CPDRegistration.scala b/src/main/scala/gingr/other/algorithms/CPDRegistration.scala
index 0bdcb8d..4fcc8a6 100644
--- a/src/main/scala/gingr/other/algorithms/CPDRegistration.scala
+++ b/src/main/scala/gingr/other/algorithms/CPDRegistration.scala
@@ -19,17 +19,17 @@ package gingr.other.algorithms
 
 import gingr.other.algorithms.cpd.CPDFactory
 import scalismo.common.{DiscreteDomain, DiscreteField, DomainWarp, Vectorizer}
-import scalismo.geometry.{Point, _3D}
+import scalismo.geometry.{_3D, Point}
 
 class RigidCPDRegistration[DDomain[_3D] <: DiscreteDomain[_3D]](
-    template: DDomain[_3D],
-    lambda: Double = 2,
-    beta: Double = 2,
-    w: Double = 0,
-    max_iterations: Int = 100
+  template: DDomain[_3D],
+  lambda: Double = 2,
+  beta: Double = 2,
+  w: Double = 0,
+  max_iterations: Int = 100
 )(implicit
-    warper: DomainWarp[_3D, DDomain],
-    vectorizer: Vectorizer[Point[_3D]]
+  warper: DomainWarp[_3D, DDomain],
+  vectorizer: Vectorizer[Point[_3D]]
 ) {
   val cpd = new CPDFactory(template.pointSet.points.toSeq, lambda, beta, w)
 
@@ -37,7 +37,7 @@ class RigidCPDRegistration[DDomain[_3D] <: DiscreteDomain[_3D]](
 
   def register(target: DDomain[_3D]): DDomain[_3D] = {
     val registrationTask = registrationMethod(target.pointSet.points.toSeq)
-    val registration     = registrationTask.Registration(max_iterations)
+    val registration = registrationTask.Registration(max_iterations)
     val warpField =
       DiscreteField(template, template.pointSet.points.toIndexedSeq.zip(registration).map { case (a, b) => b - a })
     warper.transformWithField(template, warpField)
@@ -45,28 +45,28 @@ class RigidCPDRegistration[DDomain[_3D] <: DiscreteDomain[_3D]](
 }
 
 class NonRigidCPDRegistration[DDomain[_3D] <: DiscreteDomain[_3D]](
-    template: DDomain[_3D],
-    lambda: Double = 2,
-    beta: Double = 2,
-    w: Double = 0,
-    max_iterations: Int = 100
+  template: DDomain[_3D],
+  lambda: Double = 2,
+  beta: Double = 2,
+  w: Double = 0,
+  max_iterations: Int = 100
 )(implicit
-    warper: DomainWarp[_3D, DDomain],
-    vectorizer: Vectorizer[Point[_3D]]
+  warper: DomainWarp[_3D, DDomain],
+  vectorizer: Vectorizer[Point[_3D]]
 ) extends RigidCPDRegistration[DDomain](template, lambda, beta, w, max_iterations) {
 
   override def registrationMethod(targetPoints: Seq[Point[_3D]]) = cpd.registerNonRigidly(targetPoints)
 }
 
 class AffineCPDRegistration[DDomain[_3D] <: DiscreteDomain[_3D]](
-    template: DDomain[_3D],
-    lambda: Double = 2,
-    beta: Double = 2,
-    w: Double = 0,
-    max_iterations: Int = 100
+  template: DDomain[_3D],
+  lambda: Double = 2,
+  beta: Double = 2,
+  w: Double = 0,
+  max_iterations: Int = 100
 )(implicit
-    warper: DomainWarp[_3D, DDomain],
-    vectorizer: Vectorizer[Point[_3D]]
+  warper: DomainWarp[_3D, DDomain],
+  vectorizer: Vectorizer[Point[_3D]]
 ) extends RigidCPDRegistration[DDomain](template, lambda, beta, w, max_iterations) {
 
   override def registrationMethod(targetPoints: Seq[Point[_3D]]) = cpd.registerAffine(targetPoints)
diff --git a/src/main/scala/gingr/other/algorithms/RigidICPRegistration.scala b/src/main/scala/gingr/other/algorithms/RigidICPRegistration.scala
index 61a68d5..885a5f9 100644
--- a/src/main/scala/gingr/other/algorithms/RigidICPRegistration.scala
+++ b/src/main/scala/gingr/other/algorithms/RigidICPRegistration.scala
@@ -20,15 +20,15 @@ package gingr.other.algorithms
 import gingr.other.algorithms.icp.ICPFactory
 import gingr.other.utils.Registrator
 import scalismo.common.{DiscreteDomain, DiscreteField, DomainWarp, UnstructuredPoints, Vectorizer}
-import scalismo.geometry.{Point, _3D}
+import scalismo.geometry.{_3D, Point}
 
 class RigidICPRegistration[DDomain[_3D] <: DiscreteDomain[_3D]](
-    template: DDomain[_3D],
-    max_iterations: Int = 100
+  template: DDomain[_3D],
+  max_iterations: Int = 100
 )(implicit
-    warper: DomainWarp[_3D, DDomain],
-    vectorizer: Vectorizer[Point[_3D]],
-    registration: Registrator
+  warper: DomainWarp[_3D, DDomain],
+  vectorizer: Vectorizer[Point[_3D]],
+  registration: Registrator
 ) {
   val icp = new ICPFactory(UnstructuredPoints(template.pointSet.points.toIndexedSeq))
 
@@ -36,7 +36,7 @@ class RigidICPRegistration[DDomain[_3D] <: DiscreteDomain[_3D]](
 
   def register(target: DDomain[_3D]): DDomain[_3D] = {
     val registrationTask = registrationMethod(UnstructuredPoints(target.pointSet.points.toIndexedSeq))
-    val registration     = registrationTask.Registration(max_iterations)
+    val registration = registrationTask.Registration(max_iterations)
     val warpField = DiscreteField(
       target,
       target.pointSet.points.toIndexedSeq.zip(registration.points.toIndexedSeq).map { case (a, b) => b - a }
diff --git a/src/main/scala/gingr/other/algorithms/cpd/AffineCPD.scala b/src/main/scala/gingr/other/algorithms/cpd/AffineCPD.scala
index 96b4c75..f5e805b 100644
--- a/src/main/scala/gingr/other/algorithms/cpd/AffineCPD.scala
+++ b/src/main/scala/gingr/other/algorithms/cpd/AffineCPD.scala
@@ -17,31 +17,31 @@
 
 package gingr.other.algorithms.cpd
 
-import breeze.linalg.{Axis, DenseMatrix, DenseVector, InjectNumericOps, diag, inv, sum, trace}
+import breeze.linalg.{diag, inv, sum, trace, Axis, DenseMatrix, DenseVector, InjectNumericOps}
 import scalismo.common.Vectorizer
-import scalismo.geometry.{Point, _3D}
+import scalismo.geometry.{_3D, Point}
 
 /*
  Implementation of Point Set Registration: Coherent Point Drift (CPD) - Affine transformation algorithm
  Paper: https://arxiv.org/pdf/0905.2635.pdf
  */
 private[cpd] class AffineCPD(
-    override val targetPoints: Seq[Point[_3D]],
-    override val cpd: CPDFactory
+  override val targetPoints: Seq[Point[_3D]],
+  override val cpd: CPDFactory
 )(implicit
-    vectorizer: Vectorizer[Point[_3D]]
+  vectorizer: Vectorizer[Point[_3D]]
 ) extends RigidCPD(targetPoints, cpd) {
   import cpd._
   override def Maximization(
-      X: DenseMatrix[Double],
-      Y: DenseMatrix[Double],
-      P: DenseMatrix[Double],
-      sigma2: Double
+    X: DenseMatrix[Double],
+    Y: DenseMatrix[Double],
+    P: DenseMatrix[Double],
+    sigma2: Double
   ): (DenseMatrix[Double], Double) = {
     // Update transform
     val P1: DenseVector[Double] = sum(P, Axis._1)
-    val Pt1                     = sum(P, Axis._0)
-    val Np                      = sum(P1)
+    val Pt1 = sum(P, Axis._0)
+    val Np = sum(P1)
 
     val muX = 1.0 / Np * X.t * P.t * DenseVector.ones[Double](P.rows)
     val muY = 1.0 / Np * Y.t * P1
@@ -52,10 +52,10 @@ private[cpd] class AffineCPD(
     val B = Xhat.t * P.t * Yhat * inv(Yhat.t * diag(P1) * Yhat)
     val t = muX - B * muY
 
-    val s1            = trace(Xhat.t * diag(Pt1) * Xhat)
-    val s2            = trace(Xhat.t * P.t * Yhat * B.t)
+    val s1 = trace(Xhat.t * diag(Pt1) * Xhat)
+    val s2 = trace(Xhat.t * P.t * Yhat * B.t)
     val updatedSigma2 = 1 / (Np * dim) * (s1 - s2)
-    val TY            = Y * B.t + DenseVector.ones[Double](M) * t.t
+    val TY = Y * B.t + DenseVector.ones[Double](M) * t.t
 
     (TY, updatedSigma2)
   }
diff --git a/src/main/scala/gingr/other/algorithms/cpd/BCPD.scala b/src/main/scala/gingr/other/algorithms/cpd/BCPD.scala
index fa44b13..d442636 100644
--- a/src/main/scala/gingr/other/algorithms/cpd/BCPD.scala
+++ b/src/main/scala/gingr/other/algorithms/cpd/BCPD.scala
@@ -18,21 +18,21 @@
 package gingr.other.algorithms.cpd
 
 import gingr.api.registration.utils.PointSequenceConverter
-import breeze.linalg.{Axis, DenseMatrix, DenseVector, det, diag, kron, pinv, sum, svd, tile, trace}
+import breeze.linalg.{det, diag, kron, pinv, sum, svd, tile, trace, Axis, DenseMatrix, DenseVector}
 import breeze.numerics.{abs, digamma}
 import scalismo.common.Vectorizer
-import scalismo.geometry.{Point, _3D}
+import scalismo.geometry.{_3D, Point}
 import scalismo.kernels.PDKernel
 import scalismo.statisticalmodel.MultivariateNormalDistribution
 
 // Similarity transformation parameters
 case class similarityTransformationParameters(
-    sigma: DenseMatrix[Double],
-    s: Double,
-    R: DenseMatrix[Double],
-    t: DenseVector[Double],
-    sigma2: Double,
-    alpha: DenseVector[Double]
+  sigma: DenseMatrix[Double],
+  s: Double,
+  R: DenseMatrix[Double],
+  t: DenseVector[Double],
+  sigma2: Double,
+  alpha: DenseVector[Double]
 )
 
 /*
@@ -40,42 +40,43 @@ case class similarityTransformationParameters(
  Paper: https://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=8985307
  */
 class BCPD(
-    val templatePoints: Seq[Point[_3D]],
-    val targetPoints: Seq[Point[_3D]],
-    val w: Double,      // Outlier, [0,1]
-    val lambda: Double, // Noise scaling, R+
-    val gamma: Double,  // Initial noise scaling, R+
-    val k: Double,
-    val kernel: PDKernel[_3D] // Positive semi-def kernel
+  val templatePoints: Seq[Point[_3D]],
+  val targetPoints: Seq[Point[_3D]],
+  val w: Double, // Outlier, [0,1]
+  val lambda: Double, // Noise scaling, R+
+  val gamma: Double, // Initial noise scaling, R+
+  val k: Double,
+  val kernel: PDKernel[_3D] // Positive semi-def kernel
 )(implicit
-    val vectorizer: Vectorizer[Point[_3D]]
+  val vectorizer: Vectorizer[Point[_3D]]
 ) {
   require(0.0 <= w && w <= 1.0)
   require(lambda > 0)
   require(gamma > 0)
 
-  val M: Int                          = templatePoints.length // num of reference points
-  val N: Int                          = targetPoints.length   // num of target points
-  val dim: Int                        = vectorizer.dim        // dimension
-  val G: DenseMatrix[Double]          = initializeKernelMatrixG(templatePoints, kernel)
-  val Ginv: DenseMatrix[Double]       = pinv(G)
-  val GinvLambda: DenseMatrix[Double] = lambda * Ginv
-  val X: DenseVector[Double]          = PointSequenceConverter.toVector(targetPoints)
-  val Y: DenseVector[Double]          = PointSequenceConverter.toVector(templatePoints)
+  val M: Int = templatePoints.length // num of reference points
+  val N: Int = targetPoints.length // num of target points
+  val dim: Int = vectorizer.dim // dimension
+  val G: DenseMatrix[Double] = initializeKernelMatrixG(templatePoints, kernel)
+  val Ginv: DenseMatrix[Double] = pinv(G)
+  val GinvLambda: DenseMatrix[Double] = Ginv.*(lambda)
+  val X: DenseVector[Double] = PointSequenceConverter.toVector(targetPoints)
+  val Y: DenseVector[Double] = PointSequenceConverter.toVector(templatePoints)
 
   // Helper unit matrices
   val D1Vec: DenseMatrix[Double] = DenseVector.ones[Double](dim).toDenseMatrix
-  val Dmat: DenseMatrix[Double]  = DenseMatrix.eye[Double](dim)
-
-  /** Initialize G matrix - formula in paper fig. 4
-    *
-    * @param points
-    *   // * @param kernel
-    * @return
-    */
+  val Dmat: DenseMatrix[Double] = DenseMatrix.eye[Double](dim)
+
+  /**
+   * Initialize G matrix - formula in paper fig. 4
+   *
+   * @param points
+   *   // * @param kernel
+   * @return
+   */
   private def initializeKernelMatrixG(
-      points: Seq[Point[_3D]],
-      kernel: PDKernel[_3D]
+    points: Seq[Point[_3D]],
+    kernel: PDKernel[_3D]
   ): DenseMatrix[Double] = {
     val G: DenseMatrix[Double] = DenseMatrix.zeros[Double](M, M)
     (0 until M).map { i =>
@@ -86,15 +87,16 @@ class BCPD(
     G
   }
 
-  /** Initialize sigma2 - formula in paper fig. 4
-    *
-    * @param Y
-    * @param X
-    * @return
-    */
+  /**
+   * Initialize sigma2 - formula in paper fig. 4
+   *
+   * @param Y
+   * @param X
+   * @return
+   */
   private def computeSigma2init(Y: Seq[Point[_3D]], X: Seq[Point[_3D]]): Double = {
-    val M   = Y.length
-    val N   = X.length
+    val M = Y.length
+    val N = X.length
     val dim = vectorizer.dim
     val s: Double = (0 until M).flatMap { m =>
       (0 until N).map { n =>
@@ -119,39 +121,47 @@ class BCPD(
       alpha = DenseVector.ones[Double](M) / M.toDouble
     )
 
-    val fit = (0 until max_iteration).foldLeft((templatePoints, simParsInit)) { (it, i) =>
-      val currentPars = it._2
+    var fit = (templatePoints, simParsInit)
+    var i = 0
+    var converged = false
+
+    while (i < max_iteration && !converged) {
+      val currentPars = fit._2
       println(s"CPD, iteration: ${i}, variance: ${currentPars.sigma2}")
-      val iter    = Iteration(it._1, it._2)
-      val TY      = iter._1
+      val iter = Iteration(fit._1, fit._2)
+      val TY = iter._1
       val newPars = iter._2
-      val diff    = abs(newPars.sigma2 - currentPars.sigma2)
+      val diff = math.abs(newPars.sigma2 - currentPars.sigma2)
       if (diff < tolerance) {
         println("Converged")
         printTransformation(iter._2)
-        return TY
+        converged = true
       } else {
-        iter
+        fit = iter
+        i += 1
       }
     }
+
     printTransformation(fit._2)
     fit._1
   }
 
   // TranslationAfterScalingAfterRotation
   private def vectorTransform(v: DenseVector[Double], pars: similarityTransformationParameters): DenseVector[Double] = {
-    pars.s * (kron(DenseMatrix.eye[Double](M), pars.R) * v) +
+    val tmp = (kron(DenseMatrix.eye[Double](M), pars.R) * v) +
       (kron(DenseVector.ones[Double](M).toDenseMatrix, pars.t.toDenseMatrix).toDenseVector)
+    tmp.map(p => p * pars.s)
   }
 
   // RotationAfterScalingAfterTranslation
   private def vectorInvTransform(
-      v: DenseVector[Double],
-      pars: similarityTransformationParameters
+    v: DenseVector[Double],
+    pars: similarityTransformationParameters
   ): DenseVector[Double] = {
     val sinv = 1.0 / pars.s
+    val tmp = (v + (kron(DenseVector.ones[Double](M).toDenseMatrix, pars.t.toDenseMatrix * (-1.0)).toDenseVector))
     (kron(DenseMatrix.eye[Double](M), pinv(pars.R))) *
-      (sinv * (v + (kron(DenseVector.ones[Double](M).toDenseMatrix, pars.t.toDenseMatrix * (-1.0)).toDenseVector)))
+      (tmp.map(p => p * sinv))
   }
 
   private def computeP(yPoints: Seq[Point[_3D]], pars: similarityTransformationParameters): DenseMatrix[Double] = {
@@ -159,85 +169,86 @@ class BCPD(
     (0 until M).map { m =>
       val mvnd =
         MultivariateNormalDistribution(vectorizer.vectorize(yPoints(m)), DenseMatrix.eye[Double](dim) * pars.sigma2)
-      val e = math.exp(-pars.s / (2 * pars.sigma2) * trace(pars.sigma(m, m) * DenseMatrix.eye[Double](dim)))
+      val e = math.exp(-pars.s / (2 * pars.sigma2) * trace(DenseMatrix.eye[Double](dim).*(pars.sigma(m, m))))
       (0 until N).map { n =>
         Phi(m, n) = mvnd.pdf(vectorizer.vectorize(targetPoints(n))) * e * pars.alpha(m)
       }
     }
     val Pinit = Phi.copy * (1 - w)
     // TODO: Fix outlier distribution (section 4.3.4)
-    val c      = w * 1.0 / N.toDouble // * pout(x_n) see 4.3.4 (1/V)
+    val c = w * 1.0 / N.toDouble // * pout(x_n) see 4.3.4 (1/V)
     val denRow = DenseMatrix(sum(Pinit, Axis._0).t) * (1 - w) + c
-    val den    = tile(denRow, M, 1)
+    val den = tile(denRow, M, 1)
 
     Pinit /:/ den
   }
 
-  /** One iteration of BCPD
-    *
-    * @param YhatPoints
-    *   // template points (M)
-    * @param pars
-    *   // Similarity transformation parameters
-    * @return
-    */
+  /**
+   * One iteration of BCPD
+   *
+   * @param YhatPoints
+   *   // template points (M)
+   * @param pars
+   *   // Similarity transformation parameters
+   * @return
+   */
   def Iteration(
-      YhatPoints: Seq[Point[_3D]],
-      pars: similarityTransformationParameters
+    YhatPoints: Seq[Point[_3D]],
+    pars: similarityTransformationParameters
   ): (Seq[Point[_3D]], similarityTransformationParameters) = {
     // Update P and related terms
     val P = computeP(YhatPoints, pars)
 
-    val v    = sum(P, Axis._1)        // R^M Estimated number of target points matched with each source point
-    val v_   = sum(P, Axis._0).t.copy // R^N Posterior probability that x_n is a non-outlier
+    val v = sum(P, Axis._1) // R^M Estimated number of target points matched with each source point
+    val v_ = sum(P, Axis._0).t.copy // R^N Posterior probability that x_n is a non-outlier
     val Nhat = sum(v_)
 
-    val Pkron    = kron(P, Dmat)
-    val vkron    = kron(v.toDenseMatrix, D1Vec).toDenseVector
-    val v_kron   = kron(v_.toDenseMatrix, D1Vec).toDenseVector
-    val xhat     = pinv(diag(vkron)) * Pkron * X
+    val Pkron = kron(P, Dmat)
+    val vkron = kron(v.toDenseMatrix, D1Vec).toDenseVector
+    val v_kron = kron(v_.toDenseMatrix, D1Vec).toDenseVector
+    val xhat = pinv(diag(vkron)) * Pkron * X
     val xhatTinv = vectorInvTransform(xhat, pars)
 
     // Update Local deformations
     val s2divsigma2 = (math.pow(pars.s, 2) / pars.sigma2)
-    val SigmaInv    = GinvLambda + diag(v) * s2divsigma2
-    val Sigma       = pinv(SigmaInv)
-    val SigmaKron   = kron(Sigma, Dmat)
-    val vhat        = s2divsigma2 * SigmaKron * diag(vkron) * (xhatTinv - Y)
-    val uhat        = Y + vhat
+    val SigmaInv = GinvLambda + diag(v) * s2divsigma2
+    val Sigma = pinv(SigmaInv)
+    val SigmaKron = kron(Sigma, Dmat)
+    val vhat = SigmaKron.*(s2divsigma2) * diag(vkron) * (xhatTinv - Y)
+    val uhat = Y + vhat
 
     val alpha = v.map(f => math.exp(digamma(k + f) - digamma(k * M + Nhat)))
 
     // Update Similarity transform
-    val xMean     = sum((0 until M).map(m => v(m) * xhat(m * dim until m * dim + dim))) / Nhat
-    val uMean     = sum((0 until M).map(m => v(m) * uhat(m * dim until m * dim + dim))) / Nhat
+    val xMean = sum((0 until M).map(m => xhat(m * dim until m * dim + dim).*(v(m)))) / Nhat
+    val uMean = sum((0 until M).map(m => uhat(m * dim until m * dim + dim).*(v(m)))) / Nhat
     val sigma2bar = sum((0 until M).map(m => v(m) * Sigma(m, m))) / Nhat
 
     val Sxu = sum((0 until M).map { m =>
       val xm = xhat(m * dim until m * dim + dim)
       val um = uhat(m * dim until m * dim + dim)
-      v(m) * (xm - xMean) * (um - uMean).t
+      (xm - xMean).*(v(m)) * (um - uMean).t
     }) / Nhat
 
     val Suu = sum((0 until M).map { m =>
       val um = uhat(m * dim until m * dim + dim)
-      v(m) * (um - uMean) * (um - uMean).t + DenseMatrix.eye[Double](dim) * sigma2bar
+      (um - uMean).*(v(m)) * (um - uMean).t + DenseMatrix.eye[Double](dim) * sigma2bar
     }) / Nhat
 
     val svd.SVD(phi, _, psiT) = svd(Sxu)
-    val diagphipsi            = DenseVector.ones[Double](dim)
+    val diagphipsi = DenseVector.ones[Double](dim)
     diagphipsi(dim - 1) = det(phi * psiT)
 
     val R = phi * diag(diagphipsi) * psiT
     val s = trace(R * Sxu) / trace(Suu)
-    val t = xMean - s * R * uMean
+    val t = xMean - R.*(s) * uMean
 
     val newYhat = vectorTransform(Y + vhat, pars)
 
     val sXX = X.t * diag(v_kron) * X
     val sXY = X.t * Pkron.t * newYhat
     val sYY = newYhat.t * diag(vkron) * newYhat
-    val sC  = pars.sigma2 * sigma2bar
+    val sC = pars.sigma2 * sigma2bar
     val newSigma2 =
       (sXX - 2 * sXY + sYY + sC) / (Nhat * dim) // TODO: Should sC be added to all or included in the parenthesis as currently?
 
diff --git a/src/main/scala/gingr/other/algorithms/cpd/CPDFactory.scala b/src/main/scala/gingr/other/algorithms/cpd/CPDFactory.scala
index 201adb5..8fc5399 100644
--- a/src/main/scala/gingr/other/algorithms/cpd/CPDFactory.scala
+++ b/src/main/scala/gingr/other/algorithms/cpd/CPDFactory.scala
@@ -20,7 +20,7 @@ package gingr.other.algorithms.cpd
 import gingr.api.registration.utils.PointSequenceConverter
 import breeze.linalg.DenseMatrix
 import scalismo.common.Vectorizer
-import scalismo.geometry.{Point, _3D}
+import scalismo.geometry.{_3D, Point}
 
 /*
  Implementation of Point Set Registration: Coherent Point Drift
@@ -28,33 +28,34 @@ import scalismo.geometry.{Point, _3D}
  A python implementation already exists from where parts of the implementation is from: https://github.com/siavashk/pycpd
  */
 class CPDFactory(
-    val templatePoints: Seq[Point[_3D]],
-    val lambda: Double = 2,
-    val beta: Double = 2,
-    val w: Double = 0
+  val templatePoints: Seq[Point[_3D]],
+  val lambda: Double = 2,
+  val beta: Double = 2,
+  val w: Double = 0
 )(implicit
-    val vectorizer: Vectorizer[Point[_3D]]
+  val vectorizer: Vectorizer[Point[_3D]]
 ) {
-  val M: Int                        = templatePoints.length // num of reference points
-  val dim: Int                      = vectorizer.dim        // dimension
-  val G: DenseMatrix[Double]        = initializeKernelMatrixG(templatePoints, beta)
+  val M: Int = templatePoints.length // num of reference points
+  val dim: Int = vectorizer.dim // dimension
+  val G: DenseMatrix[Double] = initializeKernelMatrixG(templatePoints, beta)
   val template: DenseMatrix[Double] = PointSequenceConverter.toMatrix(templatePoints)
 
   require(0.0 <= w && w <= 1.0)
   require(beta > 0)
   require(lambda > 0)
 
-  /** Initialize G matrix - formula in paper fig. 4
-    *
-    * @param points
-    * @param beta
-    * @return
-    */
+  /**
+   * Initialize G matrix - formula in paper fig. 4
+   *
+   * @param points
+   * @param beta
+   * @return
+   */
   private def initializeKernelMatrixG(
-      points: Seq[Point[_3D]],
-      beta: Double
+    points: Seq[Point[_3D]],
+    beta: Double
   ): DenseMatrix[Double] = {
-    val M                      = points.length
+    val M = points.length
     val G: DenseMatrix[Double] = DenseMatrix.zeros[Double](M, M)
     (0 until M).map { i =>
       (0 until M).map { j =>
diff --git a/src/main/scala/gingr/other/algorithms/cpd/NonRigidCPD.scala b/src/main/scala/gingr/other/algorithms/cpd/NonRigidCPD.scala
index 59c0cb0..42c3020 100644
--- a/src/main/scala/gingr/other/algorithms/cpd/NonRigidCPD.scala
+++ b/src/main/scala/gingr/other/algorithms/cpd/NonRigidCPD.scala
@@ -17,19 +17,19 @@
 
 package gingr.other.algorithms.cpd
 
-import breeze.linalg.{Axis, DenseMatrix, DenseVector, diag, inv, sum}
+import breeze.linalg.{diag, inv, sum, Axis, DenseMatrix, DenseVector}
 import scalismo.common.Vectorizer
-import scalismo.geometry.{Point, _3D}
+import scalismo.geometry.{_3D, Point}
 
 /*
  Implementation of Point Set Registration: Coherent Point Drift (CPD) - Non-Rigid algorithm
  Paper: https://arxiv.org/pdf/0905.2635.pdf
  */
 private[cpd] class NonRigidCPD(
-    override val targetPoints: Seq[Point[_3D]],
-    override val cpd: CPDFactory
+  override val targetPoints: Seq[Point[_3D]],
+  override val cpd: CPDFactory
 )(implicit
-    vectorizer: Vectorizer[Point[_3D]]
+  vectorizer: Vectorizer[Point[_3D]]
 ) extends RigidCPD(targetPoints, cpd) {
   import cpd._
 
@@ -43,27 +43,26 @@ private[cpd] class NonRigidCPD(
   }
 
   override def Maximization(
-      X: DenseMatrix[Double],
-      Y: DenseMatrix[Double],
-      P: DenseMatrix[Double],
-      sigma2: Double
+    X: DenseMatrix[Double],
+    Y: DenseMatrix[Double],
+    P: DenseMatrix[Double],
+    sigma2: Double
   ): (DenseMatrix[Double], Double) = {
     // Update transform
     val P1: DenseVector[Double] = sum(P, Axis._1)
-    val Pt1                     = sum(P, Axis._0)
-    val Np                      = sum(P1)
+    val Pt1 = sum(P, Axis._0)
+    val Np = sum(P1)
 
-    val myG       = G
+    val myG = G
     val diagP1inv = inv(diag(P1))
-    val PX        = P * X
-
-    val A: DenseMatrix[Double] = G + lambda * sigma2 * diagP1inv
+    val PX = P * X
+    val A: DenseMatrix[Double] = G + diagP1inv.*(lambda * sigma2)
     val B: DenseMatrix[Double] = diagP1inv * PX - Y
 
     val W = A \ B
     // Update Point Cloud
     val deform = myG * W
-    val TY     = Y + deform
+    val TY = Y + deform
 
     // Update variance
     /*
diff --git a/src/main/scala/gingr/other/algorithms/cpd/RigidCPD.scala b/src/main/scala/gingr/other/algorithms/cpd/RigidCPD.scala
index 708260f..73889b4 100644
--- a/src/main/scala/gingr/other/algorithms/cpd/RigidCPD.scala
+++ b/src/main/scala/gingr/other/algorithms/cpd/RigidCPD.scala
@@ -18,33 +18,34 @@
 package gingr.other.algorithms.cpd
 
 import gingr.api.registration.utils.PointSequenceConverter
-import breeze.linalg.{Axis, DenseMatrix, DenseVector, det, diag, norm, svd, tile, trace, *}
+import breeze.linalg.{det, diag, norm, svd, tile, trace, Axis, DenseMatrix, DenseVector, *}
 import breeze.numerics.{abs, pow}
 import scalismo.common.Vectorizer
-import scalismo.geometry.{Point, _3D}
+import scalismo.geometry.{_3D, Point}
 
 /*
  Implementation of Point Set Registration: Coherent Point Drift (CPD) - Rigid algorithm
  Paper: https://arxiv.org/pdf/0905.2635.pdf
  */
 private[cpd] class RigidCPD(
-    val targetPoints: Seq[Point[_3D]],
-    val cpd: CPDFactory
+  val targetPoints: Seq[Point[_3D]],
+  val cpd: CPDFactory
 )(implicit
-    val vectorizer: Vectorizer[Point[_3D]]
+  val vectorizer: Vectorizer[Point[_3D]]
 ) {
   import cpd._
-  val N: Int                      = targetPoints.length
+  val N: Int = targetPoints.length
   val target: DenseMatrix[Double] = PointSequenceConverter.toMatrix(targetPoints)(vectorizer)
 
-  /** Initialize sigma2 - formula in paper fig. 4
-    * @param Y
-    * @param X
-    * @return
-    */
+  /**
+   * Initialize sigma2 - formula in paper fig. 4
+   * @param Y
+   * @param X
+   * @return
+   */
   private def initializeGaussianKernel(Y: Seq[Point[_3D]], X: Seq[Point[_3D]]): Double = {
-    val M   = Y.length
-    val N   = X.length
+    val M = Y.length
+    val N = X.length
     val dim = vectorizer.dim
     val s: Double = (0 until M).flatMap { m =>
       (0 until N).map { n =>
@@ -57,20 +58,27 @@ private[cpd] class RigidCPD(
   def Registration(max_iteration: Int, tolerance: Double = 0.001): Seq[Point[_3D]] = {
     val sigmaInit = initializeGaussianKernel(templatePoints, targetPoints)
 
-    val fit = (0 until max_iteration).foldLeft((cpd.template, sigmaInit)) { (it, i) =>
-      val currentSigma2 = it._2
+    var fit = (cpd.template, sigmaInit)
+    var i = 0
+    var converged = false
+
+    while (i < max_iteration && !converged) {
+      val currentSigma2 = fit._2
       println(s"CPD, iteration: ${i}, variance: ${currentSigma2}")
-      val iter      = Iteration(target, it._1, it._2)
-      val TY        = iter._1
+      val iter = Iteration(target, fit._1, fit._2)
+      val TY = iter._1
       val newSigma2 = iter._2
-      val diff      = abs(newSigma2 - currentSigma2)
+      val diff = abs(newSigma2 - currentSigma2)
       if (diff < tolerance) {
         println("Converged")
-        return PointSequenceConverter.toPointSequence(TY)(vectorizer)
+        fit = (TY, newSigma2)
+        converged = true
       } else {
-        iter
+        fit = (TY, newSigma2)
+        i += 1
       }
     }
+
     PointSequenceConverter.toPointSequence(fit._1)(vectorizer)
   }
 
@@ -89,23 +97,23 @@ private[cpd] class RigidCPD(
       }
       vec
     }
-    val c      = w / (1 - w) * math.pow((2.0 * math.Pi * sigma2), dim.toDouble / 2.0) * (M.toDouble / N.toDouble)
+    val c = w / (1 - w) * math.pow((2.0 * math.Pi * sigma2), dim.toDouble / 2.0) * (M.toDouble / N.toDouble)
     val denRow = DenseMatrix(breeze.linalg.sum(P, Axis._0).t)
-    val den    = tile(denRow, M, 1) + c
+    val den = tile(denRow, M, 1) + c
 
     P /:/ den
   }
 
   def Maximization(
-      X: DenseMatrix[Double],
-      Y: DenseMatrix[Double],
-      P: DenseMatrix[Double],
-      sigma2: Double
+    X: DenseMatrix[Double],
+    Y: DenseMatrix[Double],
+    P: DenseMatrix[Double],
+    sigma2: Double
   ): (DenseMatrix[Double], Double) = {
     // Update transform
     val P1: DenseVector[Double] = breeze.linalg.sum(P, Axis._1)
-    val Pt1                     = breeze.linalg.sum(P, Axis._0)
-    val Np                      = breeze.linalg.sum(P1)
+    val Pt1 = breeze.linalg.sum(P, Axis._0)
+    val Np = breeze.linalg.sum(P1)
 
     val muX = 1.0 / Np * X.t * P.t * DenseVector.ones[Double](P.rows)
     val muY = 1.0 / Np * Y.t * P1
@@ -113,18 +121,18 @@ private[cpd] class RigidCPD(
     val Xhat = X - DenseVector.ones[Double](N) * muX.t
     val Yhat = Y - DenseVector.ones[Double](M) * muY.t
 
-    val A                = Xhat.t * P.t * Yhat
+    val A = Xhat.t * P.t * Yhat
     val svd.SVD(u, _, v) = svd(A)
-    val C                = DenseVector.ones[Double](dim)
+    val C = DenseVector.ones[Double](dim)
     C(dim - 1) = det(u * v.t)
 
-    val R             = u * diag(C) * v
-    val s             = trace(A.t * R) / trace(Yhat.t * diag(P1) * Yhat)
-    val s1            = trace(Xhat.t * diag(Pt1) * Xhat)
-    val s2            = s * trace(A.t * R)
-    val t             = muX - s * R * muY
+    val R = u * diag(C) * v
+    val s = trace(A.t * R) / trace(Yhat.t * diag(P1) * Yhat)
+    val s1 = trace(Xhat.t * diag(Pt1) * Xhat)
+    val s2 = s * trace(A.t * R)
+    val t = muX - s * R * muY
     val updatedSigma2 = 1 / (Np * dim) * (s1 - s2)
-    val TY            = s * Y * R.t + DenseVector.ones[Double](M) * t.t
+    val TY = s * Y * R.t + DenseVector.ones[Double](M) * t.t
 
     (TY, updatedSigma2)
   }
diff --git a/src/main/scala/gingr/other/algorithms/icp/ICPFactory.scala b/src/main/scala/gingr/other/algorithms/icp/ICPFactory.scala
index 37f4eff..0e28b6c 100644
--- a/src/main/scala/gingr/other/algorithms/icp/ICPFactory.scala
+++ b/src/main/scala/gingr/other/algorithms/icp/ICPFactory.scala
@@ -19,7 +19,7 @@ package gingr.other.algorithms.icp
 
 import gingr.other.utils.Registrator
 import scalismo.common.{UnstructuredPoints, Vectorizer}
-import scalismo.geometry.{Point, _3D}
+import scalismo.geometry.{_3D, Point}
 
 /*
  Implementation of Point Set Registration: Iterative closest points todo: check and update rest of the comment
@@ -27,10 +27,10 @@ import scalismo.geometry.{Point, _3D}
  A python implementation already exists from where parts of the implementation is from: https://github.com/siavashk/pycpd
  */
 class ICPFactory(
-    val templatePoints: UnstructuredPoints[_3D]
+  val templatePoints: UnstructuredPoints[_3D]
 )(implicit
-    val vectorizer: Vectorizer[Point[_3D]],
-    registrator: Registrator
+  val vectorizer: Vectorizer[Point[_3D]],
+  registrator: Registrator
 ) {
   def registerRigidly(targetPoints: UnstructuredPoints[_3D]): RigidICP = {
     new RigidICP(targetPoints, this)
diff --git a/src/main/scala/gingr/other/algorithms/icp/NonRigidOptimalStepICP.scala b/src/main/scala/gingr/other/algorithms/icp/NonRigidOptimalStepICP.scala
index 330f717..a5974b4 100644
--- a/src/main/scala/gingr/other/algorithms/icp/NonRigidOptimalStepICP.scala
+++ b/src/main/scala/gingr/other/algorithms/icp/NonRigidOptimalStepICP.scala
@@ -18,10 +18,10 @@
 package gingr.other.algorithms.icp
 
 import gingr.api.registration.utils.{NonRigidClosestPointRegistrator, PointSequenceConverter}
-import breeze.linalg.{CSCMatrix, SparseVector, diag}
+import breeze.linalg.{diag, CSCMatrix, SparseVector}
 import gingr.other.utils.CSCHelper
 import scalismo.common.{DomainWarp, PointId, UnstructuredPoints, Vectorizer}
-import scalismo.geometry.{Landmark, Point, _3D}
+import scalismo.geometry.{_3D, Landmark, Point}
 import scalismo.mesh.{TriangleCell, TriangleMesh}
 
 /*
@@ -29,17 +29,17 @@ import scalismo.mesh.{TriangleCell, TriangleMesh}
  */
 
 abstract class NonRigidOptimalStepICP(
-    templateMesh: TriangleMesh[_3D],
-    targetMesh: TriangleMesh[_3D],
-    templateLandmarks: Seq[Landmark[_3D]],
-    targetLandmarks: Seq[Landmark[_3D]],
-    gamma: Double = 1.0
+  templateMesh: TriangleMesh[_3D],
+  targetMesh: TriangleMesh[_3D],
+  templateLandmarks: Seq[Landmark[_3D]],
+  targetLandmarks: Seq[Landmark[_3D]],
+  gamma: Double = 1.0
 )(implicit
-    vectorizer: Vectorizer[Point[_3D]],
-    warper: DomainWarp[_3D, TriangleMesh]
+  vectorizer: Vectorizer[Point[_3D]],
+  warper: DomainWarp[_3D, TriangleMesh]
 ) {
   require(gamma >= 0)
-  val n: Int   = templateMesh.pointSet.numberOfPoints // Number of nodes
+  val n: Int = templateMesh.pointSet.numberOfPoints // Number of nodes
   val dim: Int = vectorizer.dim
 
   private val commonLmNames = templateLandmarks.map(_.id) intersect targetLandmarks.map(_.id)
@@ -53,7 +53,7 @@ abstract class NonRigidOptimalStepICP(
       .map(lm => targetMesh.pointSet.findClosestPoint(lm.point).point)
   val UL: CSCMatrix[Double] = CSCHelper.DenseMatrix2CSCMatrix(PointSequenceConverter.toMatrix(lmPointsOnTarget))
 
-  private val edges   = trianglesToEdges(templateMesh.triangulation.triangles)
+  private val edges = trianglesToEdges(templateMesh.triangulation.triangles)
   val numOfEdges: Int = edges.length
 
   val M: CSCMatrix[Double] = InitializeMatrixM(edges)
@@ -84,30 +84,32 @@ abstract class NonRigidOptimalStepICP(
   }
 
   def Registration(
-      max_iteration: Int,
-      tolerance: Double = 0.001,
-      alpha: Seq[Double] = defaultAlpha,
-      beta: Seq[Double] = defaultBeta
+    max_iteration: Int,
+    tolerance: Double = 0.001,
+    alpha: Seq[Double] = defaultAlpha,
+    beta: Seq[Double] = defaultBeta
   ): TriangleMesh[_3D] = {
     require(alpha.length == beta.length)
 
     val fit = (alpha.zip(beta)).zipWithIndex.foldLeft(templateMesh) { (temp, config) =>
       val (a, b) = config._1
-      val j      = config._2
+      val j = config._2
 
-      val innerFit = (0 until max_iteration).foldLeft((temp, Double.PositiveInfinity)) { (it, i) =>
-        val iter    = Iteration(it._1, targetMesh, a, b)
-        val TY      = iter._1
+      var innerFit = (temp, Double.PositiveInfinity)
+      var i = 0
+      while (i < max_iteration && innerFit._2 >= tolerance) {
+        val iter = Iteration(innerFit._1, targetMesh, a, b)
+        val TY = iter._1
         val newDist = iter._2
         println(
           s"ICP, iteration: ${j * max_iteration + i}/${max_iteration * alpha.length}, alpha: ${a}, beta: ${b}, average distance to target: ${newDist}"
         )
+        innerFit = (TY, newDist)
         if (newDist < tolerance) {
           println("Converged")
-          return TY
-        } else {
-          (iter._1, iter._2)
+          innerFit = (TY, newDist)
         }
+        i += 1
       }
       innerFit._1
     }
@@ -115,21 +117,21 @@ abstract class NonRigidOptimalStepICP(
   }
 
   def getClosestPoints(
-      template: TriangleMesh[_3D],
-      target: TriangleMesh[_3D]
+    template: TriangleMesh[_3D],
+    target: TriangleMesh[_3D]
   ): (Seq[Point[_3D]], Seq[Double], Double) = {
     val (corr, dist) =
       NonRigidClosestPointRegistrator.ClosestPointTriangleMesh3D.closestPointCorrespondence(template, target)
     val cp = corr.map(_._2)
-    val w  = corr.map(_._3)
+    val w = corr.map(_._3)
     (cp, w, dist)
   }
 
   def Iteration(
-      template: TriangleMesh[_3D],
-      target: TriangleMesh[_3D],
-      alpha: Double,
-      beta: Double
+    template: TriangleMesh[_3D],
+    target: TriangleMesh[_3D],
+    alpha: Double,
+    beta: Double
   ): (TriangleMesh[_3D], Double, IndexedSeq[Point[_3D]])
 }
 
@@ -137,22 +139,22 @@ abstract class NonRigidOptimalStepICP(
  Implementation of "N-ICP-T" from the paper: "Optimal Step Nonrigid ICP Algorithms for Surface Registration"
  */
 class NonRigidOptimalStepICP_T(
-    templateMesh: TriangleMesh[_3D],
-    targetMesh: TriangleMesh[_3D],
-    templateLandmarks: Seq[Landmark[_3D]],
-    targetLandmarks: Seq[Landmark[_3D]],
-    gamma: Double = 1.0
+  templateMesh: TriangleMesh[_3D],
+  targetMesh: TriangleMesh[_3D],
+  templateLandmarks: Seq[Landmark[_3D]],
+  targetLandmarks: Seq[Landmark[_3D]],
+  gamma: Double = 1.0
 )(implicit
-    vectorizer: Vectorizer[Point[_3D]],
-    warper: DomainWarp[_3D, TriangleMesh]
+  vectorizer: Vectorizer[Point[_3D]],
+  warper: DomainWarp[_3D, TriangleMesh]
 ) extends NonRigidOptimalStepICP(templateMesh, targetMesh, templateLandmarks, targetLandmarks, gamma) {
   val B1: CSCMatrix[Double] = CSCMatrix.zeros[Double](numOfEdges, dim)
 
   override def Iteration(
-      template: TriangleMesh[_3D],
-      target: TriangleMesh[_3D],
-      alpha: Double,
-      beta: Double
+    template: TriangleMesh[_3D],
+    target: TriangleMesh[_3D],
+    alpha: Double,
+    beta: Double
   ): (TriangleMesh[_3D], Double, IndexedSeq[Point[_3D]]) = {
     require(alpha >= 0.0)
     require(beta >= 0.0)
@@ -162,11 +164,11 @@ class NonRigidOptimalStepICP_T(
     val W = diag(SparseVector(w: _*))
 
     val lmPointsOnTemplate = lmIdsOnTemplate.map(id => template.pointSet.point(id))
-    val VL                 = CSCHelper.DenseMatrix2CSCMatrix(PointSequenceConverter.toMatrix(lmPointsOnTemplate))
+    val VL = CSCHelper.DenseMatrix2CSCMatrix(PointSequenceConverter.toMatrix(lmPointsOnTemplate))
 
     val U = CSCHelper.DenseMatrix2CSCMatrix(PointSequenceConverter.toMatrix(cp))
     // N-ICP-T
-    val V                     = CSCHelper.DenseMatrix2CSCMatrix(PointSequenceConverter.toMatrix(template.pointSet.points.toSeq))
+    val V = CSCHelper.DenseMatrix2CSCMatrix(PointSequenceConverter.toMatrix(template.pointSet.points.toSeq))
     val A1: CSCMatrix[Double] = M * alpha
     val A2: CSCMatrix[Double] = W * CSCHelper.eye(template.pointSet.numberOfPoints)
     val A3: CSCMatrix[Double] = CSCMatrix.zeros[Double](lmPointsOnTemplate.length, M.cols)
@@ -193,17 +195,17 @@ class NonRigidOptimalStepICP_T(
  */
 
 class NonRigidOptimalStepICP_A(
-    templateMesh: TriangleMesh[_3D],
-    targetMesh: TriangleMesh[_3D],
-    templateLandmarks: Seq[Landmark[_3D]],
-    targetLandmarks: Seq[Landmark[_3D]],
-    gamma: Double = 1.0
+  templateMesh: TriangleMesh[_3D],
+  targetMesh: TriangleMesh[_3D],
+  templateLandmarks: Seq[Landmark[_3D]],
+  targetLandmarks: Seq[Landmark[_3D]],
+  gamma: Double = 1.0
 )(implicit
-    vectorizer: Vectorizer[Point[_3D]],
-    warper: DomainWarp[_3D, TriangleMesh]
+  vectorizer: Vectorizer[Point[_3D]],
+  warper: DomainWarp[_3D, TriangleMesh]
 ) extends NonRigidOptimalStepICP(templateMesh, targetMesh, templateLandmarks, targetLandmarks, gamma) {
 
-  private val G                         = diag(SparseVector(1, 1, 1, gamma))
+  private val G = diag(SparseVector(1, 1, 1, gamma))
   private val kronMG: CSCMatrix[Double] = CSCHelper.kroneckerProduct(M, G)
 
   val B1: CSCMatrix[Double] = CSCMatrix.zeros[Double](4 * numOfEdges, dim)
@@ -226,7 +228,7 @@ class NonRigidOptimalStepICP_A(
   }
 
   private def ComputeMatrixDL(template: TriangleMesh[_3D]): CSCMatrix[Double] = {
-    val nPoints              = lmIdsOnTemplate.length
+    val nPoints = lmIdsOnTemplate.length
     val D: CSCMatrix[Double] = CSCMatrix.zeros[Double](nPoints, 4 * n)
     lmIdsOnTemplate.zipWithIndex.foreach { case (id, i) =>
       val p = template.pointSet.point(id).toArray :+ 1.0
@@ -239,10 +241,10 @@ class NonRigidOptimalStepICP_A(
   }
 
   override def Iteration(
-      template: TriangleMesh[_3D],
-      target: TriangleMesh[_3D],
-      alpha: Double,
-      beta: Double
+    template: TriangleMesh[_3D],
+    target: TriangleMesh[_3D],
+    alpha: Double,
+    beta: Double
   ): (TriangleMesh[_3D], Double, IndexedSeq[Point[_3D]]) = {
     require(alpha >= 0.0)
     require(beta >= 0.0)
diff --git a/src/main/scala/gingr/other/algorithms/icp/RigidICP.scala b/src/main/scala/gingr/other/algorithms/icp/RigidICP.scala
index 255160f..5543eed 100644
--- a/src/main/scala/gingr/other/algorithms/icp/RigidICP.scala
+++ b/src/main/scala/gingr/other/algorithms/icp/RigidICP.scala
@@ -20,44 +20,51 @@ package gingr.other.algorithms.icp
 import breeze.numerics.abs
 import gingr.other.utils.Registrator
 import scalismo.common.{PointId, UnstructuredPoints, Vectorizer}
-import scalismo.geometry.{Point, _3D}
+import scalismo.geometry.{_3D, Point}
 
 private[icp] class RigidICP(
-    val targetPoints: UnstructuredPoints[_3D],
-    val icp: ICPFactory
+  val targetPoints: UnstructuredPoints[_3D],
+  val icp: ICPFactory
 )(implicit
-    val vectorizer: Vectorizer[Point[_3D]],
-    registrator: Registrator
+  val vectorizer: Vectorizer[Point[_3D]],
+  registrator: Registrator
 ) {
   require(vectorizer.dim == 2 || vectorizer.dim == 3)
 
   def Registration(max_iteration: Int, tolerance: Double = 0.001): UnstructuredPoints[_3D] = {
     val sigmaInit = 0.0
 
-    val fit = (0 until max_iteration).foldLeft((icp.templatePoints, sigmaInit)) { (it, i) =>
-      val iter     = Iteration(it._1, targetPoints)
+    var fit = (icp.templatePoints, sigmaInit)
+    var i = 0
+    var converged = false
+
+    while (i < max_iteration && !converged) {
+      val iter = Iteration(fit._1, targetPoints)
       val distance = iter._2
-      println(s"ICP, iteration: ${i}, distance: ${distance}")
-      val TY   = iter._1
-      val diff = abs(distance - it._2)
+      println(s"ICP, iteration: $i, distance: $distance")
+      val TY = iter._1
+      val diff = abs(distance - fit._2)
       if (diff < tolerance) {
         println("Converged")
-        return TY
+        fit = (TY, distance)
+        converged = true
       } else {
-        iter
+        fit = iter
       }
+      i += 1
     }
+
     fit._1
   }
 
   private def attributeCorrespondences(
-      template: UnstructuredPoints[_3D],
-      target: UnstructuredPoints[_3D]
+    template: UnstructuredPoints[_3D],
+    target: UnstructuredPoints[_3D]
   ): (Seq[(Point[_3D], Point[_3D])], Double) = {
-    val ptIds    = template.pointIds.toIndexedSeq
+    val ptIds = template.pointIds.toIndexedSeq
     var distance = 0.0
     val corr = ptIds.map { (id: PointId) =>
-      val pt                  = template.point(id)
+      val pt = template.point(id)
       val closestPointOnMesh2 = target.findClosestPoint(pt).point
       distance += (pt - closestPointOnMesh2).norm
       (pt, closestPointOnMesh2)
@@ -66,11 +73,11 @@ private[icp] class RigidICP(
   }
 
   def Iteration(
-      template: UnstructuredPoints[_3D],
-      target: UnstructuredPoints[_3D]
+    template: UnstructuredPoints[_3D],
+    target: UnstructuredPoints[_3D]
   ): (UnstructuredPoints[_3D], Double) = {
     val (correspondences, distance) = attributeCorrespondences(template, target)
-    val registeredPoints            = registrator.register(correspondences)
+    val registeredPoints = registrator.register(correspondences)
     (registeredPoints, distance)
   }
 
diff --git a/src/main/scala/gingr/other/utils/CSCHelper.scala b/src/main/scala/gingr/other/utils/CSCHelper.scala
index b6d732f..0faecd0 100644
--- a/src/main/scala/gingr/other/utils/CSCHelper.scala
+++ b/src/main/scala/gingr/other/utils/CSCHelper.scala
@@ -38,8 +38,8 @@ object CSCHelper {
     require(matrices.forall(m => m.cols == matrices(0).cols), "Not all matrices have the same number of columns")
     val numRows = matrices.foldLeft(0)(_ + _.rows)
     val numCols = matrices(0).cols
-    val res     = CSCMatrix.zeros[Double](numRows, numCols)
-    var offset  = 0
+    val res = CSCMatrix.zeros[Double](numRows, numCols)
+    var offset = 0
     for (m <- matrices) {
       var i = 0
       while (i < m.cols) {
diff --git a/src/main/scala/gingr/other/utils/PoseRegistrator.scala b/src/main/scala/gingr/other/utils/PoseRegistrator.scala
index 2e0dae1..938377c 100644
--- a/src/main/scala/gingr/other/utils/PoseRegistrator.scala
+++ b/src/main/scala/gingr/other/utils/PoseRegistrator.scala
@@ -18,7 +18,7 @@
 package gingr.other.utils
 
 import scalismo.common.UnstructuredPoints
-import scalismo.geometry.{Point, _3D}
+import scalismo.geometry.{_3D, Point}
 import scalismo.registration.LandmarkRegistration
 
 trait Registrator {
diff --git a/src/main/scala/gingr/simple/GingrInterface.scala b/src/main/scala/gingr/simple/GingrInterface.scala
index d064131..d7e6ff4 100644
--- a/src/main/scala/gingr/simple/GingrInterface.scala
+++ b/src/main/scala/gingr/simple/GingrInterface.scala
@@ -10,7 +10,7 @@ import gingr.api.registration.config.{
   IcpRegistrationState
 }
 import gingr.api.sampling.evaluators.{EvaluationMode, ModelToTargetEvaluation}
-import scalismo.geometry.{Landmark, _3D}
+import scalismo.geometry.{_3D, Landmark}
 import scalismo.mesh.TriangleMesh
 import scalismo.statisticalmodel.PointDistributionModel
 import scalismo.transformations.TranslationAfterRotation
@@ -19,15 +19,15 @@ import scalismo.utils.Random
 import java.io.File
 
 case class GingrInterface(
-    model: PointDistributionModel[_3D, TriangleMesh],
-    target: TriangleMesh[_3D],
-    initialModelParameterTransform: Option[TranslationAfterRotation[_3D]] = None,
-    modelLandmarks: Option[Seq[Landmark[_3D]]] = None,
-    targetLandmarks: Option[Seq[Landmark[_3D]]] = None,
-    evaluatorUncertainty: Double = 1.0,
-    evaluatedPoints: Option[Int] = None,
-    evaluationMode: EvaluationMode = ModelToTargetEvaluation,
-    logFileFittingParameters: Option[File] = None
+  model: PointDistributionModel[_3D, TriangleMesh],
+  target: TriangleMesh[_3D],
+  initialModelParameterTransform: Option[TranslationAfterRotation[_3D]] = None,
+  modelLandmarks: Option[Seq[Landmark[_3D]]] = None,
+  targetLandmarks: Option[Seq[Landmark[_3D]]] = None,
+  evaluatorUncertainty: Double = 1.0,
+  evaluatedPoints: Option[Int] = None,
+  evaluationMode: EvaluationMode = ModelToTargetEvaluation,
+  logFileFittingParameters: Option[File] = None
 )(implicit rnd: Random) {
   def CPD(config: CpdConfiguration): SimpleRegistrator[CpdRegistrationState, CpdConfiguration, CpdRegistration] = {
     val algorithm = new CpdRegistration()
diff --git a/src/main/scala/gingr/simple/SimpleModels.scala b/src/main/scala/gingr/simple/SimpleModels.scala
index 3ded7b7..0a605d5 100644
--- a/src/main/scala/gingr/simple/SimpleModels.scala
+++ b/src/main/scala/gingr/simple/SimpleModels.scala
@@ -27,36 +27,36 @@ sealed trait WhichKernel {
   val printpars: String
 }
 case class InvLapKernel(scaling: Double) extends WhichKernel {
-  override val name: String      = "InvLap"
+  override val name: String = "InvLap"
   override val printpars: String = scaling.toString
 }
 case class InvLapDotKernel(scaling: Double, gamma: Double) extends WhichKernel {
-  override val name: String      = "InvLapDot"
+  override val name: String = "InvLapDot"
   override val printpars: String = scaling.toString + "_" + gamma.toString
 }
 case class GaussKernel(scaling: Double, sigma: Double) extends WhichKernel {
-  override val name: String      = "Gauss"
+  override val name: String = "Gauss"
   override val printpars: String = scaling.toString + "_" + sigma.toString
 }
 case class GaussMixKernel() extends WhichKernel {
-  override val name: String      = "GaussMix"
+  override val name: String = "GaussMix"
   override val printpars: String = ""
 }
 case class GaussDotKernel(scaling: Double, sigma: Double) extends WhichKernel {
-  override val name: String      = "GaussDot"
+  override val name: String = "GaussDot"
   override val printpars: String = scaling.toString + "_" + sigma.toString
 }
 case class GaussMirrorKernel(scaling: Double, sigma: Double) extends WhichKernel {
-  override val name: String      = "GaussMirror"
+  override val name: String = "GaussMirror"
   override val printpars: String = scaling.toString + "_" + sigma.toString
 }
 
 object SimpleTriangleModels3D {
 
   def create(
-      reference: TriangleMesh[_3D],
-      kernelSelect: WhichKernel,
-      relativeTolerance: Double = 0.01
+    reference: TriangleMesh[_3D],
+    kernelSelect: WhichKernel,
+    relativeTolerance: Double = 0.01
   ): PointDistributionModel[_3D, TriangleMesh] = {
     val model = kernelSelect match {
       case InvLapKernel(s) =>
diff --git a/src/test/scala/DummyTest.scala.scala b/src/test/scala/DummyTest.scala.scala
new file mode 100644
index 0000000..3c56e1c
--- /dev/null
+++ b/src/test/scala/DummyTest.scala.scala
@@ -0,0 +1,3 @@
+package gingr
+
+@main def Test: Unit = assert(1 > 0)