graph.cc

二维射线追踪地震层析成像

/*
 * graph.cc - graph solver implementation
 *
 * Jun Korenaga, MIT/WHOI
 * January 1999
 */

#include <iostream>
#include <algorithm>
#include <cstdlib>
#include <ctime>
#include "smesh.h"
#include "graph.h"

GraphSolver2d::GraphSolver2d(const SlownessMesh2d& m, int xorder, int zorder)
    : smesh(m), nnodes(m.numNodes()),
      is_solved(false), is_refl_solved(false), solve_refl_downward(false),
      fs_xorder(xorder), fs_zorder(zorder),
      cmp(ttime), B(cmp),
      cmp_down(ttime_down), B_down(cmp_down),
      cmp_up(ttime_up), B_up(cmp_up),
      local_search_radius(10.0), xmin(smesh.xmin()), xmax(smesh.xmax()),
      refine_if_long(false)
{
    prev_node.resize(nnodes);
    ttime.resize(nnodes);
}

void GraphSolver2d::limitRange(double x1, double x2)
{
    if (x1>x2) error("GraphSolver2d::limitRange - invalid input");
    xmin=x1-local_search_radius;
    xmax=x2+local_search_radius;
}

void GraphSolver2d::delimitRange()
{
    xmin = smesh.xmin();
    xmax = smesh.xmax();
}

void GraphSolver2d::refineIfLong(double x)
{
    refine_if_long = true;
    crit_len = x;
}

void GraphSolver2d::solve(const Point2d& s)
{
    double ttime_inf = 1e30;

    src = s;
    if (smesh.inWater(src) || smesh.inAir(src)){
	error("GraphSolver2d::solve - currently unsupported source configuration detected.\n");
    }
    int Nsrc = smesh.nearest(src);
    
    // initialization
    B.resize(0);
    prev_node(Nsrc) = Nsrc;
    ttime(Nsrc) = 0.0;
    B.push(Nsrc); 
	
    C.resize(0);
    for (int i=1; i<=nnodes; i++){
	if (i!=Nsrc){
	    double x = smesh.nodePos(i).x();
	    if (x>=xmin && x<=xmax){
		C.push_back(i);
		ttime(i) = ttime_inf;
		prev_node(i) = Nsrc;
	    }
	}
    }

    while(B.size() != 0){
	// find a node with minimum traveltime from B,
	// (and transfer it to A)
	int N0=B.top();
	B.pop();

	// form a forward star
	ForwardStar2d fstar(smesh.nodeIndex(N0), fs_xorder, fs_zorder);

	// loop over [FS and B] to update traveltime
	heapBrowser pB=B.begin();
	while(pB!=B.end()){
	    if (fstar.isIn(smesh.nodeIndex(*pB))){
		double tmp_ttime=ttime(N0)+smesh.calc_ttime(N0,*pB);
		double orig_ttime=ttime(*pB);
 		if (orig_ttime>tmp_ttime){
		    ttime(*pB) = tmp_ttime;
		    prev_node(*pB) = N0;
		    B.promote(pB);
		}
	    }
	    pB++;  
	} 
	
	// loop over [FS and C] to calculate traveltime
	// and transfer them to B
	nodeIterator pC=C.begin();
	nodeIterator epC;
	while(pC!=C.end()){
	    if (fstar.isIn(smesh.nodeIndex(*pC))){
		ttime(*pC) = ttime(N0)+smesh.calc_ttime(N0,*pC);
		prev_node(*pC) = N0;
		B.push(*pC);

		epC = pC;
		pC++;
		C.erase(epC);
	    }else{ 
		pC++;
	    }
	}
    }
    is_solved = true;
}

void GraphSolver2d::do_refl_downward(){ solve_refl_downward = true; }

