itkfemimagemetricload.txx

DTMK软件开发包,此为开源软件,是一款很好的医学图像开发资源.

  for( unsigned int k = 0; k < ImageDimension; k++ )
    { 
    //Set the size of the image region
    parameters[k]= InVec[k+ImageDimension]; // this gives the translation by the vector field 
    rindex[k] =(long)(InVec[k]+InVec[k+ImageDimension]+0.5);  // where the piece of reference image currently lines up under the above translation
    tindex[k]= (long)(InVec[k]+0.5)-(long)m_MetricRadius[k]/2;  // position in reference image
    int hibordercheck=(int)tindex[k]+(int)m_MetricRadius[k]-(int)m_TarSize[k];
    int lobordercheck=(int)tindex[k]-(int)m_MetricRadius[k];
    if (hibordercheck > 0) regionRadius[k]=m_MetricRadius[k]-(long)hibordercheck-1;
    else if (lobordercheck < 0) regionRadius[k]=m_MetricRadius[k]+(long)lobordercheck;
    else regionRadius[k]=m_MetricRadius[k];  
    tindex[k]= (long)(InVec[k]+0.5)-(long)regionRadius[k]/2;  // position in reference image
    }

// Set the associated region

  requestedRegion.SetSize(regionRadius);
  requestedRegion.SetIndex(tindex);

  m_TarImage->SetRequestedRegion(requestedRegion);  
  m_Metric->SetFixedImageRegion( m_TarImage->GetRequestedRegion() );

//--------------------------------------------------------
// Get metric values

  typename MetricBaseType::MeasureType     measure=0.0;
  try
    { 
    measure=m_Metric->GetValue( parameters);
    }
  catch( ... )
    {
    // do nothing we dont care if the metric lies outside the image sometimes
    //std::cerr << e << std::endl;
    }
      
 
  return (Float) measure;
}


template<class TMoving,class TFixed>
typename ImageMetricLoad<TMoving , TFixed>::VectorType 
ImageMetricLoad<TMoving , TFixed>::MetricFiniteDiff
( VectorType  Gpos,
  VectorType  Gsol ) 
{

  typename MetricBaseType::MeasureType     measure;
  ParametersType parameters( ImageDimension );
  typename FixedType::RegionType requestedRegion;
  typename FixedType::IndexType tindex;
  FixedRadiusType regionRadius;

  VectorType OutVec;
  OutVec.set_size(ImageDimension);

  for( unsigned int k = 0; k < ImageDimension; k++ )
    { 
    parameters[k]= Gsol[k]; // this gives the translation by the vector field 
    tindex[k]= (long)(Gpos[k]+0.5)-(long)m_MetricRadius[k]/2;  // position in reference image
    if (tindex[k] > m_TarSize[k]-1 || tindex[k] < 0) tindex[k]=(long)(Gpos[k]+0.5);
    int hibordercheck=(int)tindex[k]+(int)m_MetricRadius[k]-(int)m_TarSize[k];
    int lobordercheck=(int)tindex[k]-(int)m_MetricRadius[k];
    if (hibordercheck >= 0) regionRadius[k]=m_MetricRadius[k]-(long)hibordercheck-1;
    else if (lobordercheck < 0) regionRadius[k]=m_MetricRadius[k]+(long)lobordercheck;
    else regionRadius[k]=m_MetricRadius[k];  
    tindex[k]= (long)(Gpos[k]+0.5)-(long)regionRadius[k]/2;  // position in reference image
    }
  
  unsigned int row;
  typename ImageType::IndexType difIndex[ImageDimension][2];
  
  typename MetricBaseType::MeasureType   dPixL,dPixR;
  for(row=0; row< ImageDimension;row++)
    {
    difIndex[row][0]=tindex;
    difIndex[row][1]=tindex;
    if (tindex[row] < m_TarSize[row]-1)
      {
      difIndex[row][0][row]=tindex[row]+1;
      }
    if (tindex[row] > 0 )
      {
      difIndex[row][1][row]=tindex[row]-1;
      }
    try
      { 
      requestedRegion.SetIndex(difIndex[row][1]);
      requestedRegion.SetSize(regionRadius);
      m_TarImage->SetRequestedRegion(requestedRegion);  
      m_Metric->SetFixedImageRegion( m_TarImage->GetRequestedRegion() );
      dPixL=m_Metric->GetValue( parameters);
      }
    catch( ... )
      {
      dPixL=0.0;
      } 
    try
      { 
      requestedRegion.SetIndex(difIndex[row][0]);
      requestedRegion.SetSize(regionRadius);
      m_TarImage->SetRequestedRegion(requestedRegion);  
      m_Metric->SetFixedImageRegion( m_TarImage->GetRequestedRegion() );
      dPixR=m_Metric->GetValue( parameters);
      }
    catch( ... )
      {
      dPixR=0.0;
      }
    
    OutVec[row]=dPixL-dPixR;
    }
  return OutVec;
}


