Extended IDL Help

This page was created by the IDL library routine mk_html_help. For more information on this routine, refer to the IDL Online Help Navigator or type:

     ? mk_html_help

at the IDL command line prompt.

Last modified: Mon Sep 12 07:11:34 2005.


List of Routines


Routine Descriptions

AREA_SEARCH

[Next Routine] [List of Routines]
Find the point where A is maximum in the neighborhood
of (px,py). Search the region [px-dx:px+dx,py-dy:py+dy].

Return the location in A coordinates.

Harvey Rhody
January 25, 2005
For use in matching points between images.

(See C:\RSI\User\HarveyLib\area_search.pro)


ARRAYINTPRINT

[Previous Routine] [Next Routine] [List of Routines]
PROCEDURE ArrayIntPrint,A produces a Latex array
that is printed on the screen and can be cut and pasted
into a Latex document as a TeX field. Can be used wherever
a tabular field is desired. Use a text editor to add elements
such as \hline and so on if you need dividing lines. The
entries are printed as integers.

USAGE
ArrayPrint,A
Prints results on the screen.

KEYWORDS:
LINES  Set /Lines to have the array printed with row and column lines.
SWP_Table  Set this keyword to print out the table header and footer
           so the result can be a floating table in a SWP document
EQUATION   Set Equation to construct an equation to be pasted into
           a SWP document. Equation and SWP_TAble are mutually
           exclusive.

Version 1.0 Harvey Rhody  June 5, 2000

(See C:\RSI\User\HarveyLib\arrayintprint.pro)


ARRAYPRINT

[Previous Routine] [Next Routine] [List of Routines]
PROCEDURE Arrayprint,A produces a Latex array
that is printed on the screen and can be cut and pasted
into a Latex document as a TeX field. Can be used wherever
a tabular field is desired. Use a text editor to add elements
such as \hline and so on if you need dividing lines.

USAGE
ArrayPrint,A
Prints results on the screen.

KEYWORDS: None

Version 1.0 Harvey Rhody  June 5, 2000

(See C:\RSI\User\HarveyLib\arrayprint.pro)


ARRAYSTRPRINT

[Previous Routine] [Next Routine] [List of Routines]
PROCEDURE ArrayStrPrint,A produces a Latex array
from a string array. It inserts & characters between
intems and terminates each row with \\. The result
is printed on the screen and can be cut and pasted
into a Latex document as a TeX field. Can be used wherever
a tabular field is desired. Use a text editor to add elements
such as \hline and so on if you need dividing lines. The
entries are printed as integers.

USAGE
ArrayStrPrint,A
Prints results on the screen.

KEYWORDS:
LINES  Set /Lines to have the array printed with row and column lines.

Harvey Rhody  April 21, 2001

(See C:\RSI\User\HarveyLib\arraystrprint.pro)


CHOL

[Previous Routine] [Next Routine] [List of Routines]
L=chol(S) returns the lower triangular matrix such that
S=LL'. S must be positive definite.

(See C:\RSI\User\HarveyLib\chol.pro)


CLIP

[Previous Routine] [Next Routine] [List of Routines]
S=CLIP(A,T) returns S=A except small (absolute) values are zeroed.
The values are clipped at the threshold T. The default for
T is 1E-6. Hence, S=CLIP(A) zeros everything smaller than 1e-6.

H. Rhody, July, 2000

(See C:\RSI\User\HarveyLib\clip.pro)


CLOSEALL

[Previous Routine] [Next Routine] [List of Routines]
Procedure CLOSEALL closes all open windows.

 KEYWORDS
 Set the keyword /LIST to have a list of the indices
 of the windows printed to the screen.

(See C:\RSI\User\HarveyLib\Closeall.pro)


COORDINATE_GRID

[Previous Routine] [Next Routine] [List of Routines]
coordinate_grid,X,Y,N,M makes a grid of (x,y) coordinates of size NxM,
where N is the number of columns and M is the number of rows.

coordinate_grid,X,Y,N,M,/CENTER produces X and Y with (0,0) at position
(N/2,M/2).

H. Rhody
October 31, 2002

(See C:\RSI\User\HarveyLib\coordinate_grid.pro)


COS_ARRAY

[Previous Routine] [Next Routine] [List of Routines]
COS_ARRAY,N,M,A,W,U,V,R constructs NxM array A containing
the functions
fab(x,y)=cos(2*pi*(ax/N+by/M)) as rows for
a over the range 0,1,2...N-1 and
b over the range 0,1,2,...M-1. The SVD of A is then
calculated and returned in W,U,V. The rank R
is computed by the number of singular values
that are greater than 1e-5
See SVDC.

H. Rhody, July, 2000

(See C:\RSI\User\HarveyLib\cos_array.pro)


CUMSUM

[Previous Routine] [Next Routine] [List of Routines]
 s=CUMSUM(v) is a vector such that s[0]=v[0], s[1]=v[0]+v[1],..
 so that s is the cumulative sum over v. If v is an array
 with multiple dimensions, then s has the same shape as v.

 H. Rhody
 February 1999

(See C:\RSI\User\HarveyLib\cumsum.pro)


DCT

[Previous Routine] [Next Routine] [List of Routines]
Y=dct(X) is the discrete cosine transform of a vector
X of length N or of a square array X of size NxN.

Z=dct(Y,/INVERSE) computes the inverse

If you are going to do the same transform or the inverse
many times then it is inefficient to compute C each time.
You can save C the first time by using keyword T.
Y1=DCT(X1,T=C) constructs and returns the array C
while doing the forward transform. Now, with a new
array X2
Y2=DCT(X2,C) uses array C without computing it. To do
the inverse simply
Z2=DCT(Y2,C,/INVERSE) etc.

H. Rhody
April 20, 2005

(See C:\RSI\User\HarveyLib\dct.pro)


DIAG

[Previous Routine] [Next Routine] [List of Routines]
A=diag(v) constructs a diagonal matrix with the vector
v as the diagonal.

A=diag(v,n) constructs a diagonal matrix of size nxn
using the first n elements of v if nn_elements(v).

Note: Kept for old programs. It is better to use the built-in
DIAG_MATRIX program in IDL.

Harvey Rhody
May 18, 2002

(See C:\RSI\User\HarveyLib\diag.pro)


DIAGONAL

[Previous Routine] [Next Routine] [List of Routines]
r=diagonal(A,d) returns diagonal "d" of array A. If d GE 0
then the diagonal that starts at position index d along
the top row of the array is returned. If "d" is negative
then the diagonal that starts on the side of the array
at row |d| is returned.

NOTE: Kept to support old code. It is better to use the
DIAG_MATRIX function in IDL.

Harvey Rhody
January, 2003
Used as a midterm exam question for
SIMG-726.

(See C:\RSI\User\HarveyLib\diagonal.pro)


DIR_CHAR

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
    dir_char

 PURPOSE:
    Get the platform specific directory separation character

 CATEGORY:
   Directories, Files

 CALLING SEQUENCE:
  result=dir_char()
 
 INPUTS:
  none

 OPTIONAL INPUTS:
  none
      
 KEYWORD PARAMETERS:
  none
 
 OUTPUTS:
  returns the dir separation character
 
 OPTIONAL OUTPUTS:

 COMMON BLOCKS:

 SIDE EFFECTS:

 RESTRICTIONS:

 PROCEDURE:

 EXAMPLE:
   print, dir_char()

 MODIFICATION HISTORY:
   Dominik Brunner, ETHZ, Switzerland, first implementation

(See C:\RSI\User\HarveyLib\dir_char.pro)


DISC

[Previous Routine] [Next Routine] [List of Routines]
s=disc(n) is an nxn array that contains the distance of each point
from the center.

USAGE:
s=disc(7)
print,s
      4.24264      3.60555      3.16228      3.00000      3.16228      3.60555      4.24264
      3.60555      2.82843      2.23607      2.00000      2.23607      2.82843      3.60555
      3.16228      2.23607      1.41421      1.00000      1.41421      2.23607      3.16228
      3.00000      2.00000      1.00000     0.000000      1.00000      2.00000      3.00000
      3.16228      2.23607      1.41421      1.00000      1.41421      2.23607      3.16228
      3.60555      2.82843      2.23607      2.00000      2.23607      2.82843      3.60555
      4.24264      3.60555      3.16228      3.00000      3.16228      3.60555      4.24264

H. Rhody
Sept. 7, 2005

(See C:\RSI\User\HarveyLib\disc.pro)


DISP_IMAGE

[Previous Routine] [Next Routine] [List of Routines]
DISP_IMAGE,A opens a window of the correct size for an image
array A and displays the image. If A is a RGB image then the
display is in color. It does not set the palette or the decomposed
state.

