28

Sep

2016

The Amazing Race!

Author: Zachary Norman

It wasn't so long ago that the IDL-Python bridge was introduced to IDL 8.5. It was with this new version, that I got my first experience programming with Python and testing the IDL-Python bridge. Through the past year it has been exciting to see the new changes and improvements that have become a part of the bridge.

Some of these new features include:

-Better error catching with the IDL-Python bridge
-Enhanced Jupyter notebook that allows for the development of full IDL programs and Python code in the same environment
-Improved support for variables passing back and forth

With all the time I have spent with Python, I have always wondered what some of the advantages are between Python and IDL. One thing that I have commonly heard several engineers say was that IDL was much faster than Python. For this blog, I decided to put that to the test and see how Python and IDL really compared to one another.

Before talking about the test, I do just want to explain things a bit about how it was set up and some potential caveats about the processing times that will be shown. With the tests I created, I did my best to choose tests that were comparable between IDL and Python. Since I'm no expert at Python, there very well may have been other methods that could be faster than what I will show. Most of the pieces I included in the test are things I found easily by doing a web search - meaning that most of the approaches I used were the most common programming methods that people are likely using. This shows how much faster IDL might be than a general program than something that someone might write in Python.


The test:

Here is what was actually tested between IDL and Python with an array of [10000,10000] or 10000*10000 elements

-Array creation time
-Type conversion times
-Index array creation times (i.e. [0,1,2,3,4...,n-1])
-Incrementing array values of all elements by 1
-Complex math expression with array (exact equation: sqrt(sin(arr*arr)))
-Single threaded for IDL and multithreaded
-Array element access times (i.e. setting y = arr[i])
-Simple image processing filter times (filters: sobel, roberts, prewitt)


The results:

Average array creation time (seconds):
    Python : 0.213000 +/- 0.00953933
    IDL    : 0.0936666 +/- 0.0155028
Total time (seconds):
    Python : 0.639000
    IDL    : 0.281000
Python/IDL time ratio: 2.27402

Average array data type conversion time (seconds):
     Python : 0.171333 +/- 0.0155028
    IDL    : 0.0730000 +/- 0.00866031
Total time (seconds):
     Python : 0.514000
     IDL    : 0.219000
Python/IDL time ratio: 2.34703


Average index array creation time (seconds):
     Python : 0.229000 +/- 0.00866031
     IDL    : 0.124667 +/- 0.0160104
Total time (seconds):
    Python : 0.687000
    IDL    : 0.374000
Python/IDL time ratio: 1.83690


Average increasing array value time (seconds):
     Python : 0.0933333 +/- 0.000577446
     IDL    : 0.0313334 +/- 0.000577377
Total time (seconds):
    Python : 0.280000
    IDL    : 0.0940001
Python/IDL time ratio: 2.97872


Average complex math statements (1 thread) time (seconds):
    Python : 6.36967 +/- 0.0645319
    IDL    : 8.34667 +/- 0.0155028
Total time (seconds):
    Python : 19.1090
    IDL    : 25.0400
Python/IDL time ratio: 0.763139


Average complex math statements (8 thread) time (seconds):
    Python : 6.34400 +/- 0.0321871
    IDL    : 1.93933 +/- 0.00923762
Total time (seconds):
    Python : 19.0320
    IDL    : 5.81800
Python/IDL time ratio: 3.27123


Average loop through array element time (seconds):
    Python : 11.5290 +/- NaN
    IDL    : 3.29100 +/- NaN
Total time (seconds):
    Python : 11.5290
    IDL    : 3.29100
Python/IDL time ratio: 3.50319


Average image processing routines time (seconds):
    Python : 15.3660 +/- 0.0829635
    IDL    : 1.39900 +/- 0.0238955
Total time (seconds):
    Python : 46.0980
    IDL    : 4.19700
Python/IDL time ratio: 10.9836



Conclusion:

In short, IDL significantly outperformed Python will all the speed tests apart from the complex math statement. However, IDL has access to built in multithreading for large arrays and, with multithreading enabled, IDL outperforms Python significantly when using all available cores.

Below is the IDL code used to compare the processing speed of IDL and Python. To use it you will need a few Python modules which can be found at the beginning of the procedure "python_test". 

IDL Code:
pro printTimes, pythontimes, idltimes, TIMED = timed
  compile_opt idl2
  
  ;check if we have a name for the process to print with the mean
  if ~keyword_set(timed) then begin
    add = ' time (seconds):'
  endif else begin
    add = timed + ' time (seconds):'
  endelse
  
  print, 'Average ' + add 
  print, '  Python : ' + strtrim(pythontimes.mean(),2) + ' +/- ' + strtrim(stddev(pythontimes),2)
  print, '  IDL    : ' + strtrim(idltimes.mean(),2) + ' +/- ' + strtrim(stddev(idltimes),2)
  print, 'Total time (seconds):'
  print, '  Python : ' + strtrim(total(pythontimes),2)
  print, '  IDL    : ' + strtrim(total(idltimes),2)
  print, 'Python/IDL time ratio: ' + strtrim(total(pythontimes)/total(idltimes),2)
  
end

pro python_test
  compile_opt idl2
  
  
  ;initialize Python
  >>> 'import numpy as np'
  >>> 'import time'
  >>> 'import math'
  >>> 'from idlpy import *'
  >>> 'from skimage.filters import roberts, sobel, prewitt'
  
  ;number of times we want to run each test
  nloops = 3
  
  ;array dimensions for iterator tests
  dims = 10000
  ;>>>'dims = ' + strtrim(dims,2)
  python.dims = dims
  
  ;array dimensions for filter tests
  dims_filter = 10000
  python.dims_filter = dims_filter
  
  ;initialize arrays to hol dpython times
  pythontimes = fltarr(nloops)
  idltimes = fltarr(nloops)
  
  ;test array creation in Python and IDL
  for i=0,nloops-1 do begin
    tStart = systime(/seconds)
    >>> 'arr = np.ones((dims,dims))'
    tEndPython = systime(/seconds)
    arr = fltarr(dims,dims)
    tEndIDL = systime(/seconds)
    
    ;save the times
    pythontimes[i] = tEndPython - tStart
    idltimes[i] = tEndIDL - tEndPython
  endfor
  print
  printTimes, pythontimes, idltimes, TIMED = 'array creation'

  ;check type conversion times
  for i=0,nloops-1 do begin
    tStart = systime(/seconds)
    >>> 'arr2 = arr.astype(int)'
    tEndPython = systime(/seconds)
    arr2 = fix(arr)
    tEndIDL = systime(/seconds)
    
    ;save the times
    pythontimes[i] = tEndPython - tStart
    idltimes[i] = tEndIDL - tEndPython
  endfor
  print
  printTimes, pythontimes, idltimes, TIMED = 'array data type conversion'

  ;check index array creation times
  for i=0,nloops-1 do begin
    tStart = systime(/seconds)
    >>> 'arr = np.arange(0, dims*dims, dtype=long)'
    tEndPython = systime(/seconds)
    arr = lindgen(dims*dims)
    tEndIDL = systime(/seconds)

    ;save the times
    pythontimes[i] = tEndPython - tStart
    idltimes[i] = tEndIDL - tEndPython
  endfor
  print
  printTimes, pythontimes, idltimes, TIMED = 'index array creation'

  ;check adding values to array
  for i=0,nloops-1 do begin
    tStart = systime(/seconds)
    >>> 'arr += 1'
    tEndPython = systime(/seconds)
    arr += 1
    tEndIDL = systime(/seconds)

    ;save the times
    pythontimes[i] = tEndPython - tStart
    idltimes[i] = tEndIDL - tEndPython
  endfor
  print
  printTimes, pythontimes, idltimes, TIMED = 'increasing array value'

  ;check complex math expressions with a single thread
  pref_set, 'IDL_CPU_TPOOL_NTHREADS', 1, /commit
  for i=0,nloops-1 do begin
    tStart = systime(/seconds)
    >>> 'y = np.sin(arr*arr)**.5'
    tEndPython = systime(/seconds)
    y = sqrt(sin(arr*arr))
    tEndIDL = systime(/seconds)

    ;save the times
    pythontimes[i] = tEndPython - tStart
    idltimes[i] = tEndIDL - tEndPython
  endfor
  print
  printTimes, pythontimes, idltimes, TIMED = 'complex math statements (' + strtrim(!CPU.TPOOL_NTHREADS,2) + ' thread)'
  pref_set, 'IDL_CPU_TPOOL_NTHREADS', /default, /commit
  
  ;check complex math expressions with all threads
  for i=0,nloops-1 do begin
    tStart = systime(/seconds)
    >>> 'y = np.sin(arr*arr)**.5'
    tEndPython = systime(/seconds)
    y = sqrt(sin(arr*arr))
    tEndIDL = systime(/seconds)

    ;save the times
    pythontimes[i] = tEndPython - tStart
    idltimes[i] = tEndIDL - tEndPython
  endfor
  print
  printTimes, pythontimes, idltimes, TIMED = 'complex math statements (' + strtrim(!CPU.TPOOL_NTHREADS,2) + ' thread)'

  ;check array element access times
  nhere = 1
  for i=0,nhere-1 do begin
    tStart = systime(/seconds)
    >>>'for x in np.nditer(arr):\n    y=x'
    tEndPython = systime(/seconds)
    foreach x, arr do y = x
    tEndIDL = systime(/seconds)

    ;save the times
    pythontimes[i] = tEndPython - tStart
    idltimes[i] = tEndIDL - tEndPython
  endfor
  print
  printTimes, pythontimes[0:nhere-1], idltimes[0:nhere-1], TIMED = 'loop through array element'

  ;check times for image processing
  im_dat = lonarr(dims_filter, dims_filter)
  ;set a square in the middle of the image to 1
  im_dat[.4*dims_filter:.6*dims_filter,.4*dims_filter:.6*dims_filter] = 1
  ;send the array to Python as well
  python.im_dat = im_dat
  
  for i=0,nloops-1 do begin
    tStart = systime(/seconds)
    >>> 'edge_sobel = sobel(im_dat)'
    >>> 'edge_roberts = roberts(im_dat)'
    >>> 'edge_prewitt = prewitt(im_dat)'
    tEndPython = systime(/seconds)
    edge_sobel = sobel(im_dat)
    edge_roberts = roberts(im_dat)
    edge_prewitt = prewitt(im_dat)
    tEndIDL = systime(/seconds)

    ;save the times
    pythontimes[i] = tEndPython - tStart
    idltimes[i] = tEndIDL - tEndPython
  endfor
  print
  printTimes, pythontimes, idltimes, TIMED = 'image processing routines'

  stop
end

Comments (0) Number of views (172) Article rating: 5.0

20

Sep

2016

Using Attributes to Improve Image Classification Accuracy

Author: Jason Wolfe

In this article I will show an example of how you can improve the accuracy of a supervised classification by considering different attributes.

Supervised classification typically involves a multi-band image and ground-truth data for training. The classification is thus restricted to spectral information only; for example, radiance or reflectance values in various wavelengths. However, adding different types of data to the source image gives a classifier more information to consider when assigning pixels to classes.

Attributes are unique characteristics that can help distinguish between different objects in an image. Examples include elevation, texture, and saturation. In ENVI you can create a georeferenced layer stack  image where each band is a separate attribute. Then you can classify the attribute image using training data.

Creating an attribute image is an example of data fusion in remote sensing. This is the process of combining data from multiple sources to produce a dataset that contains more detailed information than each of the individual sources.

Think about what attributes will best distinguish between the different classes. Spectral indices are easy to create and can help identify different features. Vegetation indices are a good addition to an attribute image if different vegetation types will be classified.

Here are some other attributes to experiment with:

  • Texture measures (occurrence and co-occurrence)
  • Synthetic aperture rader (SAR) data
  • Color transforms (hue, saturation, intensity, lightness)

Including elevation data with your source imagery helps to identify features with noticeable height such as buildings and trees. If point clouds are available for the study area, you can use ENVI LiDAR to create digital surface model (DSM) and digital elevation model (DEM) images. Use Band Math to subtract the DEM from the DSM to create a height image, where pixels represent absolute heights in meters. In the following example, the brightest pixels represent trees. You can also see building outlines.

Example height image at 1-meter spatial resolution

Locating coincident point cloud data for your study area (preferably around the same date) can be a challenge. Luckily I found the perfect set of coincident data for a classification experiment.

Example

The National Ecological Observatory Network (NEON) provides field measurements and airborne remote sensing datasets for terrestrial and aquatic study sites throughout the world. Their Airborne Observation Platform (AOP) carries an imaging spectrometer with 428 narrow spectral bands extending from 380 to 2510 nanometers with a spectral sampling of 5 nanometers. Also onboard is a full-waveform LiDAR sensor and a high-resolution red/green/blue (RGB) camera. For a given observation site, the hyperspectral data, LiDAR data, and high-resolution orthophotography are available for the same sampling period.

I acquired a sample of NEON data near Grand Junction, Colorado from July 2013. My goal was to create an urban land-use classification map using an attribute image and training data. For experimentation and to reduce processing time, I only extracted the RGB and near-infrared bands from the hyperspectral dataset and created a multispectral image with these bands. I used ENVI LiDAR to extract a DEM and DSM from the point clouds. Then I created a height image whose pixels lined up exactly with the multispectral image.

I created an Enhanced Vegetation Index (EVI) image using the Spectral Indices tool. Finally, I combined the RGB/NIR images, the relative height image, and the EVI image into a single attribute image using the Layer Stacking tool.

Next, I used the Region of Interest (ROI) tool to collect samples of pixels from the multispectral image that I knew represented five dominant classes in the scene: Water, Asphalt, Concrete, Grass, and Trees. I used a NEON orthophoto to help verify the different land-cover types.

I ran seven different supervised classifiers with the multispectral image, then again with the attribute image. Here are some examples of Maximum Likelihood classifier results:

Maximum Likelihood classification result from the NEON multispectral image

Notice how the classifier assigned some Building pixels to Asphalt throughout the image. 

Here is an improved result using the attribute image. Those pixels are correctly classified as Building now.

Maximum Likelihood classification result using the attribute image

A confusion matrix and accuracy metrics can help verify the accuracy of the classification.

Confusion matrix calculated from the attribute image classification

The following table shows how the Overall Accuracy value is higher with the attribute image when using different supervised classifiers:

  Multispectral image  Attribute Image 
 Mahalanobis Distance  72.5  83.8
 Minimum Distance  57.69  95.22
 Maximum Likelihood  91.3  98.91
 Parallelepiped  57.18  95.8
 Spectral Angle Mapper  54.97  61.79
 Spectral Information Divergence  62.1  66.93
 Support Vector Machine  85.03  99.16

The accuracy of a supervised classification depends on the quality of your training data as well as a good selection of attributes. In some cases, too many attributes added to a multispectral image can make the classification result worse, so you should experiment with what works best for your study area. Also, some classifiers work better than others when considering different spatial and spectral attributes. Finally, you may need to normalize the different data layers in the attribute image if their pixel values have a wide variation.

New tutorials will be available in ENVI 5.4 for creating attribute images and for defining training data for supervised classification.

Resources

National Ecological Observatory Network. 2016. Data accessed on 28 July 2016. Available on-line at http://www.neonscience.org from National Ecological Observatory Network, Boulder, CO, USA.

Comments (0) Number of views (396) Article rating: 5.0

Categories: ENVI Blog | Imagery Speaks

Tags:

13

Sep

2016

Managing FMV with Jagwire and Passive Analytics

Author: Albert Spratley

The rapid growth of unmanned aerial vehicles (UAVs) and payloads has resulted in an ever growing deluge of data that has to be archived, sorted, analyzed, and distributed to consumers across the defense, agriculture, and utility markets. In many cases, especially in the case of full motion video (FMV), a single flight can result in several hours of data that has to be viewed and analyzed. Often only a small fraction of that data is useful for analysis purposes. For larger UAV fleets, with multiple, simultaneous missions, substantial resources are required to perform the analysis. The resources required to analyze these data products increases cost proportionally to the amount of data collected.

For systems that have adopted the use of properly formatted metadata, we can attempt to filter this glut of data by analyzing patterns and attempting to infer some operator intent based on domain knowledge. For example, identifying temporal “pauses” for the sensor center field of view may indicate an area or point of interest for further analysis. Circular patterns in the sensor center field of view could indicate the inspection of a building, object, or structure of significance. Smooth “pans” during the video or “sweeping” motions across the ground can infer a collection aimed at covering an area on the ground.

