[Numpy-discussion] scipy image processing memory leak in python 2.7?

[Joseph McGlinchy]

Hi numpy list!

I am trying to do some image processing on a number of images, 72 to be specific. I am seeing the python memory usage continually increase. I will say this is being done also using arcpy so that COULD be an issue, but want to rule out scipy first since I am only using arcpy to grab a set of polygons and convert them to a binary image, and to save images to disk. I pass the raster to my numpy/scipy processing. The code is organized as:


1)      Main loop iterating over features

2)      Function which does image processing using scipy

3)      Class which implements some of the image processing features

I placed 2 sleep calls in the code.. one at the end of the image processing function and one at the end of each iteration in the main loop. I am seeing a 20MB memory increase each iteration while looking at my processes in Task Manager:

Iteration 1: 204,312K and 209,908K, respectively
Iteration 2: 233,728K and 230,192K, respectively
Iteration 3: 255,676K and 250,600K, respectively

Any ideas?

Thanks so much!

Joe


1)      Main loop
#for each feature class, convert to raster by "Value" and reclass to 1 and 0
temp_raster = os.path.join(scratchWorkspace,"temp_1.tif")
temp_raster_elev = os.path.join(scratchWorkspace,"temp_elev.tif")
temp_reclass = os.path.join(scratchWorkspace,"temp_reclass.tif")
remap = arcpy.sa.RemapValue([[1,1], ["NoData",0]])


for i,fc in enumerate(potWaterFeatNames):


    #try to delete the temp rasters in case they exist
    try:
        #delete temp rasters
        arcpy.Delete_management(temp_raster)
        arcpy.Delete_management(temp_raster_elev)
        arcpy.Delete_management(temp_reclass)
    except e:
        arcpy.AddMessage(e.message())

    #clear in memory workspace"
    arcpy.Delete_management("in_memory/temp_table")

    #print "on feature {} of {}".format(i+1, len(potWaterFeatNames))
    arcpy.AddMessage("on feature {} of {}".format(i+1, len(potWaterFeatNames)))
    arcpy.AddMessage("fc = {}".format(fc))


    arcpy.PolygonToRaster_conversion(fc, "Value", temp_raster, "", "", 1)
    outReclass = arcpy.sa.Reclassify(temp_raster,"Value",remap)
    outReclass.save(temp_reclass)

    del outReclass

    #call function to process connected components
    try:
        proc_im = processBinaryImage(temp_reclass, scratchWorkspace, sm_pixels)
    except Exception as e:
        print "FAILED! scipy image processing"
        arcpy.AddMessage(e.message())

    #convert processed raster to polygons
    try:
        out_fc = os.path.join(outWorkspace,fc + "_cleaned")
        arcpy.RasterToPolygon_conversion(proc_im, out_fc)
        arcpy.Delete_management(proc_im)
    except Exception as e:
        print "FAILED! converting cleaned binary raster to polygon"
        arcpy.AddMessage(e.message())



    #delete zero value polyons, gridcode = 0
    try:
        uc = arcpy.UpdateCursor(out_fc, "gridcode=0",arcpy.Describe(out_fc).spatialReference)

        #loop through 0 rows and delete
        for row in uc:
            uc.deleteRow(row)

        del uc
    except Exception as e:
        print "FAILED! deleting rows from table"
        arcpy.AddMessage(e.message())


    #check that number of polygons with gridcode = 1 is greater than 0
    count = arcpy.GetCount_management(out_fc)

    #add elevation field back in
    if (count > 0):
        arcpy.PolygonToRaster_conversion(fc, "z_val", temp_raster_elev, "", "", 1)
        arcpy.sa.ZonalStatisticsAsTable(out_fc, "Id", temp_raster_elev, "in_memory/temp_table","#","MEAN")
        arcpy.JoinField_management(out_fc, "Id", "in_memory/temp_table", "Id", ["MEAN"])
    else:
        arcpy.Delete_management(out_fc)

    #delete temp rasters
    arcpy.Delete_management(temp_raster)
    arcpy.Delete_management(temp_raster_elev)
    arcpy.Delete_management(temp_reclass)

    #python garbage collection
    #collected = gc.collect()
    #print "Garbage collector: collected %d objects." % (collected)

    print "sleeping for 10 seconds"
   time.sleep(10)



2)      Image processing function
def processBinaryImage(imageName, save_dir, sm_pixels):


    fname = os.path.basename(imageName)
    imRaster = arcpy.Raster(imageName)

    #Grab AGS info from image
    #use describe module to grab info
    descData = arcpy.Describe(imRaster)
    cellSize = descData.meanCellHeight
    extent = descData.Extent
    spatialReference = descData.spatialReference
    pnt = arcpy.Point(extent.XMin, extent.YMin)

    del imRaster

    converter = ConvertRasterNumpy();
    imArray = converter.rasterToNumpyArray(imageName)
    imMin = np.min(imArray)
    imMax = np.max(imArray)
    print imMin, imMax
    arcpy.AddMessage("image min: " + str(imMin))
    arcpy.AddMessage("image max: " + str(imMax));



    #other flags
    show_image = False #verbose
    save_flag = False #under verbose
    gray = False #threshold but keep gray levels
    make_binary = False #flag if grayscale is true, to binarize final image
    save_AGS_rasters = True




    #create filter object
    filter = thresholdMedianFilter()
    lowValue = 0
    tImage = filter.thresholdImage(imArray, lowValue, gray)




    #median filter image
    filtImage1 = filter.medianFiltImage(tImage,1)

    #create structuring element for morphological operations
    sElem = np.array([[0., 1., 0.],
                      [1., 1., 1.],
                      [0., 1., 0.]])

    #open filtered image
    gray = False
    mImage = filter.mOpenImage(filtImage1, sElem, gray)



    #set list to store info change
    num_its = 100
    the_change = np.zeros(num_its)
    for it in range(100):
        prev = mImage
        filtImage = filter.medianFiltImage(mImage,1)
        mImage = filter.mOpenImage(filtImage, sElem, gray)

        #calculate difference
        (m_info, z_info) = filter.calcInfoChange(prev, mImage)
        the_change[it] = z_info
        del filtImage
        del prev

        #if the change is less than 5% of the initial change, exit
        cur_perc = the_change[it]/the_change[0]*100
        if cur_perc < 5.0:
            print "exiting filter on iteration " + str(it)
            print "change is less than 5% (this iteration: " + str(cur_perc) + ")"
            break


    ############################################################################
    # now we have a binary mask. Let's find and label the connected components #
    ############################################################################

    #clear some space
    del tImage

    del filtImage1
    del m_info

    label_im_init, nb_labels = ndimage.label(mImage)



    #Compute size, mean_value, etc. of each region:
    label_im = label_im_init.copy()
    del label_im_init

    sizes = ndimage.sum(mImage, label_im, range(nb_labels + 1))
    mean_vals = ndimage.sum(imArray, label_im, range(1, nb_labels + 1))

    #clean up small components
   mask_size = sizes < sm_pixels
    remove_pixel = mask_size[label_im]
    label_im[remove_pixel] = 0;
    labels = np.unique(label_im)
    label_im = np.searchsorted(labels, label_im)

    #make label image to a binary image, and convert to arcgis raster
    label_im[label_im > 0] = 1
    label_im = np.array(label_im, dtype = 'float32')
    print label_im.dtype

    saveit = False
    if ~saveit:
        outRaster = save_dir + "\\" + fname[:-4] + "_CC_" + str(sm_pixels) + ".tif"
        temp = arcpy.NumPyArrayToRaster(label_im,pnt,cellSize,cellSize)
        arcpy.DefineProjection_management(temp, spatialReference)
        arcpy.CopyRaster_management(temp, outRaster, "DEFAULTS","0","","","","8_BIT_UNSIGNED")


    #clear more space
    del mImage
    del nb_labels
    del sizes
    del mean_vals
    del mask_size
    del remove_pixel
    del label_im
    del labels
    del temp

    del the_change
    del sElem
    del filter
    del imArray

    print 'sleeping'
    time.sleep(20)

    return outRaster



3)      Image processing class

class thresholdMedianFilter():

    def thresholdImage(self, imArray, thresh, binary = True):
        """ threshold the image. values equal or below thresh are 0, above are 1"""

        tImage = imArray
        tImage[tImage <= thresh] = 0
        if binary:
            tImage[tImage > thresh] = 1
        return tImage


    def medianFiltImage(self,imArray,n=1):
        """median filter the image n amount of times. a single time is the default"""

        for n in range(n):
            prev = imArray
            imArray = signal.medfilt2d(prev)
            del prev

        return imArray

    def mOpenImage(self,imArray, sElem, gray = False):
        """ morphological opening """
        #Mnp = np.array(morph.binary_erosion(imArray, sElem))
        #Mnp = np.array(morph.binary_dilation(Mnp, sElem))
        if gray:
            imArray1 = np.array(morph.grey_dilation(imArray, structure = sElem))
            imArray2 = np.array(morph.grey_erosion(imArray1, structure = sElem), dtype = 'float32')
            del imArray1

            return imArray2
        else:
            Mnp1 = np.array(morph.binary_dilation(imArray, sElem))
            Mnp2 = np.array(morph.binary_erosion(Mnp1, sElem), dtype = 'float32')
            del Mnp1

            return Mnp2


    def calcInfoChange(self,imArray1, imArray2):
        """calculate entropy of an image"""
        diff = imArray1 - imArray2
        m_norm = sum(abs(diff)) #Manhattan Norm
        z_norm = norm(diff.ravel(), 0) #Zero norm
        del diff

        return (m_norm, z_norm)



Joe McGlinchy | Imagery Scientist
Database Services - 3D and Imagery Team
ESRI || 380 New York St. || Redlands, CA 92373 || USA
T 909-793-2853, ext. 4783

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20140128/8ea3e9db/attachment.html>