Flattening of a 2D surface

[mfloris]

Thank you so much (I'm going crazy trying to convert the formula!!)

And thanks for CloudCompare, it's an amazing piece of software. Here at work we use 3D software products worth thousands of euros but we always install CC too because in many occasions it's just better and simpler ;)

[daniel]

In fact the equation letters were already properly handled by CC.

I simply added a shift component on the center of gravity (hopefully in the right direction). Can you test the latest beta version?

[mfloris]

Can you test the latest beta version?

I'm on it! ;)

So far:

2.6.1 :

Code: Select all

[09:31:42] [doActionFitQuadric] Quadric equation: z = 0.191049 + 3.54679e-05 * x + 0.00128381 * y + -0.00411162 * x^2 + 0.00332104 * x.y + -0.000473885 * y^2

2.6.2.beta:

Code: Select all

[09:31:28] [doActionFitQuadric] Quadric equation: z = 0.467229 + 3.54679e-05 * x + 0.00128381 * y + -0.00411162 * x^2 + 0.00332104 * x.y + -0.000473885 * y^2

the "x.y" bug is still present (if it is a bug...)

There is little to no difference in the two produced clouds since the only thing that is changing is the constant parameter, that is the "new" cloud is only slightly higher but still doesn't fit.

[mfloris]

Full explanation of what I'm doing:

1) load scanned cloud of a distorted plane
2) fit a plane to the cloud
3) orient the cloud so that the plane is orthogonal to Z
4) fit a 2.5D quadric to the reoriented cloud
5) save the vertices of the quadric to a cloud named "vertices.txt"

run the following VBScript

Code: Select all

Set fso = CreateObject("Scripting.FileSystemObject")

Set objShell = CreateObject("Wscript.Shell")

strPath = objShell.CurrentDirectory &"\"

Set folder = fso.GetFolder(strPath)

Set files = folder.Files

COGX = 0.8245
COGY = 28.636
COGZ = 73.2375

For each folderIdx In files
    
	filename = folderIdx.Name
	name = fso.getbasename(filename)
	extension = fso.GetExtensionName(filename)
	reconstructed = name & "_reconstructed.asc"

	If ( StrComp(extension,"txt",vbTextCompare) = 0 and not fso.FileExists(reconstructed) ) Then

		Set listFile = fso.OpenTextFile(filename)
			
		Set outputFile = fso.CreateTextFile(reconstructed, True)
		
		
			do while not listFile.AtEndOfStream 

				coordinates = split(listFile.ReadLine())
				
				newZ = Quadric(coordinates(0) - COGX, coordinates(1) - COGY) + COGZ
				
				outputFile.WriteLine(coordinates(0) & " " & coordinates(1) & " " & newZ)
			loop
		
			outputFile.Close
		
		End If
	
Next

Function Quadric(ByVal x, ByVal y)
	Dim z

	z = 0.467229 + 3.54679e-05 * x + 0.00128381 * y + -0.00411162 * x^2 + 0.00332104 * x * y + -0.000473885 * y^2
 
 	Quadric = z
End Function

COGX, COGY and COGZ are not che actual COG, they are merely the center of the bounding box of the original input cloud at point (3)

Maybe that's the problem

[daniel]

Forget about what I said yesterday! I totally missed the fact that we are already projecting the input points on their best fit plane prior to compute the best quadric. The result is that the equation is expressed in a local coordinate system that has nothing to do with the original one.

I've updated CC so as to display the matrix that will let you project the points between the two coordinate systems. And by the way, now that we do this I guess the dimensions/equation will always be z = f(x,y) (as the local coordinate system is meant to minimize the spanning along Z).

This time I tested everything ;)

cc_quadric_projection.jpg (64.11 KiB) Viewed 5272 times

And here is my pseudo code:

Code: Select all

//the 4x4 transformation displayed in the Console
const ccGLMatrix& trans = quadric->getTransformation();
//the inverse transformation
ccGLMatrix invTrans = trans.inverse();

//for each point
for (unsigned i=0; i<newCloud->size(); ++i)
{
	//original point
	CCVector3* P = const_cast<CCVector3*>(newCloud->getPoint(i));
	//point in the local coordinate system
	CCVector3 Q = invTrans * (*P);
	//project the point on the quadric
	Q.z = eq[0] + eq[1]*Q.x + eq[2]*Q.y + eq[3]*Q.x*Q.x + eq[4]*Q.x*Q.y + eq[5]*Q.y*Q.y;
	//back project to the original coordinate system
	*P = trans * Q;
}

Sorry for the confusion.

P.S.: you don't need it anymore, but the translation of the transformation matrix (the 4th column) is the cloud gravity center ;)

[daniel]

Oups, I just realized that I put a wrong version of CC online yesterday. That's fixed now ;)

[mfloris]

Ok, I can see more information now... I just have to figure out what to do with it! :)

[matt_w]

Hi guys,

I've stumbled on this thread whilst trying to do the same thing, but I'm afraid that you lost me about halfway through with the technicalities of it all.

Put simply, I'd like to detrend my point cloud using the quadric (below) before running some basic stats on it, e.g. detrended standard deviation of elevations, etc. I've already fitted a standard plane and oriented the data so that this plane fits the x,y dimension of the data (as advised), but I'm struggling to figure out how to then use the quadric equation/transformation to manipulate the point cloud to essentially remove the quadric curve signal.

Would you be able to explain the steps in a simple way, please?

Many thanks, keep up the good work.

Capture.JPG (53.96 KiB) Viewed 5172 times

[daniel]

Well, you may try this trick:
- compute the distances between your point cloud and the quadric (Tools > Distances > Cloud to Mesh)
- then convert the output 'C2M Distances' scalar field into the Z coordinate of the cloud (Edit > Scalar fields > Set SF as coordinates)

THERE'S A BIG APPROXIMATION CONCERNING X AND Y (they should be displaced as well). If your quadric is too curvy, then the distortion will be high.

It will also be slightly approximated because the quadric mesh is tessellated. But it's less severe and you can increase the quadric 'drawing precision' before computing the distances to increase the accuracy (depending on how important it is for your application).

P.S.: once done, you can remove the scalar field (or simply restore the RGB colors as the default color source).

[matt_w]

Thanks Daniel - I just tried something similar to what you suggest before I saw your reply. Good to know I've been thinking along the right lines :-)