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
  • Questions like these are much easier to answer if you have some example code– SpangenFeb 13 at 17:28
  • I've added code. However, I'm not looking for subtle refinements on my code but instead just a direction to take for improvements– SeanFeb 13 at 17:35
  • Thanks Sean. I understand. I'm just trying to increase your chances of getting a worthwhile answer– SpangenFeb 13 at 17:36
  • What is 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

Most people use minc for this kind of thing. For your problem, you need mincresample:


It's mostly a command-line tool, but there do seem to be some py bindings.

    Your Answer


    By posting your answer, you agree to the privacy policy and terms of service.

    Not the answer you're looking for? Browse other questions tagged or ask your own question.