unibas-gravis / scalismo Goto Github PK

Given the recurrent problems with stl, we started relying more on other formats such as vtk and ply (which is supported in scalismo-faces).

It would be great if we could also support PLY directly in Scalismo (without the texture features of course).

Currently, tests are broken when one or more spaces are contained in the absolute path to the project folder. The returned path of the resource then contains the escaped sequence "%20" instead of spaces.

A fix could be to use a URLDecoder to transform the paths such the "%20" are decoded with spaces again.

This is due to the RandomSVD approximation behaving differently at each call. A temporary solution on the course branch was to seed the random generator.

The findClosestPoint function on an image domain does not work correctly on images, whose direction matrix is negative. The reason for this is that the internal function pointToContinuousIndex does not take the direction matrix into account when computing the continuous indices.

This function would be very useful to use when computing a distance map from a mesh as it would give a more precise value than the current implementation relying on the closest vertex.

However, distance maps are usually iterated over in parallel when performing registration, hence the issue.

The function

PivotedCholesky.computeApproximateEigGeneric

takes a parameter D. It is from the outside not clear what this parameter does and why
it is needed. Either it has to be documented, or (if possible) hidden from the interface.

A very common use case is to transform landmarks with a rigid transformation (e.g. when doing rigid alignment). Due to the uncertainty of the landmark, which has to be transformed along with it, the transform operation is non-trivial. It would be nice if the landmark class would have a transform method to do that.

Long running jobs that do a lot of MeshIO run out of memory at some point.
The problem is most likely that memory is leaked when using VTK.

I have downloaded basel face database. Please guide me how to download "01_MorphableModel.mat" file. it is not available on basel website.
i am waiting for your kind reply.
thanks

The following code

 case VTK_CHAR | VTK_SIGNED_CHAR =>
         val in = arrayVTK.asInstanceOf[vtkCharArray].GetJavaArray()

results in a classCastException as avtkSignedCharArray cannot be cast into
a vtkCharArray. A solution would be to treat the cases VTK_CHAR and VTK_SIGNED_CHAR separately.

This is due to the fact that a Rotation could have been created from a Rotation space around any given center, it does not make sense to return parameters as one would not know on which Rotation Space to use them.

The reference mesh of a DataCollection is set to the Procrustes mean when performing a GPA.

Hence the following code returns false

val dc :  DataCollection = ???
val dcAfterGPA = DataCollection.gpa(dc)
dc.reference == dcAfterGPA.reference // false

I would argue that this is unintuitive and that only the transformations should change.

One way of creating a MultivariateNormalDistribution is by specifying the mean and three orthonormal directions (with associated variance) (see https://github.com/unibas-gravis/scalismo/blob/master/src/main/scala/scalismo/statisticalmodel/MultivariateNormalDistribution.scala#L159)

The computation does not compute the correct covariance matrix.

The default Pivoted Cholesky stopping criterion (PivotedCholesky.RelativeTolerance(0)) can experience a NotConvergedException when using a very small dataset.

The following code will in most of the runs throw an error when using the default stopping criteria. If the stopping criteria instead is set to be numberOfEigenfunctions = data-1, then it consistently works.

val model = StatismoIO.readStatismoMeshModel(new File("model2017-1_bfm_nomouth-2500.h5")).get

val ref = model.referenceMesh

val data = (1 to 3).map(f => model.sample())

println(s"data length: ${data.length}")
val dc = DataCollection.fromMeshSequence(referenceMesh = ref, registeredMeshes = data)._1.get
val dcGpa = DataCollection.gpa(dc)
StatisticalMeshModel.createUsingPCA(dcGpa, NumberOfEigenfunctions.apply(data.length-1)).get
println("pca 1 computed")
StatisticalMeshModel.createUsingPCA(dcGpa).get // NotConvergedException
println("pca 2 computed")

Full error message (might need to run multiple times to get the error message - depends on the samples that are drawn from the BFM)

Exception in thread "main" breeze.linalg.NotConvergedException: 
	at breeze.linalg.svd$.breeze$linalg$svd$$doSVD_Double(svd.scala:110)
	at breeze.linalg.svd$Svd_DM_Impl$.apply(svd.scala:40)
	at breeze.linalg.svd$Svd_DM_Impl$.apply(svd.scala:39)
	at breeze.generic.UFunc.apply(UFunc.scala:48)
	at breeze.generic.UFunc.apply$(UFunc.scala:46)
	at breeze.linalg.svd$.apply(svd.scala:23)
	at scalismo.numerics.PivotedCholesky$.computeApproximateEigGeneric(PivotedCholesky.scala:193)
	at scalismo.numerics.PivotedCholesky$.computeApproximateEig(PivotedCholesky.scala:216)
	at scalismo.statisticalmodel.DiscreteLowRankGaussianProcess$.createUsingPCA(DiscreteLowRankGaussianProcess.scala:457)
	at scalismo.statisticalmodel.StatisticalMeshModel$.createUsingPCA(StatisticalMeshModel.scala:315)
	at scalismo.statisticalmodel.StatisticalMeshModel$.createUsingPCA(StatisticalMeshModel.scala:300)
	at apps.bfm.PCAtest$.main(PCAtest.scala:36)
	at apps.bfm.PCAtest.main(PCAtest.scala)

I have a set of meshes with correspondence, but as an artefact of the registration process, a number of the vertices are overlapping in each mesh.
When I import the meshes with DataCollection.fromMeshDirectory the overlapping vertices are removed, and I end up with a different number of points in each mesh, and consequently there is no longer any correspondence?
What is the correct way to solve it?
It would be nice to have a option not to eliminate vertices in the mesh import functions.

As an intermediate result in the Nyström method, an Eigendecomposition of the kernel matrix, defined on a subset of the domain points has to be computed. For efficiency reasons, this is currently done using the Pivoted Cholesky Decomposition. In the current implementation, the accuracy is not specified as a relative tolerance, but only a fixed number of basis vectors are computed. This can lead to bad approximations of the true eigenvectors.

A fix is to use a relative tolerance in the computation.

diceCoefficient from MeshMetrics.scala should be modified, at least, with a grid sampler that does not rely on a random number of samples.

On a Win7 x64 system, running JRE 1.8.0, I get very frequent crashes when trying out the first tutorials.

Steps:
1.) Double-click launch_windows.bat to start scalismoLab
2.) Open Documents > Fast Track > Rigid Alignment
3.) Select the first piece of code in the tutorial window (val mesh : ... ) and hit Shift-Enter to execute
4.) Wait for a few seconds, moving the mouse back and forth between the three windows --> crash

Attached are two log files
hs_err_pid8796.log
hs_err_pid11008.log

In ASM fitting findBestMatchingPointAtPoint can return 0 points if the filter values maxIntensityStddev and maxShapeStddev are too low.

if (minIntensityDist < config.maxIntensityStddev && shapeDistForPt < config.maxShapeStddev) Some(minPt) else None

This leads to an error when 0 training data is used to get the coefficients of the shape model.

val coeffs = asm.shapeModel.coefficients(refPtIdsWithTargetPt, sigma2 = 1e-6)

Add the robust Huber-loss function to the library. The loss function is less sensitive to outlier data and is often used with noisy data.

As of now SSM only allows for rigid transformations. The suggestion is to extend to also allow for similarity transform. At first, perform a rigid transform and then a scaling transform for the basis vectors according to the computed scaling factor between original and target landmarks (as done in similarity3DLandmarkRegistration).
Reference mesh and mean deformation field should also be updated.

The method MappedSurfaceProperty.triangulation produces an AbstractMethodError at runtime

current implementation: override val triangulation: TriangleList = values.triangulation
probable fix: override def triangulation: TriangleList = values.triangulation

(replace val by def)

to reproduce:

def produceAME(pointProperty: SurfacePointProperty): Unit = {
    pointProperty.map{x => x}.triangulation
}

The interpolation of a DiscreteGaussianProcess throws an exception, on pretty much any input.
A test should be added to ensure that it functions correctly.

The class DiscreteScalarField takes as a type parameter a subclass of DiscreteDomain. Clearly, a DiscreteScalarField[_, DD, _], which is defined on a DiscreteDomain DD which is a subclass of a DiscreteDomain D, can be replaced anywhere where a DiscreteScalarField[_, D, _].
Hence it should be covariant. (This is already correctly specified for DiscreteField).

currently Landmark IO methods are generic and rely on type parameters : readLandmarksJson[_3D]
On the other hand, image reading methods are named differently for different dimensionality: read3DScalarImage

