;FUNCTION GetPixelRayCam
;  Purpose
;    Get the pixel rays for all pixels in a camera, in camera coordinates.
;  Input
;    cam       - Camera number. 0=PX, 1=PY, 2=MX, 3=MY
;    caminfo   - camera info structure as created by LoadCamInfo()
;    /grid     - If not set (default) treats each ray as the center of a pixel. The
;                outside rays will be set in from the edge of the FOV by 1/2 pixel.
;                If set, treats each ray as a node on a grid covering the entire FOV.
;                The outer rays on all four sides will be exactly on the edge of 
;                the FOV.
;  Return
;    a 2D grid of vectors (3D array, [col,row,vec_component]). Each is
;    a normalized vector pointing in the direction a particular pixel looks,
;    in Camera coordinates.
;    Camera frame has axis pointing at +z, increasing col pointing at +y
;    Increasing row pointing at -x
;
;  Distortion model: Each vector is pointed at
;    [kx*tan(thetax)+poly(d),ky*tan(thetay)+poly(d),1]
;   where
;    thetax is angle along horizontal from camera axis to pixel column
;    thetay is angle along vertical   from camera axis to pixel row
;    kx and ky are the appropriate scale factors for each axis
;    d=sqrt(thetax^2+thetay^2) in pixels
;    poly(theta) is some polynomial correction factor
function Getraypixel,Cam,caminfo,ray,grid=grid
  dist=[0,0,0,0,0] ;[constant,Linear,quadratic,cubic,4th order,etc...]
                   ;This is a null correction - IE camera is perfectly
                   ; represented by tan(theta)
  cols=(*(caminfo.cols))[cam];
  rows=(*(caminfo.rows))[cam];

  if ~keyword_set(grid) then begin
    CamFovPerPixX=(*(caminfo.thfovx))[cam]/((cols)/2d);
    CamFovPerPixY=(*(caminfo.thfovy))[cam]/((rows)/2d);
  end else begin
    CamFovPerPixX=(*(caminfo.thfovx))[cam]/((cols-1)/2d);
    CamFovPerPixY=(*(caminfo.thfovy))[cam]/((rows-1)/2d);
  end

  resolve_grid,ray,x=x,y=y,z=z

  x=x/z
  y=y/z
; z=z/z ==1  

  ;This is correct, there is an axis swap from camera to image
  i=y/CamFovPerPixX+cols/2d -0.5d
  j=x/CamFovPerPixY+rows/2d -0.5d
  return,compose_grid(i,j)
end