4

I am trying to work with a simulation of Brain Tumor growth and I was fortunate to get a very great example at http://community.wolfram.com/. However, I need clarifications on some of the code. I want to know what it does and how it works. Here are the piece that I do not understand:

1.

diffcoeff = ListInterpolation[ImageData[img3], InterpolationOrder -> 3]

2.

 boundaries = {-y, y - 1, -x, x - 1};
 Ω = ImplicitRegion[And @@ (# <= 0 & /@ boundaries), {x, y}];

3.

 sols = NDSolveValue[{{Div[
   1./500.*(diffcoeff[798.*x, 654*y])^4* Grad[u[t, x, y], {x, y}], {x, y}] - D[u[t, x, y], t] + 0.025*u[t, x, y] == NeumannValue[0., x >= 1. || x <= 0. || y <= 0. || y >= 1.]}, 
   {u[0, x, y] == Exp[-1000. ((x - 0.6)^2 + (y - 0.6)^2)]}}, 
   u, 
   {x, y}∈ Ω, 
   {t, 0, 20}, 
   Method -> {"FiniteElement", 
     "MeshOptions" -> {
        "BoundaryMeshGenerator" -> "Continuation", 
        MaxCellMeasure -> 0.002
        }
      }
   ]

I really would like to understand each element of these code pieces (and probably the rationale behind their use). Your kind and expert contributions would be greatly appreciated. Since the web page does not load anymore, I have the pdf copy I earlier saved, here!.

Community
  • 1
Dean
  • 655
  • 5
  • 9

1 Answers1

8

  1. This takes the image's grayscale data and interpolates it with a piecewise bicubic spline function with C^1-continuity. It will be used as diffusivity for the partial differential equation in 3. Interpolation means that diffcoeff[i,j] for integer i and j has the value of the pixel ImageData[img3][[i,j]].

  2. This defines the unit square as region. The partial differential equation in 3. will be solved on this domain.

  3. This solves a one-component reaction-diffusion equation with a Gaussian as initial condition and subject to homogeneous Neumann boundary conditions by means of the finite element method.

With respect to the meaning of the numbers 798 and 654 in diffcoeff[798.*x, 654*y]: I guess that Dimensions[ImageData[img3]] equals {798,654}. The image is supposed to be mapped onto the unit square, so points in the unit square have to be scaled before feeding them to diffcoeff. Whether it is a good idea to use different scalings for the two coordinate directions and why this is supposed to be a good model for cancer growth is beyond me.

For details of each of the employed built-in symbols, please refer to its documentation.

Henrik Schumacher
  • 106,770
  • 7
  • 179
  • 309
  • 1
    Frankly, exposition of the second point in this code is overly complicated. It might make some sense if Ω has been defined to be, say, a hexagon in some other case or that this code partially originates from pre-v10 era, but in the case of a unit square... you can just say Ω = Rectangle[]. (RegionEqual[With[{boundaries = {-y, y - 1, -x, x - 1}}, ImplicitRegion[And @@ (# <= 0 & /@ boundaries), {x, y}]], Rectangle[]] evaluates to True.) – kirma Sep 02 '18 at 03:55
  • 1
    Indeed. I also wondered about that as well as how the Neumann conditons are declared. – Henrik Schumacher Sep 02 '18 at 06:23
  • 2
    If you were thinking of using True in the Neumann condition this may work in this case but it would not work for a geometry that had internal boundaries. It is more general to explicitly specify the boundaries on which a BC is active, similar to this – user21 Sep 03 '18 at 05:30
  • Thank you guys for your excellent comments. For 1, I still do not understand the physiological basis of the Interpolation function. I want to know the relationship between image data and diffusivity. For 2, I am almost okay. For 3, the basis of feeding 798 and 654 into the reaction diffusion equation is still beyond me. I need to be able to explain the mathematics of the codes! – Dean Sep 04 '18 at 03:05