From 35f8ff18d48c2a84f3a3e7ce9bc99bec4da954fe Mon Sep 17 00:00:00 2001 From: Proton local user Date: Mon, 4 Nov 2024 23:12:16 +0100 Subject: [PATCH] tril performance now Tril done only for pixels above threshold. there is a mistake beacause the image is not moving consistely in the right direction when apply small shifts in sequence --- ...nJacobsRayCastInterpolateImageFunction.hxx | 133 +++++++++--------- 1 file changed, 70 insertions(+), 63 deletions(-) diff --git a/reg23Topograms/itkDTRrecon/itkgSiddonJacobsRayCastInterpolateImageFunction.hxx b/reg23Topograms/itkDTRrecon/itkgSiddonJacobsRayCastInterpolateImageFunction.hxx index 1c99d61..7c60caa 100644 --- a/reg23Topograms/itkDTRrecon/itkgSiddonJacobsRayCastInterpolateImageFunction.hxx +++ b/reg23Topograms/itkDTRrecon/itkgSiddonJacobsRayCastInterpolateImageFunction.hxx @@ -128,7 +128,7 @@ gSiddonJacobsRayCastInterpolateImageFunction::Evaluate(c float alphaX, alphaY, alphaZ, alphaCmin, alphaCminPrev; float alphaUx, alphaUy, alphaUz; float alphaIntersectionUp[3], alphaIntersectionDown[3]; - float d12, value; + float d12, value,valuetril; float firstIntersectionIndex[3]; int firstIntersectionIndexUp[3], firstIntersectionIndexDown[3]; int iU, jU, kU; @@ -468,71 +468,78 @@ gSiddonJacobsRayCastInterpolateImageFunction::Evaluate(c (cIndex[2] < static_cast(sizeCT[2]))) { // Calculate entry and exit points using alphaCmin and alphaCminPrev - - PointType entryPoint; - entryPoint[0] = SourceWorld[0] + alphaCminPrev * rayVector[0] ; - entryPoint[1] = SourceWorld[1] + alphaCminPrev * rayVector[1] ; - entryPoint[2] = SourceWorld[2] + alphaCminPrev * rayVector[2] ; - - PointType exitPoint; - exitPoint[0] = SourceWorld[0] + alphaCmin * rayVector[0] ; - exitPoint[1] = SourceWorld[1] + alphaCmin * rayVector[1] ; - exitPoint[2] = SourceWorld[2] + alphaCmin * rayVector[2] ; - -// PointType entryPoint = SourceWorld + alphaCminPrev * rayVector; -// PointType exitPoint = SourceWorld + alphaCmin * rayVector; - PointType midpoint; - midpoint[0]= (entryPoint[0] + exitPoint[0]) * 0.5; - midpoint[1]= (entryPoint[1] + exitPoint[1]) * 0.5; - midpoint[2]= (entryPoint[2] + exitPoint[2]) * 0.5; - -// = (entryPoint + exitPoint) * 0.5; - - - // Calculate fractional parts for interpolation at the midpoint - double x_frac = midpoint[0] - std::floor(midpoint[0]); - double y_frac = midpoint[1] - std::floor(midpoint[1]); - double z_frac = midpoint[2] - std::floor(midpoint[2]); - - // Perform boundary checks for trilinear interpolation - IndexType baseIndex = cIndex; - bool within_bounds = (baseIndex[0] >= 0 && baseIndex[0] < sizeCT[0] - 1) && - (baseIndex[1] >= 0 && baseIndex[1] < sizeCT[1] - 1) && - (baseIndex[2] >= 0 && baseIndex[2] < sizeCT[2] - 1); - - if (within_bounds) - { - // Fetch intensities from neighboring voxels for trilinear interpolation - float c000 = static_cast(inputPtr->GetPixel(baseIndex)); - baseIndex[0] += 1; float c100 = static_cast(inputPtr->GetPixel(baseIndex)); - baseIndex[1] += 1; float c110 = static_cast(inputPtr->GetPixel(baseIndex)); - baseIndex[0] -= 1; float c010 = static_cast(inputPtr->GetPixel(baseIndex)); - baseIndex[2] += 1; float c011 = static_cast(inputPtr->GetPixel(baseIndex)); - baseIndex[0] += 1; float c111 = static_cast(inputPtr->GetPixel(baseIndex)); - baseIndex[1] -= 1; float c101 = static_cast(inputPtr->GetPixel(baseIndex)); - baseIndex[0] -= 1; float c001 = static_cast(inputPtr->GetPixel(baseIndex)); - - // Perform trilinear interpolation at the midpoint - float c00 = c000 * (1 - x_frac) + c100 * x_frac; - float c01 = c001 * (1 - x_frac) + c101 * x_frac; - float c10 = c010 * (1 - x_frac) + c110 * x_frac; - float c11 = c011 * (1 - x_frac) + c111 * x_frac; - - float c0 = c00 * (1 - y_frac) + c10 * y_frac; - float c1 = c01 * (1 - y_frac) + c11 * y_frac; - - value = c0 * (1 - z_frac) + c1 * z_frac; - } - else - { - // Fallback to nearest voxel intensity if out-of-bounds - value = static_cast(inputPtr->GetPixel(cIndex)); - } - + value = static_cast(inputPtr->GetPixel(cIndex)); // Accumulate the interpolated intensity along the ray path if (value > m_Threshold) /* Ignore voxels whose intensities are below the threshold. */ { - d12 += (alphaCmin - alphaCminPrev) * (value - m_Threshold); + + PointType entryPoint; + entryPoint[0] = SourceWorld[0] + alphaCminPrev * rayVector[0] ; + entryPoint[1] = SourceWorld[1] + alphaCminPrev * rayVector[1] ; + entryPoint[2] = SourceWorld[2] + alphaCminPrev * rayVector[2] ; + + PointType exitPoint; + exitPoint[0] = SourceWorld[0] + alphaCmin * rayVector[0] ; + exitPoint[1] = SourceWorld[1] + alphaCmin * rayVector[1] ; + exitPoint[2] = SourceWorld[2] + alphaCmin * rayVector[2] ; + + // PointType entryPoint = SourceWorld + alphaCminPrev * rayVector; + // PointType exitPoint = SourceWorld + alphaCmin * rayVector; + PointType midpoint; + midpoint[0]= (entryPoint[0] + exitPoint[0]) * 0.5; + midpoint[1]= (entryPoint[1] + exitPoint[1]) * 0.5; + midpoint[2]= (entryPoint[2] + exitPoint[2]) * 0.5; + + // = (entryPoint + exitPoint) * 0.5; + + + // Calculate fractional parts for interpolation at the midpoint + double x_frac = midpoint[0] - std::floor(midpoint[0]); + double y_frac = midpoint[1] - std::floor(midpoint[1]); + double z_frac = midpoint[2] - std::floor(midpoint[2]); + + // Perform boundary checks for trilinear interpolation + IndexType baseIndex = cIndex; + bool within_bounds = (baseIndex[0] >= 0 && baseIndex[0] < sizeCT[0] - 1) && + (baseIndex[1] >= 0 && baseIndex[1] < sizeCT[1] - 1) && + (baseIndex[2] >= 0 && baseIndex[2] < sizeCT[2] - 1); + + if (within_bounds) + { + // Fetch intensities from neighboring voxels for trilinear interpolation + float c000 = static_cast(inputPtr->GetPixel(baseIndex)); + baseIndex[0] += 1; + float c100 = static_cast(inputPtr->GetPixel(baseIndex)); + baseIndex[1] += 1; + float c110 = static_cast(inputPtr->GetPixel(baseIndex)); + baseIndex[0] -= 1; + float c010 = static_cast(inputPtr->GetPixel(baseIndex)); + baseIndex[2] += 1; + float c011 = static_cast(inputPtr->GetPixel(baseIndex)); + baseIndex[0] += 1; + float c111 = static_cast(inputPtr->GetPixel(baseIndex)); + baseIndex[1] -= 1; + float c101 = static_cast(inputPtr->GetPixel(baseIndex)); + baseIndex[0] -= 1; + float c001 = static_cast(inputPtr->GetPixel(baseIndex)); + + // Perform trilinear interpolation at the midpoint + float c00 = c000 * (1 - x_frac) + c100 * x_frac; + float c01 = c001 * (1 - x_frac) + c101 * x_frac; + float c10 = c010 * (1 - x_frac) + c110 * x_frac; + float c11 = c011 * (1 - x_frac) + c111 * x_frac; + + float c0 = c00 * (1 - y_frac) + c10 * y_frac; + float c1 = c01 * (1 - y_frac) + c11 * y_frac; + + valuetril = c0 * (1 - z_frac) + c1 * z_frac; + + d12 += (alphaCmin - alphaCminPrev) * (valuetril - m_Threshold); + + } else { + d12 += (alphaCmin - alphaCminPrev) * (value - m_Threshold); + } + } }