ESCI497Z: UAS for Environmental Research and Monitoring

Image Segmentation of UAS Imagery

 

Objective: The objective of this exercise is to explore the use of a new software package to conduct a segmentation and classification of a portion of our Glacier Springs image.

 

Introduction: When using moderate resolution imagery (such as LANDSAT imagery with a pixel size of 30 m), pixel-based classification works quite well. Pixel-based classification involves classifying each pixel based solely on its spectral properties without considering the characteristics of any of the surrounding pixel. With very high resolution (VHR) imagery (pixel size of ~1 m down to a few centimeters), pixel-based classification does not work. Instead, the emerging consensus is to perform image segmentation as a pre-processing step prior to classification. Image segmentation involve joining adjacent pixels with similar spectral properties into “image segments.” The image segments can be of irregular shape and size. Image segments are then classified on the basis of the mean and variance of the spectral properties of the pixels in each image segment. The spatial properties (shape and size) of each image segment can also be used as part of the classification process.

For today’s lab we will use a variety of software packages to carry out both the image segmentation and classification processes.

 

Data: The image was acquired on Sept 28, 2016 using our 3DR Solo quadcopter carrying a Canon S100 camera that has been modified to cover the near IR portion of the spectrum. The full image is available here: J:\Saldata\ESCI-497\glacier_sp.  This folder includes both a “red edge” version of the image (that will be our primary focus) and a true color version (acquired on Sept 18) that will be useful for interpretation. The GCPs for this area were in the State Plane coordinate system with coordinates in units of feet (ug!!!) so the pixel size for both images in 0.1 feet. Weird units, I know.

I’ve clipped a small subset from this larger Glacier Springs image and I’ve converted it to UTM coordinates. The pixel size for this image is 3 cm. This is the image that we will be using for our analysis. This image (rededge_ndi2_utm.tif) includes 4 bands (aka “layers”). Bands 1, 2, 3 are the red edge, green and blue, respectively. Band 4 is a “normalized difference index” that was calculated using:

NDI = (red edge – blue)/(red edge + blue)

This is an index to photosynthetic rate. More about this later.

Go ahead and move this entire J:\Saldata\ESCI-497\glacier_sp folder to your workspace on C:\temp I my case, I will be working on C:\temp\wallin

The folder also includes a shapefile (mid_poly) that you will use to cut a subset out of the full “red edge” image.

Clipping a Polygon: Open Arc and add both images to your table of contents along with the “mid_poly” shapefile.

Creating the Polygon:  You don’t need to do this part but here is how I created the shapefile. Go to Customize-Toolbars-Draw. From the Drawing toolbar, go to  Draw a rectangle. After defining the rectangle, go to Drawing-Convert Graphics to Features, and give your new shapefile a name. (Again, I’ve done this part for you and created a shapefile called “mid_poly”).

Clip the image: As we did last week, you now want to clip bot the “red edge” and the true color images using the mid_poly shapefile. T do this, open the Arc Toolbox and go to Data Management tools-Raster-Raster Processing-Clip

From the Clip dialog box, your Input raster is your red edge image, the Output extent is the mid_poly shapefile. It is critically important that you check the use Input Features for Clipping Geometry box! This is what forces your input layer to be clipped to the exact boundary of our mid_poly shapefile. Specify an Output Raster Dataset and specify a .tif output file format by just adding a .tif suffix. Click OK to run the clipping procedure. The output will be added to your project.

NOTE; Weird discovery: All tiffs are not created equally….apparently. Very long story but I discovered that bringing the tiff created by Arc directly into Montiverdi results in the image displaying upside down!@$####! After hours of beating my head against the computer, I’ve discovered that, if I bring this tiff into the software package ENVI and THEN export it as a tiff, using a different filename, THIS tiff is displayed properly in Montiverdi. Go figure.

So, I’ve clipped the image for you and converted it to UTM with a pixel size of 3 cm. This subset is called rededge3_utm.tif.

 

 

 

Image Segmentation with Montiverdi and the Orfeo Toolbox:

This is an open source image processing package. I’m still on the steep end of the learning curve with this package but it seems to have some nice features. So here goes….

Open Montiverdi, Then go to File-Open Images. Open your red edge subset image. It should look like this:

 

Image segmentation with the Orfeo Toolbox (OTB):

Go to View-OTB Application Browser. This brings up a window with a bunch of options. Scroll down until you see the options for the “Large-scale Mean-shift Segmentation (LSMSS):

 

Here is a link with some additional information about the steps that we will be running: https://www.orfeo-toolbox.org/CookBook/recipes/improc.html#large-scale-mean-shift-lsms-segmentation

 

