PoissonRecon --in bunny.points.ply --voxel bunny.voxel.ply --depth 10

My understanding is that the line above outputs bunny.voxel.ply file that contains 3D array (voxel grid), where each element is a value of the implicit function.

I tried reading bunny.voxel.ply with pcl library, but getting read error... Even though with the same code I am able to read the mesh file bunny.ply produced by

PoissonRecon --in bunny.points.ply --out bunny.ply --depth 10

So my question is how can I read that 3D array from bunny.voxel.ply?

The version 9.0 binaries for PoissonRecon and SSDRecon from:

http://www.cs.jhu.edu/~misha/Code/PoissonRecon/Version9.0

crash on Windows 10 with an Illegal Instruction. After fixing issue #14 locally I'm able to produce a debug build that does not crash and runs as expected. However, my local release builds are crashing as well. The exception is thrown at line 241 in Ply.h.

Problem traced to _GetSampleDepthAndWeight() in MultiGridOctreeData.WeightedSamples.inl

Suggested fix:
if(oldWeight>0&&newWeight>0){
double denom=log( newWeight / oldWeight ) ;
if(denom!=0){
depth = Real( _Depth( temp ) + log( newWeight ) / denom);
}
} else {
depth = Real( _Depth( temp ));
}

Hello,

I have some colored point clouds and I use this code to generate the mesh using the PoissonRecon and SurfaceTrimming codes.

I also have some material ids in the point cloud. For example each point in the cloud has, position(x,y,z), color(r,g,b), normals(nx,ny,nz) and Material ID(m_id). After the reconstruction I surely lose the material ID. Could you please tell me what would be the best way to keep the material id even after I construct the mesh?

It would be very useful because, I can get the material id while the data is in image level and definitely can use a mapping to pass this information till point cloud. If I could pass it though mesh generation it would allow me having material maps in UV.

Thanks,
Ahmad

I saw that there were a lot of changes a day ago and pulled the latest version. But the code does not compile at this moment. The error I get is:

g++ -o Bin/Linux/PoissonRecon Bin/Linux/PlyFile.o Bin/Linux/PoissonRecon.o -lgomp -lstdc++ -lz -lpng -ljpeg -O3 -g
g++ -o Bin/Linux/SSDRecon Bin/Linux/PlyFile.o Bin/Linux/SSDRecon.o -lgomp -lstdc++ -lz -lpng -ljpeg -O3 -g
g++ -o Bin/Linux/SurfaceTrimmer Bin/Linux/PlyFile.o Bin/Linux/SurfaceTrimmer.o -lgomp -lstdc++ -lz -lpng -ljpeg -O3 -g
g++ -o Bin/Linux/EDTInHeat Bin/Linux/PlyFile.o Bin/Linux/EDTInHeat.o -lgomp -lstdc++ -lz -lpng -ljpeg -O3 -g
mkdir -p Bin/Linux/
g++ -c -o Bin/Linux/ImageStitching.o -fopenmp -Wno-deprecated -Wno-write-strings -std=c++11 -O3 -DRELEASE -funroll-loops -ffast-math -g -I. Src/ImageStitching.cpp
In file included from Src/PNG.h:38:0,
                 from Src/Image.h:110,
                 from Src/ImageStitching.cpp:41:
Src/PNG.inl: In constructor ‘PNGWriter::PNGWriter(const char*, unsigned int, unsigned int, unsigned int, unsigned int)’:
Src/PNG.inl:115:40: error: ‘Z_BEST_SPEED’ was not declared in this scope
  png_set_compression_level( _png_ptr , Z_BEST_SPEED );
                                        ^~~~~~~~~~~~
Src/PNG.inl: In member function ‘virtual unsigned int PNGWriter::nextRows(const unsigned char*, unsigned int)’:
Src/PNG.inl:147:131: error: invalid use of incomplete type ‘png_struct {aka struct png_struct_def}’
  for( unsigned int r=0 ; r<rowNum ; r++ ) png_write_row( _png_ptr , (png_bytep)( rows + r * 3 * sizeof( unsigned char ) * _png_ptr->width ) );
                                                                                                                                   ^~
In file included from Src/PNG.h:7:0,
                 from Src/Image.h:110,
                 from Src/ImageStitching.cpp:41:
/usr/include/png.h:470:16: note: forward declaration of ‘png_struct {aka struct png_struct_def}’
 typedef struct png_struct_def png_struct;
                ^~~~~~~~~~~~~~
Makefile:149: recipe for target 'Bin/Linux/ImageStitching.o' failed
make: *** [Bin/Linux/ImageStitching.o] Error 1

We've been loving this library but recently realized that our code is executing in an environment where disk space is shared across multiple requests. For compliance purposes, we need to either use encrypted files on disk (would definitely require a fork), switch how we're passing information to PoissonRecon, or something else

One promising option we can think of is to pass in files via stdin and receive the result via stdout (leaving stderr to report any issues). There's at least 1 rabbit hole down the avenue but aside from that it seems straightforward:

• We currently sniff filenames to determine how to parse the data

  • PoissonRecon/Src/PoissonRecon.cpp

    Lines 478 to 481 in 16375a7

    	sampleData = new std::vector< ProjectiveData< Point3D< Real > , Real > >();
    	if ( !strcasecmp( ext , "bnpts" ) ) pointStream = new BinaryOrientedPointStreamWithData< Real , Point3D< Real > , float , Point3D< unsigned char > >( In.value );
    	else if( !strcasecmp( ext , "ply" ) ) pointStream = new PLYOrientedPointStreamWithData< Real , Point3D< Real > >( In.value , ColorInfo< Real >::PlyProperties , 6 , ColorInfo< Real >::ValidPlyProperties );
    	else pointStream = new ASCIIOrientedPointStreamWithData< Real , Point3D< Real > >( In.value , ColorInfo< Real >::ReadASCII );

  • Most likely easily resolved if we add another flag for how to parse input (and probably how to output as well)

We'd be more than happy to help out with implementing a solution to this issue =)

When the Depth and others parameters is certain, pointweight is respectively 1,2,4,8,16,the corresponding CPU time is 73.53s，73.18s.72.63s， 77,78s and 77.56s. Shouldn't CPU time like a logarithmic increase?Can you help me to answer this questions，thanks.

