std::bad_alloc with vcgIsolated

[stla]

Hello,

I have a very large mesh3d, with ~60 millions triangles. When I run vcgIsolated(mymesh, diameter=0) on this mesh, an error is generated: std::bad_alloc. Do you know a workaround, or an alternative ? (I took a look at RvtkStatismo, but I saw nothing about the connected components).

[zarquon42b]

If RAM is sufficient, I do not know. Can you send me the mesh?

[stla]

Hello,

Thank you for your prompt reply. If I find my Dropbox credentials, I can send you a link to a rds file. I can also post here the code which generates this mesh (it requires a Rcpp file). What do you prefer ?

[zarquon42b]

Code is fine. Did you check your RAM?

[stla]

Sorry for the delay, my PC has crashed. I have 8Gb of RAM and perhaps this is not enough. Even vcgClean generates a std::bad_alloc. I'm preparing the code and I post it within 5 minutes.

[zarquon42b]

I just tested it by subdividing an existing mesh. And with 12Mio Vertices, vcgIsolated takes up 7.7GB of RAM. So this will definitively be the bottleneck here.

[stla]

Ok, so there's no way ? I've noticed that "Split in connected components" in Meshlab is highly slower than Rvcg::vcgIsolated(mymesh, diameter = 0).

Here is the code if you want to try. The cpp file first:

// file mandelbulb.cpp
#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
double mandelbulb0(NumericVector x0y0z0){
  double x0 = x0y0z0[0];
  double y0 = x0y0z0[1];
  double z0 = x0y0z0[2];
  double x = x0;
  double y = y0;
  double z = z0;
  double r2, theta, phi, r8;
  for(int i = 0; i < 24; i++){
    r2 = x*x + y*y + z*z;
    if(r2 > 4){
      break;
    }
    theta = 8*atan2(sqrt(x*x+y*y),z);
    phi = 8*atan2(y,x);
    r8 = r2*r2*r2*r2;
    x = r8*cos(phi)*sin(theta) + x0;
    y = r8*sin(phi)*sin(theta) + y0;
    z = r8*cos(theta) + z0;
  }
  return sqrt(r2);
}

// [[Rcpp::export]]
NumericVector mandelbulb(double m, double M, unsigned n) {
  NumericVector out(n*n*n);
  double h = (M-m)/(n-1);
  NumericVector xyz(3);
  double x,y,z;
  unsigned l = 0;
  for(unsigned i=0; i<n; i++){
    xyz[0] = m + i*h;
    for(unsigned j=0; j<n; j++){
      xyz[1] = m + j*h;
      for(unsigned k=0; k<n; k++){
        xyz[2] = m + k*h;
        if(x*x+y*y+z*z > 4){
          out[l] = R_PosInf;
        }else{
          out[l] = mandelbulb0(xyz);
        }
        l += 1;
      }
    }
  }
  return out;
}

// test
/*** R
mandelbulb0(c(1,1,1))
mandelbulb(-1.1, 1.1, 3L)
*/

Now the R code:

# required packages : Rvcg, RvtkStatismo, neuroim, data.table

# set working directory
setwd("......")

# sourve the CPP file
Rcpp::sourceCpp("mandelbulb.cpp")

# make voxel ####
n <- 1024L
v <- mandelbulb(-1.2, 1.2, n)

# write nifti file ####
# I can't use neuroim::writeVolume because My RAM is not sufficient.
# Here is a workaround.
writeNifti <- function(voxel, n, fileName, dataType = NULL){
  vol <- neuroim::BrainVolume(numeric(2*2*2),
                               neuroim::BrainSpace(c(2,2,2), c(1,1,1)))
  vol@source@metaInfo@Dim <- vol@space@Dim <- c(n,n,n)
  hdr <- neuroim:::as.nifti.header(vol, vol@source@metaInfo)
  if (!is.null(dataType) && dataType != vol@source@metaInfo@dataType) {
    hdr$datatype <- neuroim:::.getDataCode(dataType)
    hdr$dataStorage <- neuroim:::.getDataStorage(hdr$datatype)
    hdr$bitpix <- neuroim:::.getDataSize(dataType) * 8
  }
  else {
    dataType <- vol@source@metaInfo@dataType
    hdr$datatype <- neuroim:::.getDataCode(dataType)
    hdr$dataStorage <- neuroim:::.getDataStorage(hdr$datatype)
    hdr$bitpix <- neuroim:::.getDataSize(dataType) * 8
  }
  conn <- if (substr(fileName, nchar(fileName) - 2, nchar(fileName)) == 
              ".gz") {
    gzfile(fileName, open = "wb")
  }
  else {
    file(fileName, open = "wb")
  }
  neuroim:::writeNIfTIHeader(hdr, conn, close = FALSE)
  writer <- neuroim::BinaryWriter(conn, hdr$voxOffset, dataType, hdr$bitpix/8, 
                         .Platform$endian)
  for(i in 1:ncol(voxel)){
    print(i)
    neuroim::writeElements(writer, voxel[,i])
  }
  neuroim::close(writer)
}

data.table::setattr(v, "dim", c(n*n*n/4L,4L))
writeNifti(v, n, "mandelbulb1024.nii.gz")

# v is not used anymore
rm("v")

# make isosurface ####
mesh <- RvtkStatismo::vtkTriangulate("mandelbulb1024.nii.gz", 2)

# make normals ####
mesh <- Rvcg::vcgUpdateNormals(mesh, pointcloud = c(ncol(mesh$vb),0))

# faces are inverted ####
mesh <- Morpho::invertFaces(mesh)

# keep the big connected component only ####
mesh <- Rvcg::vcgIsolated(mesh, diameter = 0)
mesh$remvert <- NULL

With n <- 512L, everything works fine. After getting the mesh I save it to a PLY file and I do "Surface reconstruction: Poisson screening" in Meshlab, then I isolate the biggest connected component, and I get a nice Mandelbulb. I've also managed to get it with n <- 768L.

[stla]

Nice, I've managed to isolate the big connected component with Meshlab (Remove isolated pieces, facenum: 100000). :+1: