Re: skfmm

[Jason Furtney]

John,

I had a closer look at the paper your linked me. I was not aware of
the fast sweeping method so thanks for that. It looks like, at least
for 2D distance functions, that this method should be faster than the
fast marching method.

As a test, I implemented the first-order version of this algorithm
which is described in this paper:
http://www.math.uci.edu/~zhao/homepage/research_files/FSM.pdf

The code is here:
https://gist.github.com/jkfurtney/890a1db0f68a36b5720d19ca4352afac

For 2D, this gives the same results as skfmm.distance(phi,order=1) but
possibly faster?

For distance functions only 2^n (where n is the number of spatial
dimensions) iterations over the array are needed so it should be
faster. For higher dimensions or for solving the Eikonal equation
maybe there is less advantage? Do you have any experience with this?

Regardless, the same initialization procedure is used for FMM and FSM.
So, improving the initialization is still a good idea.

Don't worry about changing the C++, if you can implement a method in
Python we can convert it to C++ or Cython later.

Also, have a look at this (now merged) PR:
https://github.com/scikit-fmm/scikit-fmm/pull/9
You may want to contact this person as they also seem to be using
scikit-fmm for topological optimization.

Jason


On Mon, Nov 14, 2016 at 3:23 PM, Neiferd, David John
<david....@wright.edu> wrote:
> Thanks for the additional info, I'll look into the 2nd order and higher
> order hopefully towards the end of this week. If it looks like it'll be
> useful for my type of problems, I might try adapting your code with the WENO
> method to see how that compares. My C/C++ skills are a lot weaker than
> Python since I don't use it on a regular basis, but I should be able to
> manage. I'll let you know how it goes.
>
> - John
> ________________________________
> From: Jason Furtney <jkfu...@gmail.com>
> Sent: Monday, November 14, 2016 4:08:10 PM
> To: Neiferd, David John
> Subject: Re: skfmm
>
> Thanks for the link to that paper, I will check it out. Maybe we
> should implement this method?
>
> By default scikit-fmm uses the second-order point update described in
> section 6.5 of the Sethian book. There is an optional keyword argument
> for using first-order point updates; this was done for backwards
> compatibility with another project.
>
> One weakness of the scikit-fmm implementation is that linear
> interpolation is used to establish the distance from the zero contour
> of phi to the initially "frozen" points (the points adjacent to the
> zero level-set). This is only first-order accurate and these errors
> get propagated by the marching algorithm. The "new-init" branch I have
> been working on addresses this problem.
>
> I have implemented the second-order initialization described in Salac,
> 2011 "The augmented fast marching method for level set
> reinitialization."
> https://urldefense.proofpoint.com/v2/url?u=http-3A__arxiv.org_abs_1111.6903&d=CwIFaQ&c=3buyMx9JlH1z22L_G5pM28wz_Ru6WjhVHwo-vpeS0Gk&r=G9mD5x8jN0uGw2T7G2liTJno5thoanUHnRxJNjFQvuc&m=CRpbOVK7HrBHPmYHvTNiO2mt2RiEJ0rApOvndAll_qI&s=Fj0wbIifhkTLuMDoylgz75yjrhPIsilCHiBTr_bBChA&e=
> . Specifically what
> is described in Section 2.2. This gives more accurate results. There
> are a few corner cases to resolve. I do not want to add SciPy as a
> dependency to scikit-fmm so I need to come up with a replacement for
> scipy.optimize.fsolve which this branch uses. I have done the same
> initialization for 3D but this needs more work.
>
> Thanks,
> Jason
>
> On Mon, Nov 14, 2016 at 2:11 PM, Neiferd, David John
> <david....@wright.edu> wrote:
>> I can clarify the documentation and if I setup a basic enough example I'll
>> include that also. Out of curiosity, do you have a reference I can look
>> at
>> for the second order method you are implementing? I haven't seen many
>> uses
>> of 2nd order methods, at least in topology literature, most of them use
>> higher order WENO schemes, i.e.
>> https://urldefense.proofpoint.com/v2/url?u=http-3A__www3.nd.edu_-7Eyzhang10_ZZQ-5FJSC.pdf&d=CwIFaQ&c=3buyMx9JlH1z22L_G5pM28wz_Ru6WjhVHwo-vpeS0Gk&r=G9mD5x8jN0uGw2T7G2liTJno5thoanUHnRxJNjFQvuc&m=CRpbOVK7HrBHPmYHvTNiO2mt2RiEJ0rApOvndAll_qI&s=B_k-DNNhbSb9SUNsw9lhvMHlDgx-mDycdu14f50gPaA&e=
>> .
>
>>
>>
>> ________________________________
>> From: Jason Furtney <jkfu...@gmail.com>
>> Sent: Friday, November 11, 2016 10:08:10 AM
>> To: Neiferd, David John
>> Subject: Re: skfmm
>>
>> Thanks and yes your description is correct. If you have some spare
>> effort it would be great if you could clarify the documentation and/or
>> provide another extension velocities example to help the next person
>> who comes along with the same question. Let me know how you get on
>> with the second-order initialization branch, I am happy to help move
>> this along.
>>
>> On Fri, Nov 11, 2016 at 9:00 AM, John Neiferd <david....@wright.edu>
>> wrote:
>>> Thanks for the swift reply I really appreciate it. So if I understand
>>> you
>>> correctly, as long as the points along the zero contour in the "speed"
>>> array
>>> are correct, the rest of the values in the "speed" array are arbitrary,
>>> correct?
>>>
>>> Also, thanks for the link to the in development branch, I'll check that
>>> out.
>>> Most of my experience has been in density based level set methods, but
>>> I'm
>>> moving into level-set based topology methods now as it seems more easily
>>> applicable to multidisciplinary system. I still have a lot to learn
>>> about
>>> the methods, but skfmm is a great starting point.
>>>
>>> - John
>>>
>>>
>>>
>>> On 11/11/2016 09:49 AM, Jason Furtney wrote:
>>>>
>>>> Dear John,
>>>>
>>>> Thanks for your email and for your interest in scikit-fmm.
>>>>
>>>> As I understand it, the speed input to the extsnsion_velocities
>>>> function is the velocity of the interface in the local normal
>>>> direction only on the zero contour of phi. The value of the speed
>>>> field everywhere else is discarded. The extension velocity is defined
>>>> in the rest of the domain by carrying the speeds defined on the zero
>>>> contour out into the domain.
>>>>
>>>> If you have not already check out "Level Set Methods and Fast Marching
>>>> Methods" by James Sethian as there should be a better description
>>>> there. The method implemented in scikit-fmm is described in Chapter
>>>> 11.
>>>>
>>>> I am (slowly) working on a branch which uses a second order accurate
>>>> method to initialize the zero level set. This may be helpful for your
>>>> application. The branch only works for 2d, regular grids and I am not
>>>> sure if the extension velocities are currently working in this branch.
>>>>
>>>>
>>>>
>>>> https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_scikit-2Dfmm_scikit-2Dfmm_tree_new-2Dinit&d=CwIBaQ&c=3buyMx9JlH1z22L_G5pM28wz_Ru6WjhVHwo-vpeS0Gk&r=G9mD5x8jN0uGw2T7G2liTJno5thoanUHnRxJNjFQvuc&m=HOL7JhZnzd3yDK8dt663iUcCI8IptYF8G1UNcAiDtzE&s=Lop8c41I7JJBs0ocgcdK62VUKpIuPVkxNX2qO7SgxLY&e=
>>>>
>>>> I hope that helps, let me know if you have any other questions.
>>>>
>>>> Jason
>>>>
>>>> On Thu, Nov 10, 2016 at 7:33 PM, John Neiferd <david....@wright.edu>
>>>> wrote:
>>>>>
>>>>> Hello Dr. Furtney,
>>>>>
>>>>> I had a question regarding your skfmm module. I'm trying to use the
>>>>> level-set method in some of my topology optimization research and I
>>>>> need
>>>>> to
>>>>> compute the extension velocities over the entire domain. I have my
>>>>> level-set function, phi, defined. Additionally, I have my velocities
>>>>> define
>>>>> at the zero-level contour. However, I don't quite understand what the
>>>>> "speed" input is that you require. I've looked through your
>>>>> documentation
>>>>> and the extension_velocities example, but still don't quite understand
>>>>> it.
>>>>> The speed is the magnitude of the velocity at the zero contour correct,
>>>>> but
>>>>> how can I compute the speed at the remaining points in the entire
>>>>> domain
>>>>> if
>>>>> the velocity is not know at those points?
>>>>>
>>>>> Thanks,
>>>>> John N.
>>>
>>>