DISP_IMAGE,A,wid displays A in window ID wid.

Keywords are passed through to the TV command.

(See C:\RSI\User\HarveyLib\disp_image.pro)


DISTANCE_POINT_TO_LINE

[Previous Routine] [Next Routine] [List of Routines]
d=distance_point_to_line(A,B,C) finds the
Euclidean distance from a line described by two
points A and B to a point C. All of the points
are n-dimensional vectors, where n>1. If A=B then
the distance between two points A and C is measured.

KEYWORD
D=distance_point_to_line(A,B,C,P=u)
returns the point P on the line closest to the point C.

H. Rhody
July, 2004

(See C:\RSI\User\HarveyLib\distance_point_to_line.pro)


DOW

[Previous Routine] [Next Routine] [List of Routines]
d=DOW(mon,day,year) or d=DOW(JulianDay) returns the day of the week
for the specified date. The input can be input either in the form of
mon,day,year or as the integer representing the Julian date. A vector
of dates can be given in the Julian integer form.

KEYWORDS:
NUM Set /NUM if you want the day to be returned as an integer
rather than as a string. 0='Sun', 1='Mon', etc.

EXAMPLES
print,DOW(7,4,1776) prints Thu
print,DOW(2452276) prints Tue
print,DOW(julday(12,30,1994)+Indgen(5)) prints Fri Sat Sun Mon Tue
print,DOW(julday(12,30,1994)+Indgen(5),/NUM) prints 5  6  0  1  2

SEE Also JULDAY, CALDAT, SYSTIME

HISTORY
Written by H. Rhody
January 3, 2002

(See C:\RSI\User\HarveyLib\dow.pro)


ENTROPY

[Previous Routine] [Next Routine] [List of Routines]
H=entropy(p) returns the entropy associated with the
probability array p. If two arrays, p and s, are provided then
the returned value is total(p*alog(s))/alog(2). The
calculation is done over the nonzero values of the
array in the logarithm.

(See C:\RSI\User\HarveyLib\entropy.pro)


EOL_CHAR

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
    eol_char

 PURPOSE:
    Get the platform specific end-of-line character
;
 CALLING SEQUENCE:
  result=eol_char()

 INPUTS:
  none

 OPTIONAL INPUTS:
  none

 KEYWORD PARAMETERS:
  none

 OUTPUTS:
  returns the eol character

 OPTIONAL OUTPUTS:

 COMMON BLOCKS:

 SIDE EFFECTS:

 RESTRICTIONS:

 PROCEDURE:

 EXAMPLE:
   print,'Fred is a dog.'+eol_char()+'Everybody likes Fred.'

 MODIFICATION HISTORY:
   Based on the program dir_char by
   Dominik Brunner, ETHZ, Switzerland, first implementation
   in his HipHop distribution.

   Modified to create eol_char() by Harvey Rhody, February 2002.

(See C:\RSI\User\HarveyLib\eol_char.pro)


FILTER

[Previous Routine] [Next Routine] [List of Routines]
 y=filter(b,a,x) where a and b are the ARMA filter functions
 and x is a number sequence. The filter output is generated by

 y(n)=b(0)x(n)+b(1)x(n-1)+...+b(Q)x(n-Q)-a(1)y(n-1)-a(2)y(n-2)-...
										  -a(P)y(n-P)
 This construction assumes that a(0)=1.

 INPUTS
 a=[1,a[1],...,a[P]] "autoregressive" coefficients
 b=[b[0]...b[Q]], "moving average" coefficients
 x=[x[0]...x[N-1]] input sequence of length N.

 The element a[0] is required but ignored. The value is set
 to a[0]=1 by the program. But don't omit it because then
 the program would, in effect, shorten the a vector.

 Pure Moving Average Filter: a=1 and b=[b[0]...b[Q]]
 Pure Autoregressive Filter: b=b0 and a=[1,a[1],...,a[P]]

 OUTPUTS
 The sequence y=[y[0]...y[N-1]].

 KEYWORD PARAMETERS
 XSTATE - set this parameter to a set of values that will
          serve as the past input values. This is a vector
          of length Q or less that represents the "past" values
          of the input.

 YSTATE - A vector of up to P past output values used to
          initialize the system.

 HISTORY
 Harvey E. Rhody  January, 1997, Initial version of program.
                  July, 2000, Added XSTATE and YSTATE keywords.

(See C:\RSI\User\HarveyLib\Filter.pro)


FILTER_SS

[Previous Routine] [Next Routine] [List of Routines]
y=filter_ss(b,a,u) uses a state-space model to compute the
output sequence of a filter that is described by the
rational function Y(z)=(B(z)/A(z))U(z).

See documentation to filter.

KEYWORDS
y=filter_ss(b,a,u,STATES=x) returns an array in which column
n is the state at step n. The first column of x is the
initial state. y(n)=H##x(n+1) since the state is updated
prior to computing the output.

X0=x0 is the initial state, which must be a column vector
of length N, where N is the order of polynomial a.

This implementation assumes N=n_elements(a) GE
M=n_elements(b).

H. Rhody
February 4, 2005
April 24, 2005 Revised to properly handle b longer than a

(See C:\RSI\User\HarveyLib\filter_ss.pro)


FIR_DESIGN

[Previous Routine] [Next Routine] [List of Routines]
h=FIR_DESIGN(S,M) computes the FIR coefficients for a filter of length M
that will approximate the frequency response S.

S is a vector of length N/2. The frequency response is constructed
to make S symmetrical about the origin.

H. Rhody
May, 2002

(See C:\RSI\User\HarveyLib\fir_design.pro)


FLIPUD

[Previous Routine] [Next Routine] [List of Routines]
Flip a 2D array A vertically.

(See C:\RSI\User\HarveyLib\flipUD.pro)


