;+
; $Id: plane_project.pro,v 1.11 2012/11/14 18:14:26 nathan Exp $
; NAME:
;	PLANE_PROJECT
;
; PURPOSE:
;	Deproject SECCHI/LASCO images
;
; CATEGORY:
;	Data analysis, WCS
;
; CALLING SEQUENCE:
;	plane_project, img, coord, new_img
;
; INPUTS:
;	img:		data array to be projected (2-D for one image or 3-D with dimensions <n_images, xsize, ysize>)
;	coords:		helioprojective coordinate map as returned by WCS_GET_COORD: ARR (2, X, Y) 
;
; KEYWORD PARAMETERS:
; Input:
;	INTERP:		to fill gaps in projection image--set as integer number of iterations (1 recommended)
;	SPHER:		whether projection is true spherical or planar
;	XSIZE/YSIZE:	desired pixel dimensions of returned images
;	LIMIT:		cutoff angles for longitude and optionally latitude; anything outside LIMIT is ignored
; Output:
;	SCEN:		FLTARR (2) variable to store pixel location of sun center in the projection image
;	RANGE:		FLTARR (2) variable to store range of projection image in solar radii	
;	PIXEL_MAP:	INTARR (XOUT, YOUT, 2) variable to store map from input image pixels to output image pixels
;   	NEW_COORD   	FLTARR (2, XOUT, YOUT) HPC coordinates of output array
;   	NEW_CDELT   	FLOAT 	Platescale of output array 
;   	NEWCEN=     	FLTARR (2) pixel/IDL coordinate of sun center computed from new_coord
;
; OUTPUTS:
;	new_img:	FLTARR (XOUT, YOUT) projection of img--may have different dimensions as dictated by plate-scale
;
; MODIFICATION HISTORY:
;	Written: Adam Herbst, 2007/06/13
; $Log: plane_project.pro,v $
; Revision 1.11  2012/11/14 18:14:26  nathan
; in debug donot make 2 windows
;
; Revision 1.10  2010/10/21 20:57:05  nathan
; fix bugs
;
; Revision 1.9  2010/10/21 15:49:41  nathan
; add temporary keyword NEWCEN=
;
; Revision 1.8  2010/10/20 21:00:31  nathan
; new method of computing new_cdelt and add debug plot
;
; Revision 1.7  2010/09/29 20:41:25  nathan
; Added NEW_CDELT output
;
; Revision 1.6  2010/09/28 15:08:49  nathan
; update comments
;
; Revision 1.5  2008/06/26 19:34:18  nathan
; update debug behavior call to drawcoordgrid
;
; Revision 1.4  2008/06/25 19:14:13  nathan
; Added /debug features; do not use range var for 2 different things; do not
; undefine radial,azim; radial and azim remain 2-d arrays; new_coord gets
; init value of -180 to indicate missing values
;
; Revision 1.3  2008/02/08 15:17:26  nathan
; fix units of solar radius
;
; Revision 1.2  2008/02/06 18:57:52  nathan
; added DSUN_OBS keyword to plane_project.pro
;
;-

PRO plane_project, img, coords, new_img, NEW_COORD=new_coord, LIMIT=limit, DSUN_OBS=dsun_obs, $
	SCEN=scen, RANGE=range, PIXEL_MAP=pixel_map, XSIZE=xsize, YSIZE=ysize, SPHER=spher, INTERP=interp, $
	DEBUG=debug, NEW_CDELT=new_cdelt, NEWCEN=newcen


IF(~keyword_set(PIXEL_MAP)) THEN BEGIN

;************** case 1: full calculation ******************

IF(~keyword_set(LIMIT)) THEN limit = 100.
IF(n_elements(limit) EQ 1) THEN limit = [limit, 100.]
spher = keyword_set(SPHER)

s = size(img)
IF(s[0] EQ 2) THEN BEGIN
	img = reform(img, 1, s[1], s[2])
	s = size(img)
ENDIF

pixel_map = replicate(-1, s[2], s[3], 2)

c = size(coords)
IF(c[2] NE s[2] OR c[3] NE s[3]) THEN coords = CONGRID(coords, 2, s[2], s[3])

x_ang = reform(coords[0,*,*])
y_ang = reform(coords[1,*,*])

wused = WHERE(x_ang GT -limit[0] AND x_ang LT limit[0] AND y_ang GT -limit[1] AND y_ang LT limit[1])

;stop

used = fltarr(2, n_elements(wused))
used[0,*] = x_ang[wused]
used[1,*] = y_ang[wused]
used = used * !pi/180.

IF(~spher) THEN BEGIN
	x_proj = tan(used[0,*])
	y_proj = sin(used[1,*]) / cos(used[0,*])
ENDIF ELSE BEGIN
print,'% PLANE_PROJECT: spherical'
	;radial = acos(cos(used[0,*])*cos(used[1,*]))
	radial = acos(cos(x_ang* !pi/180.)*cos(y_ang* !pi/180.))
	;azim = atan(sin(used[1,*]), cos(used[1,*])*sin(used[0,*]))
	azim = atan(sin(y_ang* !pi/180.), cos(y_ang* !pi/180.)*sin(x_ang* !pi/180.))
	x_proj = radial * cos(azim)
	y_proj = radial * sin(azim)
	rad = minmax(radial)
	;undefine, radial
	;undefine, azim
ENDELSE

x_ext = minmax(x_proj)
y_ext = minmax(y_proj)

x_diff = x_ext[1]-x_ext[0]
y_diff = y_ext[1]-y_ext[0]

x_proj = x_proj - x_ext[0]
y_proj = y_proj - y_ext[0] 

;resize to give both dimensions same platescale
IF(~keyword_set(XSIZE) OR ~keyword_set(YSIZE)) THEN BEGIN
	xsize = s[2] & ysize = s[3]
	print,'projecting: x-scale / y-scale = ', (xsize/x_diff)/(ysize/y_diff)
	IF (xsize/x_diff GT ysize/y_diff) THEN $
	    ysize = ceil(y_diff * xsize/x_diff) ELSE $
	    xsize = ceil(x_diff * ysize/y_diff)
ENDIF

IF keyword_set(DSUN_OBS) THEN dsunm=dsun_obs ELSE dsunm=1.5e11	; default to 1AU
Dsun=dsunm/6.95508E8	; distance to sun divided by solar radius in meters
;Dsun = 215. ;solar radii
IF (spher) THEN range = x_ext * 180./!pi $ ; range = rad * 180./!pi $ ;
ELSE range = x_ext * Dsun

scen = [-x_ext[0] * xsize / x_diff, -y_ext[0] * ysize / y_diff]

new_img = fltarr(s[1],xsize,ysize)
new_used = intarr(xsize,ysize)
new_coord = make_array(2,xsize,ysize,value=-200d)    ; lt -180 so I can omit unfilled pixels in plot

;stop

FOR i=0L, n_elements(wused)-1 DO BEGIN
	old_x = wused[i] MOD s[2]
	old_y = wused[i] / s[2]
	pix_x = x_proj[wused[i]] * (xsize-1) / x_diff
	pix_y = y_proj[wused[i]] * (ysize-1) / y_diff
	IF(pix_x GE xsize) THEN BEGIN
		print, 'over x'
		pix_x = xsize-1
	ENDIF
	IF(pix_y GE ysize) THEN BEGIN
		print, 'over y'
		pix_y = ysize-1
	ENDIF
	IF(pix_x LT 0) THEN BEGIN
		print, 'under x'
		pix_x = 0
	ENDIF
	IF(pix_y LT 0) THEN BEGIN
		print, 'under y'
		pix_y = 0
	ENDIF
	new_img[*, pix_x, pix_y] = new_img[*, pix_x, pix_y] + img[*, old_x, old_y]
	new_used[pix_x, pix_y] = new_used[pix_x, pix_y] + 1
	pixel_map[old_x, old_y, 0] = pix_x
	pixel_map[old_x, old_y, 1] = pix_y
	new_coord[*, pix_x, pix_y] = coords[*, old_x, old_y]
	IF keyword_set(DEBUG) THEN BEGIN
	    ;print,old_x,pix_x,old_y,pix_y,format='(f8.2,f8.2,f8.2,f8.2)'
	    ;wait,1
	ENDIF
