fakeimage = np.zeros((128,128))
fakeimage[31:96,31:96] = 10.0

"""
Explanation: Create a fake image with a step change:
End of explanation
"""


pl.subplot(111)
pl.imshow(fakeimage)
pl.show()

"""
Explanation: Plot the fake shape to see how it looks:
End of explanation
"""


fakeimage+=np.random.normal(0,0.02,size=(128,128))

"""
Explanation: Add some Gaussian noise. I'm adding noise with a mean of zero and a standard deviation of 0.01. This function creates one realisation of this noise and adds it to our fake image:
End of explanation
"""


pl.subplot(111)
pl.imshow(fakeimage)
pl.imsave("image.png",fakeimage)
pl.show()

"""
Explanation: Plot the noisy shape to see how it looks. (The noise is small compared to the amplitude of the structure so it doesn't look very different)
End of explanation
"""


img= cv2.imread('image.png',0)

# Use Otsu filtering to determine thresholds for finding the edge:
high_thresh, thresh_im = cv2.threshold(img, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
low_thresh = 0.9*high_thresh

# find the edges using Canny method:
edges = cv2.Canny(img,low_thresh,high_thresh)

# plot the output:
pl.plot(1),pl.imshow(edges,cmap = 'gray')
pl.show()

"""
Explanation: Then we detect an edge. 
I'm using Canny edge detection to find the edge of the shadow: https://docs.opencv.org/3.3.1/da/d22/tutorial_py_canny.html
I haven't tuned it very finely, so it's not particularly brilliant... it sounds like you algorithm works better, so you should replace this bit.
End of explanation
"""


n_real = 100 # number of realisations

for i in range(0,n_real):
    
    # same underlying pattern:
    fakeimage = np.zeros((128,128))
    fakeimage[31:96,31:96] = 10.0
    
    # add a new noise realisation:
    fakeimage+=np.random.normal(0,0.02,size=(128,128))
    
    # save it to an image:
    pl.imsave("image_"+str(i)+".png",fakeimage)

"""
Explanation: To make multiple noise realisations we need to make a loop. This will output $N_{\rm real}$ images, each with the same underlying pattern but a different realisation of noise from the same Gaussian distribution.
End of explanation
"""


all_edges = np.zeros(edges.shape)
for i in range(0,n_real):

    img= cv2.imread('image_'+str(i)+'.png',0)

    # Use Otsu filtering to determine thresholds for finding the edge:
    high_thresh, thresh_im = cv2.threshold(img, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
    low_thresh = 0.9*high_thresh

    # find the edges using Canny method:
    edges = cv2.Canny(img,low_thresh,high_thresh)
    
    all_edges+=edges

"""
Explanation: Then we fit the edges for each of the realisations sparately and combine the output:
End of explanation
"""


# plot the output:
pl.plot(1),pl.imshow(all_edges,cmap = 'gray')
pl.show()

"""
Explanation: The combined output shows the distribution of fitted values:
End of explanation
"""