Is there any version of the code which contains enouth comments?
I'm sorry that I have difficulty in reading the code. So I desire to know whether there are more comments about the code.
Thanks!

Hi. In the Poisson surface reconstruction paper 2006, you mentioned block Gauss-Seidel solver to solve the issue of large memory requirement of solving the Lx = b equation (size of sparse matrix L could be 125 * octree node numbers in this depth d).

To address this issue, we augment our solver with a block Gauss-Seidel solver. That is, we decompose the d-th dimensional space into overlapping regions and solve the restriction of Ld to these different regions, projecting the local solutions back into the d-dimensional space and updating the
residuals. By choosing the number of regions to be a function of the depth d, we ensure that the size of the matrix used by the solver never exceeds a desired memory threshold.

I wonder how to implement this block Gauss-Seidel solver, how to decompose the large matrix L into small ones and solve them one by one ? Could you kindly give some guidance ? I don't find much information or code on the net.

I'm quite sure that the output PLY file of PoissonRecon included the normals. Now the output file doesn't contain the points normals. What could have happened?
Can you confirm that the output file should contain the normals?
Thank you.

Hi,

I downloaded source code in zip from here:
http://www.cs.jhu.edu/~misha/Code/PoissonRecon/Version9.011/

I built it in Visual Studio 2017.
Build was successfull.

I see .exe files, but when I run those three exe files they are closed immediatelly.

How should I run the app correctly?
I have a pointcloud and just want to reconstruct a mesh from it.

Hi.
I'm trying to run the last version (binary distibution of v8, 64 bit) on a Windows 10 system. I have installed the VS2013 redist package. The application crashes on startup (even with no arguments).
Instead the 32 bit version seems to work, if I run it with no arguments it prints the usage.
There are some dependencies needed?

I tried to use the newest version of the PoissonRecon but it simply doesn't execute.
I compiled under Win10 VS2017 with Windows SDK 10.0.10586.0.

I only changed the Windows SDK because I don't have the default one and I downloaded the test data from the code website.

The sentences of this pic means : The PossionRecon stops work.

I also tried the older version and everythign works.

Hi,
I'm wondering why faces list at ply output is in binary format. Is it set by default? Can it be saved as ascii?

Input file: http://s000.tinyupload.com/?file_id=66724387703494627365
Output file: http://s000.tinyupload.com/?file_id=08871875982492369292

--in C:\dcmtest\dicom_222906.ply --depth 10 --out C:\dcmtest\output_dicom_222906.ply --verbose --pointWeight 0

Can you mark some stable version? I'm trying to use it in my master's degree project :)

Hello,

I'm trying to collect values inside the execute function of surfacetrimmer.cpp

In version 9 it was like

   for (size_t i = 0; i < vertices.size(); i++)
    {
      vertz.push_back(vertices[i].point.coords[0]);
      vertz.push_back(vertices[i].point.coords[1]);
      vertz.push_back(vertices[i].point.coords[2]);

      clrz.push_back(vertices[i].color[0]);
      clrz.push_back(vertices[i].color[1]);
      clrz.push_back(vertices[i].color[2]);
    }

In version 10.00 after including roughness as you instructed this part changed to

for (size_t i = 0; i < vertices.size(); i++)
    {
      vertz.push_back(vertices[i].point.coords[0]);
      vertz.push_back(vertices[i].point.coords[1]);
      vertz.push_back(vertices[i].point.coords[2]);

      clrz.push_back(vertices[i].color()[0]);
      clrz.push_back(vertices[i].color()[1]);
      clrz.push_back(vertices[i].color()[2]);

      roughness.push_back(vertices[i].roughness());
    }

But now in 10.01 it does not work the same way I believe and if I print out this:

    for (size_t i = 0; i < vertices.size(); i++)
    {
       std::cout << vertices[i].point.coords[0] << " "
                 << vertices[i].point.coords[1] << " "
                 << vertices[i].point.coords[2] << " "
                 << vertices[i].point.coords[3] << " "
                 << vertices[i].point.coords[4] << " "
                 << vertices[i].point.coords[5] << " "
                 << vertices[i].point.coords[6] << " "
                 << vertices[i].point.coords[7] << " "
                 << vertices[i].point.coords[8] << " "
                 << vertices[i].point.coords[9] << " "
                 << vertices[i].point.coords[10] << " "
                 << vertices[i].point.coords[11] << " "
                 << vertices[i].point.coords[12] << " "
                 << vertices[i].point.coords[13] << std::endl;

      }

I see the values of position, color, normals etc. but I'm not sure how does the coords[i] get it's sequence from. It looks like it prints points first then roughness then normals etc. Will this sequence be the same if I do not use a property as i.e. normals? How does it work?

Thanks,
Ahmad

Hi,
In PlyFile.cpp: ply_type_sizes,

1. float 32 is set to be 8 bytes
2. an entry for float64 is missing.

Hello there,

I have a question about the PoissonRecon command-line program. I'm interested in sampling the output voxel grid (i.e., output by the --voxel flag) over the output triangle mesh. In other words, I'd like to color each vertex of the output triangle mesh according to the corresponding value in the voxel grid.

I am aware that the output voxel grid sampled over the output triangle mesh should always be 0, because the triangle mesh is exactly the 0 level set of the voxel grid. But I still want to be able to sample the voxel grid in this way, because I ultimately want to want to compute derived quantities over the voxel grid, and sample these derived quantities over the triangle mesh.

Sampling the voxel grid in this way requires care, because the voxel indices need to be mapped into the coordinate system of the output triangle mesh.

In my code, I assume that the axis-aligned bounding box (AABB) of the voxel grid is exactly 1.1x the size of the AABB of the input point cloud along each dimension (i.e., I am using the default value for the --scale parameter). I also assume that the AABB of the voxel grid has the same center as the AABB of the input point cloud. Together, these assumptions uniquely determine the mapping from triangle mesh coordinates to voxel indices.

These assumptions also seem to imply that the AABB of the output triangle mesh should be at most 1.1x the size of the AABB of the input point cloud along each dimension. But I am finding that this isn't the case. I am getting a larger-than-expected AABB for the output triangle mesh. In my particular case, the size of the of the output triangle mesh's AABB along each dimension, relative to the input point cloud's AABB is [ 1.02068782 1.16858709 1.34913611]. Note that the relative size of the output triangle mesh's AABB is noticeably larger than the maximum expected relative size of 1.1x along the last two dimensions.

Is the voxel grid's AABB indeed 1.1x the size of the input point cloud's AABB? If so, where are these output vertices coming from that are so far outside the AABB of the voxel grid? I understand that PoissonRecon implements marching cubes, but as far as I am aware, marching cubes only produces triangles that are interior to a voxel grid (i.e., marching cubes does not extrapolate).

This behavior is potentially problematic, because it suggests that I'm not aligning my voxel grid and my output triangle mesh correctly.

I'm using PoissonRecon version 9.0 on Mac, and I'm invoking it as follows:

./PoissonRecon --in pset_mvs_si.ply --out poisson_surface_si.ply --voxel poisson_function_si.bin --depth 9 --color 16 --density --verbose

In case it's helpful, I've uploaded my input point cloud and output triangle mesh below.

Input point cloud: https://www.dropbox.com/s/ybhjsj3iejiix7c/pset_mvs_si.ply?dl=0
Output triangle mesh: https://www.dropbox.com/s/8oqe0jtkxujohqi/poisson_surface_si.ply?dl=0

Cheers,
Mike

