Finding the Nth ordered element in a large array

Author: Atle Borsholm

A common task when working with large arrays is to find the Nth array value in the ordered array. This can be useful for finding the Nth smallest or largest pixel value, as well as for statistical analysis of floating point data samples, (i.e. find the 95% percentile or similar). The shortest IDL code for finding the Nth value in the ordered sequence is only 2 lines of actual code. Here is a short function that accomplishes this:



; Returns the Nth number in theordered sequence


function ordinal_1, array, N

 compile_opt idl2,logical_predicate

 s = sort(array)

 return, array[s[N]]



However, because sort is an expensive computation, it runs fairly slow, especially, when the array gets larger. I did the following time test.

IDL> a = total(read_image(filepath('ohare.jpg',subdir=['examples','data'])),1)

IDL> tic & x = ordinal_1(a, 123456) & toc & print, x

% Time elapsed: 3.7200000 seconds.


In my case it took 3.72 seconds to find the 123456th smallest array element. The MEDIAN function in IDL, returns the central element in the ordered sequence without doing a full sorting. It is much faster than sorting, because it doesn't have to keep track of all elements and their ordered positions. It only cares about element N/2 in the ordered array. In the following example, repeated calls to MEDIAN and reducing the array size in half every iteration, is used to find the Nth element in the ordered sequence. The code is much longer than the code above, but it does end up running faster:


; Returns the Nth number in theordered sequence.


; Uses repeated median.


function ordinal_2, array, N

 compile_opt idl2,logical_predicate

 na = n_elements(array)

 type =size(array, /type)

 target_index = N

 tmp = arg_present(array) ? array : temporary(array)

 ntmp = na

 while ntmp ne target_index do begin

   ntmp = n_elements(tmp)

   val = fix(median(tmp), type=type)

   if target_index gt ntmp/2 then begin

     tmp = tmp[where(tmp gt val, count)]

     target_index -= ntmp-count

   endif else if target_index lt ntmp+1/2 then begin

     tmp = tmp[where(tmp lt val, count)]

   endif else break

   if target_index lt 0 then break

   if target_index ge count then break

   if target_index eq 0 then begin

     val = min(tmp)



   if target_index eq count-1 then begin

     val = max(tmp)




 return, val


This is the same time test as with the short code:

IDL> tic & x = ordinal_2(a, 123456) & toc & print, x

% Time elapsed: 0.57999992 seconds.



As can be seen here, the time saving is significant, it goes from 3.72 to 0.58 seconds, and as the array grows larger, the savings can get more significant. This function works for numeric data types such as floating point and integer arrays.

Comments (2) Number of views (205) Article rating: No rating

Categories: IDL Data Point






Dimensional Analysis, Unit Conversions, and Physical Constants

Author: Benjamin Foreback

Physical Constants

When I was in college, I memorized most of the common physical constants, such as Plank's Constant or Earth's mass. To be honest, today I don't remember most of these without looking them up. Although this isn't too difficult to do in the era of the internet where I can copy and paste the search result of "Plank's Constant" from my browser into the IDL editor, IDL makes it much easier by storing many of these physical constants in the !CONST system variable. For example, if I am writing code that calculates gravitational forces using the Gravitation Constant, I can simply use !CONST.G. More convenient than just saving me a few seconds of typing a long number as a variable, however, is the fact that IDL automatically stores these values as double precision, ensuring maximum accuracy each time they are used.

Unit Conversion in IDL

Another helpful tool when developing code that involves processing dimensional data is IDLUNIT, which can perform unit conversion, calculate simple mathematical expressions, and much more. Similar to !CONST, all values and calculations are kept in double precision for maximum accuracy.

Here is an example of how you may convert PSI to Pascals:

PsiToPa = (IDLUnit('psi -> Pa')).Quantity

PaValue = psiValue * psiToPa

IDLUNIT even allows a user to add custom units, which is helpful when using non-standard units. For example, how many Olympic sized swimming pools does it take to fill Sydney Harbour? Let's define a unit of 'pool' as 50 x 25 x 2 meters = 2500 m3, and let's define Sydney Harbour as 494 Gigaliters (source: Wikipedia: List of unusual units of measurement).

IDLUnit.AddUnit, 'pool', '2500 m^3', PLURAL='pools'

IDLUnit.AddUnit, 'Sydharb', '494 Gl', PLURAL='Sydharbs'

IDLUnit('1 Sydharb -> pools')

