Commit 7aeb1912 authored by Emmanuel Promayon's avatar Emmanuel Promayon

NEW first version of translation value computation

parent 38e44f64
......@@ -40,7 +40,9 @@
#include <vtkProperty.h>
#include <vtkMatrix4x4.h>
// Maths
#include <cmath>
#include <QVector3D>
namespace camitk {
......@@ -71,16 +73,29 @@ namespace camitk {
template<typename T>
vtkSmartPointer<vtkMatrix4x4> Multiply4x4(T a, T b) {
vtkSmartPointer<vtkMatrix4x4> c = vtkSmartPointer<vtkMatrix4x4>::New();
vtkMatrix4x4::Multiply4x4(a, b, c);
return c;
vtkSmartPointer<vtkMatrix4x4> c = vtkSmartPointer<vtkMatrix4x4>::New();
vtkMatrix4x4::Multiply4x4(a, b, c);
return c;
}
template<typename T, typename... Args>
vtkSmartPointer<vtkMatrix4x4> Multiply4x4(T a, T b, Args... args) {
return Multiply4x4(a, Multiply4x4(b, args...));
return Multiply4x4(a, Multiply4x4(b, args...));
}
QVector3D linePlaneIntersectionPoint(QVector3D lineVector, QVector3D linePoint, QVector3D planeNormal, QVector3D planePoint) {
double lineDPPlane = QVector3D::dotProduct(lineVector, planeNormal);
if (lineDPPlane < 1e-10) {
// line and plane are parallel (hopefully should never happen!)
return planePoint;
}
else {
QVector3D planeToPoint = linePoint - planePoint;
double planeToPointDPPlane = QVector3D::dotProduct(planeToPoint, planeNormal);
double k = lineDPPlane / planeToPointDPPlane;
return linePoint - k * lineVector;
}
}
// -------------------- constructor --------------------
ArbitrarySingleImageComponent::ArbitrarySingleImageComponent(Component* parentComponent, const QString& name, vtkSmartPointer<vtkWindowLevelLookupTable> lut)
......@@ -167,15 +182,11 @@ void ArbitrarySingleImageComponent::setTransformRotation(double angleX, double a
R_P->DeepCopy(transfo_R_P->GetMatrix());
// Cumulative rotation = rotate(...)
// vtkMatrix4x4::Multiply4x4(T_P2C, R_P, intermediateRotation1);
// displayMatrix4x4(intermediateRotation1)
// vtkMatrix4x4::Multiply4x4(intermediateRotation1, T_C2P, intermediateRotation2);
// displayMatrix4x4(intermediateRotation2)
// vtkMatrix4x4::Multiply4x4(intermediateRotation2, T_P2L, checkRotation);
// vtkSmartPointer<vtkMatrix4x4> checkRotation = Multiply4x4(T_P2C, R_P, T_C2P, T_P2L);
//checkRotation = by(by(by(by(T_P2C, R_P), T_P2L_rotation_only_inverse), T_C2P), T_P2L);
vtkSmartPointer<vtkMatrix4x4> checkRotation = Multiply4x4(T_P2C, R_P, T_P2L_rotation_only_inverse, T_C2P, T_P2L);
if (checkCenter(checkRotation)) {
CAMITK_INFO_ALT(QString("Ok to rotate of (%1,%2,%3)").arg(angleX).arg(angleY).arg(angleZ))
getTransform()->GetMatrix()->DeepCopy(checkRotation);
......@@ -254,10 +265,85 @@ void ArbitrarySingleImageComponent::setTransformTranslation(double x, double y,
}
}
// -------------------- updateTranslationExtremity --------------------
void ArbitrarySingleImageComponent::updateTranslationExtremity() {
double C_L[4] = {dimensions[0]* spacing[0] / 2.0 - spacing[0] / 2.0,
dimensions[1]* spacing[1] / 2.0 - spacing[1] / 2.0,
0.0,
1.0
};
vtkSmartPointer<vtkMatrix4x4> T_P2L = vtkSmartPointer<vtkMatrix4x4>::New();
T_P2L->DeepCopy(getTransform()->GetMatrix());
double C_P[4];
T_P2L->MultiplyPoint(C_L, C_P);
// z vector in local coordinate system
double V_L[4] = { 0.0, 0.0, 1.0, 1.0};
// z vector in global coordinate system
double V_P[4];
T_P2L->MultiplyPoint(V_L, V_P);
// compute intersection of line C_P + k V_P with the image
QVector3D lineVector(V_P[0], V_P[1], V_P[2]);
QVector3D startPoint(C_P[0], C_P[1], C_P[2]);
QVector3D topCorner(0.0, 0.0, 0.0);
QVector3D bottomCorner(dimensions[0]*spacing[0], dimensions[1]*spacing[1], dimensions[2]*spacing[2]);
// check intersection to front/back plane
QVector3D intersection = linePlaneIntersectionPoint(lineVector, startPoint, QVector3D(0.0, 0.0, 1.0), topCorner);
if (intersection != topCorner && pointInsideVolume(intersection)) {
cPlus = intersection;
cMinus = linePlaneIntersectionPoint(lineVector, startPoint, QVector3D(0.0, 0.0, -1.0), bottomCorner);
}
else {
// check intersection to top/bottom plane
intersection = linePlaneIntersectionPoint(lineVector, startPoint, QVector3D(0.0, 1.0, 0.0), topCorner);
if (intersection != topCorner && pointInsideVolume(intersection)) {
cPlus = intersection;
cMinus = linePlaneIntersectionPoint(lineVector, startPoint, QVector3D(0.0, -1.0, 0.0), bottomCorner);
}
else {
// intersection is with left/right plane
cPlus = linePlaneIntersectionPoint(lineVector, startPoint, QVector3D(1.0, 0.0, 0.0), topCorner);
cMinus = linePlaneIntersectionPoint(lineVector, startPoint, QVector3D(-1.0, 0.0, 0.0), bottomCorner);
}
}
// TODO reorder cPlus/cMinus if needed
// C+ = C + k z
// if k < 0 inverse C+ and C-
CAMITK_INFO_ALT(QString("C- = (%1,%2,%3) C+ = (%4,%5,%6)").arg(cMinus.x()).arg(cMinus.y()).arg(cMinus.z()).arg(cPlus.x()).arg(cPlus.y()).arg(cPlus.z()))
}
// -------------------- getTranslationInVolume --------------------
double ArbitrarySingleImageComponent::getTranslationInVolume() {
// TODO
return 0.0;
double ArbitrarySingleImageComponent::getTranslationInVolume() {
updateTranslationExtremity();
// compute center in P reference frame
double C_L[4] = {dimensions[0]* spacing[0] / 2.0 - spacing[0] / 2.0,
dimensions[1]* spacing[1] / 2.0 - spacing[1] / 2.0,
0.0,
1.0
};
vtkSmartPointer<vtkMatrix4x4> T_P2L = vtkSmartPointer<vtkMatrix4x4>::New();
T_P2L->DeepCopy(getTransform()->GetMatrix());
double C_P[4];
T_P2L->MultiplyPoint(C_L, C_P);
QVector3D C(C_P[0], C_P[1], C_P[2]);
// C- = 0.0
// C+ = 1.0
// C-C = k C-C+
// k = C-C / C-C+
double k = QVector3D((C - cMinus) / (cPlus - cMinus)).length();
CAMITK_INFO_ALT(QString("k=%1").arg(k))
return k;
}
// -------------------- checkCenter --------------------
......@@ -268,18 +354,33 @@ bool ArbitrarySingleImageComponent::checkCenter(vtkSmartPointer<vtkMatrix4x4> tr
double centerInImage[4];
transform->MultiplyPoint(centerInFrame, centerInImage);
bool inside = (centerInImage[0] >= 0 && centerInImage[0] < dimensions[0] * spacing[0] &&
centerInImage[1] >= 0 && centerInImage[1] < dimensions[1] * spacing[1] &&
centerInImage[2] >= 0 && centerInImage[2] < dimensions[2] * spacing[2]);
bool inside = pointInsideVolume(QVector3D(centerInImage[0], centerInImage[1], centerInImage[2]));
CAMITK_INFO(QString("inside=%7 p=(%1,%2,%3) pt=(%4,%5,%6) dim=(%8,%9,%10)").arg(centerInFrame[0]).arg(centerInFrame[1]).arg(centerInFrame[2]).arg(centerInImage[0]).arg(centerInImage[1]).arg(centerInImage[2])
.arg(inside).arg(dimensions[0]*spacing[0]).arg(dimensions[1]*spacing[1]).arg(dimensions[2]*spacing[2])
)
CAMITK_INFO(QString("inside=%7 p=(%1,%2,%3) pt=(%4,%5,%6) dim=(%8,%9,%10)")
.arg(centerInFrame[0]).arg(centerInFrame[1]).arg(centerInFrame[2])
.arg(centerInImage[0]).arg(centerInImage[1]).arg(centerInImage[2])
.arg(inside)
.arg(dimensions[0]*spacing[0]).arg(dimensions[1]*spacing[1]).arg(dimensions[2]*spacing[2]))
// positive check only is displacement will keep the center inside the box
return inside;
}
// -------------------- pointInsideVolume --------------------
bool ArbitrarySingleImageComponent::pointInsideVolume(QVector3D p) {
return (p.x() >= 0 && p.x() <= dimensions[0] * spacing[0]
&& p.y() >= 0 && p.y() <= dimensions[1] * spacing[1]
&& p.z() >= 0 && p.z() <= dimensions[2] * spacing[2]);
}
// -------------------- rotate --------------------
void ArbitrarySingleImageComponent::rotate(double aroundX, double aroundY, double aroundZ) {
CAMITK_WARNING("Not implemented yet")
......
......@@ -26,21 +26,11 @@
#ifndef ARBITRARYSINGLEIMAGEVOLUMECOMPONENT_H
#define ARBITRARYSINGLEIMAGEVOLUMECOMPONENT_H
// -- Core stuff
#include "SingleImageComponent.h"
#include "Slice.h"
// -- QT stuff classes
class QMenu;
// -- VTK stuff
#include <vtkImageReslice.h>
#include <vtkWindowLevelLookupTable.h>
#include <vtkImageChangeInformation.h>
// -- VTK stuff classes
class vtkImageClip;
#include <QVector3D>
namespace camitk {
/**
......@@ -133,11 +123,30 @@ private:
/// check if the center of the frame transformed usin the given matrix stays inside the initial image volume
bool checkCenter(vtkSmartPointer<vtkMatrix4x4>);
/// return true only in point is inside the image volume
bool pointInsideVolume(QVector3D);
/// update cPlus and cMinus
void updateTranslationExtremity();
/// dimension of the whole image (kept here for simplifying code)
int* dimensions;
/// spacing of the image (kept here for simplifying code)
double *spacing;
/// point of intersection of the line passing at the center of the image
/// in the z+ direction, i.e. (0,0,1) in the local frame.
/// This is the extreme possible point to translate C in the z+ direction
/// (required to compute getTranslationInVolume and set new translation)
QVector3D cPlus;
/// point of intersection of the line passing at the center of the image
/// in the z- direction, i.e. (0,0,-1) in the local frame.
/// This is the extreme possible point to translate C in the z- direction
/// (required to compute getTranslationInVolume and set new translation)
QVector3D cMinus;
};
}
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment