; Copyright (c) 1998-2000 A.P. Hitchcock  All rights reserved
;+
;NAME:
;	AX_PROF_RAD
;
;LAST CHANGED: ----------------------------------- 12-sep02
;
;PURPOSE:
; This procedure generates radial profiles from images.
;
;CATEGORY:
;	AXIS: image analysis
;
;CALLING SEQUENCE:
;	AX_PROF_RAD [, axis = axis]
;
;CALLED FROM AXIS:
;	->Images->Profiles->Radial
;
; INPUTS:
;	Uses data in currently displayed buffer
;
; KEYWORDS
;	AXIS - if set, indicates called from AXIS
;
; OUTPUT
; (1) average diameter profile spectrum (put into user selected buffer)
; (2) linearized diametric profile image (put into user selected buffer)
; further analysis of the radial profile can be done using the linescan routines
;
;COMMON BLOCKS:
;	@AXIS_COM	standard set of common blocks
;
; RESTRICTIONS
; Must be called from AXIS
;
;PROCEDURE:
;	The user identifies a diameter across the structure
; of interest, and selects the sampling mode: Linear (1) or Circular (2).
; The routine samples the image as requested
; and generates the avarage diametric profile as well as an angular map (0-180o)
; Optional symmetrization and alignment routines are included.
;
;MODIFICATION HISTORY:
; (26-dec-99 aph) first isolated from AXIS_C (written Dec-99)
; (31-dec-99 aph) AXIS standard documentation
; (23-feb-00 aph) add groupID to call to buffer
; (19-dec-00 cgz) substituted ax_white_color_index with ax_plot_axes_color_index
; (12-sep-02 aph) use correct groupID for buffer !
;-

PRO AX_PROF_RAD, axis = axis
@axis_com
on_error,2

HANDLE_VALUE, Data(CurBuf), tmp
IF n_tags(tmp) NE 0 THEN begin
  IF tmp.t EQ '2d' then begin
	SetGraf, 'MainImg'
   	WIDGET_CONTROL, Uprompt, SET_VALUE='Radial profile: select one edge'
  	CURSOR, x1, y1, /UP, /Data
  	WIDGET_CONTROL, Uprompt, SET_VALUE='Radial profile: select opposite edge'
  	CURSOR, xo, yo, /UP, /Data
  	xc = x1 + (xo-x1)/2.
  	yc = y1 + (yo-y1)/2.
  	print
  	print, 'Radial profile'
  	text = string(format='("Centre: ",f8.2,"  ",f8.2)',xc, yc)
  	print, text
  	centreInd = intarr(2)
  	centreInd = DIndex(xc,yc,tmp)
    radius = sqrt((xo-xc)*(xo-xc) + (yo-yc)*(yo-yc))
    circumf = 6.29*radius
    pixel = abs(tmp.x(centreInd(0) + 1) - tmp.x(centreInd(0)))
    text = string(format="('Radius (um) = ',f6.2, '    circumference (um) = ', f6.2,'    pixel (um) = ',f5.3)", $
     radius, circumf, pixel)
	print, text
    n_points = fix(2*radius/pixel)
    n_profiles = fix(circumf/pixel)
    angle = fix(180/n_profiles)
    if angle EQ 0 then angle = 1  ; force to a minimum of 1 deg.
    angle = get_num(Prompt = 'angle step (deg)',val = angle, group = Axis_ID)
    n_profiles = 180./angle
    print, 'Rad_map: (points, angles) ', n_points, fix(n_profiles)
	widget_control,/hourglass
	rad_map = fltarr(n_profiles, n_points)
	rad_meth = get_num(Prompt='Lines(1), circles(2)',Val = 1,group = axis_ID)
; ------------- diameter lines ------------------------
	if rad_meth EQ 1 then begin
		line = fltarr(n_points)
		angl = 3.14*angle/180.
		for j = 0, n_profiles-1 do begin	; extract diameter profile at angle y(i)
	  		x1 = xc - radius*cos(j*angl)
	  		y1 = yc - radius*sin(j*angl)
	  		x2 = xc + radius*cos(j*angl)
	  		y2 = yc + radius*sin(j*angl)
			dx = (x2 - x1)/n_points
			dy = (y2 - y1)/n_points
			plots,x1,y1,/data
			for i = 0, n_points-1 DO BEGIN
				xt = x1 + i*dx
				yt = y1 + i*dy
				ind = DIndex(xt,yt,tmp)  ;Get index of closest data point
				line(i) = tmp.d(ind(0),ind(1))
			endfor
			plots,x2, y2,/data,/continue, color = ax_plot_axes_color_index		;ax_white_color_index
			rad_map(j,*) = line
		endfor
	endif else begin
	; ---------- circles with diameter sampling
		n_points = fix(n_points/2)
		n_profiles = fix(2*n_profiles)
		rad_map = fltarr(n_profiles, n_points)
		plots,xc,yc,/data
		ind = DIndex(xc,yc,tmp)
		for i = 0,n_profiles-1 do rad_map(i,0) = tmp.d(ind(0),ind(1))
		for i = 1, n_points-1 do begin
			rt = i*radius/(n_points-1)
			np_circ = fix(6.29*rt/pixel)
			if np_circ GT n_profiles then np_circ = n_profiles
			lin_tmp = fltarr(np_circ)
			x_tmp = findgen(np_circ)*360./np_circ  ; angular pos in deg.
			SetGraf, 'MainImg'
			for j = 0,np_circ-1 do begin
				ang = 3.14*x_tmp(j)/180.
				xt = xc + rt*cos(ang)
				yt = yc + rt*sin(ang)
				ind = DIndex(xt,yt,tmp)
				lin_tmp(j) = tmp.d(ind(0),ind(1))
				plots, xt, yt,/data,/continue, color = ax_plot_axes_color_index	;ax_white_color_index
			endfor
			rad_map(*,i) = interpol(lin_tmp, n_profiles)
		endfor
; -------- convert from radial to diameter display
		n_profiles = n_profiles/2
		n_points = n_points*2
		t=fltarr(n_profiles,n_points)
		for i = 0, n_points/2-1 do $
		t(*,n_points/2-1-i) = rad_map(n_profiles:2*n_profiles-1,i)
		for i = 0, n_points/2-1 do $
			t(*,n_points/2+i) = rad_map(0:n_profiles-1,i)
		rad_map = t
	endelse

; --------------- symmetrize diameters ---------
	test = get_num(Prompt='Symmetrize diamters ? (0 = No, 1 = Yes)',Val = 1,group = axis_ID)
	if test NE 0 then begin
		twidth = intarr(n_profiles)
		for i = 0, n_profiles-1 do begin
			t = reform(rad_map(i,*))
			tmin = min(t,max=tmax)
			tmid = (tmax - tmin)/2 + tmin
			if t(n_points/2) LT tmid then begin
				tobj = where(t LT tmid, count)	; signal is LOWER inside region
			endif else begin
				tobj = where(t GT tmid, count)	; signal is HIGHER inside region
			endelse
			tshift = 0
			if count NE 0 then begin
				nlow = tobj(0)
				nhi  = tobj(n_elements(tobj)-1)
				nmid = fix((nhi - nlow )/2) + nlow
				twidth(i) = nhi - nlow
				tshift = nmid - n_points/2
	;			print, n_points, nlow,nmid, nhi
			endif
			t = shift(t,-1*tshift)
			rad_map(i,*) = t
		endfor
; ------------ optional force to constant width ----
		test = get_num(Prompt='Force constant width ? (0 = No, 1 = Yes)',Val = 1,group = axis_ID)
		if test NE 0 then begin
			width = min(twidth(where(twidth GT 0)))
			print, 'forcing diameter to ', pixel*width
			WIDGET_CONTROL, /hourglass
			xtmp = pixel*findgen(n_points) - radius
			for i = 0, n_profiles-1 do begin
				t = reform(rad_map(i,*))
				shrink = (float(width)/float(twidth(i)))
	;			print, 'line ',i,' shrunk by ', shrink
				xt = shrink*xtmp
				tt = interpol(t,xt,xtmp)
				rad_map(i,*) = tt
			endfor
		endif
	endif

; ----------- generate output
	y = pixel*findgen(n_points) - radius
	x = angle*findgen(n_profiles)
	xl = 'Angle (o)'
	yl = 'Distance (um)'
	WIDGET_CONTROL, Uprompt, SET_VALUE='Buffer for average profile'
 	CurBuf = buffer(group = axis_ID)	; generate and save the diameter average
 	if CurBuf eq -1 then return
 	d = fltarr(n_points)
 	for i = 0, n_points-1 do begin
 		d(i) = median(rad_map(*,i))
 	endfor
	s = {t:'1d', x:y, d:d, xl:yl, dn:d, dl: '<D> ' + tmp.dl}
   	HANDLE_VALUE, Data(CurBuf), s, /set
	Label(CurBuf) = s.dl
	PlotBuf,CurBuf
; 	CurBuf = 9	; generate and save the radial average
; 	r = fltarr(fix(n_points/2))
; 	x_offset = n_points/2
; 	for i = x_offset, n_points-1 do begin
; 		r(i) = (d(i) + d(i-x_offset))/2.
; 	endfor
;	s = {t:'1d', x:x, d:d, xl:xl, dn:d, dl: '<R> ' + tmp.dl}
;  	HANDLE_VALUE, Data(CurBuf), s, /set
;	Label(CurBuf) = s.dl
;	PlotBuf,CurBuf
	WIDGET_CONTROL, Uprompt, SET_VALUE='Buffer for radial profile map'
 	CurBuf = buffer(group = axis_ID)	; save and display the radial map
 	if CurBuf eq -1 then return
	s = {t:'2d', x:x, y:y, d:rad_map, xl:xl, yl:yl, dl: 'Rmap ' + tmp.dl}
    HANDLE_VALUE, Data(CurBuf), s, /set
	Label(CurBuf) = s.dl
	PlotBuf,CurBuf
  endif else WIDGET_CONTROL, Uprompt, SET_VALUE='Profiles: only for images'
ENDIF

END
