<commit_before>/*
 *  qisplit.cpp
 *
 *  Copyright (c) 2015 <NAME>.
 *
 *  This Source Code Form is subject to the terms of the Mozilla Public
 *  License, v. 2.0. If a copy of the MPL was not distributed with this
 *  file, You can obtain one at http://mozilla.org/MPL/2.0/.
 *
 */

#include <getopt.h>
#include <string>
#include <sstream>
#include <iostream>
#include <iomanip>

#include "Util.h"
#include "Types.h"

#include "itkRescaleIntensityImageFilter.h"
#include "itkOtsuThresholdImageFilter.h"
#include "itkConnectedComponentImageFilter.h"
#include "itkLabelShapeKeepNObjectsImageFilter.h"
#include "itkRelabelComponentImageFilter.h"
#include "itkGrayscaleFillholeImageFilter.h"
#include "itkLabelStatisticsImageFilter.h"
#include "itkLabelImageToLabelMapFilter.h"
#include "itkLabelMapMaskImageFilter.h"
#include "itkExtractImageFilter.h"
#include "itkImageMomentsCalculator.h"
#include "itkEuler3DTransform.h"
#include "itkTransformFileWriter.h"
#include "itkResampleImageFilter.h"
#include "itkLinearInterpolateImageFunction.h"
#include "itkNearestNeighborInterpolateImageFunction.h"
#include "itkShrinkImageFilter.h"
#include "itkSmoothingRecursiveGaussianImageFilter.h"
#include "itkMeanSquaresImageToImageMetric.h"
#include "itkRegularStepGradientDescentOptimizer.h"
#include "itkImageRegistrationMethod.h"

using namespace std;

typedef unsigned int TLabel;
typedef itk::Image<TLabel, 3> TLabelImage;

typename TLabelImage::Pointer ThresholdMask(const QI::ImageF::Pointer &img, const float thresh) {
    typedef itk::BinaryThresholdImageFilter<QI::ImageF, TLabelImage> TThreshFilter;
    auto threshold = TThreshFilter::New();
    threshold->SetInput(img);
    threshold->SetLowerThreshold(thresh);
    threshold->SetUpperThreshold(numeric_limits<float>::infinity());
    threshold->SetInsideValue(1);
    threshold->SetOutsideValue(0);
    threshold->Update();
    typename TLabelImage::Pointer mask = threshold->GetOutput();
    mask->DisconnectPipeline();
    return mask;
}

typename TLabelImage::Pointer OtsuMask(const QI::ImageF::Pointer &img) {
    auto otsuFilter = itk::OtsuThresholdImageFilter<QI::ImageF, TLabelImage>::New();
    otsuFilter->SetInput(img);
    otsuFilter->SetOutsideValue(1);
    otsuFilter->SetInsideValue(0);
    otsuFilter->Update();
    typename TLabelImage::Pointer mask = otsuFilter->GetOutput();
    mask->DisconnectPipeline();
    return mask;
}

typename TLabelImage::Pointer FindLabels(const TLabelImage::Pointer &mask, const int size_threshold, int &keep) {
    auto CC = itk::ConnectedComponentImageFilter<TLabelImage, TLabelImage>::New();
    auto relabel = itk::RelabelComponentImageFilter<TLabelImage, TLabelImage>::New();
    CC->SetInput(mask);
    relabel->SetInput(CC->GetOutput());
    relabel->Update();
    // Relabel sorts on size by default, so now work out how many make the size threshold
    auto label_sizes = relabel->GetSizeOfObjectsInPixels();
    if (keep > label_sizes.size())
        keep = label_sizes.size();
    for (int i = 0; i < keep; i++) {
        if (label_sizes[i] < size_threshold) {
            keep = i;
            break;
        }
    }
    if (keep == 0) {
        throw(runtime_error("No labels found in mask"));
    }

    typedef itk::LabelShapeKeepNObjectsImageFilter<TLabelImage> TKeepN;
    TKeepN::Pointer keepN = TKeepN::New();
    keepN->SetInput(relabel->GetOutput());
    keepN->SetBackgroundValue(0);
    keepN->SetNumberOfObjects(keep);
    keepN->SetAttribute(TKeepN::LabelObjectType::NUMBER_OF_PIXELS);
    
    typedef itk::GrayscaleFillholeImageFilter<TLabelImage, TLabelImage> TFill;
    TFill::Pointer fill = TFill::New();
    fill->SetInput(keepN->GetOutput());
    fill->Update();
    
    typename TLabelImage::Pointer labels = fill->GetOutput();
    labels->DisconnectPipeline();
    return labels;
}

QI::ImageF::Pointer MaskWithLabel(const QI::ImageF::Pointer &image, const TLabelImage::Pointer &labels, const int l, const bool crop) {
    // convert the label image into a LabelMap
    typedef itk::LabelMap<itk::LabelObject<TLabel, 3>> TLabelMap;
    auto convert = itk::LabelImageToLabelMapFilter<TLabelImage, TLabelMap> ::New();
    convert->SetInput(labels);
    auto masker = itk::LabelMapMaskImageFilter<TLabelMap, QI::ImageF>::New();
    masker->SetInput(convert->GetOutput());
    masker->SetFeatureImage(image);
    masker->SetLabel(l);
    masker->SetBackgroundValue(0.);
    masker->SetNegated(false); // Mask outside the mask
    if (crop)
        masker->SetCrop(true);
    masker->Update();
    QI::ImageF::Pointer masked = masker->GetOutput();
    masked->DisconnectPipeline();
    return masked;
}

typedef itk::ImageMomentsCalculator<QI::ImageF> TMoments;
TMoments::VectorType GetCoG(const QI::ImageF::Pointer &img) {
    auto moments = TMoments::New();
    moments->SetImage(img);
    moments->Compute();
    return moments->GetCenterOfGravity();
}

typedef itk::Euler3DTransform<double> TRigid;
template<typename TImg, typename TInterpFunction>
typename TImg::Pointer ResampleImage(const typename TImg::Pointer &image, const typename TRigid::Pointer &tfm, const QI::ImageF::Pointer &reference = ITK_NULLPTR) {
    // Set these up to use in the loop
    typedef itk::ResampleImageFilter<TImg, TImg, double> TResampler;
    typedef TInterpFunction TInterpolator;
    typename TInterpolator::Pointer interp = TInterpolator::New();
    interp->SetInputImage(image);
    typename TResampler::Pointer resamp = TResampler::New();
    resamp->SetInput(image);
    resamp->SetInterpolator(interp);
    resamp->SetDefaultPixelValue(0.);
    resamp->SetTransform(tfm);
    if (reference)
        resamp->SetOutputParametersFromImage(reference);
    else
        resamp->SetOutputParametersFromImage(image);
    resamp->Update();
    typename TImg::Pointer rimage = resamp->GetOutput();
    rimage->DisconnectPipeline();
    return rimage;
}

