seisinversionproceeds.cpp

石油勘探专业算法

	int line_id1;
	int cdp1;
	float* values = new float[num];
	int * cdps = new int[num];
	int* lids = new int[num];
	float temp = 0.0;
	
	for(i = 0;i<num;i++) {
		curves[i]->GetLineAndCDP(line_id1,cdp1);
		GetCurveT0(hor_cur_pos1,hor_cur_pos2,line_id1,cdp1,par);
		switch(par->_type) {
			case NonRestrictLayer: {
				t0 = t01_max + (layer_idx + 1) * delta_t0;
				break;
			}
			case ParallelTopHorizon: {
				t0 = hor_cur_pos1 + (layer_idx + 1) * delta_t0;
				break;
			}
			case ParallelBottomHorizon: {
				t0 = hor_cur_pos2 - (par->_layer_num - layer_idx) * delta_t0;
				break;
			}
			case IsohedralBetweenBorders: {
				t0 = hor_cur_pos1 + (layer_idx + 1) * delta_t0 / (t02_max - t01_max ) *
					(hor_cur_pos2 - hor_cur_pos1);
				break;
			}
		}
		temp = curves[i]->GetValue(t0);
		if(temp > 0){
			values[j] = temp;
			cdps[j] = cdp1;
			lids[j] = line_id1;
			j++;
		}
	}
	if(j > 0){
		interpolate_values = new float[j];
		interpolate_cdps = new int[j];
		interpolate_lids = new int[j];
		memcpy(interpolate_values,values,j * sizeof(float));
		memcpy(interpolate_cdps,cdps,j * sizeof(int));
		memcpy(interpolate_lids,lids,j * sizeof(int));
	}
	delete[] values;
	delete[] cdps;
	delete[] lids;
	return j;
}

float SeisInversion::GetWeighValue(int lid,int cdp,float* interpolate_values,int* interpolate_cdps,int* interpolate_lids,
		int interpolate_values_num)
{
	float value = 0.0;
	float dd = 0;
	float val = -1.0;
	float temp1 = 0.0,temp2 = 0.0;

	if(interpolate_values_num == 1) return interpolate_values[0];
	for(int i = 0;i < interpolate_values_num;i++) {
		dd = sqrt((interpolate_lids[i] - lid) * (interpolate_lids[i] - lid) +
			(interpolate_cdps[i] - cdp) * (interpolate_cdps[i] - cdp));
		if(dd < 0.01) {
			return interpolate_values[i];
		}
		temp1 = temp1 + interpolate_values[i] / (float)dd;
		temp2 = temp2 + 1.0 / (float)dd;
	}
	if( (temp1 > 0)&&(temp2 > 0) )
	{
		val = temp1 / temp2;
	}
	return val;
}