Harris Geospatial Solutions has designed and prototyped algorithms capable of identifying these useful segments of video by analyzing the metadata embedded within the video stream. These “passive analytics” run in real time, during the UAV flight, and identify sub-sections of video that are far more likely to be useful in a more detailed analysis. By dynamically detecting, and setting aside these sub-clips of video, the burden of first-phase analysis can be greatly reduced, allowing the user to focus their analytical and dissemination resources on meeting the challenges of their market space rather than wading through a sea of irrelevant data.

Comments (0) Number of views (356) Article rating: No rating

9

Aug

2016

Basic feature extraction with ENVI ROIs

Author: Barrett Sather

The ENVI feature extraction module allows a user to extract certain features from an image. The first step is to segment the image in to regions. Once this is done, you can use spatial and spectral qualities of the regions to pull out targets.

For this example, let’s find the traffic lights in this image taken from Wikipedia’s page on “traffic light”. https://en.wikipedia.org/wiki/Traffic_light

The first step is to do the segmentation with ENVI’s FXSegmention task. I’ve set a merge value of 85 and a segment value of 95, which I determined to be good values for this dataset through manual inspection in ENVI.

Once these segments are created, they can be converted to ENVI ROIs. In this case, I used the ImageThresholdToROI task.

From here, now we can look at qualities of each region, and create a scoring algorithm for how likely each ROI is the feature we are looking for. I felt like looking for the color red was not viable, as then this algorithm breaks down when the light changes. So instead, because our target regions we want to find are circular, let’s look for that.

There are two things we know about circles that we can check for. We know that they are symmetrical, and we know the relationship of the area and circumference as the radius of the circle increases. To start, I’ve set up the pixel coordinates of each ROI so that the origin is the average x/y location for that region. I’ve also calculated the distance away from the origin at each pixel.

x_avg = total(addr[0,*])/n_pix

y_avg = total(addr[1,*])/n_pix

x_coords = reform(addr[0,*] - x_avg)

y_coords = reform(addr[1,*] - y_avg)

power = (x_coords^2 +y_coords^2)

To test for symmetry, take the min and max of x and y, and add them together. The closer to zero that this value is, the more symmetric the ROI is.

abs(max(x_coords)+min(x_coords))+ abs(max(y_coords)+min(y_coords))

To test for how circular the area is, we can test the slope of how the distance from the origin increases. Since we know how this relationship between area and circumference behaves, we can find what slope we are looking for.

Slope =  = 1/3

Because a high score for symmetry is zero, let’s use the same scoring system for this measure.

line = linfit(lindgen(n_elements(power)),power[sort(power)])

score = abs(line[1]-.3333)

The final step is to assign weights (I used weights of 1) and calculate an overall score. The full code for extracting traffic lights (or any other circles) can be found at the bottom of this post.

This method for feature extraction takes a while to develop and perfect, which leads me to some exciting news for those that need to quickly develop feature extraction models. There is another method that we have been developing here in house to do feature extraction called MEGA, which is available through our custom solutions group. This is a machine learning tool that takes in examples of a feature you are looking for, and then generates a heat map of where is it likely that that feature is located in an image.

Stay tuned for an in depth look at how this new method compares to classic feature extraction techniques like the one I’ve presented above.

 