And here is a link to the paper that is referenced above Michel_et_al_2015_Stable_Mean_Shift_algorithm

And here is another paper that I found to be quite useful:

Chehata_etal_2014_Object-based change detection using high-resolution multispectral images.pdf

 

The Mean-Shift Segmentation Algorithm:

Briefly, the MSS algorithm works by searching for local modes in the joint feature (spectral) and spatial domains. In geek-speak, the algorithm “searches for local modes in the joint feature and spatial domain of n+2 dimensions, where n is the number of features (spectral bands) that are added to the two spatial dimensions.” So, just like ISODATA, LSMS is identifying pixels that are close together in spectral space but ALSO close together in geographic space; the spatial dimensions. There are two parameters that control the segmentation process. The spatial radius controls distance (number of pixels) that is considered when grouping pixels into image segments. The range radius refers to the degree of spectral variability (distance in the n-dimensions of spectral space) that will be permitted in a given image segment. The range radius should be higher than the maximum spectral difference between intra-segment pixel pairs and lower than the spectral difference between segment pixels and the surrounding pixels outside the segment. The spatial radius should be close to the size of the objects of interest (Chehata et al. 2014).

I’ve done a whole bunch of segmentations with a different image (centimeter resolution image acquired from an unmanned aircraft) using a range of parameters. Just to give you a feel for the different results, here are some examples (subset of our full scene) using different combinations of spatial radius and range radius. In each of the four examples below, I used a minimum region size of 500 pixels (about 0.5 m^2).

Note that as the SR and RR increase, the resulting image segments get larger and seem to include scene components that are obviously quite different (eg, vegetation and rock). After experimenting quite a bit, I found that I was not able to obtain “clean” image segments until getting down to SR: 2 and RR: 4.

Parameter selection in this segmentation step is critical. If the segments are too large, dissimilar objects are included within each segment and withi-segment variability compromises the classification accuracy. If the segments are too small, individual objects in the scene (e.g. a single building or a single tree canopy) will be subdivided into separate image segments and the spatial properties of these objects will be lost. Methods for optimizing parameter selection are under development in the literature. At present, parameter selection is mostly done using trial and error. Here is one paper that presents a possible approach to optimize parameter selection. ((link to Ming paper)).

For today’s lab I’ve just taken a wild guess at parameter selection.

 

Go ahead and perform a segmentation (the 4 steps that follow) using SR: 2 and RR: 4 and a Minimum Regions of 500.

 

Step 1: Large-Scale Mean-Shift Smoothing (LSMSS)

Let’s start with Step 1. Read the documentation, then specify your input image (the red edge subset) and filenames for the “Filtered output” and the “Spatial Image.” (be sure to give each of these output files a .tif suffix) I’d strongly recommend that you put all of your output files in a separate folder. I called mine “segout.”

Take all of the default parameters EXCEPT, you should uncheck the Mode Search option.

.Execute.

 

Run Step 2: LSMSSegmentation

Use the “Filtered output” and the “Spatial Image” from Step 1 as input files for Step 2. Specify an output image (with a .tif suffix). Don’t set a Min Region Size. If you set a Min Region size here, these small regions will be set to the background label and will not be subject to further processing. This would lead to “holes” in your final output. We’ll deal with these small regions in a different way in the next step. Specify a Directory for temporary fields. This can be your segout folder. These temporary files are automatically deleted at the end of this step.

Click Execute but be patient. It appears that nothing is happening but all is well.

When this step finishes, your red_step2.tif will be displayed. Mine looks like this:

The blocky nature of this image is just an artifact of the way that the image segments are labeled at this point and the fact that the segmentation proceeds in “tiles” or subsets of the full image. This means that each of the segments in a given tile are labeled sequentially. The labeling starts with the image segments in the upper left tile (low label numbers; darker shading) and the labels for the segments in the lower right get higher label numbers and hence lighter shading. This issue will be addressed in later steps.

 

Run Step 3: LSMSSmall Regions Merging Application

This step can be used to filter out small segments. Segments below the threshold size are merged with the closest big enough adjacent region based on radiometry.

The “Input image” is your original red edge subset. The “segmented image” is the output from step 2. Specify an output image with a .tif suffix.

Let’s go with a Minimum Region Size of 500. With 3 cm pixels, a region composed of 500 pixels would have an area of about 0.5 m^2.

Click  Execute  and again be patient. It is not clear that processing is underway but it is.

 

Step 4: LSMSS Vectorization

This step produces a vector output of the segmented image. The output is a shapefile that with a polygon that delineates each image segment. For each image segment, the attributed file provides the mean and variance of the brightness value for each band for all pixels in the segment and the number of pixels in each segment.

Again, the “Input image” is your original red edge subset. The “Segmented Image” is the output from step 3. The output file is a vector file and it must have a .shp suffix.

Again, click Execute and be patient.

 

Examine your Image Segmentation Result in Arc

Open ArcMap and use the  button to add both your original image (rededge_ndi2_utm.tif) and the shapefile that you created in Step 4 in OTB (mine is called rendi_step4_min500.shp). Turn both on. If you have the shapefile on top, you can right-click and go to Properties-Display and set the transparency to 50% so you can “see through” the shapefile to see the image below.

 

 

From the Table of Contents, right-click on your shapefile and Open the Attribute Table. The first part of mine looks like this.

 

Each row in the table represents one image segment. So note that in my case, my segmentation consists of 18,090 image segments. The fields in this table are:

FID: is just a segment ID#.

Shape: all are polygons

Label: I have no idea what this is??

nbPixels: the # of pixels in the image segment

meanB0,.., B3: mean value in band 0 (red), 1 (green), 2 (blue), 3 (NDI) for all of the pixels in this image segment

varB0,.., B3: variance of the values in B0,..,B3 for all of the pixels in this image segment

 

All of the information in your attribute table is saved in a .dbf file with the same prefix as your shapefile. Open the .dbf file in Excel and save a copy as an Excel file. For some reason, the FID does not get saved as part of the dbf file. So in the Excel file, you will need to insert a column, label it “ID” and insert segment id# starting with 0. Also, very important, delete the “Label” column!!!!!

So now, we want to identify groups of segments in the image that are similar with respect to nbPixels, mean and variance. This is referred to as image classification. To do this, we will need to move to a different software package.

 

Image classification in R:

Open R Studio. At this point, we want to perform an unsupervised classification using a routine called K-means. This routine will identify image segments that have similar spectral (mean and variance in B0,….B3) and spatial (number of pixels in the segment) properties. Image segments with similar properties will be assigned the same group. We will refer to these groups as Spectral classes and each will have a unique spectral class ID number. At a later step we will then assign these spectral classes to Information classes which will be things like vegetation, rock, logs, water.

Here is the code that you need to run in RStudio. Just copy and paste this stuff one section at a time into the Console window. Note that the text in green that follows the # are just notes describing what each section does:

 

install.packages(rJava”)

 

install.packages("xlsx")

 

library(xlsx)

 

#next line sets your working directory change to where ever you are working

setwd("c://temp/wallin/glacier/segout")

 

#Note that you will need to change the filename and sheetname below if yours are different

#And Note that you will have to change the endRow below if you had a different number of #segments in your file

dat <- read.xlsx(file="rendi_step4_min500.xlsx",sheetName = "rendi_step4_min500",

                 startRow = 2, endRow = 10724,header=FALSE)

dat <- dat[,c(1,2,3,4,5,6,7,8,9,10)]

names(dat) <- c("ID", "nbPixels", "meanB0", "meanB1", "meanB2", "meanB3",

                  "varB0", "varB1", "varB2", "varB3")

 

#next line just prints out the first few rows of the file just to confirm that it was properly read in

head(dat)

 

#next line displays the help file for the KMeans unsupervised classification routine

?kmeans

 

#this runs the kmeans routine on our data using everything except the ID as predictor variables

#it will produce 20 spectral classes, run 100 iterations and use 25 random starting points

#the output goes to an output data matrix called srun.km

srun.km <- kmeans(dat[, -1], 20, iter.max = 100, nstart = 25)

 

#next line displays a bunch of output

print(srun.km)

 

#writes an excel file (use a different name if you prefer) containing a single columns of numbers

#representing the predicted spectral class number

write.xlsx(srun.km$cluster,file="rendiclus.xlsx",sheetName="sheet1")

 

Open the excel file that you just created above (rendiclus.xlsx). Note that it just contains a single column of numbers with a column header of “x.” Insert a column to the left of this and enter a column header; something like “ID”. Then number the rows starting with the number zero. The easiest way to do this is with the Fill-Series function:

 

You will need to scroll to the bottom of the file to figure out what Stop value to enter.

 

At this point, you are done with RStudio

 

Back in Arc…….

Back in Arc, right-click on your shapefile and Open the Attribute Table again. For each image segment in the attribute table, we now want to add the spectral class ID number that we just generated in R using the K-Means unsupervised classification routine.

In the attribute table window, click on the Table option button and go to Joins and Relates-Join. In the Join Data dialog box,

Item 1: is the field in your current attribute table that the Join will be based on; select FID

Item 2: is the Excel file that you just created in R. You will need to click on the folder icon and point to this file and sheet1$ in this file

