text:
/**
 * @file   Metrics.h
 * @author Tobias Gass <tobiasgass@gmail.com>
 * @date   Thu Mar 12 14:26:25 2015
 *
 * @brief  Collection of image comparison metrics as static functions
 *
 *
 */

#pragma once

#include "FilterUtils.hpp"
#include "TransformationUtils.h"
#include "itkRecursiveGaussianImageFilter.h"
#include <itkMeanImageFilter.h>
#include "itkMeanSquaresImageToImageMetric.h"
#include "itkMattesMutualInformationImageToImageMetric.h"
//#include <itkBoxMeanImageFilter.h>
#include "itkSubtractAbsImageFilter.h"
#include <itkAbsoluteValueDifferenceImageFilter.h>
#ifdef WITH_MIND
#include "dataCostSSC.h"
#include "dataCostLCC.h"
#endif

template<class ImageType, class InputImageType>
class MultiThreadedLocalSimilarityNCC :public itk::ImageToImageFilter < ImageType, ImageType > {
public:
  /** Standard class typedefs. */
  typedef MultiThreadedLocalSimilarityNCC             Self;
  typedef itk::ImageToImageFilter< ImageType, ImageType > Superclass;
  typedef itk::SmartPointer< Self >        Pointer;
  typedef typename ImageType::RegionType RegionType;
  typedef typename itk::ImageRegionIteratorWithIndex<ImageType> ImageIteratorType;
  typedef typename itk::ImageRegionConstIterator<InputImageType> ImageConstIteratorType;
  typedef typename ImageType::IndexType IndexType;
  typedef typename ImageType::RegionType RegionType;
  /** Method for creation through the object factory. */
  itkNewMacro(Self);

  /** Run-time type information (and related methods). */
  itkTypeMacro(MultiThreadedLocalSimilarityNCC, ImageToImageFilter);

  //at each 'pixel' of this image will the local similarity be computed.
  //the radius of the block for block matching is equal to the spacing of this image
  void SetCoarseImage(const ImageType * image){
    this->SetNthInput(0, const_cast<ImageType*>(image));
  }
  void SetImage1(const InputImageType * image) {
    this->SetNthInput(1, const_cast<InputImageType*>(image));
  }
  void SetImage2(const InputImageType * image) {
    this->SetNthInput(2, const_cast<InputImageType*>(image));
  }
  //optional mask
  //metric will only be computed for pixel where the mask is NOT NULL
  void SetMask(const InputImageType * image) {
    this->SetNthInput(3, const_cast<InputImageType*>(image));
    m_haveMask = true;
  }
  virtual void VerifyInputInformation()
  {
    //no checks so far :o
    //should check:
    //image 1 and image 2 and mask shall have same resolution and phsical overlap
    //coarse image and image1/2 overlap physically

  }

protected:
  MultiThreadedLocalSimilarityNCC()
  {
    this->SetNumberOfRequiredInputs(3);
    m_haveMask = false;

  }
  ~MultiThreadedLocalSimilarityNCC(){}
  virtual void BeforeThreadedGenerateData(){

    typename ImageType::Pointer output = this->GetOutput();
    output->SetRegions(this->GetInput(0)->GetLargestPossibleRegion());
    output->Allocate();
    output->SetSpacing(this->GetInput(0)->GetSpacing());
    output->SetOrigin(this->GetInput(0)->GetOrigin());
    output->SetDirection(this->GetInput(0)->GetDirection());
    for (int d = 0; d < ImageType::ImageDimension; ++d){
      this->m_radius[d] = (this->GetInput(0)->GetSpacing()[d] / this->GetInput(1)->GetSpacing()[d]);
    }
    m_maxSize = this->GetInput(1)->GetLargestPossibleRegion().GetSize();
  }
  virtual void ThreadedGenerateData(const RegionType & region, itk::ThreadIdType threadId)
  {

    typename ImageType::ConstPointer input = this->GetInput(0);
    typename InputImageType::ConstPointer image1 = this->GetInput(1);
    typename ImageType::Pointer output = this->GetOutput();
    ImageIteratorType outIt(output, region);
    outIt.GoToBegin();
    for (; !outIt.IsAtEnd(); ++outIt){
      typename ImageType::PointType point;
      IndexType targetIndex;
      input->TransformIndexToPhysicalPoint(outIt.GetIndex(), point);
      image1->TransformPhysicalPointToIndex(point, targetIndex);
      outIt.Set(computeLocalPotential(targetIndex));
    }
  }
protected:
  typename ImageType::SizeType m_radius;
  typename ImageType::SizeType m_maxSize;
  bool m_haveMask;
  MultiThreadedLocalSimilarityNCC(const Self &); //purposely not implemented
  void operator=(const Self &);  //purposely not implemented
  inline RegionType computeRegion(const IndexType & targetIndex){
    IndexType newIndex;
    typename ImageType::SizeType localRegionSize;
    //convert newIndex to corner of patch, assumind newIndex is the central pixel
    int numberOfPixels = 1;
    for (int d = 0; d < ImageType::ImageDimension; ++d){

      typename IndexType::IndexValueType &idx = newIndex[d];
      //get max size;
      int maxSize = this->m_radius[d];// 2 * int(this->m_coarseImageSpacing[d] / this->m_scaledTargetImage->GetSpacing()[d]);
      //LOGV(2) << VAR(maxSize) << " " << VAR(this->m_scaledRadius) << " " << VAR(this->m_radius) << std::endl;

      //move targetIndex from center of patch to top left corner of region
      idx = targetIndex[d] - maxSize;

      //extend radius to diameter
      maxSize *= 2;
      //theoretical maximum number of pixels
      numberOfPixels *= maxSize;

      //check if out of bounds and adjust region index and size accordingly
      if (idx < 0){
        int diff = newIndex[d];
        maxSize += diff;
        idx = 0;
      }
      const int extent = idx + maxSize;
      const int & maxExtent = m_maxSize[d];
      if (extent >= maxExtent){
        int diff = extent - maxExtent + 1;
        maxSize -= diff;
      }
      localRegionSize[d] = maxSize;

    }
    RegionType region;
    region.SetSize(localRegionSize);
    region.SetIndex(newIndex);
    return region;
  }
  inline virtual double computeLocalPotential(const IndexType & targetIndex){

    RegionType region = computeRegion(targetIndex);

    ImageConstIteratorType targetIt(this->GetInput(1), region);
    ImageConstIteratorType atlasIt(this->GetInput(2), region);
    targetIt.GoToBegin(); atlasIt.GoToBegin();
    ImageConstIteratorType  maskIt;
    if (m_haveMask){
      maskIt =ImageConstIteratorType(this->GetInput(3),region);
      maskIt.GoToBegin();
    }
    //return 1;
    double result;

    int insideCount = 0;
    int count = 0;
    double sff = 0.0, smm = 0.0, sfm = 0.0, sf = 0.0, sm = 0.0;
    double f, m;
    for (; !targetIt.IsAtEnd(); ++targetIt, ++atlasIt){
            ++count;
            if (this->m_haveMask){
                  float mask = maskIt.Get();
                  ++maskIt;
                  if (! (mask > 0) ) continue;
            }
            f = targetIt.Get();
            m = atlasIt.Get();
            sff += f*f;
            smm += m*m;
            sfm += f*m;
            sf += f;
            sm += m;

            insideCount += 1;

    }
    double NCC = 0;
    if (insideCount){
            sff -= (sf * sf / insideCount);
            smm -= (sm * sm / insideCount);
            sfm -= (sf * sm / insideCount);
            if (smm*sff > 0){
                  NCC = 1.0*sfm / sqrt(smm*sff);
            }
    }

    result = (1.0 - (NCC)) ;

    return result;

  }

};//MultiThreadedNCC

template<class ImageType, class InputImageType>
class MultiThreadedLocalSimilaritySSD :public MultiThreadedLocalSimilarityNCC<ImageType,InputImageType>{
public:
  /** Standard class typedefs. */
  typedef MultiThreadedLocalSimilaritySSD             Self;
  typedef MultiThreadedLocalSimilarityNCC<ImageType, InputImageType> Superclass;
  typedef itk::SmartPointer< Self >        Pointer;
  typedef typename ImageType::RegionType RegionType;
  typedef typename itk::ImageRegionIteratorWithIndex<ImageType> ImageIteratorType;
  typedef typename itk::ImageRegionConstIterator<InputImageType> ImageConstIteratorType;
  typedef typename ImageType::IndexType IndexType;
  /** Method for creation through the object factory. */
  itkNewMacro(Self);

  /** Run-time type information (and related methods). */
  itkTypeMacro(MultiThreadedLocalSimilaritySSD, MultiThreadedLocalSimilarityNCC);

protected:

  inline virtual double computeLocalPotential(const IndexType & targetIndex){
    RegionType region = computeRegion(targetIndex);

    ImageConstIteratorType targetIt(this->GetInput(1), region);
    ImageConstIteratorType atlasIt(this->GetInput(2), region);
    targetIt.GoToBegin(); atlasIt.GoToBegin();
    ImageConstIteratorType  maskIt;
    if (m_haveMask){
      maskIt = ImageConstIteratorType(this->GetInput(3), region);
      maskIt.GoToBegin();
    }
    //return 1;
    double result;

    int insideCount = 0.0;
    int count = 0;
    double f, m;
    double ssd = 0.0;
    for (; !targetIt.IsAtEnd(); ++targetIt, ++atlasIt){
      ++count;
      if (m_haveMask){
        float mask = maskIt.Get();
        ++maskIt;
        if (mask > 0) continue;
      }
      f = targetIt.Get();
      m = atlasIt.Get();
      double diff = f - m;

      ssd += diff*diff;

      insideCount += 1;

    }

    if (count){
      result = ssd / count;
    }
    return result;
  }

};//MultiThreadedSSD

///class containing static functions to compute image similarities with
template<typename InputImage, typename OutputImage = InputImage, typename InternalPrecision=double>
class Metrics{

    typedef typename InputImage::Pointer  InputImagePointer;
    typedef typename InputImage::ConstPointer  ConstInputImagePointer;

    typedef typename InputImage::PixelType InputImagePixelType;
    typedef typename InputImage::IndexType InputImageIndex;
    typedef typename InputImage::RegionType InputImageRegion;
    typedef typename InputImage::SpacingType SpacingType;

    typedef typename OutputImage::Pointer OutputImagePointer;
    typedef typename OutputImage::IndexType OutputImageIndex;
    typedef typename OutputImage::RegionType OutputImageRegion;
    typedef typename OutputImage::PixelType OutputImagePixelType;
    typedef typename OutputImage::PointType OutputImagePointType;
    static const int D=InputImage::ImageDimension;
    typedef typename TransfUtils<InputImage>::DeformationFieldType DeformationFieldType;
    typedef typename DeformationFieldType::Pointer DeformationFieldPointerType;

    typedef typename itk::Image<InternalPrecision,D> InternalImage;
    typedef typename InternalImage::Pointer InternalImagePointer;

public:

