pythonimage-processing3d-reconstruction

How to create a 3D image with series of 2D Image


I have a series of 2D tiff images of a sample; I want to create or reproduce a 3D image/volume((volume rendering) using the 2D image for 3D visualization.
I found that this link Reconstructing 3D image from 2D image has a similar problem but discusses CT reconstruction using a back projection algorithm. However, I already have a 2D view of the sample in image form.

I want to know how to reproduce a 3D image(volume rendering) from those 2D slices(Tiff image) using Python or Matlab.


Solution

  • 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();