#pragma rtGlobals=3		// Use modern global access method and strict wave access.
Menu "Macros"
	"Load tif movie", LoadMoviePanel()
End

Function LoadMoviePanel()
	initialization()
	newPanel/w = (1170,50,1350,135)/n=Panel
	
	Button LoadMovie title="Load",pos={5,10},size={70,30},fColor=(0,52224,0),proc=LoadMovie
	Button Kill title="Kill",pos={5,45},size={70,30},fColor=(20000,0,20000),proc=KillData
	EnableDisableAnalysisButton()
end
	
Function initialization()
	newdatafolder/O root:WindowsExplorer
	newdatafolder/O root:RawData
	newdatafolder/O root:Analysis
	setdatafolder root:WindowsExplorer

end

Function LoadMovie(CtrlName):ButtonControl
	string CtrlName
	initialization()
	setdatafolder root:WindowsExplorer
	string/G fileName
	string/G filePath
	string/G format = ".tif"
	string/G file_Name
	WindowsExplorer()
	Setdatafolder root:RawData
	variable/G sizeOfRows,sizeOfColumns,sizeOfFrames
	variable/G MinColorValue,MaxColorValue
	variable/G frameNum
	variable/g FirstFrame=0
	variable/g LastFrame
	variable/G framesPerSec=4
	wave singleFrame
	wave movie
	
	ImageLoad/C=-1/G/O/Q/T=tiff filePath
	file_Name = s_filename
	duplicate/O/R=[][][frameNum,frameNum] $s_filename, singleFrame
	sizeOfRows = dimsize($s_filename,0)
	sizeOfColumns = dimsize($s_filename,1)
	sizeOfFrames = dimsize($s_filename,2)
	LastFrame = sizeOfFrames-1
	rename $s_filename,movie
	EnableDisableAnalysisButton()
	
	doWindow/R/f ShowMovie
	display/w=(10,10,360,200*1.618+10) as s_filename
	newImage/F/HOST=# singleFrame
	SetAxis/A left

	wavestats/Q/M=1 root:rawData:movie
	MinColorValue=V_min
	MaxColorValue=V_max
	ModifyImage singleFrame ctab= {MinColorValue,MaxColorValue,grays,0}
	movesubwindow fnum=(50,0,200,200*1.618)
	newPanel/HOST=##/FG=(FL,FT,,FB)/W=(0,0,50,200*1.618)
	ModifyPanel frameInset=0,frameStyle=0
	
	doupdate
	button play appearance={default},fstyle=0,size={50,30},pos={5,5},title="Play",fColor=(0,30000,0),proc=Play_Movie
	setVariable FramesPerSec,appearance={default},pos={5,35},title=" ",value=framesPerSec,limits={1,100,0},size={20,18} 
	setDrawEnv fname= "Times New Roman",fsize=10
	drawText 30,48,"Frame/s"
	
	ValDisplay FrameNum_shown,pos={22,380},size={40,13}
	ValDisplay FrameNum_shown,valueBackColor=(47872,47872,47872)
	ValDisplay FrameNum_shown,limits={0,0,0},barmisc={0,1000}
	ValDisplay FrameNum_shown,value= #"root:RawData:frameNum"
	
	slider frameNum appearance={os9},limits={0,sizeOfFrames-1,1},ticks=1,live=1,fsize=9,pos={22,55},Size={30,1.618*200},side=1,vert=1,variable= frameNum,proc=MovingSlider
	setVariable FirstFrame,appearance={default},pos={5,400},title=" ",value=FirstFrame,limits={0,sizeOfFrames-1,0},size={22,13} 
	setVariable LastFrame,appearance={default},pos={5,415},title=" ",value=LastFrame,limits={FirstFrame,sizeOfFrames-1,0},size={22,13} 
	button cutMovie appearance={default},fstyle=0,size={30,30},pos={30,400},title="Cut",fColor=(15000,15000,0),proc=cutMovie
	
	Dokymograph()
end

Function EnableDisableAnalysisButton()
	wave movie = root:RawData:movie
	if (waveexists(movie) == 0)
		Button Analysis title="Analyze", pos={80,10},size={50,64},fColor=(32768,54528,65280),proc=MovieAnalysis,disable=2
		Button Clear title="Clear",pos={132,10},size={38,64},fColor=(65280,16384,35840),proc=ClearData,disable=2
	else 
		Button Analysis title="Analyze", pos={80,10},size={50,64},fColor=(32768,54528,65280),proc=MovieAnalysis,disable=0
		Button Clear title="Clear",pos={132,10},size={38,64},fColor=(65280,16384,35840),proc=ClearData,disable=0
	endif
