\documentclass[psfig,preprint]{aastex}
\begin{document}
 
\title{1d IMAGE EXPLORATION AND PROCESSING USING {\tt tv} AND ITS ILK
 \\ \today} 
\author{Carl Heiles}

%\enlargethispage{2em}
\tableofcontents

%\subsection{Examine the 1d Image on Your Terminal}                             

\newpage
This handout is oriented towards image exploration and
manipulation---i.e., image processing---rather than making nice
displays for presentation..

\section{PIXELS, COLORS, BYTES} \label{intro}

Your display screen consists of about a million little areas called {\it
  pixels}.  Each pixel can show a different color/intensity combination.
Everything on your screen---text, pictures, whatever---is displayed by
filling the appropriate pixels with the appropriate color/intensity.
Your screen dimensions might be $1280 \times 1024$ ($5 \times 4$ aspect
ratio)---found on many computers of all kinds; the somewhat larger $1440
\times 900$ ($8 \times 5$ aspect ratio); $1024 \times 768$ ($4 \times 3$
aspect ratio)---found on older computers and small laptops. These aspect
ratios come in much larger dimensions. Each pixel is small! But if you
look carefully, you can see them.

Typical CCD projectors (as used for PowerPoint) have the $4 \times
3$ format with $1024 \times 768$ pixels as their native resolution, so
you get best results, with no image cropping, when you set your video
card's resolution to this value. 

\section{DECOMPOSED COLOR (i.e., TrueColor) {\it vs.}\ COMBINED COLOR
  (i.e., PseudoColor)} \label{decomposed}

 	All colors seen by the human eye can be produced by a suitable
mixture of intensities of only three colors: red, green, and blue (RGB).
Most common displays in use today allow 256 intensities of each
color\footnote{Note that 256 is the same as $2^{8}$: it's 8 bits---or,
alternatively, 1 byte.}.  This gives a total of $256^3$
combinations---this used to be billed in the PC world as ``millions of
colors''.  This is called {\it TrueColor}.

Believe it or not, it is often desirable to degrade the full ``millions
of colors'' true scheme to a 256-color scheme; this is called {\it
  PseudoColor}. In PseudoColor, you are {\it combining} the separate RGB
colors into a single 256-element set; that's why it's called {\it
  Combined color}. The particular set of 256 colors chosen is called the
{\it colortable}; each of the 256 colors is a particular combination of
RGB.

The most common colortable is the grayscale one. In this, the
intensities of RGB are all identical. They change uniformly from zero
intensity (for entry number 0---this makes black) to full intensity (for
entry number 255---this makes white). In between, we have all shades of
gray. An image made with the grayscale colortable looks like a
black-and-white photograph.  If you choose a non-gray colortable, then
you have PseudoColor, which is often called ``false color'', which can
be helpful in highlighting certain features or achieving more
contrast. We give an example below, and other examples in our ``1d2d3d:
One, Two, and Three Dimensional Color Images'' memo.

\subsection{Combined Color: 256-Element Colortables}

In your computer's memory, an image consists of a two-dimensional array
of {\it data} values $d$, one for each pixel. These are actual
numbers. Their range is restricted: $d$ can take on 256 different
values, ranging from 0 to 255.  On your screen, the image consists of
projected light in each pixel, which we call the intensity $I$.

There is a one-to-one relationship between $d$ and $I$, so there are
also 256 different possible values for $I$. Generally speaking, this
relationship is specified by the ``color table''. The color table is
often nonlinear so as to emphasize weak or strong features. Sometimes
it's an equal mixture of red, green, and blue (grayscale); sometimes the
mixture is engineered to produce color. By ``Intensity'' ($I$), we mean
the 256 different possible combinations of light intensity and color,
one for each value of $d$. 

To use combined color in IDL, you must turn the native $256^3$ RGB
colors into a 256-element colortable, which you accomplish by turning
color decomposition off and then defining a colortable. You do the
former with the command

\verb$device, decomposed=0$

\noindent By default, IDL sets the three colors {\tt r,b,b} equal, which
gives you grayscale. If you want the system colors (like {\tt !red})
defined and you're using our IDL library, the uppermost
12 entries of the grayscale colortable define the colors\footnote{e.g.,
  number 244 is {\tt !black}, 245 is {\tt !red}, etc.; use IDL's {\tt
      tvlct} comomand with the {\tt get} keyword set to see these
    definitions.}, and you load
that colortable with the command

\verb$setcolors, /system_variables  (or setcolors, /sys )$

\noindent Alternatively, you could define a modified grayscale 
(e.g., as in equation \ref{gammamapping}) or a false-color colortable that is
not gray\footnote{IDL has a wide selection of such colortables; in
  IDL, to view them type {\tt xloadct} and to load one use {\tt
    loadct}.}.  

\subsection{Decomposed Color}

For ordinary non-image work you almost always want Decomposed Color so
you can make any color you want on the graphics output; you don't want
to be restricted to 256 colors. In TrueColor, you have 256 different
possibilities for each of the three colors RGB, so you have $256^{3}$
color possibilities. It's definitely {\it not} a 256-element colortable!
In IDL, you can turn color decomposition on with (guess what!)

\verb$device, decomposed=1$.

\noindent If you want the system colors (like {\tt !red}) defined and
you're using our IDL library, also type the command
\verb$setcolors, /system_variables  (or setcolors, /sys )$

\subsection {What Have I Done?}

\noindent You can oscillate between decomposed and combined color
at will, to your heart's content. In fact, you sometimes need to do this
when making color images in PostScript files, because images in
PostScript need decomposed color, while vector graphics (like symbols or
lines) need combined color. To see what you've told IDL, enter

\verb$ help, /device$

\noindent If you
are in TrueColor with decomposed color turned off (i.e., a 256-element
colortable) it will say ``Graphics
pixels: Combined''; if it isn't turned off it says ``Graphics pixels:
Decomposed''. It also shows the sizes of all graphics windows and tells
you which visual class you're running (TrueColor or DirectColor).  If
you are in PostScript mode, in which it writes all output to a
PostScript file instead of the screen, it will tell you so.

\section {COMBINED COLOR---256-ENTRY (8-BIT) COLOR
  TABLES}
\label{colortables}

{\it The remainder of this tutorial assumes you have turned color
  decomposition off, i.e.\ you are using combined color}\footnote{You do
  this with {\tt device, decomposed=0} and, if desired, {\tt
    setcolors,/sys}.}.  In this PseudoColor mode, the maximum number of
color/intensity combinations that can be displayed simultaneously is
256.  Therefore, images are represented by numbers that range $0
\rightarrow 255$.  For this reason, displayed images are {\it always}
represented by byte arrays.  You can display other array data types, but
IDL will convert whatever you give it to a byte array before displaying
it.  Therefore, if you display, say, an integer array (integers are two
bytes long and range from $-32768 \rightarrow 32767$), and if numerical
values in this array exceed 255, then the resulting image display will
look weird.  That's because, in converting from integer to byte, numbers
that exceed 255 will ``wrap around''.  For example, integer 255 equals
byte 255, but integer 256 equals byte 0, integer 257 equals byte 1, etc.
Below, we'll deal with these conversions in more detail.

	For now, let's restrict our attention to black/white
images---otherwise known as ``grayscale'' images.  Gray, or white, is
composed of an equal mixture of red, green, and blue, and all we deal
with is the intensity $I$.  In a grayscale image, the intensity of each
pixel is related to the data value $d$ in that pixel.  Let's think of
large intensity being white and small intensity being black; there are
256 different possible intensities, so $I$ can range from $0 \rightarrow
255$.  Similarly, the data values $d$ can range from $0 \rightarrow
255$. 

	An important concept is the relationship between $I$ and $d$. 
This is known as the {\it color table}. It specifies the mapping between
data value and color/intensity---or, for a grayscale image, the mapping
between data value and intensity. 

\subsection{ Linear Mapping, both Direct and Reversed} The simplest
mapping between data value $d$ and intensity $I$ is a linear one with

\begin{equation}
I = d
\end{equation}

\noindent In this case, a data value $d = 255$ gives white and $d = 0$
gives black.  This {\it direct mapping} is the default manner in which
images are displayed on the computer screen: there is a black background
on which the image is painted with increasing data values being
increasingly white.  However, on a piece of paper the relationship is
usually reversed, because paper is white and provides a naturally white
background.  Thus, in this {\it reversed mapping}, we want to paint the
image with increasing data values being increasingly black.  This is
also a linear mapping, but reversed:

\begin{equation} \label{reverse}
I = 255 - d
\end{equation}

\noindent {\it NOTE:} Printed images usually look {\it much} better with
the reversed mapping, because printers have a hard time giving a uniformly
black area with no streaks.  This is the {\it first} reason why printed
images should be made with a reversed mapping.  The {\it second} reason is
that in scientific journals, images with the reversed mapping are
reproduced much better. The {\it third} reason is that making the paper
black uses lots of printer toner, which is expensive.  To reverse the color
table, you can use equation \ref{reverse}. Alternatively, for a byte array
called \verb$img$, you can type \verb$tv, not(img)$ instead of
\verb$tv, img$.

\subsection{ Nonlinear Mapping} 

	The linear mapping is often not very useful because you usually
want to highlight weak features or bright features; we'll see an example
below.  The most commonly used nonlinear mapping uses a power law (this
is the photographer's ``characteristic curve'') together with a
``stretch'', which cuts off the image at dim and bright intensity
levels:

\begin{mathletters} \label{gammamapping}
\begin{equation}
I = 255 \left({d - d_{bot} \over d_{top} - d_{bot}}\right)^\gamma, \ d =
d_{bot} \rightarrow d_{top}
\end{equation}

\begin{equation}
I = 0, \ d \le d_{bot}
\end{equation}

\begin{equation}
I = 255, \ d \ge d_{top} 
\end{equation}
\end{mathletters}

\noindent In a reversed mapping, you'd substitute $(255-d)$ for $d$ in
the above equations. 

There is one other commonly used nonlinear mapping, the so-called
``histogram equalization'' technique.  In this technique, the mapping is
modified on an image-by-image basis so that, all of the 255 colors are
used in an equal number of pixels.  Read about it in IDL's documentation
on {\tt hist\_equal}.

\subsection{Why We Need Nonlinear Mapping}

{\it Never forget} that the idea is to turn the data array into
an image that conveys information to the brain.  The idea is {\it not}
to be so strictly quantitative that details of interest are obscured.

You want to bring out details of interest.  For example, for many images
of the interstellar gas you want to emphasize weak structures at the
expense of the fidelity gained by a strict proportionality between image
brightness and data value.  To this end, choose a color table and
experiment with the image transfer function.  At minimum, this involves
changing the span of the data values represented in the image and
raising the data values within that span to a power: a power less than
unity to emphasize weak features, larger than unity for strong ones.

\section{DISPLAYING AN IMAGE IN IDL}

{\it We are still assuming that you have turned color decomposition off,
  i.e.\ you are using a 256-entry color table.}  To accomplish this
feat, see \S \ref{decomposed}.

There are two ways to put the image on your computer screen (in
computerese: ``write the image onto your X window''). One is simpler and
uses Tim Robishaw's (TR's) {\tt display} procedure; it is best for
producing publication-quality images as described in our handout {\tt
  tr\_display.pdf}. Here we describe the other, which uses IDL's native
procedures; it's a bit more cumbersome but is more useful when you are
manipulating images, e.g.\ when you absolute knowledge of which pixel is
which.

\subsection{An Example: Read the Image from a FITS File}

First, generate an image.  For this example there's a nice image of the
X-ray sky {\tt rass\_c.fits}, obtained by the ROSAT satellite on the web
at: \hfill \break {\tt
  http://astro.berkeley.edu/$\sim$heiles/handouts/handouts\_images.html}
\break \noindent or \break {\tt
  http://ugastro.berkeley.edu/radio/2016/index.html} \hfill \break Copy
this file to the disk area where you are running IDL.  This file is in a
format called ``fits'' format (``flexible image transport system''),
which is the same format of many astronomical images.

To read the data file into an array called {\tt image}, it is easiest to
use the IDL procedure called ``mrdfits'', which resides in the Goddard
IDL library which, in turn, is already in your IDL path.  If you are
logged into {\tt ugastro}, all you have to do is type

{\tt image = mrdfits( '/home/global/ay121/handouts/images/rass\_c.fits', 0, hdr)}

\noindent This returns two arrays: the image array ({\tt image}) and
information about the image ({\tt headerinfo}); type {\tt print,
headerinfo} to see the header information.  Now type {\tt help, image}
and IDL will tell you that it is a $480 \times 240$ FLOAT array.

\subsection{Displaying and Examining the Image with IDL's Native
  Procedures Like {\tt tv}} \label{native}

Now we're ready to write the image to the X window. If all  you want to
do is look at it, the easy thing to do is use {\tt display, image}. But
when you want to examine the contents of individual pixels on the image,
i.e., when you want to do image processing, you want each image pixel to
occupy a single pixel on your screen. 

\subsubsection{Image Size and Window Size} \label{windowsize}

This image is $480 \times 240$. Each element will occupy a single pixel on
the screen, so we need a window of at least that size to display the whole
image. We can use a bigger window, in which case the image won't fill the
window area. If we use a smaller window, only part of the image will be
displayed.  We can create a window of the appropriate size, that is with
numbers of pixels equal to the same dimensions of the data array, by typing

{\tt window, xsize=480, ysize=240}

\noindent and then displaying it with {\tt tv}.

Maybe you'd like a bigger image in a bigger window so you can see things
more clearly. Or maybe the image is too large for your screen and you
need to make it smaller so it fits.  In either case, you need to change
the size of the image as measured in pixels.  IDL does this easily with
the {\tt rebin} function. Suppose you want to increase the size by a
factor of 2 in the horizontal and 3 in the vertical
direction\footnote{Using different factors for horizontal and vertical
  changes the aspect ratio, which is usually a bad idea; we do it here
  simply for illustration.}, i.e.\ to make an array of size $960 \times
720$.  Do this by

{\tt bigimage = rebin(image, 960, 720)}

\noindent Then create an appropriately-sized window (e.g.\ with {\tt
window, 5, xsize=960, ysize=720}); this creates a new window, numbered
5, and leaves the old ones in place.

You can also resize the image using {\tt congrid}, which works for non-integral
factors. You might say, ``Well, I'll always use {\tt congrid}---it's more
flexible''. {\it But be careful!} {\tt congrid} and {\tt rebin} handle
enlargement and ensmallment differently, and treat the edges differently.
With {\tt congrid}, you {\it almost certainly} don't want to use the
default options; look carefully at the keywords and try them out on a short
1-d array to see their effects. Usually, {\tt rebin} is better; don't use
{\tt congrid} unless you know what you're doing.
 
\subsubsection{Displaying the Image with {\tt tv}} \label{displaying}

Display the image by typing

{\tt tv, image}

\noindent and you see a gray mishmash oval.  The oval is the Aitoff
projection of the entire sky in soft X-rays.  The mishmash occurs
because the data values in {\tt image} exceed the allowable $0
\rightarrow 255$ range of a byte array, so there's lots of wrapping. You
can use the {\tt max} and {\tt min} functions (or, nicer, Goddard's {\tt
  minmax} function) to determine that the data values range from about
$-174 \rightarrow 45337$, thus far exceeding the valid range for a byte
array.

	You can scale the data so that they all fit in the allowable
byte range $0 \rightarrow 255$.  We'll first produce a byte array, which
we'll call {\tt byteimage}, from {\tt image}\dots

{\tt byteimage = bytscl(image)}

\noindent This linearly scales {\tt image}, which ranges $-174
\rightarrow 45337$, into {\tt byteimage}, ranging from $0 \rightarrow
255$.  To display this image\dots

{\tt tv, byteimage}

\section{EXPLORING THE IMAGE WITH {\tt trc} AND {\tt profiles}}

What do you see in this Image? All you see is two white dots! These two
dots are the strongest X-ray sources in the sky---the one on the left is
a point source called ``Cygnus XR-1'', and the one on the right is the
Vela supernova remnant, home of the famous ``Vela pulsar''. Where are
these sources located? Use {\tt trc} to find the pixel coordinates of
these sources.  {\tt trc} is just a wrapper that calls {\tt tr\_rdplot}
with useful settings. Type

\noindent {\tt trc, x, y}

\noindent Then move the cursor inside the plot window. Left clicks write
and record the {\tt (x,y)} values in the coordinate system you select
(i.e., {\tt data, device, normal}; default is {\tt data}). A middle or
left click will exit the program.  Notice that the crosshairs are
clipped at the axes: it's often useful to make the crosshairs extend to
the edges of the window.  This is accomplished by passing the {\tt
  /noclip} keyword, and since we can use minimum matching when calling
keywords in IDL, {\tt /noc} will be good enough:

\noindent {\tt trc, x, y, /noc}

\noindent If you don't want to store the positions of your left-clicks
in {\tt x} and {\tt y}, then you don't need to include them in the
procedure call. For more details, see the
documentation of {\tt trc} and {\tt tr\_rdplot}.

%vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
These images contain much more than just these two sources! To see more,
we need to change the contrast---change the `dynamic range'---by
invoking a nonlinear mapping of data $d$ to screen intensity $I$.  In
particular, we need to change one or more of ($d_{bot}, d_{top},
\gamma$).  How much should we expand the dynamic range? We might make a
guess and try $d_{bot}=-174$ and $d_{top}=2000$.  If that didn't give a
nice result, we could try some other values.

But we don't have to guess! IDL provides several nice ways to
interactively sample the image.  You can get a quick feel for the
interesting data range by just doing {\tt plot, image} and visually
estimating the range of interest. More instructive is to sample the
image itself using {\tt profiles}, which plots horizontal or vertical
cuts of the image. Type

{\tt profiles, byteimage}

\noindent \noindent and follow the printed instructions; it plots
profile cuts across the image, horizontally or vertically, at positions
selected by the cursor. These profile plots are for the byte image.
More useful is plotting of the {\it original} data array {\tt image}:

{\tt profiles, image}

\noindent which plots the original numbers. Notice that a new window has
been created.  Move the cursor inside the plot window containing the
displayed image. You now see the profiles in the newly created window
(labeled {\it Profiles}).  Initially, you see intensity versus the {\tt
  x} data coordinates. Left click toggles betwwen horizontal/vertical;
right click exits the program and closes the profiles window.

In our example, this plot is so compressed that it is virtually
worthless, because the plot automatically scales to the minimum and
maximum values of the array; you can get around this easily by using the
${\tt <}$ and ${\tt >}$ operators, or by setting the system variable
{\tt !y.range} (easiest way: use our {\tt ver} procedure, like this:

{\tt ver, 0, 2000}

{\tt profiles, image}

{\tt ver}

After some inspection we see that limiting the data range to $200
\rightarrow 1400$ with $\gamma=.003125$
would indeed be a good start.

%gamma=0.125/4.
%bytimg= bytscl( img^gamma, min=200.^gamma, max=1400.^gamma)
%.........................................................
\section{MANIPULATING THE IMAGE BY MANUALLY CHANGING THE COLORTABLE}
\label{static}

%\section{REDISPLAY USING THE MODIFIED COLORTABLE}

Now, having determined suitable values for $d_{bot}$, $d_{top}$, and
$\gamma$, we want to display the appropriately-scaled image.  To display
the data range and gamma mentioned just above, we again use the {\tt
  bytscl} command as above but limit the data range by typing\dots

{\tt byteimage = bytscl( image\^\ gamma, min=200\^\ gamma, max=1400\^\ gamma)} 

Now display this with

{\tt tv, byteimage}

\noindent and, instead of just the two strong sources, you see beautiful
diffuse emission with lots of structure, which is produced
by hot, $\sim 10^6$ K gas heated by the expanding shock of the supernova
remnant.

\section{ EXPLORE AND HAVE SOME FUN WITH THIS IMAGE!}

There's lots more in image processing and display.  Of course, you can
manipulate images mathematically, just as you can any other IDL variable
or array---but remember to manipulate the {\it original} array instead
of its {\it byte} counterpart.  For detailed pixel-by-pixel sampling,
try {\tt rdpix, image}.  You can make a histogram [e.g.\ {\tt histo =
    histogram(image)}] of the original image; this tells you where most
of the brightness data are concentrated.  And if you're interested in
only a {\it portion} of the image, you can select this portion with the
cursor using {\tt defroi} (``defiine region of interest'').  You can do ``histogram equalization'' with
{\tt hist\_equal}; this is a simple way to see most of the interesting
structure: 

{\tt tv, hist\_equal(image)}

\noindent You can rotate with {\tt rotate} or {\tt rot}, {\tt
  transpose}, {\tt zoom}, draw {\tt contour}s, label your images and
make coordinates using {\tt plot} (with the {\tt /noerase} keyword) and
{\tt xyouts}, etc., etc., etc.

Play around with contrast by experimenting with $d_{bot}, \, d_{top}$,
  and especially {\tt gamma}. When you increase the sensitivity to small
  numbers by using $\gamma < 1$, you see the diffuse background. See
  that huge circular structure in the middle?  That's the ``North Polar
  Spur''.  It occupies an angle of about $120^\circ$.  It's
  close---almost touching our noses!  It's caused by several dozen
  supernovae that have exploded, producing a ``superbubble''.  These
  supernovae were located in the large cluster of young stars in the
  Scorpio constellation---some of the stars you see there on a dark
  night will explode as supernovae some day, adding to the energy stored
  in the hot gas and brightening the X-ray emission.  You also can see a
  bunch of fairly weak point sources and other diffuse structures.

To learn the range of native IDL capabilities, enter IDL and type ``?'',
which brings up the online help window. Go to the bottom left and click
on ``Contents'', then ``Routines (by topic)'', and take a look at the
sections entitled ``Array Manipulation'', ``Color Table Manipulation'',
``Direct Graphics'', ``Image Processing'', ``Plotting'', ``Signal
Processing''.

\end{document}
\end
