kev/Drawer/GVision/GmtSurfaceGrid/GmtSurfaceInterp.cpp

#include "GmtSurfaceInterp.h"
#include <afx.h>
#include <afxext.h>
#include <gmt/gmt.h>
#include <netcdf.h>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <vector>
#include <string>
#include <fstream>
#include <algorithm>
#include <string>

#include "DrawModel\\BaseLib.h"
#include "DrawOperator\\DrawLib.h"
#include <sstream>

template <typename T>
void qDeleteAll(std::vector<T*>& vec)
{
	for (auto* p : vec) {
		delete p;
	}

	vec.clear();
}

// »ñÈ¡±ß½çÏß/±ß¿ò
static std::vector<CCurveEx*> get_border(CString borderFile)
{
	if (!CFindFileEx::IsFileExists(borderFile)) {
		return {};
	}
	std::shared_ptr<CXy> pXy = std::make_shared<CXy>();

	borderFile = borderFile.Trim();
	if (borderFile.GetLength() <= 0 || borderFile == "NULL") {
		return {};
	}

	if (!pXy->ReadWithExtension(borderFile)) {
		return {};
	}

	std::vector<CCurveEx*> result;

	CPtrList* pl = pXy->GetValueList();
	for (POSITION pos = pl->GetHeadPosition(); pos != nullptr; pl->GetNext(pos)) {
		COne* pOne = (COne*)pl->GetAt(pos);
		if (pOne->GetType() == DOUBLEFOX_CURVE) {
			CCurveEx* pCurve = (CCurveEx*)(pOne->GetValue());
			CCurveEx* pCurveNew = new CCurveEx();
			*pCurveNew = *pCurve;
			result.push_back(pCurveNew);
		}
	}

	return result;
}

static void CreateBorder(std::shared_ptr<CXy> pXy, std::vector<CCurveEx*>& ccurves)
{
	CLayer* pLayerBorder = pXy->FindAddLayer("±ß½ç11");
	for (auto* pCurve : ccurves) {
		POSITION posNew = pXy->AddElement(pCurve, DOUBLEFOX_CURVE);
		pXy->SetElementLayer(posNew, pLayerBorder);
	}
}

static void AddContourCurve(std::shared_ptr<CXy> pXy, vector<CCurve*>* curves, vector<CString*>* layer)
{
	for (size_t i = 0; i < curves->size(); i++)
	{
		CCurve* pCurve = curves->at(i);
		if (pCurve == NULL) continue;

		CString* pName = layer->at(i);

		CCurveEx* ce = new CCurveEx(pCurve->num);
		for (int i = 0; i < ce->num; i++)
		{
			ce->x[i] = pCurve->x[i];
			ce->y[i] = pCurve->y[i];
			ce->z[i] = pCurve->z[i];
		}
		ce->nPoint = pCurve->nPoint;
		ce->GetLocation();
		POSITION pos = NULL;
		if (pCurve->name)
			ce->SetName(pCurve->name);

		pos = pXy->AddElement(ce, DOUBLEFOX_CURVE);
		pXy->SetElementLayer(pos, *pName);
	}
}

/**
 * Éú³ÉµÈÖµÏß
 *
 * \\param pXy
 * \\param pMeshNew
 * \\param pDfg
 * \\param contourStep
 * \\param contourMarkStep
 * \\return
 */
static bool CreateContourLines(std::shared_ptr<CXy> pXy, CMesh* pMeshNew, double contourStep, int contourMarkStep)
{
	// Éú³ÉµÈÖµÏß
	if (abs(contourStep) < 1E-5) 
	{
		return false;
	}
	CString strLayerMark  = _T("Layer:\\µÈÖµÏß\\±ê×¢");
	CString strLayerOther = _T("Layer:\\µÈÖµÏß\\ÎÞ±ê×¢");

	CString curLayerName = pXy->GetCurrentLayer()->GetPathName();
	CLayer* pMarkLayer   = pXy->FindLayer(strLayerMark);
	CLayer* pOtherLayer  = pXy->FindLayer(strLayerOther);

	if (pMarkLayer == nullptr)
	{
		pMarkLayer = pXy->FindAddLayer(strLayerMark);
		pMarkLayer->HowToViewCurve = new CHowToViewCurve();
		pMarkLayer->HowToViewCurve->EnableDrawSourceCurve(FALSE);
		CCurveInName* pInName = new CCurveInName();
		CRect8 rect = pMeshNew->GetRect();
		pInName->text_h = rect.Width() / 300;
		pInName->m_size.cx = pInName->text_h * 0.06;
		pInName->color = RGB(0, 0, 0);
		pMarkLayer->HowToViewCurve->Add(pInName);
	}

	if (pOtherLayer == nullptr)
	{
		pOtherLayer = pXy->FindAddLayer(strLayerOther);
		pOtherLayer->HowToViewCurve = new CHowToViewCurve();
		pOtherLayer->HowToViewCurve->EnableDrawSourceCurve(FALSE);
		CCurveProperties* pview = new CCurveProperties();
		pview->color = RGB(0, 0, 0);
		pOtherLayer->HowToViewCurve->Add(pview);
	}

	std::vector<CCurve*> pCurves;
	std::vector<CString*> pLayers;

	pMeshNew->ContourCreate(&pCurves, &pLayers, contourStep, contourMarkStep,
		strLayerMark, strLayerOther, pMeshNew->GetDfg()->range[0], pMeshNew->GetDfg()->range[1]);

	AddContourCurve(pXy, &pCurves, &pLayers);
	qDeleteAll(pCurves);
	qDeleteAll(pLayers);
	return true;
}

