function k_means,xys,nGr

	epsilon=1e-3

	numPoints=n_elements(xys)/2

	rndIDS=randomu(seed,nGr)*numPoints

	MeansX=xys[0,rndIDS]
	MeansY=xys[1,rndIDS]


	repeat begin

		distances=fltarr(nGr,numPoints)
		for j=0,nGr-1 do begin
			distances[j,*]=sqrt((xys[0,*]-MeansX[j])^2+(xys[1,*]-MeansY[j])^2)
		endfor

		groupID=intarr(numPoints)
		for j=0,numPoints-1 do groupID[j]=(sort(distances[*,j]))[0]

		newMeansX=fltarr(nGr)
		newMeansY=fltarr(nGr)
		movement=fltarr(nGr)

		for j=0,nGr-1 do begin
			idsMios=where(groupID eq j,totalMisIDs)
			if idsMios[0] ne -1 then begin
				newMeansX[j]=total(xys[0,idsMios])/float(totalMisIDs)
				newMeansY[j]=total(xys[1,idsMios])/float(totalMisIDs)
				movement[j]=sqrt((newMeansX[j]-MeansX[j])^2+(newMeansY[j]-newMeansY[j])^2)
			endif else begin
				rndIDS=randomu(seed)*numPoints
				newMeansX[j]=xys[0,rndIDS]
				newMeansY[j]=xys[1,rndIDS]
				movement[j]=1000
			endelse
		endfor

		MeansX=newMeansX
		MeansY=newMeansY

	endrep until total(movement) le epsilon

	sortOrd=sort(newMeansX)
	newMeansX=newMeansX[sortOrd]
	newMeansY=newMeansY[sortOrd]

	resultKM={	MeansX:newMeansX,$
				MeansY:newMeansY,$
				groupID:groupID}

return,resultKM
end

;xys=transpose([[10*randomu(seed,100)+50,10*randomu(seed,100)+25,5*randomu(seed,100)+15],[10*randomu(seed,100)+2,10*randomu(seed,100)+20,5*randomu(seed,100)+15]]) & print,k_means(xys,3)