end

Function MovingSlider(CtrlName, value, event)
	string CtrlName
	variable value,event
	NVAR frameNum=root:rawData:frameNum
	wave Movie=root:rawData:Movie
	wave singleFrame=root:rawData:singleFrame
	duplicate/O/R=[][][frameNum,frameNum] movie, singleFrame
end

Function Play_Movie(CtrlName):ButtonControl
	string CtrlName
	variable value,event,i,ElapseTime
	NVAR framesPerSec = root:rawData:framesPerSec
	NVAR sizeOfFrames = root:rawData:sizeOfFrames
	NVAR frameNum = root:rawData:frameNum
	ElapseTime = 1/framesPerSec
	if(frameNum==sizeOfFrames-1)
		for(i=0;i<sizeOfFrames;i=i+1)
			frameNum = i
			MovingSlider(CtrlName, value, event)
			Doupdate/E=1
			sleep/s/c=-1 ElapseTime
			if(PressEscapeToAbort())
				break
				beep
			endif
		endfor
	else
		for(i=frameNum;i<sizeOfFrames;i=i+1)
			frameNum = i
			MovingSlider(CtrlName, value, event)
			Doupdate/E=1
			sleep/s/c=-1 ElapseTime
			if(PressEscapeToAbort())
				break
				beep
			endif
		endfor
	endif
end

Function cutMovie(CtrlName):ButtonControl
	string CtrlName
	wave movie = root:rawData:movie
	wave Kymograph = root:rawData:Kymograph
	NVAR FirstFrame = root:rawData:FirstFrame
	NVAR LastFrame = root:rawData:LastFrame
	NVAR sizeOfFrames = root:rawData:sizeOfFrames
	NVAR sizeOfRows = root:rawData:sizeOfRows
	NVAR sizeOfColumns = root:rawData:sizeOfColumns
	NVAR MinColorValue = root:rawData:MinColorValue
	NVAR MaxColorValue = root:rawData:MaxColorValue
	SVAR file_Name = root:WindowsExplorer:file_Name
	wave singleFrame = root:rawData:singleFrame
	variable i,j,k
	
	sizeOfFrames = LastFrame-FirstFrame+1
	make/R/N=(sizeOfRows,sizeOfColumns,sizeOfFrames) movie_temp
	make/R/N=(sizeOfFrames,sizeOfColumns) Kymograph_temp
	
	
	for(j=0;j<sizeOfColumns;j=j+1)
		for(k=0;k<sizeOfFrames;k=k+1)
			Kymograph_temp[k][j] = Kymograph[FirstFrame+k][j]
			for(i=0;i<sizeOfRows;i=i+1)
				movie_temp[i][j][k] = movie[i][j][FirstFrame+k]
			endfor
		endfor
	endfor
	
	//duplicate/O/R=[][][FirstFrame,LastFrame] movie,movie_temp
	//duplicate/O/R=(FirstFrame,LastFrame)[] Kymograph,Kymograph_temp
	killwaves/z movie
	killwindows(0)
	killwaves/z Kymograph
	rename movie_temp movie
	rename Kymograph_temp, Kymograph
	FirstFrame=0
	LastFrame=sizeOfFrames-1
	
	doWindow/R/f ShowMovie
	display/w=(10,10,360,200*1.618+10) as file_Name
	newImage/F/HOST=# singleFrame
	SetAxis/A left

	wavestats/Q/M=1 root:rawData:movie
	MinColorValue=V_min
	MaxColorValue=V_max
	ModifyImage singleFrame ctab= {MinColorValue,MaxColorValue,grays,0}
	movesubwindow fnum=(50,0,200,200*1.618)
	newPanel/HOST=##/FG=(FL,FT,,FB)/W=(0,0,50,200*1.618)
	ModifyPanel frameInset=0,frameStyle=0
	
	doupdate
	button play appearance={default},fstyle=0,size={50,30},pos={5,5},title="Play",fColor=(0,30000,0),proc=Play_Movie
	setVariable FramesPerSec,appearance={default},pos={5,35},title=" ",value=framesPerSec,limits={1,100,0},size={20,18} 
	setDrawEnv fname= "Times New Roman",fsize=10
	drawText 30,48,"Frame/s"
	
	ValDisplay FrameNum_shown,pos={22,380},size={40,13}
	ValDisplay FrameNum_shown,valueBackColor=(47872,47872,47872)
	ValDisplay FrameNum_shown,limits={0,0,0},barmisc={0,1000}
	ValDisplay FrameNum_shown,value= #"root:RawData:frameNum"
	
	slider frameNum appearance={os9},limits={0,sizeOfFrames-1,1},ticks=1,live=1,fsize=9,pos={22,55},Size={30,1.618*200},side=1,vert=1,variable= frameNum,proc=MovingSlider
	setVariable FirstFrame,appearance={default},pos={5,400},title=" ",value=FirstFrame,limits={0,sizeOfFrames-1,0},size={22,13} 
	setVariable LastFrame,appearance={default},pos={5,415},title=" ",value=LastFrame,limits={FirstFrame,sizeOfFrames-1,0},size={22,13} 
	button cutMovie appearance={default},fstyle=0,size={30,30},pos={30,400},title="Cut",fColor=(15000,15000,0),proc=cutMovie
	
	newImage/F/HOST=## root:rawData:Kymograph
	movesubwindow fnum=(200,0,370,10+200*1.618)
	showinfo
	
end

Function Dokymograph()
	wave movie = root:rawData:movie
	NVAR sizeOfFrames = root:rawData:sizeOfFrames
	NVAR sizeOfRows = root:rawData:sizeOfRows
	NVAR sizeOfColumns = root:rawData:sizeOfColumns
	setdatafolder root:rawData
	variable i,j,k
	make/O/N=(sizeOfFrames,sizeOfColumns) Kymograph
	wave Kymograph = root:rawData:Kymograph
	for (k=0;k<sizeOfFrames;k=k+1)
		for (j=0;j<sizeOfRows;j=j+1)
			for (i=0;i<sizeOfColumns;i=i+1)
				Kymograph[k][i] = Kymograph[k][i]+movie[j][i][k]
			endfor
		endfor
	endfor
	
	newImage/F/HOST=## root:rawData:Kymograph
	movesubwindow fnum=(200,0,370,10+200*1.618)
	showinfo
end

Function MovieAnalysis(CtrlName):ButtonControl
	string CtrlName
	wave Kymograph = root:rawData:Kymograph
//	wave Kymo_upper_zero = root:Analysis:Kymo_upper_zero
//	wave Kymo_lower_zero = root:Analysis:Kymo_lower_zero
//	wave Kymo_upper_temp, Kymo_lower_temp
	NVAR sizeOfFrames = root:rawData:sizeOfFrames
	NVAR sizeOfColumns = root:rawData:sizeOfColumns
//	NVAR bleachedIntensity = root:Analysis:bleachedIntensity
	setdatafolder root:Analysis
	variable leftBoundary,rightBoundary //set by cursor position to describe photobleached area
	make/N=(sizeOfFrames)/O Boundary
	make/N=(sizeOfFrames) BleachedIntensity
	variable i,j
//	variable Upper_sizeOfColumns,Lower_sizeOfColumns
	variable TimePerFrame = 120
	variable/G pixelSize
	variable pixelSize_local = 138.6
	variable BleachedWidth = 2
	
	//string graphName = ChildWindowList("graph0")
	
	
	setActiveSubwindow Graph0#G1
	variable XA = pcsr(A,"Graph0#G1")
	variable YA = qcsr(A,"Graph0#G1")
	variable XB = pcsr(B,"Graph0#G1")
	variable YB = qcsr(B,"Graph0#G1")
	if (XA==0 && XB==sizeOfFrames-1)
		leftBoundary = YA
		rightBoundary = YB
	elseif(XA==sizeOfFrames-1 && XB==0)
		leftBoundary = YB
		rightBoundary = YA
	else
		doalert 0,"Your cursor is not at right position"
	endif
	
	for (i=0;i<sizeOfFrames;i=i+1)
		Boundary[i]=ceil(leftBoundary + i*(rightBoundary - leftBoundary)/sizeOfFrames)
		for (j=Boundary[i]-BleachedWidth;j<=Boundary[i]+BleachedWidth;j=j+1)
			BleachedIntensity[i]=BleachedIntensity[i]+Kymograph[i][j]
		endfor
		BleachedIntensity[i] = BleachedIntensity[i]/(2*BleachedWidth+1) 
	endfor
	
	if (leftBoundary > rightBoundary)
		duplicate/O/R=[0,sizeOfFrames-1][0,leftBoundary] Kymograph, Kymograph_lower
		duplicate/O/R=[0,sizeOfFrames-1][rightBoundary,sizeOfColumns-1] Kymograph, Kymograph_upper
		
		for(i=0;i<sizeOfFrames;i=i+1)
			for (j=0;j<=leftBoundary;j=j+1)
				if (j>Boundary[i])
					Kymograph_lower[i][j]=bleachedIntensity[i]
				endif
			endfor
			
			for (j=0;j<sizeOfColumns-rightBoundary;j=j+1)
				if(j<Boundary[i]-rightBoundary)
					Kymograph_upper[i][j]=bleachedIntensity[i]
				endif
			endfor
		endfor
	else
		duplicate/O/R=[0,sizeOfFrames-1][0,rightBoundary] Kymograph, Kymograph_lower
		duplicate/O/R=[0,sizeOfFrames-1][leftBoundary,sizeOfColumns-1] Kymograph, Kymograph_upper
		for(i=0;i<sizeOfFrames;i=i+1)
			for (j=0;j<=rightBoundary;j=j+1)
				if (j>Boundary[i])
					Kymograph_lower[i][j]=bleachedIntensity[i]
				endif
			endfor
			
			for (j=0;j<sizeOfColumns-leftBoundary;j=j+1)
				if(j<Boundary[i]-leftBoundary)
					Kymograph_upper[i][j]=bleachedIntensity[i]
				endif
			endfor
		endfor
	endif
	
	Prompt timePerFrame, "Enter time per frame (s)"
	Prompt pixelSize_local,"Enter size of pixel (nm)"
	DoPrompt "More info is needed",timePerFrame,pixelSize_local
	pixelSize = pixelSize_local
	
	make/O/N=(sizeOfFrames) timeSeries
	for (i=0;i<sizeOfFrames;i=i+1)
		timeSeries[i] = i * TimePerFrame
	endfor
	
	display/w=(370,10,520,10+95*1.618) as "Upper part"
	newimage/f/HOST=# Kymograph_upper
	setAxis/A=2 bottom
	movesubwindow fnum=(0,0,150,95*1.618)
	

	
//	make/N=(sizeOfFrames) Upper_centralShift
//	for (i=0;i<sizeOfFrames;i=i+1)
//		CurveFit/Q gauss, kymograph_upper[i][*]
//		rename W_coef wavetemp
//		wave wavetemp
//		Upper_centralShift[i] = round(wavetemp[2])
//		killwaves/z wavetemp
//	endfor
	 
//	//calculate cross correlation
// 	variable xm,ym,sx,sy,sxy,d
////	variable halfseriesSize_upper = 30
////	variable maxDelay_upper = 30
//	variable halfseriesSize_lower = 10
//	variable maxDelay_lower = 10


	CrossCorrelation("Kymograph_upper","_upper",0)
	wave centralshift_upper = root:Analysis:centralshift_upper
	
	CrossCorrelation("Kymograph_upper","_upperEdge",1)
	wave centralshift_upperEdge = root:Analysis:centralshift_upperEdge
	
//	Upper_sizeOfColumns = dimsize(Kymograph_upper,1)
//	make/N=(2*maxDelay_upper+1,sizeOfFrames) Upper_CrossCor
//	
//	make/N=(2*maxDelay_upper+1) upper_delayNum
//	for(i=0;i<2*maxDelay_upper+1;i=i+1)
//		upper_delayNum[i] = i - maxDelay_upper
//	endfor
//	
//	make/N=(sizeOfFrames) Upper_centralShift_temp
//	variable upper_0 = Upper_centralShift[0]-(sizeOfColumns - upper_sizeOfColumns) - halfseriesSize_upper
//	variable upper_1 = Upper_centralShift[0]-(sizeOfColumns - upper_sizeOfColumns) + halfseriesSize_upper
//	
//	duplicate/O/R=[0,0][upper_0,upper_1] Kymograph_upper,Kymo_upper_zero
//	wavestats/Q/M=1 Kymo_upper_zero
//	xm=V_avg
//	for (j=Upper_centralShift[0]- halfseriesSize_upper;j<=Upper_centralShift[0]+ halfseriesSize_upper;j=j+1)
//		sx = sx + (Kymograph_upper[0][j] - xm)^2
//	endfor
//	
//	for (i=0;i<sizeOfFrames;i=i+1)
//		upper_0 = Upper_centralShift[i]-(sizeOfColumns - upper_sizeOfColumns)- halfseriesSize_upper
//		upper_1 = Upper_centralShift[i]-(sizeOfColumns - upper_sizeOfColumns)+ halfseriesSize_upper
//		duplicate/O/R=[i,i][upper_0,upper_1] Kymograph_upper,Kymo_upper_temp
//		wavestats/Q/M=1 Kymo_upper_temp
//		ym = V_avg
//		
//		for (d=-maxDelay_upper;d<=maxDelay_upper;d=d+1)
//			sxy=0
//			sy=0
//			for (j=Upper_centralShift[0]-halfseriesSize_upper;j<=Upper_centralShift[0]+halfseriesSize_upper;j=j+1)
//				if (j+d<0 || j+d>=Upper_sizeOfColumns)
//					continue
//				else
//					sxy = sxy + (Kymograph_upper[0][j] - xm)*(Kymograph_upper[i][j-Upper_centralShift[0]+Upper_centralShift[i]+d] - ym)
//					sy = sy + (Kymograph_upper[i][j-Upper_centralShift[0]+Upper_centralShift[i]+d] - ym)^2
//				endif
//			endfor
//			Upper_CrossCor[d+maxDelay_upper][i] = sxy/(sx*sy)^0.5
//		endfor
//	
//		CurveFit/Q gauss, Upper_CrossCor[*][i]/X=upper_delayNum
//		rename W_coef tempwave
//		wave tempwave
//		Upper_centralShift_temp[i] = tempwave[2]
//		killwaves root:Analysis:W_sigma
//		killwaves root:Analysis:tempwave
//		Upper_centralShift[i] = Upper_centralShift[i] + Upper_centralShift_temp[i]
//	endfor
//	
//	for (i=sizeOfFrames-1;i>=0;i=i-1)
//		Upper_centralShift[i] = (Upper_centralShift[i] - Upper_centralShift[0])*pixelSize
//	endfor
	
	
	display/w=(370,55+95*1.618,520,55+190*1.618) as "Lower part"
	newimage/f/HOST=# Kymograph_lower
	movesubwindow fnum=(0,0,150,95*1.618)
		
	CrossCorrelation("Kymograph_lower","_lower",0)
	wave centralShift_lower = root:Analysis:centralShift_lower
	
	CrossCorrelation("Kymograph_lower","_lowerEdge",1)
	wave centralShift_lowerEdge = root:Analysis:centralShift_lowerEdge
	
