AIPreciseAgri_analysis

Soil Segmentation Code Analysis

It is the first step to port the application on the ZCU102: learning what the algorithm does.

Index

• YCbCr_test
• HSV_test
• Histogram_test
• Contour_test
• Otsu_segmentation_test
• Segmentation_test
• Soil segmentation
• Timing analysis

Behavioural analysis

YCbCr_test.py

This script converts the input image in YCbCr color space. Also it shows how to use view_as_block of skimage.util. Block views can be incredibly useful when one wants to perform local operations on non-overlapping image patches. Skimage is a collection of algorithms for image processing.

This script:

1. opens an image with imread;
2. creates a block view of the image in blocks of 32x32x1. Be careful that smaller are the blocks and bigger and larger the resulting image will be but, on the other hand, bigger are the blocks, smaller the resulting image will be and with a worse resolution.

   view = view_as_blocks(image, (32, 32, 1))
   # print(view.shape) ------- (114, 171, 3, 32, 32, 1) that is a concatenation (new_image, shape)

   For example, if the rgb image is of 3648x5472x3 px the view_as_block image will be of 114x171x3: 3648/32=114, 5472/32=171, 3/1=3.
3. Then the image is flattened:

   flatten_view = view.reshape(view.shape[0], view.shape[1], view.shape[2], -1)  # -1 one shape dimension
   # print(flatten_view.shape) ------- (114, 171, 3, 1024) (32, 32, 1) are collapsed in 1024

4. Then, the algorithm starts. It consists of three main parts: mean(), max() and median() functions of numpy library which compute the arithmetic mean, get max value and comput the median along the specified axis.

   mean_view = np.mean(flatten_view, axis=3)
   max_view = np.max(flatten_view, axis=3)
   median_view = np.median(flatten_view, axis=3)
   # they are arrays of 114x171x3 pixels. Each pixel contains the mean, max and median value of a block respectively
   # np.mean()

5. Finally the images (mean_view, max_view and median_view) are plotted in different channels. In particular the image is converted in YCbCr color space using rgb2ycbcr from skimage.color:
   • Y channel
   • Cb channel
   • Cr channel
   • RGB channel

HSV_test.py

It does the same thing of the program above but instead of plotting in YCbCr color profile, it plots now the images in HSV color profile. All the steps are similar.

Histogram_test.py

Also in this script the steps are similar:

1. image reading with imread;

2. creates a block view of the image in blocks of 32x32x1:

   view = view_as_blocks(image, (32, 32, 1))
   # print(view.shape) ------- (114, 171, 3, 32, 32, 1) that is a concatenation (new_image, shape)

3. Then the image is flattened:

   flatten_view = view.reshape(view.shape[0], view.shape[1], view.shape[2], -1)  # -1 one shape dimension
   # print(flatten_view.shape) ------- (114, 171, 3, 1024) (32, 32, 1) are collapsed in 1024

4. mean(), max() and median() are computed.

   mean_view = np.mean(flatten_view, axis=3)
   max_view = np.max(flatten_view, axis=3)
   median_view = np.median(flatten_view, axis=3)
   # they are arrays of 114x171x3 pixels. Each pixel contains the mean, max and median value of a block respectively
   # np.mean()

5. The image is converted in HSV profile.

6. The contrast stretching is applied on value_img.

   # Contrast stretching
   img = value_img
   p2, p98 = np.percentile(img, (2, 98))
   img_rescale = exposure.rescale_intensity(img, in_range=(p2, p98))

   • percentile: compute the q-th percentile of the data along the specified axis and returns the q-th percentile(s) of the array elements. The default axis is to compute the percentile(s) along a flattened version of the array.
   • rescale_intensity: rescale the intensity of the image starting from percentile.

   Now, the contrast of the image is improved, but it isn't the best.

7. The image is equalized using histogram equalization method. Note that the bins are the equal parts in which the histogram is divided for equalization

   # Equalization
   img_eq = exposure.equalize_hist(img) # it returns the image after histogram equalization

8. The image is equalized using adaptive equalization, Contrast Limited Adaptive Histogram Equalization. An algorithm for local contrast enhancement, that uses several histograms computed over different tile regions of the image. It is therefore suitable for improving the local contrast and enhancing the definitions of edges in each region of an image.

   # Adaptive Equalization
   img_adapteq = exposure.equalize_adapthist(img, clip_limit=0.03) # clip_limit (higher values give more contrast). 

9. The results are displayed:

   • Low contrast image and its histogram
   • Contrast stretching image and its histogram
   • Histogram equalization image and its histogram
   • Adaptive equalization image and its histogram

Contour_test.py

This script is similar to the Histogram_test.py:

1. image is read, divided into blocks (32x32x1) and flattened.
2. the median view is computed and it is converted into HSV color space.
3. Then contrast stretching is applied and equalizations are computed similar to Histogram_test.py but, instead of using the brightness profile color, it is used the saturation profile color.
4. The results are displayed:
   • Low contrast image and its histogram
   • Contrast stretching image and its histogram
   • Histogram equalization image and its histogram
   • Adaptive equalization image and its histogram
5. Then the gaussian filter is applied to the contrast stretched image:

   # Applying gaussian filter, with sigma=3 (standard deviation)
   img = gaussian(img_rescale, 3)

6. Find contours of the image:

   # Find contours at a constant value of 0.8
   contours = measure.find_contours(img, 0.8)

   Returns a contourslist of (n,2)-ndarrays. Each contour is an ndarray of shape (n, 2), consisting of n (row, column) coordinates along the contour.

Otsu_segmentation_test.py

This script implements the test of the otsu segmentation algorithm:

1. the image is read, divided into blocks (32x32x1) and flattened.

2. the median view is computed and it is converted into YCbCr color space.

3. the threshold_otsu() method is applied to the cr_img:

   # Threshold value based on Otsu's method
   image = cr_img.copy()
   val = filters.threshold_otsu(image)

4. the flattened cr image is resized using resize function from skimage:

   # resize the flattened image to match the original size
   image = resize(image, (_image_.shape[0], _image_.shape[1]))

5. then the pixels of the image which have values greater than the threshold val are set to 255

     # if the image[i] < val ---> mask_color[i]=true
     mask_color = image < val 
     mask = image < val
     # mask and mask_color are binary image
     _image_[mask_color, :] = 255

6. cv2 library is used to display the images,cvtColor(), imshow(). resize():

   # Convert an image from RGB color to BGR (blue, green, red)
   cvtColor(_image, cv2.COLOR_RGB2BGR)
   # show the image rescaling it by 0.15 factor
   cv2.imshow("image", cv2.resize(cv2_image, (0, 0), fx=0.15, fy=0.15))
   
   cv2_rgb_masked = cv2.cvtColor(_image_, cv2.COLOR_RGB2BGR)
   cv2.imshow("image masked ", cv2.resize(
     cv2_rgb_masked, (0, 0), fx=0.15, fy=0.15))

Segmentation_test.py

This script implements the test of segmentation to a set of images.

1. each image of the set

   1. resized to a standard size (4800x6800)
   2. divided in blocks of 100x100
   3. flattened

2. the mean image is computed and then it is rescaled by 10 factor (480x680) and converted to a floating point format

3. The script computes Felsenszwalb’s efficient graph based image segmentation using skimage.segmentation.felzenszwalb().

   # scale sets an observation level - Higher scale means less and larger segments
   # sigma is the diameter of a Gaussian kernel, used for smoothing the image prior to segmentation
   
   def felzenszwalb_and_imshow(img):
   
      segments_fz = felzenszwalb(img, scale=100, sigma=0.5, min_size=2000)
      boundaries = mark_boundaries(img, segments_fz, color=(1, 0, 0))
   
      print(f"Felzenszwalb number of segments: {len(np.unique(segments_fz))}")
      print_img((img, boundaries, segments_fz), ("Original image", "Felzenszwalbs's method",
              "Segments from felzenszwalb"))
   
      return len(np.unique(segments_fz))  # return the number of segments

4. Segments image using k-means clustering, implemented by skimage.segmentation.slic():

   # n_segments is the approximiate number of segments
   # comactness balances color proximity and space proximity. Higher values give more weight to space proximity, making superpixel shapes more square/cubic
   def slic_and_imshow(img):
    segments_slic = slic(img, n_segments=250,
                         compactness=10, sigma=1, start_label=1)
    boundaries = mark_boundaries(img, segments_slic, color=(1, 0, 0))
    
    print(f"SLIC number of segments: {len(np.unique(segments_slic))}")
    print_img((img, boundaries, segments_slic), ("Original image", "K-means clustering (slic)", "Segments from l-means clustering"))
   
    return len(np.unique(segments_slic))