We are using PoissonRecon in AWS Lambda which only permits creating files in its /tmp directory (OS' temporary directory):

http://docs.aws.amazon.com/lambda/latest/dg/limits.html

Currently BufferedReadWriteFile writes to the current working directory. As a workaround for now, we've hardcoded the BufferedReadWriteFile template to use the /tmp directory:

https://github.com/StandardCyborg/PoissonRecon/blob/93762c30d870752220973abb8252ca3b95e88d14/Src/Geometry.cpp#L52

I believe it's possible to find the OS' temporary directory via environment variables but am rather inexperienced in C++ and wanted to check the PR would be accepted first. If I wrote a PR to switch over to the OS' temporary directory, would you be interested in it in PoissonRecon?

Hello there,

I have a quick silly question about the volume returned by SSDRecon with the --voxel option. How closely does this function approximate a distance function, assuming noise-free synthetic input data and default parameter settings?

When I examine the output volume I get using my noise-free synthetic input data, I see that the distance values appear to be "normalized", in the sense that they seem to be between -1 and 1 for my data. Based on this observation, I am assuming that the largest distance value you'd ever expect to see in the output volume is sqrt(3), because this would be the normalized distance from the origin to the most distant corner of the unit bounding box. So far, so good.

However, if I assume that the volume extends from (0,0,0) to (1,1,1), then the distance values do not seem to approximate a distance function very closely. For example, here is a slice of the distance function I obtain from my input data. I am not doing any of my own scaling here, I am simply displaying the raw values from the output volume:

The zero level set of this function approximates my surface very accurately:

But when I numerically compute the gradient magnitude of this function, it seems to be non-constant and less than 1 almost everywhere. Here is a slice of the gradient magnitude, which I computed in 3D:

Assuming that I call SSDRecon with default parameter settings, how strictly is the constraint enforced that the gradient magnitude is constant? Am I wrong to expect that the distance function returned by SSDRecon is approximately correct far away from the surface?

Hi,@mkazhdan
First,thank you for providing the PoissonRecon project.
I recently used your code and made some changes to it.I used the MFC interface to run the code,rather than the console cause I found it will take part of the memory when the end of the operation.So I have some problems.When I first choose the point-set file ,run the code and it run perfectly,however when I choose the file(the same as the first or not)and run the code,it collapsed(I did not withdraw MFC in this process).
I try to debug this code,but I could not find th reason.My compiler is visual studio 2013.Thank you.

I put together a small CMake script for building your library. If interested in adding it please find it attached.

Thank you for this great library!

CMakeLists.txt

I have two questions about the algorithm. I'd be appreciated if you could offer some information.

1. Is there a paper or something to explain the algorithm of the surface trimming method ?
2. Do you know any methods related to the GPU implementation of your screened Poisson reconstruction algorithm ?

Hi , I have not read the code much. Instead I read the papers (the original one, screened one, and the GPU implementation one by MSR Asia). I implemented all the codes mainly according to the MSR paper (in CPU codes, not GPU), which is based on your paper published in 2006. But sadly I have not got the desired mesh and cannot find the problem by far. My mesh is distorted along the boundary . I have two questions. Hope you could give me some direction.

1. How do we impose the boundary constraints (Dirichlet or Neumann) to the linear system Lx = b? The descriptions about boundary constraints just appeared in screened poisson paper. And I don't know where to implement this constraints.
2. For uniformly distributed samples, is the way to compute the implicit function value is sum(F * ϕ) ? Right now I use the following algorithm given by MSR to compute the isovalue of an arbitary point:

Listing 4 Compute Implicit Function Value ϕq for Point q
1: ϕq = 0
2: nodestack ← new stack
3: nodestack.push(proot)  // push root node (depth 0 node)
4: while nodestack is not empty
5:     o ← NodeArray[nodestack.pop()]
6:     ϕq+ = Fo(q)ϕo
7:     for i = 0 to 7   // traverse 8 children of this node
8:         t ← NodeArray[o.children[i]]
           // the node function range of each node is set to be [-width, width]
9:         if abs(q.x −t.x) < t.w and abs(q.y −t.y) < t.w and abs(q.z −t.z) <t.w then
10:           nodestack.push(o.children[i])

I use this algorithm to compute isovalue in sample points and cube vertices.
This simple algorithm just iterate the nodes whose node function cover the point and accumulate F*ϕ.
But I think since the nodes are discrete and sparse, the total number of nodes that contribute ϕ to the arbitary point is not stable and can always vary, which may lead to the incontinuity of implcit function value ϕ in the 3D domain. Let alone finer depth nodes cannot cover most vertices of coarser depth leaf nodes. By the way, since the magnitude of node fuction is scaled according to node depth ( 1 / (width^3)), the contributions of coarse depth nodes are very low, and the implicit function value may change greatly in vertices from different depth nodes. I really suspect whether this method could result in a continuous, smooth and watertight isosurface. Is this algorithm enough for computing isovalue of sample points or cube vertices? I'd appreciate a lot if you could explain it to me :)

Hello, I'd like to ask about two questions from your original paper(2006), mainly related to implement with non-uniform samples case:

1. is an evaluated value of local sampling density is formed in scalar?

WˆD(q) ≡ Σs∈SΣo∈NgbrˆD(s) αo,s*Fo(q)
(sorry for readability why it doesn't rendered :| )

• I thought that every node function Fo(q) will returns scalar since it can be separated into three functions, and the evaluation can be done by multiplying all separated function's value for each queried point's coordinate.
• For the case of αo,s, I'm understanding that it'll return 3D vector type value since it'll return three interpolation parameters like tx, ty, tz.
• Finally, It seems obvious that the resulting density WˆD(q) should return a scalar value, since it's literally a density... and later when we calculate a desired depth Depth(s.p), density from a sample point is divided by average density. (Of course the average seems to be formed in scalar too). Furthermore in isovalue selection, the inverse of density is being used, which can't be done with vector form.
• However, as I tried to derive the expression, the multiplication part αo,s*Fo(q) seemed return 3d vector since it's multiplication between 3d vector and a scalar...

Did I missed something with deriving the equation? How can I derive it correctly?

2. I wonder why sample's normal didn't used for calculating vector field for non-uniform samples case. I mean, in here:

V(q) ≡ Σs∈S (1 / WˆD (s.p)) Σo∈Ngbr(Depth(s.p))(s) (αo,sFo(q))

• Should I consider sample's normal manually when I try to solve linear system?

If more description is required, I'll explain it more.

Thanks in advance!
PilJoong

Hi, I found that the marching cube codes are not really used for mesh generation. What kind of mesh generation algorithm do you use ?

Using Visual Studio 2015 doing a debug build fails with the message:

PoissonRecon.cpp : fatal error C1128: number of sections exceeded object file format limit: compile with /bigobj
After adding the /bigobj compilation option the build succeeds and the binary can be run successfully.

I tried compiling on Linux and running it, but it never seems to complete and uses 100% of all my cores. I am trying the bunny example with a depth of 9

Hello,
I saw this update in the latest diff :

Modified the reconstruction code to facilitate interpolation of other input-sample quantities, in addition to color

Although I followed your suggestion and was able to finally integrate the roughness property to the code but if you include it in the main repo I will just stash mine. But unfortunately I could not find an updated README to follow on how to run with the additional property. Could you please update the README or let me know here on how to use this new functionality? Once I can build the 10.01 version I will be able to start testing. I'll have to make it CMake compatible after.

Thanks,
Ahmad

Hello Sir
I want to know the exact effect of --depth on the out put reconstruction and how the mechanism work.
i mean is the lower the voxel rez the less point you take into consideration when reconstructing a surface ?

Than you

Saouli

When I run my model with depth 10 (color=16), the resulting color is somewhat blurry (detail level is excellent).
Is there any way to influence the interpolation process to limit it's spatial extent?

I've attempted to make it sharper by:

1. In MultiGridOctreeData.WeightedSamples.inl : Octree< Real >::_evaluate(), I limited the evaluation to the last depth levels.
2. Changing the BSpline degree (DATA_DEGREE in MultiGridOctreeData.h). Hoping that higher order polynomials, would result in higher sharpness.
   Both were unsuccessful.

Do I have any alternative to using depth 11?

Under linux mint 18.3 bash:
$ ./PoissonRecon --in horse.npts --out horse.ply
$ ./EDTInHeat --in horse.ply --out horse.edt
then, I get the error : Segmentation fault (core dumped)

Under the windoows 7 x64, I get the same error.

I have tried with the binary downloaded from github and compiled by myself, the error is same.

Can someone help me?

Hi,

I want to link this code with my own program to do some research experiments, but I find it is hard to fill points and normals easily by passing ,like two vectors, pts, normals. Could you help release some filling data interface?

Thanks a lot

Hello,
I'm using the submodule functionality of git to include your repo in a project (MeshLabJS), and I don't want to modify your code. Also, I don't want to use the omp library. Thus, I don't define the macro _OPENMP, but in the file MultiGridOctreeData.h, row #57, the library is not enclosed in #ifdef _OPENMP, unlike the functions in rows 65-66 on the same file, and compilation fails.

Hello!
I'm using poisson mesh generation on Android devices, and i run into the following issue multiple times:

Tree set in: 2.9 (s), 0.2 (GB)
Input Points: 86490
Leaves/Nodes: 397769/454593
Memory Usage: 0.195 GB
Constraints set in: 1.3 (s), 0.2 (GB)
Memory Usage: 0.195 GB
Depth[0/8]: 8
Evaluated / Got / Solved in: 0.000 / 0.000 / 0.000 (0.195 GB)
Depth[1/8]: 64
Evaluated / Got / Solved in: 0.000 / 0.000 / 0.001 (0.195 GB)
1Depth[2/8]: 512
Evaluated / Got / Solved in: 0.000 / 0.001 / 0.002 (0.195 GB)
Depth[3/8]: 4096
Evaluated / Got / Solved in: 0.000 / 0.009 / 0.004 (0.195 GB)
Depth[4/8]: 32768
Evaluated / Got / Solved in: 0.000 / 0.073 / 0.016 (0.195 GB)
Depth[5/8]: 262144
Evaluated / Got / Solved in: 0.000 / 0.620 / 2.431 (0.195 GB)
Depth[6/8]: 54352
1Evaluated / Got / Solved in: 0.000 / 0.304 / 1.266 (0.198 GB)
Depth[7/8]: 100016
Evaluated / Got / Solved in: 0.000 / 0.354 / 0.382 (0.199 GB)
Depth[8/8]: 632
Evaluated / Got / Solved in: 0.000 / 0.030 / 0.835 (0.199 GB)
Linear system solved in: 6.5 (s), 0.2 (GB)
Memory Usage: 0.199 GB
Got average in: 0.086310
Iso-Value: 3.226770e+00
[ERROR] Failed to close loop [8: 144 36 42] | (454465): 756464301900354

(You may notice i altered the printout for memory usage, as the currently used reading of '/proc/self/stat' resulted in gibberish on Android. I'm using 'getrusage' to retrieve the 'ru_maxrss' value for consumption tracking)

What bothers me is that its non-deterministic. The issue comes up around 3 times in 10 run. I tried to resolve the issue with the following:

• Compile the code without funroll-loops and ffast-math
• Turn on lInear interpolaion
• Turn on double precision

None of the above resolved the error. Do you have any suggestion what went wrong and how to eliminate it?

When building PoissonRecon, my compiler informs me of this message, and I didn't see an issue tracking it. Are there more details on the issue behind this warning?

Hi Michael,

I saw that you have "-1" _tree node in Octree<> class,
and i was wondering what is the main purpose of additional nodes (_tree node as a container of the _spaceRoot and its 7 siblings, therefore all points from stream are contained in _spaceRoot node only). There are a lot of checks: such as

if(depth > 0)
or
if(depth > _minDepth)

and some inconvenient variables such as "_minDepth" and "_maxDepth" which differs from actual input parameter Depth by 1.

Could you give some more insight on it?

Thanks in advance.

Hi.
I'm trying to run the last version on a Windows 7 system. I have installed the VS2013 package. while running the Lucy model ,there is wrong that Failed to find property in ply file: nx. but eagle.points.ply's
property containing nx,ny,nz is ＯＫ

I would like to be able to run the Poisson code in a loop over multiple meshes, but running more than once in the same process causes the program to crash.

Specifically, line 204 of Octree.inl fails because the variable _depthAndOffset is not set.
inline int OctNode< NodeData >::depth( void ) const {return int( _depthAndOffset & DepthMask );}

What needs to be modified so that the code can be run in a loop?

hello
I'm trying to use the aRatio parameter to get rid of "blobs" around my reconstruction.
but actually I'm getting the opposite of what I intuitively thought would happened , ie the results is like keeping the blobs only.is there a way to remove the little disconnected area and keep the surface ?

commands:
SurfaceTrimmer --in MVE/Poisson-L2-node1.ply --trim 5 --out MVE/Poisson-L2-trim5

SurfaceTrimmer --in MVE/Poisson-L2-node1.ply --out MVE/Poisson-L2-trim5-aRatio1 --aRatio .5

thanks in advance

Hello Misha,
There are some questions about the code of "Screened Poisson Surface Reconstruction " . They are:
(1) Why do you normalize the input data as mapping the data into the 0-1 range when you construct Octree?
(2) I have found that the functions SetSliceIsoCorners(...) , SetSliceIsoVertices(...) and so on are executed on the slices that located between 0-0.5 range at all depth, but not 0-1. I mean corners and iso-vertices are gotten by going through a part of slices, but not all slices, and just process the nodes whose returned value of function _IsValidNode< 0 >( ...) is true, which makes me confused. Could you please tell me why set node flags? And why do not go through all slices or node to get corners and iso-vertices?
(3) The reconstructed surfaces are different when using Dirichlet and Neumann boundary conditions, why does Neumann boundary condition could result in an open surface, but the Dirichlet boundary condition leads to a closed surface?
And I have other two questions that are not related to the code, what are the differences for constructing surface between triangular mesh and quadrilateral mesh? And quadrilateral mesh could be a wider application?

Cheers,
Jack

Hello Misha,
The code :
if( !neighborKey.neighbors[depth].neighbors[1][1][2z] || !neighborKey.neighbors[depth].neighbors[1][1][2z]->children )
in the function SetSliceIsoEdges(...) ,
and the code:
if(!xValues.faceSet[ eIndices[e] ] && ( !neighborKey.neighbors[depth].neighbors[xx][yy][zz] ||
!neighborKey.neighbors[depth].neighbors[xx][yy][zz]->children ) )
in the function SetXSliceIsoEdges(...).
I do not understand why you must assess whether the neighbors are empty in the if statement. Could you please explain it?

Hello, first of, greatly appreciated to providing your implementations!

I'm trying to implement my own version of Screened Poisson Surface Reconstruction, however I'm currently stuck on the problem with implementing some notation.

In your original paper(2006), there're a bunch of notation 'q', and it seems like it means an arbitrary 3D position.
We have an example of node function , which is unit integral base function aligned at node o. This notation is also can be seen in various expressions: Kernel density, Vector field, Isosurface extraction etc (Even for solving linear systems?). Unfortunately, I couldn't find any descriptions for q.

I assumed it as sample's position, however there's already a notation s.p which has equivalent meaning. It also can't be octree node's position since there's notation o.c defines node center.

1. Can I ask for the meaning of the notation q?

ps.

2. Will there be no problem if I substitute Base function (n-th convolution of box filter: approximates to Gaussian distribution) as Unit Variance Gaussian distribution? I thought you used convoluted filter for efficiency & only wants to use in specific range (which is different point against Gaussian, since Gaussian is defined in infinite range), however I'm wondering is there any other reason you used approximated Gaussian.

Thanks for reading, and I'll looking forward for the answers! If more explanation is required, I'll post it more!

Regards,
PilJoong

make debug leads to:

mkdir -p Bin/Linux/
g++ -c -o Bin/Linux/PlyFile.o -fopenmp -Wno-deprecated -Wno-write-strings -std=c++11 -DDEBUG -g3 -I. Src/PlyFile.cpp
mkdir -p Bin/Linux/
g++ -c -o Bin/Linux/PoissonRecon.o -fopenmp -Wno-deprecated -Wno-write-strings -std=c++11 -DDEBUG -g3 -I. Src/PoissonRecon.cpp
g++ -o Bin/Linux/PoissonRecon Bin/Linux/PlyFile.o Bin/Linux/PoissonRecon.o -lgomp -lstdc++ -lz -lpng -ljpeg 
Bin/Linux/PoissonRecon.o: In function `void FEMTree<3u, float>::finalizeForMultigrid<2u, FEMTree<3u, float>::HasNormalDataFunctor<UIntPack<6u, 6u, 6u> >, SparseNodeData<Point<float, 3u>, UIntPack<6u, 6u, 6u> >, FEMTree<3u, float>::DensityEstimator<2u> >(int, FEMTree<3u, float>::HasNormalDataFunctor<UIntPack<6u, 6u, 6u> >, SparseNodeData<Point<float, 3u>, UIntPack<6u, 6u, 6u> >*, FEMTree<3u, float>::DensityEstimator<2u>*) [clone ._omp_fn.5]':
.../PoissonRecon/Src/FEMTree.inl:491: undefined reference to `StaticWindow<RegularTreeNode<3u, FEMTreeNodeData, unsigned short>*, UIntPack<5u, 5u, 5u> >::Size'
Bin/Linux/PoissonRecon.o: In function `void FEMTree<3u, float>::finalizeForMultigrid<2u, FEMTree<3u, float>::HasNormalDataFunctor<UIntPack<7u, 7u, 7u> >, SparseNodeData<Point<float, 3u>, UIntPack<7u, 7u, 7u> >, FEMTree<3u, float>::DensityEstimator<2u> >(int, FEMTree<3u, float>::HasNormalDataFunctor<UIntPack<7u, 7u, 7u> >, SparseNodeData<Point<float, 3u>, UIntPack<7u, 7u, 7u> >*, FEMTree<3u, float>::DensityEstimator<2u>*) [clone ._omp_fn.165]':
.../PoissonRecon/Src/FEMTree.inl:491: undefined reference to `StaticWindow<RegularTreeNode<3u, FEMTreeNodeData, unsigned short>*, UIntPack<5u, 5u, 5u> >::Size'
Bin/Linux/PoissonRecon.o: In function `void FEMTree<3u, float>::finalizeForMultigrid<2u, FEMTree<3u, float>::HasNormalDataFunctor<UIntPack<8u, 8u, 8u> >, SparseNodeData<Point<float, 3u>, UIntPack<8u, 8u, 8u> >, FEMTree<3u, float>::DensityEstimator<2u> >(int, FEMTree<3u, float>::HasNormalDataFunctor<UIntPack<8u, 8u, 8u> >, SparseNodeData<Point<float, 3u>, UIntPack<8u, 8u, 8u> >*, FEMTree<3u, float>::DensityEstimator<2u>*) [clone ._omp_fn.301]':
.../PoissonRecon/Src/FEMTree.inl:491: undefined reference to `StaticWindow<RegularTreeNode<3u, FEMTreeNodeData, unsigned short>*, UIntPack<5u, 5u, 5u> >::Size'
collect2: error: ld returned 1 exit status
Makefile:127: recipe for target 'Bin/Linux/PoissonRecon' failed
make: *** [Bin/Linux/PoissonRecon] Error 1

P.S. Release builds Ok.

Hi Michael,

Can you please explain the order of the nodes for a given slice before MC algorithm.
given to the poisson sphere ptCloud, I'm dumping nodes of the middle slice just before the MC (Marching Squares), having problem to understand the spatial order of sValues.mcIndices[i - sValues.sliceData.nodeOffset] hence corners indicator data.

Thanks,
david

I tested the Jaffa dataset.
very good result.

How to create label.png on ImageStitching?

Hi. I think neighbors of octree nodes should be used in the algorithm. How do you compute octree node neighbors in a fast way ? Do you save 1-ring or 2-ring neighbors of each node?

Hi. Many days ago I asked a question about how to compute node neighbors. You said that when computing neighbors in a certain depth (level), all node neighbors in lower (parent) depth should be known. I think much memory will be consumed when storing all node 2-ring (125 in total) neighbors, since you may use 2nd order B-splines, where the suppor function range is [-1.5, 1.5]. However I noticed that your program consumes really small memory and computes very fast! I'm really amazed at this memory control and speed. I wonder how could you control the memory so well ?

Hi,

I am trying to use Regard 3D which seems to use the PoissonRecon code. It shows success on everything else but it always crash when trying to create surface.
I am using Windows 8.1.
If I click 'debug', it open visual studio and shows this:

Unhandled exception at 0x00007FF673747FF2 in SurfaceTrimmer.exe: 0xC0000005: Access violation reading location 0x000000000000000C.

Not sure what other information I could provide but here is the console output:
Using the OPENMP thread interface

• Features Loading -
  0% 10 20 30 40 50 60 70 80 90 100%
  |----|----|----|----|----|----|----|----|----|----|

• Features Loading -
  0% 10 20 30 40 50 60 70 80 90 100%
  |----|----|----|----|----|----|----|----|----|----|

Track building

Track filtering

Track export to internal struct

Track stats

-- Tracks Stats --
Tracks number: 449
Images Id:
3, 4, 5, 6, 7,

TrackLength, Occurrence
2 433
3 16

A-Contrario initial pair residual: 77.0975

Bundle Adjustment statistics (approximated RMSE):
#views: 2
#poses: 2
#intrinsics: 1
#tracks: 350
#residuals: 1400
Initial RMSE: 0.906527
Final RMSE: 0.777161
Time (s): 0.0277111

SequentialSfMReconstructionEngine::ComputeResidualsMSE.
-- #Tracks: 210
-- Residual min: 0.000116665
-- Residual median: 0.223267
-- Residual max: 5.90642
-- Residual mean: 0.691327

=========================
MSE Residual InitialPair Inlier: 0.691327

-- Structure from Motion (statistics):
-- #Camera calibrated: 2 from 10 input images.
-- #Tracks, #3D points: 200

SequentialSfMReconstructionEngine::ComputeResidualsMSE.
-- #Tracks: 200
-- Residual min: 0.000116665
-- Residual median: 0.20109
-- Residual max: 3.77725
-- Residual mean: 0.60566

Histogram of residuals:

0 | 487
0.378 | 72
0.755 | 71
1.13 | 48
1.51 | 55
1.89 | 24
2.27 | 16
2.64 | 10
3.02 | 9
3.4 | 7
3.78

Compute scene structure color
0% 10 20 30 40 50 60 70 80 90 100%
|----|----|----|----|----|----|----|----|----|----|

0% 10 20 30 40 50 60 70 80 90 100%
|----|----|----|----|----|----|----|----|----|----|

Reading bundle...2 cameras -- 200 points in bundle file

2 cameras -- 200 points
Reading images: *Divide:
*
Set widths/heights...done 0 secs
done 0 secs
slimNeighborsSetLinks...done 0 secs
mergeSFM...***********resetPoints...done
Rep counts: 200 -> 35 0 secs
setScoreThresholds...done 0 secs
sRemoveImages... **
Kept: 0 1

Removed:
sRemoveImages: 2 -> 2 0 secs
slimNeighborsSetLinks...done 0 secs

Cluster sizes:
2
Adding images:
0
Image nums: 2 -> 2 -> 2
done 0 secs
1 images in vis on the average

g:\Program Files\Regard3D\pmvs\pmvs2.exe
pictureset_1\matching_1\triangulation_0\densification_0\PMVS/
option-0000

--- Summary of specified options ---

of timages: 2 (enumeration)

of oimages: 0 (enumeration)

level: 1 csize: 2
threshold: 0.7 wsize: 7
minImageNum: 3 CPU: 8
useVisData: 1 sequence: -1

Reading images: **
0 1 Harris running ...Harris running ...6381 harris done
DoG running...6350 harris done
DoG running...7715 dog done
7737 dog done
done
adding seeds
(0,0)(1,0)done
---- Initial: 0 secs ----
Total pass fail0 fail1 refinepatch: 27212 0 27212 0 0
Total pass fail0 fail1 refinepatch: 100 0 100 0 0
Expanding patches...
---- EXPANSION: 0 secs ----
Total pass fail0 fail1 refinepatch: 0 0 0 0 0
Total pass fail0 fail1 refinepatch: -nan(ind) -nan(ind) -nan(ind) -nan(ind) -nan(ind)
FilterOutside
mainbody:
Gain (ave/var): 0 0
0 -> 0 (-nan(ind)%) 0 secs
Filter Exact: **
0 -> 0 (-nan(ind)%) 0 secs
FilterNeighbor: FilterGroups: STATUS: 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0
Expanding patches...
---- EXPANSION: 0 secs ----
Total pass fail0 fail1 refinepatch: 0 0 0 0 0
Total pass fail0 fail1 refinepatch: -nan(ind) -nan(ind) -nan(ind) -nan(ind) -nan(ind)
FilterOutside
mainbody:
Gain (ave/var): 0 0
0 -> 0 (-nan(ind)%) 0 secs
Filter Exact: **
0 -> 0 (-nan(ind)%) 0 secs
FilterNeighbor: FilterGroups: STATUS: 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0
Expanding patches...
---- EXPANSION: 0 secs ----
Total pass fail0 fail1 refinepatch: 0 0 0 0 0
Total pass fail0 fail1 refinepatch: -nan(ind) -nan(ind) -nan(ind) -nan(ind) -nan(ind)
FilterOutside
mainbody:
Gain (ave/var): 0 0
0 -> 0 (-nan(ind)%) 0 secs
Filter Exact: **
0 -> 0 (-nan(ind)%) 0 secs
FilterNeighbor: FilterGroups: STATUS: 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0
---- Total: 0 secs ----
Running Screened Poisson Reconstruction (Version 9.01)
--in pictureset_1\matching_1\triangulation_0\densification_0\PMVS\models\option-0000.ply
--depth 9
--out pictureset_1\matching_1\triangulation_0\densification_0\surface_0/model_surface.ply
--verbose
--samplesPerNode 1.000000
--pointWeight 4.000000
--density
Input Points / Samples: 0 / 0

Read input into tree: 0.0 (s), 3.2 (MB) / 3.2 (MB) / 3.2 (MB)

Got kernel density: 0.0 (s), 3.2 (MB) / 3.2 (MB) / 3.2 (MB)

Got normal field: 0.0 (s), 3.2 (MB) / 3.2 (MB) / 3.2 (MB)

Finalized tree: 0.0 (s), 3.4 (MB) / 3.4 (MB) / 3.4 (MB)

Set FEM constraints: 0.0 (s), 3.6 (MB) / 3.6 (MB) / 3.6 (MB)

#Set point constraints: 0.0 (s), 3.6 (MB) / 3.6 (MB) / 3.6 (MB)
Leaf Nodes / Active Nodes / Ghost Nodes: 64 / 72 / 1
Memory Usage: 3.598 MB
Depth[0/0]: Evaluated / Got / Solved in: 0.000 / 0.000 / 0.000 (3.789 MB) Nodes: 1

Linear system solved: 0.0 (s), 3.8 (MB) / 3.8 (MB) / 3.8 (MB)

Got average: 0.0 (s), 4.0 (MB) / 4.0 (MB) / 4.0 (MB)
Iso-Value: -nan(ind)
Vertices / Polygons: 0 / 0

Got triangles: 0.0 (s), 4.1 (MB) / 4.1 (MB) / 4.1 (MB)

Total Solve: 0.0 (s), 4.1 (MB)

Using the OPENMP thread interface

• Features Loading -
  0% 10 20 30 40 50 60 70 80 90 100%
  |----|----|----|----|----|----|----|----|----|----|

• Features Loading -
  0% 10 20 30 40 50 60 70 80 90 100%
  |----|----|----|----|----|----|----|----|----|----|

Track building

Track filtering

Track export to internal struct

Track stats

-- Tracks Stats --
Tracks number: 449
Images Id:
3, 4, 5, 6, 7,

TrackLength, Occurrence
2 433
3 16

A-Contrario initial pair residual: 77.0975

Bundle Adjustment statistics (approximated RMSE):
#views: 2
#poses: 2
#intrinsics: 1
#tracks: 350
#residuals: 1400
Initial RMSE: 0.906527
Final RMSE: 0.777161
Time (s): 0.0579963

SequentialSfMReconstructionEngine::ComputeResidualsMSE.
-- #Tracks: 210
-- Residual min: 0.000116665
-- Residual median: 0.223267
-- Residual max: 5.90642
-- Residual mean: 0.691327

=========================
MSE Residual InitialPair Inlier: 0.691327

-- Structure from Motion (statistics):
-- #Camera calibrated: 2 from 10 input images.
-- #Tracks, #3D points: 200

SequentialSfMReconstructionEngine::ComputeResidualsMSE.
-- #Tracks: 200
-- Residual min: 0.000116665
-- Residual median: 0.20109
-- Residual max: 3.77725
-- Residual mean: 0.60566

Histogram of residuals:

0 | 487
0.378 | 72
0.755 | 71
1.13 | 48
1.51 | 55
1.89 | 24
2.27 | 16
2.64 | 10
3.02 | 9
3.4 | 7
3.78

Compute scene structure color
0% 10 20 30 40 50 60 70 80 90 100%
|----|----|----|----|----|----|----|----|----|----|

0% 10 20 30 40 50 60 70 80 90 100%
|----|----|----|----|----|----|----|----|----|----|

Reading bundle...2 cameras -- 200 points in bundle file

2 cameras -- 200 points
Reading images: *Divide:
*
Set widths/heights...done 0 secs
done 0 secs
slimNeighborsSetLinks...done 0 secs
mergeSFM...***********resetPoints...done
Rep counts: 200 -> 35 0 secs
setScoreThresholds...done 0 secs
sRemoveImages... **
Kept: 0 1

Removed:
sRemoveImages: 2 -> 2 0 secs
slimNeighborsSetLinks...done 0 secs

Cluster sizes:
2
Adding images:
0
Image nums: 2 -> 2 -> 2
done 0 secs
1 images in vis on the average

g:\Program Files\Regard3D\pmvs\pmvs2.exe
pictureset_1\matching_2\triangulation_0\densification_0\PMVS/
option-0000

--- Summary of specified options ---

of timages: 2 (enumeration)

of oimages: 0 (enumeration)

level: 1 csize: 2
threshold: 0.7 wsize: 7
minImageNum: 3 CPU: 8
useVisData: 1 sequence: -1

Reading images: **
0 1 Harris running ...Harris running ...6350 harris done
6381 harris done
DoG running...DoG running...7715 dog done
7737 dog done
done
adding seeds
(0,0)(1,0)done
---- Initial: 0 secs ----
Total pass fail0 fail1 refinepatch: 27280 0 27280 0 0
Total pass fail0 fail1 refinepatch: 100 0 100 0 0
Expanding patches...
---- EXPANSION: 0 secs ----
Total pass fail0 fail1 refinepatch: 0 0 0 0 0
Total pass fail0 fail1 refinepatch: -nan(ind) -nan(ind) -nan(ind) -nan(ind) -nan(ind)
FilterOutside
mainbody:
Gain (ave/var): 0 0
0 -> 0 (-nan(ind)%) 0 secs
Filter Exact: **
0 -> 0 (-nan(ind)%) 0 secs
FilterNeighbor: FilterGroups: STATUS: 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0
Expanding patches...
---- EXPANSION: 0 secs ----
Total pass fail0 fail1 refinepatch: 0 0 0 0 0
Total pass fail0 fail1 refinepatch: -nan(ind) -nan(ind) -nan(ind) -nan(ind) -nan(ind)
FilterOutside
mainbody:
Gain (ave/var): 0 0
0 -> 0 (-nan(ind)%) 0 secs
Filter Exact: **
0 -> 0 (-nan(ind)%) 0 secs
FilterNeighbor: FilterGroups: STATUS: 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0
Expanding patches...
---- EXPANSION: 0 secs ----
Total pass fail0 fail1 refinepatch: 0 0 0 0 0
Total pass fail0 fail1 refinepatch: -nan(ind) -nan(ind) -nan(ind) -nan(ind) -nan(ind)
FilterOutside
mainbody:
Gain (ave/var): 0 0
0 -> 0 (-nan(ind)%) 0 secs
Filter Exact: **
0 -> 0 (-nan(ind)%) 0 secs
FilterNeighbor: FilterGroups: STATUS: 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0
---- Total: 0 secs ----
Running Screened Poisson Reconstruction (Version 9.01)
--in pictureset_1\matching_2\triangulation_0\densification_0\PMVS\models\option-0000.ply
--depth 9
--out pictureset_1\matching_2\triangulation_0\densification_0\surface_0/model_surface.ply
--verbose
--samplesPerNode 1.000000
--pointWeight 4.000000
--density
Input Points / Samples: 0 / 0

Read input into tree: 0.0 (s), 3.2 (MB) / 3.2 (MB) / 3.2 (MB)

Got kernel density: 0.0 (s), 3.2 (MB) / 3.2 (MB) / 3.2 (MB)

Got normal field: 0.0 (s), 3.2 (MB) / 3.2 (MB) / 3.2 (MB)

Finalized tree: 0.0 (s), 3.4 (MB) / 3.4 (MB) / 3.4 (MB)

Set FEM constraints: 0.0 (s), 3.6 (MB) / 3.6 (MB) / 3.6 (MB)

#Set point constraints: 0.0 (s), 3.6 (MB) / 3.6 (MB) / 3.6 (MB)
Leaf Nodes / Active Nodes / Ghost Nodes: 64 / 72 / 1
Memory Usage: 3.598 MB
Depth[0/0]: Evaluated / Got / Solved in: 0.000 / 0.000 / 0.001 (3.789 MB) Nodes: 1

Linear system solved: 0.0 (s), 3.8 (MB) / 3.8 (MB) / 3.8 (MB)

Got average: 0.0 (s), 4.0 (MB) / 4.0 (MB) / 4.0 (MB)
Iso-Value: -nan(ind)
Vertices / Polygons: 0 / 0

Got triangles: 0.0 (s), 4.1 (MB) / 4.1 (MB) / 4.1 (MB)

Total Solve: 0.0 (s), 4.1 (MB)

Thanks,
Paul.

Hi. I have two questions.

1. In Poisson Surface Reconstruction (2006), the value of the base function is sacled by (1 / o.w^3).
   
   Do you implement this feature in screened Poisson, i.e. is the B-spline function scaled by (1 / o.w^3) in 3D or (1 / o.w) in 1D whenever it is called ?
2. Is the implicit function value of an arbitary point computed as the linear combination of the solution of related nodes ( x(p) = sum(xiBi(p) ) ? Except for B-spline function, there is no other scaling factor to adjust the weight of solution value of different depth nodes, right ?

sorry to disturb, but any clue on how to achieve this?