SegmentationProjectGenerator.py

import syglass as sy
from syglass import pyglass
import numpy as np
import cv2
import tifffile as tiff
import os

PROJECT_PATH = "D:\\SyGlass-Stuff\\SyGlassProjects\\syFind Files\\2023-10-26_P0_S1_Aldh1L1 GFP_D GFP Ald M2_20x\\2023-10-26_P0_S1_Aldh1L1 GFP_D GFP Ald M2_20x.syg"

def main(projectPath: str):
    project = sy.get_project(projectPath)
    resolution_map = project.get_resolution_map()
    max_resolution_level = len(resolution_map) - 1
    blockSize = project.get_block_size() # blockSize is z, y, x
    projectSize = pyglass.vec3(float(blockSize[2] * (2 ** max_resolution_level)), 
                               float(blockSize[1] * (2 ** max_resolution_level)), 
                               float(blockSize[0] * (project.impl.GetTreeK() / 4 ** max_resolution_level)))
    # extractor needs x, y, z
    extractor = pyglass.MaskOctreeRasterExtractor(None)
    raster = extractor.GetCustomBlock(project.impl, 0, max_resolution_level, pyglass.vec3(0, 0, 0), projectSize)
    maskData: np.ndarray = pyglass.GetRasterAsNumpyArray(raster)
    emptyMaskArray = np.zeros_like(maskData)
    individualMaskBBs = []

    ### Identify Mask Bounding Boxes
    for i in range(maskData.max()):
        # bbox stores the min and max coordinates of the mask in z, y, x
        bbox = [np.array([maskData.shape[0], maskData.shape[1], maskData.shape[2]]), np.array([0, 0, 0])]
        temp = np.where(maskData == (i + 1), maskData, emptyMaskArray)
        for z in range(temp.shape[0]):
            nonzero = cv2.findNonZero(temp[z])
            if (nonzero is not None) and (nonzero.size > 0):
                if z < bbox[0][0]:
                    bbox[0][0] = z
                if z > bbox[1][0]:
                    bbox[1][0] = z
            rectangle = cv2.boundingRect(nonzero)
            if rectangle[2]:
                if rectangle[0] < bbox[0][2]:
                    bbox[0][2] = rectangle[0]
                if rectangle[0] + rectangle[2] > bbox[1][2]:
                    bbox[1][2] = rectangle[0] + rectangle[2]
                if rectangle[1] < bbox[0][1]:
                    bbox[0][1] = rectangle[1]
                if rectangle[1] + rectangle[3] > bbox[1][1]:
                    bbox[1][1] = rectangle[1] + rectangle[3] 
        individualMaskBBs.append(bbox)
    ###
    
    parentPath = os.path.dirname(PROJECT_PATH)
    tiffPath = os.path.join(parentPath, "tiffs")
    if not os.path.isdir(tiffPath):
        os.mkdir(tiffPath)
    settings = sy.ProjectCreationSettings()
    settings.voxel_dimensions = project.get_voxel_dimensions()
    settings.voxel_unit = project.get_voxel_unit()
    successfulProjectCount = 0
    ### Create Projects from Bounding Boxes
    for mask in range(len(individualMaskBBs)):
        bb = individualMaskBBs[mask]
        print(f"Mask {mask + 1} BB: xmin={bb[0][2]} ymin={bb[0][1]} zmin={bb[0][0]} xmax={bb[1][2]} ymax={bb[1][1]} zmax={bb[1][0]}")
        # Get the project data from in the bounding box
        customBlock = project.get_custom_block(0, max_resolution_level, bb[0], bb[1] - bb[0])
        # Isolate the mask data in the bounding box
        maskBlock = maskData[bb[0][0]:bb[1][0], bb[0][1]:bb[1][1], bb[0][2]:bb[1][2]]

        # Zero the unmasked voxels in the bb
        maskedCustomBlock = np.where(maskBlock == (mask+1), customBlock.data, np.zeros_like(customBlock.data))

        newProjectPath = os.path.join(os.path.dirname(parentPath), project.get_name() + "_Masks")
        newProjectTiffPath = os.path.join(tiffPath, "Mask_" + str(mask + 1))
        if not os.path.isdir(newProjectTiffPath):
            os.mkdir(newProjectTiffPath)
        # Write each channel as a separate image because I can't figure out how to write a 4 channel tiff that is recognized in syGlass
        imagePaths = []
        for channel in range(maskedCustomBlock.shape[-1]):
            imagePath = os.path.join(newProjectTiffPath, f"Image_{channel}.tiff")
            tiff.imwrite(imagePath, maskedCustomBlock[:, :, :, channel:channel+1], bigtiff=True)
            imagePaths.append(imagePath)
        successfulProjectCount += create_project(imagePaths, newProjectPath, project.get_name() + "_Mask_" + str(mask + 1), settings)
    ###
    print(f"Successfully created {successfulProjectCount} new sub-projects")
    
""" 
imagePaths is a list of channels 
settings only overrides the voxel dimensions and voxel unit
""" 
from typing import Union
def create_project(imagePaths, projectPath: str, projectName: str, settings: Union[sy.ProjectCreationSettings, None] = None):
    try:
        project: pyglass.Project = pyglass.CreateProject(pyglass.path(projectPath), projectName)
        dataProviders = []
        for imagePath in imagePaths:
            dataProvider = pyglass.OpenTIFF(imagePath)
            if dataProvider.Valid():
                dataProviders.append(dataProvider)
        finalDataProvider: pyglass.DataProvider
        if len(dataProviders) == 1:
            finalDataProvider = dataProviders[0]
        elif len(dataProviders) == 2:
            finalDataProvider = pyglass.Combine2(dataProviders[0], dataProviders[1])
        elif len(dataProviders) == 3:
            finalDataProvider = pyglass.Combine3(dataProviders[0], dataProviders[1], dataProviders[2])
        elif len(dataProviders) == 4:
            finalDataProvider = pyglass.Combine4(dataProviders[0], dataProviders[1], dataProviders[2], dataProviders[3])
        else:
            raise ValueError("More than 4 channels found, aborting.")
        
        cd = pyglass.ConversionDriver()
        cd.SetInput(finalDataProvider)
        cd.SetOutput(project)
        cd.StartAsynchronous()
        lowPercentage = 0
        print(cd.GetPercentage())
        while cd.GetPercentage() < 100:
            if cd.GetPercentage() > (lowPercentage + 10):
                lowPercentage = cd.GetPercentage()
                print("Progress: " + str(cd.GetPercentage())[0:4] + "%")
    except Exception as e:
        print("Something went wrong...")
        print(str(e))
        return 0
    else:
        print("Success...!")
        if settings is not None:
            voxel_dimensions = np.array([settings.voxel_dimensions[2], settings.voxel_dimensions[1], settings.voxel_dimensions[0]])
            voxel_size = pyglass.vec3(*voxel_dimensions)
            project.SetVoxelSize(voxel_size)
            project.SetVoxelUnit(settings.voxel_unit)
        return 1
    
    
if __name__ == "__main__":
    main(PROJECT_PATH)