    //hellishly inefficient but "clean" implementation of local normalized cross correlation
    static inline OutputImagePointer LNCC(InputImagePointer i1,InputImagePointer i2,double sigma=1.0, double exp = 1.0){
        return LNCC( (ConstInputImagePointer)i1, (ConstInputImagePointer)i2, sigma,exp);
    }
    static inline OutputImagePointer LNCC(ConstInputImagePointer i1,ConstInputImagePointer i2,double sigma=1.0, double exp = 1.0){
        InternalImagePointer i1Cast=FilterUtils<InputImage,InternalImage>::cast(i1);
        InternalImagePointer i2Cast=FilterUtils<InputImage,InternalImage>::cast(i2);
        typedef typename itk::SmoothingRecursiveGaussianImageFilter< InternalImage, InternalImage > FilterType;
        //typedef itk::DiscreteGaussianImageFilter<InternalImage,InternalImage>  FilterType;
        typename FilterType::Pointer filter=FilterType::New();
        filter->SetSigma(sigma);
        //filter->SetVariance(sigma*sigma);
        //        filter->SetMaximumKernelWidth(sigma);

        //compute local means by concolving with gaussian
        filter->SetInput(i1Cast);
        filter->Update();
        InternalImagePointer i1Bar=ImageUtils<InternalImage>::duplicate(filter->GetOutput());
        LOGI(8,ImageUtils<InternalImage>::writeImage("i1Bar.mha",i1Bar));
        filter->SetInput(i2Cast);
        filter->Update();
        InternalImagePointer i2Bar=ImageUtils<InternalImage>::duplicate(filter->GetOutput());
        FilterUtils<InternalImage,InternalImage>::lowerThresholding(i2Bar,std::numeric_limits<InputImagePixelType>::min());
        LOGI(8,ImageUtils<InternalImage>::writeImage("i2Bar.mha",i2Bar));
        //HACK!
        //compute squares of original images
        InternalImagePointer i1Square=ImageUtils<InternalImage>::multiplyImageOutOfPlace(i1Cast,i1Cast);
        InternalImagePointer i2Square=ImageUtils<InternalImage>::multiplyImageOutOfPlace(i2Cast,i2Cast);
        LOGI(8,ImageUtils<InternalImage>::writeImage("i2Square.mha",i2Square));
        LOGI(8,ImageUtils<InternalImage>::writeImage("i1Square.mha",i1Square));

        //compute local means of squared images by convolving with gaussian kernel
        filter->SetInput(i1Square);
        filter->Update();
        InternalImagePointer i1SquareBar=ImageUtils<InternalImage>::duplicate(filter->GetOutput());
        LOGI(8,ImageUtils<InternalImage>::writeImage("i1SquareBar.mha",i1SquareBar));
        filter->SetInput(i2Square);
        filter->Update();
        InternalImagePointer i2SquareBar=ImageUtils<InternalImage>::duplicate(filter->GetOutput());
        FilterUtils<InternalImage,InternalImage>::lowerThresholding(i2SquareBar,0);
        LOGI(8,ImageUtils<InternalImage>::writeImage("i2SquareBar.mha",i2SquareBar));

        //compute squares of local means
        InternalImagePointer i1BarSquare=ImageUtils<InternalImage>::localSquare(i1Bar);
        LOGI(8,ImageUtils<InternalImage>::writeImage("i1BarSquare.mha",i1BarSquare));
        InternalImagePointer i2BarSquare=ImageUtils<InternalImage>::multiplyImageOutOfPlace(i2Bar,i2Bar);
        LOGI(8,ImageUtils<InternalImage>::writeImage("i2BarSquare.mha",i2BarSquare));

        //multiply i1 and i2 locally
        InternalImagePointer i1i2=ImageUtils<InternalImage>::multiplyImageOutOfPlace(i1Cast,i2Cast);
        LOGI(8,ImageUtils<InternalImage>::writeImage("i1i2.mha",i1i2));
        //compute local means by convolving...
        filter->SetInput(i1i2);
        filter->Update();
        InternalImagePointer i1Timesi2Bar=ImageUtils<InternalImage>::duplicate(filter->GetOutput());
        LOGI(8,ImageUtils<InternalImage>::writeImage("i1Timesi2Bar.mha",i1Timesi2Bar));

        //multiply local means
        InternalImagePointer i1BarTimesi2Bar=ImageUtils<InternalImage>::multiplyImageOutOfPlace(i1Bar,i2Bar);
        LOGI(8,ImageUtils<InternalImage>::writeImage("i1BarTimesi2Bar.mha",i1BarTimesi2Bar));

        InternalImagePointer numerator=FilterUtils<InternalImage>::substract(i1Timesi2Bar,i1BarTimesi2Bar);
        LOGI(8,ImageUtils<InternalImage>::writeImage("numerator.mha",numerator));

        InternalImagePointer varianceI1=FilterUtils<InternalImage>::substract(i1SquareBar,i1BarSquare);
        FilterUtils<InternalImage,InternalImage>::lowerThresholding(varianceI1,std::numeric_limits<InternalPrecision>::epsilon());
        LOGI(8,ImageUtils<InternalImage>::writeImage("varianceI1.mha",varianceI1));

        InternalImagePointer varianceI2=FilterUtils<InternalImage>::substract(i2SquareBar,i2BarSquare);
        FilterUtils<InternalImage,InternalImage>::lowerThresholding(varianceI2,std::numeric_limits<InternalPrecision>::epsilon());

        LOGI(8,ImageUtils<InternalImage>::writeImage("varianceI2.mha",varianceI2));

        ImageUtils<InternalImage>::sqrtImage(varianceI1);
        ImageUtils<InternalImage>::sqrtImage(varianceI2);

        InternalImagePointer denominator=ImageUtils<InternalImage>::multiplyImageOutOfPlace(varianceI1,varianceI2);

        InternalImagePointer result=ImageUtils<InternalImage>::divideImageOutOfPlace(numerator,denominator);
        //convert to weights

        itk::ImageRegionIterator<InternalImage> it(result,result->GetLargestPossibleRegion());
        it.GoToBegin();
        for (;!it.IsAtEnd();++it){
            InternalPrecision r = it.Get();
            if (r< -(InternalPrecision)1.0 ){
                r = 0;
            }else if ( r > (InternalPrecision)1.0){
                r = 0;
            }
            r = pow((r+1.0)/2,exp);

            r = max(r,std::numeric_limits<InternalPrecision>::epsilon());
            it.Set(r);

        }
        // LOGI(8,ImageUtils<InternalImage>::add(result,1.0));
        //        ImageUtils<InternalImage>::multiplyImage(result,0.5);
        //        FilterUtils<InternalImage,InternalImage>::thresholding(result,std::numeric_limits<InternalPrecision>::epsilon(),1.0);

        LOGI(10,ImageUtils<InternalImage>::writeImage("result.mha",result));

        return FilterUtils<InternalImage,OutputImage>::cast(result);

    }

    static inline OutputImagePointer LSSDNorm(InputImagePointer i1,InputImagePointer i2,double sigmaWidth=1.0, double sigmaNorm=1.0){
        InternalImagePointer i1Cast=FilterUtils<InputImage,InternalImage>::cast(i1);
        InternalImagePointer i2Cast=FilterUtils<InputImage,InternalImage>::cast(i2);
        typedef typename itk::SmoothingRecursiveGaussianImageFilter< InternalImage, InternalImage > FilterType;
        //typedef itk::DiscreteGaussianImageFilter<InternalImage,InternalImage>  FilterType;
        typename FilterType::Pointer filter=FilterType::New();
        filter->SetSigma(sigmaWidth);
        //filter->SetVariance(sigmaWidth*sigmaWidth);

        InternalImagePointer diff = FilterUtils<InternalImage>::substract(i1Cast,i2Cast);
        InternalImagePointer diffSquare = ImageUtils<InternalImage>::localSquare(diff);

        if (sigmaNorm==0.0){
            itk::ImageRegionIterator<InternalImage> it(diffSquare,diffSquare->GetLargestPossibleRegion());
            itk::ImageRegionIterator<InternalImage> it2(i1Cast,diffSquare->GetLargestPossibleRegion());
            it.GoToBegin();it2.GoToBegin();
            InternalPrecision maxVal=FilterUtils<InternalImage>::getMax(diffSquare);
            sigmaNorm=-0.005*maxVal/log(0.1);
            LOGV(6)<<VAR(sqrt(sigmaNorm))<<endl;
#if 0
            for (;!it.IsAtEnd();++it,++it2){
                if (true || it2.Get()>minVal){
                    sigmaNorm+=it.Get();
                    ++c;
                }
            }
            sigmaNorm/=c;
            LOGV(6)<<VAR(sigmaNorm)<<endl;
#endif
            //sigmaNorm=(FilterUtils<InternalImage>::getMean(diffSquare));
        }else{
            sigmaNorm*=sigmaNorm;
        }

        //compute local means by concolving with gaussian
        filter->SetInput(diffSquare);
        filter->Update();
        InternalImagePointer i1Bar=ImageUtils<InternalImage>::duplicate(filter->GetOutput());
        LOGI(8,ImageUtils<InternalImage>::writeImage("lssd.nii",i1Bar));
        typename ImageUtils<InternalImage>::ImageIteratorType meanIt(i1Bar,i1Cast->GetLargestPossibleRegion());

        InternalImagePointer result = ImageUtils<InternalImage>::createEmpty(i1Bar);

        typename ImageUtils<InternalImage>::ImageIteratorType resultIt(result,result->GetLargestPossibleRegion());

        for (resultIt.GoToBegin(), meanIt.GoToBegin(); !resultIt.IsAtEnd(); ++resultIt, ++meanIt){
            InternalPrecision d = meanIt.Get();
            resultIt.Set(exp(-0.5 * fabs(d)  / sigmaNorm) );
        }

        return FilterUtils<InternalImage,OutputImage>::cast(result);
    }

    static inline OutputImagePointer LSSD(InputImagePointer i1,InputImagePointer i2,double sigmaWidth=1.0){
        if (sigmaWidth==0.0) sigmaWidth=0.001;
        return LSSD( (ConstInputImagePointer)i1, (ConstInputImagePointer)i2, sigmaWidth);
    }
    static inline OutputImagePointer LSSD(ConstInputImagePointer i1,ConstInputImagePointer i2,double sigmaWidth=1.0){
        if (sigmaWidth==0.0) sigmaWidth=0.001;
        InternalImagePointer i1Cast=FilterUtils<InputImage,InternalImage>::cast(i1);
        InternalImagePointer i2Cast=FilterUtils<InputImage,InternalImage>::cast(i2);

        //InternalImagePointer diff = FilterUtils<InternalImage>::substract(i1Cast,i2Cast);
        //InternalImagePointer diffSquare = ImageUtils<InternalImage>::localSquare(diff);
        typename ImageUtils<InternalImage>::ImageIteratorType i1It(i1Cast,i1Cast->GetLargestPossibleRegion());
        typename ImageUtils<InternalImage>::ImageIteratorType i2It(i2Cast,i1Cast->GetLargestPossibleRegion());
        for (i2It.GoToBegin(),i1It.GoToBegin(); i1It.IsAtEnd(); ++i1It, ++i2It){
            InternalPrecision d=(i1It.Get()-i2It.Get());
            i1It.Set(d*d);
        }
#if 1
        typedef typename itk::SmoothingRecursiveGaussianImageFilter< InternalImage, InternalImage > FilterType;
        typename FilterType::Pointer filter=FilterType::New();
        filter->SetSigma(sigmaWidth);
        //compute local means by concolving with gaussian
        filter->SetInput(i1Cast);
        filter->Update();
        return (filter->GetOutput());
#elif 1
        typedef typename itk::BoxMeanImageFilter< InternalImage, InternalImage > FilterType;
        typename FilterType::Pointer filter=FilterType::New();
        typename FilterType::RadiusType r;
        r.Fill(sigmaWidth);
        filter->SetRadius(r);
        //compute local means by concolving with gaussian
        filter->SetInput(i1Cast);
        return filter->GetOutput();

#endif
    }

    static inline OutputImagePointer integralSSD(InputImagePointer i1,InputImagePointer i2){
        InternalImagePointer result=FilterUtils<InputImage,InternalImage>::cast(i1);

        //InternalImagePointer diff = FilterUtils<InternalImage>::substract(i1Cast,i2Cast);
        //InternalImagePointer diffSquare = ImageUtils<InternalImage>::localSquare(diff);
        typename ImageUtils<InternalImage>::ImageIteratorType i1It(result,result->GetLargestPossibleRegion());
        typename ImageUtils<InputImage>::ImageIteratorType i2It(i2,result->GetLargestPossibleRegion());
        InternalPrecision sum=0.0;
        for (i2It.GoToBegin(),i1It.GoToBegin(); i1It.IsAtEnd(); ++i1It, ++i2It){
            InternalPrecision d=(i1It.Get()-i2It.Get());
            sum+=d*d;
            i1It.Set(sum);
        }
        return result;
    }

    static inline OutputImagePointer integralSAD(InputImagePointer i1,InputImagePointer i2){
        InternalImagePointer result=FilterUtils<InputImage,InternalImage>::cast(i1);

        //InternalImagePointer diff = FilterUtils<InternalImage>::substract(i1Cast,i2Cast);
        //InternalImagePointer diffSquare = ImageUtils<InternalImage>::localSquare(diff);
        typename ImageUtils<InternalImage>::ImageIteratorType i1It(result,result->GetLargestPossibleRegion());
        typename ImageUtils<InputImage>::ImageIteratorType i2It(i2,result->GetLargestPossibleRegion());
        InternalPrecision sum=0.0;
        for (i2It.GoToBegin(),i1It.GoToBegin(); i1It.IsAtEnd(); ++i1It, ++i2It){
            InternalPrecision d=(i1It.Get()-i2It.Get());
            sum+=fabs(d);
            i1It.Set(sum);
        }
        return result;
    }

    static inline InternalImagePointer LSAD(InputImagePointer i1,InputImagePointer i2,double sigmaWidth=1.0){
        if (sigmaWidth==0.0) sigmaWidth=0.001;
        return LSAD( (ConstInputImagePointer)i1, (ConstInputImagePointer)i2, sigmaWidth);
    }
  static inline InternalImagePointer LSAD(ConstInputImagePointer i1,ConstInputImagePointer i2,double sigmaWidth=1.0){
        if (sigmaWidth==0.0) sigmaWidth=0.001;
        typedef typename itk::SmoothingRecursiveGaussianImageFilter< InternalImage, InternalImage > FilterType;
        //typedef itk::DiscreteGaussianImageFilter<InternalImage,InternalImage>  FilterType;
        typename FilterType::Pointer filter=FilterType::New();
        filter->SetSigma(sigmaWidth);
        //filter->SetVariance(sigmaWidth*sigmaWidth);

#if 0
        InternalImagePointer i1Cast=FilterUtils<InputImage,InternalImage>::cast(i1);
        InternalImagePointer i2Cast=FilterUtils<InputImage,InternalImage>::cast(i2);

        InternalImagePointer diff =  ImageUtils<InternalImage>::createEmpty(i1Cast);
        typename ImageUtils<InternalImage>::ImageIteratorType diffIt(diff,i1Cast->GetLargestPossibleRegion());
        typename ImageUtils<InternalImage>::ImageIteratorType i1It(i1Cast,i1Cast->GetLargestPossibleRegion());
        typename ImageUtils<InternalImage>::ImageIteratorType i2It(i2Cast,i1Cast->GetLargestPossibleRegion());

        InternalPrecision mean=0.0;
        int c=0;
        //subtract and take absolute value
        for (i2It.GoToBegin(),i1It.GoToBegin(), diffIt.GoToBegin(); !diffIt.IsAtEnd(); ++i1It, ++i2It, ++diffIt,++c){
            InternalPrecision d=fabs(i1It.Get()-i2It.Get());
            diffIt.Set(d);
            mean+=d;
        }
        LOGV(6)<<VAR(mean/c)<<endl;

#else
        //typedef typename itk::SubtractAbsImageFilter<InputImage,InputImage,InternalImage> SubFilterType;

        typedef typename itk::AbsoluteValueDifferenceImageFilter<InputImage,InputImage,InternalImage> SubFilterType;
        typename SubFilterType::Pointer subFilter=SubFilterType::New();
        subFilter->SetNumberOfThreads( 8 );
        subFilter->SetInput1(i1);
        subFilter->SetInput2(i2);
        subFilter->Update();
        InternalImagePointer diff=subFilter->GetOutput();
        return FilterUtils<InternalImage,OutputImage>::cast(diff);
#endif

        //compute local means by concolving with gaussian
        filter->SetInput(diff);
        filter->Update();
        InternalImagePointer result=filter->GetOutput();
#if 0
        InternalImagePointer i1Bar=ImageUtils<InternalImage>::duplicate(filter->GetOutput());

        InternalImagePointer result = ImageUtils<InternalImage>::createEmpty(i1Bar);

        typename ImageUtils<InternalImage>::ImageIteratorType resultIt(result,result->GetLargestPossibleRegion());

        for (resultIt.GoToBegin(), diffIt.GoToBegin(); !resultIt.IsAtEnd(); ++resultIt, ++diffIt){
            InternalPrecision d = diffIt.Get();
            resultIt.Set(fabs(d));
        }
#endif

        return FilterUtils<InternalImage,OutputImage>::cast(result);
  }

