7
Jul
2016
Author: Atle Borsholm
Resampling images is a very common operation in IDL, and it can happen both implicitly as well explicitly. Implicit resampling happens with IDLgrImage rendering. When the destination rendering area contains fewer pixels than the original image, then downsampling occurs. When the destination area is larger than the original image, then upsampling occurs. There are many options for upsampling algorithms. The simplest is a pure pixel replication to fill in the gaps. This is useful when there is a need to look closely at the original data. However, if the goal is to look for details in the scene that may be approaching the limits of the image resolution, then a more sophisticated resampling algorithm should be chosen instead. There are a few options that are commonly used. Bilinear interpolation and cubic spline interpolation are both options that are available with the CONGRID function in IDL. Lanczos and Lagrange resampling are two other options that are more computationally intensive. In the code below, I am showing an example comparing the Lanczos resampling kernel with bilinear and pixel replication. Lanczos resampling is often preferred because of its ability to preserve and even enhance local contrast, whereas bilinear tends have a blurring effect.
The results are shown here, starting with the original image, then the 4x zoomed area with pixel replication, then the 4x zoomed with bilinear interpolation, and finally the 4x zoomed with Lanczos convolution. The Lanczos convolution has the advantage of retaining good contrast while avoiding looking too pixelated.
Categories: IDL Blog | IDL Data Point
Tags:
6
Author: Jon Coursey
In a blog post I wrote a few months ago, I introduced the Instant Satellite Imagery Portal section of the Harris Geospatial Marketplace. This service offers a faster, less expensive alternative to purchasing DigitalGlobe imagery directly from a sale’s representative. However, at the time, it could still take up to 48 hours for the image to be delivered following purchase. This is no longer the case. Thanks to some savvy engineering by our technical team, the wait time for imagery download has been reduced to just a few minutes.
Often times, someone will contact us with a request for imagery, but will have a deadline that would generally be impossible to meet. And even when ordering directly from a provider, the rush delivery options could still mean up to a 48 hour delivery time. This service was created with these types of requests in mind. If someone is new to purchasing imagery, they may not know the typical delivery times for satellite captures. Additionally, the need for an image could have been realized towards the very end of a project, with too little an amount of time to actually receive one. Whatever the case, our Instant Satellite Imagery Portal hopes to offer a solution for such issues.
While in the middle of automating the Portal, it occurred to our technical team that “hey, this might just work for other products”. This idea recently came to fruition, with many of our other datasets becoming immediately downloadable following purchase. Currently, the following products have been included:
- NAIP Imagery
- CityOrtho Imagery
- GeoOrtho Imagery
- SRTM Elevation Data
- NED Elevation Data
- Harris Gap Filled Elevation Data
These datasets can be ordered via the same process on the Portal. For our Aerial Imagery and Digital Elevation Model sections, you basically have access to the immediate download option for anything with a price listed.
Although only a few products have been set up for this service so far, we are working on making this available for any data that does not require custom processing. This will make the Harris Geospatial Marketplace one of the quickest and easiest sources for downloadable data. Vector data and topographic maps will be available very soon, with hopes of adding more sources for satellite imagery and elevation data by the end of 2016. For any project facing a deadline that requires some sort of geospatial data set, the Marketplace will soon be the solution.
Categories: ENVI Blog | Imagery Speaks
30
Jun
Author: Austin Coates
Recently, I was working on a project that had some pretty strict requirements for the max size of the output .png file. This did not seem like a big deal initially, until I realized that it is extremely difficult to calculate the final size of a .png file prior to creating the image. Long story short, this led me to creating a brute force method for finding the optimal size of an output image.
As a side note, I know what you are thinking, “if it is so difficult, why don’t you just choose a couple of sizes, try those out and then create a linear model so that in the future you no longer have to deal with this?” Well I tried it and it was a no-go. As you can see from figures 1 and 2 it is not a simple linear relationship, nor is it an easily predictable relationship. The blue lines in figures 1 and 2 represent the file’s size as compared to the number of pixels in the image; whereas the red line is a simple linear regression. It is possible to see a pretty distinct deviation from the trend line.
Figure 1
Figure 2
Back to the task at hand: the brute force method. The method is relatively simple:
1. Save the image at its native size (if it is smaller than the max allowable size, you are done).
2. Reduce the original image size by one half and check the file size .
a. If it is larger than the allowable size, take one quarter of the original image dimensions, or in other words, one half of the distance between the last computed value that was too large (i.e. 50% of the initial image size) and the last smallest computed value (the smallest image size has not been computed yet, since this is the first iteration so just set it to zero).
b. If it is smaller than the allowable size, take three quarters of the original file size, or one half the distance between the last largest image dimensions computed (in this case the original image dimensions) and the last smallest image dimensions computed (i.e. one half the original).
3. Continue the pattern of taking one of the values either greater than or less than the halfway point that was just computed. Each time, make sure that the number of lines and samples is being computed using either integer or long values.
4. After each iteration, check to see if the same image dimensions have already been tested and if so, that is your final image size and you have found the optimal solution.
The example below should shed a little bit more light on how this is being performed.
;Start the application
e = ENVI()
;Open an input file
File = Filepath('qb_boulder_msi', Subdir=['data'], $
Root_Dir=e.Root_Dir)
Raster = e.OpenRaster(File)
;Get the task from the catalog ofENVITasks
Task = ENVITask('ISODATAClassification')
;Define inputs
Task.Input_Raster = Raster
;Run the task
Task.Execute
;Pull out the new ISOData raster
ISORaster = Task.OUTPUT_RASTER
;get the metadata
metadata = ISORaster.Metadata
;get the information about the new file
nb = ISORaster.nb
Raster_ns = ISORaster.ns
Raster_nl = ISORaster.nl
;set the filename for the output file
output_file = 'C:\Output\PNG_Test.png'
;Set the desired output size in bytes
output_size = 150000
;delete any old files with the same name
FILE_DELETE, output_file,/ALLOW_NONEXISTENT
;Save the image out as our initial test
;Get the task from the catalog of ENVITasks
ERTP_Task = ENVITask('ExportRasterToPNG')
ERTP_Task.INPUT_RASTER = ISORaster
;Define outputs
ERTP_Task.OUTPUT_URI = output_file
ERTP_Task.Execute
;Close the PNG Raster
ERTP_Task.Output_Raster.close
;Get the new File Size
file_size = ((file_info(output_file)).size) * 1.
;create a container for the maximum number of samples and lines
max_ns = Raster_ns
max_nl = Raster_nl
;create a container for the current number of samples and lines
nl = Raster_nl
ns = Raster_ns
;create a container for the minimum number of samples and lines
min_ns = 0
min_nl = 0
;Create a container for the 1D index value of the lower right corner
;of the image with regards to the maximum image size
loc_1d = []
;check to make sure the image is currently too large
if file_size gt output_size then begin
;check to see if the same image size has come up more than once. if so
;then stop iterating
while ((where(loc_1d eq (nl * Raster_ns +ns)))[0] eq -1) do begin
;record the image width and height in 1D
loc_1d = [loc_1d, nl * Raster_ns + ns]
;decide which side of the value should be divided in half
if (file_size gt output_size) then side = 'right'
if (file_size lt output_size) then side = 'left'
case side of
;the image is too large and must be reduced in size
'right' : begin
; determine the new size of the image (i.e.1/2 of the upper bound minus the current size)
ns = max_ns - ((max_ns - min_ns)* .5)
nl = max_nl - ((max_nl - min_nl)* .5)
; Delete the old image
; Shrink the original image
; Get the task from the catalog ofENVITasks
DRR_Task=ENVITask('DimensionsResampleRaster')
; Define inputs
DRR_Task.INPUT_RASTER =ISORaster
DRR_Task.DIMENSIONS=[ns,nl]
; Run the task
DRR_Task.Execute
ns = DRR_Task.Output_Raster.ns
nl = DRR_Task.Output_Raster.nl
ERTP_Task.INPUT_RASTER =DRR_Task.Output_Raster
; Define outputs
ERTP_Task.OUTPUT_URI =output_file
; Close the Resampled Raster
DRR_Task.Output_Raster.close
; Close the PNG raster
ERTP_Task.OUTPUT_Raster.close
; Get the new File Size
file_size = ((file_info(output_file)).size)* 1.
; Record the new dimensions as the upperbound
max_ns = ns
max_nl = nl
end
'left' : begin ; too small
; determine the new size of the image (i.e.1.5 of the upper bound minus the current size)
ns = ((max_ns - min_ns)* 1.5) + max_ns
nl = ((max_nl - min_nl)* 1.5 )+ max_nl
; Record the dimensions of the new image
dims = size(test_img,/DIMENSIONS)
ns = DRR_Task.OUTPUT_RASTER.ns
nl = DRR_Task.OUTPUT_RASTER.nl
; Save the new image and test the size
; Record the new file size
; Record the old dimension as the lowerbound
min_ns = max_ns
min_nl = max_nl
; Record the new dimension as the upperbound
endcase
endwhile
endif
;Close the intial raster
ISORaster.close
24
Author: Brian Griglak
You have to be careful when you add keywords to a routine because you want to forward them on to another routine that you invoke. It’s okay when the keywords are always present, but if they are optional output keywords, you may end up copying around data that the caller of the outermost routine never sees. Here is a short example of this:
pro layer2, FOO=foo, BAR=bar compile_opt idl2 print, 'Layer2, N_ELEMENTS(foo) = ' + StrTrim(N_ELEMENTS(foo),2) print, 'Layer2, N_ELEMENTS(bar) = ' + StrTrim(N_ELEMENTS(bar),2) print, 'Layer2, ARG_PRESENT(foo) = ' + StrTrim(ARG_PRESENT(foo),2) print, 'Layer2, ARG_PRESENT(bar) = ' + StrTrim(ARG_PRESENT(bar),2) end pro layer1, FOO=foo, BAR=bar compile_opt idl2 print, 'Layer1, N_ELEMENTS(foo) = ' + StrTrim(N_ELEMENTS(foo),2) print, 'Layer1, N_ELEMENTS(bar) = ' + StrTrim(N_ELEMENTS(bar),2) print, 'Layer1, ARG_PRESENT(foo) = ' + StrTrim(ARG_PRESENT(foo),2) print, 'Layer1, ARG_PRESENT(bar) = ' + StrTrim(ARG_PRESENT(bar),2) layer2, FOO=foo, BAR=bar end
When I call layer1 with no keywords, we see how both FOO and BAR are not present in layer1, but they are in layer2:
IDL> layer1 Layer1, N_ELEMENTS(foo) = 0 Layer1, N_ELEMENTS(bar) = 0 Layer1, ARG_PRESENT(foo) = 0 Layer1, ARG_PRESENT(bar) = 0 Layer2, N_ELEMENTS(foo) = 0 Layer2, N_ELEMENTS(bar) = 0 Layer2, ARG_PRESENT(foo) = 1 Layer2, ARG_PRESENT(bar) = 1
This is because the interpreter only looks at the way layer1 invokes layer2, and has to assume that layer1 is using FOO and BAR as output keywords. While it is certainly possible that the interpreter could inspect how the foo and bar variables are used, and try to optimize them away when it realizes that these variables aren’t used elsewhere inside layer1 and weren’t passed into it, that would have a major impact on performance. Compiled languages can do this with an optimizing linker, but interpreted languages like IDL can afford to take the time to do that every time they run a line of code.
One way to interpret the values of N_ELEMENTS() and ARG_PRESENT() is shown in this table:
If I call layer1 with FOO and BAR set to a literal and an undefined variable, we see different output:
IDL> layer1, FOO=1, BAR=b Layer1, N_ELEMENTS(foo) = 1 Layer1, N_ELEMENTS(bar) = 0 Layer1, ARG_PRESENT(foo) = 0 Layer1, ARG_PRESENT(bar) = 1 Layer2, N_ELEMENTS(foo) = 1 Layer2, N_ELEMENTS(bar) = 0 Layer2, ARG_PRESENT(foo) = 1 Layer2, ARG_PRESENT(bar) = 1
Here we see how FOO is an input to layer1, and BAR is an output. The more interesting part is in layer2, where FOO in an in/out keyword while BAR is still only an output keyword. The situation changes a little when I call layer1 with a defined variable for one of the keywords:
IDL> b2 = 2 IDL> layer1, FOO=1, BAR=b2 Layer1, N_ELEMENTS(foo) = 1 Layer1, N_ELEMENTS(bar) = 1 Layer1, ARG_PRESENT(foo) = 0 Layer1, ARG_PRESENT(bar) = 1 Layer2, N_ELEMENTS(foo) = 1 Layer2, N_ELEMENTS(bar) = 1 Layer2, ARG_PRESENT(foo) = 1 Layer2, ARG_PRESENT(bar) = 1
This time FOO is an input to layer1 and in/out to layer2, but BAR is an in/out keyword to both layer1 and layer2.
Now you may be wondering what the point of this whole analysis is. It matters when an output keyword takes a lot of time to build or space to store, and is incorrectly identified as being present. Imagine that BAR uses a gigabyte of memory, if the user calls layer1 without BAR, then layer2 will allocate that memory and return it to layer1, but it gets thrown away when layer1 returns to the caller.
How do we defensively implement the code to prevent this waste of time and space? Unfortunately the best solutions I’ve come up with are to add some logic to layer1 to conditionally forward keywords to layer2. It matters what type of keyword we are talking about here, input vs output vs in/out. Pure input keywords can be blindly forwarded, while output and in/out need the extra logic.
pro layer2, INPUT=input, OUTPUT=output, INOUT=inout compile_opt idl2 print, 'Layer2, N_ELEMENTS(input) = ' + StrTrim(N_ELEMENTS(input),2) print, 'Layer2, N_ELEMENTS(output) = ' + StrTrim(N_ELEMENTS(output),2) print, 'Layer2, N_ELEMENTS(inout) = ' + StrTrim(N_ELEMENTS(inout),2) print, 'Layer2, ARG_PRESENT(input) = ' + StrTrim(ARG_PRESENT(input),2) print, 'Layer2, ARG_PRESENT(output) = ' + StrTrim(ARG_PRESENT(output),2) print, 'Layer2, ARG_PRESENT(inout) = ' + StrTrim(ARG_PRESENT(inout),2) output = 'output' inout = ISA(inout) ? inout+1 : 'new inout' end pro layer1, INPUT=input, OUTPUT=output, INOUT=inout compile_opt idl2 print, 'Layer1, N_ELEMENTS(input) = ' + StrTrim(N_ELEMENTS(input),2) print, 'Layer1, N_ELEMENTS(output) = ' + StrTrim(N_ELEMENTS(output),2) print, 'Layer1, N_ELEMENTS(inout) = ' + StrTrim(N_ELEMENTS(inout),2) print, 'Layer1, ARG_PRESENT(input) = ' + StrTrim(ARG_PRESENT(input),2) print, 'Layer1, ARG_PRESENT(output) = ' + StrTrim(ARG_PRESENT(output),2) print, 'Layer1, ARG_PRESENT(inout) = ' + StrTrim(ARG_PRESENT(inout),2) if (ARG_PRESENT(output)) then begin if (ARG_PRESENT(inout)) then begin layer2, INPUT=input, OUTPUT=output, INOUT=inout endif else begin layer2, INPUT=input, OUTPUT=output endelse endif else begin if (ARG_PRESENT(inout)) then begin layer2, INPUT=input, INOUT=inout endif else begin layer2, INPUT=input endelse endelse end
The big drawback with this approach, as you can expect, is that it is exponential in the number of keywords and quickly becomes untenable. We need a solution that is O(n) instead of O(2n), and I’ve come up with two completely different solutions which each have their pros and cons. First the simpler solution, which uses the EXECUTE () function:
pro layer1, INPUT=input, OUTPUT=output, INOUT=inout compile_opt idl2 print, 'Layer1, N_ELEMENTS(input) = ' + StrTrim(N_ELEMENTS(input),2) print, 'Layer1, N_ELEMENTS(output) = ' + StrTrim(N_ELEMENTS(output),2) print, 'Layer1, N_ELEMENTS(inout) = ' + StrTrim(N_ELEMENTS(inout),2) print, 'Layer1, ARG_PRESENT(input) = ' + StrTrim(ARG_PRESENT(input),2) print, 'Layer1, ARG_PRESENT(output) = ' + StrTrim(ARG_PRESENT(output),2) print, 'Layer1, ARG_PRESENT(inout) = ' + StrTrim(ARG_PRESENT(inout),2) cmd = 'layer2' if (N_ELEMENTS(input)) then begin cmd += ', INPUT=input' endif if (ARG_PRESENT(output)) then begin cmd += ', OUTPUT=output' endif if (N_ELEMENTS(inout) || (ARG_PRESENT(inout)) then begin cmd += ', INOUT=inout' endif !null = EXECUTE(cmd) end
This approach builds up a command string to execute, starting with the name of the layer2 routine it invokes. It then conditionally appends the keywords to the command string, using N_ELEMENTS() and/or ARG_PRESENT(), depending on what type of keyword it is. This doesn’t add a lot of new code, but it takes a performance hit because EXECUTE () has to compile the command string and then execute it. One might point out that using CALL_PROCEDURE() is more efficient, but it doesn’t support output and in/out keywords so we have to use EXECUTE().
The other approach avoids the performance hit of EXECUTE(), but it is a bit more involved, requiring the creation of a new wrapper function and the use of SCOPE_VARFETCH(), which some people are leery to use.
function layer2Wrapper, _REF_EXTRA=extra compile_opt idl2 layer2, _EXTRA=extra if (N_ELEMENTS(extra) eq 0) then return, Hash() retVal = Hash() foreach key, extra do begin retVal[key] = SCOPE_VARFETCH(key, /REF_EXTRA) endforeach return, retVal end pro layer1, INPUT=input, OUTPUT=output, INOUT=inout compile_opt idl2 print, 'Layer1, N_ELEMENTS(input) = ' + StrTrim(N_ELEMENTS(input),2) print, 'Layer1, N_ELEMENTS(output) = ' + StrTrim(N_ELEMENTS(output),2) print, 'Layer1, N_ELEMENTS(inout) = ' + StrTrim(N_ELEMENTS(inout),2) print, 'Layer1, ARG_PRESENT(input) = ' + StrTrim(ARG_PRESENT(input),2) print, 'Layer1, ARG_PRESENT(output) = ' + StrTrim(ARG_PRESENT(output),2) print, 'Layer1, ARG_PRESENT(inout) = ' + StrTrim(ARG_PRESENT(inout),2) keywords = Hash() if (N_ELEMENTS(input)) then begin keywords['input'] = input endif if (ARG_PRESENT(output)) then begin myOutput = 0 keywords['output'] = myOutput endif if (N_ELEMENTS(inout)) then begin keywords['inout'] = inout endif else if (ARG_PRESENT(inout)) then begin myInout = 0 keywords['inout'] = myInout endif outKeys = layer2Wrapper(_EXTRA=keywords.ToStruct()) if (outKeys.HasKey('OUTPUT')) then output = outKeys['OUTPUT'] if (outKeys.HasKey('INOUT')) then inout = outKeys['INOUT'] end
This version of layer1 conditionally adds the keywords to a Hash, which is converted to a struct when calling the new wrapper function. By assigning the struct to the _EXTRA keyword during invocation, the keywords are properly mapped. One will note that the variables put into the Hash depend on whether it is an input or output keyword. For input keywords, I just used the variables assigned to those keywords in layer1’s signature, but for output keywords I had to create a local variable and put that in the Hash. The reason for this is that IDL structs cannot contain !NULL or undefined values, so I had to assign them a value. In this case I used 0, but depending on the expectations of layer2 a special “not set” value may need to be used.
The new layer2Wrapper() function is worth examining more closely. It is declared using _REF_EXTRA, so that we get pass by reference semantics instead of pass by value. That _REF_EXTRA bag is simply passed along to layer2, using the _EXTRA keyword like we’re supposed to. After calling layer2, the wrapper then uses SCOPE_VARFETCH() with the /REF_EXTRA keyword to grab each keyword’s value and put it into a Hash that is returned to its caller. It’s this SCOPE_VARFETCH() trick that necessitates the creation of this wrapper routine. Once layer2Wrapper returns that Hash to layer1, we then copy any of the existing values out of the Hash into the appropriate output keyword variables.
Both solutions may seem like a lot of work for what seems like a minor inconvenience, but there are real world cases where we could waste a lot of time and/or memory on keywords that only exist for parts of the callstack and are never specified by the original caller. A prime example of this is in ENVIRaster::GetData(). This method always returns an array of pixel values, but there is the optional PIXEL_STATE keyword that can be used to retrieve a byte array of the same dimensions that indicates which pixels are good and which are bad. Without using a solution like the ones I’ve presented, we would end up with the PIXEL_STATE array being allocated and calculated every time a user calls GetData(), whether they specified that keyword or not.
Author: Zachary Norman
Believe it or not, there is a new market for precision agriculture on our big red neighbor. However, this is not your conventional agriculture analysis, but rather an innovative use of the crop counter tool in the Precision Ag Toolkit.
Our crop of choice that we will be counting is not alive, but rather lifeless scars on the surface of Mars. Because Mars has such a thin atmosphere, about about 0.6% of Earth's mean sea level atmospheric pressure, smaller meteorites don't necessarily burn up in the atmosphere like they do in our night sky. Instead, the surface or Mars is covered in small, round craters which present ideal objects to find with the crop counter tool. These craters range in size from a few meters in diameters to kilometers in size. For this case, I'm focusing on caters that are about 40-50 meters in diameter.
The coolest thing about all of the international space missions is that most of the data collected is easy to find and freely available. For Mars, I actually decided to use a DEM (Digital Elevation Model) that was derived from stereo pairs with the HiRISE (High Resolution Imaging Science Experiment) sensor which is on board the Mars Reconnaissance Orbiter. The HiRISE sensor takes images with up to .3 meter spatial resolution, which allows for lots of detail in the images that are captured. The resulting DEM I used had .99 meters per pixel for spatial resolution. The exact DEM I decided to go with can be found here
Here is a screenshot of a portion of the image that I used to count craters:
I chose this image because it had lots of uniform, round craters. After exploring the craters a bit, I decided to look for craters between 35 and 45 meters in size. Before I could count the craters, I first needed to figure out how I wanted to count them.
For the crop counter, the best way to get accurate results is to preprocess your data and do some thresholding. The first reason for this is that thresholding helps prevent false positives from being counted. Secondly, thresholding drastically increases the speed for finding crop centers. This speed increase is because, when pixels are Nan's, the crop counter skips those pixels and decreases processing time. For crater counting, this preprocessing step just becomes a matter of how to get the craters to pop out of the image.
After some thinking, I decided to isolate the crater pits themselves I would use some smoothing, band math, and then thresholding. Here are the exact steps I took to count craters:
The kernel size for smoothing the data was 91 because we were looking for craters between 35 and 45 meters in diameter. Below you can see the smoothed DEM with a similar color table as above:
With the smoothed DEM, the crater features we are looking for should be pretty much erased. This means that taking the difference between our smoothed image and our original image will return the approximate crater depth. Because the crop counting tool works with bright images we can then use the crater depth image with the plant counter after one more step is applied. Below you will see the the difference between the DEM and the smoothed image. The green areas have a higher value (height in meters) and the red areas have a lower value.
This step is very important to reduce false positives for counting craters because our difference image is quite noisy in some places. This step eliminates smaller craters and noise in our image by saying the crater depth needs to be greater than 'x' meters. I used .75 meters for this step which really isolates the craters I was looking for. To illustrate how this helps isolate the craters we are looking for, you can see the image below. This screenshot is of the original DEM with the image difference shown after a threshold has been applied.At this point you can see that we have isolated most of the craters with a few spots of noise here and there.
Step Three: Using the Crop Counter!
Now that we have preprocessed our data and have a great image to count craters in, it's time to use the crop counter. Here are the setting I used in the crop counter to find the crater locations:
CropCountTask = ENVITask('AgCropCount') CropCountTask.INPUT_RASTER = ThresholdRaster CropCountTask.WIDTH_IN_METERS = 0 CropCountTask.PERCENT_OVERLAP = .01 CropCountTask.SMOOTH_DATA = 1 CropCountTask.SEARCH_WIDTH = [35,45,1] CropCountTask.MAX_TILE_SIZE = [400,400]
After running the task, here is what our results look like. This image is a grayscale version of the DEM with the crop locations shown as cyan circles.
As you can see, the crop counting tool does a pretty good job of finding the actual crater locations!
Base 60 encoding of positive floating point numbers in IDL
First Valley Shadow Detection
Machine Learning Training for Automatic Target Detection
Sign Up for News & Updates: Stay informed with the latest news, events, technologies and special offers.