[SOLVED] Reading and manipulating a 3D stack of MRC images

Reading and manipulating a 3D stack of MRC images

The problem I am trying to solve is a three-parter.

I am trying to read a 3D stack of MRC images (not as in: Can DM script read images with the same extension (like *.mrc) from a folder?), apply a function on each image, and output the processed image, preferably in the same order as in the original stack.

Currently, I am reading in the image stack as a file stream and simply displaying it:

file_ID = OpenFileForReading(filename)
file_stream = newStreamFromFileReference(file_ID, 1)
img :=IntegerImage("Imported", 1, 0, size_x, size_y)
ImageReadImageDataFromStream(img, file_stream, 0)

First of all, is this the best way to import a image stack? Once the stack is imported, how can I apply manipulate each image and output them?

Thanks in advance for your input.

Solution

To just "load" an .mrc (stack) file, you can use:

string path = "C:\\temp\\stack.mrc"
imageDocument doc = NewImageDocumentFromFile(path)
doc.ImageDocumentShow()

Or when you want just the image of it:

string path = "C:\\temp\\stack.mrc"
image img := NewImageDocumentFromFile(path).ImageDocumentGetRootImage()
img.ShowImage()

You can then iterate over the stack using the slice commands, f.e. like f.e.

string path = "C:\\temp\\stack.mrc"
image stack := NewImageDocumentFromFile(path).ImageDocumentGetRootImage()
stack.ShowImage()



number sx,sy,sz
stack.ImageGetDimensionSizes(sx,sy,sz)
image sliceMaxProfile := realImage("Maximum of slice",4,sz)
image maskStack := realImage("Masks",4,sx,sy,sz)

for( number i=0;i<sz;i++)
{
    image plane := stack.Slice2(0,0,i, 0,sx,1, 1,sy,1)  // plane is just a reference to the memory (no copy)
    // do any processing
    
    // f.e. store the maximum value of each plane in a profile
    if (0)
        sliceMaxProfile[i,0] = max(plane)   
        
    // or display the plane rotated by 20degree
    if (0){
        image rotPlane := rotate(plane, pi()/180 * 20)
        rotPlane.ImageSetName("Plane #"+i+"+ rotated")
        rotPlane.ShowImage()
    }
    
    // or create a new stack of tresholded plane
    if (0)
     maskStack.Slice2(0,0,i, 0,sx,1, 1,sy,1) = ( plane<mean(plane) ? 1 : 0 )    
}
sliceMaxProfile.ShowImage()
maskStack.ShowImage()