13

Oct

2016

Creating Animations and Videos in IDL

Author: Jim Pendleton

Oftentimes it is easier for us to perceive change in images if we animate the images in a video.

These comparisons could be between images of the same subject taken at different times, or images acquired contemporaneously but with different techniques or wavelengths.

We in the Harris Geospatial Solutions Custom Solutions Group frequently include video production and display tools in the applications and solutions we provide to our clients.

Let's say you've been monitoring for optical changes in the galactic cluster Abell 115, an image of which is located in your IDL examples/demo/demodata directory.

IDL> read_jpeg, filepath('abell115.jpg', subdir = ['examples','demo','demodata']), image1

IDL> i = image(image1)

Let's copy a region of interest from one place to another to simulate the discovery of a supernova.

IDL> image2[*, 250:260, 250:260] = image2[*, 265:275, 350:360]

IDL> i.setdata, image2

The change is fairly small, but you should see it centered in the image.

To make the "discovery" more apparent, we can animate the changes.

One display technique involves a fade operation in which images dissolve into one another over time. Each pixel transitions over time from one color to another. Our human visual systems are drawn to areas of change in animations. Depending on the effect that's desired for the viewer, all pixels may transition equally and in parallel, or some may transition more slowly or rapidly than others.

First, select a number of individual "frames" over which the change will occur. The number of frames defines how smoothly the transition takes place between images. A value of 2, for example, would produce a blinking effect with no interpolated transition.

IDL> nframes = 32

Next, loop over the frames, weighting the displayed image by proportional contributions between the two input images at each time step.

IDL> for f = 1., nframes do i.setdata, byte(image2*(f/nframes) + image1*(1-f/nframes))

(For the purpose of code compactness in this example, I'm looping using a floating point value but in general this is not recommended.)

Next you may wish to share your "discovery" with your colleagues in the form of a video. The IDL class IDLffVideoWrite lets you create video output files in a number of different formats.

There are three basic steps you will use to create a video. Create an object that is associated with an output file, create a video stream with the appropriate frame dimensions and frame rate for display, then add each image frame to the stream.

Get the dimensions of the image and define a frame rate, in units of frames per second.

IDL> d = size(image1, /dimensions)

IDL> fps = 30

The division of the number of frames by the frame rate will determine the elapsed time of your animation.

Create a video object. In this example, I will write a file to my temporary directory.

IDL> a = filepath(/tmp, 'abell115.avi')

IDL> v = idlffvideowrite(a)

The IDLffVideoWrite object determines the default format for the video output according to the file name extension, in this case ".avi". See the online documentation for alternative formats.

Create a video stream within the video object, specifying the X and Y dimensions of each image frame, along with the frame rate, using the IDLffVideoWrite::AddVideoStream method.

IDL> vs = v.addvideostream(d[1], d[2], fps)

Rather than animating to the graphics windows, we'll loop over the frames writing data to the video stream instead using the IDLffVideoWrite::Put method.

IDL> for f = 1., nframes do !null = v.put(vs, byte(image2*(f/nframes) + image1*(1-f/nframes)))

There is no explicit call to close the action of writing frames to the file. Instead, we signal the completion of writing the data by simply destroying the object reference.

IDL> obj_destroy, v

The video file is now closed and available for sharing with our colleagues. To test the contents, simply SPAWN the path to the video file, assuming the file name extension is known to the shell.

IDL> spawn, a, /hide, /nowait

In order to make the video more interesting, you might consider adding an annotation such as TEXT or an ARROW to your IMAGE display, then copying the display data as your frame input rather than displaying only the raw data directly.  Additionally, see the documentation for the method IMAGE::CopyWindow, used below.

You can copy and paste the following to your IDL command line to see the result.

IDL> v = idlffvideowrite(a)

IDL> vs = v.addvideostream(d[1], d[2], fps)

IDL> i = image(image1, dimensions = d[1:2])

IDL> ar = arrow([200, 250], [200, 250], color='blue', /data, /current)

IDL> t = text(150, 180, 'Eureka!', color = 'green', /data, /current)

IDL> for f = 1., nframes do begin & $

IDL>   i.setdata, byte(image2*(f/nframes) + image1*(1-f/nframes)) & $

IDL>   !null = v.put(vs, i.copywindow()) & $

IDL> endfor

IDL> obj_destroy, v

IDL> spawn, a, /nowait, /hide


The video file resulting from these steps can be found here.

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

Categories: IDL Blog | IDL Data Point

Tags:

11

Oct

2016

Could satellites be the secret to detecting water leaks?

Author: Cherie Muleh

As a Channel Manager for Harris Geospatial, I look after our distributors in Australia, Latin America, and Southeast Asia. Distributors sell our products in their respective regions, handle support, training, and also uncover new uses for our products. Esri Australia's Principal Consultant Remote Sensing and Imagery, Dipak Paudyal, is driving new insights with data, and has started investigating a new way that satellites could help detect water leaks.

Water utility companies routinely face millions of dollars in lost revenue with wasted water and leaks in their pipeline infrastructure. In a recent blog, Dipak explores if there is a way for satellite data and location analytics to help preserve water loss and also enable utility companies to better identify cracks in their system. As Dipak notes, “…water utilities that are willing to think outside of the box and investigate new technologies such as SAR imagery will be guaranteed to stay ahead of the game.”

Analyzing all of this data requires the use of specialty tools like ENVI SARscape to help users transform raw SAR data into an easy-to-interpret images for further analysis. Check out Dipak's blog and let us know what do you think? Can we help the water industry better map their resources with this type of technology?

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

6

Oct

2016

Minimizing Rounding Errors with IDL's TOTAL function

Author: Atle Borsholm

The idea for my blog topic this week came from a question from Ron Kneusel about the loss of precision using IDL’s TOTAL function. IDL functions are often implemented in a way that makes them run fast, since processing speed is important. This is certainly true for the TOTAL function. If you are not using any of the keywords and additional arguments, then it is basically adding to the sum in the input array order. The exception is if you have enough elements in the array to trigger the multi-threaded code path, which I will ignore for now.
The issue with loss of precision is therefore best seen if you start with the largest term and keep adding smaller and smaller terms. At some point, the small term gets so small compared to the large sum that the additional terms get lost in round off.
Here is an example of a sum that we know should approach 128 (with a deviation in the 52 decimal place):
IDL> data = (1-1d/128) ^ lindgen(15000)
IDL> 128 - total(data)
  5.2580162446247414e-013
So, we are getting an error in the 13th decimal place in this case. This array starts with the largest term and ends with the smallest term, so that is the worst case for the TOTAL implementation. If we reverse the order of the terms, we get:
IDL> 128 - total(reverse(data))
  5.6843418860808015e-014
The error gets 10 times smaller in this case. I decided to compare the algorithm that Ron suggested, which is called Kahan sum, with a divide-and-conquer scheme that I have previously used for execution on a massively parallel GPU (which runs fast with 10000’s of independent threads). The Kahan sum algorithm can be found on Wikipedia and the IDL code is pretty short, (but very slow):
; Kahan algorithm, from wikipedia
function KahanSum, data
  sum = 0d
  c = 0d
  for i=0, data.length-1 do begin
    y = data[i] - c
    t = sum + y
    c = (t - sum) - y
    sum = t
  endfor
  return, sum
end
The massively parallel algorithm will run much faster than the Kahan Sum, and the IDL code is listed here:
; Divide and concure total,
; algorithm lends itself to massively parallel execution
function total_mp, data
  compile_opt idl2,logical_predicate
 
  n = ishft(1ull,total(ishft(1ull,lindgen(63)) lt n_elements(data),/integer))
  pad = make_array(n, type=size(data,/type))
  pad[0] = data[*]
  while n gt 1 do begin
    pad[0:n/2-1] += pad[n/2:n-1]
    n /= 2
  endwhile
  return, pad[0]
end
Here are the result with the decreasing magnitude terms:
IDL> 128 - KahanSum(data)
  5.6843418860808015e-014
IDL> 128 - total_mp(data)
  5.6843418860808015e-014
So, in this case the error is similar to the best case sorted input to TOTAL. As an independent test, I also tried randomly ordering the terms, and both KahanSum and TOTAL_MP, are still consistent on the order of 5e-14:
IDL> 128 - KahanSum(data[order])
  5.6843418860808015e-014
IDL> 128 - total_mp(data[order])
  5.6843418860808015e-014
IDL> 128 - total(data[order])
  4.1211478674085811e-013
My conclusion is that the TOTAL_MP example is just as accurate as the Kahan sum, and has the added advantage of allowing for massively parallel execution, as opposed to the highly sequential execution needed for the Kahan algorithm.
 

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

Categories: IDL Blog | IDL Data Point

Tags:

5

Oct

2016

Choosing the right type of imagery for your project

Author: Jon Coursey

Satellite Imagery is used for a wide range of projects, but often times it is tough to figure out just which type of data is needed. In the world of remote sensing, there are many different satellite sensors, all with different resolutions and multispectral capabilities. In this blog, we will explore some basics of remote sensing, what all is available, and common applications of the data.

Remote Sensing 101

As great as a natural color image can look, it is not collected in the same way as, say, the camera on your iPhone. The process of collecting imagery is most commonly referred to as remote sensing, which basically means objects on earth are detected and classified without direct contact. This blog will focus on passive sensing, as this is what is used for natural color and multispectral images (except radar). Passive sensors gather information via radiation that is either emitted or reflected by whatever object the satellite is trying to capture. This data will be captured with whatever spectral bands a sensor can acquire (generally red, green, blue, and near infrared), then sent down to a ground station for processing. From here, the natural color or multispectral image you might have ordered would be put together.

A wide variety of sensors and band combinations are available, and the best options will depend on what you are trying to view, and the size of your project area. Below are a few situations one might face when deciding how to go about selecting the most suitable captures.

What kind of imagery do I need for a base map?

This is generally our most common request. Engineering, architecture, etc. firms will go for images to view ground conditions before building on their project areas. Essentially, they just need an accurate and high resolution picture. In these cases, natural color is the best option.

How much resolution do I need?

This will come down to how large your project area might be, and what is important that you can view with quality. Say your project requires you to obtain an image of California. For this, you wouldn’t want as high of a resolution as if you needed an image of, say, an airport. This would be an extreme amount of detail for such a large area, and would result in enormous file sizes. Additionally, higher resolution sensors typically capture smaller area sizes, so many would be required to fill such an area. And if full coverage even existed, it would be from images accumulated over the course of many months or years. Lower resolution options will give a better overview of an area of this size, and generally use less images to complete the coverage.

One great option is to obtain high resolution imagery of a project site, and lower resolution for the general area. In the images below you will see captures from both the RapidEye & GeoEye sensors. RapidEye captures at 5 meter (being the lower resolution), and GeoEye captures at 0.5 meter resolution. Both show coverage of islands of Hawaii. As you can see, RapidEye still has great detail for a larger city area, but GeoEye has the ability to see specific features within such a city. 

         Figure 1: GeoEye image (0.5 meter resolution)                   Figure 2: RapidEye image (5 meter resolution)

 

When would I need the multispectral bands?

The most common uses for these are measuring crop health and doing vegetative analysis. Most sensors available at MapMart will be equipped with four bands, which are red, green, blue and near-infrared (or RGB & NIR). An image made from these four bands within the ENVI software environment is below. For crop health, one would measure relative reflectance along the wavelengths of these bands. When crops are not healthy, one will notice a decreased reflectance, mainly in the near-infrared reflectance plateau, but also with the blue and green regions of the spectrum.

Multispectral bands can also be used to detect different characteristics on the earth’s surface. Soil, plants, asphalt, etc. all have different spectral signatures, which would be apparent through these multispectral bands. Depending on how many bands you have available, and where they are along the light spectrum, you can start to see more differences within each surface. For example, an asphalt rooftop may appear as only one color with the standard RGB & NIR bands available. If you were working with shortwave infrared (or SWIR) as well, you would start to see different colors within your image. This is because the differences in the components of asphalt are now being detected.

While many applications of this type of data exist, we have covered the most common ones we are contacted for at MapMart. Fortunately, with such a large archive of satellite imagery already existing, with the opportunity to have more collected, a solution can be found for the most unique type of geospatial project. If you need data or help figuring out what type to use for a project, contact us at MapMart.com.

Comments (0) Number of views (1877) Article rating: 4.0

Categories: ENVI Blog | Imagery Speaks

Tags:

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 (1628) Article rating: 5.0

Categories: IDL Blog | IDL Data Point

Tags:

12345678910Last

MONTHLY ARCHIVE

MOST POPULAR POSTS

AUTHORS

Authors

Authors

Authors

Authors

Authors

Authors

Authors

Authors

Authors

Authors

GUEST AUTHORS

Authors

Authors

Authors

Authors

Authors

Authors

Authors



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