Ideally, one of the approaches should be used for both Landmarks and Images; preferably, generic methods for both cases.

The next test because when the objective point is already in the tree.If I put a point not in the tree works ok

it should "find the points in a ball around excluding himself" in {

    Given("some points")
    val points = Seq(
      RNDensePoint(0, 0),
      RNDensePoint(0, 1),
      RNDensePoint(0, 2)
    )

    And("the tree with the points")
    val tree = VantagePointTree(points, RNDensePoint.euclidean)

    And("an objetive point")
    val objetive = RNDensePoint(0, 1)

    And("the ball radio")
    val radio = 4.0

    And("the expected results")
    val expected = Set(
      RNDensePoint(0, 0),
      RNDensePoint(0, 1),
      RNDensePoint(0, 2)
    )

    When("the tree find points in the ball")
    val closestPoints = tree.findEpsilonNeighbours(objetive, radio)

    Then("the nearest must be the expected")
    closestPoints.toSet shouldBe expected
  }

Method 'instance' from StatisticalMeshModel.java is a little bit slow for my needs -> it takes a couple of milliseconds to execute. Is it possible to accelerate it somehow, using parallelization and/or using hardware acceleration (GPU)?

The Average distance and Hausdorff distance functions to measure the distance between two meshes, both use the closest point in pointset function. The distance should instead be taken to the closest point on the surface.
It could happen that there is 0 distance to the surface on the other mesh, whereas the vertices are not perfectly aligned, which would show as a distance > 0.

mesh.triangulation.adjacentTrianglesForPoint() currently results in a IndexOutOfBoundsException when accessed for a pointId not belonging to any cell (same for adjacentPointsForPoint). The reason is that the pointIds used for the computation of this array are obtained by iterating over the cells only.

This comes from the assumption that the meshes read into Scalismo must be "correct", i.e. without any floating vertices. This condition is however not enforced at read time.

Here we should either enforce that meshes are correct at read time (e.g. throw an exception and request to activate a repair flag) or fix the arrays above to return an empty Seq[TriangleId] for a point not belonging to any cell.

The domain of an image is too small when the LinearImageInterpolator is used to interpolate an image. More precisely, it cuts the size of the domain by one voxel, due to the incorrect use of boundingbox for images.

This results in runtime errors when an image is interpolated and then resampled with the same image grid.

MeshMetrics.procrustesDistance does not compute the Procrustes distance correctly. It correctly computes the rigid alignment that minimizes the squared error between corresponding points. However, it then uses the method avgDistance to compute the distances, which computes the distance to the closest point and not the corresponding point.

This class is defined here:
https://github.com/unibas-gravis/scalismo-faces/blob/master/src/main/scala/scalismo/faces/mesh/Mesh.scala#L43

The move would require to also move the RGBA class. The advantage would be that we could release this change to the official scalismo-ui to visualize vertex colored meshes, without the dependency on scalismo-faces.

When using the function TriangleMesh3D.operations.closestPointOnSurface on a triangle mesh with a degenerated triangle, i.e. a triangle with three colinear points, the returned closest point is not always correct. This comes from the fact that calculations used to construct the bounding sphere based index structure and others during a query result in NaN values.

The current version of the code assumes that no degenerated triangles are present in a mesh. Especially the case of a triangle with only colinear points leads to the above-described behaviour. Such triangles are easily produced when using VTK on a volume image with the steps thresholding, contour-filtering, smoothing, decimating.

Here we should introduce additional safety even that we will trade off a bit of performance.

Integrator needs a Random object although it never actually generates the samples. There are even deterministic integration methods - I think randomness should be handled by the Sampler not the integrator. It might even be removed from the Sampler interface as it is only required by some samplers and I currently do not see the value of having it everywhere.

Affects: scalismo.numerics.Sampler and scalismo.numerics.Integrator

case class Integrator[D <: Dim: NDSpace](sampler: Sampler[D])(implicit rng: Random)

trait Sampler[D <: Dim] {
  def sample()(implicit rand: Random): IndexedSeq[(Point[D], Double)]
}

Currently a result of a closestPoint operation can be either of type ClosestPoint, which indicates that only the closest point and the distance was computed, or one of the subclasses ClosestPointIsVertex , ClosestPointOnLine or ClosestPointInTriangle, which also provide access to meta data.

