Python color.separate_stains函数代码示例

 def test_hdx_rgb_roundtrip(self):
     from skimage.color.colorconv import hdx_from_rgb, rgb_from_hdx
     img_rgb = self.img_rgb
     conv = combine_stains(separate_stains(img_rgb, hdx_from_rgb),
                           rgb_from_hdx)
     with expected_warnings(['precision loss']):
         assert_equal(img_as_ubyte(conv), img_rgb)

def rgb2he(img, normalize=False):
    """
    RGB2HE: Extracts the haematoxylin and eosin components from an RGB color.

    h,e = rgb2he(img, normalize=False)

    Args:
        img (numpy.ndarray): and RGB image; no check for actual color space is
        performed
        
        normalize (bool): should the values be linearly transformed, to ensure
        they are all between -1.0 and 1.0  

    Returns:
        tuple. Contains two intensity images as numpy.ndarray, coding for the
        haematoxylin and eosin components, respectively. The values are in the
        "H&E" space or in -1..1 (if normalize==True)
    """

    # my color separation matrices
    # (from http://www.mecourse.com/landinig/software/colour_deconvolution.zip)
    # The code commented out below was used to generate the "he1_from_rgb" matrix.
    # After the matrix was obtained, it has been hard coded, for computation 
    # efficiency:
    # rgb_from_he1 = np.array([[0.644211, 0.716556, 0.266844],
    #                         [0.092789, 0.954111, 0.283111],
    #                         [0.0, 0.0, 0.0]])
    # rgb_from_he1[2, :] = np.cross(rgb_from_he1[0, :], rgb_from_he1[1, :])
    # he1_from_rgb = linalg.inv(rgb_from_he1)
    he1_from_rgb = np.array([[ 1.73057512, -1.3257525 , -0.1577248 ],
                             [-0.19972397,  1.1187028 , -0.48055639],
                             [ 0.10589662,  0.19656106,  1.67121469]])

    img_tmp = separate_stains(img_as_ubyte(img), he1_from_rgb)
    
    # The RGB -> H&E transformation maps:
    # black (0,0,0)       |-> (-1.13450710, 0.00727017021)
    # white (255,255,255) |-> (-9.08243792, 0.05.82022531)
    # red   (255,0,0)     |-> (-9.53805685, 6.44503007)
    # green (0,255,0)     |-> (-0.164661728, -5.42507111)
    # blue  (0,0,255)     |-> (-1.64873355, -0.947216369)
    
    if normalize:
        img_tmp[:,:,0] = 2*(img_tmp[:,:,0] + 9.53805685) / (-0.164661728 + 9.53805685) - 1
        img_tmp[img_tmp[:,:,0] < -1.0, 0] = -1.0
        img_tmp[img_tmp[:,:,0] >  1.0, 0] =  1.0
        
        img_tmp[:,:,1] = 2*(img_tmp[:,:,1] + 5.42507111) / (6.44503007 + 5.42507111) - 1
        img_tmp[img_tmp[:,:,1] < -1.0, 0] = -1.0
        img_tmp[img_tmp[:,:,1] >  1.0, 0] =  1.0
        
    return img_tmp[:, :, 0], img_tmp[:, :, 1]

from skimage.segmentation import random_walker
from skimage import morphology
from skimage.filters import sobel

#ihc_rgb = imread(r'C:\Users\griffin\Desktop\MicroDeconvolution\TestingScripts\SamplePics\SK108 3554 28_1_CD3_IL10_Set 1_20X_Take 1.jpg')
ihc_rgb = imread(r'TestImage.jpg')

# Color deconvolution
# Hematoxylin(0) & DAB(2)
rgb_from_hrd = np.array([[0.65, 0.70, 0.29],
                         [0.1, 0.95, 0.95],
                         [0.27, 0.57, 0.78]])
hrd_from_rgb = linalg.inv(rgb_from_hrd)

# Stain space conversion
ihc_hrd = separate_stains(ihc_rgb, hrd_from_rgb)