void RegisterImageToReference(const QI::ImageF::Pointer &image, const QI::ImageF::Pointer &reference, TRigid::Pointer tfm, const int shrink) {
    typedef itk::SmoothingRecursiveGaussianImageFilter<QI::ImageF, QI::ImageF> TSmooth;
    typedef itk::ShrinkImageFilter<QI::ImageF, QI::ImageF> TShrink;
    typedef itk::RegularStepGradientDescentOptimizer TOpt;
    typedef itk::MeanSquaresImageToImageMetric<QI::ImageF, QI::ImageF> TMetric;
    typedef itk::ImageRegistrationMethod<QI::ImageF, QI::ImageF> TReg;
    typedef itk::LinearInterpolateImageFunction<QI::ImageF, double> TInterp;
    
    TSmooth::Pointer smooth_img = TSmooth::New();
    TSmooth::Pointer smooth_ref = TSmooth::New();
    TSmooth::SigmaArrayType sigma = image->GetSpacing();
    for (int i = 0; i < sigma.Size(); i++)
        sigma[i] = sigma[i] * shrink;
    smooth_img->SetInput(image);
    smooth_img->SetSigmaArray(sigma);
    smooth_ref->SetInput(reference);
    smooth_ref->SetSigmaArray(sigma);
    
    TShrink::Pointer shrink_img = TShrink::New();
    TShrink::Pointer shrink_ref = TShrink::New();
    shrink_img->SetInput(smooth_img->GetOutput());
    shrink_img->SetShrinkFactors(shrink);
    shrink_ref->SetInput(smooth_ref->GetOutput());
    shrink_ref->SetShrinkFactors(shrink);
    
    TMetric::Pointer metric = TMetric::New();
    TInterp::Pointer interp = TInterp::New();
    TOpt::Pointer opt = TOpt::New();
    TReg::Pointer reg = TReg::New();
    
    reg->SetMetric(metric);
    reg->SetOptimizer(opt);
    reg->SetTransform(tfm);
    reg->SetInterpolator(interp);<commit_msg>Changed to Mutual Information as the metric.<commit_after>/*
 *  qisplit.cpp
 *
 *  Copyright (c) 2015 <NAME>.
 *
 *  This Source Code Form is subject to the terms of the Mozilla Public
 *  License, v. 2.0. If a copy of the MPL was not distributed with this
 *  file, You can obtain one at http://mozilla.org/MPL/2.0/.
 *
 */

#include <getopt.h>
#include <string>
#include <sstream>
#include <iostream>
#include <iomanip>

#include "Util.h"
#include "Types.h"

#include "itkRescaleIntensityImageFilter.h"
#include "itkOtsuThresholdImageFilter.h"
#include "itkConnectedComponentImageFilter.h"
#include "itkLabelShapeKeepNObjectsImageFilter.h"
#include "itkRelabelComponentImageFilter.h"
#include "itkGrayscaleFillholeImageFilter.h"
#include "itkLabelStatisticsImageFilter.h"
#include "itkLabelImageToLabelMapFilter.h"
#include "itkLabelMapMaskImageFilter.h"
#include "itkExtractImageFilter.h"
#include "itkImageMomentsCalculator.h"
#include "itkEuler3DTransform.h"
#include "itkTransformFileWriter.h"
#include "itkResampleImageFilter.h"
#include "itkLinearInterpolateImageFunction.h"
#include "itkNearestNeighborInterpolateImageFunction.h"
#include "itkShrinkImageFilter.h"
#include "itkSmoothingRecursiveGaussianImageFilter.h"
#include "itkMattesMutualInformationImageToImageMetric.h"
#include "itkRegularStepGradientDescentOptimizer.h"
#include "itkImageRegistrationMethod.h"

using namespace std;

typedef unsigned int TLabel;
typedef itk::Image<TLabel, 3> TLabelImage;

typename TLabelImage::Pointer ThresholdMask(const QI::ImageF::Pointer &img, const float thresh) {
    typedef itk::BinaryThresholdImageFilter<QI::ImageF, TLabelImage> TThreshFilter;
    auto threshold = TThreshFilter::New();
    threshold->SetInput(img);
    threshold->SetLowerThreshold(thresh);
    threshold->SetUpperThreshold(numeric_limits<float>::infinity());
    threshold->SetInsideValue(1);
    threshold->SetOutsideValue(0);
    threshold->Update();
    typename TLabelImage::Pointer mask = threshold->GetOutput();
    mask->DisconnectPipeline();
    return mask;
}

typename TLabelImage::Pointer OtsuMask(const QI::ImageF::Pointer &img) {
    auto otsuFilter = itk::OtsuThresholdImageFilter<QI::ImageF, TLabelImage>::New();
    otsuFilter->SetInput(img);
    otsuFilter->SetOutsideValue(1);
    otsuFilter->SetInsideValue(0);
    otsuFilter->Update();
    typename TLabelImage::Pointer mask = otsuFilter->GetOutput();
    mask->DisconnectPipeline();
    return mask;
}

typename TLabelImage::Pointer FindLabels(const TLabelImage::Pointer &mask, const int size_threshold, int &keep) {
    auto CC = itk::ConnectedComponentImageFilter<TLabelImage, TLabelImage>::New();
    auto relabel = itk::RelabelComponentImageFilter<TLabelImage, TLabelImage>::New();
    CC->SetInput(mask);
    relabel->SetInput(CC->GetOutput());
    relabel->Update();
    // Relabel sorts on size by default, so now work out how many make the size threshold
    auto label_sizes = relabel->GetSizeOfObjectsInPixels();
    if (keep > label_sizes.size())
        keep = label_sizes.size();
    for (int i = 0; i < keep; i++) {
        if (label_sizes[i] < size_threshold) {
            keep = i;
            break;
        }
    }
    if (keep == 0) {
        throw(runtime_error("No labels found in mask"));
    }

    typedef itk::LabelShapeKeepNObjectsImageFilter<TLabelImage> TKeepN;
    TKeepN::Pointer keepN = TKeepN::New();
    keepN->SetInput(relabel->GetOutput());
    keepN->SetBackgroundValue(0);
    keepN->SetNumberOfObjects(keep);
    keepN->SetAttribute(TKeepN::LabelObjectType::NUMBER_OF_PIXELS);
    
    typedef itk::GrayscaleFillholeImageFilter<TLabelImage, TLabelImage> TFill;
    TFill::Pointer fill = TFill::New();
    fill->SetInput(keepN->GetOutput());
    fill->Update();
    
    typename TLabelImage::Pointer labels = fill->GetOutput();
    labels->DisconnectPipeline();
    return labels;
}

QI::ImageF::Pointer MaskWithLabel(const QI::ImageF::Pointer &image, const TLabelImage::Pointer &labels, const int l, const bool crop) {
    // convert the label image into a LabelMap
    typedef itk::LabelMap<itk::LabelObject<TLabel, 3>> TLabelMap;
    auto convert = itk::LabelImageToLabelMapFilter<TLabelImage, TLabelMap> ::New();
    convert->SetInput(labels);
    auto masker = itk::LabelMapMaskImageFilter<TLabelMap, QI::ImageF>::New();
    masker->SetInput(convert->GetOutput());
    masker->SetFeatureImage(image);
    masker->SetLabel(l);
    masker->SetBackgroundValue(0.);
    masker->SetNegated(false); // Mask outside the mask
    if (crop)
        masker->SetCrop(true);
    masker->Update();
    QI::ImageF::Pointer masked = masker->GetOutput();
    masked->DisconnectPipeline();
    return masked;
}

typedef itk::ImageMomentsCalculator<QI::ImageF> TMoments;
TMoments::VectorType GetCoG(const QI::ImageF::Pointer &img) {
    auto moments = TMoments::New();
    moments->SetImage(img);
    moments->Compute();
    return moments->GetCenterOfGravity();
}

typedef itk::Euler3DTransform<double> TRigid;
template<typename TImg, typename TInterpFunction>
typename TImg::Pointer ResampleImage(const typename TImg::Pointer &image, const typename TRigid::Pointer &tfm, const QI::ImageF::Pointer &reference = ITK_NULLPTR) {
    // Set these up to use in the loop
    typedef itk::ResampleImageFilter<TImg, TImg, double> TResampler;
    typedef TInterpFunction TInterpolator;
    typename TInterpolator::Pointer interp = TInterpolator::New();
    interp->SetInputImage(image);
    typename TResampler::Pointer resamp = TResampler::New();
    resamp->SetInput(image);
    resamp->SetInterpolator(interp);
    resamp->SetDefaultPixelValue(0.);
    resamp->SetTransform(tfm);
    if (reference)
        resamp->SetOutputParametersFromImage(reference);
    else
        resamp->SetOutputParametersFromImage(image);
    resamp->Update();
    typename TImg::Pointer rimage = resamp->GetOutput();
    rimage->DisconnectPipeline();
    return rimage;
}

void RegisterImageToReference(const QI::ImageF::Pointer &image, const QI::ImageF::Pointer &reference, TRigid::Pointer tfm, const int shrink) {
    typedef itk::SmoothingRecursiveGaussianImageFilter<QI::ImageF, QI::ImageF> TSmooth;
    typedef itk::ShrinkImageFilter<QI::ImageF, QI::ImageF> TShrink;
    typedef itk::RegularStepGradientDescentOptimizer TOpt;
    typedef itk::MattesMutualInformationImageToImageMetric<QI::ImageF, QI::ImageF> TMetric;
    typedef itk::ImageRegistrationMethod<QI::ImageF, QI::ImageF> TReg;
    typedef itk::LinearInterpolateImageFunction<QI::ImageF, double> TInterp;
    
    TSmooth::Pointer smooth_img = TSmooth::New();
    TSmooth::Pointer smooth_ref = TSmooth::New();
    TSmooth::SigmaArrayType sigma = image->GetSpacing();
    for (int i = 0; i < sigma.Size(); i++)
        sigma[i] = sigma[i] * shrink;
    smooth_img->SetInput(image);
    smooth_img->SetSigmaArray(sigma);
    smooth_ref->SetInput(reference);
    smooth_ref->SetSigmaArray(sigma);
    
    TShrink::Pointer shrink_img = TShrink::New();
    TShrink::Pointer shrink_ref = TShrink::New();
    shrink_img->SetInput(smooth_img->GetOutput());
    shrink_img->SetShrinkFactors(shrink);
    shrink_ref->SetInput(smooth_ref->GetOutput());
    shrink_ref->SetShrinkFactors(shrink);
    
    TMetric::Pointer metric = TMetric::New();
    TInterp::Pointer interp = TInterp::New();
    TOpt::Pointer opt = TOpt::New();
    TReg::Pointer reg = TReg::New();
    
    metric->SetNumberOfHistogramBins(32);
    metric->SetNumberOfSpatialSamples(10000);
    
    reg->SetMetric(metric);
    reg->SetOptimizer(opt);
    reg->SetTransform(tfm);
    reg->SetInterpolator(interp);
{'fraction_non_alphanumeric': 0.1299228312932411, 'fraction_numerical': 0.0033927621075039915, 'mean_word_length': 5.083771752326993, 'pattern_counts': {'":': 0, '<': 63, '<?xml version=': 0, '>': 197, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4707863', 'n_tokens_mistral': 5000, 'n_tokens_neox': 4820, 'n_words': 929}