pro traffic_light_fx

 compile_opt idl2

 

 file = 'C:\Blogs\FX_to_ROIs\Stevens_Creek_Blvd_traffic_light.jpg'

 e=envi()

 view = e.GetView()

 ;File was created from FX segmentation only. High Edge and Merge settings.

 t0 = ENVITask('FXSegmentation')

 raster = e.Openraster(file)

 t0.INPUT_RASTER = raster

 layer = view.CreateLayer(raster)

 t0.OUTPUT_RASTER_URI = e.GetTemporaryFilename('.dat')

 t0.MERGE_VALUE = 85.0

 t0.SEGMENT_VALUE = 90.0

 t0.Execute

 

 t1 = envitask('ImageThresholdToROI')

 t1.INPUT_RASTER = t0.OUTPUT_RASTER

 t1.OUTPUT_ROI_URI = e.GetTemporaryFilename('.xml')

 loop_max = max((t0.OUTPUT_RASTER).GetData(), MIN=loop_min)

 num_areas = loop_max-loop_min+1

 t1.ROI_NAME = 'FX_Area_' + strtrim(indgen(num_areas)+1, 2)

 t1.ROI_COLOR = transpose([[bytarr(num_areas)],[bytarr(num_areas)],[bytarr(num_areas)+255]])

 arr = indgen(num_areas) + loop_min

 t1.THRESHOLD = transpose([[arr],[arr],[intarr(num_areas)]])

 t1.execute

 

 allROIs = t1.OUTPUT_ROI

 c_scores = []

 c_loc = []

 for i=0, n_elements(allROIs)-1 do begin

   addr = (allROIs[i]).PixelAddresses(raster)

   n_pix = n_elements(addr)/2

   if n_pix gt 60 then begin

     x_avg = total(addr[0,*])/n_pix

     y_avg = total(addr[1,*])/n_pix

     x_coords = reform(addr[0,*] - x_avg)

     y_coords = reform(addr[1,*] - y_avg)

     power = (x_coords^2 + y_coords^2)

     line = linfit(lindgen(n_elements(power)),power[sort(power)])

     this_score = abs(max(x_coords) + min(x_coords)) + max(y_coords) + min(y_coords)) + abs(line[1]-.3333)

     c_scores = [c_scores, this_score]

     c_loc = [c_loc, i]

   endif

 endfor

 idx = sort(c_scores)

 sort_scores = c_scores[idx]

 sort_locs = c_loc[idx]

 for i=0, 3 do begin

   lightROI = allROIs[sort_locs[i]]

    !null = layer.AddRoi(lightROI)

  endfor

end

Comments (0) Number of views (641) Article rating: 5.0

Categories: ENVI Blog | Imagery Speaks

Tags:

3

Aug

2016

Take Your Analytics to the Data

Author: Joey Griebel

ENVI Services Engine: Change Detection

I was recently on NASA’s Earth Observatory where I was reading about the shrinking glacier in Montana’s Glacier National Park. They acquired 8 images dating from 1984 going through 2015 that focused near Lake McDonald where the Jackson and Blackfoot Glaciers were very visible in false color images. They put together a slideshow showing the time series of the extent of the glacier loss due to global warming and the changing climate. The one thing that stuck out in my mind was how time consuming it must have been to search for the data to find images that didn’t have clouds, than download the data… and after that you would still need to run your analysis or comparison.

In the past few years we have started to see a shift to where users want to run the analysis where this data resides, for instance on applications like Amazon Web Services. As we continue to implement more and more ENVI tasks on Amazon Web Services, you can truly take your analysis to the where the data is. In the Change Detection example below, one of our engineers put together a quick interface utilizing ESRI basemaps to define the area of interest. By linking to ESRI’s endpoint (Landsat.arcgis.com) you can stream in the Landsat data available for that area of interest. In this example, we can search for the area we are interested in seeing (Glacier National Park), see what data is available during different years, filter out the data based upon cloud coverage, and then apply Spectral Indices if wanted:

Once you have found the two scenes you want, you simply click change detection and the ENVI tasks run through the steps of the normal analysis and provide updates along the way:

In no time at all, you are given the results of the quick change detection analysis that shows you in Dark Red (red is what has fled the image from time 1). The blue areas shown in the result are new areas to the scene. In these scenes it looks like snowpack that hasn’t quite melted. If you take a look at the examples provided on Earth Observatory, they focus on the Jackson and Blackfoot glacier, which are the areas you see in the Dark red below:

 

This gives you an idea of how you can further the original visual comparison and create shape files to highlight the glacier loss without having to take the time to comb through data for the right set and then download it for analysis. The possibilities for applications like this are really endless as we continue to wrap ENVI functionality into ENVI tasks. ENVI Services Engine allows you to quickly and easily take the analysis to where the data is and save time on downloading, as well as utilize powerful processing tools.

Give it a try yourself here

 

Comments (0) Number of views (659) Article rating: No rating

Categories: ENVI Blog | Imagery Speaks

Tags:

12345678910 Last

MOST POPULAR POSTS

AUTHORS

Authors

Authors

Authors

Authors

Authors

Authors

Authors

Authors

Authors

GUEST AUTHORS

Authors

Authors

Authors

Authors

Authors

Authors

Authors



© 2016 Exelis Visual Information Solutions, Inc., a subsidiary of Harris Corporation