# Rescale signals so that intensity ranges from 0 to 1
# ihc_hrd[:, :, (0,1, or 2 -- is the color channel)]
def stainspace2array(ihc_xyz, channel):
    rescale = rescale_intensity(ihc_xyz[:, :, channel], out_range=(0,1))
    stain_array = np.dstack((np.zeros_like(rescale), rescale, rescale))
    grey_array = rgb2grey(stain_array)

    return grey_array


DAB_Grey_Array = stainspace2array(ihc_hrd, 2)
Hema_Grey_Array = stainspace2array(ihc_hrd, 0)
red_Grey_Array = stainspace2array(ihc_hrd, 1)

def random_walk_segmentation(input_image, output_folder):

    input_image = imread(input_image)
    ihc_hrd = separate_stains(input_image, hrd_from_rgb)

    DAB_Grey_Array = stainspace_to_2d_array(ihc_hrd, 2)
    Hema_Gray_Array = stainspace_to_2d_array(ihc_hrd, 0)
    GBIred_Gray_Array = stainspace_to_2d_array(ihc_hrd, 1)

    #Perform Random Walker, fills in positive regions
    DAB_segmentation = random_walker(DAB_Grey_Array, get_markers(DAB_Grey_Array, .3, .5), beta=130, mode='cg_mg')
    Hema_segmentation = random_walker(Hema_Gray_Array, get_markers(Hema_Gray_Array, .2, .4), beta=130, mode='cg_mg')
    GBIred_segmentation = random_walker(GBIred_Gray_Array, get_markers(GBIred_Gray_Array, .4, .5), beta=130,
                                        mode='cg_mg')

    '''Compute and Output'''
    #Compute and output percentages of pixels stained by each chromagen
    pic_dimensions = np.shape(DAB_segmentation)  # both arrays same shape
    total_pixels = pic_dimensions[0] * pic_dimensions[1]

    #Change negative pixel values from 1 -> 0, positives 2 -> 1
    subtrahend_array = np.ones_like(DAB_segmentation)
    DAB_segmentation = np.subtract(DAB_segmentation, subtrahend_array)
    Hema_segmentation = np.subtract(Hema_segmentation, subtrahend_array)
    GBIred_segmentation = np.subtract(GBIred_segmentation, subtrahend_array)

    #Count positive pixels
    DAB_pixels = np.count_nonzero(DAB_segmentation)
    Hema_pixels = np.count_nonzero(Hema_segmentation)
    red_pixels = np.count_nonzero(GBIred_segmentation)


    #Percent of image covered by positive staining
    DAB_coverage_percent = (round((DAB_pixels / total_pixels * 100), 1))
    Hema_coverage_percent = (round((Hema_pixels / total_pixels * 100), 1))

    #An overlay of the DAB and Hematoxylin segmented images, for total cellular area
    total_cell_array = np.add(DAB_segmentation, Hema_segmentation)
    #Number of pixels covered by cellular area
    total_cell_pixels = np.count_nonzero(total_cell_array)

    #Percent of image covered by cellular area (DAB OR Hematoxylin)
    total_cell_percent = (round((total_cell_pixels / total_pixels * 100), 1))

    #The percentage of DAB/CD3+ cells out of the total number of cells
    percent_pos_cells = (round((DAB_pixels / total_cell_pixels * 100), 1))


    #The percentage of the image covered by cytokines
    Red_coverage_percent = (round((red_pixels / total_pixels * 100), 1))

    red_plus_total_array = np.add(total_cell_array, GBIred_segmentation)
    red_plus_total_pixels = np.count_nonzero(red_plus_total_array)

    #The percentage of the area covered by cytokines, with non-cellular regions subtracted
    adjusted_red_coverage_percent = (round((red_pixels / red_plus_total_pixels * 100), 1))




    # Plot images
    fig, axes = plt.subplots(2, 2, figsize=(12, 11))
    ax0, ax1, ax2, ax3 = axes.ravel()

    ax0.imshow(input_image, cmap=plt.cm.gray, interpolation='nearest')
    ax0.set_title("Original")

    ax1.imshow(DAB_segmentation, cmap=plt.cm.gray, interpolation='nearest')
    ax1.set_title("DAB")

    ax2.imshow(GBIred_segmentation, cmap=plt.cm.gray)
    ax2.set_title("GBI red")

    ax3.imshow(Hema_segmentation, cmap=plt.cm.gray)
    ax3.set_title("Hematoxylin")

    for ax in axes.ravel():
        ax.axis('off')

    fig.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=None, hspace=None)


    output_filename = 'output' + time.strftime("%Y-%m-%d %H:%M:%S")
    plt.savefig(output_folder + output_filename)

    #do a save csv here, maybe delete return statement after this comment

    return output_filename#, DAB_coverage_percent, Hema_coverage_percent, total_cell_percent, percent_pos_cells, Red_coverage_percent, adjusted_red_coverage_percent

 def test_hdx_rgb_roundtrip(self):
     from skimage.color.colorconv import hdx_from_rgb, rgb_from_hdx
     img_rgb = img_as_float(self.img_rgb)
     conv = combine_stains(separate_stains(img_rgb, hdx_from_rgb),
                           rgb_from_hdx)
     assert_array_almost_equal(conv, img_rgb)