    static inline OutputImagePointer LSADNorm(InputImagePointer i1,InputImagePointer i2,double sigmaWidth=1.0, double sigmaNorm=1.0){
        if (sigmaWidth==0.0) sigmaWidth=0.001;
        InternalImagePointer i1Cast=FilterUtils<InputImage,InternalImage>::cast(i1);
        InternalImagePointer i2Cast=FilterUtils<InputImage,InternalImage>::cast(i2);
        typedef typename itk::SmoothingRecursiveGaussianImageFilter< InternalImage, InternalImage > FilterType;
        //typedef itk::DiscreteGaussianImageFilter<InternalImage,InternalImage>  FilterType;
        typename FilterType::Pointer filter=FilterType::New();
        filter->SetSigma(sigmaWidth);
        //filter->SetVariance(sigmaWidth*sigmaWidth);
        InternalImagePointer diff =  ImageUtils<InternalImage>::createEmpty(i1Cast);
        typename ImageUtils<InternalImage>::ImageIteratorType diffIt(diff,i1Cast->GetLargestPossibleRegion());
        typename ImageUtils<InternalImage>::ImageIteratorType i1It(i1Cast,i1Cast->GetLargestPossibleRegion());
        typename ImageUtils<InternalImage>::ImageIteratorType i2It(i2Cast,i1Cast->GetLargestPossibleRegion());

        InternalPrecision mean=0.0;
        int c=0;
        //subtract and take absolute value
        for (i2It.GoToBegin(),i1It.GoToBegin(), diffIt.GoToBegin(); !diffIt.IsAtEnd(); ++i1It, ++i2It, ++diffIt,++c){
            InternalPrecision d=fabs(i1It.Get()-i2It.Get());
            diffIt.Set(d);
            mean+=d;
        }
        mean/=c;
        LOGV(6)<<VAR(mean)<<endl;

        if (sigmaNorm==0.0){
            sigmaNorm=mean;
        }

        //compute local means by concolving with gaussian
        filter->SetInput(diff);
        filter->Update();
        InternalImagePointer i1Bar=ImageUtils<InternalImage>::duplicate(filter->GetOutput());
        typename ImageUtils<InternalImage>::ImageIteratorType meanIt(i1Bar,i1Cast->GetLargestPossibleRegion());

        InternalImagePointer result = ImageUtils<InternalImage>::createEmpty(i1Bar);

        typename ImageUtils<InternalImage>::ImageIteratorType resultIt(result,result->GetLargestPossibleRegion());

        for (resultIt.GoToBegin(), meanIt.GoToBegin(); !resultIt.IsAtEnd(); ++resultIt, ++meanIt){
            InternalPrecision d = meanIt.Get();
            //resultIt.Set(max(0.5,exp(-0.5 * fabs(d)  / sigmaNorm) ));
            resultIt.Set(exp(-0.5 * fabs(d)  / sigmaNorm) );
            //resultIt.Set(pow(exp(-0.5 * fabs(d)  / mean ),sigmaNorm ));
        }

        return FilterUtils<InternalImage,OutputImage>::cast(result);
    }

    static inline OutputImagePointer efficientLNCC(InputImagePointer i1,InputImagePointer i2,double sigma=1.0, double exp = 1.0){
        return efficientLNCC( (ConstInputImagePointer)i1, (ConstInputImagePointer)i2, sigma,exp);
    }
    static inline OutputImagePointer efficientLNCC(ConstInputImagePointer i1,ConstInputImagePointer i2,double sigma=1.0, double exp = 1.0){
        //if (exp == 0.0 ) exp = 1.0;
        if (sigma==0.0) sigma=0.001;
        InternalImagePointer i1Cast=FilterUtils<InputImage,InternalImage>::cast(i1);
        InternalImagePointer i2Cast=FilterUtils<InputImage,InternalImage>::cast(i2);
#define RECURSIVE

#ifdef RECURSIVE
        typedef typename itk::SmoothingRecursiveGaussianImageFilter< InternalImage, InternalImage > FilterType;
        //typedef typename itk::RecursiveGaussianImageFilter< InternalImage, InternalImage > FilterType;
#else
        typedef itk::DiscreteGaussianImageFilter<InternalImage,InternalImage>  FilterType;
        //typedef typename itk::MeanImageFilter< InternalImage, InternalImage > FilterType;

#endif
        typename FilterType::Pointer filter=FilterType::New();
#ifdef RECURSIVE
        filter->SetSigma(sigma);

        LOGV(5)<<VAR(sigma)<<endl;
#else
        filter->SetVariance(sigma*sigma);
        //filter->SetRadius(sigma);
#endif
        //        filter->SetMaximumKernelWidth(sigma);

        //compute local means by concolving with gaussian
        filter->SetInput(i1Cast);
        //filter->InPlaceOff();
        filter->Update();
        InternalImagePointer i1Bar=ImageUtils<InternalImage>::duplicate(filter->GetOutput());
        LOGI(8,ImageUtils<InternalImage>::writeImage("i1Bar.mha",i1Bar));
        filter->SetInput(i2Cast);
        filter->UpdateLargestPossibleRegion();
        InternalImagePointer i2Bar=ImageUtils<InternalImage>::duplicate(filter->GetOutput());

        //FilterUtils<InternalImage,InternalImage>::lowerThresholding(i2Bar,std::numeric_limits<InputImagePixelType>::min());
        LOGI(8,ImageUtils<InternalImage>::writeImage("i2Bar.mha",i2Bar));
        //HACK!
        //compute squares of original images

        InternalImagePointer i1Square=ImageUtils<InternalImage>::multiplyImageOutOfPlace(i1Cast,i1Cast);
        InternalImagePointer i2Square=ImageUtils<InternalImage>::multiplyImageOutOfPlace(i2Cast,i2Cast);

        LOGI(8,ImageUtils<InternalImage>::writeImage("i2Square.mha",i2Square));
        LOGI(8,ImageUtils<InternalImage>::writeImage("i1Square.mha",i1Square));

        //compute local means of squared images by convolving with gaussian kernel
        filter->SetInput(i1Square);
        filter->UpdateLargestPossibleRegion();
        InternalImagePointer i1SquareBar=ImageUtils<InternalImage>::duplicate(filter->GetOutput());
        LOGI(8,ImageUtils<InternalImage>::writeImage("i1SquareBar.mha",i1SquareBar));
        i1Square=NULL;

        filter->SetInput(i2Square);
        filter->UpdateLargestPossibleRegion();
        InternalImagePointer i2SquareBar=ImageUtils<InternalImage>::duplicate(filter->GetOutput());
        //FilterUtils<InternalImage,InternalImage>::lowerThresholding(i2SquareBar,0);
        LOGI(8,ImageUtils<InternalImage>::writeImage("i2SquareBar.mha",i2SquareBar));
        i2Square=NULL;
        //compute squares of local means
        InternalImagePointer i1BarSquare=ImageUtils<InternalImage>::localSquare(i1Bar);
        InternalImagePointer i2BarSquare=ImageUtils<InternalImage>::localSquare(i2Bar);
        //multiply i1 and i2 locally
        InternalImagePointer i1i2=ImageUtils<InternalImage>::multiplyImageOutOfPlace(i1Cast,i2Cast);
        LOGI(8,ImageUtils<InternalImage>::writeImage("i1i2.mha",i1i2));
        i1Cast=NULL;
        i2Cast=NULL;
        LOGI(8,ImageUtils<InternalImage>::writeImage("i1BarSquare.mha",i1BarSquare));
        LOGI(8,ImageUtils<InternalImage>::writeImage("i2BarSquare.mha",i2BarSquare));

        //compute local means by convolving...
        filter->SetInput(i1i2);
        filter->UpdateLargestPossibleRegion();
        InternalImagePointer i1Timesi2Bar=ImageUtils<InternalImage>::duplicate(filter->GetOutput());
        LOGI(8,ImageUtils<InternalImage>::writeImage("i1Timesi2Bar.mha",i1Timesi2Bar));
        i1i2=NULL;

        //finish
        InternalImagePointer result=ImageUtils<InternalImage>::duplicate(filter->GetOutput());
        itk::ImageRegionIterator<InternalImage> resultIt(result,result->GetLargestPossibleRegion());
        resultIt.GoToBegin();

        itk::ImageRegionIterator<InternalImage> i1BarIt(i1Bar,i1Bar->GetLargestPossibleRegion());
        i1BarIt.GoToBegin();
        itk::ImageRegionIterator<InternalImage> i2BarIt(i2Bar,i2Bar->GetLargestPossibleRegion());
        i2BarIt.GoToBegin();

        itk::ImageRegionIterator<InternalImage> i1Timesi2BarIt(i1Timesi2Bar,i1Timesi2Bar->GetLargestPossibleRegion());
        i1Timesi2BarIt.GoToBegin();

        itk::ImageRegionIterator<InternalImage> i1SquareBarIt(i1SquareBar,i1SquareBar->GetLargestPossibleRegion());
        i1SquareBarIt.GoToBegin();

        itk::ImageRegionIterator<InternalImage> i2SquareBarIt(i2SquareBar,i2SquareBar->GetLargestPossibleRegion());
        i2SquareBarIt.GoToBegin();

        for (;!resultIt.IsAtEnd();++resultIt){

            InternalPrecision i1BarV=i1BarIt.Get();
            InternalPrecision i2BarV=i2BarIt.Get();
            InternalPrecision i1BarTimesi2BarV=i1BarV*i2BarV;

            InternalPrecision i1Timesi2BarV=i1Timesi2BarIt.Get();
            InternalPrecision numeratorV=i1Timesi2BarV-i1BarTimesi2BarV;
#if 1 //def RECURSIVE
            //the recursive filtering is not very exact.
            //this is hardly noticeable when only smoothing single images, but during the computation of variances this leads to errors
            //mostly these errors happen when the variance is very low, eg in smooth regions
            //we correct this errors by thresholding, which looks fine when compared to discrete gaussian filtering.
            InternalPrecision tmp=i1SquareBarIt.Get() - i1BarV*i1BarV;

            if (tmp<0.0){
                LOGV(10)<<VAR(i1SquareBarIt.Get() - i1BarV*i1BarV)<<endl;
                tmp=0.0;
            }
            InternalPrecision varianceI1V=sqrt(tmp);
            tmp=i2SquareBarIt.Get() - i2BarV*i2BarV;
            if (tmp<0.0){
                LOGV(10)<<VAR(i2SquareBarIt.Get() - i2BarV*i2BarV)<<endl;
                tmp=0.0;
            }
            InternalPrecision varianceI2V=sqrt(tmp);
            InternalPrecision denominatorV=varianceI1V*varianceI2V;
            InternalPrecision r = (abs(denominatorV)>1000.0*std::numeric_limits<InternalPrecision>::epsilon())?numeratorV/(denominatorV):0.0;

#else
            InternalPrecision varianceI1V=sqrt((i1SquareBarIt.Get() - i1BarV*i1BarV));
            InternalPrecision varianceI2V=sqrt((i2SquareBarIt.Get() - i2BarV*i2BarV));
            InternalPrecision denominatorV=varianceI1V*varianceI2V;
            InternalPrecision r = (abs(denominatorV)>0.0)?numeratorV/(denominatorV):0.0;

#endif

            ++i1BarIt;
            ++i2BarIt;
            ++i1Timesi2BarIt;
            ++i1SquareBarIt;
            ++i2SquareBarIt;

#ifdef RECURSIVE
            if (r< -(InternalPrecision)1.0 ){
                r = -1.0;
            }else if ( r > (InternalPrecision)1.0){
                r = 1.0;
            }
#endif
            r = pow((r+1.0)/2,exp);

#ifdef SAFE
            r = max(r,std::numeric_limits<InternalPrecision>::epsilon());
#endif
            //LOGV(7)<<"LNCC-local: "<<r<<endl;
            resultIt.Set(r);

        }

        LOGI(8,ImageUtils<InternalImage>::writeImage("result.mha",result));

        return FilterUtils<InternalImage,OutputImage>::cast(result);

    }
       static inline OutputImagePointer efficientLNCCNewNorm(InputImagePointer i1,InputImagePointer i2,double sigma=1.0, double exp = 1.0){
           return efficientLNCCNewNorm(ConstInputImagePointer(i1),ConstInputImagePointer(i2), sigma, exp);
       }

