<commit_before>#include "FAST/Data/Image.hpp"
#include "FAST/Data/Segmentation.hpp"
#include "FAST/Data/PointSet.hpp"
#include "FAST/DeviceManager.hpp"
#include "FAST/Algorithms/GaussianSmoothingFilter/GaussianSmoothingFilter.hpp"
#include "FAST/Algorithms/ImageGradient/ImageGradient.hpp"
#include "FAST/Utility.hpp"
#include "FAST/SceneGraph.hpp"
#include <boost/shared_array.hpp>
#include <queue>
#include "FAST/Algorithms/UltrasoundVesselDetection/UltrasoundVesselDetection.hpp"
#include "FAST/Algorithms/ImageCropper/ImageCropper.hpp"
#include "FAST/Algorithms/ImageClassifier/ImageClassifier.hpp"

namespace fast {

ProcessObjectPort UltrasoundVesselDetection::getOutputSegmentationPort() {
    mCreateSegmentation = true;
    return getOutputPort(0);
}

UltrasoundVesselDetection::UltrasoundVesselDetection() {
    createInputPort<Image>(0);
    createOutputPort<Segmentation>(0, OUTPUT_DEPENDS_ON_INPUT, 0);
    createOpenCLProgram(std::string(FAST_SOURCE_DIR) + "Algorithms/UltrasoundVesselDetection/UltrasoundVesselDetection.cl");
    mCreateSegmentation = false;

	mClassifier = ImageClassifier::New();
	std::string modelFile = "/home/smistad/vessel_detection_model/deploy.prototxt";
	std::string trainingFile = "/home/smistad/vessel_detection_model/snapshot_iter_540.caffemodel";
	std::string meanFile = "/home/smistad/vessel_detection_model/mean.binaryproto";
	//std::string modelFile = "/home/smistad/vessel_detection_model2/deploy.prototxt";
	//std::string trainingFile = "/home/smistad/vessel_detection_model2/snapshot_iter_960.caffemodel";
	//std::string meanFile = "/home/smistad/vessel_detection_model2/mean.binaryproto";
	mClassifier->loadModel(modelFile, trainingFile, meanFile);

}

struct Candidate {
	float score;
	VesselCrossSection::pointer crossSection;
};

class CandidateComparison {
	public:
		bool operator() (const Candidate& lhs, const Candidate& rhs) const {
			return lhs.score < rhs.score;
		}
};

void UltrasoundVesselDetection::execute() {
    Image::pointer input = getStaticInputData<Image>();
    if(input->getDimensions() != 2) {
        throw Exception("The UltrasoundVesselDetection algorithm is only for 2D");
    }

    mRuntimeManager->startRegularTimer("ellipse fitting");
    // Create kernel
    OpenCLDevice::pointer device = getMainDevice();
    cl::Program program = getOpenCLProgram(device);
    cl::Kernel kernel(program, "vesselDetection");

    // Run GaussianSmoothing on input
    GaussianSmoothingFilter::pointer filter = GaussianSmoothingFilter::New();
    filter->setInputData(input);
    filter->setMaskSize(5);
    filter->setStandardDeviation(2);
    filter->update();
    Image::pointer smoothedImage = filter->getOutputData<Image>();

    // Run ImageGradient on input
    ImageGradient::pointer imageGradient = ImageGradient::New();
    imageGradient->setInputConnection(filter->getOutputPort());
    imageGradient->update();
    Image::pointer gradients = imageGradient->getOutputData<Image>();
    OpenCLImageAccess::pointer inputImageAccess = input->getOpenCLImageAccess(ACCESS_READ, device);
    OpenCLImageAccess::pointer imageAccess = smoothedImage->getOpenCLImageAccess(ACCESS_READ, device);
    OpenCLImageAccess::pointer gradientAccess = gradients->getOpenCLImageAccess(ACCESS_READ, device);

    // Create output image
    boost::shared_array<float> zeros(new float[input->getWidth()*input->getHeight()*4]());
    cl::Image2D result = cl::Image2D(
            device->getContext(),
            CL_MEM_WRITE_ONLY | CL_MEM_COPY_HOST_PTR,
            getOpenCLImageFormat(device, CL_MEM_OBJECT_IMAGE2D, TYPE_FLOAT, 4),
            input->getWidth(), input->getHeight(),
			0,
			zeros.get()
    );

    // Run vessel detection kernel on smoothed image and gradient
    kernel.setArg(0, *inputImageAccess->get2DImage());
    kernel.setArg(1, *gradientAccess->get2DImage());
    kernel.setArg(2, result);
    kernel.setArg(3, input->getSpacing().x());

    const float minimumDepthInMM = 5;
    const float spacing = input->getSpacing().y();
    //const float maximumDepthInMM = 20;
    const float maximumDepthInMM = input->getHeight()*spacing*0.85;
    int startPosY = round(minimumDepthInMM/spacing);
    int endPosY = round(maximumDepthInMM/spacing);

    // Only process every second pixel
    cl::NDRange globalSize(input->getWidth() / 4, (endPosY-startPosY) / 4);
    cl::NDRange kernelOffset(0, startPosY / 4);
    device->getCommandQueue().enqueueNDRangeKernel(
            kernel,
			kernelOffset,
			globalSize,
            cl::NullRange
    );

    // Get result back
    boost::shared_array<float> data(new float[input->getWidth()*input->getHeight()*4]);
    device->getCommandQueue().enqueueReadImage(
            result,
            CL_TRUE,
            createOrigoRegion(),
            createRegion(input->getWidth(),input->getHeight(),1),
            0,0,
            data.get()
    );
    mRuntimeManager->stopRegularTimer("ellipse fitting");
    mRuntimeManager->startRegularTimer("candidate selection");

	AffineTransformation::pointer transform = SceneGraph::getAffineTransformationFromData(input);
	transform->scale(input->getSpacing());

    // Find best ellipses
	std::priority_queue<Candidate, std::vector<Candidate>, CandidateComparison> candidates;
	int startPosY2 = (startPosY / 4)*4;
    for(int x = 0; x < input->getWidth(); x+=4) {
        for(int y = startPosY2; y < endPosY; y+=4) {
            int i = x + y*input->getWidth();

            if(data[i*4] > 1.5) { // If score is higher than a threshold
				float posY = floor(data[i*4+3]/input->getWidth());
				float posX = data[i*4+3]-posY*input->getWidth();
				Vector3f voxelPosition(posX, posY, 0);
				Vector3f position = transform->multiply(voxelPosition);
				VesselCrossSection::pointer crossSection = VesselCrossSection::New();
				crossSection->create(position, voxelPosition.head(2), data[i*4+1], data[i*4+2]*data[i*4+1]);
				Candidate candidate;
				candidate.score = data[i*4];
				candidate.crossSection = crossSection;
				candidates.push(candidate);

            }
        }
    }

	mCrossSections.clear();
	// Go through all candidates
	while(!candidates.empty()) {
		Candidate next = candidates.top();
		candidates.pop();

		// Check if valid
		bool invalid = false;
		for(VesselCrossSection::pointer crossSection : mCrossSections) {
			VesselCrossSectionAccess::pointer access = crossSection->getAccess(ACCESS_READ);
			Vector2f imageCenter = access->getImageCenterPosition();
			float majorRadius = access->getMajorRadius();
			VesselCrossSectionAccess::pointer access2 = next.crossSection->getAccess(ACCESS_READ);
			if((access2->getImageCenterPosition() - imageCenter).norm() < majorRadius) {
				invalid = true;
				break;
			}
		}
		if(!invalid) {
			mCrossSections.push_back(next.crossSection);
		}
	}
	std::cout << mCrossSections.size() << " candidate vessels" << std::endl;<commit_msg>Code cleanup<commit_after>#include "FAST/Data/Image.hpp"
#include "FAST/Data/Segmentation.hpp"
#include "FAST/Data/PointSet.hpp"
#include "FAST/DeviceManager.hpp"
#include "FAST/Algorithms/GaussianSmoothingFilter/GaussianSmoothingFilter.hpp"
#include "FAST/Algorithms/ImageGradient/ImageGradient.hpp"
#include "FAST/Utility.hpp"
#include "FAST/SceneGraph.hpp"
#include <boost/shared_array.hpp>
#include <queue>
#include "FAST/Algorithms/UltrasoundVesselDetection/UltrasoundVesselDetection.hpp"
#include "FAST/Algorithms/ImageCropper/ImageCropper.hpp"
#include "FAST/Algorithms/ImageClassifier/ImageClassifier.hpp"

namespace fast {

ProcessObjectPort UltrasoundVesselDetection::getOutputSegmentationPort() {
    mCreateSegmentation = true;
    return getOutputPort(0);
}

UltrasoundVesselDetection::UltrasoundVesselDetection() {
    createInputPort<Image>(0);
    createOutputPort<Segmentation>(0, OUTPUT_DEPENDS_ON_INPUT, 0);
    createOpenCLProgram(std::string(FAST_SOURCE_DIR) + "Algorithms/UltrasoundVesselDetection/UltrasoundVesselDetection.cl");
    mCreateSegmentation = false;

	mClassifier = ImageClassifier::New();
	std::string modelFile = "/home/smistad/vessel_detection_model/deploy.prototxt";
	std::string trainingFile = "/home/smistad/vessel_detection_model/snapshot_iter_540.caffemodel";
	std::string meanFile = "/home/smistad/vessel_detection_model/mean.binaryproto";
	//std::string modelFile = "/home/smistad/vessel_detection_model2/deploy.prototxt";
	//std::string trainingFile = "/home/smistad/vessel_detection_model2/snapshot_iter_960.caffemodel";
	//std::string meanFile = "/home/smistad/vessel_detection_model2/mean.binaryproto";
	mClassifier->loadModel(modelFile, trainingFile, meanFile);

}

struct Candidate {
	float score;
	VesselCrossSection::pointer crossSection;
};

class CandidateComparison {
	public:
		bool operator() (const Candidate& lhs, const Candidate& rhs) const {
			return lhs.score < rhs.score;
		}
};

void UltrasoundVesselDetection::execute() {
    Image::pointer input = getStaticInputData<Image>();
    if(input->getDimensions() != 2) {
        throw Exception("The UltrasoundVesselDetection algorithm is only for 2D");
    }

    mRuntimeManager->startRegularTimer("ellipse fitting");
    // Create kernel
    OpenCLDevice::pointer device = getMainDevice();
    cl::Program program = getOpenCLProgram(device);
    cl::Kernel kernel(program, "vesselDetection");

    // Run GaussianSmoothing on input
    GaussianSmoothingFilter::pointer filter = GaussianSmoothingFilter::New();
    filter->setInputData(input);
    filter->setMaskSize(5);
    filter->setStandardDeviation(2);
    filter->update();
    Image::pointer smoothedImage = filter->getOutputData<Image>();

    // Run ImageGradient on input
    ImageGradient::pointer imageGradient = ImageGradient::New();
    imageGradient->setInputConnection(filter->getOutputPort());
    imageGradient->update();
    Image::pointer gradients = imageGradient->getOutputData<Image>();
    OpenCLImageAccess::pointer inputImageAccess = input->getOpenCLImageAccess(ACCESS_READ, device);
    OpenCLImageAccess::pointer imageAccess = smoothedImage->getOpenCLImageAccess(ACCESS_READ, device);
    OpenCLImageAccess::pointer gradientAccess = gradients->getOpenCLImageAccess(ACCESS_READ, device);

    // Create output image
    boost::shared_array<float> zeros(new float[input->getWidth()*input->getHeight()*4]());
    cl::Image2D result = cl::Image2D(
            device->getContext(),
            CL_MEM_WRITE_ONLY | CL_MEM_COPY_HOST_PTR,
            getOpenCLImageFormat(device, CL_MEM_OBJECT_IMAGE2D, TYPE_FLOAT, 4),
            input->getWidth(), input->getHeight(),
			0,
			zeros.get()
    );

    // Run vessel detection kernel on smoothed image and gradient
    kernel.setArg(0, *inputImageAccess->get2DImage());
    kernel.setArg(1, *gradientAccess->get2DImage());
    kernel.setArg(2, result);
    kernel.setArg(3, input->getSpacing().x());

    const float minimumDepthInMM = 5;
    const float spacing = input->getSpacing().y();
    //const float maximumDepthInMM = 20;
    const float maximumDepthInMM = input->getHeight()*spacing*0.85;
    int startPosY = round(minimumDepthInMM/spacing);
    int endPosY = round(maximumDepthInMM/spacing);

    // Only process every second pixel
    cl::NDRange globalSize(input->getWidth() / 4, (endPosY-startPosY) / 4);
    cl::NDRange kernelOffset(0, startPosY / 4);
    device->getCommandQueue().enqueueNDRangeKernel(
            kernel,
			kernelOffset,
			globalSize,
            cl::NullRange
    );

    // Get result back
    boost::shared_array<float> data(new float[input->getWidth()*input->getHeight()*4]);
    device->getCommandQueue().enqueueReadImage(
            result,
            CL_TRUE,
            createOrigoRegion(),
            createRegion(input->getWidth(),input->getHeight(),1),
            0,0,
            data.get()
    );
    mRuntimeManager->stopRegularTimer("ellipse fitting");
    mRuntimeManager->startRegularTimer("candidate selection");

	AffineTransformation::pointer transform = SceneGraph::getAffineTransformationFromData(input);
	transform->scale(input->getSpacing());

    // Find best ellipses
	std::priority_queue<Candidate, std::vector<Candidate>, CandidateComparison> candidates;
	int startPosY2 = (startPosY / 4)*4;
    for(int x = 0; x < input->getWidth(); x+=4) {
        for(int y = startPosY2; y < endPosY; y+=4) {
            int i = x + y*input->getWidth();

            if(data[i*4] > 1.5) { // If score is higher than a threshold
				float posY = floor(data[i*4+3]/input->getWidth());
				float posX = data[i*4+3]-posY*input->getWidth();
				Vector3f voxelPosition(posX, posY, 0);
				Vector3f position = transform->multiply(voxelPosition);
				VesselCrossSection::pointer crossSection = VesselCrossSection::New();
				crossSection->create(position, voxelPosition.head(2), data[i*4+1], data[i*4+2]*data[i*4+1]);
				Candidate candidate;
				candidate.score = data[i*4];
				candidate.crossSection = crossSection;
				candidates.push(candidate);

            }
        }
    }

	mCrossSections.clear();
	// Go through all candidates
	while(!candidates.empty()) {
		Candidate next = candidates.top();
		candidates.pop();

		// Check if valid
		bool invalid = false;
		VesselCrossSectionAccess::pointer access2 = next.crossSection->getAccess(ACCESS_READ);
		Vector2f candidateImageCenter = access2->getImageCenterPosition();
		for(VesselCrossSection::pointer crossSection : mCrossSections) {
			VesselCrossSectionAccess::pointer access = crossSection->getAccess(ACCESS_READ);
			Vector2f imageCenter = access->getImageCenterPosition();
			float majorRadius = access->getMajorRadius();
			if((candidateImageCenter - imageCenter).norm() < majorRadius) {
				invalid = true;
				break;
			}
		}
		if(!invalid) {
			mCrossSections.push_back(next.crossSection);
		}
	}
	std::cout << mCrossSections.size() << " candidate vessels" << std::endl;
{'fraction_non_alphanumeric': 0.12964037122969838, 'fraction_numerical': 0.010368329466357308, 'mean_word_length': 4.5955375253549695, 'pattern_counts': {'":': 0, '<': 45, '<?xml version=': 0, '>': 121, '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': '1119149', 'n_tokens_mistral': 4567, 'n_tokens_neox': 4213, 'n_words': 721}