Recursive uneven surface sampling

This script recursively tiles a surface, making finer tiling where it encounters larger local area curvature. Local area curvature is defined as= (maximum angle between the surface normal in the middle of a tile and any vector from the middle of the tile to the corner of the tile) - 90.

It will keep on subdividing untill this angle is below the treshold, or if the max number of iterations is reached.

It returns a 2D-list of UV-coordinates, where every subset contains 4 UV-corners in clockwise order. In it's current version, it deletes the former surface and builds facets on that lists for visualisation.

Problems: with large initial tiling, "gaps" might appear. (see pic below)

### Recursive uneven surface tiling (based on curvature)
### Sebastiaan Leenknegt, 2012.04.25
### Please reference properly when using (parts of) this script

import rhinoscriptsyntax as rs
import math

def MakeSample(surf,uDim,vDim):
    uDomain=rs.SurfaceDomain(surf,0)
    vDomain=rs.SurfaceDomain(surf,1)
    uStep=(uDomain[1]-uDomain[0])/uDim
    vStep=(vDomain[1]-vDomain[0])/uDim
   
    sampleList=[]
    for i in range(0,uDim):
        for j in range(0,vDim):
            uValue1=uStep*i
            uValue2=uStep*(i+1)
            vValue1=vStep*j
            vValue2=vStep*(j+1)
            sampleList.append([[uValue1,vValue1],[uValue1,vValue2],[uValue2,vValue2],[uValue2,vValue1]])
    return sampleList

def SubdivideTile(points):
    point4=[(points[0][0]+points[1][0])/2,(points[0][1]+points[1][1])/2]
    point5=[(points[1][0]+points[2][0])/2,(points[1][1]+points[2][1])/2]
    point6=[(points[2][0]+points[3][0])/2,(points[2][1]+points[3][1])/2]
    point7=[(points[3][0]+points[0][0])/2,(points[3][1]+points[0][1])/2]
    point8=[(points[0][0]+points[2][0])/2,(points[0][1]+points[2][1])/2]
    tile0=[points[0],point4,point8,point7]
    tile1=[point4,points[1],point5,point8]
    tile2=[point8,point5,points[2],point6]
    tile3=[point7,point8,point6,points[3]]
    return [tile0,tile1,tile2,tile3]

def SurfaceTiling(surf,tileList,treshold,check):
    totalmaxangle=treshold
    tileListCopy=tileList[:]
    for i in range(0,len(tileList)):
        tile=tileList[i]
        maxangle=treshold
        midpoint=[(tile[0][0]+tile[2][0])/2,(tile[0][1]+tile[2][1])/2]
        realmidpoint=rs.EvaluateSurface(surf,midpoint[0],midpoint[1])
        normalmid=rs.SurfaceCurvature(surf,midpoint)[1]
        for point in tile:
            realpoint=rs.EvaluateSurface(surf,point[0],point[1])
            vectfrommid=rs.VectorCreate(realpoint,realmidpoint)
            tempangle=abs(rs.VectorAngle(normalmid,vectfrommid))
            angle=abs(tempangle-90)
            if angle>maxangle:
                maxangle=angle
        checkdist=abs(math.sqrt(math.pow(tile[0][0]-tile[2][0],2)+math.pow(tile[0][1]-tile[2][1],2)))
        if (maxangle>treshold) and (checkdist>check):
            newtiles=SubdivideTile(tile)
            tileListCopy.insert(i,newtiles[3])
            tileListCopy.insert(i,newtiles[2])
            tileListCopy.insert(i,newtiles[1])
            tileListCopy.insert(i,newtiles[0])
            tileListCopy.remove(tile)
        if maxangle>totalmaxangle:
            totalmaxangle=maxangle
    return (tileListCopy,totalmaxangle)

def Master(maxDev):
    obj=rs.GetObject("Pick Surface")
    mindist=rs.SurfaceDomain(obj,0)[0]
    sample=MakeSample(obj,5,5)
    for i in range(0,5):
        newsample=SurfaceTiling(obj,sample,maxDev,mindist/1000)
        sample=newsample[0]
        if newsample[1]<maxDev:break
    rs.EnableRedraw(False)
    for pts in sample:
        realptList=[]
        for pt in pts:
            realptList.append(rs.EvaluateSurface(obj,pt[0],pt[1]))
        rs.AddSrfPt(realptList)
    rs.DeleteObject(obj)
    return realptList

Master(2)