pro spheretouch

;given a GIF image that is black (RGB=0,0,0) and white (RGB=255,255,255)
;get the perimeter by advancing in from the four sides and stopping at first
;black contact, convert to polar co-ordinates and remove all duplicate points
;
;click on two extremes of a calibration dip,
;sweep through 360 degrees, determining maximum, next point of contact for line
;average distance between points of contact
;all in meters x 10^-6 (microns)

;David Marshall 15, 06, 2000

choice = 'c:/temp/smooth2edited'   ;filename containing sphere profile
PRINT, 'Input GIF file name:  (c:/temp/smooth2edited)'
READ, choice
IF choice eq '' THEN choice = 'c:/temp/smooth2edited'
fileName=choice
statFileName=choice + '.dat'
gifFileName=choice + '.gif'
individualFileCheck=FINDFILE(gifFileName, count=found)

IF found eq 0 THEN BEGIN
	print, 'error: ', gifFileName,' file not found'
ENDIF ELSE BEGIN
	choice='1500'
	PRINT, 'Input approximate sphere radius in microns. [1500]'
	READ, choice
	IF choice eq '' THEN choice = '1500'
	realRadius=FLOAT(choice)
	print, 'Reading from ', gifFileName
	READ_GIF, gifFileName, bigImage, R, G, B ;image to be read

	;bigImage=congrid(bigImage, 256,256) ;decrease size so computation time reduced

	bigImageDimension=SIZE(bigImage)
	bigImageX=bigImageDimension[1] ;x dimension
	bigImageY=bigImageDimension[2] ;y dimension

	TVLCT, R, G, B ;display image
	SLIDE_IMAGE, bigImage, TITLE='Original Talyrond GIF image', SLIDE_WINDOW=slideWindowNumber, TOP_ID=base;, XSIZE=bigImageX, YSIZE=bigImageY
	WSET, slideWindowNumber
	print, 'Left mouse click at one end of calibration dip'
	CURSOR, x1, y1, /DEVICE, /DOWN
	print, 'Left mouse click at the other end of calibration dip'
	CURSOR, x2, y2, /DEVICE, /DOWN
	PLOTS,[x1,x2], [y1,y2], /DEVICE, COLOR=0
	WIDGET_CONTROL, /DESTROY, base

	calDip=SQRT((x2-x1)^2+(y2-y1)^2) ;this is the 11.7 micron calibration dip in pixels
	scaleMicronperPixel=11.7/calDip ;calDip number of pixels = 11.7e-6 meters
	pixelRadius=realRadius/scaleMicronperPixel ;full radius of sphere in pixel units

	print, 'Finding perimeter'
	black=WHERE(bigImage eq 0) ;index all the black "ink" points, that is where =0

	perimeter=LONARR(2,N_ELEMENTS(black));this is way too big but perimeter will never have more points than this
	nextPoint=0 ;perimeter point index
	columnWeight=LONARR(bigImageX) ;the weight of each column (also length since each pixel weight is 1)
	rowWeight=LONARR(bigImageY) ;the weight of each row

	FOR x=0,bigImageX-1 DO BEGIN
		columnY=WHERE(bigImage[x,*] eq 0, count)
		IF count ne 0 THEN BEGIN ;something in this column
			columnWeight[x]=columnY[N_ELEMENTS(columnY)-1]-columnY[0]+1 ;the first and last Y position of a 0 in this column
			perimeter[0,nextPoint]=x ;record first and last position in this column
			perimeter[1,nextPoint]=columnY[0]
			perimeter[0,nextPoint+1]=x
			perimeter[1,nextPoint+1]=columnY[N_ELEMENTS(columnY)-1]
			nextPoint=nextPoint+2
		ENDIF
	ENDFOR
	cmX=ROUND(TOTAL(columnWeight * INDGEN(bigImageX))/TOTAL(columnWeight)) ;cmX=[sum of (M*X)]/(total mass)

	FOR y=0,bigImageY-1 DO BEGIN
		rowX=WHERE(bigImage[*,y] eq 0, count)
		IF count ne 0 THEN BEGIN
			rowWeight[y]=rowX[N_ELEMENTS(rowX)-1]-rowX[0]+1 ;the first and last X position of a 0 in this row
			perimeter[0,nextPoint]=rowX[0] ;record first and last position in this row
			perimeter[1,nextPoint]=y
			perimeter[0,nextPoint+1]=rowX[N_ELEMENTS(rowX)-1]
			perimeter[1,nextPoint+1]=y
			nextPoint=nextPoint+2
		ENDIF
	ENDFOR
	cmY=ROUND(TOTAL(rowWeight * INDGEN(bigImageY))/TOTAL(rowWeight)) ;cmY=[sum of (M*Y)]/(total mass)

	perimeter=perimeter[*,0:nextPoint-1] ;clip it to actual number of points
	perimeter[0,*]=perimeter[0,*]-cmX ;transform co-ordinates to be centred about CM
	perimeter[1,*]=perimeter[1,*]-cmY

	print,'Converting to polar co-ordinates'
	polarCoord = CV_COORD(FROM_RECT=perimeter, /TO_POLAR) ;convert to polar co-ords

	culledPolar = polarCoord(*,UNIQ(polarCoord[0,*], SORT(polarCoord[0,*]))) ;remove duplicate points


	statData=fltarr(4,N_ELEMENTS(culledPolar)) ;will hold stats
	maxLocalPixelRadius=MAX(culledPolar[1,*], maxRadiusSubScript, MIN=minLocalPixelRadius) ;need this widest contact range calc
	averageLocalPixelRadius=MEAN(culledPolar[1,*]) ;average for whole sphere based on Talyrond data
	maxPixelRadius=pixelRadius-averageLocalPixelRadius+maxLocalPixelRadius ;radius of the maximum in pixel units
	minPixelRadius=pixelRadius-averageLocalPixelRadius+minLocalPixelRadius ;radius of the minimum in pixel units
	angleOfWidestContact=ACOS(minPixelRadius/maxPixelRadius) ;gives maximum possible range to check

	LOADCT, 0
	WINDOW, /FREE, TITLE='Perimeter of Talyrond trace', XSIZE=400, YSIZE=400 ;, XSIZE=bigImageX, YSIZE=bigImageY
	PLOT, /POLAR, /ISOTROPIC, COLOR=255, /DEVICE, culledPolar[1,*], culledPolar[0,*], PSYM=3 ;do not connect points
	fullWindowNum=!D.WINDOW

