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):for j in xrange(newdarray.shape):for k in xrange(newdarray.shape):weights=Wt=0for x in range(r):if ind+x < self.currdarray.shape:v=self.currdarray[ind+x,j,k]w=self.getWeightValue(v)Wt +=wweights=weights + (w*v)newdarray[i,j,k]=sum(weights)/Wtind +=r