//		//for Kymograph_lower
//	make/N=(sizeOfFrames) lower_centralShift
//	for (i=0;i<sizeOfFrames;i=i+1)
//		CurveFit/Q gauss, kymograph_lower[i][*]
//		rename W_coef wavetemp
//		wave wavetemp
//		lower_centralShift[i] = round(wavetemp[2])
//		killwaves/z wavetemp
//	endfor	
//	
//	lower_sizeOfColumns = dimsize(Kymograph_lower,1)
//	make/N=(2*maxDelay_lower+1,sizeOfFrames) lower_CrossCor
//	
//	make/N=(2*maxDelay_lower+1) lower_delayNum
//	for(i=0;i<2*maxDelay_lower+1;i=i+1)
//		lower_delayNum[i] = i - maxDelay_lower
//	endfor
//	
//	make/N=(sizeOfFrames) lower_centralShift_temp
//	duplicate/O/R=[0,0][lower_centralShift[0]-halfseriesSize_lower,lower_centralShift[0]+halfseriesSize_lower] Kymograph_lower,Kymo_lower_zero
//	wavestats/Q/M=1 Kymo_lower_zero
//	xm=V_avg
//	sx=0
//	for (j=lower_centralShift[0]-halfseriesSize_lower;j<=lower_centralShift[0]+halfseriesSize_lower;j=j+1)
//		sx = sx + (Kymograph_lower[0][j] - xm)^2
//	endfor
//	
//	for (i=0;i<sizeOfFrames;i=i+1)
//		duplicate/O/R=[i,i][lower_centralShift[i]-halfseriesSize_lower,lower_centralShift[i]+halfseriesSize_lower] Kymograph_lower,Kymo_lower_temp
//		wavestats/Q/M=1 Kymo_lower_temp
//		ym = V_avg
//		
//		for (d=-maxDelay_lower;d<=maxDelay_lower;d=d+1)
//			sxy = 0
//			sy = 0
//			for (j=lower_centralShift[0]-halfseriesSize_lower;j<=lower_centralShift[0]+halfseriesSize_lower;j=j+1)
//				if (j+d<0 || j+d>=lower_sizeOfColumns)
//					continue
//				else
//					sxy = sxy + (Kymograph_lower[0][j] - xm)*(Kymograph_lower[i][j-lower_centralShift[0]+lower_centralShift[i]+d] - ym)
//					sy = sy + (Kymograph_lower[i][j-lower_centralShift[0]+lower_centralShift[i]+d] - ym)^2
//				endif
//			endfor
//			lower_CrossCor[d+maxDelay_lower][i] = sxy/(sx*sy)^0.5
//		endfor
//	
//		CurveFit/Q gauss, lower_CrossCor[*][i]/X=lower_delayNum
//		rename W_coef tempwave
//		wave tempwave
//		lower_centralShift_temp[i] = tempwave[2]
//		killwaves root:Analysis:W_sigma
//		killwaves root:Analysis:tempwave
//		lower_centralShift[i] = lower_centralShift[i] + lower_centralShift_temp[i]
//	endfor
//	
//	for (i=sizeOfFrames-1;i>=0;i=i-1)
//		lower_centralShift[i] = (lower_centralShift[i] - lower_centralShift[0])*pixelSize
//	endfor
	
	display/w=(530,10,800,10+90*1.618) centralShift_upper vs timeSeries
	appendtograph centralShift_lower vs timeSeries
	appendtograph centralShift_upperEdge vs timeSeries
	appendtograph centralShift_lowerEdge vs timeSeries
	ModifyGraph rgb(centralShift_lower)=(0,9472,39168)
	ModifyGraph rgb(centralShift_upper)=(0,9472,39168)
	ModifyGraph lstyle(centralShift_lower)=3
	ModifyGraph lstyle(centralShift_lowerEdge)=3
	Label bottom "Time (s)"
	Label left "Distance shifted (nm)"
