Notice to GitKraken users: A vulnerability has been found in the SSH key generation of GitKraken versions 7.6.0 to 8.0.0 (https://www.gitkraken.com/blog/weak-ssh-key-fix). If you use GitKraken and have generated a SSH key using one of these versions, please remove it both from your local workstation and from your LRZ GitLab profile.

21.10.2021, 9:00 - 11:00: Due to updates GitLab may be unavailable for some minutes between 09:00 and 11:00.

Commit 7bd28140 authored by Christian Schulte zu Berge's avatar Christian Schulte zu Berge
Browse files

added masking and clipping to SimilarityMeasure

parent 139f8257
......@@ -41,20 +41,33 @@ uniform TextureParameters3D _movingTextureParams;
uniform float _zTex;
uniform mat4 _registrationInverse;
uniform bool _applyMask;
uniform vec2 _xClampRange;
uniform vec2 _yClampRange;
void main() {
// compute lookup coordinates
vec3 referenceLookupTexCoord = vec3(ex_TexCoord.xy, _zTex);
vec4 movingLookupTexCoord = _registrationInverse * vec4(referenceLookupTexCoord, 1.0);
float sad = 0.0;
// fetch texels
float referenceValue = texture(_referenceTexture, referenceLookupTexCoord).a;
float movingValue = texture(_movingTexture, movingLookupTexCoord.xyz).a;
if ( ex_TexCoord.x >= _xClampRange.x && ex_TexCoord.x <= _xClampRange.y
&& ex_TexCoord.y >= _yClampRange.x && ex_TexCoord.y <= _yClampRange.y
&& _zTex >= _zClampRange.x && _zTex <= _zClampRange.y) {
// compute lookup coordinates
vec3 referenceLookupTexCoord = vec3(ex_TexCoord.xy, _zTex);
vec4 movingLookupTexCoord = _registrationInverse * vec4(referenceLookupTexCoord, 1.0);
// compute differences
float difference = referenceValue - movingValue;
float sad = abs(difference);
//float ssd = difference * difference;
// fetch texels
float referenceValue = texture(_referenceTexture, referenceLookupTexCoord).a;
float movingValue = 0.0;
if (!_applyMask || referenceValue > 0.0)
movingValue = texture(_movingTexture, movingLookupTexCoord.xyz).a;
// compute differences
float difference = referenceValue - movingValue;
sad = abs(difference);
//float ssd = difference * difference;
}
// write output color
out_Color = vec4(sad, sad, sad, sad);
}
......@@ -39,30 +39,45 @@ uniform TextureParameters3D _referenceTextureParams;
uniform sampler3D _movingTexture;
uniform TextureParameters3D _movingTextureParams;
//uniform mat4 _registrationMatrix;
uniform mat4 _registrationInverse;
uniform bool _applyMask = false;
uniform vec2 _xClampRange = vec2(0.0, 1.0);
uniform vec2 _yClampRange = vec2(0.0, 1.0);
uniform vec2 _zClampRange = vec2(0.0, 1.0);
void main() {
ivec2 texel = ivec2(ex_TexCoord.xy * _referenceTextureParams._size.xy);
float depth = _referenceTextureParams._size.z;
float sad = 0.0;
float ssd = 0.0;
for (float z = _referenceTextureParams._sizeRCP.z / 2.0; z < 1.0; z += _referenceTextureParams._sizeRCP.z) {
vec3 referenceLookupTexCoord = vec3(ex_TexCoord.xy, z);
vec4 movingLookupTexCoord = _registrationInverse * vec4(referenceLookupTexCoord, 1.0);
//movingLookupTexCoord.xyz /= movingLookupTexCoord.z;
if (ex_TexCoord.x >= _xClampRange.x && ex_TexCoord.x <= _xClampRange.y && ex_TexCoord.y >= _yClampRange.x && ex_TexCoord.y <= _yClampRange.y) {
float zStart = min(_referenceTextureParams._sizeRCP.z / 2.0, _zClampRange.x);
float zEnd = min(1.0, _zClampRange.y);
// if (all(greaterThanEqual(movingLookupTexCoord.xyz, vec3(0.0))) && all(lessThanEqual(movingLookupTexCoord.xyz, vec3(1.0)))) {
for (float z = zStart; z < zEnd; z += _referenceTextureParams._sizeRCP.z) {
// fetch value from reference volume
vec3 referenceLookupTexCoord = vec3(ex_TexCoord.xy, z);
float referenceValue = texture(_referenceTexture, referenceLookupTexCoord).a;
float movingValue = texture(_movingTexture, movingLookupTexCoord.xyz).a;
float difference = referenceValue - movingValue;
sad += abs(difference);
ssd += difference * difference;
// }
// apply mask if requested
if (!_applyMask || referenceValue > 0.0) {
// compute moving lookup texture coordinates
vec4 movingLookupTexCoord = _registrationInverse * vec4(referenceLookupTexCoord, 1.0);
//movingLookupTexCoord.xyz /= movingLookupTexCoord.z;
// fetch value from moving volume
float movingValue = 0.0;
if (all(greaterThanEqual(movingLookupTexCoord.xyz, vec3(0.0))) && all(lessThanEqual(movingLookupTexCoord.xyz, vec3(1.0)))) {
movingValue = texture(_movingTexture, movingLookupTexCoord.xyz).a;
}
// compute difference metrics
float difference = referenceValue - movingValue;
sad += abs(difference);
ssd += difference * difference;
}
}
}
//out_Value = toReturn;
out_Color = vec4(ssd, sad, 0.0, 1.0);
}
......@@ -60,15 +60,14 @@ namespace campvis {
void ReducerTest::init() {
AutoEvaluationPipeline::init();
_referenceReader.p_url.setValue("D:/Medical Data/SCR/Data/RegSweeps_phantom_cropped/-1S3/Volume_0.mhd");
_referenceReader.p_url.setValue("D:/Medical Data/SCR/Data/RegSweeps_phantom_cropped/-1S3median/Volume_0.mhd");
_referenceReader.p_targetImageID.setValue("Reference Image");
_referenceReader.p_targetImageID.addSharedProperty(&_sm.p_referenceId);
_movingReader.p_url.setValue("D:/Medical Data/SCR/Data/RegSweeps_phantom_cropped/-1S3/Volume_1.mhd");
_movingReader.p_url.setValue("D:/Medical Data/SCR/Data/RegSweeps_phantom_cropped/-1S3median/Volume_1.mhd");
_movingReader.p_targetImageID.setValue("Moving Image");
_movingReader.p_targetImageID.addSharedProperty(&_sm.p_movingId);
_sm.p_differenceImageId.addSharedProperty(&_ve.p_inputVolume);
_ve.p_outputImage.setValue("renderTarget");
......
......@@ -55,8 +55,12 @@ namespace campvis {
SimilarityMeasure::SimilarityMeasure()
: VisualizationProcessor(0)
, p_referenceId("ReferenceId", "Reference Image", "", DataNameProperty::READ, AbstractProcessor::VALID)
, p_referenceId("ReferenceId", "Reference Image", "", DataNameProperty::READ, AbstractProcessor::INVALID_PROPERTIES)
, p_movingId("MovingId", "Moving Image", "", DataNameProperty::READ, AbstractProcessor::VALID)
, p_clipX("clipX", "X Axis Clip Coordinates", tgt::ivec2(0), tgt::ivec2(0), tgt::ivec2(0))
, p_clipY("clipY", "Y Axis Clip Coordinates", tgt::ivec2(0), tgt::ivec2(0), tgt::ivec2(0))
, p_clipZ("clipZ", "Z Axis Clip Coordinates", tgt::ivec2(0), tgt::ivec2(0), tgt::ivec2(0))
, p_applyMask("ApplyMask", "Apply Mask", true)
, p_translation("Translation", "Moving Image Translation", tgt::vec3(0.f), tgt::vec3(-1000.f), tgt::vec3(1000.f), tgt::vec3(1.f), tgt::vec3(1.f))
, p_rotation("Rotation", "Moving Image Rotation", tgt::vec3(0.f), tgt::vec3(-tgt::PIf), tgt::vec3(tgt::PIf), tgt::vec3(.1f), tgt::vec3(2.f))
, p_viewportSize("ViewportSize", "Viewport Size", tgt::ivec2(1), tgt::ivec2(1), tgt::ivec2(1000), tgt::ivec2(1), AbstractProcessor::VALID)
......@@ -65,18 +69,31 @@ namespace campvis {
, p_performOptimization("PerformOptimization", "Perform Optimization", AbstractProcessor::INVALID_RESULT | PERFORM_OPTIMIZATION)
, p_differenceImageId("DifferenceImageId", "Difference Image", "difference", DataNameProperty::WRITE, AbstractProcessor::VALID)
, p_computeDifferenceImage("ComputeDifferenceImage", "Compute Difference Image", AbstractProcessor::INVALID_RESULT | COMPUTE_DIFFERENCE_IMAGE)
, p_forceStop("Force Stop", "Force Stop", AbstractProcessor::VALID)
, _costFunctionShader(0)
, _glr(0)
, _opt(0)
{
addProperty(&p_referenceId);
addProperty(&p_movingId);
addProperty(&p_clipX);
addProperty(&p_clipY);
addProperty(&p_clipZ);
addProperty(&p_applyMask);
addProperty(&p_translation);
addProperty(&p_rotation);
addProperty(&p_computeSimilarity);
addProperty(&p_differenceImageId);
addProperty(&p_computeDifferenceImage);
addProperty(&p_optimizer);
addProperty(&p_performOptimization);
addProperty(&p_forceStop);
p_forceStop.s_clicked.connect(this, &SimilarityMeasure::forceStop);
_viewportSizeProperty = &p_viewportSize;
}
......@@ -104,6 +121,9 @@ namespace campvis {
delete _glr;
_glr = 0;
delete _opt;
_opt = 0;
}
void SimilarityMeasure::process(DataContainer& data) {
......@@ -132,6 +152,14 @@ namespace campvis {
ScopedTypedData<ImageData> referenceImage(dc, p_referenceId.getValue());
if (referenceImage != 0) {
p_viewportSize.setValue(referenceImage->getSize().xy());
p_clipX.setMaxValue(tgt::ivec2(static_cast<int>(referenceImage->getSize().x), static_cast<int>(referenceImage->getSize().x)));
p_clipY.setMaxValue(tgt::ivec2(static_cast<int>(referenceImage->getSize().y), static_cast<int>(referenceImage->getSize().y)));
p_clipZ.setMaxValue(tgt::ivec2(static_cast<int>(referenceImage->getSize().z), static_cast<int>(referenceImage->getSize().z)));
p_clipX.setValue(tgt::ivec2(0, static_cast<int>(referenceImage->getSize().x)));
p_clipY.setValue(tgt::ivec2(0, static_cast<int>(referenceImage->getSize().y)));
p_clipZ.setValue(tgt::ivec2(0, static_cast<int>(referenceImage->getSize().z)));
}
validate(AbstractProcessor::INVALID_PROPERTIES);
......@@ -143,9 +171,9 @@ namespace campvis {
MyFuncData_t mfd = { this, referenceImage, movingImage, 0 };
nlopt::opt opt(p_optimizer.getOptionValue(), 6);
opt.set_min_objective(&SimilarityMeasure::optimizerFunc, &mfd);
opt.set_xtol_rel(1e-4);
_opt = new nlopt::opt(p_optimizer.getOptionValue(), 6);
_opt->set_min_objective(&SimilarityMeasure::optimizerFunc, &mfd);
_opt->set_xtol_rel(1e-4);
std::vector<double> x(6);
x[0] = p_translation.getValue().x;
......@@ -162,17 +190,27 @@ namespace campvis {
stepSize[3] = 0.5;
stepSize[4] = 0.5;
stepSize[5] = 0.5;
opt.set_initial_step(stepSize);
_opt->set_initial_step(stepSize);
double minf;
nlopt::result result = opt.optimize(x, minf);
if (result >= nlopt::SUCCESS) {
nlopt::result result = nlopt::SUCCESS;
try {
result = _opt->optimize(x, minf);
}
catch (std::exception& e) {
LERROR("Excpetion during optimization: " << e.what());
}
if (result >= nlopt::SUCCESS || result <= nlopt::ROUNDOFF_LIMITED) {
LDEBUG("Optimization successful, took " << mfd._count << " steps.");
p_translation.setValue(tgt::vec3(x[0], x[1], x[2]));
p_rotation.setValue(tgt::vec3(x[3], x[4], x[5]));
}
delete _opt;
_opt = 0;
validate(PERFORM_OPTIMIZATION);
}
......@@ -202,18 +240,22 @@ namespace campvis {
// bind input images
_costFunctionShader->activate();
_costFunctionShader->setUniform("_applyMask", p_applyMask.getValue());
_costFunctionShader->setUniform("_xClampRange", tgt::vec2(p_clipX.getValue()) / static_cast<float>(size.x));
_costFunctionShader->setUniform("_yClampRange", tgt::vec2(p_clipY.getValue()) / static_cast<float>(size.y));
_costFunctionShader->setUniform("_zClampRange", tgt::vec2(p_clipZ.getValue()) / static_cast<float>(size.z));
referenceImage->bind(_costFunctionShader, referenceUnit, "_referenceTexture", "_referenceTextureParams");
movingImage->bind(_costFunctionShader, movingUnit, "_movingTexture", "_movingTextureParams");
tgt::mat4 registrationMatrix = tgt::mat4::createTranslation(translation) * euleranglesToMat4(rotation);
tgt::mat4 registrationInverse;
if (! registrationMatrix.invert(registrationInverse))
tgtAssert(false, "Could not invert registration matrix. This should not happen!");
const tgt::mat4& w2t = movingImage->getParent()->getMappingInformation().getWorldToTextureMatrix();
const tgt::mat4& t2w = referenceImage->getParent()->getMappingInformation().getTextureToWorldMatrix();
registrationMatrix = w2t * registrationMatrix * t2w;
registrationInverse = w2t * registrationInverse * t2w;
tgt::mat4 registrationInverse;
if (! registrationMatrix.invert(registrationInverse))
tgtAssert(false, "Could not invert registration matrix. This should not happen!");
// render quad to compute similarity measure by shader
_costFunctionShader->setUniform("_registrationInverse", registrationInverse);
......@@ -253,16 +295,16 @@ namespace campvis {
}
tgt::mat4 SimilarityMeasure::euleranglesToMat4(const tgt::vec3& eulerAngles) {
float sinX = sin(eulerAngles.x);
float cosX = cos(eulerAngles.x);
float sinY = sin(eulerAngles.y);
float cosY = cos(eulerAngles.y);
float sinZ = sin(eulerAngles.z);
float cosZ = cos(eulerAngles.z);
float sinX = sin(eulerAngles.x);
float cosX = cos(eulerAngles.x);
float sinY = sin(eulerAngles.y);
float cosY = cos(eulerAngles.y);
float sinZ = sin(eulerAngles.z);
float cosZ = cos(eulerAngles.z);
tgt::mat4 toReturn(cosY * cosZ, cosZ * sinX * sinY - cosX * sinZ, sinX * sinZ + cosX * cosZ * sinY, 0.f,
cosY * sinZ, sinX * sinY * sinZ + cosX * cosZ, cosX * sinY * sinZ - cosZ * sinX, 0.f,
(-1) * sinY, cosY * sinX, cosX * cosY, 0.f,
cosY * sinZ, sinX * sinY * sinZ + cosX * cosZ, cosX * sinY * sinZ - cosZ * sinX, 0.f,
(-1) * sinY, cosY * sinX, cosX * cosY, 0.f,
0.f, 0.f, 0.f, 1.f);
return toReturn;
}
......@@ -291,18 +333,22 @@ namespace campvis {
// bind input images
_differenceShader->activate();
_differenceShader->setUniform("_applyMask", p_applyMask.getValue());
_differenceShader->setUniform("_xClampRange", tgt::vec2(p_clipX.getValue()) / static_cast<float>(size.x));
_differenceShader->setUniform("_yClampRange", tgt::vec2(p_clipY.getValue()) / static_cast<float>(size.y));
_differenceShader->setUniform("_zClampRange", tgt::vec2(p_clipZ.getValue()) / static_cast<float>(size.z));
referenceImage->bind(_differenceShader, referenceUnit, "_referenceTexture", "_referenceTextureParams");
movingImage->bind(_differenceShader, movingUnit, "_movingTexture", "_movingTextureParams");
tgt::mat4 registrationMatrix = tgt::mat4::createTranslation(translation) * euleranglesToMat4(rotation);
tgt::mat4 registrationInverse;
if (! registrationMatrix.invert(registrationInverse))
tgtAssert(false, "Could not invert registration matrix. This should not happen!");
const tgt::mat4& w2t = movingImage->getParent()->getMappingInformation().getWorldToTextureMatrix();
const tgt::mat4& t2w = referenceImage->getParent()->getMappingInformation().getTextureToWorldMatrix();
registrationMatrix = w2t * registrationMatrix * t2w;
registrationInverse = w2t * registrationInverse * t2w;
tgt::mat4 registrationInverse;
if (! registrationMatrix.invert(registrationInverse))
tgtAssert(false, "Could not invert registration matrix. This should not happen!");
// render quad to compute similarity measure by shader
_differenceShader->setUniform("_registrationInverse", registrationInverse);
......@@ -328,5 +374,12 @@ namespace campvis {
validate(COMPUTE_DIFFERENCE_IMAGE);
}
void SimilarityMeasure::forceStop() {
if (_opt != 0)
_opt->force_stop();
validate(PERFORM_OPTIMIZATION);
}
}
......@@ -98,6 +98,11 @@ namespace campvis {
DataNameProperty p_referenceId; ///< image ID for reference image
DataNameProperty p_movingId; ///< image ID for moving image
IVec2Property p_clipX; ///< clip coordinates for x axis
IVec2Property p_clipY; ///< clip coordinates for y axis
IVec2Property p_clipZ; ///< clip coordinates for z axis
BoolProperty p_applyMask;
Vec3Property p_translation; ///< Moving image translation
Vec3Property p_rotation; ///< Moving image rotation
......@@ -108,6 +113,7 @@ namespace campvis {
GenericOptionProperty<nlopt::algorithm> p_optimizer;
ButtonProperty p_performOptimization;
ButtonProperty p_forceStop;
private:
struct MyFuncData_t {
......@@ -122,6 +128,8 @@ namespace campvis {
void performOptimization(const ImageRepresentationGL* referenceImage, const ImageRepresentationGL* movingImage);
void forceStop();
float computeSimilarity(const ImageRepresentationGL* referenceImage, const ImageRepresentationGL* movingImage, const tgt::vec3& translation, const tgt::vec3& rotation);
void generateDifferenceImage(DataContainer* dc, const ImageRepresentationGL* referenceImage, const ImageRepresentationGL* movingImage, const tgt::vec3& translation, const tgt::vec3& rotation);
......@@ -135,6 +143,7 @@ namespace campvis {
tgt::Shader* _costFunctionShader; ///< Shader for slice rendering
tgt::Shader* _differenceShader; ///< Shader for computing the difference image
GlReduction* _glr;
nlopt::opt* _opt;
static const std::string loggerCat_;
};
......
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