ITK 计算均方差与 马特斯互信息

const unsigned int   Dimension = 2;//const
    typedef   unsigned char  PixelType;

    typedef itk::Image< PixelType, Dimension >   ImageType;

    using ImageIOType = itk::GDCMImageIO;
    ImageIOType::Pointer gdcmImageIO = ImageIOType::New();

    typedef itk::ImageFileReader< ImageType >  ReaderType;

    ReaderType::Pointer fixedReader = ReaderType::New();
    ReaderType::Pointer movingReader = ReaderType::New();

    //const char* filenameFixed = "D:\\DicomFile\\PID089844Test\\test\\CT003001.dcm";
    //const char* filenameMove = "D:\\DicomFile\\PID089844Test\\test\\CT003001.dcm";

    CString cstringfilenameFixed = _T(""),cstringfileMove = _T("");
    ((CEdit*)GetDlgItem(IDC_EDITFixImgPath))->GetWindowTextW(cstringfilenameFixed);
    ((CEdit*)GetDlgItem(IDC_EDITMovingImgPath))->GetWindowTextW(cstringfileMove);

    char filenameFixed[1024], filenameMove[1024];
    wsprintfA(filenameFixed, "%ls", cstringfilenameFixed);
    wsprintfA(filenameMove, "%ls", cstringfileMove);

    fixedReader->SetFileName(filenameFixed);
    //fixedReader->SetFileName("D:\\DicomFile\\PID089844Test\\test\\CT003001.dcm");
    fixedReader->SetImageIO(gdcmImageIO);

    movingReader->SetFileName(filenameMove);
    //movingReader->SetFileName("D:\\DicomFile\\PID089844Test\\test\\CT003001.dcm");
    movingReader->SetImageIO(gdcmImageIO);

    try
    {
        fixedReader->Update();
        movingReader->Update();
    }
    catch (itk::ExceptionObject& excep)
    {
        CString sMsg;
        sMsg.Format(_T("Reader->Update() error! \r\n"));
        PrintStatus(sMsg);
        return;
    }

    typedef itk::TranslationTransform< double, Dimension >  TransformType;

    TransformType::Pointer transform = TransformType::New();

    typedef itk::NearestNeighborInterpolateImageFunction<
        ImageType, double >  InterpolatorType;

    InterpolatorType::Pointer interpolator = InterpolatorType::New();

    transform->SetIdentity();

    ImageType::ConstPointer fixedImage = fixedReader->GetOutput();
    ImageType::ConstPointer movingImage = movingReader->GetOutput();

    if (((CButton*)GetDlgItem(IDC_RADIO_MeanSquares))->GetCheck())
    {
        typedef itk::MeanSquaresImageToImageMetric<
            ImageType, ImageType >  MetricType;
        MetricType::Pointer metric = MetricType::New();
        metric->SetTransform(transform);
        metric->SetInterpolator(interpolator);
        metric->SetFixedImage(fixedImage);
        metric->SetMovingImage(movingImage);
        metric->SetFixedImageRegion(fixedImage->GetBufferedRegion());

        try
        {
            metric->Initialize();
        }
        catch (itk::ExceptionObject& excep)
        {
            CString sMsg;
            sMsg.Format(_T("metric->Initialize() error! \r\n"));
            PrintStatus(sMsg);
            return;
        }

        MetricType::TransformParametersType displacement(Dimension);

        const int rangex = 1;
        const int rangey = 1;

        for (int dx = 0; dx < rangex; dx++)
        {
            for (int dy = 0; dy < rangey; dy++)
            {
                displacement[0] = dx;
                displacement[1] = dy;
                const double value = metric->GetValue(displacement);
                CString sMsg;
                sMsg.Format(_T("均方差值:%f   \r\n"), value);
                PrintStatus(sMsg);

            }
        }
    }
    else if (((CButton*)GetDlgItem(IDC_RADIO_MattesMutualInformation))->GetCheck())
    {

        typedef itk::MattesMutualInformationImageToImageMetric<
            ImageType, ImageType >  MetricType;

        MetricType::Pointer metric = MetricType::New();
        metric->SetTransform(transform);
        metric->SetInterpolator(interpolator);

        metric->SetFixedImage(fixedImage);
        metric->SetMovingImage(movingImage);

        metric->SetFixedImageRegion(fixedImage->GetBufferedRegion());

        try
        {
            metric->Initialize();
        }
        catch (itk::ExceptionObject& excep)
        {
            CString sMsg;
            sMsg.Format(_T("metric->Initialize() error! \r\n"));
            PrintStatus(sMsg);
            return;
        }

        MetricType::TransformParametersType displacement(Dimension);

        const int rangex = 1;
        const int rangey = 1;

        for (int dx = 0; dx < rangex; dx++)
        {
            for (int dy = 0; dy < rangey; dy++)
            {
                displacement[0] = dx;
                displacement[1] = dy;
                const double value = metric->GetValue(displacement);
                CString sMsg;
                sMsg.Format(_T("马特斯互信息:%f   \r\n"), value);
                //sMsg.Format(_T("马特斯互信息:%f  dx:%d dy:%d \r\n"), value, dx, dy);
                PrintStatus(sMsg);

            }
        }
    }

    return;