edu.mines.jtk.interp.TimeMarker2 Maven / Gradle / Ivy
Show all versions of edu-mines-jtk Show documentation
/****************************************************************************
Copyright 2009, Colorado School of Mines and others.
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
****************************************************************************/
package edu.mines.jtk.interp;
import java.util.Random;
import java.util.concurrent.*;
import java.util.concurrent.atomic.AtomicInteger;
import edu.mines.jtk.dsp.Tensors2;
import edu.mines.jtk.util.Parallel;
import static edu.mines.jtk.util.ArrayMath.*;
/**
* A time and closest-point transform for 2D anisotropic eikonal equations.
* Transforms an array of times and marks for known samples into an array
* of times and marks for all samples. Known samples are those for which
* times are zero, and times and marks for known samples are not modified.
*
* Times for unknown samples are computed by solving an anisotropic eikonal
* equation grad(t) dot W grad(t) = 1, where W denotes a positive-definite
* (velocity-squared) metric tensor field. The solution times t represent
* the traveltimes from one known sample to all unknown samples. Separate
* solution times t are computed for each known sample that has at least
* one unknown neighbor. (Such a known sample is sometimes called a
* "source point.") The output time for each sample is the minimum time
* computed in this way for all such known samples. Therefore, the times
* output for unknown samples are the traveltimes to the closest known
* samples, where "closest" here means least time, not least distance.
*
* As eikonal solutions t are computed for each known sample, the mark
* for that known sample is used to mark all unknown samples for which
* the solution time is smaller than the minimum time computed so far.
* If marks for known samples are distinct, then output marks for unknown
* samples indicate which known sample is closest. In this way output
* marks can represent a sampled Voronoi diagram, albeit one that has
* been generalized by replacing distance with time.
*
* This transform uses an iterative sweeping method to solve for times.
* Iterations are similar to those described by Jeong and Whitaker (2007).
* Computational complexity is O(M log K), where M is the number of
* unknown (missing) samples and K is the number of known samples.
* @author Dave Hale, Colorado School of Mines
* @version 2009.01.06
*/
public class TimeMarker2 {
/**
* Type of concurrency used by this transform. Default is PARALLEL.
*/
public enum Concurrency {
PARALLELX,
PARALLEL,
SERIAL
}
/**
* Constructs a time marker for the specified tensor field.
* @param n1 number of samples in 1st dimension.
* @param n2 number of samples in 2nd dimension.
* @param tensors velocity-squared tensors.
*/
public TimeMarker2(int n1, int n2, Tensors2 tensors) {
init(n1,n2,tensors);
}
/**
* Sets the tensors used by this time marker.
* @param tensors the tensors.
*/
public void setTensors(Tensors2 tensors) {
_tensors = tensors;
}
/**
* Sets the type of concurrency used to solve for times.
* The default concurrency is parallel.
* @param concurrency the type of concurrency.
*/
public void setConcurrency(Concurrency concurrency) {
_concurrency = concurrency;
}
/**
* Transforms the specified array of times and marks.
* Known samples are those for which times are zero, and times
* and marks for these known samples are used to compute times
* and marks for unknown samples.
* @param times input/output array of times.
* @param marks input/output array of marks.
*/
public void apply(float[][] times, int[][] marks) {
// Initialize all unknown times to infinity.
for (int i2=0; i2<_n2; ++i2) {
for (int i1=0; i1<_n1; ++i1) {
if (times[i2][i1]!=0.0f)
times[i2][i1] = INFINITY;
}
}
// Indices of known samples in random order.
short[][] kk = indexKnownSamples(times);
short[] k1 = kk[0];
short[] k2 = kk[1];
shuffle(k1,k2);
int nk = k1.length;
// Array for the eikonal solution times.
float[][] t = new float[_n2][_n1];
// Active list of samples used to compute times.
ActiveList al = new ActiveList();
// For all known samples, ...
for (int ik=0; ik_a.length)
growTo(2*(_n+al._n));
int n = al._n;
for (int i=0; i.
private static class ShortStack {
void push(int k) {
if (_n==_a.length) {
short[] a = new short[2*_n];
System.arraycopy(_a,0,a,0,_n);
_a = a;
}
_a[_n++] = (short)k;
}
short pop() {
return _a[--_n];
}
int size() {
return _n;
}
void clear() {
_n = 0;
}
boolean isEmpty() {
return _n==0;
}
short[] array() {
short[] a = new short[_n];
System.arraycopy(_a,0,a,0,_n);
return a;
}
private int _n = 0;
private short[] _a = new short[2048];
}
/*
* Returns arrays of indices of known samples with times zero.
* Includes only known samples adjacent to at least one unknown sample.
* (Does not include known samples surrounded by other known samples.)
*/
private short[][] indexKnownSamples(float[][] times) {
ShortStack ss1 = new ShortStack();
ShortStack ss2 = new ShortStack();
for (int i2=0; i2<_n2; ++i2) {
for (int i1=0; i1<_n1; ++i1) {
if (times[i2][i1]==0.0f) {
for (int k=0; k<4; ++k) {
int j1 = i1+K1[k]; if (j1<0 || j1>=_n1) continue;
int j2 = i2+K2[k]; if (j2<0 || j2>=_n2) continue;
if (times[j2][j1]!=0.0f) {
ss1.push(i1);
ss2.push(i2);
break;
}
}
}
}
}
short[] i1 = ss1.array();
short[] i2 = ss2.array();
return new short[][]{i1,i2};
}
/*
* Randomly shuffles the specified arrays of indices.
*/
private static void shuffle(short[] i1, short[] i2) {
int n = i1.length;
Random r = new Random(314159); // constant seed for consistency
short ii;
for (int i=n-1; i>0; --i) {
int j = r.nextInt(i+1);
ii = i1[i]; i1[i] = i1[j]; i1[j] = ii;
ii = i2[i]; i2[i] = i2[j]; i2[j] = ii;
}
}
/*
* Solves for times by sequentially processing each sample in active list.
*/
private void solveSerial(
ActiveList al,
float[][] t, int m,
float[][] times, int[][] marks)
{
float[] d = new float[3];
ActiveList bl = new ActiveList();
//int ntotal = 0;
while (!al.isEmpty()) {
//al.shuffle(); // demonstrate that solution depends on order
int n = al.size();
//ntotal += n;
for (int i=0; i-->