I am trying to run some statistics on the contour outline of imaged zooplankton specimens from the LOKI system (https://www.researchgate.net/publication/236931182_Imaging_of_plankton_specimens_with_the_lightframe_on-sight_keyspecies_investigation_LOKI_system ).
The LOKI system extracts Regions Of Interest (ROI) directly in the underwater unit. Thus, each file generally contains the image of an individual specimen. It can be generally assumed, that the largest object in the image is the ROI to evaluate.
To get the contour information I use some functions of the EBImage package (see some sample code below). This works fine for most images.
Nevertheless, in some images the background is heterogeneous and a little brighter. Thus, the background may partially contribute as “false positive” to the ROI. In other cases parts of the background become structures of their own, being sometimes even larger than the ROI (see the two upper rows of the attached panel imge).
I already tried a couple of things like Otsu’s thresholding method or the filter2 functions. I also tried to first subtract a certain value and finally spread the range (multiplication or power function). So far I was not able to find a universal and simple approach.
The ROI is generally part of the brightest structures in the image. But specimen appendages may vary in brightness and be even close to the background. As these structures are also important the background cannot be simply darkened, as it would also loose some of these structures (see last row of the panel image).
Thus, an interesting approach might be to start from the brightest structures and to identify pixels above a local, adaptive threshold; something like a modified Hoshen-Koppelman algorithm.
Now is my question, whether here is someone who may give a hint, an approach or even some lines of code.
As I finally want to process >100k images it is of course advantageous to keep the solution as (computational) cheap as possible. However, to have a good solution is the preference.
Below I attached some quick code to get an idea.
Also included are some images for testing…
Best regards,
Jan
# short test on selected contour outlines of zooplankton specimens
# Dr. Jan Schulz, ICBM, University of Oldenburg
# data retrieved from a LOKI device (DOI: 10.2971/jeos.2010.10017s)
# file date: 30.09.2019
library ("EBImage")
#-------------------------------------------------------------------------------
# create a small fake matrix easily to distinguish from real images
# and to ensure proper functionality when an images could not be loaded correctly
# in subsequent image processing it will be used to emulate a smal 7*7 pixel image,
# with black background and a 3*3 pixels white dot in its centre
fakeMatrix = matrix (0, nrow = 7, ncol = 7)
fakeMatrix [3:5, 3:5] = 1
# initialise the greyIMG matrix with the tiny fake matrix
# just for security reasons, when import fails the routines calculate a very small
# image that can be identified pretty easily
greyIMG = fakeMatrix
aFileList = list ( "D:/Daten/HE434Test/0067_HE 434-067/Haul 1/LOKI_10001.01/Pictures/2014.10.20 11 25/20141020 112511 943 000000 2037 0102.png",
"D:/Daten/HE434Test/0067_HE 434-067/Haul 1/LOKI_10001.01/Pictures/2014.10.20 11 24/20141020 112403 993 000000 2073 0456.bmp",
"D:/Daten/HE434Test/0067_HE 434-067/Haul 1/LOKI_10001.01/Pictures/2014.10.20 11 24/20141020 112434 218 000000 0108 0933.bmp")
# create a file for the plot output
png (paste ("D:/Daten/HE434Test/0067_HE 434-067/Haul 1/Test", ".png"),
width = 21,
height = 29.7,
units = 'cm',
res = 310)
par (mfrow = c (3, 2), # display different plots as an 3x2 array
oma = c (4, 2, 3, 2), # define outer margins as lines of text
mar = c (5, 4, 2, 2))
for (i in 1:length (aFileList)) {
aFileFullPath = aFileList [[i]] [1]
# handle BMP images
if (toupper (tools::file_ext (aFileFullPath)) == "BMP")
{
#read the BMP image
aSourceImage = read.bmp (aFileFullPath)
# use the number of colour in the attributes to re-scale values to [0..1]
aSourceImage = aSourceImage [ , ] / max (attr (aSourceImage, "header")$ncolors)
greyIMG = channel (aSourceImage, "luminance")
}
#handle PNG images
if (toupper (tools::file_ext (aFileFullPath)) == "PNG")
{
# read the PNG file
aSourceImage = readPNG (aFileFullPath)
aSourceImage = channel (aSourceImage, "luminance")
greyIMG = aSourceImage [ , , 1]
}
#-------------------------- improve image and create a binary silhouette image
# do some image enhancing
# adding a value -> increase brightness
# multiplication -> increasing contrast
# power -> gamma correction
# greater as -> filter for pixels being bright enough
if (toupper (tools::file_ext (aFileFullPath)) == "BMP")
{
#enhancedgreyIMG = greyIMG [ , ] ^ 0.8
enhancedgreyIMG = (greyIMG [ , ]- min (greyIMG))*8
}
if (toupper (tools::file_ext (aFileFullPath)) == "PNG")
{
#enhancedgreyIMG = ((greyIMG [ , ] + 0.007) * 16) ^ 4
#enhancedgreyIMG = greyIMG [ , ] ^ 0.8
enhancedgreyIMG = ((greyIMG [ , ] - min (greyIMG))*8 )
}
# now find the objects in the image
IMGmask = fillHull (enhancedgreyIMG)
# use a small disk to close gaps in the structures and
# finally create the primary binary matrix from values still >0.3
# this matrix contains all objects, including dizzy background objects
IMGmask = closing (IMGmask, makeBrush (9, shape = 'diamond')) > 0.3
# now find all individual objects being not connected and
# give them an individual ID, thus IMGmask becomes an integer matrix
IMGmask = bwlabel (IMGmask)
# identify contours of all the objects having an ID
IMGcontour = ocontour (IMGmask)
# now get the ID of longest contour
# It is assumed that the largest object in the image is the
# focus object extracted by the LOKI system and we expect
# that this is our object to investigate
IDofLongestContour = which.max (lengths (IMGcontour))
# now we switch back to binary mask
# just containing the largest object
IMGmask = (IMGmask == IDofLongestContour)
# get the dimension of the image for convenient reasons
tmp_IMGwidth = dim (greyIMG) [2]
tmp_IMGheigth = dim (greyIMG) [1]
# create teh first plot with the original image
plot (x = c (1, tmp_IMGwidth),
y = c (1, tmp_IMGheigth),
main = "Original image (distorted aspect ratio)",
xlab = "Horizontal pixel",
ylab = "Vertical pixels",
type = "n")
rasterImage (aSourceImage/max (aSourceImage) , 1, 1, tmp_IMGwidth, tmp_IMGheigth)
plot (x = c (1, tmp_IMGwidth),
y = c (1, tmp_IMGheigth),
main = "Processed image and extracted contour",
xlab = "Horizontal pixel",
ylab = "Vertical pixels",
type = "n")
rasterImage (enhancedgreyIMG/max (enhancedgreyIMG) , 1, 1, tmp_IMGwidth, tmp_IMGheigth)
points ( IMGcontour[[IDofLongestContour]] [ , 2] ,
-IMGcontour[[IDofLongestContour]] [ , 1] + tmp_IMGheigth,
pch = 20,
cex = 0.1,
col = "red")
}
dev.off ()