I find myself drawing bounding boxes around things a lot. I don’t know why I do it so much, but for whatever reason I do, and as of late I wanted to up my bounding box game. In the past, I have simply used the global min and max in both the x and y directions to get the coordinates to form the bounding box; however, this is not always the most elegant solution. For example, when my data follows somewhat of a linear trend, I am left with ample white space not filled by any valuable information

Figure 1: Simple Bounding Box

Figure 2: Minimum Area Bounding Box

This got me thinking, why am I not simply drawing a bounding box around only the data? Sounds great, right? The only problem was I had no idea how to do this. Luckily, there is this thing called the internet and it has vast stores of information and ideas to pull from. I found a very elegant solution by Jesse Buesking on stackoverflow.com which was posted on November 9, 2015. The solution was written in Python which I converted to IDL. My goal in posting this is to show an awesome way to draw a bounding box and also an example of translating between IDL and Python.

function bounding_box, pts = pts, plot_results = plot_results
 compile_opt IDL2
;Get the x and y coordinates
xs = pts[0,*]
ys = pts[1,*]
 
;Find the bounding points
Triangulate, xs, ys, triangles, hull, CONNECTIVITY=CONNECTIVITY


;order hull points in a [2,n] array    
 hull_points = [[xs[hull]]##1,[ys[hull]]##1]
;calculate edge angles
edges = hull_points[*,1:-1] - hull_points[*,0:-2]
angles = atan(edges[1,*], edges[0,*])

pi2 = !DPI/2.
 
angles = abs(angles - floor(angles / pi2) * pi2)
angles = angles[sort(angles)]
angles = angles[UNIQ(angles)]


;find rotation matrices 
rotations = transpose([[cos(angles)],[cos(angles-pi2)],[cos(angles+pi2)],[cos(angles)]])
rotations = REFORM(rotations, [2,2,n_elements(angles)])
 
;apply rotations to the hull 
rot_points = fltarr( n_elements(hull_points)/2, 2, n_elements(angles))
size_rot = size(rotations)

for group = 0 , size_rot[3]-1 do begin   
  for row = 0 , size_rot[2]-1 do begin
    rot_points[*,row,group] = TRANSPOSE(rotations[*,row,group]) # hull_points
  endfor
endfor

;find the bounding points
min_x =  min(rot_points[*,0,*],DIMENSION=1, /NAN)
max_x =  max(rot_points[*,0,*],DIMENSION=1, /NAN)
min_y =  min(rot_points[*,1,*],DIMENSION=1, /NAN)
max_y =  max(rot_points[*,1,*],DIMENSION=1, /NAN)

;find the box with the best area
areas = (max_x - min_x) * (max_y - min_y)
min_val = min(areas, best_idx)

;return the best box
x1 = max_x[best_idx]
x2 = min_x[best_idx]
y1 = max_y[best_idx]
y2 = min_y[best_idx]
r = rotations[*,*,best_idx]

rval = fltarr(2,4)
rval[*,0] = TRANSPOSE(TRANSPOSE([x1, y2]) # transpose(r))
rval[*,1] = TRANSPOSE(TRANSPOSE([x2, y2]) # transpose(r))
rval[*,2] = TRANSPOSE(TRANSPOSE([x2, y1]) # transpose(r))
rval[*,3] = TRANSPOSE(TRANSPOSE([x1, y1]) # transpose(r))
 
;display results  
if KEYWORD_SET(plot_results) then begin
  p = SCATTERPLOT(xs,ys, SYM_COLOR='Red', SYM_FILL_COLOR='Red', SYM_FILLED=1,$
    XRANGE=[min(rval[0,*])-1,max(rval[0,*])+1], YRANGE=[min(rval[1,*])-1,max(rval[1,*])+1])

  p = POLYGON(rval, /FILL_BACKGROUND, $
   FILL_COLOR="light steel blue", PATTERN_ORIENTATION=45, $
   PATTERN_SPACING=4, /DATA)
endif


return, rval


end




        

Source of original Python code : http://stackoverflow.com/questions/13542855/python-help-to-implement-an-algorithm-to-find-the-minimum-area-rectangle-for-gi/33619018#33619018 

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

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