# Rescale signals so that intensity ranges from 0 to 1
# ihc_hrd[:, :, (0,1, or 2 -- is the color channel)]
def stainspace_to_2d_array(ihc_xyz, channel):
    #rescale = rescale_intensity(ihc_xyz[:, :, channel], out_range=(0,1))
    #stain_array = np.dstack((np.zeros_like(rescale), rescale, rescale))

    #try to not reverse engineer rescale right now
    stain_array = ihc_xyz[:, :, channel]
    #plt.imshow(stain_array)
    gray_array = rgb2grey(stain_array)
    #plt.imshow(gray_array)
    return gray_array


#Stain space conversion
source_ihc_hrd = separate_stains(source_ihc_rgb, hrd_from_rgb)
target_ihc_hrd = separate_stains(target_ihc_rgb, hrd_from_rgb)


source_Hema_Gray_Array = stainspace_to_2d_array(source_ihc_hrd, 0)
source_red_Gray_Array = stainspace_to_2d_array(source_ihc_hrd, 1)
source_DAB_Gray_Array = stainspace_to_2d_array(source_ihc_hrd, 2)

target_Hema_Gray_Array = stainspace_to_2d_array(target_ihc_hrd, 0)
target_red_Gray_Array = stainspace_to_2d_array(target_ihc_hrd, 1)
target_DAB_Gray_Array = stainspace_to_2d_array(target_ihc_hrd, 2)


'''
def get_stats(input_1D_array):
    mean = np.mean(input_1D_array)

	if not f.endswith(".jpg") :
		l.remove(f)


qstain = np.array([[.26451728, .5205347, .81183386], [.9199094, .29797825, .25489032], [.28947765, .80015373, .5253158]])




for im in l :

	print im

	A = transform.rescale(io.imread("../data/" + im), 0.25)

	deconv = ski.img_as_float(color.separate_stains(A, np.linalg.inv(qstain)))






	subveins1 = \
	morphology.remove_small_objects(
		filter.threshold_adaptive(
			filter.gaussian_filter(
				deconv[:, :, 2] / deconv[:, :, 0],
				11),
			250, offset = -0.13),
		60)

    # Inflammatory focus count

    qstain = np.array([[.26451728, .5205347, .81183386], [.9199094, .29797825, .25489032], [.28947765, .80015373, .5253158]])

    deconv = ski.img_as_float(color.separate_stains(transform.rescale(A, 0.25), np.linalg.inv(qstain)))



    subveins1 = \
    morphology.remove_small_objects(
        filter.threshold_adaptive(
            filter.gaussian_filter(
                deconv[:, :, 2] / deconv[:, :, 0],
                11),
            250, offset = -0.13),
        60)

    subveins2 = \
    morphology.remove_small_objects(
        filter.threshold_adaptive(

graies = [(rgb2gray(np.average(image[labels == i], axis=0)), i) for i in xrange(cluster_num)]

graies.sort()

nuclei__ = np.zeros((w, h))

nuclei__[labels == graies[0][1]] = 1

image___ = np.copy(image)

image___[labels == graies[1][1]] = 255

image___[labels == graies[2][1]] = 255


hematoxylin_ = separate_stains(image_, conv_matrix)[:, :, 0]

hematoxylin__ = separate_stains(image__, conv_matrix)[:, :, 0]

hematoxylin___ = separate_stains(image___, conv_matrix)[:, :, 0]

# opened_ = opening(hematoxylin_, disk(int(sys.argv[2])))
#
# opened__ = opening(hematoxylin__, disk(int(sys.argv[2])))

cleared = clear_border(nuclei)

cleared_ = clear_border(nuclei_)

cleared__ = clear_border(nuclei__)

# <codecell>

A = io.imread("../data/" + l[5])
io.imshow(A)

# <codecell>

#B = exposure.adjust_sigmoid(filter.gaussian_filter(A[:, :, 0], 19), cutoff=.45, gain=15)



B = exposure.adjust_sigmoid(
    filter.gaussian_filter(
        exposure.rescale_intensity(
            color.separate_stains(
                A, 
                np.linalg.inv(qstain)), 
            out_range=(0, 1))[:, :, 1], 
        29), 
    cutoff=.35, gain=20)
io.imshow(B)

# <codecell>

#b = morphology.remove_small_objects(filter.threshold_adaptive(filter.gaussian_filter(exposure.adjust_sigmoid(A[:,:,1]), 31), 501, offset=-0.05), 2000)

C = morphology.remove_small_objects(
    filter.threshold_adaptive(B, 301, offset=-0.025), 
    4000)
#io.imshow(morphology.binary_closing(np.logical_or(morphology.binary_dilation(C, morphology.disk(11)), b), morphology.disk(31)))
io.imshow(C)

def RunScript():

    # Color deconvolution
    # Normalized optical density matrix
    # see Ruifrok AC, Johnston DA. Quantification of histological staining by color deconvolution.
    # R      G      B
    # X      X      X   Hematoxylin(0)
    # X      X      X   Red(1)
    # X      X      X   DAB(2)
    # Hematoxylin(0), Red(1), DAB(2)
    rgb_from_hrd = np.array([[0.644, 0.710, 0.285],
                             [0.0326, 0.873, 0.487],
                             [0.270, 0.562, 0.781]])
    # conv_matrix
    hrd_from_rgb = linalg.inv(rgb_from_hrd)

    # Import picture
    #ihc_rgb = imread(r'TestImage.jpg')
    ihc_rgb = imread(r'TimedRunImage.jpg')

    # Rescale signals so that intensity ranges from 0 to 1
    # ihc_hrd[:, :, (0,1, or 2 -- is the color channel)]
    def stainspace_to_2d_array(ihc_xyz, channel):
        rescale = rescale_intensity(ihc_xyz[:, :, channel], out_range=(0, 1))
        stain_array = np.dstack((np.zeros_like(rescale), rescale, rescale))
        grey_array = rgb2grey(stain_array)

        return grey_array

    # Stain space conversion
    ihc_hrd = separate_stains(ihc_rgb, hrd_from_rgb)

    DAB_Grey_Array = stainspace_to_2d_array(ihc_hrd, 2)
    Hema_Gray_Array = stainspace_to_2d_array(ihc_hrd, 0)
    permRed_Gray_Array = stainspace_to_2d_array(ihc_hrd, 1)

    # Get markers for random walk
    def get_markers(grey_array, bottom_thresh, top_thresh):
        markers = np.zeros_like(grey_array)
        markers[grey_array < bottom_thresh] = 1
        markers[grey_array > top_thresh] = 2

        return markers

    # perform Random Walker, fills in positive regions
    DAB_segmentation = random_walker(DAB_Grey_Array, get_markers(DAB_Grey_Array, .3, .5), beta=130, mode='cg')
    Hema_segmentation = random_walker(Hema_Gray_Array, get_markers(Hema_Gray_Array, .2, .4), beta=130, mode='cg')
    permRed_segmentation = random_walker(permRed_Gray_Array, get_markers(permRed_Gray_Array, .4, .5), beta=130, mode='cg')

    """PRINTING OUTPUT"""

    print(20 *'-')
    print(' ')

    '''Compute and Output'''
    # Compute and output percentages of pixels stained by each chromagen
    pic_dimensions = np.shape(DAB_segmentation)  # both arrays same shape
    total_pixels = pic_dimensions[0] * pic_dimensions[1]

    # change negative pixel values from 1 -> 0, positives 2 -> 1
    subtrahend_array = np.ones_like(DAB_segmentation)
    DAB_segmentation = np.subtract(DAB_segmentation, subtrahend_array)
    Hema_segmentation = np.subtract(Hema_segmentation, subtrahend_array)
    permRed_segmentation = np.subtract(permRed_segmentation, subtrahend_array)

    # count positive pixels
    DAB_pixels = np.count_nonzero(DAB_segmentation)
    Hema_pixels = np.count_nonzero(Hema_segmentation)
    red_pixels = np.count_nonzero(permRed_segmentation)

    DAB_coverage_percent = (round((DAB_pixels / total_pixels * 100), 1))
    print("The percentage of the image covered by DAB is: " + str(DAB_coverage_percent) + "%")

    Hema_coverage_percent = (round((Hema_pixels / total_pixels * 100), 1))
    print("The percentage of the image covered by Hematoxylin is: " + str(Hema_coverage_percent) + "%")

    total_cell_array = np.add(DAB_segmentation, Hema_segmentation)
    total_cell_pixels = np.count_nonzero(total_cell_array)

    total_cell_percent = (round((total_cell_pixels / total_pixels * 100), 1))
    print("The percentage of the image covered by DAB & Hematoxylin is: " + str(total_cell_percent) + "%")

    percent_pos_cells = (round((DAB_pixels / total_cell_pixels * 100), 1))
    print("The percentage of CD3+ cells out of the total number of cells is: " + str(percent_pos_cells) + "%")

    PercentPos = percent_pos_cells

    """
    if PercentPos >= 81:
        print('Proportion Score: 5')
        proportion_score = 5
    elif PercentPos >= 61:
        print('Proportion Score: 4')
        proportion_score = 4
    elif PercentPos >= 41:
        print('Proportion Score: 3')
        proportion_score = 3
    elif PercentPos >= 21:
        print('Proportion Score: 2')
        proportion_score = 2
#.........这里部分代码省略.........

#Color deconvolution
#DAB and perm red(1)
rgb_from_drx = np.array([[0.270, 0.562, 0.781],
                         [0.0326, 0.873, 0.487],
                         [0.0, 0.0, 0.0]])
rgb_from_drx[2, :] = np.cross(rgb_from_drx[0, :], rgb_from_drx[1, :])
drx_from_rgb = linalg.inv(rgb_from_drx)

#Import picture
ihc_rgb = data.imread(r'TestImage.jpg')


#Stain space conversion
#ihc_hax = separate_stains(ihc_rgb, hax_from_rgb)
ihc_drx = separate_stains(ihc_rgb, drx_from_rgb)


#Rescale signals? - might help, might not
#[:, :, 012 color]
permred_rescale = rescale_intensity(ihc_drx[:, :, 1], out_range=(0, 1))
permred_array = np.dstack((np.zeros_like(permred_rescale), permred_rescale, permred_rescale))


#Blob detection
image2d = rgb2grey(permred_array)


blobs_dog = blob_dog(image2d, min_sigma=1, max_sigma=40, threshold=.5, overlap=0.1)
blobs_dog[:, 2] = blobs_dog[:, 2] * sqrt(2)

 def setUp(self):
     self.rgb = misc.imread(
         os.path.join(curdir, "../immunopy/image/hdab256.tif"))
     self.hdx = color.separate_stains(self.rgb, color.hdx_from_rgb)
     self.hem = self.hdx[:,:,0]
     self.dab = self.hdx[:,:,1]