[Neuroimaging] [DIPY] propagator anisotropy estimation using MAP(L)MRI
Ping-Hong Yeh
pinghongyeh at gmail.com
Thu Mar 29 12:02:23 EDT 2018
Hi Rutger,
We have some bad PA maps created using default settings, and I would like
to hear your opinions on improving the fitting.
Attached are the screenshots of PA_GCV, norm_laplacian, L_opt and
PA_laplacian_weighted0.2 maps.
I am currently running the fitting using 0.05 for the minimum bound of the
GCV, but I am not sure if that would help.
In order to do comparisons between controls and disease population, we need
to make sure that the same fitting parameters are applied for the MAPMRI
fitting for avoiding any biases. Do you have suggestions regarding this
matter?
Thank you.
Ping
On Tue, Jan 23, 2018 at 7:42 AM, Rutger Fick <fick.rutger at gmail.com> wrote:
> Hi Ping,
>
> Salt and pepper noise is not a good sign (I just didn't see it so much on
> the second set of slices you sent). To spot badly estimated voxels is
> typically pretty easy - RTOP and many others can have negative or huge
> values, which typically come from oscillations in the signs extrapolation.
> You can often see these as bright spots in the laplacian norm.
>
> If you go through the data and see that salt and pepper noise corresponds
> to unusually high norms, Increasing the laplacian minimum weight in the
> code as i told you wil usually resolve it (or fixing it to a value like
> 0.05 or 0.1 or something, see what works without overdoing it).
>
> Best,
> Rutger
>
>
>
>
> On 23 Jan 2018 03:06, "Ping-Hong Yeh" <pinghongyeh at gmail.com> wrote:
>
> Hi Rutger,
>
> Thank you very much for the detailed reply.
>
> I guess i do not need to worry about those salt-pepper dots?
>
> Would you recommend output laplacian norm and laplacian_weighted maps and
> go through the images for each data set? Any tips for realizing something
> really goes wrong when looking at the propagator anisotropy map?
>
> Best,
>
> Ping
>
>
> On Jan 22, 2018 6:55 PM, "Rutger Fick" <fick.rutger at gmail.com> wrote:
>
>> Hi Ping,
>>
>> In my experience, badly estimated voxels typically have super high
>> laplacian norm and very low estimated laplacian weight (lopt).
>> Looking at these results I would say things actually look pretty good!
>>
>> Getting the best results is always tricky finding a balance of optimally
>> regularizing: not fitting the noise but also not over-regularizing, which
>> is why the GCV option is nice.
>> But, in rare cases it does mess up. So, if you want to give the GCV a bit
>> less freedom to go low (to be on the safe side) you can increase the
>> minimum bound of the GCV optimization in line 2272 of the code.
>>
>> There's many ways to speed up the code I gave you if you want to put in
>> the effort ;-) Using parallel processing is not standardly implemented in
>> dipy, but maybe you can hack it somehow.
>> You can also set the laplacian_weight = 0.1 or something to avoid GCV,
>> but it won't make a huge difference. I only ever used this code to do
>> research - so speed was not really a concern.
>>
>> Anyway, hope this all helped! Let me know if everything works out,
>>
>> Best,
>> Rutger
>>
>> On 19 January 2018 at 22:03, Ping-Hong Yeh <pinghongyeh at gmail.com> wrote:
>>
>>> Hi Rutger,
>>>
>>> Attached please find the snapshot of norm_of_laplacian_signal, lopt,
>>> and pa maps of the same data set i used earlier.
>>>
>>> BTW, is there a way to speed up the mapmri_pa processing? Will the
>>> OpenMP help?
>>>
>>> Thank you,
>>>
>>> ping
>>>
>>> On Fri, Jan 19, 2018 at 1:25 PM, Rutger Fick <fick.rutger at gmail.com>
>>> wrote:
>>>
>>>> Hi Ping,
>>>>
>>>> So far, so good.
>>>> In my opinion the TORTOISE PA reconstruction looks a bit
>>>> flat/overregularized - but then again I don't know what kind of
>>>> regularization they implemented for themselves.
>>>> The PA of the implementation I gave you seems to give more consistent
>>>> contrast for different tissue configurations - which is a good - but looks
>>>> like it under-regularizes in some individual voxels (the salt-pepper noise
>>>> in the PA/RTOP).
>>>>
>>>> To check if this is the case, can you show me the
>>>> mapfit_L.norm_of_laplacian_signal() and mapfit_L.lopt maps?
>>>>
>>>> Rutger
>>>>
>>>>
>>>>
>>>>
>>>> On 19 January 2018 at 17:43, Ping-Hong Yeh <pinghongyeh at gmail.com>
>>>> wrote:
>>>>
>>>>> Hi Rutger,
>>>>>
>>>>> Just give you an update of the results (see the attached snapshots)
>>>>> using GCV weighted and Laplacian regularization for MAPMRI
>>>>> estimation.
>>>>>
>>>>> The other PA mapping was calculated using TORTOISE. I have also
>>>>> attached RTOP mapping calculated using DIPY with and without GCV
>>>>> weighted and Laplacian regularization.
>>>>>
>>>>> Comparing to the TORTOISE, PA values in the one using GCV weighted
>>>>> and Laplacian regularization method are relatively smaller,
>>>>> particularly over the regions with the less dense white matter.
>>>>>
>>>>> For RTOP images, I am not sure whether GCV weighted and Laplacian
>>>>> regularization method outperforms the one without using GCV weighted
>>>>> and Laplacian regularization.
>>>>>
>>>>> Any comments?
>>>>> Thank you,
>>>>>
>>>>> ping
>>>>>
>>>>> On Wed, Jan 17, 2018 at 7:48 PM, Rutger Fick <fick.rutger at gmail.com>
>>>>> wrote:
>>>>>
>>>>>> Hi Ping,
>>>>>>
>>>>>> If it's still running and gave only that error then probably it was
>>>>>> just a single voxel that failed and the rest is working. However, I
>>>>>> recommend you first try to fit a smaller dataset (just a few voxels or a
>>>>>> single slice) just to check the results make sense.
>>>>>>
>>>>>> I should mention that the code I gave you is slower than Dipy's
>>>>>> public version for reasons I won't get into, so don't worry if you have to
>>>>>> wait longer than before.
>>>>>>
>>>>>> Best,
>>>>>> Rutger
>>>>>>
>>>>>> On 18 Jan 2018 00:58, "Ping-Hong Yeh" <pinghongyeh at gmail.com> wrote:
>>>>>>
>>>>>>> Hi Rutger,
>>>>>>>
>>>>>>> Thanks again for the prompt reply.
>>>>>>>
>>>>>>> Adding "mask" to mapmri have fixed the error; however, another error
>>>>>>> shows up,
>>>>>>>
>>>>>>> mapfit_L = map_model_L.fit(data,mask=data[..., 0]>0)
>>>>>>> dipy/core/geometry.py:129: RuntimeWarning: invalid value encountered
>>>>>>> in true_divide
>>>>>>> theta = np.arccos(z / r)
>>>>>>> dipy/reconst/mapmri_pa.py:364: UserWarning: The MAPMRI positivity
>>>>>>> constraint depends on CVXOPT (http://cvxopt.org/). CVXOPT is
>>>>>>> licensed under the GPL (see: http://cvxopt.org/copyright.html) and
>>>>>>> you may be subject to this license when using the positivity constraint.
>>>>>>> warn(w_s)
>>>>>>> dipy/reconst/mapmri_pa.py:413: UserWarning: Optimization did not
>>>>>>> find a solution
>>>>>>> warn('Optimization did not find a solution')
>>>>>>> Error: Couldn't find per display information
>>>>>>>
>>>>>>>
>>>>>>> It is still running though. Should i stop the running?
>>>>>>>
>>>>>>> Thank you.
>>>>>>> ping
>>>>>>>
>>>>>>> On Tue, Jan 16, 2018 at 7:18 PM, Rutger Fick <fick.rutger at gmail.com>
>>>>>>> wrote:
>>>>>>>
>>>>>>>> Hi Ping,
>>>>>>>>
>>>>>>>> Reading the error messages, it looks like you're fitting a masked
>>>>>>>> voxel. The following error:
>>>>>>>>
>>>>>>>> /Library/Python/2.7/site-packages/dipy/reconst/mapmri_pa.py:389:
>>>>>>>> RuntimeWarning: invalid value encountered in divide
>>>>>>>> data = np.asarray(data / data[self.gtab.b0s_mask].mean())
>>>>>>>>
>>>>>>>> says you're dividing by either zero or NaN, which means your b0
>>>>>>>> value of that voxel was zero (or you had no b0 values possibly). Note that
>>>>>>>> mapmri needs at least one b0 measurement.
>>>>>>>> I recommend you check if it works when you fit a voxel that you
>>>>>>>> know for sure is in white matter. If it works, you can do something like
>>>>>>>> map_model_L.fit(data, mask=data[..., 0]>0) to use a mask that only
>>>>>>>> fits if the first measured DWI is positive (assuming your first measurement
>>>>>>>> is a b0).
>>>>>>>>
>>>>>>>> Best,
>>>>>>>> Rutger
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>> On 16 January 2018 at 23:46, Ping-Hong Yeh <pinghongyeh at gmail.com>
>>>>>>>> wrote:
>>>>>>>>
>>>>>>>>> Hi Rutger,
>>>>>>>>>
>>>>>>>>> I got an error running the map_model.fit using mapmri_pa. Here is
>>>>>>>>> the scripts i used,
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> map_model_L = mapmri_pa.MapmriModel(gtab,
>>>>>>>>> radial_order=radial_order,
>>>>>>>>> laplacian_regularization=True, #
>>>>>>>>> this regularization enhances reproducibility of estimated q-space indices
>>>>>>>>> by imposing smoothness
>>>>>>>>> laplacian_weighting="GCV", # this
>>>>>>>>> makes it use generalized cross-validation to find the best regularization
>>>>>>>>> weight
>>>>>>>>> positivity_constraint=True) # this
>>>>>>>>> ensures the estimated PDF is positive
>>>>>>>>>
>>>>>>>>> mapfit_L = map_model_L.fit(data)
>>>>>>>>>
>>>>>>>>> , and the error message,
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> /Library/Python/2.7/site-packages/dipy/core/geometry.py:129:
>>>>>>>>> RuntimeWarning: invalid value encountered in true_divide
>>>>>>>>> theta = np.arccos(z / r)
>>>>>>>>> /Library/Python/2.7/site-packages/dipy/reconst/mapmri_pa.py:364:
>>>>>>>>> UserWarning: The MAPMRI positivity constraint depends on CVXOPT (http:
>>>>>>>>> xopt.org/). CVXOPT is licensed under the GPL (see:
>>>>>>>>> http://cvxopt.org/copyright.html) and you may be subject to this
>>>>>>>>> license when using positivity constraint.
>>>>>>>>> warn(w_s)
>>>>>>>>> /Library/Python/2.7/site-packages/dipy/reconst/mapmri_pa.py:389:
>>>>>>>>> RuntimeWarning: invalid value encountered in divide
>>>>>>>>> data = np.asarray(data / data[self.gtab.b0s_mask].mean())
>>>>>>>>> /Library/Python/2.7/site-packages/dipy/reconst/mapmri_pa.py:413:
>>>>>>>>> UserWarning: Optimization did not find a solution
>>>>>>>>> warn('Optimization did not find a solution')
>>>>>>>>> /Library/Python/2.7/site-packages/dipy/reconst/mapmri_pa.py:444:
>>>>>>>>> UserWarning: Optimization did not find a solution
>>>>>>>>> warn('Optimization did not find a solution')
>>>>>>>>> Traceback (most recent call last):
>>>>>>>>> File "<stdin>", line 1, in <module>
>>>>>>>>> File "/Library/Python/2.7/site-pack
>>>>>>>>> ages/dipy/reconst/multi_voxel.py", line 33, in new_fit
>>>>>>>>> fit_array[ijk] = single_voxel_fit(self, data[ijk])
>>>>>>>>> File "/Library/Python/2.7/site-pack
>>>>>>>>> ages/dipy/reconst/mapmri_pa.py", line 465, in fit
>>>>>>>>> coef_iso = coef_iso / sum(coef_iso * self.Bm_iso)
>>>>>>>>> UnboundLocalError: local variable 'coef_iso' referenced before
>>>>>>>>> assignment
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> Any suggestions?
>>>>>>>>>
>>>>>>>>> Thank you.
>>>>>>>>>
>>>>>>>>> ping
>>>>>>>>>
>>>>>>>>> On Fri, Jan 12, 2018 at 6:24 PM, Rutger Fick <
>>>>>>>>> fick.rutger at gmail.com> wrote:
>>>>>>>>>
>>>>>>>>>> Hi Ping,
>>>>>>>>>>
>>>>>>>>>> Attached is the mapmri code that also has PA, just put it in the
>>>>>>>>>> dipy/reconst/ folder (where also the current mapmri.py file is) and run
>>>>>>>>>> "python setup.py install" from dipy's main folder. That should make it
>>>>>>>>>> usable in the same way as the current mapmri module.
>>>>>>>>>> Note that its based on an old implementation that still works
>>>>>>>>>> with the "cvxopt" optimizer package, so you'll have to install cvxopt to
>>>>>>>>>> make it run.
>>>>>>>>>>
>>>>>>>>>> I recommend you use the model with both laplacian regularization
>>>>>>>>>> and positivity constraint, this give the best results in general.
>>>>>>>>>>
>>>>>>>>>> from dipy.reconst import mapmri_pa
>>>>>>>>>> mapmod = mapmri_pa.MapmriModel(gtab,
>>>>>>>>>> laplacian_regularization=True, #
>>>>>>>>>> this regularization enhances reproducibility of estimated q-space indices
>>>>>>>>>> by imposing smoothness
>>>>>>>>>> laplacian_weighting="GCV", # this
>>>>>>>>>> makes it use generalized cross-validation to find the best regularization
>>>>>>>>>> weight
>>>>>>>>>> positivity_constraint=True) #
>>>>>>>>>> this ensures the estimated PDF is positive
>>>>>>>>>> mapfit = mapmod.fit(data)
>>>>>>>>>> pa = mapfit.pa()
>>>>>>>>>>
>>>>>>>>>> Aside from the original MAPMRI citation for Ozarslan et al.
>>>>>>>>>> (2013), note that the relevant citation for dipy's laplacian-regularized
>>>>>>>>>> MAP-MRI implementation is [1].
>>>>>>>>>> [1] Fick, Rutger HJ, et al. "MAPL: Tissue microstructure
>>>>>>>>>> estimation using Laplacian-regularized MAP-MRI and its application to HCP
>>>>>>>>>> data." *NeuroImage* 134 (2016): 365-385.
>>>>>>>>>>
>>>>>>>>>> Hope it helps and let me know if you need anything else,
>>>>>>>>>> Rutger
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> On 12 January 2018 at 21:48, Ping-Hong Yeh <pinghongyeh at gmail.com
>>>>>>>>>> > wrote:
>>>>>>>>>>
>>>>>>>>>>> Hi Roger,
>>>>>>>>>>>
>>>>>>>>>>> Thanks for the prompt reply.
>>>>>>>>>>> May I have the code for estimating PA?
>>>>>>>>>>>
>>>>>>>>>>> Ping
>>>>>>>>>>>
>>>>>>>>>>> On Jan 12, 2018 3:21 PM, "Rutger Fick" <fick.rutger at gmail.com>
>>>>>>>>>>> wrote:
>>>>>>>>>>>
>>>>>>>>>>>> Hi Ping,
>>>>>>>>>>>>
>>>>>>>>>>>> MAPL is just a name for using laplacian-regularized MAP-MRI. If
>>>>>>>>>>>> you're using the dipy mapmri implementation then you're using MAPL by
>>>>>>>>>>>> default.
>>>>>>>>>>>> From a fitted mapmri model you can estimate overall
>>>>>>>>>>>> non-gaussianity using fitted_model.ng(), and parallel and perpendicular
>>>>>>>>>>>> non-Gaussianity using ng_parallel() and ng_perpendic
>>>>>>>>>>>> perpendicularular().
>>>>>>>>>>>> Propagator Anisotropic is not included in the current dipy
>>>>>>>>>>>> implementation. However, I do have a personal version of dipy's mapmri
>>>>>>>>>>>> implementation that includes it, if you're interested.
>>>>>>>>>>>>
>>>>>>>>>>>> Best,
>>>>>>>>>>>> Rutger
>>>>>>>>>>>>
>>>>>>>>>>>> On 12 January 2018 at 16:49, Ping-Hong Yeh <
>>>>>>>>>>>> pinghongyeh at gmail.com> wrote:
>>>>>>>>>>>>
>>>>>>>>>>>>> Hi DIPY users,
>>>>>>>>>>>>>
>>>>>>>>>>>>> I would like to know the way of estimating non-Gaussian and
>>>>>>>>>>>>> PA, mentioned in the Avram et al. “Clinical feasibility of
>>>>>>>>>>>>> using mean apparent propagator (MAP) MRI to characterize brain tissue
>>>>>>>>>>>>> microstructure” paper, using MAPMRI or MAPL model.
>>>>>>>>>>>>>
>>>>>>>>>>>>> Thank you.
>>>>>>>>>>>>>
>>>>>>>>>>>>> Ping
>>>>>>>>>>>>>
>>>>>>>>>>>>> _______________________________________________
>>>>>>>>>>>>> Neuroimaging mailing list
>>>>>>>>>>>>> Neuroimaging at python.org
>>>>>>>>>>>>> https://mail.python.org/mailman/listinfo/neuroimaging
>>>>>>>>>>>>>
>>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>> _______________________________________________
>>>>>>>>>>>> Neuroimaging mailing list
>>>>>>>>>>>> Neuroimaging at python.org
>>>>>>>>>>>> https://mail.python.org/mailman/listinfo/neuroimaging
>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>> _______________________________________________
>>>>>>>>>>> Neuroimaging mailing list
>>>>>>>>>>> Neuroimaging at python.org
>>>>>>>>>>> https://mail.python.org/mailman/listinfo/neuroimaging
>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> _______________________________________________
>>>>>>>>>> Neuroimaging mailing list
>>>>>>>>>> Neuroimaging at python.org
>>>>>>>>>> https://mail.python.org/mailman/listinfo/neuroimaging
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>
>>>>>>>>> _______________________________________________
>>>>>>>>> Neuroimaging mailing list
>>>>>>>>> Neuroimaging at python.org
>>>>>>>>> https://mail.python.org/mailman/listinfo/neuroimaging
>>>>>>>>>
>>>>>>>>>
>>>>>>>>
>>>>>>>> _______________________________________________
>>>>>>>> Neuroimaging mailing list
>>>>>>>> Neuroimaging at python.org
>>>>>>>> https://mail.python.org/mailman/listinfo/neuroimaging
>>>>>>>>
>>>>>>>>
>>>>>>>
>>>>>>> _______________________________________________
>>>>>>> Neuroimaging mailing list
>>>>>>> Neuroimaging at python.org
>>>>>>> https://mail.python.org/mailman/listinfo/neuroimaging
>>>>>>>
>>>>>>>
>>>>>> _______________________________________________
>>>>>> Neuroimaging mailing list
>>>>>> Neuroimaging at python.org
>>>>>> https://mail.python.org/mailman/listinfo/neuroimaging
>>>>>>
>>>>>>
>>>>>
>>>>> _______________________________________________
>>>>> Neuroimaging mailing list
>>>>> Neuroimaging at python.org
>>>>> https://mail.python.org/mailman/listinfo/neuroimaging
>>>>>
>>>>>
>>>>
>>>> _______________________________________________
>>>> Neuroimaging mailing list
>>>> Neuroimaging at python.org
>>>> https://mail.python.org/mailman/listinfo/neuroimaging
>>>>
>>>>
>>>
>>> _______________________________________________
>>> Neuroimaging mailing list
>>> Neuroimaging at python.org
>>> https://mail.python.org/mailman/listinfo/neuroimaging
>>>
>>>
>>
>> _______________________________________________
>> Neuroimaging mailing list
>> Neuroimaging at python.org
>> https://mail.python.org/mailman/listinfo/neuroimaging
>>
>>
> _______________________________________________
> Neuroimaging mailing list
> Neuroimaging at python.org
> https://mail.python.org/mailman/listinfo/neuroimaging
>
>
>
> _______________________________________________
> Neuroimaging mailing list
> Neuroimaging at python.org
> https://mail.python.org/mailman/listinfo/neuroimaging
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/neuroimaging/attachments/20180329/b1481044/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: NCNC3308_pa_lapW0.2.png
Type: image/png
Size: 139310 bytes
Desc: not available
URL: <http://mail.python.org/pipermail/neuroimaging/attachments/20180329/b1481044/attachment-0004.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: NCNC3308_lopt.png
Type: image/png
Size: 142237 bytes
Desc: not available
URL: <http://mail.python.org/pipermail/neuroimaging/attachments/20180329/b1481044/attachment-0005.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: NCNC3308_norm_laplacian.png
Type: image/png
Size: 138960 bytes
Desc: not available
URL: <http://mail.python.org/pipermail/neuroimaging/attachments/20180329/b1481044/attachment-0006.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: NCNC3308_pa.png
Type: image/png
Size: 154492 bytes
Desc: not available
URL: <http://mail.python.org/pipermail/neuroimaging/attachments/20180329/b1481044/attachment-0007.png>
More information about the Neuroimaging
mailing list