LI Xu's World: ML

Showing posts with label ML. Show all posts

Monday, June 10, 2024

ML: The code of K-Nearest Neighbors

% Created by LI Xu
% Version 1.0
% May 31, 2024

% If you have any question about this code,
% please do not hesitate to contact me via E-mail: 
% jeremy456@163.com

% Blog:
% http://blog.sciencenet.cn/u/lixujeremy
% http://lixuworld.blogspot.com/


clear;
clc;

timebegin=tic;
cur_data=date;
cur_time=fix(clock);
str1=sprintf('%s %.2d:%.2d:%.2d', cur_data, cur_time(4), cur_time(5), cur_time(6)); 
fprintf('Time Begin: ');
fprintf(str1);
fprintf('\n');


% settings********************************************
% BackValue
BackVal=0;
fieldname='id';
cropname='canola';
% ***************************************************

% Source Directory
SouDir='input';
% Destination Directory
DesDir='output';

% All images
files=dir(fullfile(SouDir, "*.tif"));


txtpath=fullfile(DesDir, ['note.txt']);
fid=fopen(txtpath, 'w', 'n', 'US-ASCII');

% Loop
for ii=1:numel(files)
    filename=files(ii).name;
    filepath=fullfile(SouDir, filename);
    [~, fname, ext]=fileparts(filename);
    cr_folder=fullfile(DesDir, fname);
    strname=strsplit(fname, '_');
    strname=strname(end-2:end);
    strname=strjoin(strname, '_');

    if ~isfolder(cr_folder)
        mkdir(cr_folder);
    else
        cmd_rmdir(cr_folder);
        mkdir(cr_folder);
    end


    rows=['shapefiles\BL4rows_'];
    plots=['shapefiles\weeds_grassy_'];
    wfplots=['shapefiles\wf_grassy_'];
    samples=['shapefiles\BL_'];

    rows=[rows, strname, '.shp'];
    plots=[plots, strname, '.shp'];
    wfplots=[wfplots, strname, '.shp'];
    samples=[samples, strname, '.tif'];

    disp(['[', num2str(ii), '\', num2str(numel(files)), ']~', filename]);
    fprintf(fid, '%s\r\n', ['[', num2str(ii), '\', num2str(numel(files)), ']~', filename]);

    % Clip the BL/G community from the image
    simagepath=fullfile(cr_folder, ['BL_', strname, '.tif']);
    GenClip(rows, filepath, 0, simagepath);
    % Convert to the ordinary image
    tifpath=GenOrdTifImage(simagepath);

    % Create the Mask
    maskpaths=GenCompMark(simagepath, samples);









end


fprintf('Time Begin: ');
fprintf(str1);
fprintf('\n');


cur_data=date;
cur_time=fix(clock);
str2=sprintf('%s %.2d:%.2d:%.2d', cur_data, cur_time(4), cur_time(5), cur_time(6)); 
fprintf('Time End: ');
disp(str2);
timespan=toc(timebegin);
fprintf('Time Span: %.4f s\n', timespan);

disp('*******************************************************************');

function tifpath=GenOrdTifImage(inputpath)


    [srcdir, fname, ~]=fileparts(inputpath);
    tifpath=fullfile(srcdir, [fname, '_ordinary.tif']);


    image=imread(inputpath);


    imwrite(image, tifpath);

end



function GenClip(oneshp, filepath, BackVal, otpath)

    [~, fname, ~]=fileparts(oneshp);
    strcmd=['gdalwarp -of GTiff -cutline ', oneshp, ' -cl ', fname, ' -crop_to_cutline '];
    strcmd=[strcmd, '-dstnodata ' num2str(BackVal),' ', filepath, ' ', otpath];
    [~, cmdout]=system(strcmd);

end

function maskpaths=GenCompMark(filepath, samples)


    xlspath='classes.sets.xlsx';
    % Color plate
    uniValues=readcell(xlspath, 'Sheet', 'colorplate');
    uniValues(1, :)=[];
    uniValues(:, 1)=[];


    [soudir, fname, ~]=fileparts(filepath);
    maskpaths=fullfile(soudir, [fname, '_mask.tif']);

    sampimage=imread(samples);


    try
        [image, geo]=readgeoraster(filepath);
        try
            info=geotiffinfo(filepath);
        catch
            info=georasterinfo(filepath);
        end
    catch
        image=imread(filepath);
    end


    sampimage=GenFalse(sampimage);
    %% Use nearest neighbor classifier
    mask=GenNearestNeighborClass(image, sampimage);



    % rendering the mask
    outMat=zeros(size(mask, 1), size(mask, 2), 3);
    outMat=RenderUniValues(mask, uniValues, outMat);

    showmat=[image; outMat];
    % imshow(showmat);

    knn_euclidean=fullfile(soudir, [fname, '_knn_euclidean.png']);
    imwrite(showmat, knn_euclidean);
    

    try
        geotiffwrite(maskpaths, uint8(mask), geo);
    catch
        strcmd=['gdalinfo ', filepath];
        [~, cmdout]=system(strcmd);
        epsg=extractBetween(cmdout, "EPSG:"," got from GeoTIFF keys");
        epsg=epsg{1};
        % geotiffwrite(otpath, uint8(class_imag), geo, 'GeoKeyDirectoryTag', info.GeoTIFFTags.GeoKeyDirectoryTag);
        geotiffwrite(maskpaths, uint8(mask), geo, 'CoordRefSysCode', ['EPSG:', epsg]);
    end


    




end


function output=GenFalse(input)
    output=[];
    [rows, cols, ~]=size(input);
    values=unique(input(:));
    values(values>=100)=[];


    sample_regions=false([rows, cols, numel(values)]);


    % Loop to assign the matrixs
    for ii=1:numel(values)
        val=values(ii);
        index=find(input==val);
        band=sample_regions(:, :, ii);
        band(index)=1;
        sample_regions(:, :, ii)=band;

    end

    output=sample_regions;

end




function mask=GenNearestNeighborClass(image, sampimage)

    [rows, cols, ~]=size(image);

    %  Enhance the image**************************
    ycbcr=rgb2ycbcr(image);
    ycbcr(:, :, 1)=0;
    ycbcr=imadjust(ycbcr, stretchlim(ycbcr), []);

    % imwrite(ycbcr, 'ycbcr.tif');
    % *****************************************


    % classes={'soil', 'canola', 'soybean'};
    % nClasses=numel(classes);
    nClasses=size(sampimage, 3);
    % sample_regions=false([rows, cols, nClasses]);
    sample_regions=sampimage;

    mask=[];
    % select each sample region
    % figure;
    % imshow(image);
    % f=figure;
    % for ii=1:nClasses
    %     set(f, 'name', ['Select sample region for ', classes{ii}]);
    %     sample_regions(:, :, ii)=roipoly(image);
    % end
    % 
    % close(f);


    % Convert RGB to L*a*b colorspace
    lab=rgb2lab(image);
    % Calcualate the mean 'a*' and 'b*' value for each ROI area
    % a=lab(:, :, 2);
    % b=lab(:, :, 3);

    a=ycbcr(:, :, 2);
    b=ycbcr(:, :, 3);
    color_markers=repmat(0, [nClasses, 2]);

    for count=1:nClasses
        color_markers(count, 1)=mean2(a(sample_regions(:, :, count)));
        color_markers(count, 2)=mean2(b(sample_regions(:, :, count)));
    end


    % https://www.youtube.com/watch?v=3hEvcyCJNRc&list=PLEo-jHOqGNyUWoCSD3l3V-FjX9PnHvx5n&index=33
    % Classify each pixel using the nearest neighbor rule
    % Each class marker now has an 'a*' and 'b*' value.
    % You can classify each pixel in the |lab_x| image by calculating the
    % Euclidean distance bewteen that pixel and each marker. The smallest
    % distance will tell you that the pixel most closely matched that
    % marker. For example, if the distance between a pixel and the read
    % color marker is the smallest, then the pixel would be labeled as a
    % red pixel.


   color_labels=0:nClasses-1;
   a=double(a);
   b=double(b);
   distance=repmat(0, [size(a), nClasses]);

   % Perform classification
   for count=1:nClasses
       distance(:, :, count)=((a-color_markers(count, 1)).^2+...
           (b-color_markers(count, 2)).^2).^0.5;

   end

   % The other formulas as follows:
   % https://www.saedsayad.com/k_nearest_neighbors.htm

   [value, label]=min(distance, [], 3);
   label=color_labels(label);

   % clear value distance

   colors=[0, 0, 0; 0, 255, 0; 255, 0, 0];
   y=zeros(size(image));
   l=double(label)+1;

   for m=1:rows
       for n=1:cols
           y(m, n, :)=colors(l(m, n), :);

       end
   end



   mask=uint8(label);

end

Monday, February 5, 2024

ML: A simple SAM program from input image to output mask, saving as .tif format

# https://github.com/facebookresearch/segment-anything/issues/221

import cv2, os
import matplotlib.pyplot as plt
sam_checkpoint='D:/PyTest/kkk/sam_vit_l_0b3195.pth'
model_type="vit_l"
from segment_anything import SamAutomaticMaskGenerator, sam_model_registry
sam = sam_model_registry[model_type](checkpoint=sam_checkpoint)

inpath='D:/PyTest/kkk/dog.jpg'
img_arr=cv2.imread(inpath)
img_arr=cv2.cvtColor(img_arr,cv2.COLOR_BGR2RGB)

mask_generator=SamAutomaticMaskGenerator(sam)
# mask_generator = SamAutomaticMaskGenerator(
#     model=sam,
#     points_per_side=32,
#     pred_iou_thresh=0.86,
#     stability_score_thresh=0.92,
#     crop_n_layers=1,
#     crop_n_points_downscale_factor=2,
# #     # min_mask_region_area=100,  # Requires open-cv to run post-processing
# )
predictor=mask_generator.generate(img_arr)

# Choose the first mask
# mask=predictor[0]['segmentation']
# # Remove background by turn it to white
# img_arr[mask==False]=[255, 255, 255]


newimg = img_arr[:, :, 0] * 0
for ii in range(len(predictor)):
    # print(ii)
    mask=predictor[ii]['segmentation']
    # newimg = img_arr[:, :, 0] * 0
    newimg[mask == True] = ii+1
    # filename = os.path.join('D:/PyTest/kkk/export',str(ii+1)+'.tif')
    # cv2.imwrite(filename, newimg)

# plt.imshow(img_arr)
# plt.axis('off')
# plt.show()
filename='D:/PyTest/kkk/dog_new.tif'
cv2.imwrite(filename, newimg)

ML: SAM changed the background to white or other colors

# https://github.com/facebookresearch/segment-anything/issues/221
import cv2
import matplotlib.pyplot as plt
sam_checkpoint='D:/PyTest/kkk/sam_vit_l_0b3195.pth'
model_type="vit_l"
from segment_anything import SamAutomaticMaskGenerator, sam_model_registry
sam = sam_model_registry[model_type](checkpoint=sam_checkpoint)

inpath='D:/PyTest/kkk/girl.png'
img_arr=cv2.imread(inpath)
img_arr=cv2.cvtColor(img_arr,cv2.COLOR_BGR2RGB)

mask_generator=SamAutomaticMaskGenerator(sam)
predictor=mask_generator.generate(img_arr)

# Choose the first mask
mask=predictor[0]['segmentation']
# Remove background by turn it to white
img_arr[mask==False]=[255, 255, 255]

# plt.imshow(img_arr)
# plt.axis('off')
# plt.show()
filename='D:/PyTest/kkk/girl_new.png'
img_arr=cv2.cvtColor(img_arr,  cv2.COLOR_BGR2RGB)
cv2.imwrite(filename, img_arr)

Thursday, February 1, 2024

ML: Deep Learning for Image Segmentation with TensorFlow

# https://www.analyticsvidhya.com/blog/2023/04/deep-learning-for-image-segmentation-with-tensorflow/#

import cv2
import os
import numpy as np
import tensorflow as tf
from tensorflow import keras
import matplotlib.pyplot as plt
import matplotlib as mpl
from tqdm import tqdm
from sklearn.model_selection import train_test_split


# a list to collect paths of 1000 images
image_path=[]
for root, dirs, files in os.walk('D:/PyTest/kkk/png_images'):
    # iterate over 1000 images
    for file in files:
        path=os.path.join(root, file)
        image_path.append(path)
print(len(image_path))

# a list to collect paths of 1000 masks
mask_path=[]
for root, dirs, files in os.walk('D:/PyTest/kkk/png_masks'):
    # iterate over 1000 masks
    for file in files:
        # obtain the path
        path=os.path.join(root, file)
        # add path to the list
        mask_path.append(path)
print(len(mask_path))

# Create a list to store images
images=[]
for path in tqdm(image_path):
    # read file
    file=tf.io.read_file(path)
    # decode png file into a tensor
    image=tf.image.decode_png(file, channels=3, dtype=tf.uint8)
    # append to the list
    images.append(image)


# Create a list to store masks
masks=[]
for path in tqdm(mask_path):
    file=tf.io.read_file(path)
    # decode png file into a tensor
    mask=tf.image.decode_png(file, channels=1, dtype=tf.uint8)
    masks.append(mask)



def resize_image(image):
    image=tf.cast(image, tf.float32)
    image=image/255.0
    # resize image
    image=tf.image.resize(image, (128, 128))
    return image

def resize_mask(mask):
    mask=tf.image.resize(mask, (128, 128))
    mask=tf.cast(mask, tf.uint8)
    return mask

X=[resize_image(i) for i in images]
y=[resize_mask(m) for m in masks]
print(len(X))
print(len(y))


# split data into 80/20 ratio
train_X, val_X, train_y, val_y=train_test_split(X, y,
                                                test_size=0.2, random_state=0)
# Develop tf Dataset objects
train_X=tf.data.Dataset.from_tensor_slices(train_X)
val_X=tf.data.Dataset.from_tensor_slices(val_X)

train_y=tf.data.Dataset.from_tensor_slices(train_y)
val_y=tf.data.Dataset.from_tensor_slices(val_y)

# verify the shapes and data types
train_X.element_spec, train_y.element_spec, val_X.element_spec, val_y.element_spec

# adjust brightness of image
# don't alter in mask
def brightness(img, mask):
    img=tf.image.adjust_brightness(img, 0.1)
    return img, mask

def gamma(img, mask):
    img=tf.image.adjust_gamma(img, 0.1)
    return img, mask

def hue(img, mask):
    img=tf.image.adjust_hue(img, -0.1)
    return img, mask

def crop(img, mask):
    img=tf.image.central_crop(img, 0.7)
    img=tf.image.resize(img, (128,128))
    mask=tf.image.central_crop(mask, 0.7)
    mask=tf.image.resize(mask, (128,128))
    # cast to integers as they are class numbers
    mask=tf.cast(mask, tf.uint8)
    return img, mask

def flip_hori(img, mask):
    img=tf.image.flip_left_right(img)
    mask=tf.image.flip_up_down(mask)
    return img, mask

def flip_vert(img, mask):
    img=tf.image.flip_up_down(img)
    mask=tf.image.flip_up_down(mask)
    return img, mask

# rotate both image and mask identically
def rotate(img, mask):
    img=tf.image.rot90(img)
    mask=tf.image.rot90(mask)
    return img,mask

# zip images and masks
train=tf.data.Dataset.zip((train_X, train_y))
val=tf.data.Dataset.zip((val_X, val_y))


# perform augmentation on train data only
a=train.map(brightness)
b=train.map(gamma)
c=train.map(hue)
d=train.map(crop)
e=train.map(flip_hori)
f=train.map(flip_vert)
g=train.map(rotate)

# concatenate every new augmented sets
train=train.concatenate(a)
train=train.concatenate(b)
train=train.concatenate(c)
train=train.concatenate(d)
train=train.concatenate(e)
train=train.concatenate(f)

# Setting the batch size
BATCH=64

AT=tf.data.AUTOTUNE

# Buffer size
BUFFER=1000

STEPS_PER_EPOCH=800//BATCH
VALIDATION_STEPS=200//BATCH

train=train.cache().shuffle(BUFFER).batch(BATCH).repeat()
train=train.prefetch(buffer_size=AT)
val=val.batch(BATCH)


# Use pre-trained DenseNet21 without head
base=keras.applications.DenseNet121(input_shape=[128, 128, 3],
                                    include_top=False, weights='imagenet')

skip_names=[
    'conv1_relu',
    'pool2_relu',
    'pool3_relu',
    'pool4_relu',
    'relu'
]


skip_outputs=[base.get_layer(name).output for name in skip_names]
# Building the downstack with the above layers.
# We use the pre-trained model as much, without any fine-tuning
downstack=keras.Model(inputs=base.input, outputs=skip_outputs)
# freeze the downstack layers
downstack.trainable=False


from tensorflow_examples.models.pix2pix import pix2pix
upstack=[
    pix2pix.upsample(512, 3),
    pix2pix.upsample(256, 3),
    pix2pix.upsample(128, 3),
    pix2pix.upsample(64, 3)
]

# define the input layer
inputs=keras.layers.Input(shape=[128, 128, 3])
# downsample
down=downstack(inputs)
out=down[-1]

# prepare skip connection
skips=reversed(down[:-1])

# upsample with skip-connections
for up,skip in zip(upstack, skips):
    out=up(out)
    out=keras.layers.Concatenate()([out,skip])


# define the final transpose conv layer
out=keras.layers.Conv2DTranspose(
    59, 3, strides=2, padding='same',
)(out)

unet=keras.Model(inputs=inputs, outputs=out)

def Compile_Model():
    unet.compile(loss=keras.losses.SparseCategoricalCrossentropy(from_logits=True),
                 optimizer=keras.optimizers.RMSprop(learning_rate=0.001),
                 metrics=['accuracy'])

Compile_Model()

# training and fine-tuning
hist_1=unet.fit(
    train,
    validation_data=val,
    steps_per_epoch=STEPS_PER_EPOCH,
    validation_steps=VALIDATION_STEPS,
    epochs=20,
    verbose=2
)

# select a validation data batch
img, mask=next(iter(val))
# make prediction
pred=unet.predict(img)
plt.figure(figsize=(20,28))

NORM = mpl.colors.Normalize(vmin=0, vmax=58)


k=0
for i in pred:
    plt.subplot(4, 3, 1+k*3)
    i=tf.argmax(i, axis=-1)
    plt.imshow(i, cmap='jet', norm=NORM)
    plt.axis('off')
    plt.title('Prediction')

    # plot the ground truth mask
    plt.subplot(4, 3, 2+k*3)
    plt.imshow(mask[k], cmap='jet', norm=NORM)
    plt.axis('off')
    plt.title('Ground Truth')

    # plot the actual image
    plt.subplot(4, 3, 3+k*3)
    plt.imshow(img[k])
    plt.axis('off')
    plt.title('Actual Image')
    k=k+1
    if k==4:
        break
plt.suptitle('Prediction After 20 Epochs (No Fine-tuning)', color='red',
             size=20)
# plt.show()


downstack.trainable=True
# compile again
Compile_Model()
# train from epoch 20 to 40
hist_2=unet.fit(
    train,
    validation_data=val,
    steps_per_epoch=STEPS_PER_EPOCH,
    epochs=40,
    initial_epoch=20,
    verbose=2
)

# select a validation data batch
img, mask=next(iter(val))
# make prediction
pred=unet.predict(img)
plt.figure(figsize=(20, 30))

k=0
for i in pred:
    plt.subplot(4, 3, 1+k*3)
    i=tf.argmax(i, axis=-1)
    plt.imshow(i, cmap='jet', norm=NORM)
    plt.axis('off')
    plt.title('Prediction')

    plt.subplot(4, 3, 2+k*3)
    plt.imshow(mask[k], cmap='jet', norm=NORM)
    plt.axis('off')
    plt.title('Ground Truth')

    plt.subplot(4, 3, 3+k*3)
    plt.imshow(img[k])
    plt.axis('off')
    plt.title('Actual Image')
    k=k+1

    if k==4: break

plt.suptitle('Predition After 40 Epochs (By Fine-tuning from 21th Epoch)', color='red', size=20)
plt.show()

Wednesday, January 31, 2024

ML: SVM vs. Random Forest for image segmentation

# https://github.com/bnsreenu/python_for_microscopists/blob/master/068b-ML_06_04_TRAIN_ML_segmentation_All_filters_RForest_SVM.py


import numpy as np
import cv2
import pandas as pd

img=cv2.imread("D://Mask//IMG_0686.JPG")
img=cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)

img2=img.reshape(-1)
df=pd.DataFrame()
df['Original Image']=img2

# Generate Gabor features
# To count numbers up in order to give Gabor features a label in the data frame
num=1
kernels=[]

for theta in range(2):  # Define number of thetas
    theta=theta/4.*np.pi
    for sigma in (1, 3):
        for lamda in np.arange(0, np.pi, np.pi/4):
            for gamma in (0.05, 0.5):
                gabor_label='Gabor'+str(num)
                ksize=9
                kernel=cv2.getGaborKernel((ksize, ksize),sigma,theta,lamda,gamma, 0, ktype=cv2.CV_32F)
                kernels.append(kernel)
                fimg=cv2.filter2D(img2, cv2.CV_8UC3, kernel)
                filtered_img=fimg.reshape(-1)
                df[gabor_label]=filtered_img
                print(gabor_label, ': theta=', theta, ': sigma, ', sigma, ': lamda=', lamda, ':gamma=', gamma)
                num=num+1


# Canny Edge
edges=cv2.Canny(img, 100, 200)  # Image, min and max values
edges1=edges.reshape(-1)
df['Canny Edge']=edges1

from skimage.filters import roberts, sobel, scharr, prewitt


# Robert Edge
edge_roberts=roberts(img)
edge_roberts1=edge_roberts.reshape(-1)
df['Roberts']=edge_roberts1


# Sobel
edge_sobel=sobel(img)
edge_sobel1=edge_sobel.reshape(-1)
df['Sobel']=edge_sobel1

# Scharr
edge_scharr=scharr(img)
edge_scharr1=edge_scharr.reshape(-1)
df['Scharr']=edge_scharr1

# Prewitt
edge_prewitt=prewitt(img)
edge_prewitt1=edge_prewitt.reshape(-1)
df['Prewitt']=edge_prewitt1

# Gaussian with sigma=3
from scipy import ndimage as nd
gaussian_img=nd.gaussian_filter(img, sigma=3)
gaussian_img1=gaussian_img.reshape(-1)
df['Gaussian s3']=gaussian_img1

# Gaussian with sigma=7
gaussian_img2=nd.gaussian_filter(img, sigma=7)
gaussian_img3=gaussian_img2.reshape(-1)
df['Gaussian s7']=gaussian_img3

# Median with sigma=3
median_img=nd.median_filter(img, size=3)
median_img1=median_img.reshape(-1)
df['Median s3']=median_img1


# Variance with size=3
variance_img=nd.generic_filter(img, np.var, size=3)
variance_img1=variance_img.reshape(-1)
df['Variance s3']=variance_img1

# Now, add a column in the data frame for the labels
# For this, we need to import the labeled image
labeled_img=cv2.imread("D://Mask//mask.tif")
# Remember that you can load an image with partial labels
# But, drop the rows with unlabeled data

labeled_img=cv2.cvtColor(labeled_img, cv2.COLOR_BGR2GRAY)
labeled_img1=labeled_img.reshape(-1)
df['Labels']=labeled_img1

print(df.head())



# Define the dependent variable that needs to be predicted (labels)
Y=df["Labels"].values
# Define the independent variables
X=df.drop(labels=["Labels"], axis=1)

# Split data into train and test to verify accuracy after fitting the model
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, Y, test_size=0.4, random_state=20)