void SeisInversion::InterpolateOperation()
{
	time_t clock;

	time(&clock);
	char*time_str = NULL;
	time_str = ctime(&clock);
	cout<<"Interpolate Operation Begin Time: "<<time_str<<endl;
	GetHorizonSurfaceData();

	int i = 0,j = 0,k = 0,m = 0;
	float delta_t0 = 0.0;
	wURHorizonRestrictPars* hor_sestrict_par = NULL;
	int num = 0;
	wDBInversionCurveInSurvey** interpolate_curves = NULL;
	int interpolate_curves_num = 0;

	hor_sestrict_par = _seis_tpt->GetHorizonPars(num);
	_seis_tpt->Obtain3DTemplate()->ObtainSeis3D()->GetLineRange(_slid_proceed,_elid_proceed,_lid_inv);
	_seis_tpt->Obtain3DTemplate()->ObtainSeis3D()->GetCDPRange(_scdp_proceed,_ecdp_proceed,_cdp_inv);

	wDBTempHorDataVolume* hor_data_volume1 = NULL;
	hor_data_volume1 = new wDBTempHorDataVolume(_survey_id,_slid_proceed,_elid_proceed,_lid_inv,
			_scdp_proceed,_ecdp_proceed,_cdp_inv,"t0");
	wDBTempHorDataVolume* hor_data_volume2 = NULL;
	hor_data_volume2 = new wDBTempHorDataVolume(_survey_id,_slid_proceed,_elid_proceed,_lid_inv,
			_scdp_proceed,_ecdp_proceed,_cdp_inv,"v0");

	float* t0_tmp = NULL;
	if(_slice_t0_tmp) delete[] _slice_t0_tmp;
	_slice_t0_tmp = NULL;
	float* v0_tmp = NULL;
	if(_slice_v0_tmp) delete[] _slice_v0_tmp;
	_slice_v0_tmp = NULL;

	SeisInversionProgressDialog* dlg = NULL;
	dlg = new SeisInversionProgressDialog(this);
	dlg->Manage();

	int line_id_max_pos = 0,cdp_max_pos = 0;
	float t01_max = 0.0;
	float t02_max = 0.0;
	float pos1,pos2,pos3,pos4;

	for(i = 0;i<num;i++) {
		_hor_surface_data1_num1 = 0;
		_hor_surface_data1_num2 = 0;
		if(_hor_surface_data1) delete[] _hor_surface_data1;
		_hor_surface_data1 = NULL;
		_hor_surface_data2_num1 = 0;
		_hor_surface_data2_num2 = 0;
		if(_hor_surface_data2) delete[] _hor_surface_data2;
		_hor_surface_data2 = NULL;

		delta_t0 = GetDeltaT0AndDataFromDisk(i,num,&hor_sestrict_par[i],
			line_id_max_pos,cdp_max_pos,t01_max,t02_max);

		int lid_num = (_elid_proceed - _slid_proceed)/_lid_inv + 1;
		int cdp_num = (_ecdp_proceed - _scdp_proceed)/_cdp_inv + 1;
		int slice_t0_tmp_num = ( (_elid_proceed - _slid_proceed)/_lid_inv + 1)*
					( (_ecdp_proceed - _scdp_proceed)/_cdp_inv + 1);
		if(!_slice_t0_tmp) {
				_slice_t0_tmp = new float[slice_t0_tmp_num];
		}
		if(!_slice_v0_tmp) {
				_slice_v0_tmp = new float[slice_t0_tmp_num];
		}

		interpolate_curves = GetInterpolateCurvesForLayer(i,&hor_sestrict_par[i],interpolate_curves_num);

		pos1 = 100.0 * (float)(i)/(float)num;
		dlg->SetCurrentValues((int)pos1,-1,-1,-1);
		for(j = 0;j<hor_sestrict_par[i]._layer_num + 1;j++) {
			pos2 = 100.0 * (float)(j)/(float)(hor_sestrict_par[i]._layer_num + 1);
			dlg->SetCurrentValues((int)pos1,(int)pos2,-1,-1);
			t0_tmp = _slice_t0_tmp;
			v0_tmp = _slice_v0_tmp;
			float *interpolate_values = NULL;
			int *interpolate_cdps = NULL;
			int *interpolate_lids = NULL;
			int interpolate_values_num = 0;
			interpolate_values_num = GetInterpolateValue(interpolate_curves,interpolate_curves_num,i,
				&hor_sestrict_par[i],j,
				delta_t0,cdp_max_pos,t01_max,t02_max,
				interpolate_values,interpolate_cdps,interpolate_lids);
			for(k = 0;k < lid_num;k++) {
				pos3 = 100.0 * (float)(k + 1)/(float)lid_num;
				dlg->SetCurrentValues((int)pos1,(int)pos2,(int)pos3,-1);
				t0_tmp = _slice_t0_tmp + k * cdp_num;
				v0_tmp = _slice_v0_tmp + k * cdp_num;
				for(m = 0;m<cdp_num;m++){
					int line_id = 0,cdp = 0;
					float t0 = 0.0;
					float value = 0.0;
					GetLineIDAndCDPForInterpolate(line_id,cdp,t0,&hor_sestrict_par[i],
						delta_t0,line_id_max_pos,cdp_max_pos,t01_max,t02_max,j,k,m);

					if(t0 > 0) {
						value = GetWeighValue(line_id,cdp,interpolate_values,interpolate_cdps,
								interpolate_lids,interpolate_values_num);
					} else {
						value = -1.0;
					}
					if( (m > 0)&&(t0 > 0)&&(value < 0) ){
						*t0_tmp = *(t0_tmp - 1);
						*v0_tmp = *(v0_tmp - 1);
					} else {
						*t0_tmp = t0;
						*v0_tmp = value;
					}
					t0_tmp++;
					v0_tmp++;
				}//for(m = 0;m<_hor_surface_data2_num2;m++){
			}//for(k = 0;k<_hor_surface_data2_num1;k++){
			int status = hor_data_volume1->WriteHorData(_slice_t0_tmp);
			status = hor_data_volume2->WriteHorData(_slice_v0_tmp);
			delete[] interpolate_values;
			delete[] interpolate_cdps;
			delete[] interpolate_lids;
		}//for(j = 0;j<hor_sestrict_par[i]._layer_num;j++) {
		for(int mmm = 0;mmm < interpolate_curves_num;mmm++) delete interpolate_curves[mmm];
		delete[] interpolate_curves;
	}//for(i = 0;i<num;i++) {

	time(&clock);
	time_str = ctime(&clock);
	cout<<"Interpolate End Time: "<<time_str<<endl;

	if(_slice_t0_tmp) delete[] _slice_t0_tmp;
	_slice_t0_tmp = NULL;
	if(_slice_v0_tmp) delete[] _slice_v0_tmp;
	_slice_v0_tmp = NULL;

	if(_hor_surface_data1) delete[] _hor_surface_data1;
	_hor_surface_data1 = NULL;
	if(_hor_surface_data2) delete[] _hor_surface_data2;
	_hor_surface_data2 = NULL;
	int lid_num = (_elid_proceed - _slid_proceed)/_lid_inv + 1;
	int cdp_num = (_ecdp_proceed - _scdp_proceed)/_cdp_inv + 1;

	int status = 0;
	int layer_data_num = hor_data_volume1->GetNumOfHorizons();
	float** inter_data_t0 = new float*[layer_data_num];
	float** inter_data_v0 = new float*[layer_data_num];
	for(i = 0;i<layer_data_num;i++) inter_data_t0[i] = NULL;
	for(i = 0;i<layer_data_num;i++) inter_data_v0[i] = NULL;
	if(_seis_attr_pro) delete _seis_attr_pro;
	wDBSeis3DVolumeRange* seis_3d_range = NULL;
	seis_3d_range = _seis_tpt->Obtain3DTemplate()->ObtainSeis3D();
	_seis_attr_pro = new wPRSeisAttrsBase();
	_seis_attr_pro->SetOrigDataName("Amp");
	_seis_attr_pro->SetSurveyID(_survey_id);
	_seis_attr_pro->SetDomain(_data_domain);
	_seis_attr_pro->SetResultDataName("Lfv");
	_seis_attr_pro->SetProcessType(Seis3DVolumeProcessing);
	_seis_attr_pro->Set3DRange(seis_3d_range);
	_seis_attr_pro->SetCreateFlag(SeisLineVolume);
	wDBSeis3DVolumeData* seis3d_volume_data = NULL;
	seis3d_volume_data = _seis_attr_pro->Create3DVolume("Lfv");

	for(i = 0;i<lid_num;i++) {
		pos4 = 100.0 * (float)(i + 1)/(float)lid_num;
		dlg->SetCurrentValues(100,100,100,(int)pos4);
		int line_no = _slid_proceed + _lid_inv * i;
		for(j = 0;j<layer_data_num;j++) {
			if(inter_data_t0[j]) delete inter_data_t0[j];
			inter_data_t0[j] = hor_data_volume1->ReadHorData(j + 1,line_no,status);
			if(inter_data_v0[j]) delete inter_data_v0[j];
			inter_data_v0[j] = hor_data_volume2->ReadHorData(j + 1,line_no,status);
		}
		float* inter_data_t0_0 = new float[layer_data_num];
		float* inter_data_v0_0 = new float[layer_data_num];
		float* volume_data = new float[cdp_num * _t0s_num];
		for(j = 0;j<cdp_num;j++) {
			int count = 0;
			for(k = 0;k<layer_data_num;k++) {
				if( (inter_data_t0[k][j] > 0)&&(inter_data_v0[k][j] > 0) ) {
					inter_data_t0_0[count] = inter_data_t0[k][j];
					inter_data_v0_0[count] = inter_data_v0[k][j];
					if( (count > 0)&&(inter_data_t0_0[count] > inter_data_t0_0[count - 1]) ) count++;
					if(count == 0) count++;
				}
			}
			float* result;
			int result_num;
			if(count == 0) {
//				if(j > 0) memcpy(volume_data + j * _t0s_num,
//					volume_data + (j - 1) * _t0s_num,_t0s_num* sizeof(float));
				for(int ppk = 0;ppk<_t0s_num;ppk++){
					*(volume_data + j * _t0s_num + ppk) = 0.0;
				}
				continue;
			}
			result = GetLinearInterpolation(inter_data_t0_0,inter_data_v0_0,count,
				(float)_st0,(float)_et0,(float)_dt0,result_num);
			memcpy(volume_data + j * _t0s_num,result,result_num * sizeof(float));
			delete[] result;
		}
		status = seis3d_volume_data->WriteInLine(line_no,volume_data);
		delete[] volume_data;
		delete[] inter_data_t0_0;
		delete[] inter_data_v0_0;
	}
	for(i = 0;i<layer_data_num;i++) delete[] inter_data_t0[i];
	delete[] inter_data_t0;
	for(i = 0;i<layer_data_num;i++) delete[] inter_data_v0[i];
	delete[] inter_data_v0;
	delete seis3d_volume_data;


	delete hor_data_volume1;
	delete hor_data_volume2;
	delete dlg;
	time(&clock);
	time_str = ctime(&clock);
	cout<<"End Time: "<<time_str<<endl;
}