;	get an initial starting point - the first angle/radius pair in data
	angleFirstContactPoint=culledPolar[0,0] ;start at the beginning of angle/radius data(-PI)
	maxLocalPixelRadius=culledPolar[1,0] ;radius of first point

	index=0
	REPEAT BEGIN
		maxPixelRadius=double(pixelRadius-averageLocalPixelRadius+maxLocalPixelRadius) ;radius of the first contact point in pixel units
		angleTangentContactPoint=angleFirstContactPoint+angleOfWidestContact ;angle of furthest point (on horizon is very smooth)
		possibleContactRange=double(culledPolar[*,WHERE((culledPolar[0,*] gt angleFirstContactPoint) AND (culledPolar[0,*] lt angleTangentContactPoint))]) ;array to of possible points to look at

		pt2ptLength=SQRT(maxPixelRadius^2+(pixelRadius-averageLocalPixelRadius+possibleContactRange[1,*])^2- $
			2*maxPixelRadius*(pixelRadius-averageLocalPixelRadius+possibleContactRange[1,*])*COS(possibleContactRange[0,*]-angleFirstContactPoint));cosine law - length from point to point
		sinBeta=(pixelRadius-averageLocalPixelRadius+possibleContactRange[1,*])*SIN(possibleContactRange[0,*]-angleFirstContactPoint)/pt2ptLength ;sine of beta
		sinBeta=sinBeta < 1.0 ;sometimes this is > 1.0 (due to round off?) clamp to 1.0
		beta=ASIN(sinBeta)
		sinCorrection=WHERE((pixelRadius-averageLocalPixelRadius+possibleContactRange[1,*])^2-pt2ptLength^2-maxPixelRadius^2 gt 0.0, count) ;assign proper angle if beta >90 degrees
		IF count gt 0 THEN beta[sinCorrection]=!PI-beta[sinCorrection] ;the ones to correct
		maxBeta=MAX(beta, secondContactSubScript) ;find largest beta thereby finding next peak

		angleSecondContactPoint=possibleContactRange[0,secondContactSubScript] ;this is where second contact point is
		secondLocalContactRadius=possibleContactRange[1,secondContactSubScript] ;its radius
		secondContactRadius=pixelRadius-averageLocalPixelRadius+secondLocalContactRadius
		angleBetweenContactPoints=angleSecondContactPoint-angleFirstContactPoint
		actualContactRangeSubScripts=WHERE((culledPolar[0,*] ge angleFirstContactPoint) AND (culledPolar[0,*] le angleSecondContactPoint));pull out the data in the range of possible contact
		actualContactRange=double(culledPolar[*,actualContactRangeSubScripts])
		bufferedActualContactRangeSubScriptsStart=(ActualContactRangeSubScripts[0]-2) > 0 ;for investigative plotting - includes previous and following 2 points
		bufferedActualContactRangeSubScriptsEnd=(ActualContactRangeSubScripts[N_ELEMENTS(ActualContactRangeSubScripts)-1]+2) < N_ELEMENTS(culledPolar[0,*]) ;for plotting purposes

