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