( algorithm ) Tension 3DMESH
Here is a function I wrote same time ago that calculates the minimal mesh inside of four arbitrary curves in space.
You have to select the left (red), right (green), front (blue) and back (yellow) curve. The code will create a NxM mesh.
To use it for your case you just have to generate the bounding curves from your data i.e. as AcDb3dPolyline or AcDbSpline.
void cmdMinSurf()
{
ads_name ename;
ads_point pickpoint;
int stat;
LPCWSTR name[2][2] = { { L"left", L"right"}, { L"front", L"back"} };
TCHAR msg[256];
AcDbObjectId idCurve[2][2];
Acad::ErrorStatus es;
int i, j, count;
static int N = 40, M = 50; // Grid size. Made static for quick experiments.
// I'm using AcGeVector3d instead of AcGePoint3d because so I can simply add and multiply.
// p[0][M-1] - p[1][M-1] - ... p[N-2][M-1] - p[N-1][M-1]
// +-----------------------------+
// p[0][M-2] | p[1][M-2] - ... p[N-2][M-2] | p[N-1][M-2]
// ... . ... ... .
// ... . ... ... .
// p[0][1] | p[1][1] - ... p[N-2][1] | p[N-1][1]
// +-----------------------------+
// p[0][0] - p[1][0] - ... p[N-2][0] - p[N-1][0]
//
// p[0][j] and p[N-1][j] are the fixed points on the given curves in X-directions
// p[i][0] and p[i][M-1] are the fixed points on the given curves in Y-directions
// p[1..N-2][1..M-2] are the points inside of the grid
std::vector< std::vector<AcGeVector3d > > p; // works like AcGeVector3d p[N][M] - but dynamic size
p.resize(N);
for (auto& arr : p)
arr.resize(M);
AcGePoint3d pt;
AcGePoint3d startpoints[2][2], endpoints[2][2];
for (i=0; i<2; ++i) // 0: left/right 1:front/back
{
for (j = 0; j < 2; ++j) // 0: left/front 1:right/back
{
swprintf_s(msg, L"\nSelect %s border curve: ", name[i][j]);
stat = acedEntSel(msg, ename, pickpoint);
if (stat != RTNORM)
return;
es = acdbGetObjectId(idCurve[i][j], ename);
if (!es)
{
AcDbCurve *curve;
if ((es = acdbOpenObject(curve, idCurve[i][j], AcDb::kForRead)) == Acad::eOk)
{
curve->getStartPoint(startpoints[i][j]);
curve->getEndPoint(endpoints[i][j]);
curve->close();
}
}
}
}
std::vector<double> distances;
distances.resize(4);
bool reverse[2][2];
//startpoints[2][2], endpoints[2][2];
// The "normal" case is: curves run from left to right and from front to back
// we must make sure that the l/r curves and the f/b curves don't run in opposite directions
distances[0] = startpoints[0][0].distanceTo(startpoints[1][0]); // should be 0 (or minimal)
distances[1] = startpoints[0][0].distanceTo(endpoints[1][0]);
distances[2] = endpoints[0][0].distanceTo(startpoints[1][0]);
distances[3] = endpoints[0][0].distanceTo(endpoints[1][0]);
auto itMin = std::min_element(distances.begin(), distances.end());
size_t ixLF = itMin - distances.begin();
reverse[1][0] = !!(ixLF & 1);
reverse[0][0] = !!(ixLF & 2);
distances[0] = endpoints[1][1].distanceTo(endpoints[0][1]);
distances[1] = endpoints[1][1].distanceTo(startpoints[0][1]);
distances[2] = startpoints[1][1].distanceTo(endpoints[0][1]);
distances[3] = startpoints[1][1].distanceTo(startpoints[0][1]); // should be 0 (or minimal)
itMin = std::min_element(distances.begin(), distances.end());
size_t ixRB = itMin - distances.begin();
reverse[0][1] = !!(ixRB & 1);
reverse[1][1] = !!(ixRB & 2);
for (i = 0; i < 2; ++i) // 0: left/right 1:front/back
{
for (j = 0; j < 2; ++j) // 0: left/front 1:right/back
{
AcDbCurve* curve;
if ((es = acdbOpenObject(curve, idCurve[i][j], AcDb::kForRead)) == Acad::eOk)
{
bool r = reverse[i][j];
double par, parStart, parEnd, dPar;
es = curve->getStartParam(parStart);
es = curve->getEndParam(parEnd);
count = i ? M : N;
dPar = (parEnd - parStart) / (count - 1);
for (int k = 0; k < count; ++k)
{
par = parStart + k * dPar;
es = curve->getPointAtParam(par, pt);
int kr = r ? count-k-1 : k; // reverse or not?
if (i)
p[j ? N - 1 : 0][kr] = pt.asVector();
else
p[kr][j ? M - 1 : 0] = pt.asVector();
}
curve->close();
}
}
}
// Roughly initialize grid
double _N = 1.0 / N, _M = 1.0 / M;
for (i = 1; i < N - 1; ++i)
{
for (j = 1; j < M - 1; ++j)
{
double fn = i * _N;
double fm = j * _M;
p[i][j] = (
(fm * p[i][0] + (1.0 - fm) * p[i][M-1])
+ (fn * p[0][j] + (1.0 - fn) * p[N-1][j])
) / 2.0;
}
}
// Vary grid points until there is no more relevant change
double deltaMax = 0.1;
double delta;
int innerGridCount = (N - 2) * (M - 2);
AcGeVector3d pNew;
do {
delta = 0.0;
for (i = 1; i < N - 1; ++i)
for (j = 1; j < M - 1; ++j)
{
// Calculate a new point p[i][j] by the "average" from p[i-1][j], p[i+1][j], p[i][j-1] and p[i][j+1].
// This simulates p[i][j] being "dragged" to the point with the smallest distance to these neighbour points.
pNew = (p[i - 1][j] + p[i + 1][j] + p[i][j - 1] + p[i][j + 1]) / 4;
delta += (p[i][j] - pNew).lengthSqrd();
p[i][j] = pNew;
}
delta = delta / innerGridCount;
} while (delta > deltaMax);
//AcDbPolyFaceMesh, AcDbSubDMesh
AcGePoint3dArray totalArray(N*M);
for (i = 0; i < N; ++i)
for (j = 0; j < M; ++j)
totalArray.append(asPnt3d(asDblArray(p[i][j])));
AcDbDatabase* pDB = acdbHostApplicationServices()->workingDatabase();
AcDbPolygonMesh* mesh = new AcDbPolygonMesh(AcDb::kSimpleMesh, N, M, totalArray, false, false);
mesh->setDatabaseDefaults(pDB);
AcDbObjectId meshId;
postToDb(pDB, mesh, meshId);
}
Acad::ErrorStatus postToDb(AcDbDatabase *db, AcDbEntity* ent, AcDbObjectId& objId)
{
Acad::ErrorStatus es;
AcDbBlockTable* pBlockTable;
AcDbBlockTableRecord* pSpaceRecord;
if ((es = db->getSymbolTable(pBlockTable, AcDb::kForRead)) != Acad::eOk)
return es;
if ((es = pBlockTable->getAt(ACDB_MODEL_SPACE, pSpaceRecord, AcDb::kForWrite)) != Acad::eOk)
return es;
pBlockTable->close();
es = pSpaceRecord->appendAcDbEntity(objId, ent);
pSpaceRecord->close();
ent->close();
return es;
}
You can change the values of N and M without problems. I made them static to be able to "play around" with them easily while debugging.
I realized that you don't need to create curves to use the code.
You can skip the selection of curves and the following calculations that initialize the "grid border points".
p[0][M-1] ... p[N-1][M-1]
... ...
p[0][0] ... p[N-1][0]
Just initialize them yourself with your mesh border points.
Than you can start with the code after comment // Roughly initialize grid
I don't know what you mean by "equilateral length". Which lenghts shall be equal?
My code generated the white mesh on the left from the arcs and the lines in test.dwg as the "minimum tension surface":
@1127204185 wrote:
N and M remain unchanged, and there is no parameter that can the shape of the surface. For example X=1 Y=1 X=1 Y=2 X=1 Y=3
The final shape of the surface is determined by the "border frame points" of the mesh. The mesh points at the borders are fixed at their positions but all other points are allowed to move in 3D to a place with lower "tension".
I added a buildable sample with functions that allow to change M and N. It implements two commands:
MyCommand / MC = My version that lets you select four curves. The code is modified, so that it is possible to adjust M and N:
// Command MyCommand / MC
void cmdMinSurf()
{
AcDbObjectIdArray curves;
if (!SelectBorderCurveArray(curves))
return;
bool reverse[4];
if (!BuildCurveLoop(curves, reverse))
return;
std::vector< std::vector<AcGeVector3d > > p;
int M = 40, N = 50;
if (!initializeGridFrame(p, M, N, curves, reverse))
return;
double deltaMax = 0.1;
if (!MinimizeSurf(p, deltaMax))
return;
AcDbObjectId meshId;
ShowMesh(p, meshId);
}Note than M and N are basically the array dimensions of p.
Much of my code deals with initializing the "border frame" of p from the selected curves. But you don't have to deal with this. To make things more clear I added another command:
YourCommand / YC
// Command YourCommand / YC
void cmdYourMinSurf()
{
std::vector< std::vector<AcGeVector3d > > p;
int M = 40, N = 50;
double deltaMax = 0.1;
YourGeometryData data;
if (!YourInitializeGridFrame(p, M, N, data))
return;
if (!MinimizeSurf(p, deltaMax))
return;
AcDbObjectId meshId;
ShowMesh(p, meshId);
}Now it is your task to implement the function YourInitializeGridFrame(p, M, N, data) using your own geometry data to shape the surface. I can't do this for you. The project contains a sample implementation that should show you what to do.
