Commit 020c7497 authored by Martinez-Carvajal German's avatar Martinez-Carvajal German
Browse files

New script Connect CR (connected regions) !!!!!

Yes it sounds weird...

FIJI will compute the connected regions (CR) but
you may have RAM limitations. So you will have
to split your original image in different parts.

This split will connect this parts as long as you
have left one slice of "recouvrement" (covering zone)
when you splitted the original image.

The programm deals fine with all the possibilities
of connections due to the splitting of the image
and will reassign new labels for every pore.
parent 68c5cf95
# -*- coding: utf-8 -*-
"""
Created on Wed Sep 11 16:40:02 2019
@author: german.martinez-carvajal
"""
import tifffile
from os import chdir
import numpy as np
def sector_mask(shape,centre,radius,angle_range):
"""
CODE FROM STACK OVER FLOW
Return a boolean mask for a circular sector. The start/stop angles in
`angle_range` should be given in clockwise order.
"""
x,y = np.ogrid[:shape[0],:shape[1]]
cx,cy = centre
tmin,tmax = np.deg2rad(angle_range)
# ensure stop angle > start angle
if tmax < tmin:
tmax += 2*np.pi
# convert cartesian --> polar coordinates
r2 = (x-cx)*(x-cx) + (y-cy)*(y-cy)
theta = np.arctan2(x-cx,y-cy) - tmin
# wrap angles between 0 and 2*pi
theta %= (2*np.pi)
# circular mask
circmask = r2 <= radius*radius
# angular mask
anglemask = theta <= (tmax-tmin)
return circmask*anglemask
def eliminating_circle_artifact(image):
print("eliminating_circle_artifact ...")
# CREATING CIRCULAR MASK
# The cylinder is supposed to be the biggest "circle" inscribed in a squared image
z,x,y = np.shape(image)[0], np.shape(image)[1], np.shape(image)[2]
center = int(x/2), int(y/2)
radius = int((x/2))-10
print("radius", radius)
#Selecting circular ROI (region of interest)
mask_2D = sector_mask((x,y),center,radius,(0,360 ))
mask_3D = np.zeros(image.shape, dtype = np.bool)
for j in range(z):
mask_3D[j] = mask_2D
image[np.logical_not(mask_3D)] = 0
return image
def connect_CR(path_read, file_top, file_bot, path_save, sample_name):
chdir(path_read)
print("reading ... ", path_read, file_top)
top = tifffile.imread(file_top)
print("reading ... ", path_read, file_bot)
bot = tifffile.imread(file_bot)
a = len(np.unique(top))
b = len(np.unique(bot))
# checking if there will be enough labels
assert a < 2**16-1
assert b < 2**16-1
assert a+b < 2**16-1
print("eliminating circle artifact")
top = eliminating_circle_artifact(top)
bot = eliminating_circle_artifact(bot)
chdir(path_save)
tifffile.imsave('bot.tif', bot)
tifffile.imsave('top.tif', top)
a = len(np.unique(top))
b = len(np.unique(bot))
print("relabeling bot to eliminate equal labels ...")
k = max(np.unique(top))
print("k = " , k)
mask = (bot != 0)
bot = bot + k
bot = bot*mask
# connecting slices
cs_top = top[-1]
cs_bot = bot[0]
print("checking recouvrement...")
final = (cs_top==0)
start = (cs_bot==0)
assert (np.alltrue(final == start))
print("getting labels from pores at the connecting slice...")
list_pores_top = np.unique(cs_top)
list_pores_bot = np.unique(cs_bot)
print("list pores top :",list_pores_top )
print("list pores bot :",list_pores_bot )
print("checking if the number is equal for both parts: bot and top")
print("Number of pores (top), (bottom): ", len(list_pores_top), len(list_pores_bot))
if len(list_pores_top) != len(list_pores_bot):
print("THEY ARE NOT THE SAME !!!!!!")
else:
print("they are the same")
z1, y1, x1 = np.shape(top)
z2, y2, x2 = np.shape(bot)
assert ((y1==y2) and (x1==x2))
len_x = x1
len_y = y1
print("getting all possible couples")
couples = set()
for x in range(len_x):
for y in range(len_y):
couple = (cs_top[y,x], cs_bot[y,x])
couples.add(couple)
assert ~((couple[0] == 0) and couple[1] != 0)
assert ~((couple[0] != 0) and couple[1] == 0)
print('reducing couples by analyzing connections')
# couples reduced
couples_r = list()
flag_skip_label_top = False
for label_top in list_pores_top:
# do not take into account label 0
if label_top == 0:
continue
# skipping the step if label_top has already been connected
if len(couples_r)!= 0:
for i in range(len(couples_r)):
if label_top in couples_r[i][0]:
#print("skipping Label Top:", label_top )
flag_skip_label_top = True
if flag_skip_label_top == True:
flag_skip_label_top = False
continue
print ("!!!!!! Analyzing Label Top: ", label_top, "!!!!!!!")
set_label_top = set()
set_label_bot = set()
set_label_top.add(label_top)
flag_iteration_connections = True
while flag_iteration_connections:
#print('increasing list of labels bot...')
for l_top in set_label_top:
#print(l_top)
for couple in couples:
if couple[0]== l_top:
#print(couple)
set_label_bot.add(couple[1])
#print(sorted(set_label_top) , sorted(set_label_bot))
#print('increasing list of labels top ...')
set_label_top_new = set()
for label_bot in set_label_bot:
#print(label_bot)
for couple in couples:
if couple[1] == label_bot:
#print(couple)
set_label_top_new.add(couple[0])
#print(sorted(set_label_top_new) , sorted(set_label_bot))
if set_label_top_new != set_label_top:
print("repeating")
flag_iteration_connections = True
set_label_top = set_label_top_new
else:
flag_iteration_connections = False
couples_r.append((sorted(set_label_top) , sorted(set_label_bot)))
print("len of final couple list:", len(couples_r))
print("Relabelling each part (bot and top)...")
top_d = 0
bot_d = 0
for i in range(len(couples_r)):
top_d -= 1
bot_d -= 1
labels_top = couples_r[i][0]
#print("labels top: \n", labels_top)
labels_bot = couples_r[i][1]
#print("labels bot: \n", labels_bot)
min_top = min(labels_top)
min_bot = min(labels_bot)
#print(min_top, min_bot)
for label_top in labels_top:
mask = (top == label_top)
top[mask] = min_top
top_d += 1
for label_bot in labels_bot:
mask = (bot == label_bot)
bot[mask] = min_bot
bot_d += 1
print("checking if the number is equal for both parts: bot and top...")
cs_top = top[-1]
cs_bot = bot[0]
list_pores_top = np.unique(cs_top)
list_pores_bot = np.unique(cs_bot)
assert(len(list_pores_top)==len(list_pores_bot))
print("getting all possible couples")
couples = set()
for x in range(len_x):
for y in range(len_y):
couple = (cs_top[y,x], cs_bot[y,x])
couples.add(couple)
assert ~((couple[0] == 0) and couple[1] != 0)
assert ~((couple[0] != 0) and couple[1] == 0)
assert (len(list_pores_top) == len(couples))
couples = couples.difference({(0,0)})
print("couples = \n", sorted(couples))
print("Assiging the same label from top to bottom where there are connections...")
for couple in couples:
if couple == (0,0):
continue
assert(couple[1] > k)
mask = (bot == couple[1])
bot[mask] = couple[0]
print("Maximal cheking ...")
cs_top = top[-1]
cs_bot = bot[0]
for x in range(len_x):
for y in range(len_y):
assert cs_top[y,x] == cs_bot[y,x]
print("Connecting...")
connected_image = np.concatenate((top, bot[1:]))
pores = np.unique(connected_image)
print(len(pores))
pores = pores[ pores != 0]
print("Final checking...")
r = (a-1)-top_d + (b-1)-bot_d -len(couples)
print(len(pores), a, b, top_d, bot_d, len(couples), r)
assert(len(pores) == r )
print("Final relabel...")
pores = np.sort(np.unique(connected_image))
final_label = 0
for pore in pores:
mask = (connected_image == pore)
connected_image[mask] = final_label
final_label +=1
chdir(path_save)
tifffile.imsave("CR_{}.tif".format(sample_name) , connected_image)
print('Connection OK !')
path_read = "/home/german.martinez-carvajal/Bureau/These/Connecting_CR/test"
path_save = "/home/german.martinez-carvajal/Bureau/These/Connecting_CR/test"
file_top = 'small_top.tif'
file_bot = 'small_bot.tif'
sample_name = "small"
connect_CR(path_read, file_top, file_bot, path_save, sample_name)
\ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment