Tetrahedral meshing from micro-CT image stack

[JIAZHENG HU]

Dear Prof. Fang,

       First of all, Happy Lunar New Year to you and your crew! I appreciate your help on "hacking" from 3D to 2D meshing (my previous post "From 2D image to 2D triangular mesh"). I have been playing around with it recently, and it works out just great!

       Recently, I moved to the next step, which is volumetric meshing using micro-CT image stack (.tif). The script is printed below (Direct_volmesh.m in the shared folder). 

%------------Image stack 2 volumetric mesh---------------------

% Direct meshing from image stack 

fprintf(1,'loading segmented concrete image...\n');

for i=1:700                                    

  concrete(:,:,i)=imread('UCS_4.tif',i);

end

concrete=uint8(concrete);

% "UCS_4.tif" (shared in google drive) is the original 16bit stack image, it has not been segmented

% "UCS_4_segmented.tif" (shared in google drive) is pre-segmented using my script "color_change.m" (also in the shared folder)

% check a random slice for image properties

single = concrete(:,:,172);

% from previous work (color_change.m)---> UCS_4_segmented.tif havs segmented CT image where:

% pixel value = 1 ---> pores in concrete

% pixel value = 2 ---> cement paste in concrete

% pixel value = 3 ---> aggregate in concrete

% call cgalmesher to mesh the segmented volume

fprintf(1,'meshing the segmented concrete ...\n');

[node,elem,face]=v2m(concrete,[],10,20,'cgalmesh');

% plot

figure

Crt=plottetra(node(:,1:3),elem);        

axis equal;

%----------------END-------------------

        When I used the original image stack, called "UCS_4.tif", the meshing function worked out well. however, it seemed that the algorithm was having little trouble segmenting the 16bit image, because the 5th column of the output variable "elem"  had only 1 value. Therefore, I wrote another script (color_change.m in the shared folder) to apply global thresholdings to segment all of the CT images into 3 channels. Then, I used ImageJ to collect individual images into one tiff file.

        Unfortunately, problems came across when I threw the pre-segmented tiff file (UCS_4_semented.tif) into the above script. I also had a hard time interpreting the error message (text file in the shared folder). I am using WINx64, R2018b.

        I put all the relevant files in a shared google folder (those tiff files are kinda large), and the link is printed underneath:

        https://drive.google.com/open?id=1SJEe1UyET5h_fQELGleMfINXr30rlBi-

Sincerely,

Jiazheng

  

[Qianqian Fang]

On 1/28/19 3:11 PM, JIAZHENG HU wrote:

Dear Prof. Fang,

       First of all, Happy Lunar New Year to you and your crew! I appreciate your help on "hacking" from 3D to 2D meshing (my previous post "From 2D image to 2D triangular mesh"). I have been playing around with it recently, and it works out just great!

       Recently, I moved to the next step, which is volumetric meshing using micro-CT image stack (.tif). The script is printed below (Direct_volmesh.m in the shared folder). 

...

% call cgalmesher to mesh the segmented volume

fprintf(1,'meshing the segmented concrete ...\n');

[node,elem,face]=v2m(concrete,[],10,20,'cgalmesh');

...

        When I used the original image stack, called "UCS_4.tif", the meshing function worked out well. however, it seemed that the algorithm was having little trouble segmenting the 16bit image, because the 5th column of the output variable "elem"  had only 1 value. Therefore, I wrote another script (color_change.m in the shared folder) to apply global thresholdings to segment all of the CT images into 3 channels. Then, I used ImageJ to collect individual images into one tiff file.

        Unfortunately, problems came across when I threw the pre-segmented tiff file (UCS_4_semented.tif) into the above script. I also had a hard time interpreting the error message (text file in the shared folder). I am using WINx64, R2018b.

hi Jiazheng

the default release of iso2mesh only included 32bit Windows binaries for cgalmesh.exe. As a result, it can not handle large volumes - anything requires 2GB or more memory to process will result in a crash.

you can either downsample your input volume and try again, or download the 64bit cgal*.exe binaries in the win64 branch (provided by another contributor)

https://github.com/fangq/iso2mesh/tree/win64

