graph.cc

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

    }
    tmp_path.push_back(src);

    if (refine_if_long){
	list<Point2d>::iterator pt=tmp_path.begin();
	Point2d prev_p = *pt;
	pt++; // skip the first point
	while(pt!=tmp_path.end()){
	    double seg_len = prev_p.distance(*pt);
	    if (seg_len>crit_len){
		int ndiv = int(seg_len/crit_len+1.0);
		double frac = 1.0/ndiv;
		for (int j=1; j<ndiv; j++){
		    double ratio = frac*j;
		    Point2d new_p = (1.0-ratio)*prev_p+ratio*(*pt);
		    tmp_path.insert(pt,new_p);
		}
	    }
	    prev_p = *pt++;
	}
    }

    path.resize(tmp_path.size());
    list<Point2d>::const_iterator pt=tmp_path.begin();
    int i=1;
    while(pt!=tmp_path.end()){
	path(i++) = *pt++;
    }
}

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

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

    int cur = smesh.nearest(rcv);
    findRayEntryPoint(rcv,cur,ttime);

    Point2d entryp=smesh.nodePos(cur);
    tmp_path.push_back(entryp);
    tmp_path.push_back(entryp);
    tmp_path.push_back(entryp);
    i0 = 2; i1 = 4;
    
    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+entryp));
    }
    tmp_path.push_back(src);

    if (refine_if_long){
	list<Point2d>::iterator pt=tmp_path.begin();
	pt++; pt++; pt++; // now at the entry point
	Point2d prev_p = *pt;
	pt++; // skip the entry point
	while(pt!=tmp_path.end()){
	    double seg_len = prev_p.distance(*pt);
	    if (seg_len>crit_len){
		int ndiv = int(seg_len/crit_len+1.0);
		double frac = 1.0/ndiv;
		for (int j=1; j<ndiv; j++){
		    double ratio = frac*j;
		    Point2d new_p = (1.0-ratio)*prev_p+ratio*(*pt);
		    tmp_path.insert(pt,new_p);
		}
	    }
	    prev_p = *pt++;
	}
    }
    
    path.resize(tmp_path.size());
    list<Point2d>::const_iterator pt=tmp_path.begin();
    int i=1;
    while(pt!=tmp_path.end()){
	path(i++) = *pt++;
    }
}

void GraphSolver2d::pickReflPath(const Point2d& rcv, Array1d<Point2d>& path,
				 int& ir0, int& ir1) const
{
    if (!is_refl_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_up(cur)) != -1){
	Point2d p = smesh.nodePos(prev);
	tmp_path.push_back(p); nadd++; 
	cur = prev;
    }
    double reflx = smesh.nodePos(cur).x();
    double reflz = itfp->z(reflx);
    Point2d onrefl(reflx,reflz);

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

    tmp_path.push_back(onrefl); ir0 = tmp_path.size();
    tmp_path.push_back(onrefl);
    tmp_path.push_back(onrefl); ir1 = tmp_path.size();
    
    nadd=0;
    while ((prev=prev_node_down(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));
    }
    tmp_path.push_back(src);

    if (refine_if_long){
	int old_i=1, nadd_pre=0;
	list<Point2d>::iterator pt=tmp_path.begin();
	Point2d prev_p = *pt;
	pt++; old_i++; // skip the first point
	while(pt!=tmp_path.end()){
	    double seg_len = prev_p.distance(*pt);
	    if (seg_len>crit_len){
		int ndiv = int(seg_len/crit_len+1.0);
		double frac = 1.0/ndiv;
		for (int j=1; j<ndiv; j++){
		    double ratio = frac*j;
		    Point2d new_p = (1.0-ratio)*prev_p+ratio*(*pt);
		    tmp_path.insert(pt,new_p);
		    if (old_i<=ir0) nadd_pre++;
		}
	    }
	    prev_p = *pt++; old_i++;
	}
	ir0 += nadd_pre;
	ir1 += nadd_pre;
    }
    
    path.resize(tmp_path.size());
    list<Point2d>::const_iterator pt=tmp_path.begin();
    int i=1;
    while(pt!=tmp_path.end()){
	path(i++) = *pt++;
    }
}

void GraphSolver2d::pickReflPathThruWater(const Point2d& rcv, Array1d<Point2d>& path,
					  int& i0, int& i1, int& ir0, int& ir1) const
{
    if (!is_refl_solved) error("GraphSolver2d::not yet solved");

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

    int cur = smesh.nearest(rcv);
    findRayEntryPoint(rcv,cur,ttime_up);

    Point2d entryp=smesh.nodePos(cur);
    tmp_path.push_back(entryp);
    tmp_path.push_back(entryp);
    tmp_path.push_back(entryp);
    i0 = 2; i1 = 4;
    
    int prev;
    int nadd=0;
    while ((prev=prev_node_up(cur)) != -1){
	Point2d p = smesh.nodePos(prev);
	tmp_path.push_back(p); nadd++; 
	cur = prev;
    }
    double reflx = smesh.nodePos(cur).x();
    double reflz = itfp->z(reflx);
    Point2d onrefl(reflx,reflz);

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

    tmp_path.push_back(onrefl); ir0 = tmp_path.size();
    tmp_path.push_back(onrefl);
    tmp_path.push_back(onrefl); ir1 = tmp_path.size();
    
    nadd=0;
    while ((prev=prev_node_down(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));
    }
    tmp_path.push_back(src);

    if (refine_if_long){
	int nadd_pre=0;
	list<Point2d>::iterator pt=tmp_path.begin();
	pt++; pt++; pt++; // now at the entry point
	Point2d prev_p = *pt;
	pt++; // skip the entry point
	int old_i = 5;
	while(pt!=tmp_path.end()){
	    double seg_len = prev_p.distance(*pt);
	    if (seg_len>crit_len){
		int ndiv = int(seg_len/crit_len+1.0);
		double frac = 1.0/ndiv;
		for (int j=1; j<ndiv; j++){
		    double ratio = frac*j;
		    Point2d new_p = (1.0-ratio)*prev_p+ratio*(*pt);
		    tmp_path.insert(pt,new_p);
		    if (old_i<=ir0) nadd_pre++;
		}
	    }
	    prev_p = *pt++; old_i++;
	}
	ir0 += nadd_pre;
	ir1 += nadd_pre;
    }

    path.resize(tmp_path.size());
    list<Point2d>::const_iterator pt=tmp_path.begin();
    int i=1;
    while(pt!=tmp_path.end()){
	path(i++) = *pt++;
    }
}

void GraphSolver2d::findRayEntryPoint(const Point2d& rcv, int& cur, const Array1d<double>& nodetime) const
{
    Point2d p0 = smesh.nodePos(cur);
    double pwater = smesh.atWater();
    double ttime0 = pwater*rcv.distance(p0)+nodetime(cur);

    // local grid search for best ray entry point (+/- 10km radius)
    Point2d pmin(rcv.x()-local_search_radius,rcv.y());
    Point2d pmax(rcv.x()+local_search_radius,rcv.y()); 
    Index2d indmin = smesh.nodeIndex(smesh.nearest(pmin)); 
    Index2d indmax = smesh.nodeIndex(smesh.nearest(pmax)); 
    for (int i=indmin.i()+1; i<indmax.i(); i++){
	int inode = smesh.nodeIndex(i,1);
	Point2d p = smesh.nodePos(inode);
	double tmp_ttime = pwater*rcv.distance(p)+nodetime(inode);
	if (tmp_ttime < ttime0){
	    ttime0 = tmp_ttime;
	    cur = inode; 
	}
    }
}

void GraphSolver2d::printPath(ostream& os) const
{
    for (int cur=1; cur<=nnodes; cur++){
	int prev = prev_node(cur);
	Point2d cur_p = smesh.nodePos(cur);
	Point2d prev_p = smesh.nodePos(prev);
	    os << ">\n"
	       << cur_p.x() << " "  << cur_p.y() << "\n"
	       << prev_p.x() << " " << prev_p.y() << '\n';
    }
}

//
// forward star
//
ForwardStar2d::ForwardStar2d(const Index2d& id, int ix, int iz)
    : orig(id)
{
    if (ix>0 && iz>0){ xorder=ix; zorder=iz; }
    else{ error("ForwardStar2d::non-positive order detected."); }
}

bool ForwardStar2d::isIn(const Index2d& node) const
{
    int diff;
    int di = (diff=orig.i()-node.i()) > 0 ? diff : -diff;
    int dk = (diff=orig.k()-node.k()) > 0 ? diff : -diff;

    if (di>xorder || dk>zorder) return false;
    if (di==1 || dk==1) return true;
    if (di==0 || dk==0) return false; // NB: (0,1),(1,0) are already returned true.

    // now di && dk must be >=2 
    int mod = dk>di ? dk%di : di%dk;
    return mod==0 ? false : true;
}