kev/Drawer/GVision/FaciesOutline/FaciesOutlineDetection/GSurface.cpp

//#include "stdafx.h"
#include "GSurface.h"
//#include "BaseElements.h"
#include <vector>
#include <numeric>
#include <algorithm>
using namespace std;

GSurface::GSurface()
	:m_u(nullptr)
{
	memset(m_num, 0, sizeof(int) * 2);
	memset(m_delt, 0, sizeof(double) * 2);
	memset(m_P0, 0, sizeof(double) * 2);
	memset(m_range, 0, sizeof(double) * 2);
}


GSurface::~GSurface()
{
	Clear();
}

bool GSurface::ReadDfg(FILE* fr)
{


	
	int numx = 0, numy= 0;
	double xmin, ymin, dx, dy, zmin, zmax;
	

	fseek(fr, 3, SEEK_CUR);   fread((char *)(&numx), 4, 1, fr);
	fseek(fr,3, SEEK_CUR);   fread((char *)(&numy), 4, 1, fr);
	fseek(fr,6, SEEK_CUR);   fread((char *)(&xmin), 8, 1, fr);
	fseek(fr,6, SEEK_CUR);   fread((char *)(&ymin), 8, 1, fr);
	fseek(fr,6, SEEK_CUR);   fread((char *)(&dx), 8, 1, fr);
	fseek(fr,6, SEEK_CUR);   fread((char *)(&dy), 8, 1, fr);
	fseek(fr,6, SEEK_CUR);   fread((char *)(&zmin), 8, 1, fr);
	fseek(fr,6, SEEK_CUR);   fread((char *)(&zmax), 8, 1, fr);

	char format = 8;
	fread(&format, 1, 1, fr);
	//fr.Read(&format, 1);

	Create(numx, numy, xmin, ymin, dx, dy);
	m_range[0] = zmin;
	m_range[1] = zmax;

	int i, nn;
	nn = numx*numy;


	switch (format)
	{
	case 2:		//short
	{

		short* v = new short[numx];
		for (int j = 0; j < numy; ++j)
		{
		
			fread(v, 2 * numx,1,fr);
			for (i = 0; i < numx; ++i)
			{
				m_u[j*numx + i] = v[i];
			}
		}
		delete v;
	}
	
	break;
	case 4:		//float
	{

		float* v = new float[numx];
		for (int j = 0; j < numy; ++j)
		{
			fread(v, 4 * numx, 1, fr);
		
			for (i = 0; i < numx; ++i)
			{
				m_u[j*numx + i] = v[i];
			}
		}
		delete v;
	}

	break;
	default:	//double
		fread(m_u, sizeof(double)*nn,1,fr);	//2531272

	}





	return true;

}
bool GSurface::ReadDfg(const char* filename)
{
	FILE* fr = fopen(filename, "rb");
	if (nullptr == fr)
		return false;

	bool bflag = ReadDfg(fr);

	fclose(fr);

	return bflag;


}

bool GSurface::WriteDfg(FILE* fw)
{


	double flen;
	//unsigned long ll;
	//unsigned short ss;

	////////////////////////////////////////////////////////////////////
	int numx, numy;
	double xmin, ymin, xmax, ymax, zmin, zmax, dx, dy;
	xmin = m_P0[0];
	ymin = m_P0[1];
	xmax = m_P0[0] + m_num[0] * m_delt[0];
	ymax = m_P0[1] + m_num[1] * m_delt[1];
	dx = m_delt[0];
	dy = m_delt[1];
	numx = m_num[0];
	numy = m_num[1];
	zmin = 1;



	fwrite("m=\n", 1, 3, fw);	fwrite((char *)(&numx), 4, 1, fw);
	fwrite("n=\n", 1, 3, fw);	fwrite((char *)(&numy), 4, 1, fw);
	fwrite("xmin=\n", 1, 6, fw);	fwrite((char *)(&xmin), 8, 1, fw);
	fwrite("ymin=\n", 1, 6, fw);	fwrite((char *)(&ymin), 8, 1, fw);
	fwrite("  dx=\n", 1, 6, fw);	fwrite((char *)(&dx), 8, 1, fw);
	fwrite("  dy=\n", 1, 6, fw);	fwrite((char *)(&dy), 8, 1, fw);
	int i;


	zmin = m_range[0];
	zmax = m_range[1];
	fwrite("zmin=\n", 1, 6, fw); fwrite((char *)(&zmin), 8, 1, fw);
	fwrite("zmax=\n", 1, 6, fw); fwrite((char *)(&zmax), 8, 1, fw);

	char nDataType = 8;
	fwrite(&nDataType, 1, 1, fw);

	//////////////////////////////////////////
	flen = numx*numy;
	switch (nDataType)
	{
	case 2:		//short
	{
		short v;
		for (int j = 0; j < numy; ++j)
		{
			
			for (i = 0; i < numx; ++i)
			{
				v = (short)m_u[j*numx + i];
				fwrite((char *)(&v), 2, 1, fw);
			}
		}
	}

	break;
	case 4:		//float
	{
		float v;
		for (int j = 0; j < numy; ++j)
		{
			
			for (i = 0; i < numx; ++i)
			{
				v = (float)m_u[j*numx + i];
				fwrite((char *)(&v), 4, 1, fw);
			}
		}
	}

	break;
	case 8:		//double
	default:
		for (int j = 0; j < numy; ++j)
		{
			
			fwrite((char *)(m_u + j*numx), 8, numx, fw);
	
		}

	}

	return 1;
}
bool GSurface::WriteDfg(const char* filename)
{
	FILE* fw = fopen(filename, "wb");
	if (nullptr == fw)
		return false;

	bool bflag = WriteDfg(fw);

	fclose(fw);

	return bflag;


}