       static inline OutputImagePointer efficientLNCCNewNorm(ConstInputImagePointer i1,ConstInputImagePointer i2,double sigma=1.0, double exp = 1.0){
        if (exp == 0.0 ) exp = 1.0;
        if (sigma==0.0) sigma=0.001;
        InternalImagePointer i1Cast=FilterUtils<InputImage,InternalImage>::cast(i1);
        InternalImagePointer i2Cast=FilterUtils<InputImage,InternalImage>::cast(i2);
        //#define RECURSIVE

#ifdef RECURSIVE
        typedef typename itk::SmoothingRecursiveGaussianImageFilter< InternalImage, InternalImage > FilterType;
        //typedef typename itk::RecursiveGaussianImageFilter< InternalImage, InternalImage > FilterType;
#else
        typedef itk::DiscreteGaussianImageFilter<InternalImage,InternalImage>  FilterType;
        //typedef typename itk::MeanImageFilter< InternalImage, InternalImage > FilterType;

#endif
        typename FilterType::Pointer filter=FilterType::New();
#ifdef RECURSIVE
        filter->SetSigma(sigma);

        LOGV(5)<<VAR(sigma)<<endl;
#else
        filter->SetVariance(sigma*sigma);
        //filter->SetRadius(sigma);
#endif
        //        filter->SetMaximumKernelWidth(sigma);

        //compute local means by concolving with gaussian
        filter->SetInput(i1Cast);
        //filter->InPlaceOff();
        filter->Update();
        InternalImagePointer i1Bar=ImageUtils<InternalImage>::duplicate(filter->GetOutput());
        LOGI(8,ImageUtils<InternalImage>::writeImage("i1Bar.mha",i1Bar));
        filter->SetInput(i2Cast);
        filter->UpdateLargestPossibleRegion();
        InternalImagePointer i2Bar=ImageUtils<InternalImage>::duplicate(filter->GetOutput());

        //FilterUtils<InternalImage,InternalImage>::lowerThresholding(i2Bar,std::numeric_limits<InputImagePixelType>::min());
        LOGI(8,ImageUtils<InternalImage>::writeImage("i2Bar.mha",i2Bar));
        //HACK!
        //compute squares of original images

        InternalImagePointer i1Square=ImageUtils<InternalImage>::multiplyImageOutOfPlace(i1Cast,i1Cast);
        InternalImagePointer i2Square=ImageUtils<InternalImage>::multiplyImageOutOfPlace(i2Cast,i2Cast);

        LOGI(8,ImageUtils<InternalImage>::writeImage("i2Square.mha",i2Square));
        LOGI(8,ImageUtils<InternalImage>::writeImage("i1Square.mha",i1Square));

        //compute local means of squared images by convolving with gaussian kernel
        filter->SetInput(i1Square);
        filter->UpdateLargestPossibleRegion();
        InternalImagePointer i1SquareBar=ImageUtils<InternalImage>::duplicate(filter->GetOutput());
        LOGI(8,ImageUtils<InternalImage>::writeImage("i1SquareBar.mha",i1SquareBar));
        i1Square=NULL;

        filter->SetInput(i2Square);
        filter->UpdateLargestPossibleRegion();
        InternalImagePointer i2SquareBar=ImageUtils<InternalImage>::duplicate(filter->GetOutput());
        //FilterUtils<InternalImage,InternalImage>::lowerThresholding(i2SquareBar,0);
        LOGI(8,ImageUtils<InternalImage>::writeImage("i2SquareBar.mha",i2SquareBar));
        i2Square=NULL;
        //compute squares of local means
        InternalImagePointer i1BarSquare=ImageUtils<InternalImage>::localSquare(i1Bar);
        InternalImagePointer i2BarSquare=ImageUtils<InternalImage>::localSquare(i2Bar);
        //multiply i1 and i2 locally
        InternalImagePointer i1i2=ImageUtils<InternalImage>::multiplyImageOutOfPlace(i1Cast,i2Cast);
        LOGI(8,ImageUtils<InternalImage>::writeImage("i1i2.mha",i1i2));
        i1Cast=NULL;
        i2Cast=NULL;
        LOGI(8,ImageUtils<InternalImage>::writeImage("i1BarSquare.mha",i1BarSquare));
        LOGI(8,ImageUtils<InternalImage>::writeImage("i2BarSquare.mha",i2BarSquare));

        //compute local means by convolving...
        filter->SetInput(i1i2);
        filter->UpdateLargestPossibleRegion();
        InternalImagePointer i1Timesi2Bar=ImageUtils<InternalImage>::duplicate(filter->GetOutput());
        LOGI(8,ImageUtils<InternalImage>::writeImage("i1Timesi2Bar.mha",i1Timesi2Bar));
        i1i2=NULL;

        //finish
        InternalImagePointer result=ImageUtils<InternalImage>::duplicate(filter->GetOutput());
        itk::ImageRegionIterator<InternalImage> resultIt(result,result->GetLargestPossibleRegion());
        resultIt.GoToBegin();

        itk::ImageRegionIterator<InternalImage> i1BarIt(i1Bar,i1Bar->GetLargestPossibleRegion());
        i1BarIt.GoToBegin();
        itk::ImageRegionIterator<InternalImage> i2BarIt(i2Bar,i2Bar->GetLargestPossibleRegion());
        i2BarIt.GoToBegin();

        itk::ImageRegionIterator<InternalImage> i1Timesi2BarIt(i1Timesi2Bar,i1Timesi2Bar->GetLargestPossibleRegion());
        i1Timesi2BarIt.GoToBegin();

        itk::ImageRegionIterator<InternalImage> i1SquareBarIt(i1SquareBar,i1SquareBar->GetLargestPossibleRegion());
        i1SquareBarIt.GoToBegin();

        itk::ImageRegionIterator<InternalImage> i2SquareBarIt(i2SquareBar,i2SquareBar->GetLargestPossibleRegion());
        i2SquareBarIt.GoToBegin();

        for (;!resultIt.IsAtEnd();++resultIt){

            InternalPrecision i1BarV=i1BarIt.Get();
            InternalPrecision i2BarV=i2BarIt.Get();
            InternalPrecision i1BarTimesi2BarV=i1BarV*i2BarV;

            InternalPrecision i1Timesi2BarV=i1Timesi2BarIt.Get();
            InternalPrecision numeratorV=i1Timesi2BarV-i1BarTimesi2BarV;
#if 1 //def RECURSIVE
            //the recursive filtering is not very exact.
            //this is hardly noticeable when only smoothing single images, but during the computation of variances this leads to errors
            //mostly these errors happen when the variance is very low, eg in smooth regions
            //we correct this errors by thresholding, which looks fine when compared to discrete gaussian filtering.
            InternalPrecision tmp=i1SquareBarIt.Get() - i1BarV*i1BarV;

            if (tmp<0){
                LOGV(10)<<VAR(i1SquareBarIt.Get() - i1BarV*i1BarV)<<endl;
                tmp=0;
            }
            InternalPrecision varianceI1V=sqrt(tmp);
            tmp=i2SquareBarIt.Get() - i2BarV*i2BarV;
            if (tmp<0){
                LOGV(10)<<VAR(i2SquareBarIt.Get() - i2BarV*i2BarV)<<endl;
                tmp=0;
            }
            InternalPrecision varianceI2V=sqrt(tmp);
            InternalPrecision denominatorV=varianceI1V*varianceI2V;
            InternalPrecision r = (abs(denominatorV)!=0.0)?numeratorV/(denominatorV):0.0;

#else
            InternalPrecision varianceI1V=sqrt((i1SquareBarIt.Get() - i1BarV*i1BarV));
            InternalPrecision varianceI2V=sqrt((i2SquareBarIt.Get() - i2BarV*i2BarV));
            InternalPrecision denominatorV=varianceI1V*varianceI2V;
            InternalPrecision r = (abs(denominatorV)>0.0)?numeratorV/(denominatorV):0.0;

#endif

            ++i1BarIt;
            ++i2BarIt;
            ++i1Timesi2BarIt;
            ++i1SquareBarIt;
            ++i2SquareBarIt;

#ifdef RECURSIVE
            if (r< -(InternalPrecision)1.0 ){
                r = 0;
            }else if ( r > (InternalPrecision)1.0){
                r = 0;
            }
#endif
            //r = pow(fabs(r),exp);
            //r = max((1.0/(min(1.0-pow(fabs(r),exp),0.00001))),0.0);
            r = (1.0/
                 (1.0
                  -pow(min(0.9999,fabs(r)),exp)
                  )
                 );

#ifdef SAFE
            r = max(r,std::numeric_limits<InternalPrecision>::epsilon());
#endif
            resultIt.Set(r);

        }

        LOGI(8,ImageUtils<InternalImage>::writeImage("result.mha",result));

        return FilterUtils<InternalImage,OutputImage>::cast(result);

    }

    static inline OutputImagePointer coarseLNCC(InputImagePointer i1,InputImagePointer i2, InputImagePointer coarseImg,double sigma=1.0, double exp = 1.0){
        return coarseLNCC( (ConstInputImagePointer)i1, (ConstInputImagePointer)i2, coarseImg,sigma,exp);
    }
    static inline OutputImagePointer coarseLNCC(ConstInputImagePointer i1,ConstInputImagePointer i2, InputImagePointer coarseImg,double sigma=1.0, double exp = 1.0){
        if (exp == 0.0 ) exp = 1.0;
        if (sigma==0.0) sigma=0.001;
        InternalImagePointer i1Cast=FilterUtils<InputImage,InternalImage>::cast(i1);
        InternalImagePointer i2Cast=FilterUtils<InputImage,InternalImage>::cast(i2);
        typedef typename itk::SmoothingRecursiveGaussianImageFilter< InternalImage, InternalImage > FilterType;
        //typedef itk::DiscreteGaussianImageFilter<InternalImage,InternalImage>  FilterType;
        typename FilterType::Pointer filter=FilterType::New();

        filter->SetSigma(sigma);
        //filter->SetVariance(sigma*sigma);
        //        filter->SetMaximumKernelWidth(sigma);

        //compute local means by concolving with gaussian
        filter->SetInput(i1Cast);
        filter->Update();
        InternalImagePointer i1Bar=ImageUtils<InternalImage>::duplicate(filter->GetOutput());
        LOGI(8,ImageUtils<InternalImage>::writeImage("i1Bar.mha",i1Bar));
        filter->SetInput(i2Cast);
        filter->Update();
        InternalImagePointer i2Bar=ImageUtils<InternalImage>::duplicate(filter->GetOutput());

        //FilterUtils<InternalImage,InternalImage>::lowerThresholding(i2Bar,std::numeric_limits<InputImagePixelType>::min());
        LOGI(8,ImageUtils<InternalImage>::writeImage("i2Bar.mha",i2Bar));
        //HACK!
        //compute squares of original images

        InternalImagePointer i1Square=ImageUtils<InternalImage>::multiplyImageOutOfPlace(i1Cast,i1Cast);
        InternalImagePointer i2Square=ImageUtils<InternalImage>::multiplyImageOutOfPlace(i2Cast,i2Cast);
        LOGI(8,ImageUtils<InternalImage>::writeImage("i2Square.mha",i2Square));
        LOGI(8,ImageUtils<InternalImage>::writeImage("i1Square.mha",i1Square));

        //compute local means of squared images by convolving with gaussian kernel
        filter->SetInput(i1Square);
        filter->Update();
        InternalImagePointer i1SquareBar=ImageUtils<InternalImage>::duplicate(filter->GetOutput());
        LOGI(8,ImageUtils<InternalImage>::writeImage("i1SquareBar.mha",i1SquareBar));

        filter->SetInput(i2Square);
        filter->Update();
        InternalImagePointer i2SquareBar=ImageUtils<InternalImage>::duplicate(filter->GetOutput());
        FilterUtils<InternalImage,InternalImage>::lowerThresholding(i2SquareBar,0);
        LOGI(8,ImageUtils<InternalImage>::writeImage("i2SquareBar.mha",i2SquareBar));

        //compute squares of local means
        InternalImagePointer i1BarSquare=ImageUtils<InternalImage>::localSquare(i1Bar);
        InternalImagePointer i2BarSquare=ImageUtils<InternalImage>::multiplyImageOutOfPlace(i2Bar,i2Bar);
        //multiply i1 and i2 locally
        InternalImagePointer i1i2=ImageUtils<InternalImage>::multiplyImageOutOfPlace(i1Cast,i2Cast);
        LOGI(8,ImageUtils<InternalImage>::writeImage("i1i2.mha",i1i2));

        LOGI(8,ImageUtils<InternalImage>::writeImage("i1BarSquare.mha",i1BarSquare));
        LOGI(8,ImageUtils<InternalImage>::writeImage("i2BarSquare.mha",i2BarSquare));

        //compute local means by convolving...
        filter->SetInput(i1i2);
        filter->Update();
        InternalImagePointer i1Timesi2Bar=ImageUtils<InternalImage>::duplicate(filter->GetOutput());
        LOGI(8,ImageUtils<InternalImage>::writeImage("i1Timesi2Bar.mha",i1Timesi2Bar));

        //finish
        InternalImagePointer result=ImageUtils<InternalImage>::duplicate(filter->GetOutput());
        itk::ImageRegionIterator<InternalImage> resultIt(result,result->GetLargestPossibleRegion());
        resultIt.GoToBegin();

        itk::ImageRegionIterator<InternalImage> i1BarIt(i1Bar,i1Bar->GetLargestPossibleRegion());
        i1BarIt.GoToBegin();
        itk::ImageRegionIterator<InternalImage> i2BarIt(i2Bar,i2Bar->GetLargestPossibleRegion());
        i2BarIt.GoToBegin();

        itk::ImageRegionIterator<InternalImage> i1Timesi2BarIt(i1Timesi2Bar,i1Timesi2Bar->GetLargestPossibleRegion());
        i1Timesi2BarIt.GoToBegin();

        itk::ImageRegionIterator<InternalImage> i1SquareBarIt(i1SquareBar,i1SquareBar->GetLargestPossibleRegion());
        i1SquareBarIt.GoToBegin();

        itk::ImageRegionIterator<InternalImage> i2SquareBarIt(i2SquareBar,i2SquareBar->GetLargestPossibleRegion());
        i2SquareBarIt.GoToBegin();

        for (;!resultIt.IsAtEnd();++resultIt){

            InternalPrecision i1BarV=i1BarIt.Get();
            InternalPrecision i2BarV=i2BarIt.Get();
            InternalPrecision i1BarTimesi2BarV=i1BarV*i2BarV;

            InternalPrecision i1Timesi2BarV=i1Timesi2BarIt.Get();
            InternalPrecision numeratorV=i1Timesi2BarV-i1BarTimesi2BarV;

            InternalPrecision varianceI1V=sqrt(max(0.0,i1SquareBarIt.Get() - i1BarV*i1BarV));
            InternalPrecision varianceI2V=sqrt(max(0.0,i2SquareBarIt.Get() - i2BarV*i2BarV));
            InternalPrecision denominatorV=varianceI1V*varianceI2V;

            InternalPrecision r = numeratorV/(denominatorV+0.00000000001);

            ++i1BarIt;
            ++i2BarIt;
            ++i1Timesi2BarIt;
            ++i1SquareBarIt;
            ++i2SquareBarIt;

            if (r< -(InternalPrecision)1.0 ){
                r = 0;
            }else if ( r > (InternalPrecision)1.0){
                r = 0;
            }
            r = pow((r+1.0)/2,exp);

            r = max(r,std::numeric_limits<InternalPrecision>::epsilon());

            resultIt.Set(r);

        }

        LOGI(8,ImageUtils<InternalImage>::writeImage("result.mha",result));

        return FilterUtils<InternalImage,OutputImage>::cast(result);

    }