However, these types should all be distinct and ClosestPoint should not be the base class, but rather a separate type, in order to make it possible to understand from the type signature, which information is available.

DiscreteDomain calculates its hashCode based on the points iterator:
DiscreteDomain, L73

override def hashCode() = points.hashCode()

The hash value then depends on the concrete iterator instance rather than the actual point sequence.

Consequence: Meshes (which use UnstructuredPointsDomain) do not hash properly

• Equality uses points.toSeq to compare - It is thus not affected by this bug

Our util.Random object creates a seeded breeze RandBasis internally. The current call creates identical random number generators for each thread. In a concurrent application with multiple threads you easily end up with threads which all generate the same random numbers (same seed for generators).

Manual seeding is now deprecated. It is possible to use RandBasis.withSeed(seed) which ensures initialization of each thread with a different seed. This at least ensures proper behavior in the unseeded case and only leads to non-reproducible results if you seed it. It will not constantly regenerate the same "random" numbers.

General remark: Seeded random generators do not go at all with concurrent programs with an ill-defined execution order! In these situations make sure you explicitly create random generators per job/batch. Single, synchronized random sources will not work (ill-defined order of execution) and also the typical per-thread-abstraction is almost always not what you expect, e.g. in a ThreadPool threads are allocated randomly (an reused) to execute your jobs.

I often use samples of models, and from time to time I use them in a par-environment
However this leads to ugly behavior and repeated values for randomly generated numbers

the issue arises from the both different ways to draw random values:
scalismo/src/main/scala/scalismo/utils/Random.scala

The following code exploits the Issue:

object TestRandom extends App{
scalismo.initialize()
implicit val rand = Random(1986)
val dim = 5
(0 until 10).par.foreach(f => {
val standardNormal = Gaussian(0, 1)(rand.breezeRandBasis)
println("breeze: " + standardNormal.sample(dim))
println("scalaRandom:" + IndexedSeq.fill(dim)(rand.scalaRandom.nextGaussian()))
}
)
}

The result on a CPU with 8 threads is the following

breeze: Vector(0.46599886070617236, -0.12885325103464149, 1.058985097854918, -0.42287179689623533, 0.7199570539297094)
breeze: Vector(0.46599886070617236, -0.12885325103464149, 1.058985097854918, -0.42287179689623533, 0.7199570539297094)
breeze: Vector(0.46599886070617236, -0.12885325103464149, 1.058985097854918, -0.42287179689623533, 0.7199570539297094)
breeze: Vector(0.46599886070617236, -0.12885325103464149, 1.058985097854918, -0.42287179689623533, 0.7199570539297094)
breeze: Vector(0.46599886070617236, -0.12885325103464149, 1.058985097854918, -0.42287179689623533, 0.7199570539297094)
breeze: Vector(0.46599886070617236, -0.12885325103464149, 1.058985097854918, -0.42287179689623533, 0.7199570539297094)
breeze: Vector(0.46599886070617236, -0.12885325103464149, 1.058985097854918, -0.42287179689623533, 0.7199570539297094)
breeze: Vector(0.46599886070617236, -0.12885325103464149, 1.058985097854918, -0.42287179689623533, 0.7199570539297094)
scalaRandom:Vector(2.6161687929921436, 0.2360950134905143, 0.8244545674853196, 0.06492437537385155, 0.20930382403878728)
scalaRandom:Vector(0.177437296826372, 0.48657608758097104, -0.8857920497973689, 0.473885303306742, 1.5259638060371277)
scalaRandom:Vector(1.5126100549074337, -0.14147546744354175, -1.9873098137265353, -0.28301772753378446, 1.2807364396457321)
scalaRandom:Vector(-0.8651463600754313, 1.0063388340603403, 2.1126032645499624, 1.229623373736994, 0.4149855473947657)
scalaRandom:Vector(0.394497947362624, -0.7393869056632871, -0.17787934329181332, 1.6869857048905423, 1.0229810189041688)
scalaRandom:Vector(-0.6280144016292583, 0.3647319891335585, -2.2806044448504994, 1.402314541890555, -0.6547768767613916)
scalaRandom:Vector(0.13926731439585482, 0.6065878188194251, 1.6987171740770581, 0.07523230373380176, -1.6504160243355648)
breeze: Vector(1.2244046440044116, -0.24603691367390496, 0.03245754960911975, 1.0163216358427813, -0.9333993201074842)
breeze: Vector(1.2244046440044116, -0.24603691367390496, 0.03245754960911975, 1.0163216358427813, -0.9333993201074842)
scalaRandom:Vector(-1.7852408613435788, -1.1608923517856675, 0.5384125740849208, 0.13609386079673452, -1.5371141011242655)
scalaRandom:Vector(0.5405371986279452, -1.2386742393392802, -0.26588409904826177, -0.4826293117366927, 0.1533976096574605)
scalaRandom:Vector(-0.6517619731602057, 0.22941031281788035, -0.23238720213806432, 2.1542631616243217, -0.7935287454787473)

breezeRandBasis therefore does not procude what you would like to have (in most applications). It seems that every thread starts from the same seed.

behind breezerandbasis there is a scalismo function:
implicit val breezeRandBasis: RandBasis = new RandBasis(new ThreadLocalRandomGenerator(new MersenneTwister(seed)))

perhaps the ThreadLocalRandomGenerator is the problem?

The current implementation allows to get reproducible results but on the cost of all parallel stuff you want to do?

scalismo uses this at many spots in the code

Perhaps the best solution would be to define if you want to have a global random generator or a threadlocal-one when building the random object from a seed?
implicit val rand = Random(1986, 'global')

While the README mentions the online course about shape modelling, it does not directly link to the PMM website under https://gravis.dmi.unibas.ch/PMM/.
Adding such a link would be a good way to attract more visitors (especially those coming from search engines) to the tutorials on the PMM website.

The objects and methods to read/write statistical mesh models in Scalismo are named after Statismo. This can be confusing for the following reasons:

1. For people unfamiliar with Statismo will at best be seen as an weird abbreviation of statistical model
2. Once we add io methods for other types of statistical models apart from mesh models (such as e.g. deformation models), the name statismo does not distinguish these conepts.

I propose instead to name the IO object StatisticalModelIO and name the methods
readStatisticalMeshModel, readStatisticalDeformationModel, etc.

For the transition period, we can keep the existing IO methods and mark them as deprecated.

Why is the lossfunction value in the Mean Huber Loss normalized by (1 + v * v)?

This differs from the classical definition of the Huber Loss https://en.wikipedia.org/wiki/Huber_loss

and creates a discontinuity since (v * v / 2) / (1 + v * v) and (delta * (Math.abs(v) - delta / 2)) have different values for v = delta.

(Also, in practice, the current implementation does not seem to work for me, while changing it to the classical definition does).

Currently there is no way of going from Point[_3D] to Index[_3D] on a discrete image domain.
The best way would be that DiescreteImageDomain implements the findClosestPoint method of SpatiallyIndexedDiscreteDomain.

The interface of DiscreteGaussianProcess is incomplete. We should have an interpolate method, as we have it for the DiscreteLowRankGaussianProcess. Also marginal and other functions are not yet fully supported.

Update the VTK library to the newest version. This will allow running headless mode - thereby avoiding X forwarding.

It would be useful if one could combine some additional attributes, such as the subject's age, weight, height, gender, or even some properties of the shape itself, e.g. length or width of the bone, together with the shape model.

An example use case of such a combined model would be to complete a fractured bone by conditioning the model, not only on the observed healthy part, but additionally on the patient's attributes: age, weight, length of the contralateral bone , etc., which should improve the completion results.

The method presented in this paper https://link.springer.com/chapter/10.1007/978-3-319-64870-5_5 could for example be used in the implementation (generalized to all shapes, not just faces).

DiscreteScalarImage features a resampling method. Currently this method simply interpolates the image and resamples it on a new grid. However, it does not do a smoothing in case the image is subsampled.

The DiscreteGaussianProcess is a finite view obtained by sampling a finite number of points from a GaussianProcess. This is only sound if the sampling is done uniformly, and with sufficiently many points such that the characteristics of the GP can be represented. However, the current implementation does not make sure that this is the case.

Problems that appear include:

• The basis functions are not necessarily orthogonal
• The probabilities (e.g. logpdf, pdf) of the distribution is dependent on the sampling