int GSurface::Create(int numx, int numy, double x0, double y0, double dx, double dy)
{
	Clear();
	m_num[0] = numx;
	m_num[1] = numy;
	m_P0[0] = x0;
	m_P0[1] = y0;
	m_delt[0] = dx;
	m_delt[1] = dy;

	int N = numx*numy;
	m_u = new double[N];
	for (int i = 0; i < N; i++)
		m_u[i] = 1e31;

	return 1;
}
int GSurface::Create(int numx, int numy, double x0, double y0, double dx, double dy,double* values)
{
	Clear();
	m_num[0] = numx;
	m_num[1] = numy;
	m_P0[0] = x0;
	m_P0[1] = y0;
	m_delt[0] = dx;
	m_delt[1] = dy;

	int N = numx * numy;
	m_u = new double[N];
	for (int i = 0; i < N; i++)
		m_u[i] = values[i];

	this->CalcRange();
	return 1;
}
void GSurface::Clear(void) //<2F><><EFBFBD><EFBFBD>
{
	memset(m_num, 0, sizeof(int) * 2);
	memset(m_delt, 0, sizeof(double) * 2);
	memset(m_P0, 0, sizeof(double) * 2);
	memset(m_range, 0, sizeof(double) * 2);
	if (m_u)
	{
		delete[]m_u;
		m_u = nullptr;
	}
}

double GSurface::X(int i)
{
	return m_P0[0] + m_delt[0] * i;
}
double GSurface::Y(int j)
{
	return m_P0[1] + m_delt[1] * j;
}

int GSurface::XNum(void)
{
	return m_num[0];
}
int GSurface::YNum(void)
{
	return m_num[1];
}

double GSurface::Z(int i, int j)
{
	return m_u[j*XNum() + i];
}

double GSurface::DeltX()
{
	return m_delt[0];
}
double GSurface::DeltY()
{
    return m_delt[1];
}

double *GSurface::GetRange()
{
    return m_range;
}

void GSurface::CalcRange(void)
{

	int N = XNum()*YNum();
	double zmin, zmax;
	zmin = zmax = m_u[0];
	double u = 0;
	for (int i = 1; i < N; i++)
	{
		u = m_u[i];
		if (u > 1e20)
			continue;
        zmin =  std::min(zmin, u);
		zmax = max(zmax, u);
	}

	m_range[0] = zmin;
	m_range[1] = zmax;
}


bool GSurface::IsInRange(double z0)
{
	if (z0 > m_range[1] + 1e-6)
		return false;
	if (z0 < m_range[0] - 1e-6)
		return false;
	return true;
}

////ֻȡ<D6BB>ϲ<EFBFBD><CFB2>ⲿ<EFBFBD>Ľ<EFBFBD><C4BD><EFBFBD><EFBFBD><EFBFBD>ֵ
//double GSurface::Value(double x0, double y0, GPline* png)
//{

//	int i, j, k;
//	i = (int)((x0 - m_P0[0]) / m_delt[0]);
//	if (i<0 || i>m_num[0] - 1)
//		return 1e31;