FSC_COLOR

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
       FSC_COLOR

 PURPOSE:

       The purpose of this function is to obtain drawing colors
       by name and in a device-decomposition independent way. The
       color names and values may be read in as a file, or 88
       color names and values are supplied from the program. These
       were obtained from the file rgb.txt, found on most X-Window
       distributions. Representative colors were chose from across
       the color spectrum. To see a list of colors available, type:
       Print, FSC_Color(/Names).

 AUTHOR:

       FANNING SOFTWARE CONSULTING:
       David Fanning, Ph.D.
       1645 Sheely Drive
       Fort Collins, CO 80526 USA
       Phone: 970-221-0438
       E-mail: davidf@dfanning.com
       Coyote's Guide to IDL Programming: http://www.dfanning.com

 CATEGORY:

       Graphics, Color Specification.

 CALLING SEQUENCE:

       color = FSC_COLOR(theColor, theColorIndex)

 NORMAL CALLING SEQUENCE FOR DEVICE-INDEPENDENT COLOR:

       axisColor = FSC_COLOR("Green", !D.Table_Size-2)
       backColor = FSC_COLOR("Charcoal", !D.Table_Size-3)
       dataColor = FSC_COLOR("Yellow", !D.Table_Size-4)
       Plot, Findgen(11), Color=axisColor, Background=backColor, /NoData
       OPlot, Findgen(11), Color=dataColor

 OPTIONAL INPUT PARAMETERS:

       TheColor: A string with the "name" of the color. To see a list
           of the color names available set the NAMES keyword. Colors available
           are these:

                Almond   Antique White      Aquamarine           Beige          Bisque           Black
                  Blue     Blue Violet           Brown       Burlywood        Charcoal      Chartreuse
             Chocolate           Coral        Cornsilk            Cyan  Dark Goldenrod       Dark Gray
            Dark Green      Dark Khaki     Dark Orchid     Dark Salmon       Deep Pink     Dodger Blue
             Firebrick    Forest Green            Gold       Goldenrod            Gray           Green
          Green Yellow        Honeydew        Hot Pink      Indian Red           Ivory           Khaki
              Lavender      Lawn Green     Light Coral      Light Cyan      Light Gray    Light Salmon
          Light Yellow      Lime Green           Linen         Magenta          Maroon     Medium Gray
         Medium Orchid        Moccasin            Navy           Olive      Olive Drab          Orange
            Orange Red          Orchid  Pale Goldenrod      Pale Green          Papaya            Peru
                  Pink            Plum     Powder Blue          Purple             Red            Rose
            Rosy Brown      Royal Blue    Saddle Brown          Salmon     Sandy Brown       Sea Green
              Seashell          Sienna        Sky Blue      Slate Gray            Snow    Spring Green
            Steel Blue             Tan         Thistle          Tomato       Turquoise          Violet
            Violet Red           Wheat           White          Yellow

           The color WHITE is used if this parameter is absent. To see a list
           of the color names available in the program, type this:

              Print, FSC_COLOR(/Names)

       TheColorIndex: The color table index where the specified color is loaded.
           The color table index parameter should always be used if you wish to
           obtain a color value in a color-decomposition-independent way in your
           code. See the NORMAL CALLING SEQUENCE for details.

 RETURN VALUE:

       The value that is returned by FSC_COLOR depends upon the keywords
       used to call it and on the version of IDL you are using. In general,
       the return value will be either a color index number where the specified
       color is loaded by the program, or a 24-bit color value that can be
       decomposed into the specified color on true-color systems.

       If you are running IDL 5.2 or higher, the program will determine which
       return value to use, based on the color decomposition state at the time
       the program is called. If you are running a version of IDL before IDL 5.2,
       then the program will return the color index number. This behavior can
       be overruled in all versions of IDL by setting the DECOMPOSED keyword.
       If this keyword is 0, the program always returns a color index number. If
       the keyword is 1, the program always returns a 24-bit color value.

       If the TRIPLE keyword is set, the program always returns the color triple,
       no matter what the current decomposition state or the value of the DECOMPOSED
       keyword. Normally, the color triple is returned as a 1 by 3 column vector.
       This is appropriate for loading into a color index with TVLCT:

          IDL> TVLCT, FSC_Color('Yellow', /Triple), !P.Color

       But sometimes (e.g, in object graphics applications) you want the color
       returned as a row vector. In this case, you should set the ROW keyword
       as well as the TRIPLE keyword:

          viewobj= Obj_New('IDLgrView', Color=FSC_Color('charcoal', /Triple, /Row))

       If the ALLCOLORS keyword is used, then instead of a single value, modified
       as described above, then all the color values are returned in an array. In
       other words, the return value will be either an NCOLORS-element vector of color
       table index numbers, an NCOLORS-element vector of 24-bit color values, or
       an NCOLORS-by-3 array of color triples.

       If the NAMES keyword is set, the program returns a vector of
       color names known to the program.

 INPUT KEYWORD PARAMETERS:

       ALLCOLORS: Set this keyword to return indices, or 24-bit values, or color
              triples, for all the known colors, instead of for a single color.

       DECOMPOSED: Set this keyword to 0 or 1 to force the return value to be
              a color table index or a 24-bit color value, respectively.

       FILENAME: The string name of an ASCII file that can be opened to read in
              color values and color names. There should be one color per row
              in the file. Please be sure there are no blank lines in the file.
              The format of each row should be:

                  redValue  greenValue  blueValue  colorName

              Color values should be between 0 and 255. Any kind of white-space
              separation (blank characters, commas, or tabs) are allowed. The color
              name should be a string, but it should NOT be in quotes. A typical
              entry into the file would look like this:

                  255   255   0   Yellow

       NAMES: If this keyword is set, the return value of the function is
              a ncolors-element string array containing the names of the colors.
              These names would be appropriate, for example, in building
              a list widget with the names of the colors. If the NAMES
              keyword is set, the COLOR and INDEX parameters are ignored.

                 listID = Widget_List(baseID, Value=GetColor(/Names), YSize=16)

       ROW:   If this keyword is set, the return value of the function when the TRIPLE
              keyword is set is returned as a row vector, rather than as the default
              column vector. This is required, for example, when you are trying to
              use the return value to set the color for object graphics objects. This
              keyword is completely ignored, except when used in combination with the
              TRIPLE keyword.

       SELECTCOLOR: Set this keyword if you would like to select the color name with
              the PICKCOLORNAME program. Selecting this keyword automaticallys sets
              the INDEX positional parameter. If this keyword is used, any keywords
              appropriate for PICKCOLORNAME can also be used. If this keyword is used,
              the first positional parameter can be either a color name or the color
              table index number. The program will figure out what you want.

       TRIPLE: Setting this keyword will force the return value of the function to
              *always* be a color triple, regardless of color decomposition state or
              visual depth of the machine. The value will be a three-element column
              vector unless the ROW keyword is also set.

       In addition, any keyword parameter appropriate for PICKCOLORNAME can be used.
       These include BOTTOM, COLUMNS, GROUP_LEADER, INDEX, and TITLE.

 OUTPUT KEYWORD PARAMETERS:

       CANCEL: This keyword is always set to 0, unless that SELECTCOLOR keyword is used.
              Then it will correspond to the value of the CANCEL output keyword in PICKCOLORNAME.

       COLORSTRUCTURE: This output keyword (if set to a named variable) will return a
              structure in which the fields will be the known color names (without spaces)
              and the values of the fields will be either color table index numbers or
              24-bit color values.

       NCOLORS: The number of colors recognized by the program. It will be 88 by default.

 COMMON BLOCKS:
       None.

 SIDE EFFECTS:
       None.

 ADDITIONAL PROGRAMS REQUIRED:

   PICKCOLORNAME: This file can be found in the Coyote Library:

             http://www.dfanning.com/programs/pickcolorname.pro

 EXAMPLE:

       To get drawing colors in a device-decomposed independent way:

           axisColor = FSC_COLOR("Green", !D.Table_Size-2)
           backColor = FSC_COLOR("Charcoal", !D.Table_Size-3)
           dataColor = FSC_COLOR("Yellow", !D.Table_Size-4)
           Plot, Findgen(11), Color=axisColor, Background=backColor, /NoData
           OPlot, Findgen(11), Color=dataColor

       To set the viewport color in object graphics:

           theView = Obj_New('IDLgrView', Color=FSC_Color('Charcoal', /Triple))

       To change the viewport color later:

           theView->SetProperty, Color=FSC_Color('Antique White', /Triple)

 MODIFICATION HISTORY:
       Written by: David Fanning, 19 October 2000. Based on previous
          GetColor program.
       Fixed a problem with loading colors with TVLCT on a PRINTER device. 13 Mar 2001. DWF.
       Added the ROW keyword. 30 March 2001. DWF.
       Added the PICKCOLORNAME code to the file, since I keep forgetting to
          give it to people. 15 August 2001. DWF.

(See C:\RSI\User\HarveyLib\fsc_color.pro)


GAUS

[Previous Routine] [Next Routine] [List of Routines]
 FUNCTION GAUS,X
 Computes the standard GAUS function for any list of values represented
 by the vector x. See Gaskill, Page 47.

 USAGE:
 x=FINDGEN(1001)/20-25 CREATE A VECTOR ON THE INTERVAL [-25,25]
 y=GAUS((x-8)/5) Construct a function that is shifted and stretched

 HISTORY
 Written by Harvey Rhody, September 1997

(See C:\RSI\User\HarveyLib\GAUS.PRO)


GRAMSCHMIDT

[Previous Routine] [Next Routine] [List of Routines]
Q=gramSchmidt(A,R) returns Q containing the Gram-Schmidt
orthogonalization of the columns of A. An upper triangular
matrix R such that A=Q##R is returned in the second argument.

H. Rhody July, 2000

(See C:\RSI\User\HarveyLib\gramschmidt.pro)


HAAR_MATRIX

[Previous Routine] [Next Routine] [List of Routines]
Compute the Haar transform matrix H of size NxN
using Gonzalez and Woods equation 7.1-26, page 361.

N must be a power of 2.

The Haar transform of an NxN array F is computed by
T=H##F##H.

(See C:\RSI\User\HarveyLib\haar_matrix.pro)


HDR_INFO

[Previous Routine] [Next Routine] [List of Routines]
info=hdr_info(fname) returns a structure with a field and
value for each line of the file that contains an = sign.

EXAMPLE
Suppose that file 'fred.hdr' contains the lines

ENVI
DESCRIPTION = {SWIR image}
SAMPLES = 640
LINES = 510

Then info=hdr_info('fred.hdr') would produce  a structure
info with info.DESCRIPTION='{SWIR image}', info.SAMPLES='640'
and info.LINES=510.

The program only pays attention to lines that contain = .
Everything after = is the value.

H. Rhody
June 28, 2004

(See C:\RSI\User\HarveyLib\hdr_info.pro)


HR_LDP

[Previous Routine] [Next Routine] [List of Routines]
algorithm 23.27 to minimize ||z|| subject
to Gz>=h. If the inequality constraints
are consistent then phi=TRUE and a value will
be returned for vector z. The value of
phi=total(abs(z)) is returned. A very small
value of phi indicates a possible inconsistent
constraint set or a failure of the algorithm to
converge to a good answer.

H. Rhody
May, 2004

(See C:\RSI\User\HarveyLib\HR_LDP.pro)


HUFFMAN

