...
 
Commits (15)
......@@ -11,6 +11,7 @@ THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLI
#define PTM_CONSTANTS_H
#include <cmath>
#include <cstdlib>
//------------------------------------
// definitions
......@@ -214,5 +215,14 @@ const double ptm_template_graphene[PTM_NUM_POINTS_GRAPHENE][3] = {
{ -18./11+3*sqrt(3)/11, 0, 0 },
};
typedef struct
{
int num;
int8_t correspondences[PTM_MAX_INPUT_POINTS];
size_t atom_indices[PTM_MAX_INPUT_POINTS];
int32_t numbers[PTM_MAX_INPUT_POINTS];
double points[PTM_MAX_INPUT_POINTS][3];
} ptm_atomicenv_t;
#endif
......@@ -25,10 +25,10 @@ extern "C" {
int ptm_index( ptm_local_handle_t local_handle,
size_t atom_index, int (get_neighbours)(void* vdata, size_t _unused_lammps_variable, size_t atom_index, int num, int* ordering, size_t* nbr_indices, int32_t* numbers, double (*nbr_pos)[3]), void* nbrlist,
size_t atom_index, int (get_neighbours)(void* vdata, size_t _unused_lammps_variable, size_t atom_index, int num, ptm_atomicenv_t* env), void* nbrlist,
int32_t flags, bool output_conventional_orientation, //inputs
int32_t* p_type, int32_t* p_alloy_type, double* p_scale, double* p_rmsd, double* q, double* F, double* F_res, double* U, double* P, double* p_interatomic_distance, double* p_lattice_constant,
int* p_best_template_index, const double (**p_best_template)[3], int8_t* output_indices); //outputs
int* p_best_template_index, const double (**p_best_template)[3], ptm_atomicenv_t* output_env); //outputs
int ptm_remap_template( int type, bool output_conventional_orientation, int input_template_index, double* qtarget, double* q,
......
......@@ -179,9 +179,11 @@ int ptm_remap_template( int type, bool output_conventional_orientation, int inpu
int8_t temp[PTM_MAX_POINTS];
memset(temp, -1, PTM_MAX_POINTS * sizeof(int8_t));
int ret = ptm_undo_conventional_orientation(type, input_template_index, q, mapping);
if (ret != 0)
return -1;
if (input_template_index != 0) {
int ret = ptm_undo_conventional_orientation(type, input_template_index, q, mapping);
if (ret != 0)
return -1;
}
int bi = 0;
if (qtarget != NULL)
......@@ -238,14 +240,14 @@ int ptm_remap_template( int type, bool output_conventional_orientation, int inpu
return template_index;
}
static void output_data(ptm::result_t *res, ptm::atomicenv_t* env,
static void output_data(ptm::result_t *res, ptm_atomicenv_t* env,
bool output_conventional_orientation, int32_t *p_type,
int32_t *p_alloy_type, double *p_scale, double *p_rmsd,
double *q, double *F, double *F_res, double *U,
double *P, double *p_interatomic_distance,
double *p_lattice_constant,
int* p_best_template_index, const double (**p_best_template)[3],
int8_t *output_indices) {
ptm_atomicenv_t* output_env) {
const ptm::refdata_t *ref = res->ref_struct;
if (ref == NULL)
return;
......@@ -300,9 +302,14 @@ static void output_data(ptm::result_t *res, ptm::atomicenv_t* env,
ptm::polar_decomposition_3x3(F, false, U, P);
}
if (output_indices != NULL)
for (int i = 0; i < ref->num_nbrs + 1; i++)
output_indices[i] = env->ordering[res->mapping[i]];
if (output_env != NULL) {
for (int i = 0; i < ref->num_nbrs + 1; i++) {
output_env->correspondences[i] = env->correspondences[res->mapping[i]];
output_env->atom_indices[i] = env->atom_indices[res->mapping[i]];
memcpy(output_env->points[i], env->points[res->mapping[i]], 3 * sizeof(double));
}
}
double interatomic_distance = calculate_interatomic_distance(ref->type, res->scale);
double lattice_constant = calculate_lattice_constant(ref->type, interatomic_distance);
......@@ -321,14 +328,13 @@ static void output_data(ptm::result_t *res, ptm::atomicenv_t* env,
extern bool ptm_initialized;
int ptm_index(ptm_local_handle_t local_handle,
size_t atom_index, int (get_neighbours)(void* vdata, size_t _unused_lammps_variable, size_t atom_index, int num, int* ordering, size_t* nbr_indices, int32_t* numbers, double (*nbr_pos)[3]), void* nbrlist,
size_t atom_index, int (get_neighbours)(void* vdata, size_t _unused_lammps_variable, size_t atom_index, int num, ptm_atomicenv_t* env), void* nbrlist,
int32_t flags,
bool output_conventional_orientation, int32_t *p_type,
int32_t *p_alloy_type, double *p_scale, double *p_rmsd, double *q,
double *F, double *F_res, double *U, double *P,
double *p_interatomic_distance, double *p_lattice_constant,
int* p_best_template_index, const double (**p_best_template)[3],
int8_t *output_indices)
int* p_best_template_index, const double (**p_best_template)[3], ptm_atomicenv_t* output_env)
{
int ret = 0;
assert(ptm_initialized);
......@@ -340,15 +346,12 @@ int ptm_index(ptm_local_handle_t local_handle,
res.ref_struct = NULL;
res.rmsd = INFINITY;
if (output_indices != NULL)
memset(output_indices, -1, PTM_MAX_INPUT_POINTS * sizeof(int8_t));
*p_type = PTM_MATCH_NONE;
if (p_alloy_type != NULL)
*p_alloy_type = PTM_ALLOY_NONE;
//------------------------------------------------------------
ptm::atomicenv_t env, dmn_env, grp_env;
ptm_atomicenv_t env, dmn_env, grp_env;
ptm::convexhull_t ch;
double ch_points[PTM_MAX_INPUT_POINTS][3];
......@@ -361,7 +364,7 @@ int ptm_index(ptm_local_handle_t local_handle,
if (flags & PTM_CHECK_BCC)
min_points = PTM_NUM_POINTS_BCC;
int num_points = get_neighbours(nbrlist, -1, atom_index, PTM_MAX_INPUT_POINTS, env.ordering, env.nbr_indices, env.numbers, env.points);
int num_points = get_neighbours(nbrlist, -1, atom_index, PTM_MAX_INPUT_POINTS, &env);
if (num_points < min_points)
return -1;
......@@ -405,7 +408,7 @@ int ptm_index(ptm_local_handle_t local_handle,
if (res.ref_struct == NULL)
return PTM_NO_ERROR;
ptm::atomicenv_t* res_env = &env;
ptm_atomicenv_t* res_env = &env;
if (res.ref_struct->type == PTM_MATCH_DCUB || res.ref_struct->type == PTM_MATCH_DHEX)
res_env = &dmn_env;
else if (res.ref_struct->type == PTM_MATCH_GRAPHENE)
......@@ -413,7 +416,7 @@ int ptm_index(ptm_local_handle_t local_handle,
output_data( &res, res_env, output_conventional_orientation, p_type, p_alloy_type, p_scale,
p_rmsd, q, F, F_res, U, P, p_interatomic_distance,
p_lattice_constant, p_best_template_index, p_best_template, output_indices);
p_lattice_constant, p_best_template_index, p_best_template, output_env);
return PTM_NO_ERROR;
}
......
......@@ -28,10 +28,10 @@ typedef struct
{
int rank;
int inner;
int ordering;
int correspondences;
size_t atom_index;
int32_t number;
double offset[3];
double delta[3];
} atomorder_t;
static bool atomorder_compare(atomorder_t const& a, atomorder_t const& b)
......@@ -42,39 +42,39 @@ static bool atomorder_compare(atomorder_t const& a, atomorder_t const& b)
#define MAX_INNER 4
int calculate_two_shell_neighbour_ordering( int num_inner, int num_outer,
size_t atom_index, int (get_neighbours)(void* vdata, size_t _unused_lammps_variable, size_t atom_index, int num, int* ordering, size_t* nbr_indices, int32_t* numbers, double (*nbr_pos)[3]), void* nbrlist,
ptm::atomicenv_t* output)
size_t atom_index, int (get_neighbours)(void* vdata, size_t _unused_lammps_variable, size_t atom_index, int num, ptm_atomicenv_t* env), void* nbrlist,
ptm_atomicenv_t* output)
{
assert(num_inner <= MAX_INNER);
ptm::atomicenv_t central;
int num_input_points = get_neighbours(nbrlist, -1, atom_index, PTM_MAX_INPUT_POINTS, central.ordering, central.nbr_indices, central.numbers, central.points);
ptm_atomicenv_t central;
int num_input_points = get_neighbours(nbrlist, -1, atom_index, PTM_MAX_INPUT_POINTS, &central);
if (num_input_points < num_inner + 1)
return -1;
std::unordered_set<size_t> claimed;
for (int i=0;i<num_inner+1;i++)
{
output->ordering[i] = central.ordering[i];
output->nbr_indices[i] = central.nbr_indices[i];
output->correspondences[i] = central.correspondences[i];
output->atom_indices[i] = central.atom_indices[i];
output->numbers[i] = central.numbers[i];
memcpy(output->points[i], central.points[i], 3 * sizeof(double));
claimed.insert(central.nbr_indices[i]);
claimed.insert(central.atom_indices[i]);
}
int num_inserted = 0;
atomorder_t data[MAX_INNER * PTM_MAX_INPUT_POINTS];
for (int i=0;i<num_inner;i++)
{
ptm::atomicenv_t inner;
int num_points = get_neighbours(nbrlist, -1, central.nbr_indices[1 + i], PTM_MAX_INPUT_POINTS, inner.ordering, inner.nbr_indices, inner.numbers, inner.points);
ptm_atomicenv_t inner;
int num_points = get_neighbours(nbrlist, -1, central.atom_indices[1 + i], PTM_MAX_INPUT_POINTS, &inner);
if (num_points < num_inner + 1)
return -1;
for (int j=0;j<num_points;j++)
{
size_t key = inner.nbr_indices[j];
size_t key = inner.atom_indices[j];
bool already_claimed = claimed.find(key) != claimed.end();
if (already_claimed)
......@@ -82,13 +82,13 @@ int calculate_two_shell_neighbour_ordering( int num_inner, int num_outer,
data[num_inserted].inner = i;
data[num_inserted].rank = j;
data[num_inserted].ordering = inner.ordering[j];
data[num_inserted].atom_index = inner.nbr_indices[j];
data[num_inserted].correspondences = inner.correspondences[j];
data[num_inserted].atom_index = inner.atom_indices[j];
data[num_inserted].number = inner.numbers[j];
memcpy(data[num_inserted].offset, inner.points[j], 3 * sizeof(double));
memcpy(data[num_inserted].delta, inner.points[j], 3 * sizeof(double));
for (int k=0;k<3;k++)
data[num_inserted].offset[k] += central.points[1 + i][k];
data[num_inserted].delta[k] += central.points[1 + i][k];
num_inserted++;
}
......@@ -107,11 +107,11 @@ int calculate_two_shell_neighbour_ordering( int num_inner, int num_outer,
if (counts[inner] >= num_outer || already_claimed)
continue;
output->ordering[1 + num_inner + num_outer * inner + counts[inner]] = data[i].ordering;
output->correspondences[1 + num_inner + num_outer * inner + counts[inner]] = data[i].correspondences;
output->nbr_indices[1 + num_inner + num_outer * inner + counts[inner]] = nbr_atom_index;
output->atom_indices[1 + num_inner + num_outer * inner + counts[inner]] = nbr_atom_index;
output->numbers[1 + num_inner + num_outer * inner + counts[inner]] = data[i].number;
memcpy(output->points[1 + num_inner + num_outer * inner + counts[inner]], &data[i].offset, 3 * sizeof(double));
memcpy(output->points[1 + num_inner + num_outer * inner + counts[inner]], &data[i].delta, 3 * sizeof(double));
claimed.insert(nbr_atom_index);
counts[inner]++;
......
......@@ -14,18 +14,9 @@ THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLI
namespace ptm {
typedef struct
{
int ordering[PTM_MAX_INPUT_POINTS];
size_t nbr_indices[PTM_MAX_INPUT_POINTS];
int32_t numbers[PTM_MAX_INPUT_POINTS];
double points[PTM_MAX_INPUT_POINTS][3];
} atomicenv_t;
int calculate_two_shell_neighbour_ordering( int num_inner, int num_outer,
size_t atom_index, int (get_neighbours)(void* vdata, size_t _unused_lammps_variable, size_t atom_index, int num, int* ordering, size_t* nbr_indices, int32_t* numbers, double (*nbr_pos)[3]), void* nbrlist,
ptm::atomicenv_t* output);
size_t atom_index, int (get_neighbours)(void* vdata, size_t _unused_lammps_variable, size_t atom_index, int num, ptm_atomicenv_t* env), void* nbrlist,
ptm_atomicenv_t* output);
}
#endif
......
......@@ -390,5 +390,40 @@ double quat_misorientation(double* q1, double* q2)
return acos(quat_quick_misorientation(q1, q2));
}
double quat_quick_disorientation_cubic(double* q0, double* q1)
{
double qrot[4];
double qinv[4] = {q0[0], -q0[1], -q0[2], -q0[3]};
quat_rot(qinv, q1, qrot);
rotate_quaternion_into_cubic_fundamental_zone(qrot);
double t = qrot[0];
t = std::min(1., std::max(-1., t));
return 2 * t * t - 1;
}
double quat_disorientation_cubic(double* q0, double* q1)
{
return acos(quat_quick_disorientation_cubic(q0, q1));
}
double quat_quick_disorientation_hcp_conventional(double* q0, double* q1)
{
double qrot[4];
double qinv[4] = {q0[0], -q0[1], -q0[2], -q0[3]};
quat_rot(qinv, q1, qrot);
rotate_quaternion_into_hcp_conventional_fundamental_zone(qrot);
double t = qrot[0];
t = std::min(1., std::max(-1., t));
return 2 * t * t - 1;
}
double quat_disorientation_hcp_conventional(double* q0, double* q1)
{
return acos(quat_quick_disorientation_hcp_conventional(q0, q1));
}
}
......@@ -170,6 +170,9 @@ int map_quaternion_hcp(double* q, int i);
int map_quaternion_hcp_conventional(double* q, int i);
int map_quaternion_diamond_hexagonal(double* q, int i);
double quat_disorientation_cubic(double* q0, double* q1);
double quat_disorientation_hcp_conventional(double* q0, double* q1);
}
#endif
......
......@@ -192,6 +192,16 @@ public:
/// \return A reference to \c this quaternion, which has been changed.
QuaternionT& operator/=(T s) { x() /= s; y() /= s; z() /= s; w() /= s; return *this; }
/// \brief Component-wise increment operator.
/// \param q The quaternion to add to this quaternion.
/// \return A reference to \c this quaternion, which has been changed.
QuaternionT& operator+=(const QuaternionT& q) { x() += q.x(); y() += q.y(); z() += q.z(); w() += q.w(); return *this; }
/// \brief Component-wise decrement operator.
/// \param q The quaternion to subtract from this quaternion.
/// \return A reference to \c this quaternion, which has been changed.
QuaternionT& operator-=(const QuaternionT& q) { x() -= q.x(); y() -= q.y(); z() -= q.z(); w() -= q.w(); return *this; }
/// \brief Computes the scalar product of two quaternions.
Q_DECL_CONSTEXPR T dot(const QuaternionT& b) const { return x()*b.x() + y()*b.y() + z()*b.z() + w()*b.w(); }
......
......@@ -45,6 +45,8 @@ SET(SourceFiles
modifier/dxa/DislocationTracer.cpp
modifier/elasticstrain/ElasticStrainModifier.cpp
modifier/elasticstrain/ElasticStrainEngine.cpp
modifier/grains/GrainSegmentationModifier.cpp
modifier/grains/GrainSegmentationEngine.cpp
modifier/microstructure/SimplifyMicrostructureModifier.cpp
util/DelaunayTessellation.cpp
)
......
......@@ -32,6 +32,7 @@ OVITO_STANDARD_PLUGIN(CrystalAnalysisGui
modifier/DislocationAnalysisModifierEditor.cpp
modifier/ElasticStrainModifierEditor.cpp
modifier/SimplifyMicrostructureModifierEditor.cpp
modifier/GrainSegmentationModifierEditor.cpp
util/DislocationInspectionApplet.cpp
PLUGIN_DEPENDENCIES CrystalAnalysis ParticlesGui
GUI_PLUGIN
......
///////////////////////////////////////////////////////////////////////////////
//
// Copyright (2018) Alexander Stukowski
//
// This file is part of OVITO (Open Visualization Tool).
//
// OVITO is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// OVITO is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
//
///////////////////////////////////////////////////////////////////////////////
#include <plugins/crystalanalysis/CrystalAnalysis.h>
#include <plugins/crystalanalysis/modifier/grains/GrainSegmentationEngine.h>
#include <plugins/crystalanalysis/modifier/grains/GrainSegmentationModifier.h>
#include <plugins/particles/gui/modifier/analysis/StructureListParameterUI.h>
#include <gui/properties/FloatParameterUI.h>
#include <gui/properties/IntegerParameterUI.h>
#include <gui/properties/BooleanParameterUI.h>
#include <gui/utilities/concurrent/ProgressDialog.h>
#include <core/dataset/DataSetContainer.h>
#include "GrainSegmentationModifierEditor.h"
#include <3rdparty/qwt/qwt_plot_zoneitem.h>
namespace Ovito { namespace Plugins { namespace CrystalAnalysis {
IMPLEMENT_OVITO_CLASS(GrainSegmentationModifierEditor);
SET_OVITO_OBJECT_EDITOR(GrainSegmentationModifier, GrainSegmentationModifierEditor);
/******************************************************************************
* Sets up the UI widgets of the editor.
******************************************************************************/
void GrainSegmentationModifierEditor::createUI(const RolloutInsertionParameters& rolloutParams)
{
// Create the rollout.
QWidget* rollout = createRollout(tr("Grain segmentation"), rolloutParams);
QVBoxLayout* layout = new QVBoxLayout(rollout);
layout->setContentsMargins(4,4,4,4);
layout->setSpacing(6);
QGroupBox* paramsBox = new QGroupBox(tr("Parameters"));
layout->addWidget(paramsBox);
QGridLayout* sublayout2 = new QGridLayout(paramsBox);
sublayout2->setContentsMargins(4,4,4,4);
sublayout2->setSpacing(4);
sublayout2->setColumnStretch(1, 1);
FloatParameterUI* rmsdCutoffUI = new FloatParameterUI(this, PROPERTY_FIELD(GrainSegmentationModifier::rmsdCutoff));
sublayout2->addWidget(rmsdCutoffUI->label(), 0, 0);
sublayout2->addLayout(rmsdCutoffUI->createFieldLayout(), 0, 1);
FloatParameterUI* mergingThresholdUI = new FloatParameterUI(this, PROPERTY_FIELD(GrainSegmentationModifier::mergingThreshold));
sublayout2->addWidget(mergingThresholdUI->label(), 1, 0);
sublayout2->addLayout(mergingThresholdUI->createFieldLayout(), 1, 1);
IntegerParameterUI* minGrainAtomCountUI = new IntegerParameterUI(this, PROPERTY_FIELD(GrainSegmentationModifier::minGrainAtomCount));
sublayout2->addWidget(minGrainAtomCountUI->label(), 2, 0);
sublayout2->addLayout(minGrainAtomCountUI->createFieldLayout(), 2, 1);
QGroupBox* debuggingParamsBox = new QGroupBox(tr("Debugging options"));
layout->addWidget(debuggingParamsBox);
sublayout2 = new QGridLayout(debuggingParamsBox);
sublayout2->setContentsMargins(4,4,4,4);
sublayout2->setSpacing(4);
sublayout2->setColumnStretch(1, 1);
// Orphan atom adoption
BooleanParameterUI* orphanAdoptionUI = new BooleanParameterUI(this, PROPERTY_FIELD(GrainSegmentationModifier::orphanAdoption));
sublayout2->addWidget(orphanAdoptionUI->checkBox(), 0, 0, 1, 2);
BooleanParameterUI* outputBondsUI = new BooleanParameterUI(this, PROPERTY_FIELD(GrainSegmentationModifier::outputBonds));
sublayout2->addWidget(outputBondsUI->checkBox(), 1, 0, 1, 2);
#if 0
QPushButton* grainTrackingButton = new QPushButton(tr("Perform grain tracking..."));
layout->addWidget(grainTrackingButton);
connect(grainTrackingButton, &QPushButton::clicked, this, &GrainSegmentationModifierEditor::onPerformGrainTracking);
#endif
// Status label.
layout->addWidget(statusLabel());
// Structure list.
StructureListParameterUI* structureTypesPUI = new StructureListParameterUI(this, true);
layout->addSpacing(10);
layout->addWidget(new QLabel(tr("Structure types:")));
layout->addWidget(structureTypesPUI->tableWidget());
// Create plot widget for RMSD distribution.
_rmsdPlotWidget = new DataSeriesPlotWidget();
_rmsdPlotWidget->setMinimumHeight(200);
_rmsdPlotWidget->setMaximumHeight(200);
_rmsdRangeIndicator = new QwtPlotZoneItem();
_rmsdRangeIndicator->setOrientation(Qt::Vertical);
_rmsdRangeIndicator->setZ(1);
_rmsdRangeIndicator->attach(_rmsdPlotWidget);
_rmsdRangeIndicator->hide();
layout->addSpacing(10);
layout->addWidget(_rmsdPlotWidget);
connect(this, &GrainSegmentationModifierEditor::contentsReplaced, this, &GrainSegmentationModifierEditor::plotHistogram);
}
/******************************************************************************
* This method is called when a reference target changes.
******************************************************************************/
bool GrainSegmentationModifierEditor::referenceEvent(RefTarget* source, const ReferenceEvent& event)
{
if(source == modifierApplication() && event.type() == ReferenceEvent::PipelineCacheUpdated) {
plotHistogramLater(this);
}
return ModifierPropertiesEditor::referenceEvent(source, event);
}
/******************************************************************************
* Replots the histogram computed by the modifier.
******************************************************************************/
void GrainSegmentationModifierEditor::plotHistogram()
{
GrainSegmentationModifier* modifier = static_object_cast<GrainSegmentationModifier>(editObject());
if(modifier && modifier->rmsdCutoff() > 0) {
_rmsdRangeIndicator->setInterval(0, modifier->rmsdCutoff());
_rmsdRangeIndicator->show();
}
else {
_rmsdRangeIndicator->hide();
}
if(modifierApplication()) {
// Request the modifier's pipeline output.
const PipelineFlowState& state = getModifierOutput();
// Look up the data series in the modifier's pipeline output.
_rmsdPlotWidget->setSeries(state.getObjectBy<DataSeriesObject>(modifierApplication(), QStringLiteral("grains-rmsd")));
}
else {
_rmsdPlotWidget->reset();
}
}
/******************************************************************************
* Is called when the user clicks the 'Perform grain tracking' button.
******************************************************************************/
void GrainSegmentationModifierEditor::onPerformGrainTracking()
{
#if 0
GrainSegmentationModifier* modifier = static_object_cast<GrainSegmentationModifier>(editObject());
GrainSegmentationModifierApplication* modApp = dynamic_object_cast<GrainSegmentationModifierApplication>(someModifierApplication());
if(!modifier || !modApp) return;
undoableTransaction(tr("Grain tracking"), [this,modifier,modApp]() {
ProgressDialog progressDialog(container(), modifier->dataset()->container()->taskManager(), tr("Grain tracking"));
modifier->trackGrains(progressDialog.taskManager(), modApp);
});
#endif
GrainSegmentationModifier* modifier = static_object_cast<GrainSegmentationModifier>(editObject());
if(!modifier) return;
undoableTransaction(tr("Grain tracking"), [this,modifier]() {
ProgressDialog progressDialog(container(), modifier->dataset()->container()->taskManager(), tr("Grain tracking"));
});
}
} // End of namespace
} // End of namespace
} // End of namespace
///////////////////////////////////////////////////////////////////////////////
//
// Copyright (2018) Alexander Stukowski
//
// This file is part of OVITO (Open Visualization Tool).
//
// OVITO is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// OVITO is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
//
///////////////////////////////////////////////////////////////////////////////
#pragma once
#include <plugins/crystalanalysis/CrystalAnalysis.h>
#include <plugins/stdobj/gui/widgets/DataSeriesPlotWidget.h>
#include <gui/properties/ModifierPropertiesEditor.h>
#include <core/utilities/DeferredMethodInvocation.h>
class QwtPlotZoneItem;
namespace Ovito { namespace Plugins { namespace CrystalAnalysis {
/**
* Properties editor for the GrainSegmentationModifier class.
*/
class GrainSegmentationModifierEditor : public ModifierPropertiesEditor
{
Q_OBJECT
OVITO_CLASS(GrainSegmentationModifierEditor)
public:
/// Default constructor.
Q_INVOKABLE GrainSegmentationModifierEditor() {}
protected Q_SLOTS:
/// Replots the histogram computed by the modifier.
void plotHistogram();
/// Is called when the user clicks the 'Perform grain tracking' button.
void onPerformGrainTracking();
protected:
/// Creates the user interface controls for the editor.
virtual void createUI(const RolloutInsertionParameters& rolloutParams) override;
/// This method is called when a reference target changes.
virtual bool referenceEvent(RefTarget* source, const ReferenceEvent& event) override;
private:
/// The graph widget to display the RMSD histogram.
DataSeriesPlotWidget* _rmsdPlotWidget;
/// Marks the RMSD cutoff in the histogram plot.
QwtPlotZoneItem* _rmsdRangeIndicator;
/// For deferred invocation of the plot repaint function.
DeferredMethodInvocation<GrainSegmentationModifierEditor, &GrainSegmentationModifierEditor::plotHistogram> plotHistogramLater;
};
} // End of namespace
} // End of namespace
} // End of namespace
///////////////////////////////////////////////////////////////////////////////
//
// Copyright (2018) Alexander Stukowski
//
// This file is part of OVITO (Open Visualization Tool).
//
// OVITO is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// OVITO is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
//
///////////////////////////////////////////////////////////////////////////////
#pragma once
#include <plugins/crystalanalysis/CrystalAnalysis.h>
#include <plugins/particles/modifier/analysis/StructureIdentificationModifier.h>
#include <plugins/particles/objects/ParticlesObject.h>
#include <plugins/particles/objects/BondsObject.h>
#include <plugins/particles/modifier/analysis/ptm/PTMAlgorithm.h>
#include <boost/optional/optional.hpp>
namespace Ovito { namespace Plugins { namespace CrystalAnalysis {
/*
* Computation engine of the GrainSegmentationModifier, which decomposes a polycrystalline microstructure into individual grains.
*/
class GrainSegmentationEngine : public StructureIdentificationModifier::StructureIdentificationEngine
{
public:
/// Constructor.
GrainSegmentationEngine(
ParticleOrderingFingerprint fingerprint, ConstPropertyPtr positions, const SimulationCell& simCell,
const QVector<bool>& typesToIdentify, ConstPropertyPtr selection,
FloatType rmsdCutoff, FloatType misorientationThreshold,
int minGrainAtomCount, bool orphanAdoption, bool outputBonds);
/// Performs the computation.
virtual void perform() override;
/// This method is called by the system to free memory and release any working data after the
/// computation has been successfully completed.
virtual void cleanup() override {
_neighborLists.reset();
decltype(_distanceSortedAtoms){}.swap(_distanceSortedAtoms);
decltype(_clusterOrientations){}.swap(_clusterOrientations);
decltype(_clusterSizes){}.swap(_clusterSizes);
if(!_outputBonds) {
decltype(_latticeNeighborBonds){}.swap(_latticeNeighborBonds);
_neighborDisorientationAngles.reset();
}
StructureIdentificationEngine::cleanup();
}
/// Injects the computed results into the data pipeline.
virtual void emitResults(TimePoint time, ModifierApplication* modApp, PipelineFlowState& state) override;
/// Returns the per-atom RMSD values computed by the PTM algorithm.
const PropertyPtr& rmsd() const { return _rmsd; }
/// Returns the RMSD value range of the histogram.
FloatType rmsdHistogramRange() const { return _rmsdHistogramRange; }
/// Returns the histogram of computed RMSD values.
const PropertyPtr& rmsdHistogram() const { return _rmsdHistogram; }
/// Returns the computed per-particle lattice orientations.
const PropertyPtr& orientations() const { return _orientations; }
/// Returns the particle to cluster assignment.
const PropertyPtr& atomClusters() const { return _atomClusters; }
/// Returns the bonds generated between neighboring lattice atoms.
const std::vector<Bond>& latticeNeighborBonds() const { return _latticeNeighborBonds; }
/// Returns the computed disorientation angles between neighboring lattice atoms.
const PropertyPtr& neighborDisorientationAngles() const { return _neighborDisorientationAngles; }
private:
/// Performs the PTM algorithm. Determines the local structure type and the local lattice orientation.
bool identifyAtomicStructures();
/// Builds the graph of neighbor atoms and calculates the misorientation angle for each graph edge (i.e. bond).
bool buildNeighborGraph();
/// Merges adjacent clusters with similar lattice orientations.
bool mergeSuperclusters();
/// Computes the average lattice orientation of each cluster.
bool calculateAverageClusterOrientations();
/// Builds grains by iterative region merging
bool regionMerging();
/// Randomizes cluster IDs for testing purposes (giving more color contrast).
bool randomizeClusterIDs();
/// Merges any orphan atoms into the closest cluster.
bool mergeOrphanAtoms();
private:
/// Supercluster IDs
std::vector<size_t> _atomSuperclusters;
/// Counts the number of superclusters
size_t _numSuperclusters = 0;
/// Stores the number of atoms in each supercluster.
std::vector<size_t> _superclusterSizes;
/// The cutoff parameter used by the PTM algorithm.
FloatType _rmsdCutoff;
/// The merging criterion threshold.
FloatType _mergingThreshold;
/// The minimum number of crystalline atoms per grain.
int _minGrainAtomCount;
/// Whether to adopt orphan atoms.
int _orphanAdoption;
/// Stores the list of neighbors of each lattice atom.
PropertyPtr _neighborLists;
/// Stores indices of crystalline atoms sorted by value of distance transform.
std::vector<size_t> _distanceSortedAtoms;
/// Counts the number of clusters
size_t _numClusters = 0;
/// Stores the average lattice orientation of each cluster.
std::vector<Quaternion> _clusterOrientations;
/// Stores the number of atoms in each cluster.
std::vector<qlonglong> _clusterSizes;
/// The per-atom RMSD values computed by the PTM algorithm.
const PropertyPtr _rmsd;
/// Histogram of the RMSD values computed by the PTM algorithm.
PropertyPtr _rmsdHistogram;
/// The value range of the RMSD histogram.
FloatType _rmsdHistogramRange;
/// The computed per-particle lattice orientations.
PropertyPtr _orientations;
/// The particle to cluster assignment.
PropertyPtr _atomClusters;
/// Flag controlling the output of generated bonds to the data pipeline.
bool _outputBonds;
/// The bonds generated between neighboring lattice atoms.
std::vector<Bond> _latticeNeighborBonds;
/// The computed disorientation angles between neighboring lattice atoms.
PropertyPtr _neighborDisorientationAngles;
FloatType _misorientationThreshold = 4 / FloatType(180) * FLOATTYPE_PI;
};
} // End of namespace
} // End of namespace
} // End of namespace
///////////////////////////////////////////////////////////////////////////////
//
// Copyright (2018) Alexander Stukowski
//
// This file is part of OVITO (Open Visualization Tool).
//
// OVITO is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// OVITO is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
//
///////////////////////////////////////////////////////////////////////////////
#pragma once
#include <plugins/crystalanalysis/CrystalAnalysis.h>
#include <plugins/particles/modifier/analysis/StructureIdentificationModifier.h>
#include <plugins/particles/objects/BondsVis.h>
#include <plugins/particles/objects/ParticlesObject.h>
#include <plugins/stdobj/series/DataSeriesObject.h>
namespace Ovito { namespace Plugins { namespace CrystalAnalysis {
class GrainSegmentationEngine; // defined in GrainSegmentationEngine.h
/*
* Decomposes a polycrystalline microstructure into individual grains.
*/
class OVITO_CRYSTALANALYSIS_EXPORT GrainSegmentationModifier : public StructureIdentificationModifier
{
Q_OBJECT
OVITO_CLASS(GrainSegmentationModifier)
Q_CLASSINFO("DisplayName", "Grain segmentation");
Q_CLASSINFO("ModifierCategory", "Analysis");
public:
/// Constructor.
Q_INVOKABLE GrainSegmentationModifier(DataSet* dataset);
/// Performs grain tracking over the whole simulation trajectory.
bool trackGrains(TaskManager& taskManager, ModifierApplication* modApp);
protected:
/// Creates a computation engine that finds the grains in a single frame.
std::shared_ptr<GrainSegmentationEngine> createSegmentationEngine(TimePoint time, ModifierApplication* modApp, const PipelineFlowState& input);
/// Creates a computation engine that will compute the modifier's results.
virtual Future<ComputeEnginePtr> createEngine(TimePoint time, ModifierApplication* modApp, const PipelineFlowState& input) override;
private:
/// The RMSD cutoff for the PTM algorithm.
DECLARE_MODIFIABLE_PROPERTY_FIELD_FLAGS(FloatType, rmsdCutoff, setRmsdCutoff, PROPERTY_FIELD_MEMORIZE);
/// Controls the amount of noise allowed inside a grain.
DECLARE_MODIFIABLE_PROPERTY_FIELD(FloatType, mergingThreshold, setMergingThreshold);
/// The minimum number of crystalline atoms per grain.
DECLARE_MODIFIABLE_PROPERTY_FIELD(int, minGrainAtomCount, setMinGrainAtomCount);
/// Controls whether to adopt orphan atoms
DECLARE_MODIFIABLE_PROPERTY_FIELD_FLAGS(bool, orphanAdoption, setOrphanAdoption, PROPERTY_FIELD_MEMORIZE);
/// Controls whether only selected particles should be taken into account.
DECLARE_MODIFIABLE_PROPERTY_FIELD(bool, onlySelectedParticles, setOnlySelectedParticles);
/// The visual element for rendering the bonds created by the modifier.
DECLARE_MODIFIABLE_REFERENCE_FIELD_FLAGS(BondsVis, bondsVis, setBondsVis, PROPERTY_FIELD_DONT_PROPAGATE_MESSAGES | PROPERTY_FIELD_MEMORIZE);
/// Controls the output of bonds by the modifier.
DECLARE_MODIFIABLE_PROPERTY_FIELD(bool, outputBonds, setOutputBonds);
};
} // End of namespace
} // End of namespace
} // End of namespace
......@@ -36,7 +36,7 @@ StructureListParameterUI::StructureListParameterUI(PropertiesEditor* parentEdito
: RefTargetListParameterUI(parentEditor, PROPERTY_FIELD(StructureIdentificationModifier::structureTypes)),
_showCheckBoxes(showCheckBoxes)
{
connect(tableWidget(230), &QTableWidget::doubleClicked, this, &StructureListParameterUI::onDoubleClickStructureType);
connect(tableWidget(300), &QTableWidget::doubleClicked, this, &StructureListParameterUI::onDoubleClickStructureType);
connect(parentEditor, &PropertiesEditor::contentsReplaced, tableWidget(), &QTableView::resizeRowsToContents);
tableWidget()->setAutoScroll(false);
......
......@@ -22,7 +22,6 @@
#include <plugins/particles/Particles.h>
#include "PTMAlgorithm.h"
#include <3rdparty/ptm/ptm_functions.h>
namespace Ovito { namespace Particles { OVITO_BEGIN_INLINE_NAMESPACE(Util)
......@@ -62,7 +61,7 @@ typedef struct
} ptmnbrdata_t;
static int get_neighbours(void* vdata, size_t _unused_lammps_variable, size_t atom_index, int num_requested, int* ordering, size_t* nbr_indices, int32_t* numbers, double (*nbr_pos)[3])
static int get_neighbours(void* vdata, size_t _unused_lammps_variable, size_t atom_index, int num_requested, ptm_atomicenv_t* env)
{
ptmnbrdata_t* nbrdata = (ptmnbrdata_t*)vdata;
const NearestNeighborFinder* neighFinder = nbrdata->neighFinder;
......@@ -79,31 +78,34 @@ static int get_neighbours(void* vdata, size_t _unused_lammps_variable, size_t at
ptm_index_to_permutation(numNeighbors, precachedNeighbors[atom_index], permutation);
// Bring neighbor coordinates into a form suitable for the PTM library.
ordering[0] = 0;
nbr_indices[0] = atom_index;
nbr_pos[0][0] = nbr_pos[0][1] = nbr_pos[0][2] = 0;
env->correspondences[0] = 0;
env->atom_indices[0] = atom_index;
env->points[0][0] = 0;
env->points[0][1] = 0;
env->points[0][2] = 0;
for(int i = 0; i < numNeighbors; i++) {
//ordering[index] = permutation[i] + 1;
nbr_indices[i+1] = neighQuery.results()[permutation[i]].index;
nbr_pos[i+1][0] = neighQuery.results()[permutation[i]].delta.x();
nbr_pos[i+1][1] = neighQuery.results()[permutation[i]].delta.y();
nbr_pos[i+1][2] = neighQuery.results()[permutation[i]].delta.z();
//env->correspondences[index] = permutation[i] + 1;
env->atom_indices[i+1] = neighQuery.results()[permutation[i]].index;
env->points[i+1][0] = neighQuery.results()[permutation[i]].delta.x();
env->points[i+1][1] = neighQuery.results()[permutation[i]].delta.y();
env->points[i+1][2] = neighQuery.results()[permutation[i]].delta.z();
}
// Build list of particle types for ordering identification.
if(particleTypes != nullptr) {
numbers[0] = particleTypes->getInt(atom_index);
env->numbers[0] = particleTypes->getInt(atom_index);
for(int i = 0; i < numNeighbors; i++) {
numbers[i+1] = particleTypes->getInt(neighQuery.results()[permutation[i]].index);
env->numbers[i+1] = particleTypes->getInt(neighQuery.results()[permutation[i]].index);
}
}
else {
for(int i = 0; i < numNeighbors + 1; i++) {
numbers[i] = 0;
env->numbers[i] = 0;
}
}
env->num = numNeighbors + 1;
return numNeighbors + 1;
}
......@@ -159,7 +161,7 @@ PTMAlgorithm::StructureType PTMAlgorithm::Kernel::identifyStructure(size_t parti
nullptr,
&_bestTemplateIndex,
&_bestTemplate,
_correspondences);
&_env);
// Convert PTM classification back to our own scheme.
if(type == PTM_MATCH_NONE || (_algo._rmsdCutoff != 0 && _rmsd > _algo._rmsdCutoff)) {
......@@ -187,18 +189,20 @@ PTMAlgorithm::StructureType PTMAlgorithm::Kernel::identifyStructure(size_t parti
}
}
#if 0
if (_structureType != OTHER && qtarget != nullptr) {
//arrange orientation in PTM format
double qtarget_ptm[4] = { qtarget->w(), qtarget->x(), qtarget->y(), qtarget->z() };
double disorientation = 0; //TODO: this is probably not needed
int template_index = ptm_remap_template(type, true, _bestTemplateIndex, qtarget_ptm, _q, &disorientation, _correspondences, &_bestTemplate);
int template_index = ptm_remap_template(type, true, _bestTemplateIndex, qtarget_ptm, _q, &disorientation, _env.correspondences, &_bestTemplate);
if (template_index < 0)
return _structureType;
_bestTemplateIndex = template_index;
}
#endif
return _structureType;
......@@ -294,7 +298,7 @@ const NearestNeighborFinder::Neighbor& PTMAlgorithm::Kernel::getNeighborInfo(int
{
OVITO_ASSERT(_structureType != OTHER);
OVITO_ASSERT(index >= 0 && index < numStructureNeighbors());
int mappedIndex = _correspondences[index + 1] - 1;
int mappedIndex = _env.correspondences[index + 1] - 1;
OVITO_ASSERT(mappedIndex >= 0 && mappedIndex < results().size());
return results()[mappedIndex];
}
......
......@@ -25,6 +25,7 @@
#include <plugins/particles/Particles.h>
#include <plugins/stdobj/properties/PropertyStorage.h>
#include <plugins/particles/util/NearestNeighborFinder.h>
#include <3rdparty/ptm/ptm_functions.h>
extern "C" {
typedef struct ptm_local_handle* ptm_local_handle_t;
......@@ -51,9 +52,9 @@ public:
HCP, //< Hexagonal close-packed
BCC, //< Body-centered cubic
ICO, //< Icosahedral structure
SC, //< Simple cubic structure
CUBIC_DIAMOND, //< Cubic diamond structure
HEX_DIAMOND, //< Hexagonal diamond structure
SC, //< Simple cubic structure
CUBIC_DIAMOND, //< Cubic diamond structure
HEX_DIAMOND, //< Hexagonal diamond structure
GRAPHENE, //< Graphene structure
NUM_STRUCTURE_TYPES //< This counts the number of defined structure types.
......@@ -199,6 +200,9 @@ public:
/// identified for the current particle.
const Vector_3<double>& getIdealNeighborVector(int index) const;
//TODO: don't leave this public
ptm_atomicenv_t _env;
private:
/// Reference to the parent algorithm object.
const PTMAlgorithm& _algo;
......@@ -216,8 +220,8 @@ public:
int32_t _orderingType = ORDERING_NONE;
int _bestTemplateIndex;
const double (*_bestTemplate)[3] = nullptr;
int8_t _correspondences[MAX_INPUT_NEIGHBORS+1];
std::vector<uint64_t> _cachedNeighbors;
//int8_t _correspondences[MAX_INPUT_NEIGHBORS+1];
std::vector<uint64_t> _cachedNeighbors;
};
private:
......