Saving an Image matrix with imageio for opening using ImageJ

Question

Hey so I am trying to save an images that are in the form (10, 1926, 1926, 9). That is 10 z-stacks of 1926x1926 pixel images over 9 channels. However, when I try to save the processed images using the tiffile package and open them up in imageJ they seem open up differently the the preprocessed images ie. the images will be 10x1926 pixels or it doesn't differentiate the z-stacks from the channel resulting in 90 stacks instead of 10 stacks with 9 channels. Should I be using another file format or something? My code is below

referenceImage = io.imread("C:\\Users\\antho\\OneDrive\\Documents\\NerettiLab\\Montage-000-

000_reduced\\img000.ome.tif")
dapi_channel = referenceImage[:,:,:,0]
bead_channel = referenceImage[:,:,:,1]
CT_channel = referenceImage[:,:,:,2]

referenceImage2 = io.imread("C:\\Users\\antho\\OneDrive\\Documents\\NerettiLab\\Montage-000-000_reduced\\img001.ome.tif")
dapi_channel2 = referenceImage2[:,:,:,0]
bead_channel2 = referenceImage2[:,:,:,1]
CT_channel2 = referenceImage2[:,:,:,2]


corrected_image = np.empty_like(referenceImage)
corrected_image[:,:,:,0] = dapi_channel 
corrected_image[:,:,:,1] = bead_channel 
corrected_image[:,:,:,2] = CT_channel 

for file in range(0,len(files)):
    
    referenceImage = io.imread(myCTpath + files[0])
    dapi_channel = referenceImage[:,:,:,0]
    bead_channel = referenceImage[:,:,:,1]
    CT_channel = referenceImage[:,:,:,2]

    referenceImage2 = io.imread(myCTpath + files[file])
    
    dapi_channel2 = referenceImage2[:,:,:,0]
    bead_channel2 = referenceImage2[:,:,:,1]
    CT_channel2 = referenceImage2[:,:,:,2]
    
    images_list = []
    for img in range(0,10): 
        image1_gray = bead_channel[img]
        # Normalize pixel values to the range [0, 1]
        normalized_image = (image1_gray - image1_gray.min()) / (image1_gray.max() - image1_gray.min())

        # Scale the normalized values to the range [0, 255]
        scaled_image1 = (normalized_image * 255).astype(np.uint8)
        image2_gray = bead_channel2[img]

        # Normalize pixel values to the range [0, 1]
        normalized_image = (image2_gray - image2_gray.min()) / (image2_gray.max() - image2_gray.min())

        # Scale the normalized values to the range [0, 255]
        scaled_image2 = (normalized_image * 255).astype(np.uint8)
        shift, error, diffphase = phase_cross_correlation(scaled_image1, scaled_image2)
        shift_x = int(shift[0])
        shift_y = int(shift[1])
        shifted_image2 = np.empty_like(scaled_image2)
        #print(scaled_image2.shape)
        shifted_image2 = np.roll(scaled_image2, (shift_x, shift_y), axis=(0, 1))
        images_list.append(shifted_image2)
    
    new_images_stacked = np.stack(images_list, axis=0)
    corrected_image = np.concatenate((corrected_image, new_images_stacked[..., np.newaxis]), axis=-1)


new_images_stacked = np.stack(images_list, axis=0)

corrected_image = np.concatenate((corrected_image, new_images_stacked[..., np.newaxis]), axis=-1)


# Choose the path and filename for the saved image
output_file = 'output_image.tif'

# Save the 4D image as a 3D stack using tifffile
tiff.imsave(myCTpath+output_file, image_rearranged)

# Check if the image was saved successfully
if isfile(myCTpath+output_file):
    print(f"Image saved as {output_file}")
else:
    print("Image saving failed")