template<class TMoving,class TFixed>
typename ImageMetricLoad<TMoving , TFixed>::VectorType 
ImageMetricLoad<TMoving , TFixed>::GetPolynomialFitToMetric
( VectorType  Gpos,
  VectorType  Gsol ) 
{

//discrete orthogonal polynomial fitting
//see p.394-403 haralick computer and robot vision
//
//here, use chebyshev polynomials for fitting a plane to the data
//
//f(x,y,z) = a0 + a1*x + a2*y + a3*z
//
  ParametersType parameters( ImageDimension );
  typename FixedType::RegionType requestedRegion;
  typename FixedType::IndexType tindex;
  FixedRadiusType regionRadius;

  typename ImageType::IndexType temp;

  VectorType chebycoefs; // gradient direction
  chebycoefs.set_size(ImageDimension);
  double chebycoefs0=0.0;  // the constant term
  double datatotal=0.0;
  double a0norm=1.0;
  double a1norm=1.0/2.0;
  
  double met, ind1,ind2;
  double inds[3]; inds[0]=-1.0;  inds[1]=0.0;  inds[2]=1.0;

  for( unsigned int k = 0; k < ImageDimension; k++ )
    { 
    a0norm/=3.0;
    if (k < ImageDimension-1) a1norm/=3.0;
    chebycoefs[k]=0.0;
    parameters[k]= Gsol[k]; // this gives the translation by the vector field 
    tindex[k]= (long)(Gpos[k]+0.5)-(long)m_MetricRadius[k]/2;  // position in reference image
    if (tindex[k] > m_TarSize[k]-1 || tindex[k] < 0) tindex[k]=(long)(Gpos[k]+0.5);
    int hibordercheck=(int)tindex[k]+(int)m_MetricRadius[k]-(int)m_TarSize[k];
    int lobordercheck=(int)tindex[k]-(int)m_MetricRadius[k];
    if (hibordercheck >= 0) regionRadius[k]=m_MetricRadius[k]-(long)hibordercheck-1;
    else if (lobordercheck < 0) regionRadius[k]=m_MetricRadius[k]+(long)lobordercheck;
    else regionRadius[k]=m_MetricRadius[k];
    tindex[k]= (long)(Gpos[k]+0.5)-(long)regionRadius[k]/2;  // position in reference image
    }
  

  if (ImageDimension==2)
    {
    
    double measure[3][3];
    for(int row=-1; row< 2; row++)
      {
      for(int col=-1; col< 2; col++)
        {
      
        temp[0]=tindex[0]+(long)row;
        temp[1]=tindex[1]+(long)col;
        
        for (unsigned int i=0; i<ImageDimension; i++)
          {
          if (temp[i] > m_TarSize[i]-1)
            {
            temp[i]=m_TarSize[i]-1;
            }
          else if (temp[i] < 0 )
            {
            temp[i]=0;
            }
          }

        requestedRegion.SetIndex(temp);
        requestedRegion.SetSize(regionRadius);
        m_TarImage->SetRequestedRegion(requestedRegion);  
        m_Metric->SetFixedImageRegion( m_TarImage->GetRequestedRegion() );
        measure[row+1][col+1]=0.0;
        
        try
          { 
          measure[row+1][col+1]=m_Metric->GetValue( parameters);
          }
        catch( ... )
          {
          }
 
    
        datatotal+=measure[row+1][col+1];
        }
      }
    for( unsigned int cb1 = 0; cb1 < 3; cb1++ )
      { 
      for( unsigned int cb2 = 0; cb2 < 3; cb2++ ) 
        {
        met=measure[cb1][cb2];
        ind1=inds[cb1]*a1norm;
        ind2=inds[cb2]*a1norm;
        chebycoefs[0]+=met*ind1;
        chebycoefs[1]+=met*ind2;
        }
      }
    }
  else if (ImageDimension == 3)
    {

    double measure3D[3][3][3];
    for(int row=-1; row< 2; row++)
      {
      for(int col=-1; col< 2; col++)
        {
        for(int z=-1; z< 2; z++)
          {

          temp[0]=tindex[0]+(long)row;
          temp[1]=tindex[1]+(long)col;
          temp[2]=tindex[2]+(long)z;

          for (unsigned int i=0; i<ImageDimension; i++)
            {
            if (temp[i] > m_TarSize[i]-1)
              {
              temp[i]=m_TarSize[i]-1;
              }
            else if (temp[i] < 0 )
              {
              temp[i]=0;
              }
            }

          requestedRegion.SetIndex(temp);
          requestedRegion.SetSize(regionRadius);
          m_TarImage->SetRequestedRegion(requestedRegion);  
          m_Metric->SetFixedImageRegion( m_TarImage->GetRequestedRegion() );
          measure3D[row+1][col+1][z+1]=0.0;
   
          try
            { 
            measure3D[row+1][col+1][z+1]=m_Metric->GetValue( parameters);
            }
          catch( ... )
            {
            }
 
    
          datatotal+=measure3D[row+1][col+1][z+1];
          }
        }
      }
    for( unsigned int cb1 = 0; cb1 < 2; cb1++ )
      {
      for( unsigned int cb2 = 0; cb2 < 2; cb2++ )
        { 
        for( unsigned int cb3 = 0; cb3 < 2; cb3++ ) 
          {
          chebycoefs[0]+=measure3D[cb1][cb2][cb3]*inds[cb1]*a1norm;
          chebycoefs[1]+=measure3D[cb1][cb2][cb3]*inds[cb2]*a1norm;
          chebycoefs[2]+=measure3D[cb1][cb2][cb3]*inds[cb3]*a1norm;
          }
        }
      }
    }
  
  chebycoefs0=a0norm*datatotal;
//  std::cout << " cb " << chebycoefs << std::endl;
  return chebycoefs;

}


template<class TMoving,class TFixed> 
int ImageMetricLoad<TMoving,TFixed>::CLID()
{
  static const int CLID_ = FEMOF::Register( ImageMetricLoad::NewB,(std::string("ImageMetricLoad(")
                                                                   +typeid(TMoving).name()+","+typeid(TFixed).name()+")").c_str());
  return CLID_;
}


template<class TMoving,class TFixed> 
const int ImageMetricLoad<TMoving,TFixed>::DummyCLID=ImageMetricLoad<TMoving,TFixed>::CLID();


} // end namespace fem
} // end namespace itk

#endif