end

Function CrossCorrelation(KymographName,NameAddOn,mode)
 // mode defines whether targets bright band or dark boundary for cross-correlation, 0 is bright band, 1 is dark boundary 
 //The length of segment used for cross-correlation is 2*halfseriesSize+1 
 //maxDelay defines the max number of pixel one segment can move away from the other during correlation
	string KymographName,NameAddOn
	variable mode
	variable/G halfseriesSize
	variable/G maxDelay

	if(strlen(WaveList(KymographName,";",""))==0)
		Doalert 0,"Wrong wave name input, man"
		return -1
	endif
	
	string pathOfKymo = "root:Analysis:" + KymographName 
	string centralShift = "centralShift" + NameAddOn
	string crossCor = "crossCor"+ NameAddOn
	string delayNum = "delayNum" + NameAddOn
	string centralShift_temp = "centralShift_temp" + NameAddOn
	string halfseriesSizeName = "halfseriesSize" + NameAddOn
	string maxDelayName = "maxDelay" + NameAddOn
	string centralshift_band = "centralshift_band" + NameAddOn
	wave Kymowave = $pathOfKymo
	 
	variable xm,ym,sx,sy,sxy,d,i,j
	variable sizeOfColumns_part = dimsize(Kymowave,1)
	NVAR sizeOfFrames = root:RawData:sizeOfFrames
	NVAR sizeOfColumns = root:RawData:sizeOfColumns
	NVAR pixelSize = root:Analysis:pixelSize
	variable temp_0,temp_1
	
	variable/G offset
	if (dimOffset(Kymowave,1)>0)
		offset = sizeOfColumns - sizeOfColumns_part
	elseif (dimOffset(Kymowave,1)==0)
		offset = 0
	endif
	
	
	make/O/N=(sizeOfFrames) tempwave1 //centralShift
	
	if (mode==0) //detect bright band
		make/O/N=(sizeOfFrames) tempwave0
		for (i=0;i<sizeOfFrames;i=i+1)
			CurveFit/Q gauss, Kymowave[i][*]
			rename W_coef wavetemp1
			wave wavetemp1
			tempwave0[i] = wavetemp1[3]
			tempwave1[i] = round(wavetemp1[2])
			killwaves/z wavetemp1
		endfor
		wavestats/Q tempwave0
		
		halfseriesSize = ceil(V_avg - V_adev)
		if (halfseriesSize < 10)
			halfseriesSize = 10
		elseif (halfseriesSize > 30)
			halfseriesSize = 30
		endif
		maxdelay = halfseriesSize
		killwaves/z tempwave0
		duplicate tempwave1 $centralshift_band
		appendToGraph/C=(0,0,65000) $centralshift_band
		temp_0 = tempwave1[0]-offset - halfseriesSize
		temp_1 = tempwave1[0]-offset + halfseriesSize
	elseif (mode==1) //Detect edge
		FindThreeStates(KymographName)
		string pathOfEdge1 = "root:Analysis:" + KymographName + "Edge1"
		string pathOfEdge2 = "root:Analysis:" + KymographName + "Edge2" 
		wave Edge1 = $pathOfEdge1
		wave Edge2 = $pathOfEdge2
		wave tempwave1 = root:Analysis:tempwave1
		AppendToGraph Edge1
		AppendToGraph Edge2 
		halfseriesSize = 5
		maxDelay = 5
		temp_0 = tempwave1[0] - halfseriesSize
		temp_1 = tempwave1[0] + halfseriesSize
	else 
		Doalert 0, "we only have mode 0 and 1..."
	endif
	
	make/N=(2*maxDelay+1,sizeOfFrames) tempwave2 //$crossCor
	
	make/N=(2*maxDelay+1) tempwave3 //$delayNum
	for(i=0;i<2*maxDelay+1;i=i+1)
		tempwave3[i] = i - maxDelay
	endfor
	
	
	make/N=(sizeOfFrames) tempwave4 //$centralShift_temp
	
	duplicate/O/R=[0,0][temp_0,temp_1] Kymowave,Kymo_zero
	wavestats/Q/M=1 Kymo_zero
	xm=V_avg
	
	for (j=temp_0;j<=temp_1;j=j+1)
		sx = sx + (Kymowave[0][j] - xm)^2
	endfor
	
	variable temp_2, temp_3
	for (i=0;i<sizeOfFrames;i=i+1)
		if (mode==0)
			temp_2 = tempwave1[i]-offset- halfseriesSize
			temp_3 = tempwave1[i]-offset+ halfseriesSize
		elseif (mode==1)
			temp_2 = tempwave1[i]- halfseriesSize
			temp_3 = tempwave1[i]+ halfseriesSize
		endif
		
		duplicate/O/R=[i,i][temp_2,temp_3] Kymowave,Kymowave_temp
		wavestats/Q/M=1 Kymowave_temp
		ym = V_avg
		
		for (d=-maxDelay;d<=maxDelay;d=d+1)
			sxy=0
			sy=0
			for (j=temp_0;j<=temp_1;j=j+1)
				if (j-tempwave1[0]+tempwave1[i]+d<0 || j-tempwave1[0]+tempwave1[i]+d>sizeOfColumns_part-1)
					continue
				else
					sxy = sxy + (Kymowave[0][j] - xm)*(Kymowave[i][j-tempwave1[0]+tempwave1[i]+d] - ym)
					sy = sy + (Kymowave[i][j-tempwave1[0]+tempwave1[i]+d] - ym)^2
				endif
			endfor
			tempwave2[d+maxDelay][i] = sxy/(sx*sy)^0.5
		endfor
	
		CurveFit/Q gauss, tempwave2[*][i]/X=tempwave3
		rename W_coef wavetemp
		wave wavetemp
		tempwave4[i] = wavetemp[2]
		killwaves root:Analysis:W_sigma
		killwaves root:Analysis:wavetemp
		tempwave1[i] = tempwave1[i] + tempwave4[i]
	endfor
	

	for (i=sizeOfFrames-1;i>=0;i=i-1)
		tempwave1[i] = (tempwave1[i] - tempwave1[0])*pixelSize
	endfor
	rename tempwave1 $centralShift
	rename halfseriesSize $halfseriesSizeName
	rename maxDelay $maxDelayName
	killwaves/z tempwave2
	killwaves/z tempwave3
	killwaves/z tempwave4 
	killwaves/z Kymo_zero
	Killwaves/z Kymowave_temp