/**
 * ÅÐ¶Ï dx dY ×ø±êÊÇ·ñÔÚÆäÖÐÈÎÒâÒ»¸öÇøÓò
 *
 * \\param vec ÇøÓò¼¯ºÏ£¬ËäÈ»ÊÇÕÛÏßÀ࣬µ«ÕâÀïÊÇÒ»¸ö¸ö±ÕºÏµÄÇøÓò
 * \\param dX
 * \\param dY
 * \\return
 */
static bool any_contains(std::vector<CCurveEx*>& vec, double dX, double dY)
{
	return std::any_of(vec.begin(), vec.end(), [dX, dY](CCurveEx* pBorder) {
		return (pBorder->IsInside(dX, dY));
		});
}

/**
 * ½«Íø¸ñÐÅÏ¢±£´æµ½ÎļþÖУ¬¸ù¾ÝÎļþºó׺Ãû×Ô¶¯Éú³É¶ÔÓ¦¸ñʽ
 *
 * \\param pXy
 * \\param gridFile 
 * \\param outputFile  Òª±£´æµÄÎļþÃû
 * \\return
 */
static bool SaveFile(CString gridFile, CString outputFile, double contourStep, int contourMarkStep, CString borderFile)
{
	std::shared_ptr<CXy> pXy = std::make_shared<CXy>();
	pXy->FromGridAuto(gridFile);

	CLayer *pLayer = pXy->FindAddLayer("±³¾°");
	CPositionList pointList;
	int index = pXy->GetElement(DOUBLEFOX_MESH, pointList);
	if (pointList.GetSize() == 0)
	{
		return false;
	}

	COne* pNewOne = pXy->GetAt(pointList.GetHead());
	pNewOne->SetLayer(pLayer);
	CMesh* pMeshNew = (CMesh*)pNewOne->value;
	CreateContourLines(pXy, pMeshNew, contourStep, contourMarkStep);
	
	pXy->SaveAsWithExtension(outputFile);
	return true;
}

bool WriteDSAA(const char* netcdf_file, const char* output_grd)
{
	int ncid, varid, xvar, yvar;
	if (nc_open(netcdf_file, NC_NOWRITE, &ncid) != NC_NOERR)
	{
		fprintf(stderr, "Failed to open NetCDF file: %s\n", netcdf_file);
		return false;
	}

	nc_inq_varid(ncid, "z", &varid);
	nc_inq_varid(ncid, "x", &xvar);
	nc_inq_varid(ncid, "y", &yvar);

	size_t x_len, y_len;
	nc_inq_dimlen(ncid, 0, &x_len);
	nc_inq_dimlen(ncid, 1, &y_len);

	std::vector<float> z(x_len * y_len);
	std::vector<double> xs(x_len);
	std::vector<double> ys(y_len);

	nc_get_var_float(ncid, varid, z.data());
	nc_get_var_double(ncid, xvar, xs.data());
	nc_get_var_double(ncid, yvar, ys.data());

	nc_close(ncid);

	float zmin = 1e30f, zmax = -1e30f;
	for (float val : z)
	{
		if (!std::isnan(val)) 
		{
			zmin = min(zmin, val);
			zmax = max(zmax, val);
		}
	}

	std::ofstream out(output_grd);
	out << "DSAA\n";
	out << x_len << " " << y_len << "\n";
	out << xs.front() << " " << xs.back() << "\n";
	out << ys.front() << " " << ys.back() << "\n";
	out << zmin << " " << zmax << "\n";

	for (size_t j = 0; j < y_len; ++j)
	{
		for (size_t i = 0; i < x_len; ++i) 
		{
			size_t idx = j * x_len + i;
			float v = z[idx];
			if (std::isnan(v)) out << "1.000000E+30 ";
			else out << v << " ";
			if ((i + 1) % 10 == 0) out << "\n";
		}
		out << "\n";
	}
	out.close();
	return true;
}

//bool ConvertFaultFile(const char* input_file, const char* output_file) 
//{
//	std::ifstream infile(input_file);
//	if (!infile.is_open())
//	{
//		fprintf(stderr, "Failed to open fault file: %s\n", input_file);
//		return false;
//	}
//
//	std::ofstream outfile(output_file);
//	if (!outfile.is_open()) 
//	{
//		fprintf(stderr, "Failed to create converted fault file: %s\n", output_file);
//		return false;
//	}
//
//	double x, y;
//	int id, last_id = -1;
//	while (infile >> x >> y >> id)
//	{
//		if (id != last_id)
//		{
//			outfile << "> -Z\n";
//			last_id = id;
//		}
//		outfile << x << " " << y << "\n";
//	}
//
//	infile.close();
//	outfile.close();
//	return true;
//}

/**
 * ½«dfd¸ñʽÎļþת»»³ÉgmtÄÜʶ±ðµÄÊý¾Ý
 *
 * \\param input_file  .dfd
 * \\param output_file .dat
 * \\param layer_name  Layer M ¶Ï²ã / Layer M ±ß½ç
 * \\return
 */
bool ConvertDfdFile(const char* input_file, const char* output_file, const char* layer_name)
{
	std::ofstream outfile(output_file);
	if (!outfile.is_open()) {
		fprintf(stderr, "Failed to create converted file: %s\n", output_file);
		return false;
	}
	CString cstrBorderFile(input_file);
	std::vector<CCurveEx*> borderList = get_border(cstrBorderFile);
	for (const auto& curve : borderList) {
		if (!curve) continue;

		// ÿ¸öÇúÏßдһ¸ö Pline
		outfile << "> -Z\n";
		for (size_t i = 0; i < curve->num; ++i) {
			outfile << curve->x[i] << " " << curve->y[i] << "\n";
		}
	}

	outfile.close();
	return true;
}


bool ExtractBounds(const char* xyz_input, double bounds[4]) 
{
	char cmd[256];
	snprintf(cmd, sizeof(cmd), "SplineInterpolation info %s -C", xyz_input);

	FILE* pipe = _popen(cmd, "r");
	if (!pipe)
	{
		fprintf(stderr, "Failed to run SplineInterpolation info.\n");
		return false;
	}

	if (fscanf_s(pipe, "%lf %lf %lf %lf", &bounds[0], &bounds[1], &bounds[2], &bounds[3]) != 4) 
	{
		fprintf(stderr, "Failed to parse SplineInterpolation info output.\n");
		_pclose(pipe);
		return false;
	}

	_pclose(pipe);
	return true;
}

bool GmtSurfaceInterp(const char* xyz_input, const char* grid_output, double step, double tension,
	const char* fault_file, const char* boundary_file, const char* outputfile, double contourStep, int contourMarkStep,
	double XMin, double YMin, double XMax, double YMax)
{
	if (step <= 0 || tension < 0)
		return false;

	double clip = 0.001;

	void* session = GMT_Create_Session("surface_interp", GMT_SESSION_EXTERNAL, 0, NULL);
	if (!session) {
		fprintf(stderr, "Failed to create GMT session.\n");
		return false;
	}

	bool hasFault = fault_file && strlen(fault_file) > 0;
	bool hasBoundary = boundary_file && strlen(boundary_file) > 0;

	// ת»»¶Ï²ã dfd
	if (hasFault) {
		if (!ConvertDfdFile(fault_file, "converted_fault.dat", "Layer M ¶Ï²ã")) {
			fprintf(stderr, "Failed to convert fault file.\n");
			GMT_Destroy_Session(session);
			return false;
		}
	}

	// ת»»±ß½ç dfd
	if (hasBoundary) {
		if (!ConvertDfdFile(boundary_file, "converted_boundary.dat", "Layer M ±ß½ç")) {
			fprintf(stderr, "Failed to convert boundary file.\n");
			GMT_Destroy_Session(session);
			return false;
		}
	}

	GMT_DATASET* data = (GMT_DATASET*)GMT_Read_Data(session, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_PLP,
		GMT_READ_NORMAL, NULL, xyz_input, NULL);
	if (!data) {
		fprintf(stderr, "Failed to read input file: %s\n", xyz_input);
		GMT_Destroy_Session(session);
		return false;
	}

	char input_vf[GMT_VF_LEN] = { 0 };
	if (GMT_Open_VirtualFile(session, GMT_IS_DATASET, GMT_IS_PLP, GMT_IN, data, input_vf) != GMT_NOERROR) {
		fprintf(stderr, "Failed to register virtual file.\n");
		GMT_Destroy_Session(session);
		return false;
	}

	double bounds[4];
	bounds[0] = XMin;
	bounds[1] = XMax;
	bounds[2] = YMin;
	bounds[3] = YMax;

	if (XMin == 0 && XMax == 0)
	{
		//Èç¹ûĬÈϱ߽çûÓУ¬ÔòʹÓÃgmt»ñÈ¡
		if (!ExtractBounds(xyz_input, bounds)) {
			fprintf(stderr, "Failed to extract region from dataset.\n");
			GMT_Destroy_Session(session);
			return false;
		}
	}

	int nx = static_cast<int>((bounds[1] - bounds[0]) / step + 0.5);
	int ny = static_cast<int>((bounds[3] - bounds[2]) / step + 0.5);
	double adj_x_max = bounds[0] + nx * step;
	double adj_y_max = bounds[2] + ny * step;

	char tmp_nc[] = "tmp_output.nc";

	char args[512];
	snprintf(args, sizeof(args),
		"-R%.8f/%.8f/%.8f/%.8f -I%.8f/%.8f -T%.3f -C%.4f -G%s -Vn %s",
		bounds[0], adj_x_max, bounds[2], adj_y_max,
		step, step, tension, clip, tmp_nc, input_vf);

	if (GMT_Call_Module(session, "surface", GMT_MODULE_CMD, args) != GMT_NOERROR) {
		fprintf(stderr, "Surface interpolation failed.\n");
		GMT_Destroy_Session(session);
		return false;
	}

	if (hasFault) {
		char mask_cmd[512];
		snprintf(mask_cmd, sizeof(mask_cmd),
			"converted_fault.dat -R%.8f/%.8f/%.8f/%.8f -I%.8f/%.8f -N1/1/NaN -Gfault_mask.nc -Vn",
			bounds[0], adj_x_max, bounds[2], adj_y_max, step, step);

		if (GMT_Call_Module(session, "grdmask", GMT_MODULE_CMD, mask_cmd) != GMT_NOERROR) {
			fprintf(stderr, "Failed to generate mask from fault file.\n");
			GMT_Destroy_Session(session);
			return false;
		}
	}

	if (hasBoundary) {
		char boundary_mask_cmd[512];
		snprintf(boundary_mask_cmd, sizeof(boundary_mask_cmd),
			"converted_boundary.dat -R%.8f/%.8f/%.8f/%.8f -I%.8f/%.8f -NNaN/1 -Gboundary_mask.nc -Vn",
			bounds[0], adj_x_max, bounds[2], adj_y_max, step, step);

		if (GMT_Call_Module(session, "grdmask", GMT_MODULE_CMD, boundary_mask_cmd) != GMT_NOERROR) {
			fprintf(stderr, "Failed to generate mask from boundary file.\n");
			GMT_Destroy_Session(session);
			return false;
		}
	}

	if (hasFault) {
		char math_cmd[256];
		snprintf(math_cmd, sizeof(math_cmd), "fault_mask.nc %s MUL = temp.nc", tmp_nc);

		if (GMT_Call_Module(session, "grdmath", GMT_MODULE_CMD, math_cmd) != GMT_NOERROR) {
			fprintf(stderr, "Failed to apply fault mask.\n");
			GMT_Destroy_Session(session);
			return false;
		}
		strcpy_s(tmp_nc, "temp.nc");
	}

	if (hasBoundary) {
		char math_cmd2[256];
		snprintf(math_cmd2, sizeof(math_cmd2), "boundary_mask.nc %s MUL = clipped.nc", tmp_nc);

		if (GMT_Call_Module(session, "grdmath", GMT_MODULE_CMD, math_cmd2) != GMT_NOERROR) {
			fprintf(stderr, "Failed to apply boundary mask.\n");
			GMT_Destroy_Session(session);
			return false;
		}
		strcpy_s(tmp_nc, "clipped.nc");
	}

	GMT_Destroy_Session(session);

	if (!WriteDSAA(tmp_nc, grid_output)) {
		fprintf(stderr, "Failed to write .grd output file.\n");
		return false;
	}

	if (hasFault) {
		std::remove("converted_fault.dat");
		std::remove("fault_mask.nc");
		std::remove("temp.nc");
	}

	if (hasBoundary) {
		std::remove("converted_boundary.dat");
		std::remove("boundary_mask.nc");
		std::remove("clipped.nc");
	}

	std::remove("tmp_output.nc");

	SaveFile(grid_output, outputfile, contourStep, contourMarkStep, boundary_file);
	std::remove(grid_output);

	return true;
}