ENDFOR

;
; Get platescale of new_img
;
newn=new_coord
newn[where(new_coord LT -190)]=!values.d_nan
; find y-value of center (IDL coords)
ycens=dblarr(xsize-20)
FOR i=0,xsize-20-1 DO ycens[i]=interp4x(0,newn[1,i+10,*],minsam=20)
ycens=ycens[where(ycens GT 0)]
ycen=avg(ycens)
IF keyword_set(DEBUG) THEN BEGIN
    help,ycens,ycen
    window,1,xsi=xsize,ysi=ysize
    print,'Drawing helioprojective new_coord grid '
    drawcoordgrid,coords=newn,c_charsize=1.5
endif

; now find x-value
xcen=interp4x(0,newn[0,*,ycen],debug=debug,sample=(xsize-100), slope=m)
new_cdelt=m
newcen=[xcen,ycen]

ENDIF ELSE BEGIN

;******************* case 2: use pixel map from previous calculation *****************


s = size(img)
IF(s[0] EQ 2) THEN BEGIN
	img = reform(img, 1, s[1], s[2])
	s = size(img)
ENDIF

new_img = fltarr(s[1],xsize,ysize)
new_used = intarr(xsize,ysize)

FOR a=0,s[2]-1 DO $
	FOR b=0,s[3]-1 DO BEGIN
		pix_x = pixel_map[a,b,0] & pix_y = pixel_map[a,b,1]
		IF(pix_x GE 0) THEN BEGIN
			new_img[*,pix_x,pix_y] = new_img[*,pix_x,pix_y] + img[*,a,b]
			new_used[pix_x, pix_y] = new_used[pix_x, pix_y] + 1
		ENDIF
	ENDFOR

ENDELSE
;******************* end case 2                                      *****************

;******************* interpolate to fill unused pixels in the new image **************
iter = 0
empty = 0
iter_limit = keyword_set(INTERP) ? interp : 0
tmp_new_used = new_used
IF(iter_limit GT 0) THEN message, 'Interpolating to fill projection image', /inf

WHILE(empty[0] GE 0 AND iter LT iter_limit) DO BEGIN
	empty = WHERE(tmp_new_used EQ 0)
	nempty = n_elements(empty)
	FOR i=0L, nempty-1 DO BEGIN
		pix = [empty[i] MOD xsize, empty[i] / xsize]
		xmin = -(1<pix[0]) & xmax = 1<(xsize-1-pix[0]) & ymin = -(1<pix[1]) & ymax = 1<(ysize-1-pix[1])
		FOR a=xmin,xmax DO FOR b=ymin,ymax DO BEGIN
			IF(a EQ 0 AND b EQ 0) THEN CONTINUE
			neighbor = [pix[0] + a, pix[1] + b]
			IF(new_used[neighbor[0], neighbor[1]] GT 0) THEN BEGIN
				new_img[*, pix[0], pix[1]] = new_img[*, pix[0], pix[1]] + new_img[*, neighbor[0], neighbor[1]]
				tmp_new_used[pix[0], pix[1]] = tmp_new_used[pix[0], pix[1]] + new_used[neighbor[0], neighbor[1]]
			ENDIF
		ENDFOR
	ENDFOR
	iter = iter + 1
;	print, 'Replaced ', count, ' of ', nempty, ' points on iteration ', iter
ENDWHILE

used = intarr(s[1], xsize, ysize)
FOR i=0,s[1]-1 DO used[i,*,*] = tmp_new_used > 1
new_img = reform(new_img / used)

IF keyword_set(DEBUG) and datatype(NEW_COORD) NE 'UND' THEN BEGIN
    help,used,x_diff,y_diff,dsun,new_img,img,new_cdelt
    print,xcen,ycen
    print,scen
    ;window,2,xsi=xsize,ysi=ysize
    tv,new_img[0,*,*]
    newrth=hpc2hpr(new_coord,/degr)
    ; newrth is [2,pa,elong]
    print,'Drawing helioprojective radial coordinate grid from hpc2hpr on new image.'
    drawcoordgrid,coords=newrth,c_charsize=1.5
    print,'done'
    wait,2
ENDIF

END

;************************ end current version *********************************