Item 3: is the field in this Excel file to base the join on. This will be the ID number that you inserted into the first column of this file.

Click OK. Arc will chug for quite awhile as it complete the link for all of the 10,000+ image segments. When it finishes, you will see two more columns added on the right side of your attribute table, “ID” and “x”.Note that you just created a dynamic link between your Attribute Table and the Excel file. The information from the Excel file is not saved into the Attribute table at this point…..and this is OK.

 

Create a new raster file:

At this point, we want to create a new raster file with the same cell size as the original starting image (rededge_ndi2_utm.tif) but with every pixel having the spectral class number that we just added to the Attribute Table above. To do this, go the index and search for Feature to raster. In this dialog box, enter:

Note that the default Output Cell Size is about 0.44 meters. This is because, way back when we created the image segmentation, we specified a minimum region size of 500 cells, which is just under 0.5 meters. So Arc wants to try to create a raster with a cell size that is about the size of the smallest image segment. We don’t want this! We want a cell size of 0.03 m (3 cm) which is the cell size of our original image. Make sure that you enter this properly. Then hit OK. Arc will chug for a bit, then add this new raster to your Table of Contents. Take a look at it. Mine looks like this:

a

Export to ENVI:

Next we need to move to another software package. Before we do, we need to export your new raster to a format that ENVI can read. Right-click on your new raster and go to Data-Export Data. In the Export Raster Data dialog, select a Location and a Name. It is also very important to select a Format which you should set to ENVI. Then Save.

Open ENVI:

From the Start button on your computer, go to All Programs-ENVI-Tools-ENVI5.3 Classic (64-bit). From the ENVI Classic toolbar, go to Open Image File. Open both your spectral class file (mine is called rendi_class1.dat) and also your original image (rededge_NDI2_utm.tif). This will bring up your Available Bands List dialog:

Right-click on your class file and go to Edit header. Note that the file type is currently “ENVI Standard.” Change this to ENVI Classification. When you do, you will be asked how many classes are in your file. You have 20. Entering this value brings up the Class Color Map Editing dialog. Click OK to close this. We’ll come back to this shortly. Click OK in the Header Info dialog to close it.

 

Display your images:

In the Available Bands List, go to No Display-New Display. Click on RGB Color, then, one at a time, click on Band1, Band2, Band3 in your red edge image, then Load RGB. You should see something like this:

Move the Available Bands List to the side, go to Display #1-New Display, then click on Gray Scale, then Band1 in your spectral class image, then Load Band.  You should now see something like this:

 

For each display, you have an Image window (upper left), a Scroll window (lower left) and a Zoom window (lower right). The red box in the Scroll window is the section displayed in the Image window. The red box in the Image window is the section displayed in the Zoom window. You can grab the little red box in either the Scroll or Image window to change the display.

 

Link the displays:

In either Image window, go to Tools-Link-Link Displays. Then OK. Now, as you move around either display, the other display will move to show the same location.

 

Overlay Classification:

Now the fun begins. We need to assign spectral classes to information classes. I’m not sure how many information classes we can break out but at least let’s try for these:

Vegetation

Rocks

Logs

Water

Shadow

 

To do this, in Image window Display #2 (the one with the spectral classes displayed) go to Overlay-Classification. Select our spectral class layer, click OK. This brings up the Interactive Class Tool dialog. Also, from the Image window (either one) go to Tools-Cursor Location Value. You should now have something like this:

 

 

I have moved the scroll window to a forested area in the lower right corner of the image. Note that my cross-hair is on a yellow blog that seems to correspond to some vegetation. The Inquire Cursor window indicates that this is class#4. So, you can go the Interactive Class Tool dialog and go to Options-Edit Class Color Names. In the Class Color Map Editing dialog, click on Class#4. Under Class Name: you can then add to the label by calling this something like “vegetation”. Note that at this point, I’d strongly suggest that you retain the “Class #4” part of the name. You can also change the color assigned to this class. For now, change it to black.

 

Click OK. Back in the Interactive Class Tool window, you will now see that Class #4 has your “vegetation” label. You can now click the “ON” box on and off and you will see a whole bunch of image segments toggle back and forth between the original color and the new color (black) that you have assigned. What you are seeing is every image segment in the image that is spectral class#4 turning on and off. As you can see, there are a large number of image segments in spectral class#4. The same is true for each of the other spectral classes.

 

At this point, you can proceed to assign each of your spectral classes to one of our five information classes:

 

Vegetation

Rocks

Logs

Water

Shadow

 

After doing so, you will need to save your color scheme and class assignments. You can do so in the Interactive Class Tool window by going to File-Save Changes to file.