cxSliceComputer.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 "cxSliceComputer.h"
#include "cxDefinitions.h"
#include <math.h>

namespace cx
{

SlicePlane::SlicePlane(const Vector3D& c_, const Vector3D& i_, const Vector3D& j_) :
        c(c_), i(i_), j(j_)
{
}

std::ostream& operator<<(std::ostream& s, const SlicePlane& val)
{
        s << "center   : " << val.c << std::endl;
        s << "i_vector : " << val.i << std::endl;
        s << "j_vector : " << val.j << std::endl;
        return s;
}

bool similar(const SlicePlane& a, const SlicePlane& b)
{
        return similar(a.c, b.c) && similar(a.i, b.i) && similar(a.j, b.j);
}

SliceComputer::SliceComputer() :
        mClinicalApplication(mdNEUROLOGICAL),
        mOrientType(otORTHOGONAL),
        mPlaneType(ptAXIAL),
        mFollowType(ftFIXED_CENTER),
        mFixedCenter(Vector3D(0,0,0)),
        m_rMt(Transform3D::Identity()),
        mToolOffset(0.0),
        mUseGravity(false),
        mGravityDirection(Vector3D(0,0,-1)) ,
        mUseViewOffset(false),
        mViewportHeight(1),
        mViewOffset(0.5)
{
}

SliceComputer::~SliceComputer()
{
}

void SliceComputer::initializeFromPlane(PLANE_TYPE plane, bool useGravity, const Vector3D& gravityDir, bool useViewOffset, double viewportHeight, double toolViewOffset, CLINICAL_VIEW application)
{
        setPlaneType(plane);
        mClinicalApplication = application;

        if (plane == ptSAGITTAL || plane == ptCORONAL || plane == ptAXIAL )
        {
                setOrientationType(otORTHOGONAL);
                setFollowType(ftFIXED_CENTER);
        }
    else if (plane == ptANYPLANE || plane==ptRADIALPLANE || plane==ptSIDEPLANE)
        {
                setOrientationType(otOBLIQUE);
                setFollowType(ftFOLLOW_TOOL);

                setGravity(useGravity, gravityDir);
                setToolViewOffset(useViewOffset, viewportHeight, toolViewOffset); // TODO finish this one
        }
    else if (plane==ptTOOLSIDEPLANE)
    {
        setOrientationType(otOBLIQUE);
        setFollowType(ftFIXED_CENTER);
        setGravity(useGravity, gravityDir);
    }
}

void SliceComputer::setClinicalApplication(CLINICAL_VIEW application)
{
        mClinicalApplication = application;
}

ORIENTATION_TYPE SliceComputer::getOrientationType() const
{
        return mOrientType;
}

PLANE_TYPE SliceComputer::getPlaneType() const
{
        return mPlaneType;
}

FOLLOW_TYPE SliceComputer::getFollowType() const
{
    return mFollowType;
}

Transform3D SliceComputer::getToolPosition() const
{
        return m_rMt;
}

void SliceComputer::setToolPosition(const Transform3D& rMt) 
{ 
        m_rMt = rMt; 
}

void SliceComputer::setOrientationType(ORIENTATION_TYPE val) 
{ 
        mOrientType = val; 
}

void SliceComputer::setPlaneType(PLANE_TYPE val) 
{ 
        mPlaneType = val; 
}

void SliceComputer::setFixedCenter(const Vector3D& center) 
{ 
        mFixedCenter = center; 
}

void SliceComputer::setFollowType(FOLLOW_TYPE val) 
{ 
        mFollowType = val; 
}

void SliceComputer::setGravity(bool use, const Vector3D& dir) 
{ 
        mUseGravity = use; 
        mGravityDirection = dir; 
}

void SliceComputer::setToolOffset(double val) 
{ 
        mToolOffset = val; 
}

void SliceComputer::setToolViewOffset(bool use, double viewportHeight, double viewOffset)
{
        mUseViewOffset = use; 
        mViewportHeight = viewportHeight; 
        mViewOffset = viewOffset;
}

void SliceComputer::setToolViewportHeight(double viewportHeight)
{
        mViewportHeight = viewportHeight; 
}

SlicePlane SliceComputer::getPlane()  const
{
        std::pair<Vector3D,Vector3D> basis = generateBasisVectors();
        SlicePlane plane;
        plane.i = basis.first;
        plane.j = basis.second;
        plane.c = Vector3D(0,0,mToolOffset);

        // transform position from tool to reference space
        plane.c = m_rMt.coord(plane.c);
        // transform orientation from tool to reference space for the oblique case only
        if (mOrientType==otOBLIQUE)
        {
                plane.i = m_rMt.vector(plane.i);
                plane.j = m_rMt.vector(plane.j);
        }

    if (mPlaneType == ptTOOLSIDEPLANE)
    {
        plane = this->orientToGravityAroundToolZAxisAndAlongTheOperatingTable(plane);
    }
    else
    {
        plane = orientToGravity(plane);
    }

        // try to to this also for oblique views, IF the ftFIXED_CENTER is set.
        // use special acs centermod algo
        plane.c = generateFixedIJCenter(mFixedCenter, plane.c, plane.i, plane.j);

        // set center so that it is a fixed distance from viewport top
        plane = applyViewOffset(plane);

        return plane; 
}

SlicePlane SliceComputer::applyViewOffset(const SlicePlane& base) const
{
        if (!mUseViewOffset)
        {
                return base;
        }

        double centerOffset = this->getViewOffsetAbsoluteFromCenter();

        SlicePlane retval = base;

        // limit by tooloffset
        Vector3D toolOffsetCenter = m_rMt.coord(Vector3D(0,0,mToolOffset));
        Vector3D newCenter = toolOffsetCenter + centerOffset * base.j;
        double toolOffsetDistance = dot(newCenter - base.c, base.j);

        // limit by tooltip
        Vector3D toolCenter = m_rMt.coord(Vector3D(0,0,0));
        newCenter = toolCenter - centerOffset * base.j;
        double toolDistance = dot(newCenter - base.c, base.j);

        // select a dist and apply
        double usedDistance = std::min(toolOffsetDistance, toolDistance);
        retval.c = base.c + usedDistance * base.j; // extract j-component of newCenter

        return retval;
}

double SliceComputer::getViewOffsetAbsoluteFromCenter() const
{
        if (mPlaneType==ptRADIALPLANE)
                return 0; // position in the center

        return mViewportHeight*(0.5-mViewOffset);
}

std::pair<Vector3D,Vector3D> SliceComputer::generateBasisVectors() const
{
        switch (mClinicalApplication)
        {
        case mdRADIOLOGICAL:
                return this->generateBasisVectorsRadiology();
        case mdNEUROLOGICAL:
        default:
                return this->generateBasisVectorsNeurology();
        }
}

std::pair<Vector3D,Vector3D> SliceComputer::generateBasisVectorsNeurology() const
{
        switch (mPlaneType)
        {
        // use left-right ordering for axial/coronal
        case ptAXIAL:       return std::make_pair(Vector3D(-1, 0, 0), Vector3D( 0,-1, 0));
        case ptCORONAL:     return std::make_pair(Vector3D(-1, 0, 0), Vector3D( 0, 0, 1));
        case ptSAGITTAL:    return std::make_pair(Vector3D( 0, 1, 0), Vector3D( 0, 0, 1));

        // use planes corresponding to the cx Tool definitions
        case ptANYPLANE:    return std::make_pair(Vector3D( 0,-1, 0), Vector3D( 0, 0,-1));
        case ptSIDEPLANE:   return std::make_pair(Vector3D(-1, 0, 0), Vector3D( 0, 0,-1));
        case ptRADIALPLANE: return std::make_pair(Vector3D( 0,-1, 0), Vector3D(-1, 0, 0));
    case ptTOOLSIDEPLANE: return std::make_pair(Vector3D(-1, 0, 0), Vector3D( 0, 0,-1)); //SIDE
        default:
                throw std::exception();
        }
}

std::pair<Vector3D,Vector3D> SliceComputer::generateBasisVectorsRadiology() const
{
        switch (mPlaneType)
        {
        // use right-left ordering for axial/coronal
        case ptAXIAL:       return std::make_pair(Vector3D( 1, 0, 0), Vector3D( 0,-1, 0));
        case ptCORONAL:     return std::make_pair(Vector3D( 1, 0, 0), Vector3D( 0, 0, 1));
        case ptSAGITTAL:    return std::make_pair(Vector3D( 0, 1, 0), Vector3D( 0, 0, 1));

        // use planes corresponding to the cx Tool definitions
        case ptANYPLANE:    return std::make_pair(Vector3D( 0,-1, 0), Vector3D( 0, 0,-1));
        case ptSIDEPLANE:   return std::make_pair(Vector3D(-1, 0, 0), Vector3D( 0, 0,-1));
        case ptRADIALPLANE: return std::make_pair(Vector3D( 0,-1, 0), Vector3D(-1, 0, 0));
    case ptTOOLSIDEPLANE: return std::make_pair(Vector3D(-1, 0, 0), Vector3D( 0, 0,-1)); //SIDE
        default:
                throw std::exception();
        }
}

Vector3D SliceComputer::generateFixedIJCenter(const Vector3D& center_r, const Vector3D& cross_r, const Vector3D& i, const Vector3D& j) const
{
        if (mFollowType==ftFIXED_CENTER)
        {
                // r is REF, s is SLICE
                Transform3D M_rs = createTransformIJC(i, j, Vector3D(0,0,0)); // transform from data to slice, zero center.
                Transform3D M_sr = M_rs.inv();
                Vector3D center_s = M_sr.coord(center_r);
                Vector3D cross_s = M_sr.coord(cross_r);
                // in SLICE space, use {xy} values from center and {z} value from cross.
                Vector3D q_s(center_s[0], center_s[1], cross_s[2]);
                Vector3D q_r = M_rs.coord(q_s);
                return q_r;
        }
        return cross_r;
}

SlicePlane SliceComputer::orientToGravity(const SlicePlane& base) const
{
        if (!mUseGravity)
        {
                return base;
        }

        SlicePlane retval = base;
        const Vector3D k = cross(base.i, base.j); // plane normal. Constant
        Vector3D up;
        up = -mGravityDirection; // normal case

        // weight of nongravity, 0=<w=<1, 1 means dont use gravity
        double w_n = dot(up, k); 
        w_n = w_n*w_n; // square to keep stability near normal use.

        Vector3D i_g = cross(up, k); //  |i_g| varies from 0 to 1 depending on 1-w_n
    Vector3D i_n = base.i; // |i_n|==1 //It seems to me that i_n might need a change to e.g. i_n = -base.j, to be good in the singularity situation. Look into that if using this method later. jone, 20160712

        // set i vector to a weighted mean of the two definitions
        // can also experiment with a tanh function or simply a linear interpolation
        //
        // Note: i_g is already small here if w_n is small, this will increase
        // the effect of i_n. Investigate.
        //
        retval.i = i_g*(1.0-w_n) + i_n*w_n; 
        retval.i = retval.i.normal(); // |i|==1 
        retval.j = cross(k, retval.i);

    return retval;
}

SlicePlane SliceComputer::orientToGravityAroundToolZAxisAndAlongTheOperatingTable(const SlicePlane &base) const
{
    if (!mUseGravity)
    {
        return base;
    }

    SlicePlane retval = base;
    Vector3D up = -mGravityDirection;

    Vector3D k_neg = -cross(base.i, base.j);
    Vector3D toolVector = m_rMt.vector(Vector3D(0, 0, 1));
    Vector3D i_perpendicular = cross(up, toolVector).normal();

        double w_n = this->getWeightForAngularDifference(up, toolVector);

        i_perpendicular = i_perpendicular*(1.0-w_n) + k_neg*w_n;

    Vector3D i_mark = cross(i_perpendicular, up).normal();
    Vector3D j_mark = up;

    retval.i = i_mark;
    retval.j = j_mark;

    return retval;
}

double SliceComputer::getWeightForAngularDifference(Vector3D a, Vector3D b) const
{
        double w_n = dot(a, b);
        w_n = fabs(w_n);
        // w_n = 0 : normal case
        // w_n = 1 : singularity

        double cutoff = sqrt(3.0)/2.0; // cutoff at 30*, i.e. use only toolvector up to that angle between up and tool
        if (w_n<cutoff)
                w_n = 0;
        else
                w_n = (w_n-cutoff)/(1.0-cutoff);
        return w_n;
}


} // namespace cx