def test_compare_with_cmorph_erode():
    # compare the result of maximum filter with erode

    image = (np.random.random((100, 100)) * 256).astype(np.uint8)
    out = np.empty_like(image)
    mask = np.ones(image.shape, dtype=np.uint8)

    for r in range(1, 20, 1):
        elem = np.ones((r, r), dtype=np.uint8)
        rank.minimum(image=image, selem=elem, out=out, mask=mask)
        cm = cmorph.erode(image=image, selem=elem)
        assert_array_equal(out, cm)

def test_percentile_min():
    # check that percentile p0 = 0 is identical to local min
    img = data.camera()
    img16 = img.astype(np.uint16)
    selem = disk(15)
    # check for 8bit
    img_p0 = rank.percentile(img, selem=selem, p0=0)
    img_min = rank.minimum(img, selem=selem)
    assert_array_equal(img_p0, img_min)
    # check for 16bit
    img_p0 = rank.percentile(img16, selem=selem, p0=0)
    img_min = rank.minimum(img16, selem=selem)
    assert_array_equal(img_p0, img_min)

def test_16bit():
    image = np.zeros((21, 21), dtype=np.uint16)
    selem = np.ones((3, 3), dtype=np.uint8)

    for bitdepth in range(17):
        value = 2 ** bitdepth - 1
        image[10, 10] = value
        assert rank.minimum(image, selem)[10, 10] == 0
        assert rank.maximum(image, selem)[10, 10] == value
        assert rank.mean(image, selem)[10, 10] == int(value / selem.size)

def test_smallest_selem16():
    # check that min, max and mean returns identity if structuring element
    # contains only central pixel

    image = np.zeros((5, 5), dtype=np.uint16)
    out = np.zeros_like(image)
    mask = np.ones_like(image, dtype=np.uint8)
    image[2, 2] = 255
    image[2, 3] = 128
    image[1, 2] = 16

    elem = np.array([[1]], dtype=np.uint8)
    rank.mean(image=image, selem=elem, out=out, mask=mask,
              shift_x=0, shift_y=0)
    assert_array_equal(image, out)
    rank.minimum(image=image, selem=elem, out=out, mask=mask,
                 shift_x=0, shift_y=0)
    assert_array_equal(image, out)
    rank.maximum(image=image, selem=elem, out=out, mask=mask,
                 shift_x=0, shift_y=0)
    assert_array_equal(image, out)

    def run(self, workspace):
 
        image = workspace.image_set.get_image(self.image_name.value)
        nuclei_image = image.pixel_data[:,:]
        image_collection = []
      
        #
        #Get the global Threshold with Otsu algorithm and smooth nuclei image
        #
#         nuclei_smoothed = self.smooth_image(image_collection[3][0], image.mask, 1)

        global_threshold_nuclei = otsu(nuclei_image, min_threshold=0, max_threshold=1)
        print "the threshold compute by the Otsu algorythm is %f" % global_threshold_nuclei      
        
        
        #
        #Binary thee "DAPI" Image (Nuclei) and labelelize the nuclei
        #

        binary_nuclei = (nuclei_image >= global_threshold_nuclei)
        labeled_nuclei, object_count = scipy.ndimage.label(binary_nuclei, np.ones((3,3), bool))
        print "the image got %d detected" % object_count
        
        #
        #Fill the hole and delete object witch touch the border. 
        #labeled_nuclei is modify after the function
        #Filter small object and split object
        #
        labeled_nuclei = fill_labeled_holes(labeled_nuclei)        
        labeled_nuclei = self.filter_on_border(labeled_nuclei)
        labeled_nuclei = self.filter_on_size(labeled_nuclei, object_count)         
        labeled_nuclei = self.split_object(labeled_nuclei)
        
        #
        #Edge detection of nuclei image and object are more separated 
        #
        labeled_nuclei_canny = skf.sobel(labeled_nuclei)       
        labeled_nuclei[labeled_nuclei_canny > 0] = 0
        labeled_nuclei = skr.minimum(labeled_nuclei.astype(np.uint16), skm.disk(3))
        
        image_collection.append((nuclei_image, "Original"))
        image_collection.append((labeled_nuclei, "Labelized image"))
        workspace.display_data.image_collection = image_collection

        #
        #Create a new object which will be add to the workspace
        #
        objects = cellprofiler.objects.Objects()
        objects.segmented = labeled_nuclei
        objects.parent_image = nuclei_image
        
        workspace.object_set.add_objects(objects, self.object_name.value)

def test_empty_selem():
    # check that min, max and mean returns zeros if structuring element is empty

    image = np.zeros((5, 5), dtype=np.uint16)
    out = np.zeros_like(image)
    mask = np.ones_like(image, dtype=np.uint8)
    res = np.zeros_like(image)
    image[2, 2] = 255
    image[2, 3] = 128
    image[1, 2] = 16

    elem = np.array([[0, 0, 0], [0, 0, 0]], dtype=np.uint8)

    rank.mean(image=image, selem=elem, out=out, mask=mask,
              shift_x=0, shift_y=0)
    assert_array_equal(res, out)
    rank.minimum(image=image, selem=elem, out=out, mask=mask,
                 shift_x=0, shift_y=0)
    assert_array_equal(res, out)
    rank.maximum(image=image, selem=elem, out=out, mask=mask,
                 shift_x=0, shift_y=0)
    assert_array_equal(res, out)