Qianqian

        I put all the relevant files in a shared google folder (those tiff files are kinda large), and the link is printed underneath:

        https://drive.google.com/open?id=1SJEe1UyET5h_fQELGleMfINXr30rlBi-

Sincerely,

Jiazheng

  

--
You received this message because you are subscribed to the Google Groups "iso2mesh-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to iso2mesh-user...@googlegroups.com.
To post to this group, send email to iso2mes...@googlegroups.com.
Visit this group at https://groups.google.com/group/iso2mesh-users.
For more options, visit https://groups.google.com/d/optout.

[JIAZHENG HU]

Thanks Qianqian, I have tried my script on volume reduced .tif, it works out:) thanks for you help! 

To unsubscribe from this group and stop receiving emails from it, send an email to iso2mesh-use...@googlegroups.com.
To post to this group, send email to iso2me...@googlegroups.com.

[iso2meshNewbie]

Dear Prof. Fang,

I am getting error in meshing a simple binary tiff stack (segmented to 0s and 1s):

% clear memory

clear variables;

% change the current folder 

cd 'C:\Users\Stack'

%------------Image stack 2 volumetric mesh---------------------

% Direct meshing from image stack 

fprintf(1,'loading segmented image stack...\n');

for i=1:2988

  mydata(:,:,i)=imread(sprintf('strut_35.26_2_tomo_%04d.tif', i));

end

mydata=uint8(mydata);

fprintf(1,'loading completed...\n');

% mesh the segmented volume

fprintf(1,'meshing the segmented image stack...\n');

[node,elem,face]=v2m(mydata,0.5,5,100);

% plot

figure

plotmesh(node,face);

%----------------END-------------------

Error message:

Error using fgets

Invalid file identifier. Use fopen to generate a valid file identifier.

Error in fgetl (line 32)

[tline,lt] = fgets(fid);

Error in readoff (line 23)

line=fgetl(fid);

Error in vol2restrictedtri (line 54)

[node,elem]=readoff(mwpath('post_extract.off'));