//	rename tempwave2 $crossCor
//	rename tempwave3 $delayNum
//	rename tempwave4 $centralShift_temp
End

Function FindThreeStates(KymographName)
  //this function serves to find two boundary points to segregate each slice of Kymograph into three parts 
  	string kymographName 
 	string pathOfKymo = "root:Analysis:" + KymographName 
 	string Edge1name = KymographName + "Edge1"
 	string Edge2name = KymographName + "Edge2"
	wave Kymowave = $pathOfKymo
	variable sizeOfColumns_part = dimsize(Kymowave,1)
	NVAR sizeOfFrames = root:RawData:sizeOfFrames
	NVAR offset = root:Analysis:offset
	variable i,j,k,tempvar1,tempvar2
	make/O/N=(sizeOfColumns_part,sizeOfColumns_part,4,sizeOfFrames) bestFitOfSteps//this matrix stores the best fitting parameter for each row of wave "column_addup_part"
	//first layer is sum of chi-square, second layer is average of the first segment
	//third layer is average of the second segment, forth layer is average of the third segment
	//rows num equals the point of the first step, columns num is the point of the second step 
	bestFitOfSteps = Nan
	make/O/N=(sizeOfFrames) Edge1
	make/O/N=(sizeOfFrames) Edge2
	make/O/N=(sizeOfFrames) tempwave1 
	
	for (i=0;i<sizeOfFrames;i=i+1)
		for (j=0;j<sizeOfColumns_part-2;j=j+1)
			imagestats/G={i,i,0,j} Kymowave
			tempvar1 = (V_npnts*V_adev)^2/V_avg
			tempvar2 = V_avg
			for (k=j+1;k<sizeOfColumns_part-1;k=k+1)
				bestFitOfSteps[j][k][0][i] = tempvar1
				bestFitOfSteps[j][k][1][i] = tempvar2
				imagestats/G={i,i,j+1,k} Kymowave
				bestFitOfSteps[j][k][0][i] = bestFitOfSteps[j][k][0][i] + (V_npnts*V_adev)^2/V_avg
				bestFitOfSteps[j][k][2][i] = V_avg
				imagestats/G={i,i,k+1,sizeOfColumns_part-1} Kymowave 
				bestFitOfSteps[j][k][0][i] = bestFitOfSteps[j][k][0][i] + (V_npnts*V_adev)^2/V_avg
				bestFitOfSteps[j][k][3][i] = V_avg
			endfor
		endfor
		duplicate/O/R=[][][0,0][i,i] bestFitOfSteps bestFitOfSteps_slide 
		imagestats/M=1 bestFitOfSteps_slide
		Edge1[i] = V_minRowLoc + offset
		Edge2[i] = V_minColLoc + offset
		
		if (dimOffset(Kymowave,1)==0)
			tempwave1[i] = V_minColLoc +5
		elseif ((dimOffset(Kymowave,1)>0))
			tempwave1[i] = V_minRowLoc
		endif
	endfor
	
	rename Edge1 $Edge1name
	rename Edge2 $Edge2name
