fattori_g/reg23Topograms
0
0
Fork 0
Code Issues 1 Pull Requests Actions Packages Projects Releases Wiki Activity

tril performance

Browse Source 
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
This commit is contained in:
Proton local user
2024-11-04 23:12:16 +01:00
parent 8e56950e24
commit 35f8ff18d4
1 changed files with 70 additions and 63 deletions
Download Patch File Download Diff File Expand all files Collapse all files

133
reg23Topograms/itkDTRrecon/itkgSiddonJacobsRayCastInterpolateImageFunction.hxx
View File

@ -128,7 +128,7 @@ gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::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<TInputImage, TCoordRep>::Evaluate(c
(cIndex[2] < static_cast<IndexValueType>(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<float>(inputPtr->GetPixel(baseIndex));
baseIndex[0] += 1; float c100 = static_cast<float>(inputPtr->GetPixel(baseIndex));
baseIndex[1] += 1; float c110 = static_cast<float>(inputPtr->GetPixel(baseIndex));
baseIndex[0] -= 1; float c010 = static_cast<float>(inputPtr->GetPixel(baseIndex));
baseIndex[2] += 1; float c011 = static_cast<float>(inputPtr->GetPixel(baseIndex));
baseIndex[0] += 1; float c111 = static_cast<float>(inputPtr->GetPixel(baseIndex));
baseIndex[1] -= 1; float c101 = static_cast<float>(inputPtr->GetPixel(baseIndex));
baseIndex[0] -= 1; float c001 = static_cast<float>(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<float>(inputPtr->GetPixel(cIndex));
}
value = static_cast<float>(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<float>(inputPtr->GetPixel(baseIndex));
baseIndex[0] += 1;
float c100 = static_cast<float>(inputPtr->GetPixel(baseIndex));
baseIndex[1] += 1;
float c110 = static_cast<float>(inputPtr->GetPixel(baseIndex));
baseIndex[0] -= 1;
float c010 = static_cast<float>(inputPtr->GetPixel(baseIndex));
baseIndex[2] += 1;
float c011 = static_cast<float>(inputPtr->GetPixel(baseIndex));
baseIndex[0] += 1;
float c111 = static_cast<float>(inputPtr->GetPixel(baseIndex));
baseIndex[1] -= 1;
float c101 = static_cast<float>(inputPtr->GetPixel(baseIndex));
baseIndex[0] -= 1;
float c001 = static_cast<float>(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);
}
}
}