4

I have series of 2D tiff images of a sample, I want to create or reproduce 3D image/volume using those 2d image for 3D visualization.
I found this link Reconstructing 3D image from 2D image have similar question but It discussed about CT reconstruction using backprojection algorithm. But I already have 2D view of sample in image form.

I want to know how can I reproduce 3D image from those 2D slices(Tiff image) using python or Matlab.

akshahi17
  • 129
  • 1
  • 2
  • 11
  • Read the data into a container that support n-dimensional data - Concatenate the 2-d data along a third axis/dimension? - [https://numpy.org/doc/stable/reference/routines.array-manipulation.html#joining-arrays](https://numpy.org/doc/stable/reference/routines.array-manipulation.html#joining-arrays) – wwii Mar 18 '21 at 22:36
  • Are those 2D slices of a volume scan (e.g. MRI / CT scan) or something else ? (if you could post an example tiff that might helpful for others to advise) – George Profenza Mar 18 '21 at 22:54
  • Using [napari](https://napari.org/), visualizing the image series might be as easy as `python -m napari *.tif` and then switch to 3D view. – cgohlke Mar 19 '21 at 01:28
  • Yes, 2D sciences are from Micro-CT volume scanned – akshahi17 Mar 20 '21 at 05:19
  • @wwii It's a kind of image stacking to form a 3D model. – akshahi17 Mar 20 '21 at 05:27

2 Answers2

6

I wanna check that this is what you're looking for before I go on a long explanation of something that could be irrelevant.

I have a series of 2d images of a tumor. I'm constructing a 3d shell from the image slices and creating a .ply file from that shell.

2D slices

enter image description here

3D Reconstruction

enter image description here

Is this the sort of thing that you're looking for?

Edit:

I downloaded the dataset and ran it through the program.

enter image description here

I set the resolution of the image to 100x100 to reduce the number of points in the .ply file. You can increase it or decrease it as you like.

Program for creating the .ply file

import cv2
import math
import numpy as np
import os

# creates a point cloud file (.ply) from numpy array
def createPointCloud(filename, arr):
    # open file and write boilerplate header
    file = open(filename, 'w');
    file.write("ply\n");
    file.write("format ascii 1.0\n");

    # count number of vertices
    num_verts = arr.shape[0];
    file.write("element vertex " + str(num_verts) + "\n");
    file.write("property float32 x\n");
    file.write("property float32 y\n");
    file.write("property float32 z\n");
    file.write("end_header\n");

    # write points
    point_count = 0;
    for point in arr:
        # progress check
        point_count += 1;
        if point_count % 1000 == 0:
            print("Point: " + str(point_count) + " of " + str(len(arr)));

        # create file string
        out_str = "";
        for axis in point:
            out_str += str(axis) + " ";
        out_str = out_str[:-1]; # dump the extra space
        out_str += "\n";
        file.write(out_str);
    file.close();


# extracts points from mask and adds to list
def addPoints(mask, points_list, depth):
    mask_points = np.where(mask == 255);
    for ind in range(len(mask_points[0])):
        # get point
        x = mask_points[1][ind];
        y = mask_points[0][ind];
        point = [x,y,depth];
        points_list.append(point);

def main():
    # tweakables
    slice_thickness = .2; # distance between slices
    xy_scale = 1; # rescale of xy distance

    # load images
    folder = "images/";
    files = os.listdir(folder);
    images = [];
    for file in files:
        if file[-4:] == ".tif":
            img = cv2.imread(folder + file, cv2.IMREAD_GRAYSCALE);
            img = cv2.resize(img, (100, 100)); # change here for more or less resolution
            images.append(img);

    # keep a blank mask
    blank_mask = np.zeros_like(images[0], np.uint8);

    # create masks
    masks = [];
    masks.append(blank_mask);
    for image in images:
        # mask
        mask = cv2.inRange(image, 0, 100);

        # show
        cv2.imshow("Mask", mask);
        cv2.waitKey(1);
        masks.append(mask);
    masks.append(blank_mask);
    cv2.destroyAllWindows();

    # go through and get points
    depth = 0;
    points = [];
    for index in range(1,len(masks)-1):
        # progress check
        print("Index: " + str(index) + " of " + str(len(masks)));

        # get three masks
        prev = masks[index - 1];
        curr = masks[index];
        after = masks[index + 1];

        # do a slice on previous
        prev_mask = np.zeros_like(curr);
        prev_mask[prev == 0] = curr[prev == 0];
        addPoints(prev_mask, points, depth);

        # # do a slice on after
        next_mask = np.zeros_like(curr);
        next_mask[after == 0] = curr[after == 0];
        addPoints(next_mask, points, depth);

        # get contour points (_, contours) in OpenCV 2.* or 4.*
        _, contours, _ = cv2.findContours(curr, cv2.RETR_TREE, cv2.CHAIN_APPROX_NONE);
        for con in contours:
            for point in con:
                p = point[0]; # contours have an extra layer of brackets
                points.append([p[0], p[1], depth]);

        # increment depth
        depth += slice_thickness;

    # rescale x,y points
    for ind in range(len(points)):
        # unpack
        x,y,z = points[ind];

        # scale
        x *= xy_scale;
        y *= xy_scale;
        points[ind] = [x,y,z];

    # convert points to numpy and dump duplicates
    points = np.array(points).astype(np.float32);
    points = np.unique(points.reshape(-1, points.shape[-1]), axis=0);
    print(points.shape);

    # save to point cloud file
    createPointCloud("test.ply", points);

if __name__ == "__main__":
    main();
Ian Chu
  • 2,924
  • 9
  • 14
  • Yes, Something like what you have mentioned. Here is an example of my tiff file [sample image][1] This is what I'm looking for [Desire output][2] [1]: https://i.stack.imgur.com/frbQy.png [2]: https://i.stack.imgur.com/9abv9.png – akshahi17 Mar 20 '21 at 05:16
  • Can you post a link to your images? – Ian Chu Mar 20 '21 at 13:17
  • Also, you show a cube as your desired output. Do you care about what's inside the cube, or is it enough to just see the textured faces of the cube? – Ian Chu Mar 20 '21 at 13:18
  • here is the link to 2D slices, https://www.digitalrocksportal.org/projects/92/origin_data/384/ – akshahi17 Mar 20 '21 at 17:04
  • Actually, each pixel contains some information in 2D that I want to see in 3D as well. – akshahi17 Mar 20 '21 at 17:07
  • Hi, @IanChu I ran the code you provided, however I couldn't open the .ply file generated. How could I open the file to view? – Nicholas TJ May 23 '21 at 08:49
  • I used Open3D (https://pypi.org/project/open3d/) to view the file. You can use anything that can open a .ply file though, it's a standard file format for 3d stuff. – Ian Chu May 24 '21 at 14:26
2

Another option is to apply open-source library MeshLib, which can be called both from C++ and python code.

Assuming we have bunch of slices like following: bunch of 2D slices

On Windows, the simplest way to install it in python 3.10 is

py -3.10 -m pip install --upgrade meshlib

And the code for constructing 3D object from a series of images will look like

import meshlib.mrmeshpy as mr
settings = mr.LoadingTiffSettings()

# load images from specified directory
settings.dir = "C:/slices"

# specifiy size of 3D image element
settings.voxelSize = mr.Vector3f(1, 1, 5)

#create voxel object from the series of images
volume = mr.loadTiffDir(settings)

#define ISO value to build surface
iso=127.0

#convert voxel object to mesh
mesh=mr.gridToMesh(volume, iso)

#save mesh to .stl file
mr.saveMesh(mesh, mr.Path("C:/slices/mesh.stl"))

Finally, result looks like 3D image

Fedor
  • 17,146
  • 13
  • 40
  • 131
Alex Koksin
  • 129
  • 1
  • 6