← Back to branch summary

~ubuntu-branches/ubuntu/precise/insighttoolkit/precise

« back to all changes in this revision

Viewing changes to Code/Numerics/FEM/itkFEMImageMetricLoad.txx

  • Committer: Bazaar Package Importer
  • Author(s): Steve M. Robbins
  • Date: 2008-05-31 12:07:29 UTC
  • mfrom: (3.1.3 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080531120729-7g7layn480le43ko
Tags: 3.6.0-3
http://bugs.debian.org/478500
http://bugs.debian.org/478940
http://bugs.debian.org/483312
* debian/patches/gccxml-workaround.patch: New.  Work around gccxml issue
  with #include_next; c.f. http://www.gccxml.org/Bug/view.php?id=7134.  
* debian/patches/gcc43.patch: include <cstring> in itkNeighbourhood.h.
  This only showed up in the tcl wrapping step.

* Above two entries fix FTBFS for GCC 4.3-based systems.
  Closes: #478500.

* debian/patches/sharedforward.patch: New.  Ensure that linux/sparc
  systems are not also configured as a SUN sparc system, which requires
  SUN header sys/isa_defs.h.  Closes: #478940, #483312.

Show diffs side-by-side

added added

removed removed

Lines of Context:
Code/Numerics/FEM/itkFEMImageMetricLoad.txx
3
3
  Program:   Insight Segmentation & Registration Toolkit (ITK)
4
4
  Module:    $RCSfile: itkFEMImageMetricLoad.txx,v $
5
5
  Language:  C++
6
 
  Date:      $Date: 2006/03/19 04:37:20 $
7
 
  Version:   $Revision: 1.31 $
 
6
  Date:      $Date: 2008-01-07 21:28:29 $
 
7
  Version:   $Revision: 1.32 $
8
8
 
9
9
=========================================================================*/
10
10
#ifndef _itkFEMImageMetricLoad_txx_
42
42
//  typename MovingType::IndexType rindex;
43
43
  // initialize the offset/vector part
44
44
  for( unsigned int k = 0; k < ImageDimension; k++ )
45
 
  { 
46
 
  //Set the size of the image region
 
45
    { 
 
46
    //Set the size of the image region
47
47
    size[k] = 1;
48
48
    tindex[k]=0;
49
 
  }
 
49
    }
50
50
 
51
51
// Set the number of integration points to zero (default for an element)
52
52
 
70
70
// internal components: Interpolator, Transform and Images
71
71
//------------------------------------------------------------
72
72
  try 
73
 
  {
 
73
    {
74
74
    m_Metric->Initialize();
75
 
  } catch( ExceptionObject & e )
76
 
  {
 
75
    }
 
76
  catch( ExceptionObject & e )
 
77
    {
77
78
    std::cout << "Metric initialization failed" << std::endl;
78
79
    std::cout << "Reason " << e.GetDescription() << std::endl;
79
 
  }
 
80
    }
80
81
 
81
82
}
82
83
 
91
92
  m_Sign=1.0;
92
93
 
93
94
  for (unsigned int i=0; i<ImageDimension; i++)
94
 
  {
 
95
    {
95
96
    m_MetricRadius[i] = 1;
96
 
  }
 
97
    }
97
98
  m_MetricGradientImage=NULL;
98
99
 
99
100
}
101
102
 
102
103
template<class TMoving,class TFixed>
103
104
typename ImageMetricLoad<TMoving , TFixed>::Float 
104
 
ImageMetricLoad<TMoving , TFixed>::EvaluateMetricGivenSolution( Element::ArrayType* el,Float step)
 