.. note::

    `skimage.dilate` and `skimage.erode` are equivalent filters (see below for
    comparison).

Here is an example of the classical morphological greylevel filters: opening,
closing and morphological gradient.

"""

from skimage.filter.rank import maximum, minimum, gradient

ima = data.camera()

closing = maximum(minimum(ima, disk(5)), disk(5))
opening = minimum(maximum(ima, disk(5)), disk(5))
grad = gradient(ima, disk(5))

# display results
fig = plt.figure(figsize=[10, 7])
plt.subplot(2, 2, 1)
plt.imshow(ima, cmap=plt.cm.gray)
plt.xlabel('original')
plt.subplot(2, 2, 2)
plt.imshow(closing, cmap=plt.cm.gray)
plt.xlabel('greylevel closing')
plt.subplot(2, 2, 3)
plt.imshow(opening, cmap=plt.cm.gray)
plt.xlabel('greylevel opening')
plt.subplot(2, 2, 4)

.. note::

    `skimage.dilate` and `skimage.erode` are equivalent filters (see below for
    comparison).

Here is an example of the classical morphological gray-level filters: opening,
closing and morphological gradient.

"""

from skimage.filter.rank import maximum, minimum, gradient

noisy_image = img_as_ubyte(data.camera())

closing = maximum(minimum(noisy_image, disk(5)), disk(5))
opening = minimum(maximum(noisy_image, disk(5)), disk(5))
grad = gradient(noisy_image, disk(5))

# display results
fig = plt.figure(figsize=[10, 7])

plt.subplot(2, 2, 1)
plt.imshow(noisy_image, cmap=plt.cm.gray)
plt.title('Original')
plt.axis('off')

plt.subplot(2, 2, 2)
plt.imshow(closing, cmap=plt.cm.gray)
plt.title('Gray-level closing')
plt.axis('off')

#.........这里部分代码省略.........
            
   if infdef_tuple:
      inf_idx = [i for i, v in enumerate(cropped_contour[:,0,:]) if v[0] == infdef_tuple[0] and v[1] == infdef_tuple[1]]
      if inf_idx:
         if right_side:
            cropped_contour = cropped_contour[:inf_idx[0],:,:]
         else:
            cropped_contour = cropped_contour[inf_idx[0]:,:,:]
         
   if right_side:
      cropped_contour = cropped_contour[cropped_contour[:,0,1] != 1]
   else:
      cropped_contour = cropped_contour[cropped_contour[:,0,0] != 1]

    # Draw the cropped contour
    #cv2.drawContours(imarraymarkup,cropped_contour,-1,(255,255,0),30)
    #cv2.drawContours(imarraymarkup,biggest_contour,-1,(255,0,0),30)

    # Fill in the cropped polygon to mask
    #cv2.fillPoly(maskarray, pts = [cropped_contour], color=(255,255,255))
   cv2.fillPoly(maskarray, pts = [cropped_contour], color=(255,255,255))
    #maskarray = ~np.all(maskarray == 0, axis=1)a
    #print maskarray
    #maskarray[~np.all(maskarray == 0, axis=2)]

    # Threshold the cropped version
    #ret_otsu,thresh_otsu = cv2.threshold(imarray,0,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU)

    # Multiply original image to the mask to get the cropped image
   imarray2 = imarray + onesarray
   imcrop_orig = cv2.bitwise_and(imarray2, maskarray)
    #cv2.drawContours(maskarray,biggest_contour,-1,(255,0,0),3)
    
   rankmin_out = rank.minimum(imcrop_orig, disk(20))
   thresh = mahotas.otsu(rankmin_out, ignore_zeros = True)

    # Draw thick black contour to eliminate the skin and nipple from the image
   cv2.drawContours(rankmin_out,cropped_contour,-1,(0,0,0),255) # 
    #cv2.drawContours(maskarray,cropped_contour,-1,(0,0,0),255) # 

    # Apply the thresholding to generate a new matrix and convert to int type
   otsubool_out = rankmin_out > thresh
   otsubinary_out = otsubool_out.astype('uint8')
   otsuint_out = otsubinary_out * 255

    # Crop out the fibroglandular tissue
    #print output2.shape, imarray2.shape, np.amax(output2), np.amax(imarray2), np.amax(maskarray)
   imcrop_fgt = cv2.bitwise_and(imarray2, otsuint_out) # both arrays are uint8 type

   segmented = maskarray > 0
   segmented = segmented.astype(int)
   segmented_sum = segmented.sum()
   otsubinary_sum = otsubinary_out.sum()
   
   density = (otsubinary_sum*100/segmented_sum).astype(int)
   
   if density < 25:
      dcat = 'Fatty'
   elif density < 50:
      dcat = 'Scattered'
   elif density < 75:
      dcat = 'Heterogenous'
   else:
      dcat = 'Extremely Dense'
        
   if right_side: