From 24b5ed7235932238311095653b06ffddd7709fb8 Mon Sep 17 00:00:00 2001 From: David Gobbi Date: Thu, 14 Jul 2016 11:10:28 -0600 Subject: [PATCH 1/3] Use vtkImagePointDataIterator in vtkImageReslice The vtkImagePointDataIterator provides a more efficient and sometimes more convenient way to iterate through image data. It also takes care of progress reporting. --- Imaging/Core/vtkImageReslice.cxx | 644 ++++++++++++++++--------------- 1 file changed, 324 insertions(+), 320 deletions(-) diff --git a/Imaging/Core/vtkImageReslice.cxx b/Imaging/Core/vtkImageReslice.cxx index ed73eacccae..3c8392ec1fc 100644 --- a/Imaging/Core/vtkImageReslice.cxx +++ b/Imaging/Core/vtkImageReslice.cxx @@ -24,6 +24,7 @@ #include "vtkStreamingDemandDrivenPipeline.h" #include "vtkTransform.h" #include "vtkPointData.h" +#include "vtkImagePointDataIterator.h" #include "vtkImageInterpolator.h" #include "vtkGarbageCollector.h" @@ -1883,50 +1884,6 @@ void vtkImageResliceRescaleScalars(F *floatData, int components, int n, } } -//---------------------------------------------------------------------------- -// helper function for clipping of the output with a stencil -int vtkResliceGetNextExtent(vtkImageStencilData *stencil, - int &r1, int &r2, int rmin, int rmax, - int yIdx, int zIdx, - void *&outPtr, void *background, - int numscalars, - void (*setpixels)(void *&out, - const void *in, - int numscalars, - int n), - int &iter) -{ - // trivial case if stencil is not set - if (!stencil) - { - if (iter++ == 0) - { - r1 = rmin; - r2 = rmax; - return 1; - } - return 0; - } - - // for clearing, start at last r2 plus 1 - int clear1 = r2 + 1; - if (iter == 0) - { // if no 'last time', start at rmin - clear1 = rmin; - } - - int rval = stencil->GetNextExtent(r1, r2, rmin, rmax, yIdx, zIdx, iter); - int clear2 = r1 - 1; - if (rval == 0) - { - clear2 = rmax; - } - - setpixels(outPtr, background, numscalars, clear2 - clear1 + 1); - - return rval; -} - //---------------------------------------------------------------------------- // This function simply clears the entire output to the background color, // for cases where the transformation places the output extent completely @@ -1937,12 +1894,6 @@ void vtkImageResliceClearExecute(vtkImageReslice *self, { void (*setpixels)(void *&out, const void *in, int numscalars, int n) = 0; - // for the progress meter - unsigned long count = 0; - unsigned long target = static_cast - ((outExt[5]-outExt[4]+1)*(outExt[3]-outExt[2]+1)/50.0); - target++; - // Get Increments to march through data vtkIdType outIncX, outIncY, outIncZ; outData->GetContinuousIncrements(outExt, outIncX, outIncY, outIncZ); @@ -1958,26 +1909,12 @@ void vtkImageResliceClearExecute(vtkImageReslice *self, vtkGetSetPixelsFunc(&setpixels, scalarType, scalarSize, numscalars, outPtr); - // Loop through output voxels - for (int idZ = outExt[4]; idZ <= outExt[5]; idZ++) + vtkImagePointDataIterator iter(outData, outExt, NULL, self, threadId); + for (; !iter.IsAtEnd(); iter.NextSpan()) { - for (int idY = outExt[2]; idY <= outExt[3]; idY++) - { - if (threadId == 0) - { // update the progress if this is the main thread - if (!(count%target)) - { - self->UpdateProgress(count/(50.0*target)); - } - count++; - } - // clear the pixels to background color and go to next row - setpixels(outPtr, background, numscalars, outExt[1]-outExt[0]+1); - outPtr = static_cast( - static_cast(outPtr) + outIncY*scalarSize); - } - outPtr = static_cast( - static_cast(outPtr) + outIncZ*scalarSize); + // clear the pixels to background color and go to next row + outPtr = iter.GetVoidPointer(outData, iter.GetId()); + setpixels(outPtr, background, numscalars, outExt[1]-outExt[0]+1); } vtkFreeBackgroundPixel(&background); @@ -2016,12 +1953,6 @@ void vtkImageResliceExecute(vtkImageReslice *self, void (*setpixels)(void *&out, const void *in, int numscalars, int n) = 0; void (*composite)(F *in, int numscalars, int n) = 0; - // for the progress meter - unsigned long count = 0; - unsigned long target = static_cast - ((outExt[5]-outExt[4]+1)*(outExt[3]-outExt[2]+1)/50.0); - target++; - // get the input stencil vtkImageStencilData *stencil = self->GetStencil(); // get the output stencil @@ -2088,9 +2019,7 @@ void vtkImageResliceExecute(vtkImageReslice *self, optimizeNearest = 1; } - // get Increments to march through data - vtkIdType outIncX, outIncY, outIncZ; - outData->GetContinuousIncrements(outExt, outIncX, outIncY, outIncZ); + // get pixel information int scalarType = outData->GetScalarType(); int scalarSize = outData->GetScalarSize(); int outComponents = outData->GetNumberOfScalarComponents(); @@ -2142,210 +2071,220 @@ void vtkImageResliceExecute(vtkImageReslice *self, vtkGetCompositeFunc(&composite, self->GetSlabMode(), self->GetSlabTrapezoidIntegration()); - // Loop through output pixels - for (int idZ = outExt[4]; idZ <= outExt[5]; idZ++) + // create some variables for when we march through the data + int idY = outExt[2] - 1; + int idZ = outExt[4] - 1; + F inPoint0[4] = { 0.0, 0.0, 0.0, 0.0 }; + F inPoint1[4] = { 0.0, 0.0, 0.0, 0.0 }; + + // create an iterator to march through the data + vtkImagePointDataIterator iter(outData, outExt, stencil, self, threadId); + char *outPtr0 = static_cast(iter.GetVoidPointer(outData)); + for (; !iter.IsAtEnd(); iter.NextSpan()) { - F inPoint0[4]; - inPoint0[0] = origin[0] + idZ*zAxis[0]; // incremental transform - inPoint0[1] = origin[1] + idZ*zAxis[1]; - inPoint0[2] = origin[2] + idZ*zAxis[2]; - inPoint0[3] = origin[3] + idZ*zAxis[3]; + int span = static_cast(iter.SpanEndId() - iter.GetId()); + outPtr = outPtr0 + iter.GetId()*scalarSize*outComponents; - for (int idY = outExt[2]; idY <= outExt[3]; idY++) + if (!iter.IsInStencil()) + { + // clear any regions that are outside the stencil + setpixels(outPtr, background, outComponents, span); + } + else { - F inPoint1[4]; - inPoint1[0] = inPoint0[0] + idY*yAxis[0]; // incremental transform - inPoint1[1] = inPoint0[1] + idY*yAxis[1]; - inPoint1[2] = inPoint0[2] + idY*yAxis[2]; - inPoint1[3] = inPoint0[3] + idY*yAxis[3]; + // get output index, and compute position in input image + int outIndex[3]; + iter.GetIndex(outIndex); - if (!threadId) + // if Z index increased, then advance position along Z axis + if (outIndex[2] > idZ) { - if (!(count%target)) - { - self->UpdateProgress(count/(50.0*target)); - } - count++; + idZ = outIndex[2]; + inPoint0[0] = origin[0] + idZ*zAxis[0]; + inPoint0[1] = origin[1] + idZ*zAxis[1]; + inPoint0[2] = origin[2] + idZ*zAxis[2]; + inPoint0[3] = origin[3] + idZ*zAxis[3]; + idY = outExt[2] - 1; } - int iter = 0; - int idXmin, idXmax; - while (vtkResliceGetNextExtent(stencil, idXmin, idXmax, - outExt[0], outExt[1], idY, idZ, - outPtr, background, outComponents, - setpixels, iter)) + // if Y index increased, then advance position along Z axis + if (outIndex[1] > idY) { - if (!optimizeNearest) - { - bool wasInBounds = 1; - bool isInBounds = 1; - int startIdX = idXmin; - int idX = idXmin; - F *tmpPtr = floatPtr; + idY = outIndex[1]; + inPoint1[0] = inPoint0[0] + idY*yAxis[0]; + inPoint1[1] = inPoint0[1] + idY*yAxis[1]; + inPoint1[2] = inPoint0[2] + idY*yAxis[2]; + inPoint1[3] = inPoint0[3] + idY*yAxis[3]; + } - while (startIdX <= idXmax) + // march through one row of the output image + int idXmin = outIndex[0]; + int idXmax = idXmin + span - 1; + + if (!optimizeNearest) + { + bool wasInBounds = 1; + bool isInBounds = 1; + int startIdX = idXmin; + int idX = idXmin; + F *tmpPtr = floatPtr; + + while (startIdX <= idXmax) + { + for (; idX <= idXmax && isInBounds == wasInBounds; idX++) { - for (; idX <= idXmax && isInBounds == wasInBounds; idX++) + F inPoint2[4]; + inPoint2[0] = inPoint1[0] + idX*xAxis[0]; + inPoint2[1] = inPoint1[1] + idX*xAxis[1]; + inPoint2[2] = inPoint1[2] + idX*xAxis[2]; + inPoint2[3] = inPoint1[3] + idX*xAxis[3]; + + F inPoint3[4]; + F *inPoint = inPoint2; + isInBounds = 0; + + int sampleCount = 0; + for (int sample = 0; sample < nsamples; sample++) { - F inPoint2[4]; - inPoint2[0] = inPoint1[0] + idX*xAxis[0]; - inPoint2[1] = inPoint1[1] + idX*xAxis[1]; - inPoint2[2] = inPoint1[2] + idX*xAxis[2]; - inPoint2[3] = inPoint1[3] + idX*xAxis[3]; - - F inPoint3[4]; - F *inPoint = inPoint2; - isInBounds = 0; - - int sampleCount = 0; - for (int sample = 0; sample < nsamples; sample++) + if (nsamples > 1) { - if (nsamples > 1) - { - double s = sample - 0.5*(nsamples - 1); - s *= slabSampleSpacing; - inPoint3[0] = inPoint2[0] + s*zAxis[0]; - inPoint3[1] = inPoint2[1] + s*zAxis[1]; - inPoint3[2] = inPoint2[2] + s*zAxis[2]; - inPoint3[3] = inPoint2[3] + s*zAxis[3]; - inPoint = inPoint3; - } - - if (perspective) - { // only do perspective if necessary - F f = 1/inPoint[3]; - inPoint[0] *= f; - inPoint[1] *= f; - inPoint[2] *= f; - } - - if (newtrans) - { // apply the AbstractTransform if there is one - vtkResliceApplyTransform(newtrans, inPoint, inOrigin, - inInvSpacing); - } - - if (interpolator->CheckBoundsIJK(inPoint)) - { - // do the interpolation - sampleCount++; - isInBounds = 1; - interpolator->InterpolateIJK(inPoint, tmpPtr); - tmpPtr += inComponents; - } + double s = sample - 0.5*(nsamples - 1); + s *= slabSampleSpacing; + inPoint3[0] = inPoint2[0] + s*zAxis[0]; + inPoint3[1] = inPoint2[1] + s*zAxis[1]; + inPoint3[2] = inPoint2[2] + s*zAxis[2]; + inPoint3[3] = inPoint2[3] + s*zAxis[3]; + inPoint = inPoint3; + } + + if (perspective) + { // only do perspective if necessary + F f = 1/inPoint[3]; + inPoint[0] *= f; + inPoint[1] *= f; + inPoint[2] *= f; + } + + if (newtrans) + { // apply the AbstractTransform if there is one + vtkResliceApplyTransform(newtrans, inPoint, inOrigin, + inInvSpacing); } - tmpPtr -= sampleCount*inComponents; - if (sampleCount > 1) + if (interpolator->CheckBoundsIJK(inPoint)) { - composite(tmpPtr, inComponents, sampleCount); + // do the interpolation + sampleCount++; + isInBounds = 1; + interpolator->InterpolateIJK(inPoint, tmpPtr); + tmpPtr += inComponents; } - tmpPtr += inComponents; + } - // set "was in" to "is in" if first pixel - wasInBounds = ((idX > idXmin) ? wasInBounds : isInBounds); + tmpPtr -= sampleCount*inComponents; + if (sampleCount > 1) + { + composite(tmpPtr, inComponents, sampleCount); } + tmpPtr += inComponents; - // write a segment to the output - int endIdX = idX - 1 - (isInBounds != wasInBounds); - int numpixels = endIdX - startIdX + 1; + // set "was in" to "is in" if first pixel + wasInBounds = ((idX > idXmin) ? wasInBounds : isInBounds); + } + + // write a segment to the output + int endIdX = idX - 1 - (isInBounds != wasInBounds); + int numpixels = endIdX - startIdX + 1; - if (wasInBounds) + if (wasInBounds) + { + if (outputStencil) { - if (outputStencil) - { - outputStencil->InsertNextExtent(startIdX, endIdX, idY, idZ); - } + outputStencil->InsertNextExtent(startIdX, endIdX, idY, idZ); + } - if (rescaleScalars) - { - vtkImageResliceRescaleScalars(floatPtr, inComponents, - idXmax - idXmin + 1, - scalarShift, scalarScale); - } + if (rescaleScalars) + { + vtkImageResliceRescaleScalars(floatPtr, inComponents, + idXmax - idXmin + 1, + scalarShift, scalarScale); + } - if (convertScalars) - { - (self->*convertScalars)(tmpPtr - inComponents*(idX-startIdX), - outPtr, - vtkTypeTraits::VTKTypeID(), - inComponents, numpixels, - startIdX, idY, idZ, threadId); - - outPtr = static_cast(static_cast(outPtr) - + numpixels*outComponents*scalarSize); - } - else - { - convertpixels(outPtr, tmpPtr - inComponents*(idX - startIdX), - outComponents, numpixels); - } + if (convertScalars) + { + (self->*convertScalars)(tmpPtr - inComponents*(idX-startIdX), + outPtr, + vtkTypeTraits::VTKTypeID(), + inComponents, numpixels, + startIdX, idY, idZ, threadId); + + outPtr = static_cast(outPtr) + + numpixels*outComponents*scalarSize; } else { - setpixels(outPtr, background, outComponents, numpixels); + convertpixels(outPtr, tmpPtr - inComponents*(idX - startIdX), + outComponents, numpixels); } - - startIdX += numpixels; - wasInBounds = isInBounds; } - } - else // optimize for nearest-neighbor interpolation - { - char *outPtrTmp = static_cast(outPtr); + else + { + setpixels(outPtr, background, outComponents, numpixels); + } - int inExtX = inExt[1] - inExt[0] + 1; - int inExtY = inExt[3] - inExt[2] + 1; - int inExtZ = inExt[5] - inExt[4] + 1; + startIdX += numpixels; + wasInBounds = isInBounds; + } + } + else // optimize for nearest-neighbor interpolation + { + char *outPtrTmp = static_cast(outPtr); - int startIdX = idXmin; - int endIdX = idXmin-1; - bool isInBounds = false; + int inExtX = inExt[1] - inExt[0] + 1; + int inExtY = inExt[3] - inExt[2] + 1; + int inExtZ = inExt[5] - inExt[4] + 1; - for (int iidX = idXmin; iidX <= idXmax; iidX++) - { - char *inPtrTmp = static_cast(background); - int bytesPerPixel = inputScalarSize*inComponents; + int startIdX = idXmin; + int endIdX = idXmin-1; + bool isInBounds = false; - F inPoint[3]; - inPoint[0] = inPoint1[0] + iidX*xAxis[0]; - inPoint[1] = inPoint1[1] + iidX*xAxis[1]; - inPoint[2] = inPoint1[2] + iidX*xAxis[2]; + for (int iidX = idXmin; iidX <= idXmax; iidX++) + { + char *inPtrTmp = static_cast(background); + int bytesPerPixel = inputScalarSize*inComponents; - int inIdX = vtkInterpolationMath::Round(inPoint[0]) - inExt[0]; - int inIdY = vtkInterpolationMath::Round(inPoint[1]) - inExt[2]; - int inIdZ = vtkInterpolationMath::Round(inPoint[2]) - inExt[4]; + F inPoint[3]; + inPoint[0] = inPoint1[0] + iidX*xAxis[0]; + inPoint[1] = inPoint1[1] + iidX*xAxis[1]; + inPoint[2] = inPoint1[2] + iidX*xAxis[2]; - if ((inIdX >= 0) & (inIdX < inExtX) & - (inIdY >= 0) & (inIdY < inExtY) & - (inIdZ >= 0) & (inIdZ < inExtZ)) - { - inPtrTmp = static_cast(inPtr) + - (inIdX*inInc[0] + inIdY*inInc[1] + inIdZ*inInc[2])* - inputScalarSize; + int inIdX = vtkInterpolationMath::Round(inPoint[0]) - inExt[0]; + int inIdY = vtkInterpolationMath::Round(inPoint[1]) - inExt[2]; + int inIdZ = vtkInterpolationMath::Round(inPoint[2]) - inExt[4]; - startIdX = (isInBounds ? startIdX : iidX); - endIdX = iidX; - isInBounds = true; - } + if ((inIdX >= 0) & (inIdX < inExtX) & + (inIdY >= 0) & (inIdY < inExtY) & + (inIdZ >= 0) & (inIdZ < inExtZ)) + { + inPtrTmp = static_cast(inPtr) + + (inIdX*inInc[0] + inIdY*inInc[1] + inIdZ*inInc[2])* + inputScalarSize; - int oc = bytesPerPixel; - do { *outPtrTmp++ = *inPtrTmp++; } while (--oc); + startIdX = (isInBounds ? startIdX : iidX); + endIdX = iidX; + isInBounds = true; } - outPtr = outPtrTmp; + int oc = bytesPerPixel; + do { *outPtrTmp++ = *inPtrTmp++; } while (--oc); + } - if (outputStencil && endIdX >= startIdX) - { - outputStencil->InsertNextExtent(startIdX, endIdX, idY, idZ); - } + if (outputStencil && endIdX >= startIdX) + { + outputStencil->InsertNextExtent(startIdX, endIdX, idY, idZ); } } - outPtr = static_cast( - static_cast(outPtr) + outIncY*scalarSize); } - outPtr = static_cast( - static_cast(outPtr) + outIncZ*scalarSize); } vtkFreeBackgroundPixel(&background); @@ -2888,115 +2827,180 @@ void vtkReslicePermuteExecute(vtkImageReslice *self, vtkAllocBackgroundPixel(&background, self->GetBackgroundColor(), scalarType, scalarSize, outComponents); - // for tracking progress - unsigned long count = 0; - unsigned long target = static_cast - ((outExt[5]-outExt[4]+1)*(outExt[3]-outExt[2]+1)/50.0); - target++; + // generate the extent we will iterate over while painting output + // voxels with input data (anything outside will be background color) + int iterExt[6]; + bool empty = false; + for (int jj = 0; jj < 6; jj += 2) + { + iterExt[jj] = clipExt[jj]; + iterExt[jj+1] = clipExt[jj+1]; + empty |= (iterExt[jj] > iterExt[jj+1]); + } + if (empty) + { + for (int jj = 0; jj < 6; jj += 2) + { + iterExt[jj] = outExt[jj]; + iterExt[jj+1] = outExt[jj]-1; + } + } + else if (nsamples > 1) + { + // adjust extent for multiple samples if slab mode + int adjust = nsamples - 1; + int maxAdjustDown = iterExt[4] - outExt[4]; + iterExt[4] -= (adjust <= maxAdjustDown ? adjust : maxAdjustDown); + int maxAdjustUp = outExt[5] - iterExt[5]; + iterExt[5] += (adjust <= maxAdjustUp ? adjust : maxAdjustUp); + } - // Loop through output pixels - for (int idZ = outExt[4]; idZ <= outExt[5]; idZ++) + // clear any leading slices + for (int idZ = outExt[4]; idZ < iterExt[4]; idZ++) { + int fullspan = outExt[1] - outExt[0] + 1; for (int idY = outExt[2]; idY <= outExt[3]; idY++) { - if (threadId == 0) - { // track progress if we are main thread - if (!(count%target)) + setpixels(outPtr, background, outComponents, fullspan); + outPtr = static_cast(outPtr) + outIncY*scalarSize; + } + outPtr = static_cast(outPtr) + outIncZ*scalarSize; + } + + if (!empty) + { + vtkImagePointDataIterator iter(outData, iterExt, stencil, self, threadId); + for (; !iter.IsAtEnd(); iter.NextSpan()) + { + // get output index + int outIndex[3]; + iter.GetIndex(outIndex); + int span = static_cast(iter.SpanEndId() - iter.GetId()); + int idXmin = outIndex[0]; + int idXmax = idXmin + span - 1; + int idY = outIndex[1]; + int idZ = outIndex[2]; + + if (idXmin == iterExt[0]) + { + // clear rows that were outside of iterExt + if (idY == iterExt[2]) { - self->UpdateProgress(count/(50.0*target)); + int fullspan = outExt[1] - outExt[0] + 1; + for (idY = outExt[2]; idY < iterExt[2]; idY++) + { + setpixels(outPtr, background, outComponents, fullspan); + outPtr = static_cast(outPtr) + outIncY*scalarSize; + } + } + // clear leading pixels + if (iterExt[0] > outExt[0]) + { + setpixels(outPtr, background, outComponents, iterExt[0] - outExt[0]); } - count++; } - // do extent check - if (idZ < clipExt[4]-(nsamples-1) || idZ > clipExt[5]+(nsamples-1) || - idY < clipExt[2] || idY > clipExt[3]) - { // just clear, we're completely outside - setpixels(outPtr, background, outComponents, outExt[1]-outExt[0]+1); + if (!iter.IsInStencil()) + { + // clear any regions that are outside the stencil + setpixels(outPtr, background, outComponents, span); } else { - // clear pixels to left of input extent - setpixels(outPtr, background, outComponents, clipExt[0]-outExt[0]); - - int iter = 0; - int idXmin, idXmax; - while (vtkResliceGetNextExtent(stencil, idXmin, idXmax, - clipExt[0], clipExt[1], idY, idZ, - outPtr, background, outComponents, - setpixels, iter)) - { - int idX = idXmin; + int idX = idXmin; - if (doConversion) + if (doConversion) + { + // these six lines are for handling incomplete slabs + int lowerSkip = clipExt[4] - idZ; + lowerSkip = (lowerSkip >= 0 ? lowerSkip : 0); + int upperSkip = idZ + (nsamples-1) - clipExt[5]; + upperSkip = (upperSkip >= 0 ? upperSkip : 0); + int idZ1 = idZ + lowerSkip; + int nsamples1 = nsamples - lowerSkip - upperSkip; + + for (int isample = 0; isample < nsamples1; isample++) { - // these six lines are for handling incomplete slabs - int lowerSkip = clipExt[4] - idZ; - lowerSkip = (lowerSkip >= 0 ? lowerSkip : 0); - int upperSkip = idZ + (nsamples-1) - clipExt[5]; - upperSkip = (upperSkip >= 0 ? upperSkip : 0); - int idZ1 = idZ + lowerSkip; - int nsamples1 = nsamples - lowerSkip - upperSkip; - - for (int isample = 0; isample < nsamples1; isample++) - { - F *tmpPtr = ((nsamples1 > 1) ? floatSumPtr : floatPtr); - interpolator->InterpolateRow( - weights, idX, idY, idZ1, tmpPtr, idXmax - idXmin + 1); - - if (composite && (nsamples1 > 1)) - { - composite(floatPtr, floatSumPtr, inComponents, - idXmax - idXmin + 1, isample, nsamples1); - } + F *tmpPtr = ((nsamples1 > 1) ? floatSumPtr : floatPtr); + interpolator->InterpolateRow( + weights, idX, idY, idZ1, tmpPtr, span); - idZ1++; - } - - if (rescaleScalars) + if (composite && (nsamples1 > 1)) { - vtkImageResliceRescaleScalars(floatPtr, inComponents, - idXmax - idXmin + 1, - scalarShift, scalarScale); + composite(floatPtr, floatSumPtr, inComponents, span, + isample, nsamples1); } - if (convertScalars) - { - (self->*convertScalars)(floatPtr, outPtr, - vtkTypeTraits::VTKTypeID(), - inComponents, idXmax - idXmin + 1, - idXmin, idY, idZ, threadId); + idZ1++; + } - outPtr = static_cast(static_cast(outPtr) - + (idXmax-idXmin+1)*outComponents*scalarSize); - } - else - { - conversion(outPtr, floatPtr, inComponents, idXmax - idXmin + 1); - } + if (rescaleScalars) + { + vtkImageResliceRescaleScalars(floatPtr, inComponents, span, + scalarShift, scalarScale); + } + + if (convertScalars) + { + (self->*convertScalars)(floatPtr, outPtr, + vtkTypeTraits::VTKTypeID(), + inComponents, span, + idXmin, idY, idZ, threadId); + + outPtr = (static_cast(outPtr) + + static_cast(span)*outComponents*scalarSize); } else { - // fast path for when no conversion is necessary - summation(outPtr, idX, idY, idZ, inComponents, - idXmax - idXmin + 1, weights); + conversion(outPtr, floatPtr, inComponents, span); } + } + else + { + // fast path for when no conversion is necessary + summation(outPtr, idX, idY, idZ, inComponents, span, weights); + } + + if (outputStencil) + { + outputStencil->InsertNextExtent(idXmin, idXmax, idY, idZ); + } + } + + if (idXmax == iterExt[1]) + { + // clear trailing pixels + if (iterExt[1] < outExt[1]) + { + setpixels(outPtr, background, outComponents, outExt[1] - iterExt[1]); + } + outPtr = static_cast(outPtr) + outIncY*scalarSize; - if (outputStencil) + // clear trailing rows + if (idY == iterExt[3]) + { + int fullspan = outExt[1] - outExt[0] + 1; + for (idY = iterExt[3]+1; idY <= outExt[3]; idY++) { - outputStencil->InsertNextExtent(idXmin, idXmax, idY, idZ); + setpixels(outPtr, background, outComponents, fullspan); + outPtr = static_cast(outPtr) + outIncY*scalarSize; } + outPtr = static_cast(outPtr) + outIncZ*scalarSize; } - - // clear pixels to right of input extent - setpixels(outPtr, background, outComponents, outExt[1] - clipExt[1]); } + } + } - outPtr = static_cast( - static_cast(outPtr) + outIncY*scalarSize); + // clear any trailing slices + for (int idZ = iterExt[5]+1; idZ <= outExt[5]; idZ++) + { + int fullspan = outExt[1] - outExt[0] + 1; + for (int idY = outExt[2]; idY <= outExt[3]; idY++) + { + setpixels(outPtr, background, outComponents, fullspan); + outPtr = static_cast(outPtr) + outIncY*scalarSize; } - outPtr = static_cast( - static_cast(outPtr) + outIncZ*scalarSize); + outPtr = static_cast(outPtr) + outIncZ*scalarSize; } vtkFreeBackgroundPixel(&background); -- GitLab From 41fdbbf64b2e9efd396f738813a8e7ee7570a624 Mon Sep 17 00:00:00 2001 From: David Gobbi Date: Fri, 15 Jul 2016 06:47:45 -0600 Subject: [PATCH 2/3] Use logical and in reslice nearest-neighbor path. For this path, the resulting speed increase can be 30% or more. --- Imaging/Core/vtkImageReslice.cxx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Imaging/Core/vtkImageReslice.cxx b/Imaging/Core/vtkImageReslice.cxx index 3c8392ec1fc..5bbd793a389 100644 --- a/Imaging/Core/vtkImageReslice.cxx +++ b/Imaging/Core/vtkImageReslice.cxx @@ -2262,9 +2262,9 @@ void vtkImageResliceExecute(vtkImageReslice *self, int inIdY = vtkInterpolationMath::Round(inPoint[1]) - inExt[2]; int inIdZ = vtkInterpolationMath::Round(inPoint[2]) - inExt[4]; - if ((inIdX >= 0) & (inIdX < inExtX) & - (inIdY >= 0) & (inIdY < inExtY) & - (inIdZ >= 0) & (inIdZ < inExtZ)) + if (inIdX >= 0 && inIdX < inExtX && + inIdY >= 0 && inIdY < inExtY && + inIdZ >= 0 && inIdZ < inExtZ) { inPtrTmp = static_cast(inPtr) + (inIdX*inInc[0] + inIdY*inInc[1] + inIdZ*inInc[2])* -- GitLab From 4cc6d5113b056ac64fe8786b25951af9f2586ac1 Mon Sep 17 00:00:00 2001 From: David Gobbi Date: Fri, 15 Jul 2016 09:45:57 -0600 Subject: [PATCH 3/3] Optimize the reslice nearest neighbor interpolation Use the efficient "setpixels" method for setting out-of-bounds regions to the background color. Also, use a switch() to unroll the inner loop for common pixels sizes, this provides a performance boost of a few %. --- Imaging/Core/vtkImageReslice.cxx | 83 +++++++++++++++++++++++++++----- 1 file changed, 71 insertions(+), 12 deletions(-) diff --git a/Imaging/Core/vtkImageReslice.cxx b/Imaging/Core/vtkImageReslice.cxx index 5bbd793a389..d16c063f1c3 100644 --- a/Imaging/Core/vtkImageReslice.cxx +++ b/Imaging/Core/vtkImageReslice.cxx @@ -2238,8 +2238,13 @@ void vtkImageResliceExecute(vtkImageReslice *self, } else // optimize for nearest-neighbor interpolation { + const char *inPtrTmp0 = static_cast(inPtr); char *outPtrTmp = static_cast(outPtr); + vtkIdType inIncX = inInc[0]*inputScalarSize; + vtkIdType inIncY = inInc[1]*inputScalarSize; + vtkIdType inIncZ = inInc[2]*inputScalarSize; + int inExtX = inExt[1] - inExt[0] + 1; int inExtY = inExt[3] - inExt[2] + 1; int inExtZ = inExt[5] - inExt[4] + 1; @@ -2247,12 +2252,10 @@ void vtkImageResliceExecute(vtkImageReslice *self, int startIdX = idXmin; int endIdX = idXmin-1; bool isInBounds = false; + int bytesPerPixel = inputScalarSize*inComponents; for (int iidX = idXmin; iidX <= idXmax; iidX++) { - char *inPtrTmp = static_cast(background); - int bytesPerPixel = inputScalarSize*inComponents; - F inPoint[3]; inPoint[0] = inPoint1[0] + iidX*xAxis[0]; inPoint[1] = inPoint1[1] + iidX*xAxis[1]; @@ -2266,19 +2269,75 @@ void vtkImageResliceExecute(vtkImageReslice *self, inIdY >= 0 && inIdY < inExtY && inIdZ >= 0 && inIdZ < inExtZ) { - inPtrTmp = static_cast(inPtr) + - (inIdX*inInc[0] + inIdY*inInc[1] + inIdZ*inInc[2])* - inputScalarSize; - - startIdX = (isInBounds ? startIdX : iidX); + if (!isInBounds) + { + // clear leading out-of-bounds pixels + startIdX = iidX; + isInBounds = true; + setpixels(outPtr, background, outComponents, startIdX-idXmin); + outPtrTmp = static_cast(outPtr); + } + // set the final index that was within input bounds endIdX = iidX; - isInBounds = true; - } - int oc = bytesPerPixel; - do { *outPtrTmp++ = *inPtrTmp++; } while (--oc); + // perform nearest-neighbor interpolation via pixel copy + const char *inPtrTmp = inPtrTmp0 + + inIdX*inIncX + inIdY*inIncY + inIdZ*inIncZ; + + // manual loop unrolling significantly boosts performance, + // and is much less bloat than templating over the type + switch (bytesPerPixel) + { + case 1: + outPtrTmp[0] = inPtrTmp[0]; + break; + case 2: + outPtrTmp[0] = inPtrTmp[0]; + outPtrTmp[1] = inPtrTmp[1]; + break; + case 3: + outPtrTmp[0] = inPtrTmp[0]; + outPtrTmp[1] = inPtrTmp[1]; + outPtrTmp[2] = inPtrTmp[2]; + break; + case 4: + outPtrTmp[0] = inPtrTmp[0]; + outPtrTmp[1] = inPtrTmp[1]; + outPtrTmp[2] = inPtrTmp[2]; + outPtrTmp[3] = inPtrTmp[3]; + break; + case 8: + outPtrTmp[0] = inPtrTmp[0]; + outPtrTmp[1] = inPtrTmp[1]; + outPtrTmp[2] = inPtrTmp[2]; + outPtrTmp[3] = inPtrTmp[3]; + outPtrTmp[4] = inPtrTmp[4]; + outPtrTmp[5] = inPtrTmp[5]; + outPtrTmp[6] = inPtrTmp[6]; + outPtrTmp[7] = inPtrTmp[7]; + break; + default: + int oc = 0; + do + { + outPtrTmp[oc] = inPtrTmp[oc]; + } + while (++oc != bytesPerPixel); + break; + } + outPtrTmp += bytesPerPixel; + } + else if (isInBounds) + { + // leaving input bounds + break; + } } + // clear trailing out-of-bounds pixels + outPtr = outPtrTmp; + setpixels(outPtr, background, outComponents, idXmax-endIdX); + if (outputStencil && endIdX >= startIdX) { outputStencil->InsertNextExtent(startIdX, endIdX, idY, idZ); -- GitLab