5. Segments image using quickshift clustering - the quickshift mode-seeking algorithm, implemented by skimage.segmentation.quickshift():

   # ratio alances color-space proximity and image-space proximity
   # kernel_size width of Gaussian kernel used in smoothing the sample density
   # max_dist cut-off point for data distances. Higher means fewer clusters
   def quick_and_imshow(img):
      segments_quick = quickshift(img, kernel_size=3, max_dist=6, ratio=0.5)
      boundaries = mark_boundaries(img, segments_quick, color=(1, 0, 0))
   
      print(f"Quickshift number of segments: {len(np.unique(segments_quick))}")
   
      print_img((img, boundaries, segments_quick), ("Original image", "Quick method", "Segments from quickshift"))
      
      return len(np.unique(segments_quick))  
      

6. Find edges in an image using the Sobel filter using skimage.filters.sobel():

   def sobel_and_imshow(img):
      segments_sobel = sobel(rgb2gray(img))
   
      print_img((img, segments_sobel), ("Original image", "Sobel filter"))
   
      return len(np.unique(segments_sobel))

7. Find watershed basins in image flooded from given markers using skimage.segmentation.watershed):

   # markers the desired number of markers 
   def watershed_and_imshow(img):
      gradient = sobel(rgb2gray(img))
      segments_watershed = watershed(gradient, markers=250, compactness=0.001)
      boundaries = mark_boundaries(img, segments_watershed, color=(1,0,0))
      
      print(f"Watershed number of segments: {len(np.unique(segments_watershed))}")
      
      print_img((img, boundaries, segments_watershed), ("Original image", "Watershed method", "Segments from watershed"))
      
      return len(np.unique(segments_watershed))
      

Soil_Segmentation.py

This script applies the otsu segmentation algorithm on a set of sample images. The output images are available in SOIL_SEGMENTED_IMAGES and SOIL_SEGMENTED_MASKS directories.

1. it divides each image in 32x32x1 blocks
2. it flattens it and compute the mean, max and median views
3. then it extracts the ycbcr profile colors and applies the otsu segmentation algorithm:

   # Threshold value based on Otsu's method
   val = filters.threshold_otsu(image)
   # resize the flattened image to match the original size
   image = resize(image, (_image_.shape[0], _image_.shape[1]))
   mask_color = image < val
   mask = image < val
   _image_[mask_color, :] = 255

4. finally, it writes the images segmented and the masks to the output directories

Timing analysis

Average latencies report of Soil_Segmentation.py. Script executed on AMD Ryzen 7 5800HS @ 3.2GHz (silent mode).

Code	Latency (ms)	Latency (ms)
Divide and flatten	48.99 ms	48994 μs
Mean filter	25.79 ms	25792 μs
Max filter	40.13 ms	40134 μs
Median filter	392.20 ms	392196 μs
rgb2ycbcr color space	0.59 ms	588 μs
Otsu threshold	0.92 ms	917 μs
Resize image	467.52 ms	467519 μs
Mask	24.64 ms	24644 μs
Masked image	308.34 ms	308343 μs
Writing images	402.91 ms	402909 μs
Total avg time elapsed	1576.34 ms	1576341 μs

Script executed on AMD Ryzen 7 5800HS @ 3.2GHz (boost mode).

Code	Latency (ms)	Latency (ms)
Divide and flatten	38.53 ms	38533 μs
Mean filter	20.77 ms	20771 μs
Max filter	31.32 ms	31324 μs
Median filter	300.19 ms	300193 μs
rgb2ycbcr color space	0.51 ms	513 μs
Otsu threshold	0.74 ms	743 μs
Resize image	373.20 ms	373198 μs
Mask	26.72 ms	26725 μs
Masked image	237.42 ms	237425 μs
Writing images	312.81 ms	312815 μs
Total avg time elapsed	1245.39 ms	1245386 μs

This latencies are the average of the elapsed times of 23 image elaborations (Sample_images). These tables show that there are two parts which can be implemented on FPGA about 10 to 20 times slower than the rest of the code and they are:

1. Median filter
2. Masked image (code snippet to mask the image)

Writing image files which consists on saving elaborated images is an exclusive task of the ARM and cannot implemented on the FPGA.