I’m using EMAN 2.22 (a convolutional neural network) to extract the inner membrane of a cell from tomograms. The output from EMAN 2.22 is quite noisy however. I’m using Otsu’s threshold to get rid of the salt and pepper noise and then converting the remainder to binary 0 or 1 (because parts which are not in the membrane are sometimes brighter than those which are). The noise is quite large clusters and the membrane is not one large connected cluster either. I need to convert the membrane output to a single line model to do further analysis but to do this I need to get rid of the noise.
I’ve tried morphological transformations using scikit-learn. Orig is after applying Otsu and converting to binary.
The best I have (to try and increase the joining of the real membrane but reducing connection to clusters) is by dilating, closing and then skeletonizing. (This is the third image in this post that has regression mapped on top).
I was doing this to try to increase the size of membrane clusters and reduce noise clusters to try to do a connected component analysis (and remove smaller clusters). My membrane and noise seem to be in similarly sized clusters, however.
(I also tried to do an active contour on the whole image, given the noise is on the inner of the membrane, but I can’t seem to get it to work well.)
I’ve split the image into tiles and then my thought was to do a Hough transform. I’ve given up on because my lines aren’t necessarily straight and it becomes harder and harder as I reduce my tile size to know where the true membrane is. I moved onto linear regression (on the dilated, closed, skeletonized tiles) to get a line through my membrane which works quite well but is still skewed by large noise.
a = blockshaped(orig, 256, 256)
fig, axs = plt.subplots(4,4, figsize=(10, 10), sharex="col", sharey="row")
fig.subplots_adjust(hspace = .5, wspace=.001)
axs = axs.ravel()
for i in range(len(a)):
a[i] = dilation(a[i])
a[i] = closing(a[i])
a[i] = skeletonize(a[i])
axs[i].imshow(a[i])
i += 1
x, y = np.where(a[2]>0)
plt.scatter(x, y)
plt.show()
from numpy import polyfit
from numpy import polyval
p1 = polyfit(x, y, 1)
p2 = polyfit(x, y, 2)
p3 = polyfit(x, y, 3)
print p1
%matplotlib
plt.scatter(x, y)
xp = np.linspace(0, 256, 100)
plt.plot(xp, polyval(p1, xp), "r-")
plt.plot(xp, polyval(p2, xp), "b--")
plt.plot(xp, polyval(p3, xp), "m:")
I’ve considered gradient descent algorithms like Adam to find a better function but this seems a bit too extreme?
Another thought was to make a knn graph and collapse branches leaving the main branch but I have no idea how to do this.
I’m sure that I’m missing a really simple thing and any advice which could point me in the right direction of another technique to consider would be fantastic.
Thank you so much!