    static inline OutputImagePointer ITKLMI(InputImagePointer i1,InputImagePointer i2,double sigma=1.0, double exp = 1.0, InputImagePointer coarseImg=NULL){
        return ITKLMI( (ConstInputImagePointer)i1, (ConstInputImagePointer)i2, sigma,exp,coarseImg);
    }
    static inline OutputImagePointer ITKLMI(ConstInputImagePointer i1,ConstInputImagePointer i2,double sigma=1.0, double exp = 1.0, InputImagePointer coarseImg=NULL){

        InternalImagePointer result;
        if (coarseImg.IsNotNull()){
            result=FilterUtils<InputImage,InternalImage>::createEmpty(coarseImg);
            //LOGV(6)<<VAR(result->GetSpacing()[0])<<" "<<VAR(i1->GetSpacing()[0])<<endl;
            //sigma=sigma*result->GetSpacing()[0]/i1->GetSpacing()[0];
        }
        else{
            result=FilterUtils<InputImage,InternalImage>::createEmpty(i1);
        }
        result->FillBuffer(0.0);
        typedef typename InputImage::RegionType RegionType;
        RegionType region;
        //LOGV(4)<<VAR(sigma)<<endl;
        typename InputImage::SizeType regionSize;
        regionSize.Fill(2*(2*sigma+1)); //needs some spacing information!?
        region.SetSize(regionSize);
        itk::ImageRegionIteratorWithIndex<InternalImage> resultIt(result,result->GetLargestPossibleRegion());

        typedef typename itk::MattesMutualInformationImageToImageMetric<InputImage,InputImage> MIType;
        typename MIType::Pointer nccMetric = MIType::New();
        typedef typename itk::LinearInterpolateImageFunction<InputImage> InterpolatorType;
        typename InterpolatorType::Pointer interpolator=InterpolatorType::New();
        interpolator->SetInputImage(i2);
        nccMetric->SetFixedImage(i1);
        nccMetric->SetMovingImage(i2);
        nccMetric->SetInterpolator(interpolator);
        typename itk::IdentityTransform< double, D >::Pointer iTrans= itk::IdentityTransform< double, D >::New();
        nccMetric->SetTransform(iTrans);
        //nccMetric->SetSubtractMean(true);

  resultIt.GoToBegin();
        for (;!resultIt.IsAtEnd();++resultIt){

            typename InputImage::IndexType idx = resultIt.GetIndex(), newIndex;
            bool fullInside=true;
            typename InputImage::PointType pt;
            //get coordinate in fine image in case we're simultaneously downsampling
            if (coarseImg.IsNotNull()){
                result->TransformIndexToPhysicalPoint(idx,pt);
                i1->TransformPhysicalPointToIndex(pt,newIndex);
            }else{
                newIndex=idx;
            }
            LOGV(9)<<VAR(idx)<<" "<<VAR(newIndex)<<" "<<VAR(sigma)<<endl;

#if 0
            //convert newIndex to corner of patch, assumind newIndex is the central pixel
            for (int d=0;d<D;++d){
                if (newIndex[d]>=i1->GetLargestPossibleRegion().GetSize()[d]-sigma){
                    fullInside=false;
                    break;
                }
                newIndex[d]-=sigma;
                if (newIndex[d]<0){
                    fullInside=false;
                    break;
                }
            }
#else

            typename InputImage::SizeType localRegionSize;
             //convert newIndex to corner of patch, assumind newIndex is the central pixel
            for (int d=0;d<D;++d){
                //make region smaller
                int maxSize=2*sigma+1;

                //size gets smaller when center pixel is outside maxRange-sigma
                int idxDifference=newIndex[d] - (i1->GetLargestPossibleRegion().GetSize()[d]-sigma -1);
                maxSize-=max(0, idxDifference);

                //size gets smaller when newIndex-sigma is smaller than zero
                idxDifference=sigma-newIndex[d];
                maxSize-=max(0, idxDifference);

                localRegionSize[d]=maxSize;
                newIndex[d]=max(0.0,newIndex[d]-sigma);

            }
            region.SetSize(localRegionSize);
#endif

            if (fullInside){
                region.SetIndex(newIndex);
                nccMetric->SetFixedImageRegion(region);
                InternalPrecision val=(1.0-nccMetric->GetValue(iTrans->GetParameters()))/2;
                LOGV(9)<<VAR(newIndex)<<" "<<VAR(pt)<<" " <<VAR(val)<<" "<<endl;
                resultIt.Set(pow(val,exp));
            }

        }
        return FilterUtils<InternalImage,OutputImage>::cast(result);

    }

    static inline OutputImagePointer ITKLNCC2(InputImagePointer i1,InputImagePointer i2,double sigma=1.0, double exp = 1.0, InputImagePointer coarseImg=NULL){
        return ITKLNCC2( (ConstInputImagePointer)i1, (ConstInputImagePointer)i2, sigma,exp,coarseImg);
    }
    static inline OutputImagePointer ITKLNCC(ConstInputImagePointer i1,ConstInputImagePointer i2,double sigma=1.0, double exp = 1.0, InputImagePointer coarseImg=NULL){

        InternalImagePointer result;
        if (coarseImg.IsNotNull()){
            result=FilterUtils<InputImage,InternalImage>::createEmpty(coarseImg);
            //LOGV(6)<<VAR(result->GetSpacing()[0])<<" "<<VAR(i1->GetSpacing()[0])<<endl;
            //sigma=sigma*result->GetSpacing()[0]/i1->GetSpacing()[0];
        }
        else{
            result=FilterUtils<InputImage,InternalImage>::createEmpty(i1);
        }
        result->FillBuffer(0.0);
        typedef typename InputImage::RegionType RegionType;
        RegionType region;
        //LOGV(4)<<VAR(sigma)<<endl;
        typename InputImage::SizeType regionSize;
        regionSize.Fill(2*(2*sigma+1)); //needs some spacing information!?
        region.SetSize(regionSize);
        itk::ImageRegionIteratorWithIndex<InternalImage> resultIt(result,result->GetLargestPossibleRegion());

        typedef typename itk::NormalizedCorrelationImageToImageMetric<InputImage,InputImage> NCCType;
        typename NCCType::Pointer nccMetric = NCCType::New();
        typedef typename itk::LinearInterpolateImageFunction<InputImage> InterpolatorType;
        typename InterpolatorType::Pointer interpolator=InterpolatorType::New();
        interpolator->SetInputImage(i2);
        nccMetric->SetFixedImage(i1);
        nccMetric->SetMovingImage(i2);
        nccMetric->SetInterpolator(interpolator);
        typename itk::IdentityTransform< double, D >::Pointer iTrans= itk::IdentityTransform< double, D >::New();
        nccMetric->SetTransform(iTrans);
        nccMetric->SetSubtractMean(true);

  resultIt.GoToBegin();
        for (;!resultIt.IsAtEnd();++resultIt){

            typename InputImage::IndexType idx = resultIt.GetIndex(), newIndex;
            bool fullInside=true;
            typename InputImage::PointType pt;
            //get coordinate in fine image in case we're simultaneously downsampling
            if (coarseImg.IsNotNull()){
                result->TransformIndexToPhysicalPoint(idx,pt);
                i1->TransformPhysicalPointToIndex(pt,newIndex);
            }else{
                newIndex=idx;
            }
            LOGV(9)<<VAR(idx)<<" "<<VAR(newIndex)<<" "<<VAR(sigma)<<endl;

#if 0
            //convert newIndex to corner of patch, assumind newIndex is the central pixel
            for (int d=0;d<D;++d){
                if (newIndex[d]>=i1->GetLargestPossibleRegion().GetSize()[d]-sigma){
                    fullInside=false;
                    break;
                }
                newIndex[d]-=sigma;
                if (newIndex[d]<0){
                    fullInside=false;
                    break;
                }
            }
#else

            typename InputImage::SizeType localRegionSize;
             //convert newIndex to corner of patch, assumind newIndex is the central pixel
            for (int d=0;d<D;++d){
                //make region smaller
                int maxSize=2*sigma+1;

                //size gets smaller when center pixel is outside maxRange-sigma
                int idxDifference=newIndex[d] - (i1->GetLargestPossibleRegion().GetSize()[d]-sigma -1);
                maxSize-=max(0, idxDifference);

                //size gets smaller when newIndex-sigma is smaller than zero
                idxDifference=sigma-newIndex[d];
                maxSize-=std::max(0, idxDifference);

                localRegionSize[d]=maxSize;
                newIndex[d]=std::max(0.0,newIndex[d]-sigma);

            }
            region.SetSize(localRegionSize);
#endif

            if (fullInside){
                region.SetIndex(newIndex);
                nccMetric->SetFixedImageRegion(region);
                InternalPrecision val=(1.0-nccMetric->GetValue(iTrans->GetParameters()))/2;
                LOGV(9)<<VAR(newIndex)<<" "<<VAR(pt)<<" " <<VAR(val)<<" "<<endl;
                resultIt.Set(pow(val,exp));
            }

        }
        return FilterUtils<InternalImage,OutputImage>::cast(result);

    }

    static OutputImagePointer localMetricAutocorrelation(InputImagePointer img1, InputImagePointer img2,double sigma,int nSamples,string metric, double expo=1.0){

        InternalImagePointer lncc=efficientLNCC(img1,img2,sigma);
        InternalImagePointer centerMetric;
        string acc="max";
        string eval="diff";
        if (metric == "lncc"){
            centerMetric = lncc;//efficientLNCC(img1,img2,sigma);
        }else if (metric == "itklncc"){
            centerMetric = ITKLNCC(img1,img2,sigma);
        }else if (metric == "lsad"){
            centerMetric = LSADNorm(img1,img2,sigma,0.0);
        }else if (metric == "lssd"){
            centerMetric = LSSDNorm(img1,img2,sigma,0.0);
        }else{
            std::cout<<"Unknown metric "<<metric<<std::endl;
            exit(0);
        }

        InternalImagePointer accumulator=ImageUtils<InternalImage>::createEmpty(centerMetric);
        accumulator->FillBuffer(0.0);

        int nTotalSamples=1;
        std::vector<InternalImagePointer> metrics(nTotalSamples,NULL);

        DeformationFieldPointerType def=TransfUtils<InputImage>::createEmpty(img2);

        int count=0;
        for (int d=0;d<D;++d){

            for (int sign=-1;sign<2;sign+=2){
                for (int sample=0;sample<nSamples;++sample){

                    InternalPrecision defMag=sign*pow(2,sample)*sigma;
                    typename TransfUtils<InputImage>::DisplacementType disp;
                    disp.Fill(0);
                    disp[d]=defMag;
                    def->FillBuffer(disp);
                    InputImagePointer warped=TransfUtils<InputImage>::warpImage(img1,def);
                    if (metric == "lncc"){
                        metrics[count] = efficientLNCC(warped,img2,sigma);
                    }else if (metric == "itklncc"){
                        metrics[count] = ITKLNCC(warped,img2,sigma);
                    }else if (metric == "lsad"){
                        metrics[count] = LSADNorm(warped,img2,sigma,0.0);
                    }else if (metric == "lssd"){
                        metrics[count] = LSSDNorm(warped,img2,sigma,0.0);
                    }else{
                        std::cout<<"Unknown metric "<<metric<<std::endl;
                        exit(0);
                    }

                    if (acc=="mean"){
                        accumulator=FilterUtils<InternalImage>::add(accumulator,metrics[count]);
                    }else if (acc=="max"){
                        FilterUtils<InternalImage>::localMax(accumulator,metrics[count]);
                    }else{
                        LOG<<"unknown accumulation method "<<acc<<endl;
                        exit(0);
                    }

                    //++count;

                }
            }
        }
        if (acc == "mean")
            ImageUtils<InternalImage>::multiplyImage(accumulator,1.0/(D*2*nSamples));

        if (eval=="diff"){
            // l= 0.5 + (\hat l - mu)/(2*\hat l)
            //InternalImagePointer diff=FilterUtils<InternalImage>::substract(centerMetric,accumulator);
            //centerMetric=FilterUtils<InternalImage>::lowerThresholding(centerMetric,1e-10);
            //centerMetric=ImageUtils<InternalImage>::divideImageOutOfPlace(diff,centerMetric);
            //ImageUtils<InternalImage>::multiplyImage(centerMetric,0.5);

            //ImageUtils<InternalImage>::add(centerMetric,0.5);
            centerMetric=ImageUtils<InternalImage>::divideImageOutOfPlace(centerMetric,accumulator);
            centerMetric=ImageUtils<InternalImage>::multiplyImageOutOfPlace(lncc,centerMetric);
            //centerMetric=FilterUtils<InternalImage>::thresholding(centerMetric,0,1);
            typename ImageUtils<InternalImage>::ImageIteratorType resultIt(centerMetric,centerMetric->GetLargestPossibleRegion());
#if 1
            if (expo !=1.0){
                for (resultIt.GoToBegin();!resultIt.IsAtEnd();++resultIt){
                    resultIt.Set(pow(resultIt.Get(),expo));
                }
            }
#endif

        }else if (eval == "exp" ){

            typename ImageUtils<InternalImage>::ImageIteratorType accIt(accumulator,accumulator->GetLargestPossibleRegion());
            typename ImageUtils<InternalImage>::ImageIteratorType resultIt(centerMetric,centerMetric->GetLargestPossibleRegion());

            accIt.GoToBegin();
            for (resultIt.GoToBegin();!resultIt.IsAtEnd();++resultIt,++accIt){
                InternalPrecision Lprime=resultIt.Get();
                InternalPrecision Lbar=accIt.Get();

                InternalPrecision weight;
                if (Lprime != 0.0)
                    weight = exp (- 0.6931 * Lbar/Lprime);
                else
                    weight = 0.0;
                resultIt.Set(weight);

            }

        }else{
            LOG<<"Unknown evaluation method "<<eval<<endl;
            exit(0);
        }
        return centerMetric;
    }

    static double nCC(InputImagePointer img1, InputImagePointer img2, DeformationFieldPointerType def){

        InputImagePointer warpedImage2=TransfUtils<InputImage>::warpImage(img2,def);
        typedef typename itk::NormalizedCorrelationImageToImageMetric<InputImage,InputImage> NCCType;
        typename NCCType::Pointer nccMetric = NCCType::New();
        typedef typename itk::LinearInterpolateImageFunction<InputImage> InterpolatorType;
        typename InterpolatorType::Pointer interpolator=InterpolatorType::New();
        interpolator->SetInputImage(warpedImage2);
        nccMetric->SetFixedImage(img1);
        nccMetric->SetMovingImage(warpedImage2);
        nccMetric->SetInterpolator(interpolator);
        typename itk::IdentityTransform< double, D >::Pointer iTrans= itk::IdentityTransform< double, D >::New();
        nccMetric->SetTransform(iTrans);
        nccMetric->SetSubtractMean(true);
        nccMetric->SetFixedImageRegion(img1->GetLargestPossibleRegion());
        return nccMetric->GetValue(iTrans->GetParameters());

    }

  static double nCC(InputImagePointer img1, InputImagePointer warpedImage2){

      typedef typename itk::NormalizedCorrelationImageToImageMetric<InputImage,InputImage> NCCType;
        typename NCCType::Pointer nccMetric = NCCType::New();
        typedef typename itk::LinearInterpolateImageFunction<InputImage> InterpolatorType;
        typename InterpolatorType::Pointer interpolator=InterpolatorType::New();
        interpolator->SetInputImage(warpedImage2);
        nccMetric->SetFixedImage(img1);
        nccMetric->SetMovingImage(warpedImage2);
        nccMetric->SetInterpolator(interpolator);
        typename itk::IdentityTransform< double, D >::Pointer iTrans= itk::IdentityTransform< double, D >::New();
        nccMetric->SetTransform(iTrans);
        nccMetric->SetSubtractMean(true);
        nccMetric->SetFixedImageRegion(img1->GetLargestPossibleRegion());
        return nccMetric->GetValue(iTrans->GetParameters());

    }

      static double msd(InputImagePointer img1, InputImagePointer warpedImage2){

           itk::ImageRegionIterator<InputImage> it1(img1,img1->GetRequestedRegion());
      itk::ImageRegionIterator<InputImage> it2(warpedImage2,warpedImage2->GetRequestedRegion());
      double result=0.0;
      int c=0;
      it1.GoToBegin();
      it2.GoToBegin();
      for (;!it1.IsAtEnd();++it1,++it2,++c){
          double val=(it1.Get()-it2.Get());
          val*=val;
          result+=val;
      }
      return result/c;

    }

  static double mad(InputImagePointer img1, InputImagePointer warpedImage2){
      itk::ImageRegionIterator<InputImage> it1(img1,img1->GetRequestedRegion());
      itk::ImageRegionIterator<InputImage> it2(warpedImage2,warpedImage2->GetRequestedRegion());
      double result=0.0;
      int c=0;
      it1.GoToBegin();
      it2.GoToBegin();
      for (;!it1.IsAtEnd();++it1,++it2,++c){
          result+=fabs(it1.Get()-it2.Get());
      }
      return result/c;

    }

       static inline InternalImagePointer CategoricalDiff(InputImagePointer i1,InputImagePointer i2,double sigmaWidth=1.0){
        if (sigmaWidth==0.0) sigmaWidth=0.001;
        return CategoricalDiff( (ConstInputImagePointer)i1, (ConstInputImagePointer)i2, sigmaWidth);
    }
    static inline InternalImagePointer CategoricalDiff(ConstInputImagePointer i1,ConstInputImagePointer i2,double sigmaWidth=1.0){
        if (sigmaWidth==0.0) sigmaWidth=0.001;
        //typedef typename itk::SmoothingRecursiveGaussianImageFilter< InternalImage, InternalImage > FilterType;
        typedef itk::DiscreteGaussianImageFilter<InternalImage,InternalImage>  FilterType;
        typename FilterType::Pointer filter=FilterType::New();
        //filter->SetSigma(sigmaWidth);
        filter->SetVariance(sigmaWidth*sigmaWidth);

        InternalImagePointer i1Cast=FilterUtils<InputImage,InternalImage>::cast(i1);
        InternalImagePointer i2Cast=FilterUtils<InputImage,InternalImage>::cast(i2);

        InternalImagePointer diff =  ImageUtils<InternalImage>::createEmpty(i1Cast);
        typename ImageUtils<InternalImage>::ImageIteratorType diffIt(diff,i1Cast->GetLargestPossibleRegion());
        typename ImageUtils<InternalImage>::ImageIteratorType i1It(i1Cast,i1Cast->GetLargestPossibleRegion());
        typename ImageUtils<InternalImage>::ImageIteratorType i2It(i2Cast,i1Cast->GetLargestPossibleRegion());

        InternalPrecision mean=0.0;
        int c=0;
        //subtract and take absolute value
        for (i2It.GoToBegin(),i1It.GoToBegin(), diffIt.GoToBegin(); !diffIt.IsAtEnd(); ++i1It, ++i2It, ++diffIt,++c){
            InternalPrecision d=fabs(i1It.Get()-i2It.Get());
            diffIt.Set(d<std::numeric_limits<InternalPrecision>::epsilon());
            mean+=d;
        }
        LOGV(6)<<VAR(mean/c)<<endl;

        //compute local means by concolving with gaussian
        filter->SetInput(diff);
        filter->Update();
        InternalImagePointer result=filter->GetOutput();

        return FilterUtils<InternalImage,OutputImage>::cast(result);
  }
    static inline InternalImagePointer CategoricalDiffNorm(InputImagePointer i1,InputImagePointer i2,double sigmaWidth=1.0, double sigmaNorm=1.0){
        if (sigmaWidth==0.0) sigmaWidth=0.001;
        return CategoricalDiffNorm( (ConstInputImagePointer)i1, (ConstInputImagePointer)i2, sigmaWidth,sigmaNorm);
    }
    static inline InternalImagePointer CategoricalDiffNorm(ConstInputImagePointer i1,ConstInputImagePointer i2,double sigmaWidth=1.0,double sigmaNorm=1.0){
        if (sigmaWidth==0.0) sigmaWidth=0.001;
        //typedef typename itk::SmoothingRecursiveGaussianImageFilter< InternalImage, InternalImage > FilterType;
        typedef itk::DiscreteGaussianImageFilter<InternalImage,InternalImage>  FilterType;
        typename FilterType::Pointer filter=FilterType::New();
        //filter->SetSigma(sigmaWidth);
        filter->SetVariance(sigmaWidth*sigmaWidth);

        InternalImagePointer i1Cast=FilterUtils<InputImage,InternalImage>::cast(i1);
        InternalImagePointer i2Cast=FilterUtils<InputImage,InternalImage>::cast(i2);

        InternalImagePointer diff =  ImageUtils<InternalImage>::createEmpty(i1Cast);
        typename ImageUtils<InternalImage>::ImageIteratorType diffIt(diff,i1Cast->GetLargestPossibleRegion());
        typename ImageUtils<InternalImage>::ImageIteratorType i1It(i1Cast,i1Cast->GetLargestPossibleRegion());
        typename ImageUtils<InternalImage>::ImageIteratorType i2It(i2Cast,i1Cast->GetLargestPossibleRegion());

        InternalPrecision mean=0.0;
        int c=0;
        //subtract and take absolute value
        for (i2It.GoToBegin(),i1It.GoToBegin(), diffIt.GoToBegin(); !diffIt.IsAtEnd(); ++i1It, ++i2It, ++diffIt,++c){
            InternalPrecision d=fabs(i1It.Get()-i2It.Get());
            diffIt.Set(d<std::numeric_limits<InternalPrecision>::epsilon());
            mean+=d;
        }
        LOGV(6)<<VAR(mean/c)<<endl;

        //compute local means by concolving with gaussian
        filter->SetInput(diff);
        filter->Update();
        InternalImagePointer result=filter->GetOutput();
        typename ImageUtils<InternalImage>::ImageIteratorType resultIt(result,result->GetLargestPossibleRegion());

        for (resultIt.GoToBegin(); !resultIt.IsAtEnd(); ++resultIt){
            InternalPrecision d = resultIt.Get();
            //resultIt.Set(max(0.5,exp(-0.5 * fabs(d)  / sigmaNorm) ));
            resultIt.Set(exp(-0.5 * fabs(d)  / sigmaNorm) );
            LOGV(7)<<VAR(d)<<" "<<resultIt.Get()<<endl;
            //resultIt.Set(pow(exp(-0.5 * fabs(d)  / mean ),sigmaNorm ));
        }

        return FilterUtils<InternalImage,OutputImage>::cast(result);
    }
#ifdef WITH_MIND
    static inline OutputImagePointer deedsMIND(InputImagePointer i1,InputImagePointer i2,double sigmaWidth=1.0,double sigmaNorm=1.0){
        InternalImagePointer i1Cast=FilterUtils<InputImage,InternalImage>::cast(i1);
        InternalImagePointer i2Cast=FilterUtils<InputImage,InternalImage>::cast(i2);
        float * i1Data=i1Cast->GetBufferPointer();
        float * i2Data=i2Cast->GetBufferPointer();
        InternalImagePointer result=ImageUtils<InternalImage>::createEmpty(i1Cast);
        float * resultData=result->GetBufferPointer();
        typename InputImage::SizeType size=i1->GetLargestPossibleRegion().GetSize();
        int hw = 0;
        int sparse=1;
        int r=int(sigmaWidth);
        dataRegSSC(resultData,i1Data,i2Data,hw,sparse,r,0,size[0],size[1],D>2?size[2]:1);
        typename ImageUtils<InternalImage>::ImageIteratorType resultIt(result,result->GetLargestPossibleRegion());
        for (resultIt.GoToBegin(); !resultIt.IsAtEnd(); ++resultIt){
            InternalPrecision d = resultIt.Get();
            resultIt.Set(InternalPrecision(1.0)-pow(min(d,InternalPrecision(1.0)),InternalPrecision(sigmaNorm)));
            //resultIt.Set(-d);
            //resultIt.Set(pow(exp(-0.5 * fabs(d)  / mean ),sigmaNorm ));
        }
        return FilterUtils<InternalImage,OutputImage>::cast(result);
    }
    static inline OutputImagePointer deedsLCC(InputImagePointer i1,InputImagePointer i2,double sigmaWidth=1.0,double sigmaNorm=1.0){
        InternalImagePointer i1Cast=FilterUtils<InputImage,InternalImage>::cast(i1);
        InternalImagePointer i2Cast=FilterUtils<InputImage,InternalImage>::cast(i2);
        float * i1Data=i1Cast->GetBufferPointer();
        float * i2Data=i2Cast->GetBufferPointer();
        InternalImagePointer result=ImageUtils<InternalImage>::createEmpty(i1Cast);
        float * resultData=result->GetBufferPointer();
        typename InputImage::SizeType size=i1->GetLargestPossibleRegion().GetSize();
        int hw = 0;
        int sparse=1;
        int r=int(sigmaWidth);
        dataRegLCC(resultData,i1Data,i2Data,hw,sparse,r,0,size[0],size[1],D>2?size[2]:1);
        typename ImageUtils<InternalImage>::ImageIteratorType resultIt(result,result->GetLargestPossibleRegion());
        for (resultIt.GoToBegin(); !resultIt.IsAtEnd(); ++resultIt){
            InternalPrecision d = resultIt.Get();
            resultIt.Set(pow(1.0-d,InternalPrecision(sigmaNorm)));
            //resultIt.Set(pow(exp(-0.5 * fabs(d)  / mean ),sigmaNorm ));
        }
        return FilterUtils<InternalImage,OutputImage>::cast(result);
    }
#endif
    static inline OutputImagePointer multiScaleLNCCAbs(InputImagePointer i1,InputImagePointer i2,double sigmaMax=1.0, double exp = 1.0){
        double sigma=1.0;
        int count=0;
        InternalImagePointer result;
        for (;sigma<=sigmaMax;sigma*=2){
            if (!count)
                result=efficientLNCCNewNorm(i1,i2,sigma,1.0);
            else
                result=FilterUtils<InternalImage>::add(result,efficientLNCCNewNorm(i1,i2,sigma,1.0));
            ++count;

        }
        typename ImageUtils<InternalImage>::ImageIteratorType resultIt(result,result->GetLargestPossibleRegion());
        for (resultIt.GoToBegin(); !resultIt.IsAtEnd(); ++resultIt){
            InternalPrecision d = resultIt.Get();
            resultIt.Set(pow(d/count,InternalPrecision(exp)));
        }
        return FilterUtils<InternalImage,OutputImage>::cast(result);

    }
};


text:
import * as npmVersioning from '../../versioning/npm';

export { resetMemCache, resetCache } from './get';
export { getReleases } from './releases';
export { getNpmrc, setNpmrc } from './npmrc';
export { id } from './common';
export const defaultVersioning = npmVersioning.id;
export const customRegistrySupport = false;

redpajama_set_name: RedPajamaGithub


text:
import { PrimaryGeneratedColumn, Column, Entity, OneToOne, OneToMany, ManyToOne, JoinColumn } from "typeorm";
import { Post } from "./Post";

@Entity("PostDetails")
export class PostDetails {

    @PrimaryGeneratedColumn()
    id: number;

    @Column({
        nullable: true
    })
    authorName: string;

    @Column({
        nullable: true
    })
    comment: string;

    @Column({
        nullable: true
    })
    metadata: string;

    @OneToMany(type => Post, post => post.details)
    posts: Post[];

}


text:
package gobblin.runtime.std;

import java.net.URI;
import java.net.URISyntaxException;

import com.google.common.base.Optional;
import com.google.common.base.Preconditions;
import com.google.common.base.Predicate;
import com.google.common.base.Predicates;

import gobblin.runtime.api.JobSpec;

import lombok.AllArgsConstructor;
import lombok.Getter;
import lombok.Setter;

/**
 * A helper class to create JobSpec predicates.
 */
@AllArgsConstructor
public class JobSpecFilter implements Predicate<JobSpec> {
  private final Optional<Predicate<URI>> uriPredicate;
  private final Optional<Predicate<String>> versionPredicate;

  @Override
  public boolean apply(JobSpec input) {
    Preconditions.checkNotNull(input);

    boolean res = true;
    if (this.uriPredicate.isPresent()) {
      res &= this.uriPredicate.get().apply(input.getUri());
    }
    if (res && this.versionPredicate.isPresent()) {
      res &= this.versionPredicate.get().apply(input.getVersion());
    }
    return res;
  }

  public static Builder builder() {
    return new Builder();
  }

  public static JobSpecFilter eqJobSpecURI(URI jobSpecURI) {
    return builder().eqURI(jobSpecURI).build();
  }

  public static JobSpecFilter eqJobSpecURI(String jobSpecURI) {
    return builder().eqURI(jobSpecURI).build();
  }

  @Getter
  @Setter
  public static class Builder {
    private Predicate<URI> uriPredicate;
    private Predicate<String> versionPredicate;

    public Builder eqURI(URI uri) {
      this.uriPredicate = Predicates.equalTo(uri);
      return this;
    }

    public Builder eqURI(String uri) {
      try {
        this.uriPredicate = Predicates.equalTo(new URI(uri));
      } catch (URISyntaxException e) {
        throw new RuntimeException("invalid URI: " + uri, e);
      }
      return this;
    }

    public Builder eqVersion(String version) {
      this.versionPredicate = Predicates.equalTo(version);
      return this;
    }

    public JobSpecFilter build() {
      return new JobSpecFilter(Optional.fromNullable(this.uriPredicate),
             Optional.fromNullable(this.versionPredicate));
    }
  }

}


text:
class N_PUBLIC TMouseCursor : public TObject
{
	TString	Image;
    TTexture* Texture;
	TPoint Hotspot;
	int Width;
	int	Height;

public:
    TMouseCursor();
    virtual ~TMouseCursor();
	virtual void Create(const char* Filename, const TPoint& AHotspot);
	virtual void Create(TString Filename, const TPoint& AHotspot);

	TTexture* GetTexture(void) const { return Texture; }

	void SetImage(TString Filename);
	TString GetImage(void) const { return Image; }

	void SetHotspot(TPoint Value) { Hotspot = Value; }
    TPoint GetHotspot(void) const { return Hotspot; }

	void SetWidth(int Value) { Width = Value; }
	int GetWidth(void) const { return Width; }

	void SetHeight(int Value) { Height = Value; }
	int GetHeight(void) const { return Height; }
};

//-------------------------------------------------------------------
#endif


text:
package com.mxgraph.util;

import java.util.ArrayList;
import java.util.List;

/**
 * Implements a 2-dimensional rectangle with double precision coordinates.
 */
public class mxUndoableEdit
{

	/**
	 * Defines the requirements for an undoable change.
	 */
	public interface mxUndoableChange
	{

		/**
		 * Undoes or redoes the change depending on its undo state.
		 */
		void execute();

	}

	/**
	 * Holds the source of the undoable edit.
	 */
	protected Object source;

	/**
	 * Holds the list of changes that make up this undoable edit.
	 */
	protected List<mxUndoableChange> changes = new ArrayList<mxUndoableChange>();

	/**
	 * Specifies this undoable edit is significant. Default is true.
	 */
	protected boolean significant = true;

	/**
	 * Specifies the state of the undoable edit.
	 */
	protected boolean undone, redone;

	/**
	 * Constructs a new undoable edit for the given source.
	 */
	public mxUndoableEdit(Object source)
	{
		this(source, true);
	}

	/**
	 * Constructs a new undoable edit for the given source.
	 */
	public mxUndoableEdit(Object source, boolean significant)
	{
		this.source = source;
		this.significant = significant;
	}

	/**
	 * Hook to notify any listeners of the changes after an undo or redo
	 * has been carried out. This implementation is empty.
	 */
	public void dispatch()
	{
		// empty
	}

	/**
	 * Hook to free resources after the edit has been removed from the command
	 * history. This implementation is empty.
	 */
	public void die()
	{
		// empty
	}

	/**
	 * @return the source
	 */
	public Object getSource()
	{
		return source;
	}

	/**
	 * @return the changes
	 */
	public List<mxUndoableChange> getChanges()
	{
		return changes;
	}

	/**
	 * @return the significant
	 */
	public boolean isSignificant()
	{
		return significant;
	}

	/**
	 * @return the undone
	 */
	public boolean isUndone()
	{
		return undone;
	}

	/**
	 * @return the redone
	 */
	public boolean isRedone()
	{
		return redone;
	}

	/**
	 * Returns true if the this edit contains no changes.
	 */
	public boolean isEmpty()
	{
		return changes.isEmpty();
	}

	/**
	 * Adds the specified change to this edit. The change is an object that is
	 * expected to either have an undo and redo, or an execute function.
	 */
	public void add(mxUndoableChange change)
	{
		changes.add(change);
	}

	/**
	 *
	 */
	public void undo()
	{
		if (!undone)
		{
			int count = changes.size();

			for (int i = count - 1; i >= 0; i--)
			{
				mxUndoableChange change = changes.get(i);
				change.execute();
			}

			undone = true;
			redone = false;
		}

		dispatch();
	}

	/**
	 *
	 */
	public void redo()
	{
		if (!redone)
		{
			int count = changes.size();

			for (int i = 0; i < count; i++)
			{
				mxUndoableChange change = changes.get(i);
				change.execute();
			}

			undone = false;
			redone = true;
		}

		dispatch();
	}

}


text:
use crossterm::style::Color;
use pueue_lib::settings::Settings;

/// Color wrapper for actual colors depending on settings
/// Using dark colors if dark_mode is enabled
pub struct Colors {
    /// red color
    red: Color,
    /// green color
    green: Color,
    /// white color
    white: Color,
    /// yellow color
    yellow: Color,
}

impl Colors {
    /// init color-scheme depending on settings
    pub const fn new(settings: &Settings) -> Self {
        if settings.client.dark_mode {
            Self {
                green: Color::DarkGreen,
                red: Color::DarkRed,
                yellow: Color::DarkYellow,
                white: Color::White,
            }
        } else {
            Self {
                green: Color::Green,
                red: Color::Red,
                yellow: Color::Yellow,
                white: Color::White,
            }
        }
    }

    /// return green color
    pub const fn green(&self) -> Color {
        self.green
    }
    /// return red color
    pub const fn red(&self) -> Color {
        self.red
    }
    /// return yellow color
    pub const fn yellow(&self) -> Color {
        self.yellow
    }
    /// return white color
    pub const fn white(&self) -> Color {
        self.white
    }
}


text:
using System.Threading;
using Microsoft.CodeAnalysis.CSharp.Extensions.ContextQuery;

namespace Microsoft.CodeAnalysis.CSharp.Completion.KeywordRecommenders
{
    internal class AsKeywordRecommender : AbstractSyntacticSingleKeywordRecommender
    {
        public AsKeywordRecommender()
            : base(SyntaxKind.AsKeyword)
        {
        }

        protected override bool IsValidContext(int position, CSharpSyntaxContext context, CancellationToken cancellationToken)
        {
            return !context.IsInNonUserCode && context.IsIsOrAsContext;
        }
    }
}


text:
<?php
/**
 * Created by IntelliJ IDEA.
 * User: mduncan
 * Date: 8/17/15
 * Time: 3:47 PM
 */

namespace Fulfillment\OMS\Api;
use Fulfillment\OMS\Models\Request\Contracts\User;

/**
 * Class UsersApi
 * @package Fulfillment\OMS\Api
 */
class UsersApi extends ApiRequestBase
{
    /**
     * Get an array of Users
     *
     * Will only returns users from the merchant authorizing the request
     *
     * @param $queryString
     * @return array|User[]
     * @throws \Exception
     */
    public function getUsers($queryString)
    {
        $json = $this->apiClient->get('users', $queryString);

        if ($this->jsonOnly) {
            $users = $json;
        } else {
            $users = $this->jsonMapper->mapArray($json, [], get_class($this->container->get(User::class)));
        }
        return $users;
    }

    /**
     * Get a user with the specified Id
     *
     * @param $userId
     * @return array|User
     * @throws \Exception
     * @throws \JsonMapper_Exception
     */
    public function getUser($userId)
    {
        $json = $this->apiClient->get('users/' . $userId);

        if ($this->jsonOnly) {
            $user = $json;
        } else {
            $user = $this->jsonMapper->map($json, $this->container->get(User::class));
        }

        return $user;
    }

    /**
     * Create a new user in the Fulfillment system
     *
     * @param array|\Fulfillment\OMS\Models\Request\User $user
     * @return array|User
     * @throws \Exception
     * @throws \Fulfillment\OMS\Exceptions\ValidationFailureException
     * @throws \JsonMapper_Exception
     */
    public function createUser($user)
    {
        if ($user instanceof \Fulfillment\OMS\Models\Request\User && $this->validateRequests) {
            $user->validate();
        }
        $json = $this->apiClient->post('users/', $user);

        if ($this->jsonOnly) {
            $returnUser = $json;
        } else {
            $returnUser = $this->jsonMapper->map($json, $this->container->get(User::class));
        }
        return $returnUser;
    }