Error in vol2surf (line 172)

          [v0,f0]=vol2restrictedtri(newimg,isovalues(i)-perturb,regions(i,:),...

Error in vol2mesh (line 63)

[no,el,regions,holes]=vol2surf(img,ix,iy,iz,opt,dofix,method,isovalues);

Error in v2m (line 18)

[node,elem,face]=vol2mesh(img,1:size(img,1),1:size(img,2),1:size(img,3),opt,maxvol,1,method,isovalues);

Error in Strut_mesh (line 18)

[node,elem,face]=v2m(mydata,0.5,5,100);

Any solution to overcome this error?

To unsubscribe from this group and stop receiving emails from it, send an email to iso2me...@googlegroups.com.
To post to this group, send email to iso2me...@googlegroups.com.

[Qianqian Fang]

On 11/25/19 3:18 AM, iso2meshNewbie wrote:

Dear Prof. Fang,

I am getting error in meshing a simple binary tiff stack (segmented to 0s and 1s):

I don't see a problem with the code, can you save your mydata array into a .mat file and send it to me offline (with a dropbox, google drive share link)?

To unsubscribe from this group and stop receiving emails from it, send an email to iso2mesh-user...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/iso2mesh-users/a8f4ea13-0dae-41d5-969d-1a41ae4ad57b%40googlegroups.com.

[Qianqian Fang]

sorry for the delay in getting back to you.

I just tested your dataset, it is just too big for most computer to handle. it crashed my desktop with 32GB memory by throwing an out of memory error.

dowsampling the volume by 2 in each dimension (i.e. a 8x reduction in volume size) made it a lot easier to mesh. see the screenshot

Qianqian

To view this discussion on the web visit https://groups.google.com/d/msgid/iso2mesh-users/51f56ad0-e208-6f69-c5d7-abd8b3be1d34%40neu.edu.

[iso2meshNewbie]

Thank you for the solution, can you please share the script here?

> I just tested your dataset, it is just too big for most computer to handle. it crashed my desktop with 32GB memory by throwing an out of memory error.

My computer has 256GB memory but it still failed.

> dowsampling the volume by 2 in each dimension (i.e. a 8x reduction in volume size) made it a lot easier to mesh. see the screenshot

Does it mean the volume is reduced or the resolution is downsampled?

[Qianqian Fang]

On 12/4/19 11:31 AM, iso2meshNewbie wrote:

Thank you for the solution, can you please share the script here?

it was highlighted in the screenshot, just one line

[no,el]=v2m(mydata(1:2:end,1:2:end,1:2:end), 0.5, 10, 100);

> I just tested your dataset, it is just too big for most computer to handle. it crashed my desktop with 32GB memory by throwing an out of memory error.

My computer has 256GB memory but it still failed.

I just tried it on one of my servers with 128GB memory, and it also crashed even though the peak memory used was about 44GB. so, I believe there must be a memory problem with the default surface extraction tool (cgalsurf). I admit that the CGAL executables we provided in the iso2mesh/bin folder were compiled using a very old CGAL library (v3.8), so such memory issue may have fixed in the later releases.

to use the newer versions of CGAL, you will have to recompile these tools on your computer using the provided source codes in iso2mesh/tools.

you need to install a few dependencies, including cmake, cgal and superlu in order to recompile, see this packaging script I made for Fedora recently

https://src.fedoraproject.org/rpms/octave-iso2mesh/blob/master/f/octave-iso2mesh.spec#_28

> dowsampling the volume by 2 in each dimension (i.e. a 8x reduction in volume size) made it a lot easier to mesh. see the screenshot

Does it mean the volume is reduced or the resolution is downsampled?

both. the volume is downsampled, and the size of the raw input volume is reduced.

I often found a common misunderstanding on what is the "desired" mesh density one should be using (and related, what should be the desired volume resolution as the input).

many researchers appear to be keen on producing high density meshes from high density volumes. this is fine as long as the memory and time needed are not excessive - but I want to say is that this is not necessarily always needed.

One must be aware that the process of creating a mesh from a volume itself is a downsampling process - you want to use a small number of nodes to accurately capture the shapes, but you don't want to go to the extreme - where your mesh/surface density is close to voxel density. In such case, the stair-shaped voxel boundary will appear in your mesh, and such feature is not real but an artifact due to discretization. Therefore, you don't want to set the opt.radbound parameter to be near 1 when you call v2m/v2s, instead, it should be something 5-10 or even larger so that the stair-shaped artifacts can be smoothen out. As a result, your mesh is already performing a downsampling by 5x or more.

another big downside of producing an extensively dense mesh is to place a huge burden to your post processing such as FEM modeling using the generated mesh.

for example, if you plan to set your opt.radbound to 10 (voxels), a faster way to do this is to downsample your volume by a factor of 2x, and then set opt.radbound to 5. This will produce mesh of similar density or quality, but the input volume is much smaller, thus much faster to process.

let me know if this makes sense to you.

Qianqian

--
You received this message because you are subscribed to the Google Groups "iso2mesh-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to iso2mesh-user...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/iso2mesh-users/8b54e6dc-f70e-4766-aba9-612e66af1b72%40googlegroups.com.

[iso2meshNewbie]

I have ran your code and still having the same error, wondering what could wrong?

To unsubscribe from this group and stop receiving emails from it, send an email to iso2me...@googlegroups.com.

[Qianqian Fang]

if you are able to get it to work by further downsampling the volume, then the other possibility is that you are using a windows machine.

if so, try downloading the win64 branch of iso2mesh, and use the exe files in the bin folder.

the default iso2mesh packages include only 32bit windows binaries, thus is limited in maximum memory size

To unsubscribe from this group and stop receiving emails from it, send an email to iso2mesh-user...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/iso2mesh-users/bd32188f-15da-42db-906a-7d68f8a5a8e2%40googlegroups.com.

[iso2meshNewbie]

I was able to mesh it and export to ABAQUS, but I have some issues:

if you look into the mat file sent earlier, there are small voids inside the material but the meshing procedure is unable to capture the fine details of voids. How should I approach in retaining the finer voids (smallest element size near voids, but relatively larger element size away from voids)?