Page MenuHomePhabricator

Doubts about global fiber tracking
Closed, ResolvedPublic

Description

Thanks for providing your implementation of global fiber tracking. When comparing your implementation with the paper by Reisert et al., MICCAI 2009, we have the following doubts:

  1. mitkGibbEnergyComputer.cpp, line 203ff: It appears to us that the computation of the internal energy has its sign flipped w.r.t,

to the definition by Reisert.

  1. mitkGibbEnergyComputer.cpp, line 146: What is the reasoning of adding m_ParticleChemicalPotential here?
  2. mitkGibbEnergyComputer.cpp, line 152: In our understanding of the derivation, the energy term should be computed as:

float energy = odfVal/m_ParticleWeight-2*modelVal+D2, where D2 is some constant.

  1. mitkEnergyComputer.cp, lines 150ff: This function does not approximate the modified Bessel function I0, although its use in

mitkGibbEnergyComputer.cpp, line 143 supposes to do so.
Thanks for your consideration. Best wishes, Frithjof

Event Timeline

I Frithjof. Thank you for looking into our implementation! A bit of double checking is always appreciated. It's been a while since I looked at the code but i will try my best :)

  1. mitkGibbEnergyComputer.cpp, line 203ff: It appears to us that the computation of the internal energy has its sign flipped w.r.t, to the definition by Reisert.

You are right that "ComputeInternalEnergyConnection" returns -E(M) instead of E(M), but in successive calculations Reisert uses exp(-E(M)) in his paper. So we directly apply the minus sign inside of "ComputeInternalEnergyConnection".

  1. mitkGibbEnergyComputer.cpp, line 146: What is the reasoning of adding m_ParticleChemicalPotential here?

Good question. If I recall correctly it influences the probability of accepting particles that do not fit well with their neighborhood. So if there are many particles that are not in line with the new particle and m_ParticleChemicalPotential is really small, this increases the external energy because bw is negative. This is realte to your question 4. We got the implementation mbesseli0 directly from Marco Reisert, since we included his code in MITK. I'm not sure how these values came to be, but mbesseli0 should correspond to I_0 in the paper. But I'm also not sure how close the implementation is to it's theoretical foundation :)
In general it also works with the raw dot product, since it boils down to the idea that parallel particles should contribute much and anti-parallel particles should contribute less signal.

To get rid of the dependency on the FRT-ODF stuff, I wanted to implement a version that directly works on the raw data instead of the Q-ball ODFs for years now, but I never got around to it.

As for point 3, I think you are right there. I will fix this. It won't actually have much of an effect since this is/can be compensated by the particle weight.

Hello Peter,

thanks for your quick and extensive answer. I do not want to bug you too much with my requests... we have our own
implementation of global tracking, and I just wanted to share some of my concerns about this implementation.

Re 1: mitkGibbEnergyComputer.cpp: in successive calculations Reisert uses exp(-E(M)) in his paper.
Yes, I know. But both internal and external energy are added with the same sign (e.g., mitkMetropolisHastingsSampler.cpp:127),
so I think that there is some snag here.

Re: 4: Thanks for your explanation.

So if there are many particles that are not in line with the new particle and m_ParticleChemicalPotential is really small,
this increases the external energy because bw is negative.

Given the definition in the paper, I do not think that I0 can ever be negative. Thinking over the integral again, I do not think
that the Bessel function is involved here. Instead of deriving a closed form, I have written a small function to integrate this
expression numerically. I find that I0 does not depend much on the angle ni, nj - and might even be approximated by a constant.

Let me know if you really want to pick this up again. Thanks again for your kind help. Frithjof

Don't worry about bugging me, we really vlaue any feedback! If you find more, don't hesitate to let me know. I'm also happy to integrate changes/new stuff. If you already implemented changes you can also attach a patch here or contribute it in one of the following ways: http://mitk.org/wiki/How_to_contribute

  1. I will look into the sign of the energy calculation again.
  1. I think you are right that it does not really reflect the intgral. I found it really difficlut in parts to relate the original implementation with the paper. Happy to test you implementation of I0!
This comment was removed by neher.

I checked the formulas again and

float energy = 2*odfVal/m_ParticleWeight-modelVal - (mbesseli0(1.0)+m_ParticleChemicalPotential);

looks correct to me. This is basically -( |Fm|² - <Fm,D>) so the signe is ok in the implementation.

I agree with your assessment that -( |Fm|² - <Fm,D>) is implemented here, but has a sign opposite to Reisert's paper.
The internal energy is computed with the same sign as in the paper. In mitkMetropolisHastingsSampler.cpp, I find that
both energies are summed up (e.g., in birth proposal). As stated before, my comments are based on a code review,
not a test, so it is likely that I am missing a point here.

The internal energy is also calculated with the opposite sign, as stated in my first comment. Here is the complete formula and how it is implemented:

Paper:

P(M|D) = exp( -E_int -E_ext ) = exp( -1/l²(  ||x_1 + a_1*l*n_1 - C ||² + ||x_2 + a_2*l*n_2 - C||²  - L) -( |Fm|² - 2<Fm,D>) ) 
= exp( -(norm1 + norm2 - m_ConnectionPotential)*m_IntStrength - ( modelVal + (mbesseli0(1.0)+m_ParticleChemicalPotential) - 2*odfVal/m_ParticleWeight)*m_ExtStrength  ) 
= exp( (m_ConnectionPotential - norm1 - norm2)*m_IntStrength + ( 2*odfVal/m_ParticleWeight - modelVal - (mbesseli0(1.0)+m_ParticleChemicalPotential) )*m_ExtStrength  )

So in the third line we get rid of the minus signs and change all the signs inside the brackets accordingly. Therefore, we add the outputs of ComputeInternalEnergy and ComputeExternalEnergy instead of subtracting them.

If I got everyting right here, this shows that the implementation fits the formula in the paper.

Yes, I agree. Both potentials have their signs flipped w.r.t. the paper. Thanks for your assessment.

kislinsk triaged this task as Normal priority.May 30 2017, 8:17 AM
kislinsk added a subscriber: kislinsk.

Is this task resolved? 😄

neher claimed this task.