[Previous Routine] [Next Routine] [List of Routines]
 C=Huffman(problist,symlist) makes a binary Huffman code

 symlist is a list of character symbols that corresponds to
 the list of probabilities.

 c=huffman(problist,symlist) returns an array of code strings
 associated with the symbols as determined by the Huffman coding
 algorithm. The code is returned in an array of structures sorted
 in order of descending probability. The structure structue tags
 are "codeword", "probability", "length" and "symbol". If the symlist is
 omitted then the symbol field is blank in each case.


KEYWORDS

 ENTROPY
  c=huffman(problist,symlist,Entropy=entropy) causes a value to be
  returned for entropy of the code.

 AVGLEN
  c=huffman(problist,symlist,Entropy=entropy,AvgLen=avglen) causes a value to be
  returned for entropy and the average length of the code.

 VERBOSE
  Setting /VERBOSE causes the code progress to be printed out and the entropy and
  average codeword length to be reported at the end.

 LATEX
  Setting /LATEX causes a printout that can be cut and pasted into a LaTeX report.

 This version will only produce a binary code.


 EXAMPLE:
 p=[.2,.3,.1,.4]
 s=['a','b','c','d']
 c=huffman(p,s)
 for n=0,3 DO print,format='(a2,2x,f3.1,1x,a4)',c[n].symbol,$
	c[n].probability,c[n].codeword

 produces the printout shown below:

 d  0.4    0
 b  0.3   10
 a  0.2  111
 c  0.1  110

 HISTORY:
 Written by H. Rhody February, 1999
 Modified March 2001 to make symlist a positional
   parameter and ENTROPY  and AVGLEN keyword parameters
   and to include keyword VERBOSE.
 Modified May 2001 to add the length tag.

(See C:\RSI\User\HarveyLib\huffman.pro)


HWT

[Previous Routine] [Next Routine] [List of Routines]
B=HWT(A) returns the Haar wavelet transform of array A.
A may be a row vector or a 2D array. The number of
columns must be even (or the result will be truncated
to an even number). The result B is a
floating point array with the same dimensions as A
provided A has an even number of columns.

If A has dimension NxM then B[0:N/2-1,*] contains the
low-frequency part of the result and B[N/2:2*(N/2)-1,*] contains
the high-frequency part.

The inverse transform is computed by A=HWT(B,/INVERSE).
The result is a floating point array. If it is important
to maintain the array type then that must be done
as an external step.

EXAMPLE 1:
A=[1,-3,0,4]
B=HWT(A)
PRINT,B
     -1.41421      2.82843      2.82843     -2.82843
C=HWT(B,/INVERSE)
PRINT,C
      1.00000     -3.00000     0.000000      4.00000

EXAMPLE 2:
A=DIST(4) & PRINT,A

      0.000000      1.00000      2.00000      1.00000
      1.00000      1.41421      2.23607      1.41421
      2.00000      2.23607      2.82843      2.23607
      1.00000      1.41421      2.23607      1.41421

Compute the HWT along the rows.
B1=HWT(A) & PRINT,B1

      0.707107      2.12132    -0.707107     0.707107
      1.70711      2.58114    -0.292893     0.581139
      2.99535      3.58114    -0.166925     0.418861
      1.70711      2.58114    -0.292893     0.581139

Transpose B1 and compute the HWT, which does the columns.
B2=HWT(TRANSPOSE(B1)) & PRINT,B2

      1.70711      3.32514    -0.707107     0.910927
      3.32514      4.35739    -0.325141     0.707107
    -0.707107    -0.325141    -0.292893    0.0890728
     0.910927     0.707107    0.0890728    -0.114748

Compute the inverse HWT of B2
C1=HWT(B2,/INVERSE) & PRINT,C1

      0.707107      1.70711      2.99535      1.70711
      2.12132      2.58114      3.58114      2.58114
    -0.707107    -0.292893    -0.166925    -0.292893
     0.707107     0.581139     0.418861     0.581139

C1 is seen to be the same as TRANSPOSE(B1). We can
recover the original array by
C=HWT(TRANSPOSE(C1),/INVERSE)

      0.000000      1.00000      2.00000      1.00000
      1.00000      1.41421      2.23607      1.41421
      2.00000      2.23607      2.82843      2.23607
      1.00000      1.41421      2.23607      1.41421

To do the HWT on both rows and columns it is necessary that
the number of rows also be an even number.

;EXAMPLE 3: Odd dimensions.
A=DIST(5) & PRINT,A
      0.000000      1.00000      2.00000      2.00000      1.00000
      1.00000      1.41421      2.23607      2.23607      1.41421
      2.00000      2.23607      2.82843      2.82843      2.23607
      2.00000      2.23607      2.82843      2.82843      2.23607
      1.00000      1.41421      2.23607      2.23607      1.41421

B1=HWT(A) & PRINT,B1

      0.707107      2.82843    -0.707107     0.000000
      1.70711      3.16228    -0.292893     0.000000
      2.99535      4.00000    -0.166925     0.000000
      2.99535      4.00000    -0.166925     0.000000
      1.70711      3.16228    -0.292893     0.000000

