How to find out which polygon contains the most points using Python OGR?

I have two shapefiles: point (railway stations) and polygon (regions). My task is to write a script using OGR module in Python that can determine in which region the biggest amount of railway stations is.

I don't have any problem loading shp and reading its content, it's just the mechanics of solving the task that I can't seem to figure out.

I thought if I iterate over the regions (kind of abusing SetAttributeFilter), get their geometry, then in another loop get geometry of the points and finally do intersect (or maybe within, I am not sure, seemed to un/work the same), I get the result but I must be doing something wrong. Funny thing is, the first result I get is correct but following on it doesn't work as expected.

You can check the crucial part of my code (until here it works fine), I am a beginner so it may be ugly.

filters = []

for code in codes:
    filter1 = str("KOD_KRAJ='" + code + "'")
    filter2 = str('"'+filter1+'"')
    filters.append(filter2)


stations = 0

for filter in filters:
    regions.SetAttributeFilter(filter)

    for i in range(0,numberRe):
        regionFeat = regions.GetNextFeature()
        regionGeometry = regionFeat.GetGeometryRef()


        for r in range(0,numberRa):
            railwayFeature = railway.GetFeature(r)
            railwayGeometry = railwayFeature.GetGeometryRef()

            if railwayGeometry.Within(regionGeometry):
                stations = stations + 1
        print filter, str(stations)

Note: KOD_KRAJ is the header in shapefile's attribute table and variable "codes" is a list containing its instances (administrative code of each region). SetAttributeFilter should be used like this ("NAME='something'") which is why the weird first loop.

I just want to get it working to the stage when I get the numbers for each region, after that it will not be hard for me to do the rest.