Skip to content

Commit

Permalink
Updated README
Browse files Browse the repository at this point in the history
  • Loading branch information
@patnr
patnr committed Dec 7, 2018
1 parent 0a48bb7 commit 09614b9
Showing 1 changed file with 6 additions and 66 deletions.
72 changes: 6 additions & 66 deletions README.md
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
The fast marching (FM) method calculates first arrival time of a seismic wave travelling through various geological layers. The results serve as a kind of preview and comparison to the results obtained by the more comprehensive full-wave model. I implemented this to help out my dad at his work in the spring of 2011. The existing implementation on [file exchange](http://www.mathworks.com/matlabcentral/fileexchange/24531-accurate-fast-marching) only allowed for grids where the space between nodes is the same in all dimensions. This is clearly impractical if, for example, you want to model a geography that is complex in the horizontal plane but uniform in the vertical dimension. Supporting this feature was the client's main reason for requesting the project.
## The fast marching (FM) method
calculates first arrival time of a seismic wave travelling through various geological layers. The results serve as a kind of preview and comparison to the results obtained by the more comprehensive full-wave model. I implemented this to help out my dad at his work in the spring of 2011. The existing implementation on [file exchange](http://www.mathworks.com/matlabcentral/fileexchange/24531-accurate-fast-marching) only allowed for grids where the space between nodes is the same in all dimensions. This is clearly impractical if, for example, you want to model a geography that is complex in the horizontal plane but uniform in the vertical dimension. Supporting this feature was the client's main reason for requesting the project.



Expand All @@ -15,20 +16,13 @@ Both 1st order and 2nd order finite difference schemes are available, and warnin


## Installation
Just download. Test with one of the `TS_` scripts.
Just download. Test with one of the `TS_` scripts. Output is described below.


Note: Multicore mode can be had using [this multicore module](http://www.mathworks.com/matlabcentral/fileexchange/13775-multicore-parallel-processing-on-multiple-cores).



Note: running the code requires a newer version of Matlab, as it relies on some newer Matlab functionality and object oriented code. I'm unsure if the plotting calls are still valid after Matlab updated their graphics interface.

Note: The C++ code uses infinity as a number to avoid use of boolean switches. This works on the GNU compiler, but I don't think it's portable to many compilers.
Note: The C++ code uses infinity as a number to avoid use of boolean switches. This works on the GNU compiler, but I don't think it ports to all compilers. Alternatively, just stick to the Matlab implementations.



## Results
## Examples

The cooler output plots are further below.

Expand All @@ -37,67 +31,37 @@ If you drop a stone into the water, rings representing wavefronts propagate outw

![circle](./Pics/circle.jpg)



* * *

Here's the sort of thing the client might use this for. Let's say the earth below us consists of 4 distinct layers, each with its own propagation speed, F(x,z), illustrated below.


![strata-speedmap](./Pics/strata-speedmap.jpg)




This time, the wavefronts will not produce circular rings, or analytic solutions that are easy to find. However, the fast marching method is not deterred. All it needs is a specification of the speed map F(x,z) and it will give you T(x,z).



![strata](./Pics/strata.jpg)








* * *

What if we're in a 3D domain? The fast marching method is easily generalized, requiring only some adaptation of the finite difference scheme and the data storage. The plot below shows 9 cross sections (i.e. z=constant) of T(x,y,z).

![crosSections_solution](./Pics/crosSections_solution.jpg)



* * *

The implementation is flexible enough to handle trivial/exceptional cases. For example, here is the output from running the 3D fast marching method on a domain with two singleton dimensions (i.e. where two of the dimensions only contain one discretization step).



![bar](./Pics/bar.jpg)



* * *

This example illustrates that one can insert impenetrable "walls" in the domain, where the speed of the wave is defined to be zero. This is in many ways an exceptional case, as it may result in division by zero or zero*inf computations, and so it requires special attention. Also, the speed map has been bombarded with random noise, creating the noisy look of the plot.




![random](./Pics/random.jpg)



* * *






Now we're getting to the cool parts. Did you know the method can be used to solve labyrinths? Here's how. Make a speed map F(x,y) where F=0 wherever there's a wall and F=1 elsewhere. Put the source point at the end (exit) of the labyrinth, and run the fast marching method. Then the wavefront will propagate through the labyrinth, including cul-de-sacs, and produce a map of T(x,y) as a result again. This tells you how long it takes to reach any point in the labyrinth from the exit, and is shown by the color in the maze below.


Expand All @@ -112,12 +76,8 @@ How do you find the path through the labyrinth from the color map? The white lin

Cooler yet, the very same method can be used to figure out the best way to, say, move your bed through the maze that is your house! All you need is to make a 3D speed map model of the house, and set F=0 wherever there's a wall. However, since the bed is not a point, but an object with a finite size, you might want to add some complexity -- increasing the dimensionality of your model to 5D. The two new dimensions will be used to specify the orientation of the bed (elevation and azimuth). Now the speed map is a function of five variables: F(x,y,z,θ,ϕ). So, say you're in a narrow hallway. Then F=1 for most values of θ and ϕ, except when the bed would crash with the walls, in which case F=0. Also, you might want to slow down (by lowering F) when you're in a tight spot, to reflect the fact that progress will be slower. In the end, you run a ray tracer backwards, which will give you the fastest/optimal way of moving the bed out of the apartment.



![maze](./Pics/maze.jpg)



* * *

How about getting from A to B as fast as possible? Same principle. Propagate a wave front from a source point, taking into account different propagation speeds.
Expand All @@ -130,40 +90,20 @@ Unless you're James Bond, the speed of your boat will be around 20 knot on the o

Starting the ray tracer from, say, LA, will lead you back to Oslo through the white path -- the shortest path.



![LA1](./Pics/LA1.jpg)


![world_LA1](./Pics/world_LA1.jpg)

Same time-distance map T, different starting point for the ray tracer.



![world_Tokyo1](./Pics/world_Tokyo1.jpg)





This time, let's open up the Bering strait, and close the Suez canal. I.e. set the speed map pixels in those regions to F=0 and F=1. Run the fast marching method again. Then the ray tracer from Tokyo. This time we see that it's shorter to go north of Russia than south of Africa, though you might get slightly cold on the way. Note that I'm not sure what kind of [world map projection](http://en.wikipedia.org/wiki/Map_projection) this is, but it probably distorts global distances.



![world_Tokyo2](./Pics/world_Tokyo2.jpg)





What if we close the Bering strait and the Suez canal? Then we gotta go by the Cape of Good Hope, like the spice traders of yore.



![world_Tokyo3](./Pics/world_Tokyo3.jpg)



* * *

These are some of the possible applications of the fast marching method. There are of course infinitely (countably?) more.

0 comments on commit 09614b9

Please sign in to comment.