105
ImageMetricLoad<TMoving , TFixed>::EvaluateMetricGivenSolution( Element::ArrayType* element,Float step)
105
106
{
106
107
  Float energy=0.0,defe=0.0; 
107
108
 
111
112
  Element::MatrixType solmat;
112
113
  Element::Float w;
113
114
 
114
 
  Element::ArrayType::iterator elt=el->begin();
 
115
  Element::ArrayType::iterator elt=element->begin();
115
116
  const unsigned int Nnodes=(*elt)->GetNumberOfNodes();
116
117
 
117
118
  solmat.set_size(Nnodes*ImageDimension,1);
118
119
 
119
 
  for(  ; elt!=el->end(); elt++) 
120
 
  {
 
120
  for(  ; elt!=element->end(); elt++) 
 
121
    {
121
122
    for(unsigned int i=0; i<m_NumberOfIntegrationPoints; i++)
122
 
    {
 
123
      {
123
124
      dynamic_cast<Element*>(&*(*elt))->GetIntegrationPointAndWeight(i,ip,w,m_NumberOfIntegrationPoints); // FIXME REMOVE WHEN ELEMENT NEW IS BASE CLASS
124
125
      shapef = (*elt)->ShapeFunctions(ip);
125
126
 
127
128
      Float detJ=(*elt)->JacobianDeterminant(ip);
128
129
        
129
130
      for(unsigned int f=0; f<ImageDimension; f++)
130
 
      {
 
131
        {
131
132
        solval=0.0;
132
133
        posval=0.0;
133
134
        for(unsigned int n=0; n<Nnodes; n++)
134
 
        {
 
135
          {
135
136
          posval+=shapef[n]*(((*elt)->GetNodeCoordinates(n))[f]);
136
137
          float nodeval=( (m_Solution)->GetSolutionValue( (*elt)->GetNode(n)->GetDegreeOfFreedom(f) , m_SolutionIndex)
137
 
            +(m_Solution)->GetSolutionValue( (*elt)->GetNode(n)->GetDegreeOfFreedom(f) , m_SolutionIndex2)*step);
 
138
                          +(m_Solution)->GetSolutionValue( (*elt)->GetNode(n)->GetDegreeOfFreedom(f) , m_SolutionIndex2)*step);
138
139
      
139
140
          solval+=shapef[n] * nodeval;   
140
141
          solmat[(n*ImageDimension)+f][0]=nodeval;
141
 
        }
 
142
          }
142
143
        InVec[f]=posval;
143
144
        InVec[f+ImageDimension]=solval;
144
 
      }
 
145
        }
145
146
 
146
147
      float tempe=0.0;
147
148
      try
148
 
      {
149
 
      tempe=vcl_fabs(GetMetric(InVec));
150
 
      }
 
149
        {
 
150
        tempe=vcl_fabs(GetMetric(InVec));
 
151
        }
151
152
      catch( itk::ExceptionObject & )
152
 
      { 
153
 
      // do nothing we dont care if the metric region is outside the image
154
 
      //std::cerr << e << std::endl;
155
 
      }
 
153
        { 
 
154
        // do nothing we dont care if the metric region is outside the image
 
155
        //std::cerr << e << std::endl;
 
156
        }
156
157
      for(unsigned int n=0; n<Nnodes; n++)
157
 
      {
 
158
        {
158
159
        itk::fem::Element::Float temp=shapef[n]*tempe*w*detJ;
159
160
        energy+=temp;
160
 
      }
161
 
    }  
 
161
        }
 
162
      }  
162
163
    
163
164
    defe+=0.0;//(double)(*elt)->GetElementDeformationEnergy( solmat );
164
 
  }
 
165
    }
165
166
   
166
167
  //std::cout << " def e " << defe << " sim e " << energy*m_Gamma << std::endl;
167
168
  return vcl_fabs((double)energy*(double)m_Gamma-(double)defe);
189
190
 
190
191
  VectorType OutVec;
191
192
  
