Image Processing: Segmentation of Optic Disc

Question

I'm new in image processing. I work in my project now to segment optic disc in retinal images. I did the ROI selection, red channel extraction, applying filter, thresholding, and then morphological operations. Here is the code and the input image:

ROI

from numpy.ma.core import ndim
def getROI(image):
    b,g,r = cv2.split(image)
    b = cv2.Blur(b,(15,15),0)
    kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE,(15,15))
    b = ndimage.grey_opening(b,structure=kernel)    
    (minVal, maxVal, minLoc, maxLoc) = cv2.minMaxLoc(b)

    x0 = int(maxLoc[0])-150
    y0 = int(maxLoc[1])-150
    x1 = int(maxLoc[0])+150
    y1 = int(maxLoc[1])+150
    
    return image[y0:y1,x0:x1]

plt.imshow(image)
plt.show()

RED CHANNEL EXTRACTION

def rgb2Red(img):
    b,g,r = cv2.split(img)
    return r
def rgb2Gray(img):
    return cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)

APPLYING FILTER

def clahe_image(image, clipLimit = 1.0, channels = 'a'):
    old_image = cv2.cvtColor(image, cv2.COLOR_BGR2RGB)
    R,G,B = cv2.split(old_image)
    
    clahe = cv2.createCLAHE(clipLimit=clipLimit, tileGridSize=(8,8))
    
    if channels is 'g': 
        G = clahe.apply(G)
    if channels is 'b': 
        B = clahe.apply(B) 
    if channels is 'r': 
        R = clahe.apply(R) 
    if channels is 'a':
        G = clahe.apply(G)
        B = clahe.apply(B)
        R = clahe.apply(R)

    clahe_image = cv2.merge((R, G, B))

    return clahe_image

def preprocess(image):
    gray_blur = cv2.GaussianBlur(image, (3,3), 0);
    gray = cv2.addWeighted(image, 1.5, gray_blur, -0.5, 0, image);
    kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE,(31,31));
    gray = ndimage.grey_closing(gray,structure=kernel);
    gray = cv2.equalizeHist(gray);
    gray = clahe_image(gray);
    return gray

THRESHOLDING

def optimal_threshold(h,t):
    h1 = h[:t]
    h2 = h[t:]
    m1 = (h1*np.arange(0,t)).sum()/h1.sum()
    m2 = (h2*np.arange(t,len(h))).sum()/h2.sum()
    t2 = int(np.round((m1+m2)/2))
    print(m1,m2,t2)
    if( t2 != t ) : return optimal_threshold(h,t2)
    return t2

h,bins = np.histogram(image,range(257))
t = optimal_threshold(h,255)
print(t)

ret, final_img = cv2.threshold(image,250,255,cv2.THRESH_BINARY)

MORPHOLOGICAL OPERATIONS

img = final_img
kernel = np.ones((3,3), np.uint8)  
img_erosion = cv2.erode(img, kernel, iterations=1)  
img_dilation = cv2.dilate(img, kernel, iterations=4)  
plt.imshow(img_dilation, cmap='gray')   

The result I get:

I want to make the result as perfect as the ground truth and the gray one like this:

Can anyone give a tips how to get optic disc segmented aa perfect as the ground truth? And what if i add vessel removal? Is it okay?