/////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////
int SeisInversion::Get2DHorizonSurfaceData()
{
	int i = 0;
	wURHorizonRestrictPars* hor_sestrict_par = NULL;
	int num = 0;
	hor_sestrict_par = _seis_tpt->GetHorizonPars(num);
	_seis_tpt->Obtain2DTemplate()->ObtainSeis2D()->GetCDPRange(_scdp,_ecdp,_cdp_inv);
	_seis_tpt->Obtain2DTemplate()->ObtainSeis2D()->GetT0DepthRange(_st0,_et0,_dt0);
	_2d_cdp_num = (_ecdp - _scdp)/_cdp_inv + 1;
	wURSeisHorizonCurve* temp_hor_curves = NULL;
	int temp_hor_curves_num = 0;
	if(num > 0) {
		wDBSeisHorizonInSurvey* h = NULL;
		h = new wDBSeisHorizonInSurvey();
		for(i = 0;i<num - 1;i++) {
			temp_hor_curves = new wURSeisHorizonCurve();
			int status = h->OpenHorInSurvey(hor_sestrict_par[i]._bottom_hor_id,
				hor_sestrict_par[i]._survey_id);
			if(status != OK_STATUS){
				delete temp_hor_curves;
				continue;
			}
			
			char *hs_id = NULL;
			hs_id = h->GetHorSurveyID();
			if(!hs_id) {
				delete temp_hor_curves;
				continue;
			}

			temp_hor_curves->SetHorSurveyID(hs_id);
			status = temp_hor_curves->SetPars(_otype,_oid,_scdp,_ecdp,_cdp_inv);
			delete hs_id;
			if(status != OK_STATUS){
				delete temp_hor_curves;
				continue;
			}
			float* data = NULL;
			int data_num = 0;
			data = temp_hor_curves->ObtainValues(data_num);
			if( (data_num == 0)||(data == NULL) ){
				delete temp_hor_curves;
				continue;
			}
			temp_hor_curves_num++;
			delete temp_hor_curves;
		}
		delete h;
	}
	if(temp_hor_curves_num != num - 1) return 100;

	if(num > 0) {
		_2d_hor_curves = new wURSeisHorizonCurve[num - 1];
		for(i = 0;i<num - 1;i++) {
			wDBSeisHorizonInSurvey* h = NULL;
			h = new wDBSeisHorizonInSurvey();
			h->OpenHorInSurvey(hor_sestrict_par[i]._bottom_hor_id,
				hor_sestrict_par[i]._survey_id);
			char *hs_id = NULL;
			hs_id = h->GetHorSurveyID();
			_2d_hor_curves[i].SetHorSurveyID(hs_id);
			int status = _2d_hor_curves[i].SetPars(_otype,_oid,_scdp,_ecdp,_cdp_inv);

			delete hs_id;
			delete h;
			_2d_buffers_num = _2d_buffers_num + hor_sestrict_par[i]._layer_num + 1;
		}

		_2d_buffers_num = _2d_buffers_num + hor_sestrict_par[num - 1]._layer_num + 1;
		_2d_buffers_t0 = new float*[_2d_buffers_num];
		_2d_buffers_v0 = new float*[_2d_buffers_num];
		for(i = 0;i<_2d_buffers_num;i++) {
			_2d_buffers_t0[i] = new float[_2d_cdp_num];
			_2d_buffers_v0[i] = new float[_2d_cdp_num];
		}
	}

	//Get Curves:
	wDBScatterDataInSurvey *scatter_data = NULL;
	scatter_data = new wDBScatterDataInSurvey();
	wURScatterDataInSection* section_data = NULL;
	section_data = _seis_tpt->ObtainScatterSection();
	char sc_id[16];
	section_data->GetScatterID(sc_id);
	scatter_data->SetScatterID(sc_id);
	
	scatter_data->OpenScatterData();
	wDBScatterCurveInSurvey** curves = NULL;
	curves = scatter_data->ListScatterCurves(_curves_in_survey_num);
	cout<<"_curves_in_survey_num = "<<_curves_in_survey_num<<endl;
	_curves_in_survey = new wDBInversionCurveInSurvey*[_curves_in_survey_num];
	int values_num = 0;
	_seis_tpt->Obtain2DTemplate()->ObtainSeis2D()->GetT0DepthRange(_st0,_et0,_dt0);