### Navigation

- index
- modules |
- PySptools 0.12.2 beta documentation »

# Smokestack

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.

Figure 1: Reveal Viewer screen of the smokestack.

First, we extract the endmembers of the whole cube (figure 2). It will be use to compare to the ROI extracted endmembers. The spectrum tagged EM2 in green is a SO2 gas signature.

Figure 2: whole cube endmembers extracted with NFINDR, in green a SO2 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.

Figure 3: binary mask used to find the second endmembers set.

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.

Figure 3: masked endmembers set, EM1 to EM7 are the SO2 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 spectrum. 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: spectra EM1 to EM4.

Figure 5: same as figure 4; from left to right: spectra EM5 to EM7.

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.

Figure 6: at left, the SAM mean image, at right, the unmixed mean image.

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 spectrum 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)

### Navigation

- index
- modules |
- PySptools 0.12.2 beta documentation »