[Insight-users] my ITK code is slow: what is wrong?

[Renaud Isabelle]

Hi,
   
  I translate a matlab procedure into C++ with ITK to make things quicker. Now, comparing both of them, I found that results are the same but my C++ version with ITK is twice slower than matlab, whereas the contrary was expected. Could someone tell me what is wrong with my code: is there things that I can improve and time that I can save?
   
  Here is what I did:
     
  constructeur()
  {
    this->input = NULL;
    this->image1 = NULL;
    this->image2 = NULL;
    this->coeff = NULL;
    this->shiftx = NULL;
    this->shifty = NULL;
   
    //filters
    extractSlice = ExtractSliceFilter::New();
    extractRegion = ExtractRegionFilter::New();
    fft1 = FFTFilterType::New();
    fft2 = FFTFilterType::New(); 
    mult = MultiplyFilter::New(); 
    ifft = IFFTFilterType::New();
   
    //release extra data
    extractSlice->GetOutput()->ReleaseDataFlagOn();
    extractRegion->GetOutput()->ReleaseDataFlagOn();
    fft1->GetOutput()->ReleaseDataFlagOn();
    fft2->GetOutput()->ReleaseDataFlagOn();
    mult->GetOutput()->ReleaseDataFlagOn();
    ifft->GetOutput()->ReleaseDataFlagOn();
    
    mult->SetInput1(fft1->GetOutput());
    mult->SetInput2(fft2->GetOutput());
    ifft->SetInput(mult->GetOutput()); 
  }
  void ComputeElastogram::SetInput( ImageType3D::ConstPointer image)
  {
  this->input = image;
  extractSlice->SetInput(this->input);
  this->calcul_deplacements();}
   
  void ComputeElastogram::calcul_deplacements() 
  {
  unsigned int k1,k2,n;
   
  //Creation des matrices de deplacement
  ImageType::IndexType start;
  start.Fill(0); 
  ImageType::SizeType sz;
  sz[0] = dim_s[0]; 
  sz[1] = dim_s[1];
  ImageType::RegionType region;
  region.SetSize( sz );
  region.SetIndex( start );
  shifty = ImageType::New();
  shifty->SetRegions( region );
  shifty->Allocate();
  shiftx = ImageType::New();
  shiftx->SetRegions( region );
  shiftx->Allocate();
   
  rowa = dim_fen[0]; cola = dim_fen[1];
  /*--------Cropped ROI-------------*/
  ImageType3D::SizeType croppedSize;
  croppedSize[0] = dim_img[0];//dim_roi[0];
  croppedSize[1] = dim_img[1];//dim_roi[1];
  croppedSize[2] = 0;
   
  
  ImageType3D::IndexType croppedStart;
  croppedStart[0] = 0;//bords_c[0];
  croppedStart[1] = 0;//bords_l[0];
  croppedStart[2] = 0; 
   
  ImageType3D::RegionType croppedRegion;
  croppedRegion.SetSize(croppedSize);
  croppedRegion.SetIndex(croppedStart);
   
  /*------------------------------------*/
  /*--------WINDOWS--------------------*/
  ImageType::SizeType extractedSize;
  extractedSize[0] =80;
  extractedSize[1] = 20;
  ImageType::IndexType extractedStart;
  extractedStart.Fill(0); 
  ImageType::RegionType desiredRegion;
  desiredRegion.SetSize(extractedSize); 
  desiredRegion.SetIndex(extractedStart); 
  /*------------------------------------*/
   
  /*--------Hanning Window--------------*/ 
  ImageType::Pointer hanning = ImageType::New();
  hanning->SetRegions( desiredRegion );
  hanning->Allocate();
  float **w = hanning_2D(dim_fen[1],dim_fen[0]);//taille de la fenetre 20lignes * 80colonnes
  ImageRegionIterator it(hanning, desiredRegion); 
  for(int i=0;i<dim_fen[1];i++)
  for(int j=0; j<dim_fen[0]; j++)
  {
    it.Set(w[i][j]);
    ++it;
  }
  free_fmatrix_2d(w);
   
  
  fhan = WindowFilter::New();
  fhan->GetOutput()->ReleaseDataFlagOn();
  fhan->SetInput1(hanning);
  fhan->SetInput2(extractRegion->GetOutput());
  /*------------------------------------*/
   
  ImageRegionIterator it1(shifty, shifty->GetLargestPossibleRegion()); 
  ImageRegionIterator it2(shiftx, shiftx->GetLargestPossibleRegion());
  
   
  croppedStart[2] = 0;
  croppedRegion.SetIndex(croppedStart); 
  extractSlice->SetExtractionRegion(croppedRegion);
  extractSlice->Update();
  ImageType2D::Pointer im0 = extractSlice->GetOutput();
  im0->DisconnectPipeline(); 
   
  croppedStart[2] = 1;
  croppedRegion.SetIndex(croppedStart); 
  extractSlice->SetExtractionRegion(croppedRegion);
  extractSlice->Modified();
  extractSlice->Update();
  ImageType2D::Pointer im1 = extractSlice->GetOutput();
  im1->DisconnectPipeline();
  
  //-------------------------------------------------------------------
   
  it1.GoToBegin();
  it2.GoToBegin(); 
   
  
  // Correlation 
  for(k1 = 0; k1 < dim_s[0]; k1++)
  {
  for(k2 = 0; k2 < dim_s[1]; k2++)
  { 
    extractedStart[0] = k1*pas_fen[0]+ bords_c[0];
    extractedStart[1] = k2*pas_fen[1]+ bords_l[0];
    desiredRegion.SetIndex(extractedStart); 
  
    extractRegion->SetRegionOfInterest(desiredRegion); 
    extractRegion->SetInput(im0); 
    fhan->Update();
    image1 = fhan->GetOutput();
    image1->DisconnectPipeline();
    extractRegion->SetInput(im1); 
    fhan->Update();
    image2 = fhan->GetOutput();
    image2->DisconnectPipeline();
  
    int *shft = shift(image1,image2);
  
    it1.Set(shft[0]);
    it2.Set(shft[1]);
    ++it1;
    ++it2;
  
  }//k1
  }//k2 
   
  timer.Stop();
  fprintf(fp, "\n Elapsed time: %f\n", timer.GetMeanTime());
  
  }
   
  int* ComputeElastogram::shift( ImageType::Pointer image1, ImageType::Pointer image2) 
  {
  int *shft = new int[2];
  this->calcul_coef(image2,image1);
   
  ImageType::IndexType index;
  index[0] = 0; index[1] = 0;
  unsigned int row, col; 
  row=col=0;
  ImageType::PixelType max = coeff->GetPixel(index); 
  
  ImageRegionIterator it(coeff, coeff->GetLargestPossibleRegion()); 
  while( !it.IsAtEnd() )
  {
    if(it.Get()>max)
    {
    max = it.Get();
    index=it.GetIndex();
    row = index[0];
    col = index[1];
  }
  ++it;
  }
  
  shft[0] = - rowa * (row > (int)(rowa/2) ) + row;
  shft[1] = - cola * (col > (int)(cola/2)) + col;}
  return shft;
  }
  void ComputeElastogram::calcul_coef( ImageType::Pointer image1, ImageType::Pointer image2) 
  { 
  // Forward FFT filter
  fft1->SetInput( image1 );
  fft2->SetInput( image2 );
  
  try{ ifft->Update(); }
  catch(itk::ExceptionObject &e){
  CString msg;
  msg.Format("Erreur %s in %s\n", e.GetDescription(), e.GetLocation());
  return;
  }
  // Iterator which traverse the image
  coeff = ifft->GetOutput();
  }
   
  Actually, I don't see what is wrong. I tried to release data and so on, but it seems to make any difference. 
   
  Please help,
   
  Isabelle


		
---------------------------------
 Appel audio GRATUIT partout dans le monde avec le nouveau Yahoo! Messenger
 Téléchargez le ici !  
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://public.kitware.com/pipermail/insight-users/attachments/20051128/36964b2b/attachment-0001.htm