cxSyntheticVolume.cpp

Go to the documentation of this file.

/*=========================================================================
This file is part of CustusX, an Image Guided Therapy Application.

Copyright (c) 2008-2014, SINTEF Department of Medical Technology
All rights reserved.

Redistribution and use in source and binary forms, with or without 
modification, are permitted provided that the following conditions are met:

1. Redistributions of source code must retain the above copyright notice, 
   this list of conditions and the following disclaimer.

2. Redistributions in binary form must reproduce the above copyright notice, 
   this list of conditions and the following disclaimer in the documentation 
   and/or other materials provided with the distribution.

3. Neither the name of the copyright holder nor the names of its contributors 
   may be used to endorse or promote products derived from this software 
   without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE 
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR 
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER 
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, 
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=========================================================================*/

#include "cxSyntheticVolume.h"
#include "vtkImageData.h"
#include "cxImage.h"
#include <cstdlib>
#include <time.h>
#include "cxTypeConversions.h"
#include <QTime>
#include "cxLogger.h"
#include "cxRegistrationTransform.h"
#include "cxVolumeHelpers.h"


double noiseValue(double noiseSigma,
                         double noiseMean)
{
        double random_value_1 = (rand()+1.0)/(RAND_MAX+1.0);
        double random_value_2 = (rand()+1.0)/(RAND_MAX+1.0);

        double random_normal = sqrt(-2*log(random_value_1)) * cos(2*M_PI*random_value_2);

        return random_normal*noiseSigma + noiseMean;
}