# Import the model we are using
# RandomForestRegressor is for regression type of problems
# For classification we use RandomForestClassifier
# Both yield similar results except for regressor the result is float
# and for classifier it is an integer

from sklearn.ensemble import RandomForestClassifier
model=RandomForestClassifier(n_estimators=100, random_state=42)


# Train the model on training data
model.fit(X_train, y_train)


# Testing the model by predicting on test data
# and Calculate the accuracy score
# First test predication on the training data itself. Should be good
prediction_test_train=model.predict(X_train)

# Test prediction on testing data
prediction_test=model.predict(X_test)


# Let us check the accuracy on test data
from sklearn import metrics

# First check the accuracy on training data. This will be higher than test data predication accuracy
print("Accuracy on training data= ", metrics.accuracy_score(y_train, prediction_test_train))
print("Accuracy= ", metrics.accuracy_score(y_test, prediction_test))

Monday, January 29, 2024

ML: Auto segmentation using multi-otsu

# Auto segmentation using multi-otsu, https://www.youtube.com/watch?v=YdhhiXDQDl4&list=PLZsOBAyNTZwYx-7GylDo3LSYpSompzsqW&index=36
# https://github.com/bnsreenu/python_for_microscopists/blob/master/115_auto_segmentation_using_multiotsu.py
import matplotlib.pyplot as plt
import numpy as np


from skimage import data, io, img_as_ubyte
from skimage.filters import threshold_multiotsu
from skimage.color import rgb2gray

# Read an image
image=io.imread("D:/Mask/IMG_0252.JPG")
image=rgb2gray(image)


# image=Image.open("D:/Mask/IMG_0252.JPG").convert('L')
# Apply multi-Otsu threshold
thresholds=threshold_multiotsu(image, classes=2)

# Digitize (segment) original image into multiple classes
# np.digitize assign values 0, 1, 2, 3, ... to pixels in each class
regions=np.digitize(image, bins=thresholds)
output=img_as_ubyte(regions)

plt.imsave("D:/Mask/otsu_segment.jpg", output)