;shiftedAngles=actualContactRange
;shiftedAngles[0,*]=shiftedAngles[0,*]-(angleFirstContactPoint+angleSecondContactPoint)/2+!PI/2
;shiftedAngles[1,*]=shiftedAngles[1,*]+pixelRadius-averageLocalPixelRadius
;cartCoord = CV_COORD(FROM_POLAR=shiftedAngles, /TO_RECT) ;convert to cart co-ords
;cartCoord=cartCoord[*,SORT(cartCoord[0,*])] ;sort on X
;centerX=(cartCoord[0,N_ELEMENTS(cartCoord[0,*])-1]+cartCoord[0,0])/2
;centerY=cartCoord[1,0]+SQRT(pixelRadius^2-((cartCoord[0,N_ELEMENTS(cartCoord[0,*])-1]-cartCoord[0,0])/2)^2)
;otherSphere=cartCoord
;otherSphere[1,*]=centerY-SQRT(pixelRadius^2-otherSphere[0,*]^2)
;separation=otherSphere[1,*]-cartCoord[1,*]
;HunchaverageContactHeight=MEAN(separation)

		localIntegration=INT_TABULATED(actualContactRange[0,*], actualContactRange[1,*])
		boxIntegration=INT_TABULATED([angleFirstContactPoint, angleSecondContactPoint],[maxLocalPixelRadius, secondLocalContactRadius])
		HunchaverageContactHeight=(boxIntegration-localIntegration)/angleBetweenContactPoints

;		averageContactHeight=(maxLocalPixelRadius+secondLocalContactRadius)/2-MEAN(actualContactRange[1,*])
;		minContactHeight=(maxLocalPixelRadius < secondLocalContactRadius)-MEAN(actualContactRange[1,*]) ;lesser of the two

		statData[0,index]=angleFirstContactPoint*!RADEG
		statData[1,index]=angleBetweenContactPoints*!RADEG
		statData[2,index]=maxBeta*!RADEG
		statData[3,index]=HunchaverageContactHeight*scaleMicronperPixel

		WSHOW, fullWindowNum, 0
		WSET, fullWindowNum
		PLOT, /POLAR, /ISOTROPIC, COLOR=255, /DEVICE, culledPolar[1,*], culledPolar[0,*], /NODATA, /NOERASE ;do not connect points
		OPLOT, /POLAR, [maxLocalPixelRadius, secondLocalContactRadius], [angleFirstContactPoint, angleSecondContactPoint], COLOR=255, PSYM=-2
		WSHOW, fullWindowNum, 1

		WINDOW, /FREE, TITLE='contact range'
		contactWindowNum=!D.WINDOW
		PLOT, /POLAR, pixelRadius-averageLocalPixelRadius+culledPolar[1,bufferedActualContactRangeSubScriptsStart:bufferedActualContactRangeSubScriptsEnd], $
			culledPolar[0,bufferedActualContactRangeSubScriptsStart:bufferedActualContactRangeSubScriptsEnd],/YNOZERO
		OPLOT, /POLAR, [maxPixelRadius, secondContactRadius], [angleFirstContactPoint, angleSecondContactPoint], COLOR=255, PSYM=-2
		PRINT,angleFirstContactPoint*!RADEG, angleBetweenContactPoints*!RADEG, maxBeta*!RADEG, $
			secondContactSubScript, HunchaverageContactHeight*scaleMicronperPixel

		index=index+1

;		choice = 'string'
;		PRINT, 'hit enter to continue, (q to quit)'
;		READ, choice
		angleFirstContactPoint=angleSecondContactPoint ;move on - second contact point becomes first now
		maxLocalPixelRadius=secondLocalContactRadius
;		IF STRMID(choice,0,1) eq 'q' THEN angleFirstContactPoint=!PI;get out of loop early
		WDELETE, contactWindowNum

	ENDREP UNTIL angleFirstContactPoint+angleOfWidestContact ge !PI ;we're back at the starting point


	statData=statData[*,0:index-1]
	meanStatData=fltarr(3)
	meanStatData[0]=MEAN(statData[1,*]) ;angleBetweenContactPoints
	meanStatData[1]=MEAN(statData[2,*]) ;maxBeta*!RADEG
	meanStatData[2]=MEAN(statData[3,*]) ;Hunch average gap separation

	print, '   angle start      between     beta        sphr sep'
	FOR index=0,N_ELEMENTS(StatData[0,*])-1 DO BEGIN
		print, statData[*,index]
	ENDFOR
	print, '   angle start      between     beta        sphr sep'
	print, 'averages - - - -',meanStatData

	choice = 'yes'
	PRINT, 'Would you like to print these data to file ',statFileName, '?    [YES]'
	READ, choice
	IF choice eq '' THEN choice = 'yes'
	choice=STRLOWCASE(STRMID(choice,0,1)) ;grab first letter
	IF choice eq 'y' THEN BEGIN
		PRINT,'Writing to file ', statFileName
		OPENW, unit, statFileName, /GET_LUN
		PRINTF, unit,  FORMAT='("averages - - - - - - ",4(F10.4,2x))',meanStatData
		PRINTF, unit, '   angle start      between     beta        sphr sep'
		FOR index=0,N_ELEMENTS(StatData[0,*])-1 DO BEGIN
			PRINTF, unit, FORMAT='(4(F10.4,2x))', statData[*,index]
		ENDFOR
		CLOSE,unit
		FREE_LUN,unit
		print, 'Done writing'
	ENDIF

ENDELSE

print, 'Program finished'
END