IDL tells us that 197,600 pools could fill Sydney Harbour.

IDL's help page for the IDLUNIT topic contains many additional examples.

Unit Converter Widget

In addition to performing unit conversions in code or at the command line, IDL has an interactive unit converter. To launch this converter, run 'idl_convert' or select 'IDL Converter' from the IDL Workbench's Macros menu. Notice that the custom units defined in the example above are included in this tool.


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

Categories: IDL Data Point





Gum Drop Plot

Author: David Starbuck

On this webpage advertising IDL graphics, there is a image that displays a plot made out of a series of colorful spheres. Based on my research, this plot was generated a long time ago most likely using IDL Object Graphics.

I wanted to generate a similar plot using the IDL 8 Graphics Functions. To do this, I used the SCATTERPLOT3D routine, the ORB object (provided with the IDL distribution but not documented), and the POLYLINE routine.  An example of the output is shown below:

The data used to produce this plot is randomly generated. Therefore, if you run the code, the output will look different each time. The code used to generate this plot is shown below:

pro gum_drop

  compile_opt idl2
  ;Generate some random data to
  x = RANDOMU(seed, 10)
  y = RANDOMU(seed, 10)
  z = RANDOMU(seed, 10)
  ;Draw an initail scatter
  ;plot with symbolsize of 1    
  scat_plot = scatterplot3d(x,y,z,RGB_TABLE=2, $
              SYM_OBJECT=orb(),$ ;use an orb object as symbol
              SYM_SIZE=1, $ ;Set symbol size to 1
              MAGNITUDE=z, $ ;change color with Z value
              /SYM_FILLED, clip=0,$ ;Fill symbols and no clipping
              xticklen=0, yticklen=0, zticklen=0, $ ;remove ticks
              xsubticklen=0, ysubticklen=0, zsubticklen=0, $ ;remove ticks
              xmajor=5, ymajor=5, $ ;Only use 5 ticks on each axis
              xrange=[0,1], yrange=[0,1],$ ;force the x and y range
              ASPECT_RATIO=1.0,$  ;Don't distort the image
              DEPTH_CUE=[0,4], $ ;Make things farther away fade
              AXIS_STYLE=2, $  ;Make the axis a box
              background_color = 'light yellow')
  ;Draw 9 more plots where the symbol size
  ;changes each plot
  for ind = 2L, 10 do begin
   z = RANDOMU(seed, 10)
   scat_plot_loop = scatterplot3d(x,y,z,$
      RGB_TABLE=2, SYM_OBJECT=orb(),SYM_SIZE=ind/2, $

  ;Generate polygons to create a grid
  ;on the Z-Y and Z-X planes
  x = [0.25,0.25,0.5,0.5,0.75,0.75]
  y = [0.999,0.999,0.999,0.999,0.999,0.999]
  z = [0.00,1.00,0.0,1.0,0.00,1.00]

  ;Connect every 2 points in polygon data
  ;with lines using the CONNECTIVITY keyword
  con = [2,0,1,2,2,3,2,4,5]

  poly0 = polyline(x,y,z,/DATA,CONNECTIVITY=con)
  temp = x
  poly1 = polyline(x,y,z,/DATA,CONNECTIVITY=con)
  temp = z
  z = y
  y = temp
  poly2 = polyline(x,y,z,/DATA,CONNECTIVITY=con)
  temp = x
  poly3 = polyline(x,y,z,/DATA,CONNECTIVITY=con)
  ;Remove the axis from the front of
  ;the plot
  ax = scat_plot.AXES

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

Categories: IDL Data Point

Tags: (New) Graphics




DebuggerHelper - A Handy Debugging Class for IDL Developers

Author: Jim Pendleton

Remember our old friend the IDL PRINT procedure?

In addition to faithfully showing us our derived output, it has also bailed us out of a lot of sticky code development problems over the years when Workbench breakpoints alone didn't suffice. Yes, sometimes bugs end up even in my code, put there by cosmic rays, random keystrokes and bouts of brain freeze.

PRINT does have its drawbacks as a debugging utility, however.  Here are four of them.

  1. You have to type "PRINT,". And then some. That's a lot of typing.  (See Implied Print.)
  2. Unless you develop applications for an audience of one, you generally need to eventually locate and comment out or remove all your debugging PRINT statements. If you choose otherwise, your live demo in front of your customer, manager, or thesis committee may show questionable comments in the IDL Console tab. Don't ask me how I know this.
  3. When your workbench or command line window is closed as a result of willful misconduct (see: bugs), you'll lose any PRINT output that could have helped you debug the problem causing the needed information to disappear in a puff of illogic.
  4. If you distribute IDL Runtime applications, as we frequently produce in the Exelis VIS Professional Services Group, the PRINT statement is of little or no value because there is no console output. You need a way to get useful information from your end user, to supplement the standard description of "IDL crashed!"

Over the course of a number of years a design pattern emerged in my application development. I found I really wanted a number of useful features

  1. It has to be easy to use and not obstruct the real work my application must perform
  2. Debugging information should display in a widget separate from the IDL Console, if desired
  3. Debugging output should be written to an auto-generated text file with a unique name, if desired
  4. The command to produce output should be simpler to type than "PRINT,"
  5. The utility should be a class that I can instantiate or from which I can subclass
  6. The utility should automatically produce time-stamping and contextual information
  7. I should be able to enable debugging without modifying my source code

Now I will hook you in with the simplicity of the resulting IDL class, one that provides all this capability but is only about 150 lines long. The code is listed at the bottom of the article.  Copy, paste and compile to play along.

Create the DebuggerHelper object and toggle on the debugging

IDL> d = debuggerhelper(/on)

This creates both a log file and a debug text widget. It can't get much easier to instantiate, thereby satisfying my first three requirements.

Generate a line of debugging output.

IDL> d.o, 'Hello bug'

Requirement 4?  Check.

I'm only showing the output in the IDL Console here, but I'm also writing information to a stand-alone text widget AND a text file. The displayed data shows the class from which the output is produced, which may be a subclass of DebuggerHelper, along with the context of the call. In this case we're at the IDL command line so this is $MAIN$.

If it were called from within an IDL function or procedure, that routine's name would replace $MAIN$, as shown later.

I can turn on TIC and TOC management to help test performance. This wasn't an original requirement, but I added it to my class while writing this article since it was incredibly simple to do so.

IDL> d.tic, 'Timer 1'
IDL> d.o, 'How long did it take me to type this?', toc = 'Timer 1'
(DEBUGGERHELPER) $MAIN$, Performance timer 1 dt=21.950000:How long did it take me to type this?
IDL> d.o, 'That was not very fast.'
(DEBUGGERHELPER) $MAIN$:That was not very fast.

I can have multiple timers defined by name to manage different code sections, if I desire.

Obviously, TIC and TOC aren't going to be as precise for high-performance testing when embedded in debugging statements since the act of watching and logging perturbs the system. However, they're still useful for cases where the operations under test take significantly more time than the debug I/O.

I can programmatically toggle my debugging on and off, which means I don't need to edit my source code, checking another item off my requirements.

IDL> d.off
IDL> d.o, 'Do not criticize my typing speed'
IDL> d.on
IDL> d.o, 'Sorry about that'
(DEBUGGERHELPER) $MAIN$:Sorry about that

The debug log file is automatically generated with a unique name, satisfying another requirement.  Let's see where it was put.

IDL> d.file

In this implementation, the log file is created relative to the location of the source .pro file (or SAVE file) by utilizing ROUTINE_FILEPATH, but you can follow your bliss when modifying the source code, or when subclassing from the original implementation.

Even after the object is destroyed or IDL exits, I can still view my debug information unless I proactively delete the file.

IDL> file = d.file
IDL> obj_destroy, d
IDL> b = strarr(file_lines(file))
IDL> openr,lun,file,/get_lun & readf, lun, b & free_lun, lun
IDL> b
[        11.310, 1403199708.570] (DEBUGGERHELPER) $MAIN$:Hello bug
[        66.141, 1403199763.401] (DEBUGGERHELPER) $MAIN$, Timer 1 dt=21.950000:How long did it take me to type this?
[        84.052, 1403199781.312] (DEBUGGERHELPER) $MAIN$:That was not very fast.
[       119.896, 1403199817.156] (DEBUGGERHELPER) $MAIN$:Sorry about that

Notice that the log file contains additional time stamp information that was not printed to the IDL Console.

I like to show two times. The first is the elapsed time in seconds since the object was created.  It's easier for me to do elapsed-time math in my head with numbers formatted this way. The second time is SYSTIME(1) output. If I have multiple debug objects in an application, I can later merge the log files and sort messages by time. For example I may want to sequence device communication protocols used in different parts of my application, or even among individual, stand-alone applications or IDL_IDLBridge processes.

I may want to ask a customer running my application to turn debugging on so that internal statements I'd built into the code will come to life. How do we fully address the final requirement, that I should not need to change the source code to enable debugging from a Runtime IDL application?

GETENV or COMMAND_LINE_ARGS offer two alternatives. For example, I may ask a client to create an environment variable in their operating system, whose name I have hardwired into the code, before starting my IDL Runtime application. During execution, I always create the debug object but I toggle the debugging statements on or off based on the environment variable setting, for example

toDebug = getenv('DEBUG_IDL_APP') eq 'YES'

d = debuggerhelper(on = toDebug, no_file = ~toDebug)

Here, if the user has not created an environment variable named DEBUG_IDL_APP or they have not populated that environment variable with the string "YES", then the debug helper object will be created, but it will be silent.

Even the presence or absence of a file placed in a particular directory can be used for this purpose, along with a corresponding FILE_TEST call in your code.

Adding the debug object to another class, or subclassing from it, will provide additional feedback. Here, we will create a DebuggerHelper object as a member variable within a new "MyObj" class.

IDL> .run
- function myobj::init
- return, 1
- end
% Compiled module: MYOBJ::INIT.
IDL> .run
- function myobj::init
- self.d = debuggerhelper(/on)
- return, 1
- end

IDL> .run
- pro myobj::yelp
- self.d.o, 'Within myobj'
- end
% Compiled module: MYOBJ::YELP.

Create an object of this class and call the method that generates debugging output

IDL> m = myobj()
IDL> m.yelp

Notice that the class and method names from which the debug output was generated are now displayed in my output as well.

Here's the source IDL for DebuggerHelper. I think it's straightforward, but leave a comment if you have questions. I'd enjoy the feedback, and I'd get a better sense if only web crawlers read our blogs.

Pro DebuggerHelper::_ConstructDebugWidget, Group_Leader
   Compile_Opt IDL2
   self._DebugTLB = Widget_Base(/Column, $
      Title = Obj_Class(self) + ' Debug', $
      Group_Leader = Group_Leader)
    DebugText = Widget_Text(self._DebugTLB, Value = '', $
      UName = 'DebugText', $
   XSize = 140, YSize = 24, /Scroll)
   Widget_Control, self._DebugTLB, /Realize
Pro DebuggerHelper::_CreateDebugFile
   Compile_Opt IDL2
   LogDir = FilePath('', $
      Root = File_DirName(Routine_Filepath()), $
      SubDir = ['logs'])
   If (~File_Test(LogDir)) then Begin
      File_MkDir, LogDir
      File_ChMod, LogDir, /A_Read, /A_Write
   self._DebugFile = FilePath(Obj_Class(self) + $
      StrTrim(ULong(SysTime(1)), 2) + '.log', $
      Root = LogDir)
   If (File_Test(self._DebugFile)) then Begin
      File_Delete, self._DebugFile
   OpenW, DebugLUN, self._DebugFile, /Get_LUN
   self._DebugLUN = DebugLUN
   File_ChMod, self._DebugFile, /A_Read, /A_Write
Pro DebuggerHelper::Cleanup
   Compile_Opt IDL2
   If (self._DebugLUN ne 0) then $
      Free_LUN, self._DebugLUN
   If (Widget_Info(self._DebugTLB, /Valid_ID)) then $
      Widget_Control, self._DebugTLB, /Destroy
Pro DebuggerHelper::DebugOutput, Output, $
   No_Print = NoPrint, $
   Up = Up, $
   Toc = Toc
   Compile_Opt IDL2
   If (~self._DebugOn) then Begin
   Elapsed = ''
   If (Toc ne !null) then Begin
      If (Size(Toc, /TName) eq 'STRING') then Begin
         Elapsed = ', ' + Toc + ' dt=' + $
            StrTrim(Toc(self._DebugClocks[Toc]), 2)
      EndIf Else If (Keyword_Set(Toc)) then begin
         Elapsed = ', dt=' + StrTrim(Toc(), 2)
   ThisClass = Obj_Class(self)
   Routines = (Scope_Traceback(/Structure)).Routine
   Result = '(' + ThisClass + ') ' + $
      Routines[N_elements(Routines) - $
      (2 + Keyword_Set(Up))] + $
      Elapsed + ':' + Output
   If (~Keyword_Set(NoPrint)) then Begin
      Print, Result, Format = '(a)'
   If (~Widget_Info(self._DebugTLB, /Valid_ID)) then Begin
   Now = SysTime(1)
   DT = '[' + $
      String(Now - self._DebugCreationTime, $
         Format = '(f14.3)') + ', ' + $
         String(Now, Format = '(f14.3)') + $
         '] '
   DebugText = Widget_Info(self._DebugTLB, $
      Find_by_UName = 'DebugText')
   DebugYSize = (Widget_Info(DebugText, /Geometry)).YSize
   Widget_Control, DebugText, Get_UValue = NLines
   If (N_elements(NLines) eq 0) then Begin
      NLines = 0L
   Line = DT + Result
   If (self._DebugLUN ne 0) then Begin
      PrintF, self._DebugLUN, Line
      Flush, self._DebugLUN
   EndIf Else Begin
      Print, Line, Format = '(a)'
   If (NLines gt self._DebugWindowMaxLines) then Begin
      Widget_Control, DebugText, Get_Value = Old
      Lines = Old[self._DebugWindowMaxLines/2:*]
      Old = 0
      NLines = N_elements(Lines)
      Widget_Control, DebugText, $
         Set_Value = [Temporary(Lines), Line]
      Widget_Control, DebugText, $
         Set_UValue = NLines + 1L
   EndIf Else Begin
      Widget_Control, DebugText, $
         Set_Value = Line, /Append
      Widget_Control, DebugText, $
         Set_UValue = ++NLines
   Widget_Control, DebugText, Set_Text_Top_Line = $
      NLines - DebugYSize + 3 > 0
Pro DebuggerHelper::O, Output, _Extra = Extra
   Compile_Opt IDL2
   self->DebugOutput, Output, /Up, _Extra = Extra
Pro DebuggerHelper::Off
   Compile_Opt IDL2
   self._DebugOn = 0
Pro DebuggerHelper::On
   Compile_Opt IDL2
   self._DebugOn = 1
Pro DebuggerHelper::Tic, Name
   Compile_Opt IDL2
   self._DebugClocks += Hash(Name, Tic(Name))
Pro DebuggerHelper::GetProperty, $
   On = On, $
   File = File
   Compile_Opt IDL2
   If (Arg_Present(On)) then On = self._DebugOn
   If (Arg_Present(File)) then $
      File = self._DebugFile
Function DebuggerHelper::Init, $
   On = On, $
   Group_Leader = Group_Leader, $
   No_File = No_File
   Compile_Opt IDL2
   self._DebugCreationTime = SysTime(1)
   self._DebugWindowMaxLines = 500
   self._DebugOn = Keyword_Set(On)
   self._DebugClocks = Hash()
   self->_ConstructDebugWidget, Group_Leader
   If (~Keyword_Set(No_File)) then $
   Return, 1
Pro DebuggerHelper__Define
   !null = {DebuggerHelper, Inherits IDL_Object, $
   _DebugOn : 0B, $
   _DebugLUN : 0L, $
   _DebugFile : '', $
   _DebugTLB : 0L, $
   _DebugCreationTime : 0D, $
   _DebugClocks : Hash(), $
   _DebugWindowMaxLines : 0L}

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

Categories: IDL Data Point





Sharing Ideas and Methods on Data

Author: Beau Legeer

When I am out of the office with customers or at tradeshows I am often asked about data, tools and algorithms. I get questions like, “what problems can I solve with this data” or “what are the best approaches towards solving this problem?”  It is natural to ask these questions since my colleagues and I at Exelis VIS work with any modality of remotely sensed data that provide information to solve a problem. We have seen all types formats that characterize data from multiple types of sensors.  

We want to lend our expertise to the community and a great place to start is the Esri International User Conference in July. We go every year, bring a booth and our experts. This year we are going to try something new – expert remote sensing consultations. During the conference we are going to accept requests for short 15-20 minute consultations on data, algorithms and other methods around solving problems with remotely sensed data.

We are excited to talk with you about the types of problems you are trying to solve and the methods you are planning to use. We are happy to provide feedback, suggestions and bring people who have worked on various forms of data to you. I am excited for this new idea and for the knowledge sharing that we can provide each other as a community of remote sensing experts. Fill out this form to request one of the 20 minute consultation spots we have available. We’ll contact you to let you know what time slot is yours. I look forward to seeing you at the UC.

Comments (1) Number of views (682) Article rating: No rating

Categories: IDL Data Point

Tags: imagery speaks

12345678910 Last




© 2014 Exelis Visual Information Solutions