192
 
  for( unsigned int k = 0; k < ImageDimension; k++ ) {
 
193
  for( unsigned int k = 0; k < ImageDimension; k++ )
 
194
    {
193
195
    if ( vnl_math_isnan(Gpos[k])  || vnl_math_isinf(Gpos[k]) ||
194
 
        vnl_math_isnan(Gsol[k])  || vnl_math_isinf(Gsol[k]) ||
 
196
         vnl_math_isnan(Gsol[k])  || vnl_math_isinf(Gsol[k]) ||
195
197
         vcl_fabs(Gpos[k]) > 1.e33  || vcl_fabs(Gsol[k]) > 1.e33  ) 
196
 
    {
 
198
      {
197
199
      OutVec.set_size(ImageDimension);  OutVec.fill(0.0);  return OutVec;
 
200
      }
198
201
    }
199
 
  }
200
202
//  OutVec=this->MetricFiniteDiff(Gpos,Gsol); // gradient direction
201
203
//  OutVec=this->GetPolynomialFitToMetric(Gpos,Gsol); // gradient direction
202
204
//  for( unsigned int k = 0; k < ImageDimension; k++ ) {
215
217
 
216
218
  int lobordercheck=0,hibordercheck=0;
217
219
  for( unsigned int k = 0; k < ImageDimension; k++ )
218
 
  { 
219
 
  //Set the size of the image region
 
220
    { 
 
221
    //Set the size of the image region
220
222
    parameters[k]= Gsol[k]; // this gives the translation by the vector field 
221
223
    rindex[k] =(long)(Gpos[k]+Gsol[k]+0.5);  // where the piece of reference image currently lines up under the above translation
222
224
    tindex[k]= (long)(Gpos[k]+0.5)-(long)m_MetricRadius[k]/2;  // position in reference image
226
228
    else if (lobordercheck < 0) regionRadius[k]=m_MetricRadius[k]+(long)lobordercheck;
227
229
    else regionRadius[k]=m_MetricRadius[k];
228
230
    tindex[k]= (long)(Gpos[k]+0.5)-(long)regionRadius[k]/2;  // position in reference image
229
 
  }
 
231
    }
230
232
 
231
233
// Set the associated region
232
234
 
243
245
  typename MetricBaseType::DerivativeType  derivative;
244
246
 
245
247
  try
246
 
  { 
 
248
    { 
247
249
    m_Metric->GetValueAndDerivative( parameters, measure, derivative );
248
 
  //  m_Metric->GetDerivative( parameters, derivative );
249
 
  }
 
250
    //  m_Metric->GetDerivative( parameters, derivative );
 
251
    }
250
252
  catch( ... )
251
 
  {
252
 
  // do nothing we don't care if the metric lies outside the image sometimes
253
 
  //std::cerr << e << std::endl;
254
 
  }
 
253
    {
 
254
    // do nothing we don't care if the metric lies outside the image sometimes
 
255
    //std::cerr << e << std::endl;
 
256
    }
255
257
 
256
258
  m_Energy+=(double)measure;
257
259
  float gmag=0.0;
258
260
  for( unsigned int k = 0; k < ImageDimension; k++ )
259
 
  {
 
261
    {
260
262
    if (lobordercheck < 0 || hibordercheck >=0 ||
261
 
       vnl_math_isnan(derivative[k])  || vnl_math_isinf(derivative[k]) ) 
262
 
    {
 
263
        vnl_math_isnan(derivative[k])  || vnl_math_isinf(derivative[k]) ) 
 
264
      {
263
265
      OutVec[k]=0.0;
264
 
    } 
 
266
      } 
265
267
    else OutVec[k]= m_Sign*m_Gamma*derivative[k];
266
268
    gmag+=OutVec[k]*OutVec[k];
267
 
  }
 
269
    }
268
270
  if (gmag==0.0) gmag=1.0;
269
 
 // NOTE : POSSIBLE THAT DERIVATIVE DIRECTION POINTS UP OR DOWN HILL!
270
 
 // IN FACT, IT SEEMS MEANSQRS AND NCC POINT IN DIFFT DIRS
 
271
  // NOTE : POSSIBLE THAT DERIVATIVE DIRECTION POINTS UP OR DOWN HILL!
 
272
  // IN FACT, IT SEEMS MEANSQRS AND NCC POINT IN DIFFT DIRS
271
273
  //std::cout   << " deriv " << derivative <<  " val " << measure << endl;
272
274
  //if (m_Temp !=0.0) 
273
275
  //return OutVec * vcl_exp(-1.*OutVec.magnitude()/m_Temp);
299
301
  FixedRadiusType regionRadius;
300
302
  VectorType OutVec(ImageDimension,0.0); // gradient direction
301
303
  //std::cout << " pos   translation " << InVec  << endl;
302
 
   // initialize the offset/vector part
 
304
  // initialize the offset/vector part
303
305
  for( unsigned int k = 0; k < ImageDimension; k++ )
304
 
  { 
305
 
  //Set the size of the image region
 
306
    { 
 
307
    //Set the size of the image region
306
308
    parameters[k]= InVec[k+ImageDimension]; // this gives the translation by the vector field 
307
309
    rindex[k] =(long)(InVec[k]+InVec[k+ImageDimension]+0.5);  // where the piece of reference image currently lines up under the above translation
308
310
    tindex[k]= (long)(InVec[k]+0.5)-(long)m_MetricRadius[k]/2;  // position in reference image
312
314
    else if (lobordercheck < 0) regionRadius[k]=m_MetricRadius[k]+(long)lobordercheck;
313
315
    else regionRadius[k]=m_MetricRadius[k];  
314
316
    tindex[k]= (long)(InVec[k]+0.5)-(long)regionRadius[k]/2;  // position in reference image
315
 
  }
 
317
    }
316
318
 
317
319
// Set the associated region
318
320
 
327
329
 
328
330
  typename MetricBaseType::MeasureType     measure=0.0;
329
331
  try
330
 
  { 
331
 
  measure=m_Metric->GetValue( parameters);
332
 
  }
 
332
    { 
 
333
    measure=m_Metric->GetValue( parameters);
 
334
    }
333
335
  catch( ... )
334
 
  {
335
 
  // do nothing we dont care if the metric lies outside the image sometimes
336
 
  //std::cerr << e << std::endl;
337
 
  }
 
336
    {
 
337
    // do nothing we dont care if the metric lies outside the image sometimes
 
338
    //std::cerr << e << std::endl;
 
339
    }
338
340
      
339
341
 
340
342
  return (Float) measure;
358
360
  OutVec.set_size(ImageDimension);
359
361
 
360
362
  for( unsigned int k = 0; k < ImageDimension; k++ )
361
 
  { 
 
363
    { 
362
364
    parameters[k]= Gsol[k]; // this gives the translation by the vector field 
363
365
    tindex[k]= (long)(Gpos[k]+0.5)-(long)m_MetricRadius[k]/2;  // position in reference image
364
366
    if (tindex[k] > m_TarSize[k]-1 || tindex[k] < 0) tindex[k]=(long)(Gpos[k]+0.5);
368
370
    else if (lobordercheck < 0) regionRadius[k]=m_MetricRadius[k]+(long)lobordercheck;
369
371
    else regionRadius[k]=m_MetricRadius[k];  
370
372
    tindex[k]= (long)(Gpos[k]+0.5)-(long)regionRadius[k]/2;  // position in reference image
371
 
  }
 
373
    }
372
374
  
373
375
  unsigned int row;
374
376
  typename ImageType::IndexType difIndex[ImageDimension][2];
375
377
  
376
378
  typename MetricBaseType::MeasureType   dPixL,dPixR;
377
 
  for(row=0; row< ImageDimension;row++){
 
379
  for(row=0; row< ImageDimension;row++)
 
380
    {
378
381
    difIndex[row][0]=tindex;
379
382
    difIndex[row][1]=tindex;
380
 
    if (tindex[row] < m_TarSize[row]-1) difIndex[row][0][row]=tindex[row]+1;
381
 
    if (tindex[row] > 0 )               difIndex[row][1][row]=tindex[row]-1;
382
 
 
 
383
    if (tindex[row] < m_TarSize[row]-1)
 
384
      {
 
385
      difIndex[row][0][row]=tindex[row]+1;
 
386
      }
 
387
    if (tindex[row] > 0 )
 
388
      {
 
389
      difIndex[row][1][row]=tindex[row]-1;
 
390
      }
383
391
    try
384
 
    { 
385
 
    requestedRegion.SetIndex(difIndex[row][1]);
386
 
    requestedRegion.SetSize(regionRadius);
387
 
    m_TarImage->SetRequestedRegion(requestedRegion);  
388
 
    m_Metric->SetFixedImageRegion( m_TarImage->GetRequestedRegion() );
389
 
    dPixL=m_Metric->GetValue( parameters);
390
 
    }
 
392
      { 
 
393
      requestedRegion.SetIndex(difIndex[row][1]);
 
394
      requestedRegion.SetSize(regionRadius);
 
395
      m_TarImage->SetRequestedRegion(requestedRegion);  
 
396
      m_Metric->SetFixedImageRegion( m_TarImage->GetRequestedRegion() );
 
397
      dPixL=m_Metric->GetValue( parameters);
 
398
      }
391
399
    catch( ... )
392
 
    {
 
400
      {
393
401
      dPixL=0.0;
394
 
    } 
 
402
      } 
395
403
    try
396
 
    { 
397
 
    requestedRegion.SetIndex(difIndex[row][0]);
398
 
    requestedRegion.SetSize(regionRadius);
399
 
    m_TarImage->SetRequestedRegion(requestedRegion);  
400
 
    m_Metric->SetFixedImageRegion( m_TarImage->GetRequestedRegion() );
401
 
    dPixR=m_Metric->GetValue( parameters);
402
 
    }
 
404
      { 
 
405
      requestedRegion.SetIndex(difIndex[row][0]);
 
406
      requestedRegion.SetSize(regionRadius);
 
407
      m_TarImage->SetRequestedRegion(requestedRegion);  
 
408
      m_Metric->SetFixedImageRegion( m_TarImage->GetRequestedRegion() );
 
409
      dPixR=m_Metric->GetValue( parameters);
 
410
      }
403
411
    catch( ... )
404
 
    {
 
412
      {
405
413
      dPixR=0.0;
406
 
    }
 
414
      }
407
415
    
408
416
    OutVec[row]=dPixL-dPixR;
409
 
  }
 
417
    }
410
418
  return OutVec;
411
419
}
412
420
 
425
433
//
426
434
//f(x,y,z) = a0 + a1*x + a2*y + a3*z
427
435
//
428
 
  typename MetricBaseType::MeasureType     measure;
429
436
  ParametersType parameters( ImageDimension );
430
437
  typename FixedType::RegionType requestedRegion;
431
438
  typename FixedType::IndexType tindex;
444
451
  double inds[3]; inds[0]=-1.0;  inds[1]=0.0;  inds[2]=1.0;
445
452
 
446
453
  for( unsigned int k = 0; k < ImageDimension; k++ )
447
 
  { 
 
454
    { 
448
455
    a0norm/=3.0;
449
456
    if (k < ImageDimension-1) a1norm/=3.0;
450
457
    chebycoefs[k]=0.0;
457
464
    else if (lobordercheck < 0) regionRadius[k]=m_MetricRadius[k]+(long)lobordercheck;
458
465
    else regionRadius[k]=m_MetricRadius[k];
459
466
    tindex[k]= (long)(Gpos[k]+0.5)-(long)regionRadius[k]/2;  // position in reference image
460
 
  }
 
467
    }
461
468
  
462
469
 
463
 
  if (ImageDimension==2){
464
 
 
465
 
  double measure[3][3];
466
 
  for(int row=-1; row< 2; row++){
467
 
  for(int col=-1; col< 2; col++){
468
 
 
469
 
    temp[0]=tindex[0]+(long)row;
470
 
    temp[1]=tindex[1]+(long)col;
471
 
 
472
 
    for (unsigned int i=0; i<ImageDimension; i++){
473
 
      if (temp[i] > m_TarSize[i]-1) temp[i]=m_TarSize[i]-1;
474
 
      else if (temp[i] < 0 ) temp[i]=0;
475
 
    }
476
 
 
477
 
    requestedRegion.SetIndex(temp);
478
 
    requestedRegion.SetSize(regionRadius);
479
 
    m_TarImage->SetRequestedRegion(requestedRegion);  
480
 
    m_Metric->SetFixedImageRegion( m_TarImage->GetRequestedRegion() );
481
 
    measure[row+1][col+1]=0.0;
482
 
   
483
 
    try
484
 
    { 
485
 
      measure[row+1][col+1]=m_Metric->GetValue( parameters);
486
 
    }
487
 
    catch( ... )
488
 
    {
489
 
    }
490
 
 
491
 
    
492
 
     datatotal+=measure[row+1][col+1];
493
 
  }}
494
 
  for( unsigned int cb1 = 0; cb1 < 3; cb1++ ) 
495
 
  for( unsigned int cb2 = 0; cb2 < 3; cb2++ ) 
496
 
  {
497
 
    met=measure[cb1][cb2];
498
 
    ind1=inds[cb1]*a1norm;
499
 
    ind2=inds[cb2]*a1norm;
500
 
    chebycoefs[0]+=met*ind1;
501
 
    chebycoefs[1]+=met*ind2;
502
 
  }
503
 
  }
504
 
  else if (ImageDimension == 3) {
505
 
 
506
 
  double measure3D[3][3][3];
507
 
  for(int row=-1; row< 2; row++){
508
 
  for(int col=-1; col< 2; col++){
509
 
  for(int z=-1; z< 2; z++){
510
 
 
511
 
    temp[0]=tindex[0]+(long)row;
512
 
    temp[1]=tindex[1]+(long)col;
513
 
    temp[2]=tindex[2]+(long)z;
514
 
 
515
 
    for (unsigned int i=0; i<ImageDimension; i++){
516
 
      if (temp[i] > m_TarSize[i]-1) temp[i]=m_TarSize[i]-1;
517
 
      else if (temp[i] < 0 ) temp[i]=0;
518
 
    }
519
 
 
520
 
    requestedRegion.SetIndex(temp);
521
 
    requestedRegion.SetSize(regionRadius);
522
 
    m_TarImage->SetRequestedRegion(requestedRegion);  
523
 
    m_Metric->SetFixedImageRegion( m_TarImage->GetRequestedRegion() );
524
 
    measure3D[row+1][col+1][z+1]=0.0;
525
 
   
526
 
    try
527
 
    { 
528
 
      measure3D[row+1][col+1][z+1]=m_Metric->GetValue( parameters);
529
 
    }
530
 
    catch( ... )
531
 
    {
532
 
    }
533
 
 
534
 
    
535
 
     datatotal+=measure3D[row+1][col+1][z+1];
536
 
  }}}
537
 
  for( unsigned int cb1 = 0; cb1 < 2; cb1++ ) 
538
 
  for( unsigned int cb2 = 0; cb2 < 2; cb2++ ) 
539
 
  for( unsigned int cb3 = 0; cb3 < 2; cb3++ ) 
540
 
  {
541
 
    chebycoefs[0]+=measure3D[cb1][cb2][cb3]*inds[cb1]*a1norm;
542
 
    chebycoefs[1]+=measure3D[cb1][cb2][cb3]*inds[cb2]*a1norm;
543
 
    chebycoefs[2]+=measure3D[cb1][cb2][cb3]*inds[cb3]*a1norm;
544
 
  }
545
 
  }
 
470
  if (ImageDimension==2)
 
471
    {
 
472
    
 
473
    double measure[3][3];
 
474
    for(int row=-1; row< 2; row++)
 
475
      {
 
476
      for(int col=-1; col< 2; col++)
 
477
        {
 
478
      
 
479
        temp[0]=tindex[0]+(long)row;
 
480
        temp[1]=tindex[1]+(long)col;
 
481
        
 
482
        for (unsigned int i=0; i<ImageDimension; i++)
 
483
          {
 
484
          if (temp[i] > m_TarSize[i]-1)
 
485
            {
 
486
            temp[i]=m_TarSize[i]-1;
 
487
            }
 
488
          else if (temp[i] < 0 )
 
489
            {
 
490
            temp[i]=0;
 
491
            }
 
492
          }
 
493
 
 
494
        requestedRegion.SetIndex(temp);
 
495
        requestedRegion.SetSize(regionRadius);
 
496
        m_TarImage->SetRequestedRegion(requestedRegion);  
 
497
        m_Metric->SetFixedImageRegion( m_TarImage->GetRequestedRegion() );
 
498
        measure[row+1][col+1]=0.0;
 
499
        
 
500
        try
 
501
          { 
 
502
          measure[row+1][col+1]=m_Metric->GetValue( parameters);
 
503
          }
 
504
        catch( ... )
 
505
          {
 
506
          }
 
507
 
 
508
    
 
509
        datatotal+=measure[row+1][col+1];
 
510
        }
 
511
      }
 
512
    for( unsigned int cb1 = 0; cb1 < 3; cb1++ )
 
513
      { 
 
514
      for( unsigned int cb2 = 0; cb2 < 3; cb2++ ) 
 
515
        {
 
516
        met=measure[cb1][cb2];
 
517
        ind1=inds[cb1]*a1norm;
 
518
        ind2=inds[cb2]*a1norm;
 
519
        chebycoefs[0]+=met*ind1;
 
520
        chebycoefs[1]+=met*ind2;
 
521
        }
 
522
      }
 
523
    }
 
524
  else if (ImageDimension == 3)
 
525
    {
 
526
 
 
527
    double measure3D[3][3][3];
 
528
    for(int row=-1; row< 2; row++)
 
529
      {
 
530
      for(int col=-1; col< 2; col++)
 
531
        {
 
532
        for(int z=-1; z< 2; z++)
 
533
          {
 
534
 
 
535
          temp[0]=tindex[0]+(long)row;
 
536
          temp[1]=tindex[1]+(long)col;
 
537
          temp[2]=tindex[2]+(long)z;
 
538
 
 
539
          for (unsigned int i=0; i<ImageDimension; i++)
 
540
            {
 
541
            if (temp[i] > m_TarSize[i]-1)
 
542
              {
 
543
              temp[i]=m_TarSize[i]-1;
 
544
              }
 
545
            else if (temp[i] < 0 )
 
546
              {
 
547
              temp[i]=0;
 
548
              }
 
549
            }
 
550
 
 
551
          requestedRegion.SetIndex(temp);
 
552
          requestedRegion.SetSize(regionRadius);
 
553
          m_TarImage->SetRequestedRegion(requestedRegion);  
 
554
          m_Metric->SetFixedImageRegion( m_TarImage->GetRequestedRegion() );
 
555
          measure3D[row+1][col+1][z+1]=0.0;
 
556
   
 
557
          try
 
558
            { 
 
559
            measure3D[row+1][col+1][z+1]=m_Metric->GetValue( parameters);
 
560
            }
 
561
          catch( ... )
 
562
            {
 
563
            }
 
564
 
 
565
    
 
566
          datatotal+=measure3D[row+1][col+1][z+1];
 
567
          }
 
568
        }
 
569
      }
 
570
    for( unsigned int cb1 = 0; cb1 < 2; cb1++ )
 
571
      {
 
572
      for( unsigned int cb2 = 0; cb2 < 2; cb2++ )
 
573
        { 
 
574
        for( unsigned int cb3 = 0; cb3 < 2; cb3++ ) 
 
575
          {
 
576
          chebycoefs[0]+=measure3D[cb1][cb2][cb3]*inds[cb1]*a1norm;
 
577
          chebycoefs[1]+=measure3D[cb1][cb2][cb3]*inds[cb2]*a1norm;
 
578
          chebycoefs[2]+=measure3D[cb1][cb2][cb3]*inds[cb3]*a1norm;
 
579
          }
 
580
        }
 
581
      }
 
582
    }
546
583
  
547
584
  chebycoefs0=a0norm*datatotal;
548
585
//  std::cout << " cb " << chebycoefs << std::endl;
555
592
int ImageMetricLoad<TMoving,TFixed>::CLID()
556
593
{
557
594
  static const int CLID_ = FEMOF::Register( ImageMetricLoad::NewB,(std::string("ImageMetricLoad(")
558
 
                +typeid(TMoving).name()+","+typeid(TFixed).name()+")").c_str());
 
595
                                                                   +typeid(TMoving).name()+","+typeid(TFixed).name()+")").c_str());
559
596
  return CLID_;
560
597
}
561
598