filament.est=function(data,eps,a,b,q.poisson,FDR)

	{

	n=dim(data)[1]

	vol=pi*eps^2											#volume of ball

	cutoff=qpois(q.poisson,n*vol)									#threshold

	####START FDR Calculation

	h.max=c()

	for (j in 1:100)

		{

		stuff=c()	

		atad=cbind(runif(n),runif(n))

		for (i in 1:n)
	
			{

			x=atad[i,1]

			y=atad[i,2]

			ball=atad[which(((atad[,1]-x)^2+(atad[,2]-y)^2)<eps^2),]

			if (dim(ball)[1]>cutoff)
	
				{
	
				axis=princomp(ball)
	
				x1=axis$l[1,1]

				y1=axis$l[2,1]

				x2=mean(ball[,1])

				y2=mean(ball[,2])

				weight=(axis$sdev[1]-axis$sdev[2])/eps^5

				stuff=rbind(stuff,c(x2,y2,atan(y1/x1),weight))

				}

			}

		#####estimate h(x), using elliptical kernel, constant weights in kernel

		m=dim(stuff)[1]

		stuff=cbind(stuff[,1:4],rep(0,m))

		atad=stuff[,1:2]

		for (i in 1:m)
	
			{

			center=stuff[i,1:2]

			theta=stuff[i,3]

			r.m=matrix(c(cos(theta),sin(theta),-sin(theta),cos(theta)),nrow=2)

			c.x=cbind(atad[,1]-center[1],atad[,2]-center[2])

			new.data=c.x%*%r.m

			x=new.data[,1]

			y=new.data[,2]

			ball=stuff[which(x^2/b^2+y^2/a^2<=1),]

			if (length(ball)>=10)

				{
	
				stuff[i,5]=sum(ball[abs(ball[,3]-theta)<.5,4])

				}

			}

		h.max[j]=max(stuff[,5])

		}

	thresh=quantile(h.max,1-FDR)	

	###END FDR Calculation

	plot(c(0,1),c(0,1),t='n')

	points(data)	

	#####calculate weights for each ball

	stuff=c()

	for (i in 1:n)
	
		{
	
		x=data[i,1]

		y=data[i,2]

		ball=data[which(((data[,1]-x)^2+(data[,2]-y)^2)<eps^2),]

		if (dim(ball)[1]>cutoff)

			{
	
			axis=princomp(ball)
	
			x1=axis$l[1,1]

			y1=axis$l[2,1]

			x2=mean(ball[,1])

			y2=mean(ball[,2])

			weight=(axis$sdev[1]-axis$sdev[2])/eps^5

			stuff=rbind(stuff,c(x2,y2,atan(y1/x1),weight))

			}

		}
	eps2=eps/3						#radius of balls on which to average theta
	
	plot.est(stuff,a,b,eps2,thresh)

	}

###############################

plot.est=function(stuff,a,b,eps,thresh)

	{

	x1=seq(0,1,.01)

	x2=seq(0,1,.01)

	n=length(x1)

	h=(1:n)%*%t(1:n)*0

	for (i in 1:n)

		{

		for (j in 1:n)

			{

			center=c(x1[i],x2[j])

			if (length(which((stuff[,1]-center[1])^2+(stuff[,2]-center[2])^2<=eps^2))>0)

				{

				theta=mean(stuff[which((stuff[,1]-center[1])^2+(stuff[,2]-center[2])^2<=eps^2),3])

				r.m=matrix(c(cos(theta),sin(theta),-sin(theta),cos(theta)),nrow=2)

				c.x=cbind(stuff[,1]-center[1],stuff[,2]-center[2])

				new.data=c.x%*%r.m

				x=new.data[,1]

				y=new.data[,2]

				ball=stuff[which(x^2/b^2+y^2/a^2<=1),]

				if (length(ball)>=10)

					{

					h[i,j]=sum(ball[abs(ball[,3]-theta)<.5,4])

					}

				}

			}

		}

	goods=c()

	for (i in 1:n)

		{

		for (j in 1:n)

			{

			if (h[i,j]>thresh)

				goods=rbind(goods,c((i-1)/(n-1),(j-1)/(n-1)))

			}

		}

	points(goods,col='blue',pch=19)

	}

################################




n=1000													#sample size

delta=.15 												#proportion of data assigned to filament

sigma=.02 												#noise st dev in both directions

data=cbind(runif(n*(1-delta)),runif(n*(1-delta))) 							#generate background noise

u=runif(n*delta,0,1)

data=rbind(data,cbind(u+rnorm(n*delta,0,sigma),sqrt(u^2+.15)+sin(u*5)/3+rnorm(n*delta,0,sigma)))	#concatonate background noise with filament

eps=.10 												#radius of ball

q.poisson=.8												#quantile for threshold for considering ball

a=.02		   											#minor axis length is 2a
	
b=.05													#major axis is 2b

FDR=.1													#false detection rate

filament.est(data,eps,a,b,q.poisson,FDR)

u=seq(0,1,.01)

lines(u,sqrt(u^2+.15)+sin(u*5)/3)