end
			
 
 
//Function BleachedArea():GraphMarquee
//	setdatafolder root:Analysis
//	variable/G startrow,endrow,startcolumn,endcolumn
//	variable/G BleachedIntensity = 0 
//	variable marqueeWidth,marqueeHeight
//	variable i,j
//	wave Kymograph = root:rawData:Kymograph
//	getMarquee left,bottom
//	startrow = floor(v_bottom)
//	endrow = ceil(v_top)
//	startcolumn = floor(v_left)
//	endcolumn = ceil(v_right)
//	if (startcolumn < 0)
//		startcolumn = 0
//	endif
//	marqueeWidth = endColumn-startColumn+1
//	marqueeHeight = endRow-startRow+1
//	
//	for (i=startcolumn;i<=endcolumn;i=i+1)
//		for (j=startrow;j<=endrow;j=j+1)
//			BleachedIntensity = BleachedIntensity + Kymograph[i][j]
//		endfor
//	endfor
//	
//	BleachedIntensity = BleachedIntensity/MarqueeWidth/MarqueeHeight
//end	
	
Function WindowsExplorer()
	setdatafolder root:WindowsExplorer
	string path = "C:Users:lc:Desktop:For LuChen Data2"
	string subpath
	SVAR fileName = root:WindowsExplorer:fileName
	SVAR filePath = root:WindowsExplorer:filePath
	SVAR format = root:WindowsExplorer:format
	string fileList
	
	Prompt path, "Where to load your data?"
	DoPrompt "Enter your path", path
	if (V_Flag) //user canceled
		return -1
	endif
	
	Do
		newpath/O/Q DataPath, path 
		fileList = indexedDir(Datapath,-1,0)
		fileList = fileList + indexedfile (Datapath,-1,format) 
		Prompt subpath "Which would you choose?" popup fileList
		DoPrompt "Choose your path" subpath
		if (V_Flag)
			return -1
		endif
	
		if (stringmatch (subpath,"!*.*"))	
			path = path+":"+subpath
		else
			break
		endif
	while(1)

	path = path+":"+subpath
	filePath = path
	fileName = ParseFilePath (3, filePath,":",0,0)
	fileName = CleanupName (fileName,0)
end

Function KillData(CtrlName):ButtonControl
	String CtrlName
	String graphName
	killwindows(0)
	
	setdatafolder root:
	Killdatafolder/z root:WindowsExplorer
	Killdatafolder/z root:RawData
	Killdatafolder/z root:Analysis
	EnableDisableAnalysisButton()
end

Function ClearData(CtrlName):ButtonControl
	string CtrlName
	dowindow/F graph0
	killwindows(1)
	Killdatafolder/z root:Analysis
	newdatafolder/O root:Analysis
end

Function killwindows(value)
	variable value
	String graphName
	Do
		graphName=winName(value,1)
		if (strlen(graphName)==0)
			break
		endif
		killwindow $graphName
	While(1)
end

Function PressEscapeToAbort()
	Variable doAbort = 0
	if (GetKeyState(0) & 32) // Is Escape key pressed now?
		doAbort = 1
	endif
	return doAbort
End