namespace cx {


ProcessedUSInputDataPtr
cxSyntheticVolume::sampleUsData(const std::vector<Transform3D>& planes_rMt,
                         const ProbeDefinition& probe,
                         const Transform3D& output_dMr,
                         const double noiseSigma,
                         const unsigned char noiseMean) const
{
        cx::ProbeSector sector;
        sector.setData(probe);

        std::vector<Transform3D> planes_rMf(planes_rMt.size());
        for (unsigned i=0; i<planes_rMt.size(); ++i)
                planes_rMf[i] = planes_rMt[i] * sector.get_tMu() * sector.get_uMv();

        Eigen::Array2f pixelSpacing = probe.getSpacing().block(0,0,2,1).cast<float>();
        Eigen::Array2i sliceDimension(probe.getSize().width(), probe.getSize().height());

        return this->sampleUsData(planes_rMf,
                                                          pixelSpacing,
                                                          sliceDimension,
                                                          output_dMr,
                                                          noiseSigma, noiseMean);
}

vtkImageDataPtr
cxSyntheticVolume::sampleUsData(const Transform3D& plane_rMt,
                         const ProbeDefinition& probe,
                         const double noiseSigma,
                         const unsigned char noiseMean) const
{
        cx::ProbeSector sector;
        sector.setData(probe);

        Transform3D rMf = plane_rMt * sector.get_tMu() * sector.get_uMv();
        Eigen::Array2f pixelSpacing = probe.getSpacing().block(0,0,2,1).cast<float>();
        Eigen::Array2i sliceDimension(probe.getSize().width(), probe.getSize().height());

        return this->sampleUsData(rMf,
                                                          pixelSpacing,
                                                          sliceDimension,
                                                          noiseSigma, noiseMean);
}

ProcessedUSInputDataPtr
cxSyntheticVolume::sampleUsData(const std::vector<Transform3D>& planes_rMf,
                                const Eigen::Array2f& pixelSpacing,
                                const Eigen::Array2i& sliceDimension,
                                                                const Transform3D& output_dMr,
                                                                const double noiseSigma, const unsigned char noiseMean) const
{
        std::vector<TimedPosition> positions;
        std::vector<vtkImageDataPtr> images;
        // For each plane
        for(std::vector<Transform3D>::const_iterator i = planes_rMf.begin();
                planes_rMf.end() != i;
                        ++i)
        {
                const Transform3D rMf = *i;

                vtkImageDataPtr us_frame;
                us_frame = this->sampleUsData(rMf, pixelSpacing, sliceDimension, noiseSigma, noiseMean);

                // Build the TimedPosition for this frame
                TimedPosition t;
                t.mTime = i - planes_rMf.begin();
                t.mPos = output_dMr*rMf;

                positions.push_back(t);
                images.push_back(us_frame);
        }

        vtkImageDataPtr mask = this->createEmptyMask(sliceDimension);
//      std::cout << "elapsed: " << time.elapsed() << std::endl;

        ProcessedUSInputDataPtr ret;
        ret.reset(new ProcessedUSInputData(images, positions, mask, "VIRTUAL_DATA", "VIRTUAL_DATA_"));
        return ret;
}

vtkImageDataPtr
cxSyntheticVolume::sampleUsData(const Transform3D& rMf,
                                                                const Eigen::Array2f& pixelSpacing,
                                                                const Eigen::Array2i& sliceDimension,
                                                                const double noiseSigma, const unsigned char noiseMean) const
{


        const Vector3D p0 = rMf.coord(Vector3D(0,0,0));
        const Vector3D e_x = rMf.vector(Vector3D(pixelSpacing[0],0,0));
        const Vector3D e_y = rMf.vector(Vector3D(0,pixelSpacing[1],0));

        vtkImageDataPtr us_frame = vtkImageDataPtr::New();
        us_frame->SetExtent(0, sliceDimension[0]-1, 0, sliceDimension[1]-1, 0, 0);
        us_frame->SetSpacing(pixelSpacing[0], pixelSpacing[1], 0.0);
        us_frame->AllocateScalars(VTK_UNSIGNED_CHAR, 1);

        unsigned char* us_data = (unsigned char*)us_frame->GetScalarPointer();
        // For each pixel on that plane
        for(unsigned int px = 0; px < sliceDimension[0]; px++)
        {
                // optimization: use transformed pixel vectors
                const Vector3D px0_vol = p0 + e_x*px;

                for(unsigned int py = 0; py < sliceDimension[1]; py++)
                {
                        const Vector3D volume_coords = px0_vol + e_y*py;

                        // Evaluate volume at that position
                        const unsigned char val = this->evaluate(volume_coords);

                        const double noise_val = noiseValue(noiseSigma, noiseMean);
                        const int noised_val = noise_val + val;
                        unsigned char final_val = this->constrainToUnsignedChar(noised_val);
                        // Store that value in the US slice
                        us_data[px + py*sliceDimension[0]] = final_val;

                }
        }

        setDeepModified(us_frame);
        return us_frame;
}

vtkImageDataPtr cxSyntheticVolume::createEmptyMask(const Eigen::Array2i& sliceDimension) const
{
        vtkImageDataPtr mask = vtkImageDataPtr::New();
        mask->SetExtent(0, sliceDimension[0]-1, 0, sliceDimension[1]-1, 0, 0);
        mask->AllocateScalars(VTK_UNSIGNED_CHAR, 1);
        unsigned char* mask_data = (unsigned char*)mask->GetScalarPointer();
        memset(mask_data, 1, sizeof(unsigned char)*sliceDimension[0]*sliceDimension[1]);
        setDeepModified(mask);
        return mask;
}

unsigned char cxSyntheticVolume::constrainToUnsignedChar(const int val) const
{
        if(val < 0)
        {
                return 0;
        }
        else if(val > 255)
        {
                return 255;
        }
        else
        {
                return (unsigned char)val;
        }
}

float cxSyntheticVolume::computeRMSError(ImagePtr vol)
{
        vtkImageDataPtr input = vol->getBaseVtkImageData();

        vtkImageDataPtr nominal = vtkImageDataPtr::New();
        nominal->DeepCopy(input);
        cx::ImagePtr nominal_img(new cx::Image("nominal", nominal));
        nominal_img->get_rMd_History()->setRegistration(vol->get_rMd());
        this->fillVolume(nominal_img);

        return calculateRMSError(input, nominal);
}

void cxSyntheticVolume::fillVolume(cx::ImagePtr vol)
{
        vtkImageDataPtr raw = vol->getBaseVtkImageData();
        cx::Transform3D rMd = vol->get_rMd();

        Eigen::Array3i dims = Eigen::Array3i(raw->GetDimensions());
        Eigen::Array3d spacing = Eigen::Array3d(raw->GetSpacing());
        unsigned char *pixels = (unsigned char*)raw->GetScalarPointer();

        for(int z = 0; z < dims[2]; ++z)
        {
                for(int y = 0; y < dims[1]; ++y)
                {
                        for(int x = 0; x < dims[0]; ++x)
                        {
                                Vector3D p_d = Vector3D(x, y, z).array()*spacing;
                                int index = x + y*dims[0] + z*dims[1]*dims[0];
                                pixels[index] = this->evaluate(rMd.coord(p_d));
                        }
                }
        }
        setDeepModified(raw);
}

double calculateRMSError(vtkImageDataPtr a, vtkImageDataPtr b)
{
        CX_ASSERT(Eigen::Array3i(a->GetDimensions()).isApprox(Eigen::Array3i(b->GetDimensions())));
        CX_ASSERT(Eigen::Array3d(a->GetSpacing()).isApprox(Eigen::Array3d(b->GetSpacing())));

        float sse = 0.0f;
        Eigen::Array3i dims = Eigen::Array3i(a->GetDimensions());
        unsigned char *pa = (unsigned char*)a->GetScalarPointer();
        unsigned char *pb = (unsigned char*)b->GetScalarPointer();

        for(int z = 0; z < dims[2]; ++z)
        {
                for(int y = 0; y < dims[1]; ++y)
                {
                        for(int x = 0; x < dims[0]; ++x)
                        {
                                int index = x + y*dims[0] + z*dims[1]*dims[0];
                                float error = pa[index] - pb[index];
                                sse += error*error;
                        }
                }
        }

        return sqrt(sse/dims.prod());
}

struct MassFunctor
{
        MassFunctor()
        {
                mSum = 0;
                mThreshold = 2;
        }

        void operator()(int x, int y, int z, unsigned char* val)
        {
                if (*val < mThreshold)
                        return;
                mSum += *val;
        }

        double getVolume()
        {
                return mSum;
        }

private:
        double mSum;
        unsigned char mThreshold;
};

struct CentroidFunctor
{
        CentroidFunctor()
        {
                mSum = 0;
                mWeight = Vector3D::Zero();
                mThreshold = 2;
        }

        void operator()(int x, int y, int z, unsigned char* val)
        {
                if (*val < mThreshold)
                        return;
                mSum += *val;
                mWeight += Vector3D(x,y,z) * (*val);
        }

        Vector3D getCentroid()
        {
                return mWeight/mSum;
        }

private:
        double mSum;
        Vector3D mWeight;
        unsigned char mThreshold;
};

template<class FUNCTOR>
void applyFunctor(cx::ImagePtr image, FUNCTOR& func)
{
        vtkImageDataPtr raw = image->getBaseVtkImageData();
        Eigen::Array3i dims = Eigen::Array3i(raw->GetDimensions());
        unsigned char *pixels = (unsigned char*)raw->GetScalarPointer();

        for(int z = 0; z < dims[2]; ++z)
        {
                for(int y = 0; y < dims[1]; ++y)
                {
                        for(int x = 0; x < dims[0]; ++x)
                        {
                                int index = x + y*dims[0] + z*dims[1]*dims[0];
                                func(x,y,z, pixels+index);
                        }
                }
        }
}

cx::Vector3D calculateCentroid(cx::ImagePtr image)
{
        vtkImageDataPtr raw = image->getBaseVtkImageData();
        cx::Transform3D rMd = image->get_rMd();
        Eigen::Array3d spacing = Eigen::Array3d(raw->GetSpacing());

        CentroidFunctor func;
        applyFunctor(image, func);
        Vector3D ci = func.getCentroid().array() * spacing;
        return rMd.coord(ci);
}

double calculateMass(cx::ImagePtr image)
{
        MassFunctor func;
        applyFunctor(image, func);
        return func.getVolume();
}

}