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

[Lorensen, William E (Research)]

Which compiler are you using? What optimization level? ITK is very sensitive to optimization. On Windows you should build RelWithDebInfo or Release. Debug builds are very slow.
 
Bill

-----Original Message-----
From: insight-users-bounces+lorensen=crd.ge.com at itk.org [mailto:insight-users-bounces+lorensen=crd.ge.com at itk.org]On Behalf Of Renaud Isabelle
Sent: Monday, November 28, 2005 5:04 PM
To: insight-users at itk.org
Subject: [Insight-users] my ITK code is slow: what is wrong? 


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 croppedSizee;
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  <http://us.rd.yahoo.com/messenger/mail_taglines/default/*http://fr.messenger.yahoo.com> le ici ! 

-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://public.kitware.com/pipermail/insight-users/attachments/20051129/e2694b67/attachment.htm