B2=HWT(TRANSPOSE(B1) & PRINT,B2

      1.70711      4.23607    -0.707107     0.000000
      4.23607      5.65685    -0.236068     0.000000
    -0.707107    -0.236068    -0.292893     0.000000
     0.000000     0.000000     0.000000     0.000000

We now have an even number of rows and columns. Now
see what we get when we transform back.

B3=HWT(B2,/INVERSE) & PRINT,B3

     0.707107      1.70711      2.99535      2.99535
      2.82843      3.16228      4.00000      4.00000
    -0.707107    -0.292893    -0.166925    -0.166925
     0.000000     0.000000     0.000000     0.000000

This is the transpose of B1 with one less column. Now
invert this.

B4=HWT(TRANSPOSE(B3),/INVERSE) & PRINT,B4

     0.000000      1.00000      2.00000      2.00000
      1.00000      1.41421      2.23607      2.23607
      2.00000      2.23607      2.82843      2.82843
      2.00000      2.23607      2.82843      2.82843

The result is the same as the original array except for the
truncation of the number of rows and columns to an even
number.


H. Rhody
November, 2002

(See C:\RSI\User\HarveyLib\hwt.pro)


INTERSECTION

[Previous Routine] [Next Routine] [List of Routines]
C=INTERSECTION(A,B) is the intersection of the sets A and B. It is assumed
that A and B have some form of integer values. That is, they
come from the set of non-negative integers.

(See C:\RSI\User\HarveyLib\intersection.pro)


KRONECKER

[Previous Routine] [Next Routine] [List of Routines]
C=kronecker(A,B) is the Kronecker product of A with B. Array
C has size na*nb X ma*mb, where A is of size na X nb and B is
of size nb X mb. C is formed by "tiling". Tile [i,j] of C is
a[i,j]B.

(See C:\RSI\User\HarveyLib\kronecker.pro)


LINE_NORMAL_FORM

[Previous Routine] [Next Routine] [List of Routines]
line_normal_form,px,py,rho,theta returns the normal
rho,theta representation of the regression line
through the points represented by px and py.

rho is in the same units as px and py and is
always non-negative.

theta is in radians and ranges from (-pi,pi)

H. Rhody
July, 2004

(See C:\RSI\User\HarveyLib\line_normal_form.pro)


LS

[Previous Routine] [Next Routine] [List of Routines]
ls lists the files in the current directory.
ls,filter causes files containing the specified filter string
       to be selected.
       Example: ls,'*.tif' selects tiff files.
       Any selection string accepted by FINDFILE can be used.

KEYWORDS
FILES=files returns the list in the string array "files".

DIRECTORY=dir causes files in directory "dir" to be listed.

COUNT=count returns the number of files found.

NOPRINT: /NOPRINT to suppress printing. (Default is to print.)

February 2002: Added NOPRINT keyword.

(See C:\RSI\User\HarveyLib\ls.pro)


LSIE

[Previous Routine] [Next Routine] [List of Routines]
Minimize ||Ex-f|| subject to Cx=d and Gx>=h

Use ortho_factor to find the C=HRK^T

The parameter phi is returned by LDP, which is the
search routine. If phi is very small, the solution may
be faulty.

H. Rhody
June, 2004

(See C:\RSI\User\HarveyLib\LSIE.PRO)


MAKEMYHTML

[Previous Routine] [Next Routine] [List of Routines]
makeMyHtml constructs HarveyHelp.html from the
imglib and HarveyLib directories and puts the
resulting HTML file on the desktop.

USAGE: makeMyHTML

(See C:\RSI\User\HarveyLib\makemyhtml.pro)


MAKE_BASIS_IMAGES

[Previous Routine] [Next Routine] [List of Routines]
B=MAKE_BASIS_IMAGES(A,n,m) constructs the (n,m) basis image
formed from the unitary matrix A. The image is the
outer product of columns n and m of the Hermetian
of A.

A unitary matrix for some common transforms may be
constructed with the function TRANSFORM_BASIS.

REFERENCE:
A. K. JAIN, Funcamentals of Digital Image Processing, p. 135.

H. Rhody, July, 2000

(See C:\RSI\User\HarveyLib\make_basis_images.pro)


MATCH_POINT

[Previous Routine] [Next Routine] [List of Routines]
R=SAMPLE_CLICK(win1,win2) uses the mouse cursor to collect matching
points from windows win1 and win2. The windows need not be the same size.

To use the process, click the left mouse in the active window. This
causes the coordinates of the point to be determined. The other window
is then made active. Clicking there collects the matching point.

To terminate the process, right-click in the current window. If the first
window is active the process returns the previously-collected pairs of
points. If the second window is active, then that point is collected
and the last pair of points is included in the return.

R is an anonymous structure with arrays P=[X1,Y1] and Q=[X2,Y2]. Display
the points with R.P and R.Q.


HISTORY:
Written by H. Rhody, April 8, 2004

(See C:\RSI\User\HarveyLib\match_point.pro)


MATRIXINTPRINT

[Previous Routine] [Next Routine] [List of Routines]
PROCEDURE MatrixIntPrint,A produces a Latex matrix
that is printed on the screen and can be cut and pasted
into a Latex document as a TeX field. Can be used wherever
a math array is desired. Use the ArrayPrint procedure to
get a tabular array.

Produces output in integer format. Use MatrixPrint for a
floating point output

USAGE
ArrayPrint,A
Prints results on the screen.

KEYWORDS: None

Version 1.0 Harvey Rhody  June 5, 2000

(See C:\RSI\User\HarveyLib\MatrixIntPrint.pro)


MATRIXPRINT

[Previous Routine] [Next Routine] [List of Routines]
PROCEDURE MatrixPrint,A produces a Latex matrix
that is printed on the screen and can be cut and pasted
into a Latex document as a TeX field. Can be used wherever
a math array is desired. Use the ArrayPrint procedure to
get a tabular array.

USAGE
ArrayPrint,A
Prints results on the screen.

KEYWORDS: None

Version 1.0 Harvey Rhody  June 5, 2000

(See C:\RSI\User\HarveyLib\matrixprint.pro)


MESHDOM

[Previous Routine] [Next Routine] [List of Routines]
 MESHDOM,x,y,XR,YR

 MESHDOM calculates arrays XR and YR that index a region with
 value pairs from (x,y). The arrays are of size Nx columns by
 Ny rows, where Nx is the length of x and Ny is the length of y.
 Each row of XR equals the vector x, and each column of YR
 equals the vector y.

 USAGE
 x=indgen(10) + 5
 y=indgen(5) - 3
 MESHDOM,x,y,XR,YR

 The results are the arrays shown below.

 ARRAY XR
 -5   -4   -3   -2   -1    0    1    2    3    4
 -5   -4   -3   -2   -1    0    1    2    3    4
 -5   -4   -3   -2   -1    0    1    2    3    4
 -5   -4   -3   -2   -1    0    1    2    3    4
 -5   -4   -3   -2   -1    0    1    2    3    4

 ARRAY YR
 -3   -3   -3   -3   -3   -3   -3   -3   -3   -3
 -2   -2   -2   -2   -2   -2   -2   -2   -2   -2
 -1   -1   -1   -1   -1   -1   -1   -1   -1   -1
  0    0    0    0    0    0    0    0    0    0
  1    1    1    1    1    1    1    1    1    1


 HISTORY
 Written by Harvey Rhody, November, 1996
 Updated September, 1997

NAME
MESHDOM

(See C:\RSI\User\HarveyLib\Meshdom.pro)


MESHGRID

[Previous Routine] [Next Routine] [List of Routines]
 PRO MESHGRID,N,M,X,Y calculates two arrays, X
 and Y that can be used to calculate and plot a 2D function.

 N and M can be either positive integers or vectors. If they are
 vectors then N is used as the rows of X and M is used as the columns of
 Y. If they are integers then the rows of X are IndGen(N) and the columns
 of Y are Indgen(M).

 Example 1
 MESHGRID,31,31,X,Y
 X=(X-15)/2.
 Y=(Y-11)/3.
 Z=EXP(-(X^2+Y^2))
 SURFACE,Z

 Example 2
 MESHGRID,[2,4,6],[3,5,7,9],X,Y
 creates two arrays of size 3x4 with the X array containg four rows
 with [2,4,6] and Y with columns [3,5,7,9]'.

 HISTORY:
 Originated by Harvey Rhody September 17, 1997.
 Revised April 2001 to accomodate vectors for N and M.

(See C:\RSI\User\HarveyLib\meshgrid.pro)


MSEMATCH

[Previous Routine] [Next Routine] [List of Routines]
p=MSEmatch(x,y) provides a vector p=[a,b] such that
mean-squared difference between x and (ay+b) is minimized.

(See C:\RSI\User\HarveyLib\MSEmatch.pro)


NORM

[Previous Routine] [Next Routine] [List of Routines]
y=norm(x) returns the Euclidean norm of a vector x

H. Rhody
July, 2000

(See C:\RSI\User\HarveyLib\norm.pro)


NORMALFORM

[Previous Routine] [Next Routine] [List of Routines]
n=NormalForm(X,Y) returns a vector n=[rho,theta] that
represents the normal form of a line that is devined
by the points represented by X and Y. The vector X
contains the x-coordinates and Y contains the y-coordinates
of the points.

HISTORY
H. Rhody, November 16, 2001--version 1.0

(See C:\RSI\User\HarveyLib\normalForm.pro)


ONES

[Previous Routine] [Next Routine] [List of Routines]
x=ones(n) constructs a row vector containing n ones.
x=ones(n,/COL) generates a column vector of n ones.

The default array type is INTEGER. However, the type
can be set to any of the types accepted by MAKE_ARRAY
by using the keyword TYPE. See Help on SIZE for a list
of type codes.

EXAMPLE
V=ones(10,/COL,TYPE=4) constructs a column of 10 floats.

H. Rhody July, 2000

(See C:\RSI\User\HarveyLib\ones.pro)


ORTHO_FACTOR

[Previous Routine] [Next Routine] [List of Routines]
rankA=Ortho_Factor(A,H,R,K,TAU=tau) returns the rank
of a matrix A and provides orthogonal matrices
H and K and an upper-triangular matrix R such that
A=H##R##Transpose(K). H,R and K are computed by the
QR algorithm.

The SVDC decomposition is available through arguments
s,U,V by calling
rankA=Ortho_Factor(A,H,R,K,s,U,V)

KEYWORD
TAU is the threshold for rank determination. Rank is
calculated by examining the singular values. If
s[k] is the last first singular value such that
s[k] LT s[0]*tau then the array has rank=k. The default
value is TAU=1e-6.

If A is a scalar then rankA=0. The SVD arrays are not
defined for this case.

If A is a row vector then rankA=1. The SVD arrays are
not defined for this case.

Returns -1 in cases not handled, such as arrays with
dimensions greater than 2D.

USES algorithm QR and IDL procedure SVDC.

H. Rhody
June, 2004

(See C:\RSI\User\HarveyLib\ortho_factor.pro)


ORTHO_RANK

[Previous Routine] [Next Routine] [List of Routines]
rankA=Ortho_Rank(A,TAU=tau) is the rank of a matrix A.
If A is a scalar then rankA=0. tau is the threshold
for zero values. Default: tau=1e-6.

Uses QR,A,Q,R and counts the number rows of R
that contain at least one value GE tau.

Arrays H,R,K such that A=H##R##Transpose(K) are
returned through the arguments if ortho_rank is called
with rankA=Ortho_Rank(A,TAU=tau,H,R,K)

H. Rhody
June, 2004

(See C:\RSI\User\HarveyLib\ortho_rank.pro)


POINTS_NEAR_LINE

[Previous Routine] [Next Routine] [List of Routines]
R=points_near_line(x0,y0,x1,y1,P,delta) or
R=points_near_line(V,pts,delta) returns the
points in the array P that are near a line described
by the endpoints [x0,y0], [x1,y1]. The distance threshold
is set by delta. The endpoints may be supplied as
the four-parameters x0,y0,x1,y1 or a vector V=[x0,y0,x1,y1].

P must be a [2,n] array of n points. R is the index list
of the points in P that are near the line. The points
may be selected by Q=P[*,R]

H. Rhody
July, 2004

(See C:\RSI\User\HarveyLib\points_near_line.pro)


POINT_SELECT

[Previous Routine] [Next Routine] [List of Routines]
R=POINT_SELECT(P,Q) uses the mouse cursor to alternately select values
from arrays P and Q that have been displayed in windows winP and winQ.

Left click in the image to collect and print the coordinates of a point.
Right-click to end the procedure. It is assumed that the windows are
sized to fit the image arrays.

Returned value R is a structure with six keys: Pv,Px,Py,Qv,Qx,Qy where
Pv and Qv are the sample values at the selected coordinates and Px,Py,Qx
Qy are the coordinates.

HISTORY:
Written by H. Rhody, January 25, 2005

Constructed by modification of sample_click.

(See C:\RSI\User\HarveyLib\point_select.pro)


POISSONPROCESS

[Previous Routine] [Next Routine] [List of Routines]
 S=PoissonProcess(t,lambda,p)
 divides the interval [0,t] into intervals of size
 deltaT=p/lambda where p is sufficiently small so that
 the Poisson assumptions are satisfied. The value of
 p may be specified by the third positional variable.
 If it is not provided, it is set to the default value
 p=0.1.

 The interval (0,t) is divided into n=t*lambda/p intervals
 and the number of events in the interval (0,k*deltaT) is
 returned in the array S. The maximum length of S is 10000.

 USAGE
 S=PoissonProcess(10,1,0.1)
 Plot,S
 FOR m=1,10 DO OPLOT,PoissonProcess(10,1,0.1)

 The above generate a number of overlayed plots of the process.

 Harvey Rhody
 March, 2000

(See C:\RSI\User\HarveyLib\PoissonProcess.pro)


POLY_INSIDE

[Previous Routine] [Next Routine] [List of Routines]
r=poly_inside(X,Y,P) determines whether a point P is
within the boundary of a polygon with vertices described
by the coordinate vectors X and Y, or on the boundary, or
outside the boundary.

H. Rhody
October 5, 2004

(See C:\RSI\User\HarveyLib\poly_inside.pro)


PSEUDO_INVERSE

[Previous Routine] [Next Routine] [List of Routines]
B=pseudo_inverse(A,[U,V,S,tol,EXTRA=extra)
computes the pseudo-inverse using SVDC.
B=pseudo_inverse(A,U,V,S) returns the SVDC values
U,V and S such that A=U##S##Transpose(V). This is
especially useful in examining the sensitivity of
the solution. Small values along the diagonal of
S can be zeroed out in computing the pseudo-inverse
to reduce numerical sensitivity.

PARAMETERS
Set TOL to a threshold value on the singular values
to be accepted in computing the inverse. Example
Ap=pseudo_inverse(A,TOL=1e-5) causes singular values
smaller than 1e-5 to be excluded from the inverse
calculation. Default: tol=1e-10

Set SORT to sort the singular values in descending
order before computing the inverse.

2004-05-13 HR Added SORT and TOL keywords and replaced
use of the special DIAG() function with DIAG_MATRIX()

Reference:
Moon and Stirling,
Mathematical Methods and Algorithms
Prentice-Hall, 2000
Section 7.4

H. Rhody  December 30, 2003

(See C:\RSI\User\HarveyLib\pseudo_inverse.pro)


QR

[Previous Routine] [Next Routine] [List of Routines]
QR,A,Q,R computes matrices Q and R such that A=QR where R is an
upper triangular matrix. Q is constructed by a sequence of
Householder transformations.

Harvey Rhody
Laboratory for Imaging Algorithms and Systems
Rochester Institute of Technology
March, 2004

(See C:\RSI\User\HarveyLib\qr.pro)


QRFAC

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
   QRFAC

 AUTHOR:
   Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
   craigm@lheamail.gsfc.nasa.gov

 PURPOSE:
   Perform QR decomposition of a rectangular matrix

 MAJOR TOPICS:
   Linear Systems

 CALLING SEQUENCE:
   QRFAC, A, R, [ IPVT, /PIVOT ]

 DESCRIPTION:

  Given an MxN matrix A (M>N), the procedure QRFAC computes the QR
  decomposition (factorization) of A.  This factorization is useful
  in least squares applications solving the equation, A ## x = B.
  Together with the procedure QRSOLV, this equation can be solved in
  a least squares sense.

  The QR factorization produces two matrices, Q and R, such that

     A = Q ## R

  where Q is orthogonal such that TRANSPOSE(Q)##Q equals the identity
  matrix, and R is upper triangular.  This procedure does not compute
  Q directly, but returns the more-compact Householder reflectors,
  which QRSOLV applies in constructing the solution.

  Pivoting can be performed by setting the PIVOT keyword.  Rows with
  the largest L2-norm are pivoted into the top positions of the
  matrix.  The permutation matrix is returned in the IPVT parameter.


 PARAMETERS:

   A - upon input, an MxN matrix to be factored.

       Upon output, the upper triangular MxN matrix of Householder
       reflectors used in reconstructing Q.  Obviously the original
       matrix A is destroyed upon output.

       Note that the dimensions of A in this routine are the
       *TRANSPOSE* of the conventional appearance in the least
       squares matrix equation.

   R - upon ouptut, an upper triangular NxN matrix

   IPVT - upon output, the permutation indices used in partial
          pivoting.  If pivoting is used, this array should be passed
          to the PIVOTS keyword of QRSOLV.  If the PIVOT keyword is
          not set, then IPVT returns an unpermuted array of indices.

 KEYWORD PARAMETERS:

   PIVOT - if set, then partial pivoting is performed, to bring the
           rows with the largest norm to the top of the matrix.

   QMATRIX - upon return, the fully explicit "Q" matrix is returned.
             This square matrix can be used to perform explicit
             matrix multiplication (although not super efficiently).
             The values returned modified in A are the Householder
             vectors, which are then used to compute QMAT.


 EXAMPLE:

  Decompose the 3x2 matrix [[9.,2.,6.],[4.,8.,7.]]
    aa = [[9.,2.,6.],[4.,8.,7.]]
    qrfac, aa, r, ipvt

     IDL> print, aa
          1.81818      0.181818      0.545455 
         XXXXXXXXX      1.90160      0.432573 
    (position marked with Xs is undefined)

  Construct the matrix Q by expanding the Householder reflectors
  returned in AA.  ( M = 3, N = 2 )  This same procedure is
  accomplished by using the QMATRIX keyword.

    ident = fltarr(m,m)  ;; Construct an identity matrix
    ident(lindgen(m),lindgen(m)) = 1

    q = ident
    for i = 0, n-1 do begin
     v = aa(*,i) & if i GT 0 then v(0:i-1) = 0  ;; extract reflector
     q = q ## (ident - 2*(v # v)/total(v * v))  ;; generate matrix
    endfor

  Verify that Q ## R returns to the original AA

     print, q(0:1,*) ## r
         9.00000      4.00000
         2.00000      8.00000
         6.00000      7.00000
     (transposed)

  See example in QRSOLV to solve a least squares problem.
   

 REFERENCES:

   More', Jorge J., "The Levenberg-Marquardt Algorithm:
     Implementation and Theory," in *Numerical Analysis*, ed. Watson,
     G. A., Lecture Notes in Mathematics 630, Springer-Verlag, 1977.

 MODIFICATION HISTORY:
   Written (taken from MPFIT), CM, Feb 2002
   Added usage message, error checking, CM 15 Mar 2002
   Corrected error in EXAMPLE, CM, 10 May 2002
   Now returns Q matrix explicitly if requested, CM, 14 Jul 2002
   Documented QMATRIX keyword, CM, 22 Jul 2002

  $Id: qrfac.pro,v 1.5 2003/12/18 09:28:43 craigm Exp $

(See C:\RSI\User\HarveyLib\qrfac.pro)


QRSOLV

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
   QRSOLV

 AUTHOR:
   Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
   craigm@lheamail.gsfc.nasa.gov

 PURPOSE:
   Solve a linear equation after performing QR factorization

 MAJOR TOPICS:
   Linear Systems

 CALLING SEQUENCE:
   X = QRSOLV(A, R, B, PIVOTS=IPVT)

 DESCRIPTION:

  The procedure QRSOLV completes the solution of a linear equation,

        A ## x = B

  after the MxN matrix has been factorized by QR decomposition.
  After being factorized once using QRFAC, the matrices can be used
  for multiple righthand sides (i.e., different B's).

  The solution technique is to first compute the factorization using
  QRFAC, which yields the orthogonal matrix Q and the upper
  triangular matrix R.  [ Actually, Q is represented by its
  Householder reflectors. ]  Then the solution vector, X, is computed
  using QRSOLV.

  If pivoting was performed in the factorization, the permutation
  vector IPVT returned by QRFAC must also be passed to QRSOLV.
  

 PARAMETERS:

   A - upon input, the factorized matrix A, returned by QRFAC.

   R - upon input, the upper diagonal matrix R, returned by QRFAC.

   B - upon input, the righthand vector B, which fits into the
       equation,  A ## x = B

   X - upon ouptut, the solution vector X, to the above linear
       equation.  For an overdetermined system, X is the least
       squares solution which minimizes TOTAL( (A ## X - B)^2 ).


 KEYWORD PARAMETERS:

   PIVOTS - upon input, the permutation matrix IPVT returned by
            QRFAC, if pivoting is to be performed.


 EXAMPLE:

  Solve the equation A ## X = B, in the least squares sense, where:

    A = [[1.0,1.0,1.0,1.0,1.0,1.0],$
         [0.6,0.8,0.5,0.8,0.7,0.9],$
         [0.2,0.3,0.1,0.4,0.3,0.4]]

  and B = [0.57E,0.69,0.5,0.7,0.6,0.8]

  qrfac, a, r, ipvt, /PIVOT
  x = qrsolv(a, r, b, PIVOTS=ipvt)

  print, x
       0.0834092     0.852273    -0.179545

 REFERENCES:

   More', Jorge J., "The Levenberg-Marquardt Algorithm:
     Implementation and Theory," in *Numerical Analysis*, ed. Watson,
     G. A., Lecture Notes in Mathematics 630, Springer-Verlag, 1977.

 MODIFICATION HISTORY:
   Written (taken from MPFIT), CM, Feb 2002
   Usage message, error checking, CM, 15 Mar 2002
   Error checking is fixed, CM, 10 May 2002

  $Id: qrsolv.pro,v 1.3 2002/05/10 18:41:31 craigm Exp $

(See C:\RSI\User\HarveyLib\qrsolv.pro)


RAMP

[Previous Routine] [Next Routine] [List of Routines]
 FUNCTION RAMP,X
 Computes the standard RAMP function for any list of values represented
 by the vector x.

 USAGE:
 x=FINDGEN(1001)/20-25 CREATE A VECTOR ON THE INTERVAL [-25,25]
 y=RAMP((x-8)/5) Construct a function that is shifted and stretched

 HISTORY
 Written by Harvey Rhody, September 1997

(See C:\RSI\User\HarveyLib\RAMP.PRO)


RANDOMQ

[Previous Routine] [Next Routine] [List of Routines]
z=RandomQ(Q,N) returns an array z of normal random vectors
with N columns which have the covariances defined by Q.

OPTIONAL PARAMETERS

SEED is a keyword specifiying the seed to the random number generator.
MEAN is a keyword specifying the desired mean vector. Default is 0.

NOTES:
Q must be a positive definite symmetric matrix.
Uses RANDOMN and CHOLDC.

Harvey Rhody
November, 2002.

(See C:\RSI\User\HarveyLib\randomq.pro)


RANK

[Previous Routine] [Next Routine] [List of Routines]
rankA=Ortho_Rank(A,TAU=tau) is the rank of a matrix A.
If A is a scalar then rankA=0. tau is the threshold
for zero values. Default: tau=1e-6.

Uses QR,A,Q,R and counts the number rows of R
that contain at least one value GE tau.

Arrays H,R,K such that A=H##R##Transpose(K) are
returned through the arguments if ortho_rank is called
with rankA=Ortho_Rank(A,TAU=tau,H,R,K)

H. Rhody
June, 2004

(See C:\RSI\User\HarveyLib\rank.pro)


RDPIX

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	RDPIX

 PURPOSE:
	Interactively display the X position, Y position, and pixel value
	of the cursor.

 CATEGORY:
	Image display.

 CALLING SEQUENCE:
	RDPIX, Image [, X0, Y0]

 INPUTS:
	Image:	The array that represents the image being displayed.  This
		array may be of any type.  Rather reading pixel values from
		the display, they are taken from this parameter, avoiding
		scaling difficulties.

 OPTIONAL INPUT PARAMETERS:
	X0, Y0:	The location of the lower-left corner of the image area on
		screen.  If these parameters are not supplied, they are
		assumed to be zero.

 OUTPUTS:
	None.

 COMMON BLOCKS:
	None.

 SIDE EFFECTS:
	The X, Y, and value of the pixel under the cursor are continuously
	displayed.

 RESTRICTIONS:
	None.

 PROCEDURE:
	Instructions are printed and the pixel values are printed as the
	cursor is moved over the image.

	Press the left or center mouse button to create a new line of output,
	saving the previous line.

	Press the right mouse button to exit the procedure.

 MODIFICATION HISTORY:
	DMS, Dec, 1987.
	Rob Montgomery (rob@hao.ucar.edu), 9/21/92;
		Correct indices for case of !order = 1

(See C:\RSI\User\HarveyLib\imgprobe.pro)


RECT

[Previous Routine] [Next Routine] [List of Routines]
 FUNCTION RECT,X
 Computes the standard RECT function for any list of values represented
 by the vector x.

 USAGE:
 x=FINDGEN(1001)/20-25 CREATE A VECTOR ON THE INTERVAL [-25,25]
 y=RECT((x-8)/5) Construct a function that is shifted and stretched

 HISTORY
 Written by Harvey Rhody, September 1997

(See C:\RSI\User\HarveyLib\RECT.PRO)


RECT2D

[Previous Routine] [Next Routine] [List of Routines]
 Z=RECT2D(X,Y)
 Computes the 2D rect function for the arrays X and Y, The
 result is unity wherever the corresponding values in the
 X and Y array are within a unit rectangle centered on the origin.

 X and Y must be arrays that are of the same size.

 USAGE
 Compute X and Y arrays using MESHDOM
 s=INDGEN(101)/10 -5
 t=INDGEN(101)/10 -5
 MESHDOM,s,t,X,Y
 Z=RECT2D(X/2,Y/3)
 SURFACE,Z,X/2,Y/3

 HISTORY
 Written by Harvey Rhody, September 1997

(See C:\RSI\User\HarveyLib\Rect2d.pro)


REFLECTION

[Previous Routine] [Next Routine] [List of Routines]
khat=REFLECTION(k,kc,ncols,nrows) is the reflection of image set k
about the origin represented by kc. The number of columns and rows in the
base image must be provided. Keywords /row and /column may be set to
reflect the image around the row or column that contains kc instead of the point
kc.

If the reflected set is empty then khat=-1 is returned.

HISTORY
Default reflection written by HR Sept. 1999  for SIMG-782
Modified to reflect about row and column by JH Oct. 2000 for SIMG-782

(See C:\RSI\User\HarveyLib\reflection.pro)


RUBBERBOX

[Previous Routine] [Next Routine] [List of Routines]
Creates a widget in which a box can be drawn using the
mouse. On exit the opposite corners of the box are
returned

RESTRICTIONS
Works only with monochrome images.

Does not rescale the brightness range.

H. Rhody
July, 2004

(See C:\RSI\User\HarveyLib\rubberbox.pro)


RUBBERLINE

[Previous Routine] [Next Routine] [List of Routines]
Creates a widget in which a line can be drawn using the
mouse. On exit the start and endpoints of the line are
returned

RESTRICTIONS
Works only with monochrome images.

Does not rescale the brightness range.

H. Rhody
July, 2004

(See C:\RSI\User\HarveyLib\rubberline.pro)


SAMPLE_CLICK

[Previous Routine] [Next Routine] [List of Routines]
P=SAMPLE_CLICK(Q) uses the mouse cursor to show values of an array Q that
has been displayed in the current window. Left click in the image to print the
coordinates and the value. Right-click to end the procedure. It is assumed
that the window is sized to fit the array image.

KEYWORDS
COORD='DEVICE','DATA','NORMAL' Default is DEVICE. If DATA or NORMAL are chosen
then the returned array P contains only the click coordinates, not the image values.

TRUE=[0,1,2,3]. Default is TRUE=0, which is a 2D image. For a 3D image
TRUE specifies the stacking coordinate. If TRUE =1,2,3 then the returned value
contains a vector of values for each pixel.

HISTORY:
Written by H. Rhody, Sept. 7, 1999

Converted to a function that returns the
coordinates and values. Nov 12, 2003 HR
COORD keyword October 29,2004
Added TRUE June 15, 2005.

(See C:\RSI\User\HarveyLib\sample_click.pro)


SGN

[Previous Routine] [Next Routine] [List of Routines]
 FUNCTION SGN,X
 Computes the standard signum function for any list of values represented
 by the vector x.

 USAGE:
 x=FINDGEN(1001)/20-25 CREATE A VECTOR ON THE INTERVAL [-25,25]
 y=SGN((x-8)/5) Construct a function that is shifted and stretched

 HISTORY
 Written by Harvey Rhody, September 1997

(See C:\RSI\User\HarveyLib\SGN.PRO)


SINC

[Previous Routine] [Next Routine] [List of Routines]
 FUNCTION SINC,X
 Computes the standard SINC function for any list of values represented
 by the vector x. See Gaskill, Page 47.

 USAGE:
 x=FINDGEN(1001)/20-25 CREATE A VECTOR ON THE INTERVAL [-25,25]
 y=SINC((x-8)/5) Construct a function that is shifted and stretched

 HISTORY
 Written by Harvey Rhody, September 1997

(See C:\RSI\User\HarveyLib\Sinc.pro)


SINC2D

[Previous Routine] [Next Routine] [List of Routines]
 Z=SINC2D(X,Y)
 Computes the 2D SINC function for the arrays X and Y,

 X and Y must be arrays that are of the same size.

 USAGE
 Compute X and Y arrays using MESHDOM
 s=FINDGEN(101)/10 -5
 t=FINDGEN(101)/10 -5
 MESHDOM,s,t,X,Y
 Z=TRI2D(X/2,Y/3)
 SURFACE,Z,X/2,Y/3

 HISTORY
 Written by Harvey Rhody, September 1997

(See C:\RSI\User\HarveyLib\Sinc2d.pro)


SOMB

[Previous Routine] [Next Routine] [List of Routines]
 Z=SOMB(X,Y)
 Computes the sombrero function for the arrays X and Y,
 See Gaskill, page 72. NOTE: Gaskill uses the radius R
 as the input. This function uses the rectulangular arrays
 X and Y. Internally, R is computed from X and Y.

 X and Y must be arrays that are of the same size.

 USAGE
 Compute X and Y arrays using MESHDOM
 s=FINDGEN(101)/10 -5
 t=FINDGEN(101)/10 -5
 MESHDOM,s,t,X,Y
 Z=SOMB(X/2,Y/3)
 SURFACE,Z,X/2,Y/3

 HISTORY
 Written by Harvey Rhody, September 1997

(See C:\RSI\User\HarveyLib\SOMB.PRO)


STEP

[Previous Routine] [Next Routine] [List of Routines]
 FUNCTION STEP,X
 Computes the standard step function for any list of values represented
 by the vector x.

 USAGE:
 x=FINDGEN(1001)/20-25 CREATE A VECTOR ON THE INTERVAL [-25,25]
 y=STEP((x-8)/5) Construct a function that is shifted and stretched

 HISTORY
 Written by Harvey Rhody, April, 2000

(See C:\RSI\User\HarveyLib\step.pro)


STRING2BINARY

[Previous Routine] [Next Routine] [List of Routines]
b=string2binary(s) returns the binary ASCII representation of the characters in string s

(See C:\RSI\User\HarveyLib\string2binary.pro)


TELEGRAPHPROCESS

[Previous Routine] [Next Routine] [List of Routines]
 S=TelegraphProcess(t,lamda,p)
 Generates a random process with X[t]=X[0] if
 the number of zero crossings in (0,t) is even and
 X[t]=-X[0] if the number is odd. The waveform is
 returned in an array of size n=lambda*t/p. If
 p is not specified then the default is p=0.1. The
 probability should be small enough that the Poisson
 assumptions are satisfied. The initial value X[0] can
 be specified with the keyword X0. If it is unspecified,
 it is randomly initialized to +1 or -1.

 Harvey Rhody
 March, 2000

(See C:\RSI\User\HarveyLib\TelegraphProcess.pro)


TRANSFORM_BASIS

[Previous Routine] [Next Routine] [List of Routines]
A=TRANSFORM_BASIS(N,NAME) returns the NxN matrix A
whose rows are the 1D basis vectors of length N defined
for the transform selected by the string NAME.

A is unitary with A##A^H=A^H##A=I, where A^H is the
conjugate transpose of A.

NAME may be one of the following (not case sensitive):
'DCT' for the discrete cosine transform
'DFT' for the DFT

H. Rhody,  July, 2000

(See C:\RSI\User\HarveyLib\transform_basis.pro)


TRI

[Previous Routine] [Next Routine] [List of Routines]
 FUNCTION TRI,X
 Computes the standard TRI function for any list of values represented
 by the vector x.

 USAGE:
 x=FINDGEN(1001)/20-25 CREATE A VECTOR ON THE INTERVAL [-25,25]
 y=TRI((x-8)/5) Construct a function that is shifted and stretched

 HISTORY
 Written by Harvey Rhody, May, 2000

(See C:\RSI\User\HarveyLib\Tri.pro)


TRIG2D

[Previous Routine] [Next Routine] [List of Routines]
z=Trig2D(n,m) returns a spatial cosine with n horizontal
and m vertical cycles. The default size is 256x256.
The size can be set to any integer value using
keywords XSIZE and YSIZE.

z=Trig2D(2,3,XSIZE=100,YSIZE=200) returns an array
of size 100 columns and 200 rows with 2 horizontal and
three vertical cycles.

EXAMPLE
The results could be computed and displayed with
n=2 & m=3 & xsize=100 & ysize=200
z=Trig2D(n,m,XSIZE=xsize,YSIZE=ysize)
txt=STRING(FORMAT='("N=",I2," M=",I2)',n,m)
WINDOW,/FREE,XSIZE=xsize,YSIZE=ysize,TITLE=txt
TVSCL,z

H. Rhody, July 2000

(See C:\RSI\User\HarveyLib\Trig2D.pro)


TRI_COEF

[Previous Routine] [Next Routine] [List of Routines]
c=tri_coef(X,Y) returns an interpolation vector of length
three that describes a linear interpolation at a known point
[X[3],Y[3]] when the values [Z[0],Z[1],Z[2]] are given at
points (X[0],Y[0]), (X[1],Y[1]), (X[2],Y[2]).

Vectors X and Y must contain four components, namely, the
coordinates of the corners of a triangle and the coordinates
of another point, presumably within the triangle.

The notion is that this program will be used to compute
coefficients to a known grid from a set of known DeLaunay
triangles. This enables us to pre-compute the triangles
and insert the Z values at the end.

H. Rhody
October 5, 2004

(See C:\RSI\User\HarveyLib\tri_coef.pro)


TRUNCATED_RAMP

[Previous Routine] [Next Routine] [List of Routines]
 FUNCTION TRUNCATED_RAMP,X
 Computes a truncated ramp function for any list of values represented
 by the vector x. The ramp is zero outside the range [-b,b]

 USAGE:
 x=FINDGEN(1001)/20-25 CREATE A VECTOR ON THE INTERVAL [-25,25]
 y=TRUNCATED_RAMP((x-8)/5,3) Construct a function that is shifted and stretched

 HISTORY
 Written by Harvey Rhody, October, 2000

(See C:\RSI\User\HarveyLib\Truncated_Ramp.pro)


UNION

[Previous Routine] [Next Routine] [List of Routines]
C=UNION(A,B) is the union of the sets A and B. It is assumed
that A and B have some form of integer values. That is, they
come from the set of non-negative integers.

(See C:\RSI\User\HarveyLib\union.pro)


WASP_WRITE_HDR

[Previous Routine] [Next Routine] [List of Routines]
Reads the header file associated with a .img file and
returns the number of columns, number of rows, and the
data type.

(See C:\RSI\User\HarveyLib\wasp_write_hdr.pro)


ZEROPAD

[Previous Routine] [Next Routine] [List of Routines]
B=zeroPad(A,N,M) returns an array B of size NxM such that B[i,j]=A[i,j]
where A is defined and B[i,j]=0 where A[i,j] is not defined. B is the same
type as A.

If A is a scalar or a row vector then zeroPad can be called with a single
dimension argument. B=zeroPad(A,N) returns a row vector of length N.

It is an error if N or M is smaller than the corresponding dimension of A.

KEYWORDS

CENTER: Set the keyword center if the array should be centered in the
field of zeros.

HISTORY
Written by H. Rhody September, 1999
Modified to include additional error checks and CENTER keyword
January 2001.

(See C:\RSI\User\HarveyLib\zeropad.pro)


ZSHIFT

[Previous Routine] [List of Routines]
B=zshift(A,dx,dy) is the same as shift(A,dx,dy) except that the
vacated pixels are filled with zeros. Works only with 2D.

(See C:\RSI\User\HarveyLib\zshift.pro)