;s=size(img)
;xsize = s[1]
;ysize = s[2]
;
;amount by which IMG is scaled relative to image associated with WCS
;xscale = float(xsize) / wcs.naxis[0]
;yscale = float(ysize) / wcs.naxis[1]
;
;scen = WCS_GET_PIXEL(wcs, [0.,0.]) * [xscale, yscale]
;
;ang = -wcs.roll_angle * !pi/180.
;new_x = ceil(abs(ysize*cos(ang)) + abs(xsize*sin(ang)))
;new_y = ceil(abs(xsize*cos(ang)) + abs(ysize*sin(ang)))
;rot_img = fltarr(new_x, new_y)
;start_x = (new_x - xsize)/2
;start_y = (new_y - ysize)/2
;rot_img[start_x:start_x + xsize-1, start_y:start_y + ysize-1] = img
;rot_img = ROT(rot_img, -wcs.roll_angle)
;
;rotT = [[cos(ang), -sin(ang)], [sin(ang), cos(ang)]] ;CW by ang
;scen = (rotT # (scen - [xsize/2.,ysize/2.]))
;
;xpt = cos(ang) LT 0.
;ypt = sin(ang) LT 0.
;pix = [[xpt*(xsize-1)/xscale, ypt*(ysize-1)/yscale], [(~xpt)*(xsize-1)/xscale, (~ypt)*(ysize-1)/yscale]]
;
;s = size(rot_img)
;xsize = s[1]
;ysize = s[2]
;
;scen = scen + [xsize/2.,ysize/2.]
;
;pixels = fltarr(2,xsize)
;
;FOR i=0, xsize-1 DO BEGIN
;	scl = float(i) / (xsize-1)
;	pixels[*,i] = (1-scl)*pix[*,0] + scl*pix[*,1]
;ENDFOR
;
;coord = WCS_GET_COORD(wcs, pixels)
;
;;stop
;
;IF(~keyword_set(LIMIT)) THEN limit = 40.
;
;start_pix = 0
;end_pix = xsize-1
;WHILE coord[0,start_pix] LE -limit DO start_pix = start_pix + 1
;WHILE coord[0,end_pix] GE limit DO end_pix = end_pix - 1
;start_ang = coord[0,start_pix]
;end_ang = coord[0,end_pix]
;rot_img = rot_img[start_pix:end_pix, *]
;coord = coord[*, start_pix:end_pix]
;scen[0] = scen[0] - (xsize - (end_pix - start_pix + 1))
;xsize = end_pix - start_pix + 1
;
;start_ang = start_ang * !pi/180.
;end_ang = end_ang * !pi/180.
;
;mid_ang = (start_ang + end_ang)/2.
;diff_ang = end_ang - start_ang
;diff_y = 2*tan(diff_ang/2.0)
;start_x = tan(start_ang)
;end_x = tan(end_ang)
;diff_x = end_x - start_x
;
;d_sun = 215.
;range = [d_sun*start_x, d_sun*end_x]
;
;scale = xsize / diff_x
;
;new_img = fltarr(xsize,ysize)
;new_used = intarr(xsize,ysize)
;
;FOR i=0,xsize-1 DO BEGIN
;;	curr_y = (i - xsize/2.) * diff_y / xsize
;;	new_x = tan(atan(curr_y) + mid_ang)
;;	IF(i EQ 1) THEN scale = 1 / (new_x - start_x)
;;	pix_x = CEIL((new_x - start_x) * scale)
;
;	lon = coord[0,i] * !pi/180.
;	pix_x = ceil((tan(lon) - start_x) * scale)
;
;	IF(pix_x LT 0) THEN pix_x = 0
;	IF(pix_x GE xsize) THEN	BEGIN
;		pix_x = xsize-1
;		print, 'Out of bounds at x = ' + STRTRIM(STRING(i),2)
;	ENDIF
;	new_img[pix_x,*] = new_img[pix_x,*] + rot_img[i,*]
;	new_used[pix_x,*] = new_used[pix_x,*] + 1
;ENDFOR
;
;scen[0] = (-start_x * scale)
;
;new_img = new_img / (new_used > 1)
;
;;stop
;
;undefine, new_used
;
;END