0


I would like to transfer my image processing from Image J (Fiji) to Python.
In Image J, I split the image into HSB, then use the Moments Auto-threshold on the B channel. On Python, I wish to have the value of the threshold t from which the segmentation is done. As I couldn't find any help or piece of code on the subject, here I am.
From the paper "moment-preserving thresholding: a new approach, Tsai" (here), I did this:

import numpy as np
import cv2
import skimage
from skimage import io

img = io.imread("C:\\Users\\Image.tif")
hsv_img = cv2.cvtColor(filt_img, cv2.COLOR_RGB2HSV)    
H, S, B = cv2.split(img) # spliting the image into HSB

B_histo = skimage.exposure.histogram (B) # making the histogram

pix_sum = B.shape[0]*B.shape[1] # calculating the sum of the pixels

#from the paper, calculating the 4 first odrers m0, m1, m2, m3 to get to p0. The name of the further variables stems from the paper. 

pj = B_histo[0] / pix_sum 
pj_z1 = np.power(B_histo[1], 1) * pj
pj_z2 = np.power(B_histo[1], 2) * pj
pj_z3 = np.power(B_histo[1], 3) * pj

m0 = np.sum(pj)
m1 = np.sum(pj_z1)
m2 = np.sum(pj_z2)
m3 = np.sum(pj_z3)

cd = (m0*m2) - (m1*m1)
c0 = ((-m2*m2) - (-m3*m1))/cd
c1 = ((m0*-m3) - (m1*-m2))/cd

z0 = 0.5 *(-c1 - (np.power(np.power(c1, 2) - 4*c0, 1/2)))
z1 = 0.5 *(-c1 + (np.power(np.power(c1, 2) - 4*c0, 1/2)))

pd = z1 - z0
p0 = (z1 - m1) / pd # p0 should be the percentage of the pixels to which the value 0 is attributed

# using cumulative histogram and comparing it to a target value by calculating the difference. When the difference is the lowest, the index indicates the value of the threshold t
cum_pix = np.cumsum(B_histo[0]) 
target_value = p0 * pix_sum
#td = cum_pix
diff = [abs(i - target_value) for i in cum_pix]
cum_pix[1]
diff[0]
t = [abs(i - target_value) for i in cum_pix].index(np.min(diff))

print(t)

Sorry for the messy code. Anyway, the value I calculate on Python is not the same as the one on Image J. Do you know where the problem might come from ? Or is there a function in Python which can retrieve the value of the moment auto-threshold ? I would greatly apprciate tips or hints to dig, thanks

PaPe
  • 1
  • 1
  • It would be helpful if you provided the input image, the **Image J** result (including any intermediate values you can extract) and your Python result. – Mark Setchell Dec 12 '18 at 17:46
  • Hi Mark, I'm sorry I didn't add any material to the post but I found an answer to my question. The problem came from skimage.exposure.histogram module which do not give automatically 256 bins for image histogram, even when specified in the arguments. I still don't know why skimage doesn't work properly on my computer but I found a way to overcome it (even if it is more messy) – PaPe Dec 17 '18 at 08:42

1 Answers1

0

Calculations for the moments and proportions to threshold were good but there was a problem with skimage histogram which shrink the bins (even if 256 bins was specified) and the way t was calculated. Here is a solution to the problem:

#Instead of skimage.exposure.histogram (B, nbins=256)   
grey_value = np.arange(256)
B_freq = cv2.calcHist([S], [0],None,[256],[0,255])
B_freq = B_freq.reshape((256,))
B_freq = np.int64(B_freq)
B_histo = (B_freq, grey_value)

pix_sum = B.shape[0]*B.shape[1] # calculating the sum of the pixels

#from the paper, calculating the 3 first odrers

pj = B_histo[0] / pix_sum 
pj_z1 = np.power(B_histo[1], 1) * pj
pj_z2 = np.power(B_histo[1], 2) * pj
pj_z3 = np.power(B_histo[1], 3) * pj

m0 = np.sum(pj)
m1 = np.sum(pj_z1)
m2 = np.sum(pj_z2)
m3 = np.sum(pj_z3)

cd = (m0*m2) - (m1*m1)
c0 = ((-m2*m2) - (-m3*m1))/cd
c1 = ((m0*-m3) - (m1*-m2))/cd


z0 = 0.5 *(-c1 - (np.power(np.power(c1, 2) - 4*c0, 1/2)))
z1 = 0.5 *(-c1 + (np.power(np.power(c1, 2) - 4*c0, 1/2)))

pd = z1 - z0
p0 = (z1 - m1) / pd # p0 should be the percentage of the pixels to which the threshold t should be done

# using cumulative histogram and comparing it to a target value by calculating the difference. When the difference is the lowest, the index indicates the value of the threshold t
cum_pix = np.cumsum(B_freq) 
target_value = p0 * pix_sum

diff = [(i - target_value) for i in cum_pix]

def find_nearest(array, value):
    array = np.asarray(array)
    idx = (np.abs(array - value)).argmin()
    return array[idx]

t = diff.index(find_nearest(diff, 0))

print(t)

So this code gives you the value of threshold of the Moments auto-threshold of Image J.

PaPe
  • 1
  • 1