This example follows the same pattern than the examples 2 and 3 but for one interesting add-on. As from the version 0.10, the two endmembers extraction algorithms, ATGP and FNINDR, can extract from a region of interest (ROI) and not all the cube, as it is actually the case. To control the ROI you have to define a binary mask. The binary mask is a 2D matrix, a True value at a position in the matrix means that you want to include the corresponding pixel in the endmembers search. The effluents exiting a smokestack are concentrated at the exit. We will define a small ROI just at this exit to catch the different signals that determine the gas.
Look at the Reveal Viewer screen at figure 1. The selected pixel, at the red crossbar, is located just above the smokestack where the gas emission is at the maximum. The two broad emission bands between 1075 and 1200 cm-1 correspond to the symmetric S=O stretching vibration of sulfur dioxide (SO2). The cube is acquired in the VLW range of infrared; between 867 to 1288 wavenumber (cm-1) (7.76 to 11.54 micrometer) and it have 165 bands. The instrument used to acquire the data comes from the Telops Company.
First, we extract the endmembers of the whole cube (figure 2). It will be use to compare to the ROI extracted endmembers. The signal tagged EM2 in green is a SO2 gas signature.
Next, we make a small ROI using SAM and the EM2 signature. After some tests, the threshold retained is 0.15. Figure 3 presents the binary mask.
With the binary mask, we extract a second endmembers set. This one is composed mainly with SO2 gas signatures presents in the effluents. Look at the figure 4, the endmembers EM1 to EM7 are the SO2 the signatures.
With this information at hand, we generate two mean images, one with SAM and the other with FCLS. This is a two steps process. First, we apply a SAM classification to the masked endmembers set excluding the background signal. The upper row at figure 4 and 5 presents the results. Next, we construct the mean image by adding the EM1 to EM7 classification images and dividing the result by 7. The same process is applied using FCLS but with an important difference. To be effective, the linear unmixing model needs a good pure pixels set. A way to give it is to replace EM2 of the full cube endmembers set with each SO2 signatures given by the masked endmembers set. For each signature replacement inside the full endmembers set, we call FCLS and keep the related abundance map. In this case, this is always the second abundance map. We present these results maps in the second row of figure 4 and 5. Like with SAM, we generate a mean image by adding each related abundance maps and dividing them by 7. You can see at figure 6 the two mean images.
Figure 4: in red, the masked endmembers set classified with SAM, in colors the same unmixed with FCLS; from left to right: signals EM1 to EM4.
Note
Transmittance is a non-linear process. The SAM and FCLS mean maps are an approximation.
The mean images are an interesting achievement. The FCLS based mean image have a better resolution than the SAM based image. However, the SAM version takes less CPU cycle.
Which one between all these SO2 signatures is a pure pixel? It is difficult to say. EM2 from the full cube endmembers set along with EM2 and EM5 from the masked endmembers set are almost identical. In addition, theirs SAM classification footprints are similar to the mean SAM footprint. Using the same reasoning with FCLS do not give so much evidence of which signal is a pure pixel.
Code follow:
"""
Plot smokestack effluents.
"""
import os
import os.path as osp
import pysptools.classifiers as cls
import matplotlib.pyplot as plt
import pysptools.util as util
import pysptools.eea as eea
import pysptools.abundance_maps as amp
import numpy as np
class Classify(object):
"""
For this problem NormXCorr works as well as SAM
SID was not tested.
"""
def __init__(self, data, E, path, threshold, suffix):
print 'Classify using SAM'
self.sam = cls.SAM()
self.sam.classify(data, E, threshold=threshold)
self.path = path
self.suffix = suffix
def get_single_map(self, idx):
return self.sam.get_single_map(idx, constrained=False)
def plot_single_map(self, idx):
self.sam.plot_single_map(self.path, idx, constrained=False, suffix=self.suffix)
def get_endmembers(data, header, q, path, mask, suffix, output=False):
print 'Endmembers extraction with NFINDR'
ee = eea.NFINDR()
U = ee.extract(data, q, maxit=5, normalize=True, ATGP_init=True, mask=mask)
if output == True:
ee.plot(path, header, suffix=suffix)
return U
def get_abundance_maps(data, U, umix_source, path, output=False):
print 'Abundance maps with FCLS'
fcls = amp.FCLS()
amap = fcls.map(data, U, normalize=True)
if output == True:
fcls.plot(path, colorMap='jet', suffix=umix_source)
return amap
def get_endmembers_sets(data, header, path):
""" Return a endmembers set for the full cube and a
endmembers set for the region of interest (ROI).
The ROI is created using a small region of the
effluents leaving near the smokestack.
"""
# Take the endmembers set for all the cube
U_full_cube = get_endmembers(data, header, 8, path, None, 'full_cube', output=True)
# A threshold of 0.15 give a good ROI
cls = Classify(data, U_full_cube, path, 0.15, 'full_cube')
# The endmember EM2 is use to define the region of interest
cls.plot_single_map(2)
# The effluents region of interest
effluents = cls.get_single_map(2)
# Create the binary mask with the effluents
mask = (effluents > 0)
# Plot the mask
plot(mask, 'gray', 'binary_mask', path)
# And use this mask to extract endmembers near the smokestack exit
U_masked = get_endmembers(data, header, 8, path, mask, 'masked', output=True)
return U_full_cube, U_masked
def classification_analysis(data, header, path, E_full_cube, E_masked):
# Classify with the masked endmembers set
cls = Classify(data, E_masked, path, 0.15, 'masked')
# and plot the results
cls.plot_single_map('all')
# Calculate the average image
# 0 to 6, the last image, number 7, is a background, we skip it
for i in range(7):
if i == 0:
gas = cls.get_single_map(i+1)
gas = gas + cls.get_single_map(i+1)
gas = gas / 7
# and plot it
plot(gas, 'spectral', 'mean_SAM', path)
def unmixing_analysis(data, header, path, E_full_cube, E_masked):
# Calculate the average image, but use a trick
# The last image, number 7, is a background, we skip it
for i in range(7):
E_full_cube[1,:] = E_masked[i,:]
amaps = get_abundance_maps(data, E_full_cube, 'masqued_{0}'.format(i+1), path, output=False)
if i == 0:
mask = amaps[:,:,1]
else:
mask = mask + amaps[:,:,1]
thresholded = (amaps[:,:,1] > 0.3) * amaps[:,:,1]
plot(thresholded, 'spectral', 'masqued_{0}'.format(i+1), path)
mask = mask / 7
thresholded = (mask > 0.3) * mask
plot(thresholded, 'spectral', 'mean_FCLS', path)
def plot(image, colormap, desc, path):
plt.ioff()
img = plt.imshow(image)
img.set_cmap(colormap)
plt.colorbar()
fout = osp.join(path, 'plot_{0}.png'.format(desc))
plt.savefig(fout)
plt.clf()
if __name__ == '__main__':
plt.ioff()
data_path = '../data1'
project_path = '../'
result_path = osp.join(project_path, 'results')
sample = 'Smokestack1.hdr'
fin = osp.join(data_path, sample)
if osp.exists(result_path) == False:
os.makedirs(result_path)
data_file = osp.join(data_path, sample)
data, header = util.load_ENVI_file(data_file)
# Telops cubes are flipped left-right
# Flipping them again restore the orientation
data = np.fliplr(data)
U_full_cube, U_masked = get_endmembers_sets(data, header, result_path)
classification_analysis(data, header, result_path, U_full_cube, U_masked)
unmixing_analysis(data, header, result_path, U_full_cube, U_masked)