Speed up scipy ndimage measurements applied on a numpy 3-d
<p>I have multiple large labeled numpy 2d array (10 000x10 000). For each label (connected cells with the same number) I want to calculate multiple measurements based on the values of another numpy 3-d array (mean, std., max, etc.). This is possible with the scipy.ndimage.labeled_comprehension tool if the 3-d numpy is converted to a 2-d. However, as the numbers of labels and the size of the arrays is rather big, the calculations take quite some time. My current code seems redundant as I am now iterating over the same labels for each 3rd dimension of the input image. I am wondering if there are ways to speed-up my code (for example by combining the three scipy.ndimage.labeled_comprehension calculations into a single calculation).</p>

<p><em>Using a test dataset of shape (4200,3000,3) and 283047 labels the calculations took 10:34 minutes</em></p>

<hr>

<p><strong>Test data</strong></p>

<pre><code>example_labels=np.array([[1, 1, 3, 3],
   [1, 2, 2, 3],
   [2, 2, 4, 4],
   [5, 5, 5, 4]])

unique_labels=np.unique(example_labels)
value_array=np.arange(48).reshape(4,4,3)
</code></pre>

<p><strong>Current code and desired output</strong></p>

<pre><code>def mean_std_measurement(x):
    xmean = x.mean()
    xstd = x.std()
    vals.append([xmean,xstd])

def calculate_measurements(labels, unique_labels, value_arr):
    global vals
    vals=[]
    ndimage.labeled_comprehension(value_array[:,:,0],labels,unique_labels,mean_std_measurement,float,-1)
    val1=np.array(vals)
    vals=[]
    ndimage.labeled_comprehension(value_array[:,:,1],labels,unique_labels,mean_std_measurement,float,-1)
    val2=np.array(vals)
    vals=[]
    ndimage.labeled_comprehension(value_array[:,:,2],labels,unique_labels,mean_std_measurement,float,-1)
    val3=np.array(vals)
    return np.column_stack((unique_labels,val1,val2,val3))

&gt;&gt;&gt; print calculate_measurements(example_labels,unique_labels,value_array)
array([[  1.        ,   5.        ,   5.09901951,   6.        ,
      5.09901951,   7.        ,   5.09901951],
   [  2.        ,  21.        ,   4.74341649,  22.        ,
      4.74341649,  23.        ,   4.74341649],
   [  3.        ,  12.        ,   6.4807407 ,  13.        ,
      6.4807407 ,  14.        ,   6.4807407 ],
   [  4.        ,  36.        ,   6.4807407 ,  37.        ,
      6.4807407 ,  38.        ,   6.4807407 ],
   [  5.        ,  39.        ,   2.44948974,  40.        ,
      2.44948974,  41.        ,   2.44948974]])
</code></pre>

<p>The list appends in mean_std_measurement is probably forming a bottleneck when dealing with large arrays. <code>labeled_comprehension</code> will return an array on its own, so a first step is to just let scipy handle the array construction under the hood. </p>

<p><code>labeled_comprehension</code> can only apply functions that output a single value---I suspect this is why you were using the list construction in the first place---but we can cheat this by having the function output a complex value. Another option would be to use structured datatypes as the output, which would be necessary if we were returning more than 2 values.</p>

<pre><code>import numpy as np
from scipy import ndimage
example_labels=np.array([[1, 1, 3, 3],
   [1, 2, 2, 3],
   [2, 2, 4, 4],
   [5, 5, 5, 4]])

unique_labels=np.unique(example_labels)
value_array=np.arange(48).reshape(4,4,3)
# return a complex number to get around labled_comprehension limitations
def mean_std_measurement(x):
    xmean = x.mean()
    xstd = x.std()
    return np.complex(xmean, xstd)

def calculate_measurements(labels, unique_labels, value_array):
    val1 = ndimage.labeled_comprehension(value_array[:,:,0],labels,unique_labels,mean_std_measurement,np.complex,-1)
    val2 = ndimage.labeled_comprehension(value_array[:,:,1],labels,unique_labels,mean_std_measurement,np.complex,-1)
    val3 = ndimage.labeled_comprehension(value_array[:,:,2],labels,unique_labels,mean_std_measurement,np.complex,-1)
    # convert the complex numbers back into reals
    return np.column_stack((np.unique(labels),
                            val1.real,val1.imag,
                            val2.real,val2.imag,
                            val3.real,val3.imag))
</code></pre>

python numpy scipy vectorization ndimage