I am working with DICOM images in python, which are series of 2-d images (512x512x160 could be typical) representing a CT scan. I am working on an algorithm which downsamples in the z-axis by an integer ratio, producing an image like 512x512x80 (or 512x512x40 etc.).

As an example, going from z=160 to z=80 means that origimg[x,y,z] and origimg[x,y,z+1] will be used to calculate the value for newimg[x,y,z]. The trick is that I have a weighting function which (linearly) defines a mapping between pixel intensity and weight. The resulting value of newimg[x,y,z] is a linear combination of origimg[x,y,z], origimg[x,y,z+1]... origimg[x,y,z+n] (where n is the downsampling ratio), but the value of the weights is dependent on the values of origimg[x,y,z] etc.

I have a naive implementation of this algorithm (no threading), but the runtime is exceedingly slow (10+ minutes), where I need it to be much faster (hopefully <15sec if not better). I am wondering if anyone would have any suggestions on how to do this in the most efficient way possible. Should I write this in C and wrap the code in Python? Will stock Python threading be enough (I am in 2.7 but could port to 3)? Is there some built-in in ITK or scikit that might help?

Edit: Including code

`r=int(newThickness/self.sliceThickness)newdarray=np.zeros([int(ceil(self.numSlices/r)), self.currMatrixSize, self.currMatrixSize])ind=0for i in xrange(newdarray.shape[0]):for j in xrange(newdarray.shape[1]):for k in xrange(newdarray.shape[2]):weights=[]Wt=0for x in range(r):if ind+x < self.currdarray.shape[0]:v=self.currdarray[ind+x,j,k]w=self.getWeightValue(v)Wt +=wweights=weights + (w*v)newdarray[i,j,k]=sum(weights)/Wtind +=r`

`self.getWeightValue`

? Is there a chance you could rewrite it to take a full 3D array as input and produce a full 3D array as output? Other than that one, I don't see a reason to loop over`i`

,`j`

, and`k`

. Numpy has indexing and element-wise operators that would do exactly what you do here. Rewrite`weights=weights + (w*v)`

as`weights +=w*v`

, and it becomes clear you only need a loop over`x`

.– Cris LuengoFeb 13 at 19:26