Commit 28e32f22 authored by luroth's avatar luroth
Browse files

Documentation with demo scripts and readme files

parent 78b3bae5
© Lukas Roth , Group of crop science, ETH Zurich
Part of: [Phenofly UAS data processing tools](../README.md)
# Image segmentation
Sample scripts to segment RAW image pixels in plant and soil.
## Background
Based on:
Roth et al. 2019 (unpublished)
## Example
Run ```ActiveLearningSegmentation/_Demo_.py```
Prerequisit:
1. RAW images
2. Image masks generated (see [Image mask generation (Standalone Agisoft Script)](ImageProjectionAgisoft/README.md))
The campaign folder should therefore contain following subfolders:
- ```RAW```
- ```image_masks```
It will generate the following folders:
- ```segmentation```: Segmented images
......@@ -539,6 +539,7 @@ class Workflow():
# load numpy array
enhanced_raw = np.load(self.path_workspace_descriptors / (image.at['image_name'] + '.npy'))
enhanced_raw = np.delete(enhanced_raw, 19, axis=2)
mask = np.load(self.path_workspace_descriptors / (image.at['image_name'] + '_mask.npy'))
# iterate over plot groups on image
......@@ -636,6 +637,8 @@ class Workflow():
# load numpy array
print("Load...", end="")
enhanced_raw = np.load(self.path_workspace_descriptors / (image.at['image_name'] + '.npy'))
enhanced_raw = np.delete(enhanced_raw, 19, axis=2)
mask = np.load(self.path_workspace_descriptors / (image.at['image_name'] + '_mask.npy'))
print(timeit.default_timer() - tic, "seconds")
......
from ActiveLearningSegmentation.Workflow import Workflow
from Common import Imagefunctions
import os
from pathlib import Path
base_dir = Path(os.path.join(os.path.dirname(os.path.realpath(__file__))))
if base_dir.name == "ActiveLearningSegmentation":
base_dir = base_dir / ".."
# Function to extract global sample id
def plot_group_label_func(df):
return df.plot_label.str[:11]
return df.plot_label
wfl = Workflow(path_RAW = './_TestData/sonyA9/RAW',
wfl = Workflow(path_RAW = base_dir / '_TestData/sonyA9/20190402/28m_M600P/RAW',
raw_extension = ".ARW",
path_masks = './_TestData/sonyA9/image_masks',
path_workspace = './_Demo_als_output',
path_masks = base_dir / '_TestData/sonyA9/20190402/28m_M600P/image_masks',
path_upload = base_dir / '_TestData/sonyA9/20190402/28m_M600P/segmentation',
path_workspace = base_dir / '_temp',
plot_group_label_func = plot_group_label_func,
minimum_group_size=40,
minimum_group_size=1,
sample_type="soil",
border_pixels=12,
read_raw_func=Imagefunctions.read_sonyA9_RAW,
enhanced_raw_func=Imagefunctions.enhance_sonyA9_rawpy,
preview_func=Imagefunctions.preview_sonyA9_rawpy)
#wfl.prepare_workspace()
#wfl.prepare_descriptors()
wfl.prepare_workspace()
#wfl.init_sample_masks()
wfl.prepare_descriptors()
#wfl.save_masks_to_workspace()
wfl.load_masks_from_workspace()
wfl.init_images()
wfl.init_new_classifier()
#path_training_coords_plants = './_TestData/sonyA9/_training_plants_coords.csv'
#path_training_coords_soil = './_TestData/sonyA9/_training_soil_coords.csv'
trainings = ["Training1", "Training2", "Training3", "Training4"]
#wfl.init_training_coords_from_files(path_training_coords_plants, path_training_coords_soil)
#wfl.init_training_data_from_coords()
for training in trainings:
path_training_data = base_dir / '_TestData' / 'sonyA9' / 'segmentation_training' / (training + '.npy')
path_training_response_data = base_dir / '_TestData' / 'sonyA9' / 'segmentation_training' / (training + "_response.npy")
wfl.load_training_data_from_file(path_training_data, path_training_response_data)
wfl.init_sample_masks()
#wfl.save_training_data_from_coords()
wfl.load_training_data_from_file('./_Demo_als_output/Training.npy', './_Demo_als_output/Training_response.npy')
wfl.train_classifier()
wfl.segment_image_cutouts()
wfl.GUI()
#wfl.segment_image_cutouts()
#wfl.GUI()
wfl.segment_images()
#wfl.segment_images()
wfl.upload_predictions()
......@@ -8,7 +8,7 @@ from pathlib import Path
import shutil
import sys
from GroundAerialCoverage import MaskRead
from Common import MaskRead
import time
def affine_transform_images(path_images, path_masks, path_transformed, nx, ny, design_label, offset= 0,
......
......@@ -2,10 +2,6 @@ import rawpy
import numpy as np
import cv2
####
# TEST
path_image = '/home/luroth/PycharmProjects/crop_growth_analysis/Scripts/Analysis/Active_Learning/samples2/RAW/_DSC7431.ARW'
def enhance_for_fft_sonyA9_rawpy(image_raw, pixel_shift):
'''Gets a rawpy.RawPy object and generates 16 bit RGB, HSV, Lab and ExG/ExR representants
......@@ -106,6 +102,8 @@ def postprocess_sonyA9_RAW(path_image, pixel_shift=0):
use_camera_wb=True, no_auto_bright=True)
# Convert XYZ to RGB space
RGB_conventional = cv2.cvtColor(XYZ_conventional, cv2.COLOR_XYZ2RGB)
RGB_conventional = RGB_conventional[pixel_shift:RGB_conventional.shape[0] - pixel_shift,
pixel_shift:RGB_conventional.shape[1] - pixel_shift, :]
return RGB_conventional
......@@ -170,15 +168,26 @@ def enhance_sonyA9_rawpy(image_raw, pixel_shift=0):
Lab_conventional_float32[:, :, 2] = (Lab_conventional_float32[:, :, 2] + 128) / 128
Lab_conventional = np.uint16(Lab_conventional_float32 * 2 ** 16)
# Calcualte ExR and ExG
# Calcualte ExR, ExG and ExB
R, G, B = cv2.split(RGB_conventional)
normalizer = np.add(np.add(R, G), B)
# Cheat: avoid div. by zero
R = np.array(R, dtype=np.float32)
G = np.array(G, dtype=np.float32)
B = np.array(B, dtype=np.float32)
normalizer = np.array(R + G + B)
# Avoid division by zero
normalizer[normalizer == 0] = 1
normalizer = normalizer
r, g, b = (R, G, B) / normalizer
# ExR + ExG
ExR = np.divide(np.subtract(1.4 * R, G), normalizer)
ExG = np.divide(np.subtract(np.subtract(2.0 * G, R), B), normalizer)
ExRG_conventional_float = np.stack([ExR, ExG], axis=2)
ExB = np.array(1.4 * b - g,
dtype=np.float32) # Mao, W., Wang, Y., Wang, Y., 2003. Real-time detection of between-row weeds using machine vision. ASAE paper number 031004. The Society for Agricultural, Food, and Biological Systems, St. Joseph, MI.
ExRG_conventional_float = np.stack([ExR, ExG, ExB], axis=2)
ExRG_conventional = np.uint16(ExRG_conventional_float * 2 ** 16)
# Concat all
......
# Prototype algorithm to transform images to better detect ears
# Author: Simon Treier
# 11. September 2019
import os
import shutil
import timeit
import imageio
import rawpy, time
import pandas as pd
import numpy as np
from Common import FFTfunctions
from Common import Imagefunctions
from skimage.feature import peak_local_max
from skimage.morphology import watershed
from scipy.ndimage.measurements import center_of_mass, label
from scipy.sparse import csr_matrix
from pathlib import Path
from matplotlib import path as geopath
from sklearn.ensemble import RandomForestClassifier
import cv2
from progressbar import ProgressBar
import multiprocessing
from multiprocessing import Process, Manager
from Common import MaskRead
def segment_ear_image(path_raw, path_processed, raw_image_name,
GSD, # Indicate ground sampling distance in mm
clf,
Min_SD, # Indicate minimal sampling distance (dimension of objects) to detect in mm
Max_SD, # Indicate maximal sampling distance (dimension of objects) to detect in mm
pixel_shift, # dim_raw - dim_rgb / 2
local_maxima_threshold,
local_maxima_min_distance,
min_intensity_watershed_area,
RGB_factor,
ExB_factor,
R_factor,
L_factor,
V_factor):
tic = time.process_time()
########################################################################################################################################################################################################################################################################################################################
pathCompl = Path(path_raw) / raw_image_name
#print(pathCompl)
image_name = raw_image_name.split(".")[0]
# Read RAW file
image_raw = rawpy.imread(str(pathCompl))
# Create dir
path_processed = Path(path_processed)
path_processed.mkdir(exist_ok=True)
# Demosaic RAW, transform to color spaces and calculate indices
a_RGB_16bitf, a_ExG, a_ExR, a_ExB, a_ExB_n, \
R, G, B, RGB_intensities, RGB_VARI, \
GCC, RCC, BCC, \
H, S, V, \
L, a, b = Imagefunctions.enhance_for_fft_sonyA9_rawpy(image_raw, pixel_shift)
# Read dimensions from imput image for calculations
rows, cols, dim = a_RGB_16bitf.shape
DimIm = np.max(a_RGB_16bitf.shape)
Dim = np.max(a_RGB_16bitf.shape)
DimPlot = Dim/2
# Build a filter to subsequently filter the FFT images
TEli, TEliInv, Tcrop_eli, X, Y = FFTfunctions.generate_filter_ellipsoid(GSD, Min_SD, Max_SD, DimPlot, a_RGB_16bitf[:,:,1].shape, rows, cols)
# Apply FFT and filter to images
RGB_intensities_trans, ExB_trans, R_trans, L_trans, V_trans = (FFTfunctions.fft_forward_back(band, Tcrop_eli) for band in (RGB_intensities, a_ExB_n, R, L, V))
sum_divisior = RGB_factor + ExB_factor + R_factor + L_factor + V_factor
sum_image = np.multiply(RGB_factor,RGB_intensities_trans) + np.multiply(ExB_factor,ExB_trans) + np.multiply(R_factor,R_trans) + np.multiply(L_factor,L_trans/100) + np.multiply(V_factor,V_trans)
sum_image = np.divide(sum_image,sum_divisior)
sum_image_trans = FFTfunctions.fft_forward_back(sum_image, Tcrop_eli)
local_maxima = peak_local_max(sum_image_trans, min_distance=local_maxima_min_distance, indices=False, threshold_rel=local_maxima_threshold)
local_maxima_labels, n_max = label(local_maxima)
# Extract coords
local_maxima_merged = np.round(center_of_mass( local_maxima, local_maxima_labels, range(1, n_max+1) ))
local_maxima_labels = csr_matrix((np.arange(local_maxima_merged.shape[0])+1, (local_maxima_merged[:, 0], local_maxima_merged[:, 1])), shape=sum_image_trans.shape).toarray()
# Do watershed with min intensity limit
watershed_labels = watershed(-sum_image, local_maxima_labels, mask=sum_image_trans>min_intensity_watershed_area)
result = np.array((watershed_labels >0)*64, dtype=np.uint8)
# Ask random forest if ear
color_value_chart_ = [(int(Y_coord), int(X_coord), R[int(Y_coord),int(X_coord)], G[int(Y_coord),int(X_coord)], B[int(Y_coord),int(X_coord)],
RGB_intensities[int(Y_coord),int(X_coord)],
H[int(Y_coord),int(X_coord)], S[int(Y_coord),int(X_coord)], V[int(Y_coord),int(X_coord)],
L[int(Y_coord),int(X_coord)], a[int(Y_coord),int(X_coord)], b[int(Y_coord),int(X_coord)],
a_ExG[int(Y_coord),int(X_coord)], a_ExR[int(Y_coord),int(X_coord)], a_ExB[int(Y_coord),int(X_coord)],
RGB_VARI[int(Y_coord),int(X_coord)],
GCC[int(Y_coord),int(X_coord)],
RCC[int(Y_coord),int(X_coord)],
BCC[int(Y_coord),int(X_coord)]) for Y_coord, X_coord in local_maxima_merged]
color_value_chart = pd.DataFrame(color_value_chart_,
columns=['Y', 'X', 'R', 'G', 'B', 'RGB_intensities', 'H', 'S', 'V', 'L', 'a', 'b', \
'ExG', 'ExR', 'ExB', 'VARI', 'GCC', 'RCC', 'BCC']
)
RF_decision = clf.predict(color_value_chart)
local_maxima_merged_RF_filtered = local_maxima_merged[RF_decision,:]
local_maxima_merged_RF = np.column_stack([local_maxima_merged, RF_decision])
local_maxima_labels_RF = csr_matrix(
(np.arange(local_maxima_merged_RF_filtered.shape[0]) + 1, (local_maxima_merged_RF_filtered[:, 0], local_maxima_merged_RF_filtered[:, 1])),
shape=sum_image_trans.shape).toarray()
watershed_labels_RF = watershed(-sum_image, local_maxima_labels_RF, mask=sum_image_trans > min_intensity_watershed_area)
result_RF = np.array((watershed_labels_RF > 0) * 128, dtype=np.uint8)
result_ = np.add(result_RF, result)
#print("Save...", end="")
imageio.imsave(path_processed / (image_name + "_watershed.tif"), result_)
pd.DataFrame(np.array(local_maxima_merged_RF, dtype=np.uint16), columns=["y", "x", "ear"]).to_csv(path_processed / (image_name + "_localmax.csv"), index=False)
tac = time.process_time()
#print(tac - tic, "seconds")
def segment_ear_image_worker(work_queue, result_queue, raw_files_folder, output_directory,
GSD, Min_SD, Max_SD, pixel_shift, local_maxima_threshold, local_maxima_min_distance,
min_intensity_watershed_area, RGB_factor, ExB_factor, R_factor, L_factor, V_factor):
# Initiate random forest
clf = RandomForestClassifier(bootstrap = False,
max_depth = 95,
max_features = 6,
min_samples_leaf = 6,
min_samples_split = 4,
n_estimators = 55,
random_state=1, n_jobs=1)
training = pd.read_csv('./EarCount/PixelvaluesOut.csv')
clf.fit(np.array(training[['Y','X','R','G','B','RGB_intensities',
'H','S','V','L','a','b','ExG','ExR','ExB','VARI','GCC','RCC','BCC']]), np.array(training['Organ']=="Ear"))
print(clf.feature_importances_)
for job in iter(work_queue.get, 'STOP'):
raw_filename = job['raw_filename_full']
segment_ear_image(raw_files_folder, output_directory, raw_filename,
GSD, clf, Min_SD, Max_SD, pixel_shift, local_maxima_threshold, local_maxima_min_distance,
min_intensity_watershed_area, RGB_factor, ExB_factor, R_factor, L_factor, V_factor)
result_queue.put(raw_filename)
def segment_ear_images_process(path_raw_files, output_directory,
GSD, # Indicate ground sampling distance in mm
Min_SD=7, # Indicate minimal sampling distance (dimension of objects) to detect in mm
Max_SD=100, # Indicate maximal sampling distance (dimension of objects) to detect in mm
pixel_shift=12, # dim_raw - dim_rgb / 2
local_maxima_threshold=0.30,
local_maxima_min_distance=2,
min_intensity_watershed_area=0.25,
RGB_factor=2,
ExB_factor=1,
R_factor=1,
L_factor=1,
V_factor=1
):
GSD = GSD * 1000
shutil.rmtree(output_directory, ignore_errors=True)
try:
os.mkdir(output_directory)
except WindowsError:
print("Dir not accessible, skip")
return None
# Job and results queue
m = Manager()
jobs = m.Queue()
results = m.Queue()
processes = []
# Progress bar counter
max_jobs = 0
count = 0
# Build up job queue by iterating over RAW files
directory = os.fsencode(path_raw_files)
for file in os.listdir(directory):
# Get filename
raw_filename_full = os.fsdecode(file)
raw_filename = os.path.splitext(raw_filename_full)[0]
if os.path.splitext(raw_filename_full)[1].lower() == ".arw" or \
os.path.splitext(raw_filename_full)[1].lower() == ".cr2":
#print(raw_filename_full, "to queue")
job = dict()
job['raw_filename_full'] = raw_filename_full
job['raw_filename'] = raw_filename
jobs.put(job)
max_jobs += 1
# Init progress bar
progress = ProgressBar(min_value=0, max_value=max_jobs)
# Start processes, number of CPU - 2 to have some left for main thread / OS
for w in range(multiprocessing.cpu_count() - 2):
p = Process(target=segment_ear_image_worker,
args=(jobs, results, path_raw_files, output_directory,
GSD, Min_SD, Max_SD, pixel_shift, local_maxima_threshold, local_maxima_min_distance,
min_intensity_watershed_area, RGB_factor, ExB_factor, R_factor, L_factor, V_factor
))
p.daemon = True
p.start()
processes.append(p)
jobs.put('STOP')
print("jobs all started")
progress.update(0)
# Get pseudo results and increment counter with it
while count < (max_jobs):
results.get()
progress.update(count)
count += 1
progress.finish()
for p in processes:
p.join()
def process_multiview_earcount_generation(path_processed_str, design_label, nx=50 * 20, ny=125, buffer=0, subplots=False,
rmdest=False, use_RF_for_smv=True):
print("Processing design", design_label)
path_processed = Path(path_processed_str)
path_seg_ear = path_processed / 'ear_segmentation'
path_masks = path_processed / 'image_masks'
path_mvEars = path_processed / 'ear_count'
if rmdest:
try:
shutil.rmtree(path_mvEars, ignore_errors=True)
except:
print("Could not remove ", path_mvEars)
path_mvEars.mkdir(exist_ok=True)
# Read all masks
gpd_masks = MaskRead.read_masks(path_masks, design_label)
gpd_masks['geometry'] = gpd_masks.buffer(buffer)
# Generate sample matrix and coords
vx, vy = np.mgrid[0:nx, 0:ny]
samplepoints_array = np.zeros((nx, ny))
samplepoints_coords = np.stack((vx.flatten(), vy.flatten()), axis=1)
samplepoints_normalized_coords = samplepoints_coords / (nx - 1, ny - 1)
# Split plot label in group and subnumber
if subplots:
gpd_masks['plot_group'] = gpd_masks.plot_label.str[0:11]
gpd_masks['plot_subnumber'] = np.int8(gpd_masks.plot_label.str[12:])
# Dicts to store sample ground coverage / aerial coverage arrays in
ear_counts = []
plot_resamples = {}
plot_image_count = {}
plot_labels = np.unique(gpd_masks['plot_label'])
for plot_label in plot_labels:
plot_resamples[plot_label] = samplepoints_array.copy().flatten()
plot_image_count[plot_label] = 0
# Process image by image
image_groups = gpd_masks.groupby("image")
for image, image_group in image_groups:
print(f'Process {image}')
# Load segmented image
path_image = path_seg_ear / (image + '_watershed.tif')
if path_image.exists():
image_segmented = imageio.imread(path_image)
ear_count_coords = np.array(pd.read_csv(path_seg_ear / (image + '_localmax.csv')), dtype=np.uint16)
# Convert to to [0, 1]
image_segmented_filtered = image_segmented.copy()
threshold = 192 if use_RF_for_smv else 64
image_segmented_filtered[image_segmented >= threshold] = 1
image_segmented_filtered[image_segmented < threshold] = 0
# Process plots on image
plot_label_groups = image_group.groupby("plot_label")
for plot_label, plot_label_group in plot_label_groups:
# Process different shapes
for _, plot in plot_label_group.iterrows():
if plot['type'] == 'canopy':
# Get postion of plot edges on image
plot_polygon = plot.at['geometry']
x, y = plot_polygon.exterior.coords.xy
plot_edges_image_coords = np.float32(np.stack([x, y], axis=1))
plot_edges_image_coords_zonalstat = np.float32(np.stack([y, x], axis=1))
# Get affine transform matrix to transform from 0:1, 0:1 space to image coords
plot_edges_normalized_coords = np.float32([[0, 0], [0, 1], [1, 1], [1, 0]])
M_normalized_to_image_coords = cv2.getAffineTransform(plot_edges_normalized_coords[0:3, :],
plot_edges_image_coords[0:3, :])
# Get position of plot raster points on image
samplepoints_image_coords = np.int16(
np.round(np.dot(
np.c_[samplepoints_normalized_coords, np.ones(samplepoints_normalized_coords.shape[0])],
M_normalized_to_image_coords.T)))
# Filter border pixels
border = 100
is_border = np.any(samplepoints_image_coords[:, 1] > image_segmented.shape[0] - 1 - border) or \
np.any(samplepoints_image_coords[:, 1] < 0 + border) or \
np.any(samplepoints_image_coords[:, 0] > image_segmented.shape[1] - 1 - border) or \
np.any(samplepoints_image_coords[:, 0] < 0 + border)
if not is_border:
# Sample oblique segmented image with samplepoints
sampled = image_segmented_filtered[
samplepoints_image_coords[:, 1], samplepoints_image_coords[:, 0]]
##DEBUG
'''
plt.figure()
plt.imshow(image_segmented_filtered)
plt.scatter(samplepoints_image_coords[:, 0], samplepoints_image_coords[:, 1])
plt.scatter(samplepoints_image_coords[50, 0], samplepoints_image_coords[50, 1])
plt.scatter(samplepoints_image_coords[100, 0], samplepoints_image_coords[100, 1])
plt.scatter(samplepoints_image_coords[200, 0], samplepoints_image_coords[200, 1])
plt.scatter(samplepoints_image_coords[600, 0], samplepoints_image_coords[600, 1])
plt.scatter(samplepoints_image_coords[125:250, 0], samplepoints_image_coords[125:250, 1])
plt.scatter(samplepoints_image_coords[1000:1125, 0], samplepoints_image_coords[1000:1125, 1])
plt.show()
plt.figure()
sampled_img = sampled.reshape((50, 125))
plt.imshow(sampled_img.T)
plt.show()
'''
plot_resamples[plot_label] += sampled
plot_image_count[plot_label] += 1
# Sample plot polygon for ear counts
# Transform coordinates to path
sample_path = geopath.Path(plot_edges_image_coords_zonalstat)
# Sample ears
ear_count_coords_orig = ear_count_coords[:,0:2]
ear_count_coords_RF = ear_count_coords[ear_count_coords[:,2]==1, 0:2]
ears_in_poly_orig = sample_path.contains_points(ear_count_coords_orig, radius=abs(1))
ear_counts_orig = np.sum(ears_in_poly_orig*1)
ears_in_poly_RF = sample_path.contains_points(ear_count_coords_RF, radius=abs(1))
ear_counts_RF = np.sum(ears_in_poly_RF*1)
#from matplotlib import pyplot as plt
#plt.imshow(image_segmented_filtered)
#plt.scatter(ear_count_coords[ears_in_poly,1], ear_count_coords[ears_in_poly,0])
plot['ear_count'] = ear_counts_orig
plot['ear_count_RF'] = ear_counts_RF
del (plot['geometry'])
ear_counts.append(plot.to_dict())
# Normalize GC_ACs with image counts
for plot_label in plot_resamples.keys():
plot_resamples[plot_label] /= plot_image_count[plot_label]
# Save to disk
try:
path_mvEars.mkdir()
except:
pass
if subplots:
plot_groups = np.unique(gpd_masks['plot_group'])
for plot_group in plot_groups:
print(f"Summarizing {plot_group} with subplots")
plots = []
for i in range(20):
plot_label = plot_group + "_" + str(i + 1)
plot = plot_resamples[plot_label].reshape(samplepoints_array.shape[0], samplepoints_array.shape[1])
plot = np.uint8(plot * 255).T
imageio.imsave(str(path_mvEars / (design_label + "_" + plot_label + ".tif")), plot)
plots.append(plot)
all_plots = np.concatenate(plots, axis=1)
imageio.imsave(str