    /**
     * Update an existing user
     *
     * @param $userId
     * @param array|\Fulfillment\OMS\Models\Request\User $user
     * @throws \Exception
     * @throws \Fulfillment\OMS\Exceptions\ValidationFailureException
     */
    public function updateUser($userId, $user)
    {
        if ($user instanceof \Fulfillment\OMS\Models\Request\User && $this->validateRequests) {
            $user->validate();
        }

        $this->apiClient->put('users/' . $userId, $user);
    }

    /**
     * Delete the user with the specified Id
     *
     * @param $userId
     * @throws \Exception
     */
    public function deleteUser($userId)
    {
        $this->apiClient->delete('users/' . $userId);
    }

    /**
     * Change a User's password
     *
     * Must provide the current password for this to work
     *
     * @param $userId
     * @param $currentPassword
     * @param $newPassword
     * @throws \Exception
     */
    public function changeUserPassword($userId, $currentPassword, $newPassword)
    {
        $this->apiClient->put('users/' . $userId . '/password', ['currentPassword' => $currentPassword, 'newPassword' => $newPassword]);
    }

    /**
     * Generate a password token for the specified user and notifies them via email associated with their account
     *
     * Note: Generating a new token will invalidate a previous token
     *
     * @param $userId
     * @throws \Exception
     */
    public function resetUserPassword($userId)
    {
        $this->apiClient->get('users/' . $userId . '/password/reset');
    }

    /**
     * Add an array of roles to a User
     *
     * @param $userId int
     * @param $rolesArray array
     * @throws \Exception
     */
    public function addRoles($userId, $rolesArray)
    {
        $this->apiClient->put('users/' . $userId, $rolesArray);
    }

    /**
     * Delete an array of roles from a User
     *
     * @param $userId int
     * @param $rolesArray array
     * @throws \Exception
     */
    public function deleteRoles($userId, $rolesArray)
    {
        $this->apiClient->delete('users/' . $userId, $rolesArray);
    }

    /**
     * Generate a password token for the user associated with this email address and notify them.
     *
     * Note: Generating a token will invalidate any previous tokens
     *
     * @param $email
     * @throws \Exception
     */
    public function generatePasswordToken($email)
    {
        $this->apiClient->post('generatePasswordToken', ['email' => $email]);
    }

    /**
     * Reset the password for a user using their password token and a new password
     *
     * @param $token string
     * @param $newPassword string
     * @throws \Exception
     */
    public function resetPassword($token, $newPassword)
    {
        $this->apiClient->post('resetPassword', ['token' => $token, 'newPassword' => $newPassword]);
    }
}


text:
module DbCharmer
  module Sharding
    autoload :Connection, 'db_charmer/sharding/connection'
    autoload :StubConnection, 'db_charmer/sharding/stub_connection'
    autoload :Method, 'db_charmer/sharding/method'

    @@sharded_connections = {}

    def self.register_connection(config)
      name = config[:name] or raise ArgumentError, "No :name in connection!"
      @@sharded_connections[name] = DbCharmer::Sharding::Connection.new(config)
    end

    def self.sharded_connection(name)
      @@sharded_connections[name] or raise ArgumentError, "Invalid sharded connection name!"
    end
  end
end


text:
<?php
namespace app\components;

use app\Events\AddPlayerEvent;
use app\Events\GetPlayersEvent;
use yii\base\Component;
use yii\base\Exception;

class Event extends Component
{
    const EVENT_ADD_PLAYER = 'addPlayer';

    public function addPlayer(array $playerParams)
    {
        $event = new AddPlayerEvent;
        $event->data = $playerParams;
        try {
            return $event->run();
        } catch (Exception $exception) {
            throw new Exception('Could not insert player');
        }
    }

    public function getAllPlayers()
    {
        $event = new GetPlayersEvent();
        return $event->run();
    }
}


text:
package net.sf.jabref.gui.groups;

import net.sf.jabref.logic.l10n.Localization;
import net.sf.jabref.model.groups.ExplicitGroup;
import net.sf.jabref.model.groups.KeywordGroup;
import net.sf.jabref.model.groups.SearchGroup;
import net.sf.jabref.model.strings.StringUtil;

public class GroupDescriptions {

    public static String getDescriptionForPreview(String field, String expr, boolean caseSensitive, boolean regExp) {
        String header = regExp ? Localization.lang("This group contains entries whose %0 field contains the regular expression %1",
                field, StringUtil.quoteForHTML(expr))
                : Localization.lang("This group contains entries whose %0 field contains the keyword %1",
                field, StringUtil.quoteForHTML(expr));
        String caseSensitiveText = caseSensitive ? Localization.lang("case sensitive") :
            Localization.lang("case insensitive");
        String footer = regExp ?
                Localization.lang("Entries cannot be manually assigned to or removed from this group.")
                : Localization.lang(
                "Additionally, entries whose %0 field does not contain "
                        + "%1 can be assigned manually to this group by selecting them "
                        + "then using either drag and drop or the context menu. "
                        + "This process adds the term %1 to "
                        + "each entry's %0 field. "
                        + "Entries can be removed manually from this group by selecting them "
                        + "then using the context menu. "
                        + "This process removes the term %1 from "
                        + "each entry's %0 field.",
                field, StringUtil.quoteForHTML(expr));
        return String.format("%s (%s). %s", header, caseSensitiveText, footer);
    }

    public static String getShortDescriptionKeywordGroup(KeywordGroup keywordGroup, boolean showDynamic) {
        StringBuilder sb = new StringBuilder();
        sb.append("");
        if (showDynamic) {
            sb.append("").append(StringUtil.quoteForHTML(keywordGroup.getName())).append("");
        } else {
            sb.append(StringUtil.quoteForHTML(keywordGroup.getName()));
        }
        sb.append(" - ");
        sb.append(Localization.lang("dynamic group"));
        sb.append(" ");
        sb.append(keywordGroup.getSearchField());
        sb.append(" ");
        sb.append(Localization.lang("contains"));
        sb.append(" ");
        sb.append(StringUtil.quoteForHTML(keywordGroup.getSearchExpression()));
        sb.append(")");
        switch (keywordGroup.getHierarchicalContext()) {
        case INCLUDING:
            sb.append(", ").append(Localization.lang("includes subgroups"));
            break;
        case REFINING:
            sb.append(", ").append(Localization.lang("refines supergroup"));
            break;
        default:
            break;
        }
        return sb.toString();

    }

    public static String getDescriptionForPreview() {
        return Localization.lang("This group contains entries based on manual assignment. "
                + "Entries can be assigned to this group by selecting them "
                + "then using either drag and drop or the context menu. "
                + "Entries can be removed from this group by selecting them "
                + "then using the context menu.");
    }

    public static String getShortDescriptionExplicitGroup(ExplicitGroup explicitGroup) {
        StringBuilder sb = new StringBuilder();
        sb.append("").append(explicitGroup.getName()).append(" - ").append(Localization.lang("static group"));
        switch (explicitGroup.getHierarchicalContext()) {
        case INCLUDING:
            sb.append(", ").append(Localization.lang("includes subgroups"));
            break;
        case REFINING:
            sb.append(", ").append(Localization.lang("refines supergroup"));
            break;
        default:
            break;
        }
        return sb.toString();
    }

    public static String getShortDescriptionAllEntriesGroup() {
        return Localization.lang("All Entries (this group cannot be edited or removed)");
    }

    public static String getShortDescription(SearchGroup searchGroup, boolean showDynamic) {

        StringBuilder sb = new StringBuilder();
        sb.append("");
        if (showDynamic) {
            sb.append("").append(StringUtil.quoteForHTML(searchGroup.getName())).append("");
        } else {
            sb.append(StringUtil.quoteForHTML(searchGroup.getName()));
        }
        sb.append(" - ");
        sb.append(Localization.lang("dynamic group"));
        sb.append(" (");
        sb.append(Localization.lang("search expression"));
        sb.append(" ").append(StringUtil.quoteForHTML(searchGroup.getSearchExpression())).append(")");
        switch (searchGroup.getHierarchicalContext()) {
        case INCLUDING:
            sb.append(", ").append(Localization.lang("includes subgroups"));
            break;
        case REFINING:
            sb.append(", ").append(Localization.lang("refines supergroup"));
            break;
        default:
            break;
        }
        return sb.toString();
    }

}


text:
{% extends '../layouts/template.html' %}
{% block title %}Check my State Pension - {{mvpversion}}{% endblock %}
{% block serviceName %}{% endblock %}
{% block main %}

Prototype 29 persona logins

can improve forecast with gaps

      davida 123456

will get maximum before SPA - gaps not needed

    janeb 123456

already reached the maximum

    johnc 123456

    GOV.UK start page

{% endblock %}


text:
maojian
haoguanwei
renwei

# Author
all

# Reviewer
all


text:
<meta http-equiv="content-type" content="text/html;charset=UTF-8" />
<title>Algebra - I</title>

<meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">

  COMBUSTION AND SHOCKWAVES

Basic background material on ODE, PDE and fluid
mechanics, Conservation laws, Combustion modelling, Blast waves and detonations,
Hamillton-Jacobi equations. Special topics from one of the following three: (1)
Numerics of conservation laws, (2) Mathematicsal problems in low Mach number
combustions, (3) Formation and propagation of singularities on a nonlinear
wavefront in two-dimensions.

Books

    Evans. L. C., Partial Differential Equations,
    AMS, 1998, Providence, RI.

    Godlewski, E., and Raviart P. A., Numerical
    Approximations of Hyperbolic systems of Conservation laws, Springer, 1996,
    New York.

    Prasad, P., Nonlinear Hyperbolic Waves in
    Multi-dimensions, Chapman & Hall/CRC, 2001, New York.

    Smoller, J., Shock Waves and
    Reaction-Diffusion equations, Springer-Verlag, 1983, New York.

    Williams, F. A., Combustion Theory, Benjamin
    Cummings, 1985, Menlo Park.


text:
  Castineta est
       un village
       géographiquement positionné dans le département de Haute-Corse en Corse. Elle totalisait 53 habitants en 2008.

À coté de Castineta sont positionnées géographiquement les villes de
Ortiporio localisée à 5 km, 123 habitants,
Ficaja située à 5 km, 49 habitants,
Morosaglia localisée à 1 km, 1 058 habitants,
La Porta localisée à 4 km, 242 habitants,
Casabianca située à 5 km, 80 habitants,
Poggio-Marinaccio située à 5 km, 24 habitants,
 entre autres. De plus, Castineta est située à seulement 17 km de Borgo.

Le parc de logements, à Castineta, se décomposait en 2011 en trois appartements et 61 maisons soit
 un marché plutôt équilibré.

Si vous envisagez de venir habiter à Castineta, vous pourrez facilement trouver une maison à acheter.


text:
void DisplayFunc(void);
void CleanUp(void);
void ExitClean(void);
void ResizeFunc(int NewWidth, int NewHeight);
void IdleFunc(void);
void KeyboardFunc(unsigned char cChar, int nMouseX, int nMouseY);
void KeySpecialFunc(int Key, int x, int y);
void SetPerspective(float Fov);
ILboolean Setup(void);
ILboolean GenSides(void);

ILboolean bCleaned = IL_FALSE;

GLuint TexID1 = 0, TexID2 = 0, TexID3 = 0, TexID4 = 0;
ILuint ImgId;

float Angle = 0.0, TransFactor = -470.0;

#endif //TEST3D_H


text:
<?php
require dirname(__FILE__) . '/../setup.php';

Loader::model("categories", "fastnews");

/**
 * Generated by PHPUnit_SkeletonGenerator 1.2.1 on 2014-06-29 at 20:33:38.
 */
class FastnewsCategoriesTest extends PHPUnit_Framework_TestCase {

    /**
     * @var FastnewsCategories
     */
    protected $object;

    /**
     * Sets up the fixture, for example, opens a network connection.
     * This method is called before a test is executed.
     */
    protected function setUp() {
        $this->object = new FastnewsCategories();
    }

    /**
     * Tears down the fixture, for example, closes a network connection.
     * This method is called after a test is executed.
     */
    protected function tearDown() {

    }

    public function testCategoriesCreate() {
        $this->object->save();

        $this->object->title = "Title";
        $this->object->description = "Content Some";

        $this->object->save();
        $this->id = $this->object->id;

        $this->object = FastnewsCategories::findOrFail($this->id);

        $this->assertEquals("Title", $this->object->title);
        $this->assertEquals("Content Some", $this->object->description);
    }

    /**
     * Testing DB edit functions
     * @depends testCategoriesCreate()
     */
    public function testCategoriesEdit() {
        $obj = FastnewsCategories;

        $obj->title = "TT";
        $obj->description = "con";

        $obj->save();
        $id = $obj->id;

        $obj = FastnewsCategories::findOrFail($id);

        $this->assertEquals("TT", $obj->title);
        $this->assertEquals("con", $obj->description);
    }

    /**
     * @depends testCategoriesEdit()
     */
    public function testCategoriesDelete() {
        $this->object->delete();
        $obj = FastnewsCategories::find($this->id);

        $this->assertEmpty($obj);
    }

    /**
     * Test Relations
     */
    public function testCategoriesRelations() {
        $c1 = FastnewsCategories::create(array("title" => "categorie1"));
        $c2 = FastnewsCategories::create(array("title" => "categorie2"));

        $n1 = FastnewsNews::create(array("title" => "news1", "categorie" => $c1->id));
        $n2 = FastnewsNews::create(array("title" => "news1", "categorie" => $c1->id));
        $n3 = FastnewsNews::create(array("title" => "news1", "categorie" => $c2->id));

        print_r($c1->news());
    }
}