file-formatdm-script

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