//	j = (int)((y0 - m_P0[1]) / m_delt[1]);
//	if (j<0 || j>m_num[1] - 1)
//		return 1e31;

//	if (i == XNum() - 1 || j == YNum() - 1)
//		return 1e31;


//	k = i + j*m_num[0];

//	//ѡ<><D1A1><EFBFBD>ϲ<EFBFBD><CFB2>ⲿ<EFBFBD>Ľ<EFBFBD><C4BD><EFBFBD>
//	GPoint3D pt;
//	vector<GPoint3D> dstpts;

//	const int NX[4] = { i,i+1,i+1,i };
//	const int NY[4] = {j,j,j+1,j+1};
//	for (int kk = 0; kk < 4; kk++)
//	{
//		pt.x0 = this->X(NX[kk]);
//		pt.y0 = this->Y(NY[kk]);
//		if (0 != png->IsInside(pt.x0, pt.y0))
//			continue;
//		pt.z0 = this->Z(NX[kk], NY[kk]);
//		dstpts.push_back(pt);
//	}


//	if (dstpts.empty() )
//		return 1e31;

//	if (dstpts.size() < 4)
//	{
//		vector<double> dists;
//		double d;
//		pt = GPoint3D(x0, y0, 0);
//		for (int kk = 0; kk < dstpts.size(); kk++)
//		{
//			d = dstpts[kk].Distance2D(pt);
//			if (d < 1e-6)
//				return dstpts[kk].z0;
//			dists.push_back(1.0/d);
//		}
//		double vv = 0.0;
//		for (int kk = 0; kk < dstpts.size(); kk++)
//			vv += dists[kk] * dstpts[kk].z0;
		

//		vv /= double(std::accumulate(dists.begin(), dists.end(), 0.0));

//		return vv;
//	}



//	if (m_num[1] - 1 == j)
//	{
//		if (!IsInRange(m_u[k]) || !IsInRange(m_u[k + 1]))
//			return 1e101;
//	}
//	else if (
//		!IsInRange(m_u[k]) ||
//		!IsInRange(m_u[k + 1]) ||
//		!IsInRange(m_u[k + m_num[0]]) ||
//		!IsInRange(m_u[k + m_num[0] + 1]))
//		return 1e31;

//	double kx = (x0 - X(i)) / m_delt[0];
//	double ky = (y0 - Y(j)) / m_delt[1];

//	if (fabs(kx) < 1e-6 && fabs(ky) < 1e-6)
//		return Z(i, j);

//	double z1 = Z(i, j) + kx*(Z(i + 1, j) - Z(i, j));
//	double z2 = Z(i, j + 1) + kx*(Z(i + 1, j + 1) - Z(i, j + 1));

//	double z0 = z1 + ky*(z2 - z1);

//	return z0;


//}



double GSurface::Z(double x0, double y0)
{
	int i, j, k;
	i = (int)((x0 - m_P0[0]) / m_delt[0]); 
	if (i<0 || i>m_num[0] - 1)
		return 1e31;
	
	j = (int)((y0 - m_P0[1]) / m_delt[1]);
	if (j<0 || j>m_num[1] - 1)
		return 1e31;

	if (i == XNum() - 1 || j == YNum() - 1)
		return 1e31;

	
	k = i + j*m_num[0];

	
	if (m_num[1] - 1 == j) 
	{
		if (!IsInRange(m_u[k]) || !IsInRange(m_u[k + 1]))
			return 1e101;
	}
	else if (
		!IsInRange(m_u[k]) ||
		!IsInRange(m_u[k + 1]) ||
		!IsInRange(m_u[k + m_num[0]]) ||
		!IsInRange(m_u[k + m_num[0] + 1]))
		return 1e31;

	double kx = (x0 - X(i)) / m_delt[0];
	double ky = (y0 - Y(j)) / m_delt[1];

	if (fabs(kx) < 1e-6 && fabs(ky) < 1e-6)
		return Z(i, j);

	double z1 = Z(i, j) + kx*(Z(i + 1, j) - Z(i, j));
	double z2 = Z(i, j + 1) + kx*(Z(i + 1, j + 1) - Z(i, j + 1));

	double z0 = z1 + ky*(z2 - z1);

	return z0;



}

void GSurface::setZ(int i, int j, double z)
{

         m_u[j*XNum() + i] = z;

}