manipulation - Cálculo de la imagen de la variación python

python image manipulation

Question

¿Existe una manera fácil de calcular un filtro de varianza en ejecución en una imagen usando Python / NumPy / Scipy? Al ejecutar la imagen de varianza quiero decir el resultado de calcular la suma ((I - mean (I)) ^ 2) / nPixels para cada subventana I en la imagen.

Dado que las imágenes son bastante grandes (12000x12000 píxeles), quiero evitar la sobrecarga de convertir las matrices entre formatos solo para poder usar una biblioteca diferente y luego convertir de nuevo.

Creo que podría hacer esto manualmente encontrando la media usando algo así como

kernel = np.ones((winSize, winSize))/winSize**2
image_mean = scipy.ndimage.convolve(image, kernel)
diff = (image - image_mean)**2
# Calculate sum over winSize*winSize sub-images
# Subsample result

pero sería mucho mejor tener algo así como la función stdfilt de Matlab.

¿Puede alguien señalarme en la dirección de una biblioteca que tiene esta funcionalidad Y admite matrices numpy, o insinuar / proporcionar una forma de hacerlo en NumPy / SciPy?

Answer 1

Puedes usar numpy.lib.stride_tricks.as_strided para obtener una vista en ventana de tu imagen:

import numpy as np
from numpy.lib.stride_tricks import as_strided
rows, cols = 500, 500
win_rows, win_cols = 5, 5
img = np.random.rand(rows, cols)
win_img = as_strided(img, shape=(rows-win_rows+1, cols-win_cols+1,
win_rows, win_cols),
strides=img.strides*2)

Y ahora win_img[i, j] es la (win_rows, win_cols) con la esquina superior izquierda en la posición [i, j] :

>>> img[100:105, 100:105]
array([[ 0.34150754, 0.17888323, 0.67222354, 0.9020784 , 0.48826682],
[ 0.68451774, 0.14887515, 0.44892615, 0.33352743, 0.22090103],
[ 0.41114758, 0.82608407, 0.77190533, 0.42830363, 0.57300759],
[ 0.68435626, 0.94874394, 0.55238567, 0.40367885, 0.42955156],
[ 0.59359203, 0.62237553, 0.58428725, 0.58608119, 0.29157555]])
>>> win_img[100,100]
array([[ 0.34150754, 0.17888323, 0.67222354, 0.9020784 , 0.48826682],
[ 0.68451774, 0.14887515, 0.44892615, 0.33352743, 0.22090103],
[ 0.41114758, 0.82608407, 0.77190533, 0.42830363, 0.57300759],
[ 0.68435626, 0.94874394, 0.55238567, 0.40367885, 0.42955156],
[ 0.59359203, 0.62237553, 0.58428725, 0.58608119, 0.29157555]])

Sin embargo, debes tener cuidado al no convertir tu vista en ventana de la imagen en una copia en ventana: en mi ejemplo, eso requeriría 25 veces más almacenamiento. Creo que Numpy 1.7 te permite seleccionar más de un eje, por lo que podrías simplemente hacer:

>>> np.var(win_img, axis=(-1, -2))

Estoy atorado con Numpy 1.6.2, así que no puedo probar eso. La otra opción, que puede fallar con ventanas no tan grandes, sería hacer, si recuerdo mis cálculos correctamente:

>>> win_mean = np.sum(np.sum(win_img, axis=-1), axis=-1)/win_rows/win_cols
>>> win_sqr_mean = np.sum(np.sum(win_img**2, axis=-1), axis=-1)/win_rows/win_cols
>>> win_var = win_sqr_mean - win_mean**2

Y ahora win_var es una matriz de formas

>>> win_var.shape
(496, 496)

y win_var[i, j] contiene la varianza de la ventana (5, 5) con la esquina superior izquierda en [i, j] .

Answer 2

Después de un poco de optimización creamos esta función para una imagen 3D genérica:

def variance_filter( img, VAR_FILTER_SIZE ):
from numpy.lib.stride_tricks import as_strided
WIN_SIZE=(2*VAR_FILTER_SIZE)+1
if ~ VAR_FILTER_SIZE%2==1:
print ''Warning, VAR_FILTER_SIZE must be ODD Integer number ''
# hack -- this could probably be an input to the function but Alessandro is lazy
WIN_DIMS = [ WIN_SIZE, WIN_SIZE, WIN_SIZE ]
# Check that there is a 3D image input.
if len( img.shape ) != 3:
print "/t variance_filter: Are you sure that you passed me a 3D image?"
return -1
else:
DIMS = img.shape
# Set up a windowed view on the data... this will have a border removed compared to the img_in
img_strided = as_strided(img, shape=(DIMS[0]-WIN_DIMS[0]+1, DIMS[1]-WIN_DIMS[1]+1, DIMS[2]-WIN_DIMS[2]+1, WIN_DIMS[0], WIN_DIMS[1], WIN_DIMS[2] ), strides=img.strides*2)
# Calculate variance, vectorially
win_mean = numpy.sum(numpy.sum(numpy.sum(img_strided, axis=-1), axis=-1), axis=-1) / (WIN_DIMS[0]*WIN_DIMS[1]*WIN_DIMS[2])
# As per http://en.wikipedia.org/wiki/Variance, we are removing the mean from every window,
# then squaring the result.
# Casting to 64 bit float inside, because the numbers (at least for our images) get pretty big
win_var = numpy.sum(numpy.sum(numpy.sum((( img_strided.T.astype(''<f8'') - win_mean.T.astype(''<f8'') )**2).T, axis=-1), axis=-1), axis=-1) / (WIN_DIMS[0]*WIN_DIMS[1]*WIN_DIMS[2])
# Prepare an output image of the right size, in order to replace the border removed with the windowed view call
out_img = numpy.zeros( DIMS, dtype=''<f8'' )
# copy borders out...
out_img[ WIN_DIMS[0]/2:DIMS[0]-WIN_DIMS[0]+1+WIN_DIMS[0]/2, WIN_DIMS[1]/2:DIMS[1]-WIN_DIMS[1]+1+WIN_DIMS[1]/2, WIN_DIMS[2]/2:DIMS[2]-WIN_DIMS[2]+1+WIN_DIMS[2]/2, ] = win_var
# output
return out_img.astype(''>f4'')

Answer 3

Solución más simple y también más rápida: use ndimage.uniform_filter de SciPy

import numpy as np
from scipy import ndimage
rows, cols = 500, 500
win_rows, win_cols = 5, 5
img = np.random.rand(rows, cols)
win_mean = ndimage.uniform_filter(img, (win_rows, win_cols))
win_sqr_mean = ndimage.uniform_filter(img**2, (win_rows, win_cols))
win_var = win_sqr_mean - win_mean**2

El "truco de zancada" es un hermoso truco, pero 4 más lento y no tan legible. el generic_filter es 20 veces más lento que los zancadas ...

Answer 4

Puedes usar scipy.ndimage.generic_filter . No puedo realizar pruebas con matlab, pero quizás esto te brinde lo que estás buscando:

import numpy as np
import scipy.ndimage as ndimage
subs = 10 # this is the size of the (square) sub-windows
img = np.random.rand(500, 500)
img_std = ndimage.filters.generic_filter(img, np.std, size=subs)

Puede hacer las subventanas de tamaños arbitrarios usando la palabra clave footprint . Vea esta pregunta para un ejemplo.