Hello, my name is JS Gauthier. I am a artist working with scientific datasets in Saskatoon, Canada (https://vimeo.com/jsgauthier or https://jsgauthier.com/ ). I have created a series of 3D scans using X-Ray micro-tomography of several species of embryos at different developmental time-points (in collaboration with Dr. Brian Eames). I’ve been seeking ways to use Processing to visualize my 3D data to create interactive volumetric representations of these image stacks.
My hope is that this thread might aid myself and others to adopt CT scans as an artistic media using Processing rather than costly industrial solutions.
I have managed to create a basic implementation to show a low resolution instances of my data. I haven’t found a thread that discusses my particular interest in the Forum and would appreciate help or dialogue on how to advance my coding project.
This is my first post in the forum and my first time sharing code, so go easy on my folks. I have no formal training, I’ve just followed the Processing language pages and also read D. Shiffman’s “Learning Processing” (which I highly recommend).
I believe you should look at the Toxiclibs library, this sketch might be very relevant:-
/**
* MRISurface demo showing how to utilize the IsoSurface class to efficiently
* visualise volumetric data, in this case using MRI scan data obtained from the
* volvis.org repository of free datasets.
* The demo also shows how to save the generated mesh as binary STL file (or
* alternatively in OBJ format) for later use in other 3D tools/digital fabrication.
*
* I've planned further classes for the toxi.geom.volume package to easier draw
* and manipulate volumetric data.
*
* Important note: This demo is fairly memory hungry, so if you receive OutOfMemory
* errors, please make sure to give Processing as much RAM as possible (see preferences)
* and if this still fails, reduce the surface resolution from DIM=128 to 64 or 32...
*
* Controls:
* Click mouse button to toggle rendering style between shaded/wireframe.
*/
/*
* Copyright (c) 2009 Karsten Schmidt
*
* This demo & library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
*
* http://creativecommons.org/licenses/LGPL/2.1/
*
* This library is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with this library; if not, write to the Free Software
* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
*/
import toxi.geom.*;
import toxi.geom.mesh.*;
import toxi.volume.*;
import toxi.math.noise.*;
import toxi.processing.*;
int DIM=128;
float ISO_THRESHOLD = 0.1;
Vec3D SCALE=new Vec3D(DIM,DIM,DIM).scaleSelf(8);
IsoSurface surface;
Mesh3D mesh;
ToxiclibsSupport gfx;
boolean isWireframe=false;
void setup() {
size(1024,768,P3D);
gfx=new ToxiclibsSupport(this);
strokeWeight(0.5);
// convert MRI scan data into floats
// MRI data is 256 x 256 x 256 voxels @ 8bit/voxel
byte[] mriData=loadBytes("aneurism.raw.gz");
// scale factor to normalize 8bit to the 0.0 - 1.0 interval
float mriNormalize=1/255.0;
// setup lower resolution grid for IsoSurface
VolumetricSpaceArray volume=new VolumetricSpaceArray(SCALE,DIM,DIM,DIM);
float[] cloud=volume.getData();
int stride=256/DIM;
for(int z=0,idx=0; z<256; z+=stride) {
for(int y=0; y<256; y+=stride) {
int sliceIdx=y*256+z*65536;
for(int x=0; x<256; x+=stride) {
byte b=mriData[x+sliceIdx];
cloud[idx++]=(int)(b<0 ? 256+b : b)*mriNormalize;
}
}
}
long t0=System.nanoTime();
// create IsoSurface and compute surface mesh for the given iso threshold value
surface=new HashIsoSurface(volume,0.15);
mesh=surface.computeSurfaceMesh(null,ISO_THRESHOLD);
float timeTaken=(System.nanoTime()-t0)*1e-6;
println(timeTaken+"ms to compute "+mesh.getNumFaces()+" faces");
}
void draw() {
background(128);
translate(width/2,height/2,0);
rotateX(mouseY*0.01);
rotateY(mouseX*0.01);
ambientLight(48,48,48);
lightSpecular(230,230,230);
directionalLight(255,255,255,0,-0.5,-1);
specular(255,255,255);
shininess(16.0);
if (isWireframe) {
stroke(255);
noFill();
}
else {
noStroke();
fill(255);
}
gfx.mesh(mesh);
}
void normal(Vec3D v) {
normal(v.x,v.y,v.z);
}
void vertex(Vec3D v) {
vertex(v.x,v.y,v.z);
}
void mousePressed() {
isWireframe=!isWireframe;
}
Unfortunately toxi has moved on from processing, otherwise I expect he might have developed the idea a bit further.
@JSGauthier
Thank you for bringing this topic to the forum.
I would like to work with and understand micro-CT data.
If anyone with experience in this area can direct me to resources I would appreciate it.
Is there a good database for micro-CT scans to work with?
References:
:)
Side note (not intended to distract from main topic or questions above):
I had some fun modifying this code and generating 3D data and visualizations with this a while back:
It would be of interest to me to work other data such as CT data.
Thanks glv,
I appreciate the links! I am aware of several of these links and glad to discover some new ones. I’ve been working for a bit over a year (on my free time) to learn coding and Coding train is my one of my youtube favorites! I am hoping to implement some vector know how from Nature of Code into my projects in the near future.
–> PLEASE NOTE You must download the attached data (see dropbox link below) for it to run this as intended.
**My main problem is that I cannot visualize both sides of the image stack.
Could this be due to PImage() having face culling?
I implemented a PShape array and applied my PGraphics layers as textures but I still have disappearing back faces. When I rotate them the will disappear.**
One less than perfect solution was to use Hint(DISABLE_DEPTH_TEST); but this causes unusual distortions when the volume is zoomed in.
Is anyone aware of a way to disable this face culling for PShape() or image() objects?
I would like to render both front and back faces so the volume does not disappear, but can’t seem to find a way to solve this. I made a way to create boxes in the place of the pixels. (This function is active in the code I am posting as makePixBoxes()
The Dataset (a stack of 572 .png files) can be downloaded from this dropbox link and should placed the sketches \data folder.
/**
Volumetric_CT_Renderer_v001
by J-S Gauthier, 2020
Dataset can be downloaded from this dropbox link and placed in the data folder.
https://www.dropbox.com/sh/8s6r0mxmy1goss3/AADrfZ4PqVOx0AL1X49Ek01la?dl=0
CONTROLS:
Use middle mouse wheel to zoom in/out & middle mouse click to rotate the datasets.
Requirements:
-->DONE load image n of stack all Images from data folder
-->DONE load pixels
-->DONE Evaluate brightness of each pixel in pixel array
-->DONE Assign alpha 0 to each pixel above chosen brightness threshold
-->DONE Make a box for each pixel
-->DONE Display all 3D boxes
-->DONE Replace each image object with a PShape object and assign a texture to the PShape
Priority items:
Make PShapes or Image objects visible have both sides visible. (no back face culling).
Later implementations:
Make groups for each colour value
Join some of the groups
**/
PShape[] slices; //create Pshape of 4 vertices to upon which to project PImage layers as textures.
PImage[] layers; //create PImage array to contain the CT scan series in memory.
PGraphics[] altLayers; //create Pgraphics array to contain altered RGBA pixel data from original CT data as required (such as alpha transparecnies).
String path = sketchPath() + "\\data\\";
// Eventually I want to change String "\\data\\" into other String values to alter folder locations on a "per object" basis in order to show multiple species.
int lCount; // Layer count corresponds to # of files in the stack.
float zDepth; // Zdepth parameter corresponds to the Zspacing between each stacked layer. This should be determined by the # pf pixels n W & H of the data and their scaler value.
float rotY; // Stack rotation param for Y axis.
float rotX; // Stack rotation param for X (hold middleMouse button to alter this value).
int cellsize = 5; // Dimensions of each box/cell in the grid, for rendering 3D boxes as voxels.
int columns, rows; // Cols and rows * cellsize to define the # of boxes to generate.
void setup() {
size(500, 500, P3D);
surface.setResizable(true); //enables changes to the window size
String[] filenames = listFileNames(path);
//printArray(filenames);
lCount = (filenames.length-1); //this make the number of layers to the ammount of files in the data folder
layers = new PImage [lCount];
altLayers = new PGraphics[lCount];
slices = new PShape[lCount];
createSlices(lCount);
loadLayers(lCount);
makeAltLayers(lCount);
loadLayersPix();
makeSlices(); //this make Pshape rather than image objects and applies the dataset as textures. I was hoping PShapes would not have face culling.
}
void draw() {
ambientLight(255, 255, 255); // lights set to max white
//hint(DISABLE_DEPTH_TEST); //this allows on to see the back of the Pshape or images but also produces distortion in the output (not ideal).
imageMode(CENTER);
background(0);
translate(width/2, height/2, 0);
push();
if (mousePressed == true && mouseButton == CENTER) {
rotY = map(mouseX, 0, width, 0, 360);
rotX = map(mouseY, 0, height, 0, 360);
}
translate(0,0,zDepth);
rotateX(radians(rotX));
rotateY(radians(rotY));
//updateStack(); //future task: implement update position of cubes and perhaps give them a lifetime? Like particles?
makePixBoxes(); //create 3d boxes from pixel locations in the image stack.
showStack();
pop();
}
void mouseWheel(MouseEvent event) {
float e = event.getCount();
zDepth = zDepth + (e*10);
}
void loadLayers(int lCountTemp) {
for (int i = 0; i < lCountTemp; i++) {
layers[i] = loadImage("stack_0" + i + ".png"); //The String "stack_0" nameSpace should be alterable to allow loading other data sets. using this method.
layers[i].loadPixels();
}
}
void createSlices(int lCountTemp) {
for (int i = 0; i < lCountTemp; i++) {
slices[i] = createShape();
}
}
void makeAltLayers(int altLayerNum) {
for (int i = 0; i < altLayerNum; i++) {
altLayers[i] = createGraphics(layers[i].width, layers[i].height);
columns = layers[i].width / cellsize; // Calculate # of columns
rows = layers[i].height / cellsize; // Calculate # of rows
altLayers[i].beginDraw();
altLayers[i].tint(0,random(0,255),0);
altLayers[i].image(layers[i], 0,0);
altLayers[i].endDraw();
altLayers[i].loadPixels();
}
}
void loadLayersPix() {
for (int i = 0; i < altLayers.length; i++) {
altLayers[i].loadPixels();
for (int x = 0; x < altLayers[i].width; x++) {
for (int y = 0; y < altLayers[i].height; y++) {
int loc = x+y*altLayers[i].width;
float b = brightness(altLayers[i].pixels[loc]);
if (b == 255) {
altLayers[i].pixels[loc] = color(0,255,255);
}
//if (b < 255 && b >= 126) {
// altLayers[i].pixels[loc] = color(255);
//}
//else if (b >= 201 && b < 255) {
// altLayers[i].pixels[loc] = color(0,255,0, 150);
//}
//else if (b <= 127) {
// altLayers[i].pixels[loc] = color(0,0);
//}
altLayers[i].pixels[loc] = altLayers[i].pixels[loc];
}
}
altLayers[i].updatePixels();
}
}
void makeSlices() {
for (int i = 0; i < slices.length-1; i++) {
layers[i].loadPixels();
slices[i].beginShape();
slices[i].noStroke();
slices[i].texture(layers[i]);
slices[i].vertex(0, 0, 0, 0);
slices[i].vertex(width, 0, 256, 0);
slices[i].vertex(height, width, 256, 256);
slices[i].vertex(0, height, 0, 256);
slices[i].endShape(CLOSE);
}
}
String[] listFileNames(String dir) {
File file = new File(dir);
if (file.isDirectory()) {
String names[] = file.list();
return names;
} else {
// If it's not a directory
return null;
}
}
void makePixBoxes() {
for (int k = 0; k < altLayers.length-1; k+=cellsize) {
for ( int i = 0; i < columns; i++) {
for ( int j = 0; j < rows; j++) {
int x = i*cellsize + cellsize/2; // x position
int y = j*cellsize + cellsize/2; // y position
int loc = x + y*altLayers[k].width; // Pixel array location
color c = altLayers[k].pixels[loc]; // Grab the color
// Calculate a z position as a function of mouseX and pixel brightness
//float z = (mouseY / float(width)) * brightness(altLayers[k].pixels[loc]) - 20.0;
// Translate to the location, set fill and stroke, and draw the rect
if (altLayers[k].pixels[loc] == (0)) {
continue;
}
else {
pushMatrix();
translate(x, y, k);
fill(255,255,255);
stroke(127);
//rectMode(CENTER);
//rotateY(radians(mouseX));
box(cellsize, cellsize, cellsize);
popMatrix();
}
}
} }
}
void showStack() {
translate(0,0,-altLayers.length/2);
for (int i = 0; i < altLayers.length-1; i++) {
push();
translate(0,0,i);
/* i should be a float calculating the area of gfx of the z height
relationship to its net width/height to the # of images in the
stack(stack height). */
//image(altLayers[i], 0, 0);
//image(layers[i], 100, 100);
shape(slices[i], 0,0);
pop();
}
}
void preLoad() {
//Could preload images files on sepreate thread if needed. Using : requestImage()
}
Any help is greatly appreciated.
Cheers.strong text