void GraphSolver2d::solve_refl(const Point2d& s, const Interface2d& itf)
{
    if (!is_solved && !solve_refl_downward)
	error("GraphSolver2d::not ready for solve_refl()");
    
    double ttime_inf = 1e30;

    if (prev_node_down.size() != nnodes) prev_node_down.resize(nnodes);
    if (ttime_down.size() != nnodes) ttime_down.resize(nnodes);
    if (prev_node_up.size() != nnodes) prev_node_up.resize(nnodes);
    if (ttime_up.size() != nnodes) ttime_up.resize(nnodes);
    itfp = &itf;
    Array1d<int> itfnodes;
    smesh.nearest(itf,itfnodes);

    if (solve_refl_downward){
	// initialization (downgoing)
	src = s;
	if (smesh.inWater(src) || smesh.inAir(src)){
	    error("GraphSolver2d::solve_refl - currently unsupported source configuration detected.\n");
	}
	int Nsrc = smesh.nearest(src);
    
	B_down.resize(0);
	prev_node_down(Nsrc) = Nsrc;
	ttime_down(Nsrc) = 0.0;
	B_down.push(Nsrc); 
	
	C.resize(0);
	for (int i=1; i<=nnodes; i++){
	    if (i!=Nsrc){
		double x = smesh.nodePos(i).x();
		double z = smesh.nodePos(i).y();
		double itfz = itf.z(x);
		if (x>=xmin && x<=xmax){
		    bool is_ok = false;
		    if (z<=itfz){
			is_ok = true;
		    }else{
			for (int n=1; n<=itfnodes.size(); n++){
			    if (itfnodes(n)==i){
				is_ok = true; break;
			    }
			}
		    }
		    if (is_ok){
			C.push_back(i);
			ttime_down(i) = ttime_inf;
			prev_node_down(i) = Nsrc;
		    }
		}
	    }
	}

	while(B_down.size() != 0){
	    // find a node with minimum traveltime from B,
	    // (and transfer it to A)
	    int N0=B_down.top();
	    B_down.pop();

	    // form a forward star
	    ForwardStar2d fstar(smesh.nodeIndex(N0), fs_xorder, fs_zorder);

	    // loop over [FS and B] to update traveltime
	    heapBrowser pB=B_down.begin();
	    while(pB!=B_down.end()){
		if (fstar.isIn(smesh.nodeIndex(*pB))){
		    double tmp_ttime=ttime_down(N0)+smesh.calc_ttime(N0,*pB);
		    double orig_ttime=ttime_down(*pB);
		    if (orig_ttime>tmp_ttime){
			ttime_down(*pB) = tmp_ttime;
			prev_node_down(*pB) = N0;
			B_down.promote(pB);
		    }
		}
		pB++;  
	    } 
	
	    // loop over [FS and C] to calculate traveltime
	    // and transfer them to B
	    nodeIterator pC=C.begin();
	    nodeIterator epC;
	    while(pC!=C.end()){
		if (fstar.isIn(smesh.nodeIndex(*pC))){
		    ttime_down(*pC) = ttime_down(N0)+smesh.calc_ttime(N0,*pC);
		    prev_node_down(*pC) = N0;
		    B_down.push(*pC);

		    epC = pC;
		    pC++;
		    C.erase(epC);
		}else{ 
		    pC++;
		}
	    }
	}
    }else{
	// use refraction's graph solution
	prev_node_down = prev_node;
	ttime_down = ttime;
    }
    
    // initialization (upgoing)
    B_up.resize(0);
    C.resize(0);
    for (int i=1; i<=nnodes; i++){
	bool is_itf = false;
	double x=smesh.nodePos(i).x();
	double z=smesh.nodePos(i).y();
	double itfz=itf.z(x);
	if (x>=xmin && x<=xmax){
	    for (int n=1; n<=itfnodes.size(); n++){
		if (itfnodes(n)==i){
		    is_itf = true;
		    prev_node_up(i) = -1;
		    ttime_up(i) = ttime_down(i);
		    B_up.push(i);
		    break;
		}
	    }
	    if (!is_itf && z<=itfz){ // take nodes above the reflector
		ttime_up(i) = ttime_inf;
		prev_node_up(i) = 0;
		C.push_back(i);
	    }
	}
    }

    while(B_up.size() != 0){
	// find a node with minimum traveltime from B,
	// (and transfer it to A)
	int N0=B_up.top();
	B_up.pop();

	// form a forward star
	ForwardStar2d fstar(smesh.nodeIndex(N0), fs_xorder, fs_zorder);

	// loop over [FS and B] to update traveltime
	heapBrowser pB=B_up.begin();
	while(pB!=B_up.end()){
	    if (fstar.isIn(smesh.nodeIndex(*pB))){
		double tmp_ttime=ttime_up(N0)+smesh.calc_ttime(N0,*pB);
		double orig_ttime=ttime_up(*pB);
 		if (orig_ttime>tmp_ttime){
		    ttime_up(*pB) = tmp_ttime;
		    prev_node_up(*pB) = N0;
		    B_up.promote(pB);
		}
	    }
	    pB++;  
	} 
	
	// loop over [FS and C] to calculate traveltime
	// and transfer them to B
	nodeIterator pC=C.begin();
	nodeIterator epC;
	while(pC!=C.end()){
	    if (fstar.isIn(smesh.nodeIndex(*pC))){
		ttime_up(*pC) = ttime_up(N0)+smesh.calc_ttime(N0,*pC);
		prev_node_up(*pC) = N0;
		B_up.push(*pC);
		epC = pC;
		pC++;
		C.erase(epC);
	    }else{ 
		pC++;
	    }
	}
    }
    
    is_refl_solved = true;
}

void GraphSolver2d::pickPath(const Point2d& rcv, Array1d<Point2d>& path) const
{
    if (!is_solved) error("GraphSolver2d::not yet solved");

    list<Point2d> tmp_path;
    tmp_path.push_back(rcv);

    int cur = smesh.nearest(rcv);
    int prev;
    int nadd=0;
    while ((prev=prev_node(cur)) != cur){
	Point2d p = smesh.nodePos(prev);
	tmp_path.push_back(p); nadd++;
	cur = prev;
    }

    if (nadd>0){
	tmp_path.pop_back();
    }
    if (nadd==1){
	// insert mid point for later bending refinement
	tmp_path.push_back(0.5*(src+rcv));