iProjectorJoseph.cc

Go to the documentation of this file.

/*
  Implementation of class iProjectorJoseph

  - separators: done
  - doxygen: done
  - default initialization: done
  - CASTOR_DEBUG: 
  - CASTOR_VERBOSE: 
*/

#include "iProjectorJoseph.hh"
#include "sOutputManager.hh"

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

iProjectorJoseph::iProjectorJoseph() : vProjector()
{
  // This projector is not compatible with SPECT attenuation correction because
  // the voxels contributing to the line are not strictly ordered with respect to
  // their distance to point2 (due to interpolations at each plane crossed)
  m_compatibleWithSPECTAttenuationCorrection = false;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

iProjectorJoseph::~iProjectorJoseph()
{
  if ( mp_boundariesX )
  {
    delete[] mp_boundariesX;
    mp_boundariesX = nullptr;
  }

  if ( mp_boundariesY )
  {
    delete[] mp_boundariesY;
    mp_boundariesY = nullptr;
  }

  if ( mp_boundariesZ )
  {
    delete[] mp_boundariesZ;
    mp_boundariesZ = nullptr;
  }
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

int iProjectorJoseph::ReadConfigurationFile(const string& a_configurationFile)
{
  // No options for joseph
  ;
  // Normal end
  return 0;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

int iProjectorJoseph::ReadOptionsList(const string& a_optionsList)
{
  // No options for joseph
  ;
  // Normal end
  return 0;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

void iProjectorJoseph::ShowHelpSpecific()
{
  cout << "This projector is a line projector that uses linear interpolation between pixels." << endl;
  cout << "It is implemented from the following published paper:" << endl;
  cout << "P. M. Joseph, \"An improved algorithm for reprojecting rays through pixel images\", IEEE Trans. Med. Imaging, vol. 1, pp. 192-6, 1982." << endl;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

int iProjectorJoseph::CheckSpecificParameters()
{
  // Nothing to check for this projector
  ;
  // Normal end
  return 0;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

int iProjectorJoseph::InitializeSpecific()
{
  // Verbose
  if (m_verbose>=2) Cout("iProjectorJoseph::InitializeSpecific() -> Use Joseph projector" << endl);
  // Allocate and compute boundaries
  mp_boundariesX = new FLTNB[ mp_nbVox[ 0 ] ];
  for( int i = 0; i < mp_nbVox[ 0 ]; ++i ) mp_boundariesX[ i ] = -mp_halfFOV[ 0 ] + i * mp_sizeVox[ 0 ];
  mp_boundariesY = new FLTNB[ mp_nbVox[ 1 ] ];
  for( int i = 0; i < mp_nbVox[ 1 ]; ++i ) mp_boundariesY[ i ] = -mp_halfFOV[ 1 ] + i * mp_sizeVox[ 1 ];
  mp_boundariesZ = new FLTNB[ mp_nbVox[ 2 ] ];
  for( int i = 0; i < mp_nbVox[ 2 ]; ++i ) mp_boundariesZ[ i ] = -mp_halfFOV[ 2 ] + i * mp_sizeVox[ 2 ];
  // Normal end
  return 0;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

INTNB iProjectorJoseph::EstimateMaxNumberOfVoxelsPerLine()
{
  // Find the maximum number of voxels along a given dimension
  INTNB max_nb_voxels_in_dimension = mp_ImageDimensionsAndQuantification->GetNbVoxX();
  if (mp_ImageDimensionsAndQuantification->GetNbVoxY()>max_nb_voxels_in_dimension) max_nb_voxels_in_dimension = mp_ImageDimensionsAndQuantification->GetNbVoxY();
  if (mp_ImageDimensionsAndQuantification->GetNbVoxZ()>max_nb_voxels_in_dimension) max_nb_voxels_in_dimension = mp_ImageDimensionsAndQuantification->GetNbVoxZ();
  // We should have at most 4 voxels in a given plane, so multiply by 4
  max_nb_voxels_in_dimension *= 4;
  // Return the value
  return max_nb_voxels_in_dimension;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

int iProjectorJoseph::ProjectWithoutTOF(int a_direction, oProjectionLine* ap_ProjectionLine )
{
  #ifdef CASTOR_DEBUG
  if (!m_initialized)
  {
    Cerr("***** iProjectorJoseph::ProjectWithoutTOF() -> Called while not initialized !" << endl);
    Exit(EXIT_DEBUG);
  }
  #endif

  #ifdef CASTOR_VERBOSE
  if (m_verbose>=10)
  {
    string direction = "";
    if (a_direction==FORWARD) direction = "forward";
    else direction = "backward";
    Cout("iProjectorJoseph::Project without TOF -> Project line '" << ap_ProjectionLine << "' in " << direction << " direction" << endl);
  }
  #endif

  // Simpler way now, hopefully it works
  FLTNB* event1 = ap_ProjectionLine->GetPosition1();
  FLTNB* event2 = ap_ProjectionLine->GetPosition2();

  // **************************************
  // STEP 1: LOR length calculation
  // **************************************

  // Distance between point event1 and event2
  FLTNB const dX21( event2[ 0 ] - event1[ 0 ] );
  FLTNB const dY21( event2[ 1 ] - event1[ 1 ] );
  FLTNB const dZ21( event2[ 2 ] - event1[ 2 ] );

  // Square of distance
  FLTNB const dX21_2( dX21 * dX21 );
  FLTNB const dY21_2( dY21 * dY21 );
  FLTNB const dZ21_2( dZ21 * dZ21 );

  // Size of the voxels
  FLTNB const szImg_X( mp_sizeVox[ 0 ] );
  FLTNB const szImg_Y( mp_sizeVox[ 1 ] );
  FLTNB const szImg_Z( mp_sizeVox[ 2 ] );

  // Get the number of the planes
  INTNB const nPlane_X( mp_nbVox[ 0 ] );
  INTNB const nPlane_Y( mp_nbVox[ 1 ] );
  INTNB const nPlane_Z( mp_nbVox[ 2 ] );
  INTNB const nPlane_XY( nPlane_X * nPlane_Y );

  // Coordinates of the first and last planes
  FLTNB const xPlane_0 = -mp_halfFOV[ 0 ];
  FLTNB const yPlane_0 = -mp_halfFOV[ 1 ];
  FLTNB const zPlane_0 = -mp_halfFOV[ 2 ];
  FLTNB const xPlane_N = mp_halfFOV[ 0 ];
  FLTNB const yPlane_N = mp_halfFOV[ 1 ];
  FLTNB const zPlane_N = mp_halfFOV[ 2 ];

  // Parameters for the index
  INTNB index_x( 0 );
  INTNB index_y( 0 );
  INTNB index_z( 0 );

  // Parameter for the intersection point in each axis
  // X
  FLTNB intersection_x( event1[ 0 ] < xPlane_0 ?
    xPlane_0 + szImg_X / ((FLTNB)2.0) : xPlane_N - szImg_X / ((FLTNB)2.0) );

  // Y
  FLTNB intersection_y( event1[ 1 ] < yPlane_0 ?
    yPlane_0 + szImg_Y / ((FLTNB)2.0) : yPlane_N - szImg_Y / ((FLTNB)2.0) );

  // Z
  FLTNB intersection_z( event1[ 2 ] < zPlane_0 ?
    zPlane_0 + szImg_Z / ((FLTNB)2.0) : zPlane_N - szImg_Z / ((FLTNB)2.0) );

  // Min. Max. coordinate for histo.
  FLTNB const xPlane_0_min = xPlane_0 + szImg_X / ((FLTNB)2.0);
  FLTNB const yPlane_0_min = yPlane_0 + szImg_Y / ((FLTNB)2.0);
  FLTNB const zPlane_0_min = zPlane_0 + szImg_Z / ((FLTNB)2.0);

  // Min. Max security. Do not take the last slice
  FLTNB const xPlane_0_min_safety = xPlane_0 + ( szImg_X / ((FLTNB)2.0) ) + szImg_X;
  FLTNB const yPlane_0_min_safety = yPlane_0 + ( szImg_Y / ((FLTNB)2.0) ) + szImg_Y;
  FLTNB const zPlane_0_min_safety = zPlane_0 + ( szImg_Z / ((FLTNB)2.0) ) + szImg_Z;
  FLTNB const xPlane_N_max_safety = xPlane_N - ( szImg_X / ((FLTNB)2.0) ) - szImg_X;
  FLTNB const yPlane_N_max_safety = yPlane_N - ( szImg_Y / ((FLTNB)2.0) ) - szImg_Y;
  FLTNB const zPlane_N_max_safety = zPlane_N - ( szImg_Z / ((FLTNB)2.0) ) - szImg_Z;

  // Incrementation of the plane
  FLTNB const increment[] = {
    event1[ 0 ] < xPlane_0 ? szImg_X : -szImg_X,
    event1[ 1 ] < yPlane_0 ? szImg_Y : -szImg_Y,
    event1[ 2 ] < zPlane_0 ? szImg_Z : -szImg_Z
  };

  INTNB const increment_index[] = {
    event1[ 0 ] < xPlane_0 ? 1 : -1,
    event1[ 1 ] < yPlane_0 ? 1 : -1,
    event1[ 2 ] < zPlane_0 ? 1 : -1
  };

  INTNB const first_index[] = {
    event1[ 0 ] < xPlane_0 ? 0 : nPlane_X - 1,
    event1[ 1 ] < yPlane_0 ? 0 : nPlane_Y - 1,
    event1[ 2 ] < zPlane_0 ? 0 : nPlane_Z - 1
  };

  // Values of interpolation
  FLTNB x1( 0.0 ); FLTNB x2( 0.0 );
  FLTNB y1( 0.0 ); FLTNB y2( 0.0 );
  FLTNB z1( 0.0 ); FLTNB z2( 0.0 );

  // Find the principal direction
  // X direction
  if( ::fabs( dX21 ) >= ::fabs( dY21 ) && ::fabs( dX21 ) >= ::fabs( dZ21 ) )
  {
    // Compute the weight to normalize the line
    FLTNB const factor( ::fabs( dX21 ) / ::sqrt( dX21_2 + dY21_2 + dZ21_2 ) );
    FLTNB const weight( szImg_X / factor );

    // Take the increment X
    FLTNB const incrementX = increment[ 0 ];

    // Compute first intersection in Y and Z
    FLTNB const p1X_IntersectionX1 = ( intersection_x - event1[ 0 ] );
    intersection_y = ( p1X_IntersectionX1 * dY21 / dX21 ) + event1[ 1 ];
    intersection_z = ( p1X_IntersectionX1 * dZ21 / dX21 ) + event1[ 2 ];

    // Compute the increment (next intersection minus current intersection)
    FLTNB const incrementY =
      ( ( intersection_x + incrementX - event1[ 0 ] ) * dY21 / dX21 )
      + event1[ 1 ] - intersection_y;
    FLTNB const incrementZ =
      ( ( intersection_x + incrementX - event1[ 0 ] ) * dZ21 / dX21 )
      + event1[ 2 ] - intersection_z;

    // Take first index for x
    index_x = first_index[ 0 ];

    // Loop over the Y-Z planes
    for( INTNB i = 0; i < nPlane_X; ++i )
    {
      // Check if intersection outside of the image space
      if( intersection_y < yPlane_0_min_safety
        || intersection_y >= yPlane_N_max_safety
        || intersection_z < zPlane_0_min_safety
        || intersection_z >= zPlane_N_max_safety )
      {
        index_x += increment_index[ 0 ];
        intersection_y += incrementY;
        intersection_z += incrementZ;
        continue;
      }

      // Find the index in the image
      index_y = (INTNB)( ( intersection_y - yPlane_0_min ) / szImg_Y );
      index_z = (INTNB)( ( intersection_z - zPlane_0_min ) / szImg_Z );

      // Compute distance between intersection point and Y and Z boundaries
      y2 = ::fabs( mp_boundariesY[ index_y ] + szImg_Y * ((FLTNB)0.5) - intersection_y )
        / szImg_Y;
      y1 = ((FLTNB)1.0) - y2;
      z2 = ::fabs( mp_boundariesZ[ index_z ] + szImg_Z * ((FLTNB)0.5) - intersection_z )
        / szImg_Z;
      z1 = ((FLTNB)1.0) - z2;

      INTNB index_final = index_x + index_y * nPlane_X + index_z * nPlane_XY;
      ap_ProjectionLine->AddVoxel(a_direction, index_final, ( y1 * z1 ) * weight);
      ap_ProjectionLine->AddVoxel(a_direction, index_final + nPlane_X, ( y2 * z1 ) * weight);
      ap_ProjectionLine->AddVoxel(a_direction, index_final + nPlane_XY, ( y1 * z2 ) * weight);
      ap_ProjectionLine->AddVoxel(a_direction, index_final + nPlane_X + nPlane_XY, ( y2 * z2 ) * weight);

      // Increment
      index_x += increment_index[ 0 ];
      intersection_y += incrementY;
      intersection_z += incrementZ;
    }
  } // Y direction
  else if( ::fabs( dY21 ) > ::fabs( dX21 ) && ::fabs( dY21 ) >= ::fabs( dZ21 ) )
  {
    // Compute the weight to normalize the line
    FLTNB const factor( ::fabs( dY21 ) / ::sqrt( dX21_2 + dY21_2 + dZ21_2 ) );
    FLTNB const weight( szImg_Y / factor );

    // Take the increment Y
    FLTNB const incrementY = increment[ 1 ];

    // Compute first intersection in X and Z
    FLTNB const p1Y_IntersectionY1 = ( intersection_y - event1[ 1 ] );
    intersection_x = ( p1Y_IntersectionY1 * dX21 / dY21 ) + event1[ 0 ];
    intersection_z = ( p1Y_IntersectionY1 * dZ21 / dY21 ) + event1[ 2 ];

    // Compute the increment (next intersection minus current intersection)
    FLTNB const incrementX =
      ( ( intersection_y + incrementY - event1[ 1 ] ) * dX21 / dY21 )
      + event1[ 0 ] - intersection_x;
    FLTNB const incrementZ =
      ( ( intersection_y + incrementY - event1[ 1 ] ) * dZ21 / dY21 )
      + event1[ 2 ] - intersection_z;

    // Take first index for y
    index_y = first_index[ 1 ];

    // Loop over the X-Z planes
    for( INTNB i = 0; i < nPlane_Y; ++i )
    {
      // Check if intersection outside of the image space
      if( intersection_x < xPlane_0_min_safety
        || intersection_x >= xPlane_N_max_safety
        || intersection_z < zPlane_0_min_safety
        || intersection_z >= zPlane_N_max_safety )
      {
        index_y += increment_index[ 1 ];
        intersection_x += incrementX;
        intersection_z += incrementZ;
        continue;
      }

      // Find the index in the image
      index_x = (INTNB)( ( intersection_x - xPlane_0_min ) / szImg_X );
      index_z = (INTNB)( ( intersection_z - zPlane_0_min ) / szImg_Z );

      // Compute distance between intersection point and Y and Z boundaries
      x2 = ::fabs( mp_boundariesX[ index_x ] + szImg_X * ((FLTNB)0.5) - intersection_x )
        / szImg_X;
      x1 = ((FLTNB)1.0) - x2;
      z2 = ::fabs( mp_boundariesZ[ index_z ] + szImg_Z * ((FLTNB)0.5) - intersection_z )
        / szImg_Z;
      z1 = ((FLTNB)1.0) - z2;

      uint32_t index_final = index_x + index_y * nPlane_X + index_z * nPlane_XY;
      ap_ProjectionLine->AddVoxel(a_direction, index_final, ( x1 * z1 ) * weight);
      ap_ProjectionLine->AddVoxel(a_direction, index_final + 1, ( x2 * z1 ) * weight);
      ap_ProjectionLine->AddVoxel(a_direction, index_final + nPlane_XY, ( x1 * z2 ) * weight);
      ap_ProjectionLine->AddVoxel(a_direction, index_final + 1 +  nPlane_XY, ( x2 * z2 ) * weight);

      // Increment
      index_y += increment_index[ 1 ];
      intersection_x += incrementX;
      intersection_z += incrementZ;
    }
  }
  else // Z direction
  {
    // Compute the weight to normalize the line
    FLTNB const factor( ::fabs( dZ21 ) / ::sqrt( dX21_2 + dY21_2 + dZ21_2 ) );
    FLTNB const weight( szImg_Z / factor );

    // Take the increment Z
    FLTNB const incrementZ = increment[ 2 ];

    // Compute first intersection in X and Y
    FLTNB const p1Z_IntersectionZ1 = ( intersection_z - event1[ 2 ] );
    intersection_y = ( p1Z_IntersectionZ1 * dY21 / dZ21 ) + event1[ 1 ];
    intersection_x = ( p1Z_IntersectionZ1 * dX21 / dZ21 ) + event1[ 0 ];

    // Compute the increment (next intersection minus current intersection)
    FLTNB const incrementY =
      ( ( intersection_z + incrementZ - event1[ 2 ] ) * dY21 / dZ21 )
      + event1[ 1 ] - intersection_y;
    FLTNB const incrementX =
      ( ( intersection_z + incrementZ - event1[ 2 ] ) * dX21 / dZ21 )
      + event1[ 0 ] - intersection_x;

    // Take first index for z
    index_z = first_index[ 2 ];

    // Loop over the Y-Z planes
    for( INTNB i = 0; i < nPlane_Z; ++i )
    {
      // Check if intersection outside of the image space
      if( intersection_y < yPlane_0_min_safety
        || intersection_y >= yPlane_N_max_safety
        || intersection_x < xPlane_0_min_safety
        || intersection_x >= xPlane_N_max_safety )
      {
        index_z += increment_index[ 2 ];
        intersection_y += incrementY;
        intersection_x += incrementX;
        continue;
      }

      // Find the index in the image
      index_y = (INTNB)( ( intersection_y - yPlane_0_min ) / szImg_Y );
      index_x = (INTNB)( ( intersection_x - xPlane_0_min ) / szImg_X );

      // Compute distance between intersection point and Y and Z boundaries
      y2 = ::fabs( mp_boundariesY[ index_y ] + szImg_Y * ((FLTNB)0.5) - intersection_y )
        / szImg_Y;
      y1 = ((FLTNB)1.0) - y2;
      x2 = ::fabs( mp_boundariesX[ index_x ] + szImg_X * ((FLTNB)0.5) - intersection_x )
        / szImg_X;
      x1 = ((FLTNB)1.0) - x2;

      uint32_t index_final = index_x + index_y * nPlane_X + index_z * nPlane_XY;
      ap_ProjectionLine->AddVoxel(a_direction, index_final, ( x1 * y1 ) * weight);
      ap_ProjectionLine->AddVoxel(a_direction, index_final + 1, ( x2 * y1 ) * weight);
      ap_ProjectionLine->AddVoxel(a_direction, index_final + nPlane_X, ( x1 * y2 ) * weight);
      ap_ProjectionLine->AddVoxel(a_direction, index_final + nPlane_X + 1, ( x2 * y2 ) * weight);

      // Increment
      index_z += increment_index[ 2 ];
      intersection_y += incrementY;
      intersection_x += incrementX;
    }
  }
  return 0;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

int iProjectorJoseph::ProjectWithTOFPos(int a_Projector, oProjectionLine* ap_ProjectionLine)
{
  Cerr("***** iProjectorJoseph::ProjectWithTOFPos() -> Not yet implemented !" << endl);
  return 1;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

int iProjectorJoseph::ProjectWithTOFBin(int a_Projector, oProjectionLine* ap_ProjectionLine)
{
  Cerr("***** iProjectorJoseph::ProjectWithTOFBin() -> Not yet implemented !" << endl);
  return 1;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================