Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
28 changes: 14 additions & 14 deletions Estimate_W.py
Original file line number Diff line number Diff line change
Expand Up @@ -37,9 +37,9 @@ def Wfast(img,nstains,lamb,num_patches,patchsize,level,background_correction=Fal
I_percentile=90

if ydim*xdim>max_size:
print "Finding patches for W estimation:"
print("Finding patches for W estimation:")
for j in range(20):
#print "Patch Sampling Attempt:",i+1
#print("Patch Sampling Attempt:",i+1
initBias=int(math.ceil(patchsize/2)+1)
xx=np.array(range(initBias,xdim-initBias,patchsize))
yy=np.array(range(initBias,ydim-initBias,patchsize))
Expand All @@ -63,20 +63,20 @@ def Wfast(img,nstains,lamb,num_patches,patchsize,level,background_correction=Fal
break
patchsize=int(patchsize*0.95)
valid_inp=np.array(valid_inp)
print "Number of patches sampled for W estimation:", len(valid_inp)
print("Number of patches sampled for W estimation:", len(valid_inp))
else:
patch=np.asarray(img.read_region((0,0),level,(xdim,ydim)))
patch=patch[:,:,:3]
valid_inp=[]
valid_inp.append(patch)
white_pixels= patch[np.sum((patch>white_cutoff),axis=2)==3]
print "Image small enough...W estimation done using whole image"
print("Image small enough...W estimation done using whole image")

if background_correction:
print "Number of white pixels sampled",len(white_pixels)
print("Number of white pixels sampled",len(white_pixels))
if len(white_pixels)<min_num_white:
i0=np.array([255.0,255.0,255.0])
print "Not enough white pixels found, default background intensity assumed"
print("Not enough white pixels found, default background intensity assumed")
elif len(white_pixels)>0:
i0 = np.percentile(white_pixels,I_percentile,axis=0)[:3]
else:
Expand All @@ -101,7 +101,7 @@ def Wfast(img,nstains,lamb,num_patches,patchsize,level,background_correction=Fal
if WS.shape[0]==1:
Wsource=WS[0,:3,:]
else:
print "Median color basis of",len(WS),"patches used as actual color basis"
print("Median color basis of",len(WS),"patches used as actual color basis")
Wsource=np.zeros((3,nstains))
for k in range(nstains):
Wsource[:,k]=[np.median(WS[:,0,k]),np.median(WS[:,1,k]),np.median(WS[:,2,k])]
Expand All @@ -110,15 +110,15 @@ def Wfast(img,nstains,lamb,num_patches,patchsize,level,background_correction=Fal

if Wsource.sum()==0:
if patchsize*0.95<100:
print "No suitable patches found for learning W. Please relax constraints"
print("No suitable patches found for learning W. Please relax constraints")
return None #to prevent infinite recursion
else:
print "W estimation failed, matrix of all zeros found. Trying again..."
print("W estimation failed, matrix of all zeros found. Trying again...")
return Wfast(img,nstains,lamb,min(100,num_patches*1.5),int(patchsize_original*0.95),level)
else:
return Wsource,i0
else:
print "No suitable patches found for learning W. Please relax constraints"
print("No suitable patches found for learning W. Please relax constraints")
return None,None

def patch_Valid(patch,threshold):
Expand All @@ -131,15 +131,15 @@ def patch_Valid(patch,threshold):
temp = tempr*tempg*tempb
r,c = np.shape((temp))
prob= float(np.sum(temp))/float((r*c))
#print prob
#print(prob)
if prob>threshold:
return False
else:
return True

def W_sort(W):
# All sorting done such that first column is H, second column is E
# print W
# print(W)

method = 3

Expand All @@ -153,15 +153,15 @@ def W_sort(W):
# 3. Using r/b ratios of the vectors. H must have a larger value of r/b.
r_b_1 = W[0][0]/W[2][0]
r_b_2 = W[0][1]/W[2][1]
# print r_b_1, r_b_2
# print(r_b_1, r_b_2)
if r_b_1<r_b_2: #else no need to switch
W[:,[0, 1]] = W[:,[1, 0]]
elif method==4:
# 4. Using r-b values of the vectors. H must have a larger value of r-b.
# This is equivalent to comparing the ratios of e^(-r)/e^(-b)
r_b_1 = W[0][0]-W[2][0]
r_b_2 = W[0][1]-W[2][1]
# print r_b_1, r_b_2
# print(r_b_1, r_b_2)
if r_b_1<r_b_2: #else no need to switch
Wsource[:,[0, 1]] = Wsource[:,[1, 0]]

Expand Down
96 changes: 48 additions & 48 deletions Run_ColorNorm.py
Original file line number Diff line number Diff line change
Expand Up @@ -80,7 +80,7 @@ def run_batch_colornorm(filenames,nstains,lamb,output_direc,img_level,background


file_no=0
print "To be normalized:",filenames[1:],"using",filenames[0]
print("To be normalized:",filenames[1:],"using",filenames[0]
for filename in filenames:

display_separator()
Expand All @@ -101,22 +101,22 @@ def run_batch_colornorm(filenames,nstains,lamb,output_direc,img_level,background


tic=time.time()
print

print()
I = openslide.open_slide(filename)
if img_level>=I.level_count:
print "Level",img_level,"unavailable for image, proceeding with level 0"
print("Level",img_level,"unavailable for image, proceeding with level 0")
level=0
else:
level=img_level
xdim,ydim=I.level_dimensions[level]
ds=I.level_downsamples[level]

if file_no==0:
print "Target Stain Separation in progress:",filename,str(xdim)+str("x")+str(ydim)
print("Target Stain Separation in progress:",filename,str(xdim)+str("x")+str(ydim))
else:
print "Source Stain Separation in progress:",filename,str(xdim)+str("x")+str(ydim)
print "\t \t \t \t \t \t \t \t \t \t Time: 0"
print("Source Stain Separation in progress:",filename,str(xdim)+str("x")+str(ydim))
print("\t \t \t \t \t \t \t \t \t \t Time: 0")


#parameters for W estimation
Expand All @@ -127,40 +127,40 @@ def run_batch_colornorm(filenames,nstains,lamb,output_direc,img_level,background

Wi,i0 = Wfast(I,nstains,lamb,num_patches,patchsize,level,background_correction)
if i0 is None:
print "No white background detected"
print("No white background detected")
i0=i0_default

if not background_correction:
print "Background correction disabled, default background intensity assumed"
print("Background correction disabled, default background intensity assumed")
i0=i0_default

if Wi is None:
print "Color Basis Matrix Estimation failed...image normalization skipped"
print("Color Basis Matrix Estimation failed...image normalization skipped")
continue
print "W estimated",
print "\t \t \t \t \t \t Time since processing started:",round(time.time()-tic,3)
print("W estimated"),
print("\t \t \t \t \t \t Time since processing started:",round(time.time()-tic,3))
Wi=Wi.astype(np.float32)

if file_no==0:
print "Target Color Basis Matrix:"
print Wi
print "Target Color Basis Matrix Size:"
print Wi.shape
print("Target Color Basis Matrix:")
print(Wi)
print("Target Color Basis Matrix Size:")
print(Wi.shape)

Wi_target=np.transpose(Wi)
tar_i0=i0
print "Target Image Background white intensity:",i0
print("Target Image Background white intensity:",i0)
else:
print "Source Color Basis Matrix:"
print Wi
print("Source Color Basis Matrix:")
print(Wi)

print "Source Image Background white intensity:",i0
print("Source Image Background white intensity:",i0)

_max=2000
#raise valueError()
print
print()
if (xdim*ydim)<=(_max*_max):
print "Small image processing..."
print("Small image processing...")
img=np.asarray(I.read_region((0,0),level,(xdim,ydim)),dtype=np.float32)[:,:,:3]

Hiv=session1.run(Hiv1,feed_dict={Img1:img, Wis1: Wi,src_i_0:i0})
Expand All @@ -173,44 +173,44 @@ def run_batch_colornorm(filenames,nstains,lamb,output_direc,img_level,background
if file_no==0:
file_no+=1
Hta_Rmax = np.copy(H_Rmax)
print "Target H calculated",
print "\t \t \t \t \t \t \t Total Time:",round(time.time()-tic,3)
print("Target H calculated"),
print("\t \t \t \t \t \t \t Total Time:",round(time.time()-tic,3))
display_separator()
continue

print "Color Normalization in progress..."
print("Color Normalization in progress...")

norm_fac = np.divide(Hta_Rmax,H_Rmax).astype(np.float32)
session1.run(fwrite,feed_dict={shape:np.array(img.shape),Wit1: Wi_target,Hiv2:Hiv,sav_name:s,tar_i_0:tar_i0,normfac:norm_fac})

print "File written to:",s
print "\t \t \t \t \t \t \t \t \t \t Total Time:",round(time.time()-tic,3)
print("File written to:"+s)
print("\t \t \t \t \t \t \t \t \t \t Total Time:",round(time.time()-tic,3))
display_separator()

else:
_maxtf=2550#changed from initial 3000
x_max=xdim
y_max=min(max(int(_maxtf*_maxtf/x_max),1),ydim)
print "Large image processing..."
print("Large image processing...")
if file_no==0:
Hivt=np.memmap('H_target', dtype='float32', mode='w+', shape=(xdim*ydim,2))
else:
Hivs=np.memmap('H_source', dtype='float32', mode='w+', shape=(xdim*ydim,2))
sourcenorm=np.memmap('wsi', dtype='uint8', mode='w+', shape=(ydim,xdim,3))
x_tl = range(0,xdim,x_max)
y_tl = range(0,ydim,y_max)
print "WSI divided into",str(len(x_tl))+"x"+str(len(y_tl))
print("WSI divided into",str(len(x_tl))+"x"+str(len(y_tl)))
count=0
print "Patch-wise H calculation in progress..."
print("Patch-wise H calculation in progress...")
ind=0
perc=[]
for x in x_tl:
for y in y_tl:
count+=1
xx=min(x_max,xdim-x)
yy=min(y_max,ydim-y)
print "Processing:",count," patch size",str(xx)+"x"+str(yy),
print "\t \t Time since processing started:",round(time.time()-tic,3)
print("Processing:",count," patch size",str(xx)+"x"+str(yy)),
print("\t \t Time since processing started:",round(time.time()-tic,3))
img=np.asarray(I.read_region((int(ds*x),int(ds*y)),level,(xx,yy)),dtype=np.float32)[:,:,:3]

Hiv = session1.run(Hiv1,feed_dict={Img1:img, Wis1: Wi,src_i_0:i0})
Expand All @@ -233,23 +233,23 @@ def run_batch_colornorm(filenames,nstains,lamb,output_direc,img_level,background
ind+=len(Hiv)

if file_no==0:
print "Target H calculated",
print("Target H calculated"),
Hta_Rmax = np.percentile(np.array(perc),50,axis=0)
file_no+=1
del Hivt
print "\t \t \t \t \t Time since processing started:",round(time.time()-tic,3)
print("\t \t \t \t \t Time since processing started:",round(time.time()-tic,3))
ind=0
continue

print "Source H calculated",
print "\t \t \t \t \t Time since processing started:",round(time.time()-tic,3)
print("Source H calculated"),
print("\t \t \t \t \t Time since processing started:",round(time.time()-tic,3))
Hso_Rmax = np.percentile(np.array(perc),50,axis=0)
print "H Percentile calculated",
print "\t \t \t \t Time since processing started:",round(time.time()-tic,3)
print("H Percentile calculated"),
print("\t \t \t \t Time since processing started:",round(time.time()-tic,3))

_normfac=np.divide(Hta_Rmax,Hso_Rmax).astype(np.float32)

print "Color Normalization in progress..."
print("Color Normalization in progress...")
count=0
ind=0
np_max=1000
Expand All @@ -258,7 +258,7 @@ def run_batch_colornorm(filenames,nstains,lamb,output_direc,img_level,background
y_max=min(max(int(np_max*np_max/x_max),1),ydim)
x_tl = range(0,xdim,x_max)
y_tl = range(0,ydim,y_max)
print "Patch-wise color normalization in progress..."
print("Patch-wise color normalization in progress...")
total=len(x_tl)*len(y_tl)

prev_progress=0
Expand All @@ -277,21 +277,21 @@ def run_batch_colornorm(filenames,nstains,lamb,output_direc,img_level,background
ind+=pix
percent=5*int(count*20/total) #nearest 5 percent
if percent>prev_progress and percent<100:
print str(percent)+" percent complete...",
print "\t \t \t \t \t Time since processing started:",round(time.time()-tic,3)
print(str(percent)+" percent complete..."),
print("\t \t \t \t \t Time since processing started:",round(time.time()-tic,3))
prev_progress=percent
print "Color Normalization complete!",
print "\t \t \t \t Time since processing started:",round(time.time()-tic,3)
print("Color Normalization complete!"),
print("\t \t \t \t Time since processing started:",round(time.time()-tic,3))

p = time.time()-tic
s = output_direc+base_s.replace(".", "_")+" (using "+base_t.replace(".", "_")+" "+correc+").tif"
print "Saving normalized image..."
print("Saving normalized image...")
#cv2.imwrite(s,cv2.cvtColor(sourcenorm, cv2.COLOR_RGB2BGR))
pyimg = numpy2vips(sourcenorm)
pyimg.tiffsave(s, tile=True, compression='lzw',xres=5000,yres=5000, bigtiff=True,pyramid=True,Q=100)#xres and yres should be controlled to produce finer or coarser tif
del sourcenorm
print "File written to:",s
print "\t \t \t \t \t \t \t \t \t Total Time:",round(time.time()-tic,3)
print("File written to:"+s)
print("\t \t \t \t \t \t \t \t \t Total Time:",round(time.time()-tic,3))
display_separator()

file_no+=1
Expand All @@ -313,5 +313,5 @@ def run_colornorm(source_filename,target_filename,nstains,lamb,output_direc,leve


def display_separator():
print "________________________________________________________________________________________________"
print
print("________________________________________________________________________________________________")
print()
16 changes: 8 additions & 8 deletions Run_StainSep.py
Original file line number Diff line number Diff line change
Expand Up @@ -12,21 +12,21 @@

def run_stainsep(filename,nstains,lamb,output_direc="",background_correction=True):

print
print "Running stain separation on:",filename
print()
print("Running stain separation on:",filename)

level=0

I = openslide.open_slide(filename)
xdim,ydim=I.level_dimensions[level]
img=np.asarray(I.read_region((0,0),level,(xdim,ydim)))[:,:,:3]

print "Fast stain separation is running...."
print("Fast stain separation is running....")
Wi,Hi,Hiv,stains=Faststainsep(I,img,nstains,lamb,level,background_correction)

#print "\t \t \t \t \t \t Time taken:",elapsed
#print("\t \t \t \t \t \t Time taken:",elapsed

print "Color Basis Matrix:\n",Wi
print("Color Basis Matrix:\n",Wi)

fname=os.path.splitext(os.path.basename(filename))[0]
cv2.imwrite(output_direc+fname+"-0_original.png",cv2.cvtColor(img, cv2.COLOR_RGB2BGR))
Expand All @@ -38,7 +38,7 @@ def Faststainsep(I_obj,I,nstains,lamb,level,background_correction):
s=I.shape
ndimsI = len(s)
if ndimsI!=3:
print "Input Image I should be 3-dimensional!"
print("Input Image I should be 3-dimensional!")
sys.exit(0)
rows = s[0]
cols = s[1]
Expand All @@ -51,10 +51,10 @@ def Faststainsep(I_obj,I,nstains,lamb,level,background_correction):


if background_correction:
print "Background intensity:",i0
print("Background intensity:",i0)
else:
i0 = np.array([255.,255.,255.])
print "Background correction disabled, default background intensity assumed"
print("Background correction disabled, default background intensity assumed")

#Beer-Lambert tranformation
V,VforW=BLtrans(I,i0) #V=WH see in paper
Expand Down
Loading