function ideproject,dimage,point,negative=negative,missing=missing, $
                    bytscale=bytscale,cubic=cubic,nointerp=nointerp

     nx_image = n_elements(dimage(*,0,0))
     ny_image = n_elements(dimage(0,*,0))
     xcorners = [0,0,nx_image-1,nx_image-1]
     ycorners = [0,ny_image-1,0,ny_image-1]
     
     if tag_index(point,'LAT') LT 0 then begin
        pindex = tag_index(point,'POINT')
        if pindex LT 0 then message,'Point structure is incomplete'
        upoint = point.(pindex)
     endif $
     else upoint=point

     dsize = size(dimage)
     if dsize(0) EQ 3 then nimage = dsize(3) $
     else if dsize(0) EQ 2 then nimage = 1 $
     else message,'Image array does not have the correct number of dimensions'

     if nimage NE n_elements(upoint) then begin
        if n_elements(upoint) eq 1 then begin
           upoint = replicate(upoint,nimage) 
           message,/info,'WARNING: point structure has only one element, replicating to match the number of images in the data cube.'
        endif else $
           message,'The point structure must have the same number of elements as the number of images in the data cube'
     endif

     if tag_index(upoint,'B0') LT 0 OR tag_index(upoint,'P') LT 0 OR $
        tag_index(upoint,'LAT') LT 0 OR tag_index(upoint,'CMD') LT 0 then $
        message,'Point structure is incomplete'

     newpoint = upoint

     for i=0L,nimage-1L do begin

        ;message,/info,'Image '+string(i)

        ; Set up the deprojection transformation

        b0 = double(upoint(i).b0)
        p  = double(upoint(i).p)
        if tag_index(upoint,'PIX_SIZE') GE 0 AND $
           tag_index(upoint,'RADIUS') GE 0 then begin
           latlon = get_latlon(nx=nx_image,ny=ny_image,point=upoint(i))
           bc = double(latlon.latitude) * !dtor
           lc = double(latlon.cmd) * !dtor
        endif else begin
           bc = double(upoint(i).lat)
           lc = double(upoint(i).cmd)
           message,/info,'WARNING: PIX_SIZE and/or RADIUS not set.  Treating image as a plane.'
        endelse

        cxx = - SIN(b0)*SIN(p)*SIN(lc) + COS(p)*COS(lc)
        cxy = - SIN(bc)*( SIN(b0)*SIN(p)*COS(lc) + COS(p)*SIN(lc) ) $
              - COS(bc)*COS(b0)*SIN(p) 
        cyx = + SIN(b0)*COS(p)*SIN(lc) + SIN(p)*COS(lc)
        cyy = + SIN(bc)*( SIN(b0)*COS(p)*COS(lc) - SIN(p)*SIN(lc) ) $
              + COS(bc)*COS(b0)*COS(p)
 
        determ = cxx*cyy-cyx*cxy  ; Get image to helio transformation
        icxx = cyy/determ         ; The determinent is double precision
        icxy = -cxy/determ
        icyx = -cyx/determ
        icyy = cxx/determ
 
        dXdxi = icxx
        dYdxi = icyx
        dXdyi = icxy
        dYdyi = icyy
 
        if tag_index(upoint,'PIX_SIZE') LT 0 then begin
           pix_size = [1.d0,1.d0]
        endif else begin
           pix_size = double(upoint(i).pix_size)
           if n_elements(pix_size) EQ 1 then $
              pix_size = [pix_size(0),pix_size(0)]
        endelse
        mpix_size = min(pix_size)
 
        x = (lindgen(nx_image,ny_image) MOD nx_image)*(pix_size(0)/mpix_size); image pixel coordinates
        y = (lindgen(nx_image,ny_image) / nx_image)*(pix_size(1)/mpix_size) 
        b_pxd = dXdxi*x + dXdyi*y           ; Heliographic pixel coordinates
        b_pyd = dYdxi*x + dYdyi*y


        nx = fix((max(b_pxd)-min(b_pxd)+1))
        ny = fix((max(b_pyd)-min(b_pyd)+1))

        if i EQ 0 then begin 
           newimage=make_array(nx,ny,nimage,type=dsize(dsize(0)+1))
           newnx = nx
           newny = ny
        endif 
        xoffset = (newnx - nx)/2.0
        yoffset = (newny - ny)/2.0

        ; p,q will be the forward poly_2d transformation 
        if n_elements(missing) NE 1 then $
           missing =(max(dimage(*,*,i))+min(dimage(*,*,i)))/2.
        polywarp,x-xoffset,y-yoffset,b_pxd-min(b_pxd), $
                 b_pyd-min(b_pyd),3,p,q
        if keyword_set(cubic) or NOT keyword_set(nointerp) then begin
           useinterp =1+keyword_set(cubic) 
           if keyword_set(cubic) then cval = -0.5 ; otherwise leave undefined
        endif else begin
           useinterp = 0
        endelse
        newimage(*,*,i) = poly_2d(dimage(*,*,i),p,q,useinterp, $
                                  newnx,newny,missing=missing,cubic=cval)

        if keyword_set(bytscale) then begin
           ; The bytscl call ensures that the image will not use
           ; the first 15 colors which will be reserved for contours.
           newimage(*,*,i)= (15+bytscl(newimage(*,*,i),top=(!d.table_size-1-15)))
           if keyword_set(negative) then begin
              newimage(*,*,i) = max(newimage(*,*,i))   $
                                + min(newimage(*,*,i)) $
                                - newimage(*,*,i)
           endif
        endif
     endfor

     return,newimage

end