Skip to content

Commit 672ba41

Browse files
stevenlabellestevenlabelle
committed
Image specific ranges, ImageMap validation, plot/log point validation
Incremental changes to allow analysis on only subvolume of mesh with template and target of different size and position. Added separate physical dimension ranges for template and target images. Standard constraints automatically assign the same range to both for backwards compatibility. Filter based template now requires each. ImageMap validation function used for forces, energy, diff, etc. Plot/log file values check validation function and return 0 for data outside image range.
1 parent 6c7671e commit 672ba41

7 files changed

Lines changed: 114 additions & 52 deletions

‎FEWarp/FEWarpFilteredImageConstraint.cpp‎

Lines changed: 15 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -15,8 +15,11 @@ FEWarpFilteredImageConstraint::FEWarpFilteredImageConstraint(FEModel* pfem) : FE
1515
{
1616
m_pt = 0.0;
1717

18-
m_r0[0] = m_r0[1] = m_r0[2] = 0.0;
19-
m_r1[0] = m_r1[1] = m_r1[2] = 1.0;
18+
m_tr0[0] = m_tr0[1] = m_tr0[2] = 0.0;
19+
m_tr1[0] = m_tr1[1] = m_tr1[2] = 1.0;
20+
m_sr0[0] = m_sr0[1] = m_sr0[2] = 0.0;
21+
m_sr1[0] = m_sr1[1] = m_sr1[2] = 1.0;
22+
2023
}
2124

2225
//-----------------------------------------------------------------------------
@@ -31,11 +34,13 @@ bool FEWarpFilteredImageConstraint::Init()
3134
int ny = m_tmp0.height();
3235
int nz = m_tmp0.depth ();
3336

34-
vec3d r0(m_r0[0], m_r0[1], m_r0[2]);
35-
vec3d r1(m_r1[0], m_r1[1], m_r1[2]);
37+
vec3d tr0(m_tr0[0], m_tr0[1], m_tr0[2]);
38+
vec3d tr1(m_tr1[0], m_tr1[1], m_tr1[2]);
39+
vec3d sr0(m_sr0[0], m_sr0[1], m_sr0[2]);
40+
vec3d sr1(m_sr1[0], m_sr1[1], m_sr1[2]);
3641

37-
m_tmap.SetRange(r0, r1);
38-
m_smap.SetRange(r0, r1);
42+
m_tmap.SetRange(tr0, tr1);
43+
m_smap.SetRange(sr0, sr1);
3944

4045
m_tmp = m_tmp0;
4146
m_trg = m_trg0;
@@ -100,8 +105,10 @@ BEGIN_FECORE_CLASS(FEWarpSingleFilteredImageConstraint, FEWarpVolumeConstraint);
100105
ADD_PARAMETER(m_k , "penalty" );
101106
ADD_PARAMETER(m_blaugon, "laugon" );
102107
ADD_PARAMETER(m_altol , "altol" );
103-
ADD_PARAMETER(m_r0 , 3, "range_min");
104-
ADD_PARAMETER(m_r1 , 3, "range_max");
108+
ADD_PARAMETER(m_tr0 , 3, "template_range_min");
109+
ADD_PARAMETER(m_tr1 , 3, "template_range_max");
110+
ADD_PARAMETER(m_sr0, 3, "target_range_min");
111+
ADD_PARAMETER(m_sr1, 3, "target_range_max");
105112

106113
ADD_PROPERTY(m_tmpReader, "template")->SetDefaultType("raw");
107114
ADD_PROPERTY(m_trgReader, "target" )->SetDefaultType("raw");

‎FEWarp/FEWarpImageConstraint.cpp‎

Lines changed: 19 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -11,8 +11,10 @@ FEWarpImageConstraint::FEWarpImageConstraint(FEModel* pfem) : FEWarpVolumeConstr
1111
m_blur_cur = 0.0;
1212
m_blur_method = 0; // default average blur
1313

14-
m_r0[0] = m_r0[1] = m_r0[2] = 0.0;
15-
m_r1[0] = m_r1[1] = m_r1[2] = 1.0;
14+
m_tr0[0] = m_tr0[1] = m_tr0[2] = 0.0;
15+
m_tr1[0] = m_tr1[1] = m_tr1[2] = 1.0;
16+
m_sr0[0] = m_sr0[1] = m_sr0[2] = 0.0;
17+
m_sr1[0] = m_sr1[1] = m_sr1[2] = 1.0;
1618
}
1719

1820
//-----------------------------------------------------------------------------
@@ -27,11 +29,16 @@ bool FEWarpImageConstraint::Init()
2729
int ny = m_tmp0.height();
2830
int nz = m_tmp0.depth ();
2931

30-
vec3d r0(m_r0[0], m_r0[1], m_r0[2]);
31-
vec3d r1(m_r1[0], m_r1[1], m_r1[2]);
32+
vec3d tr0(m_tr0[0], m_tr0[1], m_tr0[2]);
33+
vec3d tr1(m_tr1[0], m_tr1[1], m_tr1[2]);
34+
// SL: Leave this unchanged for now by copying the template range.
35+
m_sr0[0] = m_tr0[0]; m_sr0[1] = m_tr0[1]; m_sr0[2] = m_tr0[2];
36+
m_sr1[0] = m_tr1[0]; m_sr1[1] = m_tr1[1]; m_sr1[2] = m_tr1[2];
37+
vec3d sr0(m_sr0[0], m_sr0[1], m_sr0[2]);
38+
vec3d sr1(m_sr1[0], m_sr1[1], m_sr1[2]);
3239

33-
m_tmap.SetRange(r0, r1);
34-
m_smap.SetRange(r0, r1);
40+
m_tmap.SetRange(tr0, tr1);
41+
m_smap.SetRange(sr0, sr1);
3542

3643
m_tmp = m_tmp0;
3744
m_trg = m_trg0;
@@ -76,8 +83,12 @@ BEGIN_FECORE_CLASS(FEWarpSingleImageConstraint, FEWarpConstraint);
7683
ADD_PARAMETER(m_altol , "altol" );
7784
ADD_PARAMETER(m_blur , "blur" );
7885
ADD_PARAMETER(m_blur_method , "blur_method" )->setEnums("AVERAGE\0FFT\0");
79-
ADD_PARAMETER(m_r0 , 3, "range_min");
80-
ADD_PARAMETER(m_r1 , 3, "range_max");
86+
//SL: Leave parameter name unchanged for now.
87+
ADD_PARAMETER(m_tr0 , 3, "range_min");
88+
ADD_PARAMETER(m_tr1 , 3, "range_max");
89+
//SL: Don't expose for now to leave this unchanged.
90+
//ADD_PARAMETER(m_sr0, 3, "target_range_min");
91+
//ADD_PARAMETER(m_sr1, 3, "target_range_max");
8192

8293
ADD_PROPERTY(m_tmpReader, "template")->SetDefaultType("raw");
8394
ADD_PROPERTY(m_trgReader, "target" )->SetDefaultType("raw");

‎FEWarp/FEWarpLog.cpp‎

Lines changed: 23 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -20,12 +20,14 @@ double FELogWarpTemplate::value(FEElement& el)
2020
FEMesh& mesh = GetFEModel()->GetMesh();
2121
ImageMap& tmap = m_wrp->GetTemplateMap();
2222
double T = 0.0;
23+
int n = 0;
2324
for (int i = 0; i < el.Nodes(); ++i)
2425
{
2526
vec3d r0 = mesh.Node(el.m_node[i]).m_r0;
26-
T += tmap.value(r0);
27+
T += tmap.valid(r0) ? tmap.value(r0) : 0.0;
28+
n += tmap.valid(r0) ? 1 : 0;
2729
}
28-
T /= (double)el.Nodes();
30+
T /= (double)n;
2931
return T;
3032
}
3133

@@ -35,12 +37,15 @@ double FELogWarpTarget::value(FEElement& el)
3537
FEMesh& mesh = GetFEModel()->GetMesh();
3638
ImageMap& smap = m_wrp->GetTargetMap();
3739
double S = 0.0;
40+
int n = 0;
41+
3842
for (int i = 0; i < el.Nodes(); ++i)
3943
{
4044
vec3d rt = mesh.Node(el.m_node[i]).m_rt;
41-
S += smap.value(rt);
45+
S += smap.valid(rt) ? smap.value(rt) : 0.0;
46+
n += smap.valid(rt) ? 1 : 0;
4247
}
43-
S /= (double)el.Nodes();
48+
S /= (double)n;
4449
return S;
4550
}
4651

@@ -51,15 +56,17 @@ double FELogWarpEnergy::value(FEElement& el)
5156
ImageMap& tmap = m_wrp->GetTemplateMap();
5257
ImageMap& smap = m_wrp->GetTargetMap();
5358
double e = 0.0;
59+
int n = 0;
5460
for (int i = 0; i < el.Nodes(); ++i)
5561
{
5662
vec3d r0 = mesh.Node(el.m_node[i]).m_r0;
5763
vec3d rt = mesh.Node(el.m_node[i]).m_rt;
58-
double T = tmap.value(r0);
59-
double S = smap.value(rt);
60-
e += (0.5 * (T - S) * (T - S));
64+
double T = tmap.valid(r0) ? tmap.value(r0) : 0.0;
65+
double S = tmap.valid(r0) ? smap.value(rt) : 0.0;
66+
e += (tmap.valid(r0) || smap.valid(rt)) ? (0.5 * (T - S) * (T - S)) : 0.0;
67+
n += (tmap.valid(r0) || smap.valid(rt)) ? 1 : 0;
6168
}
62-
e /= (double)el.Nodes();
69+
e /= (double)n;
6370
return e;
6471
}
6572

@@ -74,9 +81,9 @@ vec3d FELogWarpForce_::force(FEElement& el)
7481
{
7582
vec3d r0 = mesh.Node(el.m_node[i]).m_r0;
7683
vec3d rt = mesh.Node(el.m_node[i]).m_rt;
77-
double T = tmap.value(r0);
78-
double S = smap.value(rt);
79-
vec3d G = smap.gradient(rt);
84+
double T = tmap.valid(r0) ? tmap.value(r0) : 0.0;
85+
double S = smap.valid(rt) ? smap.value(rt) : 0.0;
86+
vec3d G = smap.valid(rt) ? smap.gradient(rt) : vec3d(0.0);
8087
fw += G * ((T - S));
8188
}
8289
fw /= (double)el.Nodes();
@@ -108,13 +115,15 @@ double FELogWarpDiff::value(FEElement& el)
108115
ImageMap& tmap = m_wrp->GetTemplateMap();
109116
ImageMap& smap = m_wrp->GetTargetMap();
110117
double e = 0.0;
118+
int n = 0;
111119
for (int i = 0; i < el.Nodes(); ++i)
112120
{
113121
vec3d r0 = mesh.Node(el.m_node[i]).m_r0;
114122
vec3d rt = mesh.Node(el.m_node[i]).m_rt;
115-
double T = tmap.value(r0);
116-
double S = smap.value(rt);
117-
e += S - T;
123+
double T = tmap.valid(r0) ? tmap.value(r0) : 0.0;
124+
double S = smap.valid(rt) ? smap.value(rt) : 0.0;
125+
e += (tmap.valid(r0) || smap.valid(rt)) ? S - T : 0.0;
126+
n += (tmap.valid(r0) || smap.valid(rt)) ? 1 : 0.0;
118127
}
119128
e /= (double)el.Nodes();
120129
return e;

‎FEWarp/FEWarpMultiImageConstraint.cpp‎

Lines changed: 8 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -11,8 +11,14 @@ BEGIN_FECORE_CLASS(FEWarpMultiImageConstraint, FEWarpImageConstraint);
1111
ADD_PARAMETER(m_blaugon, "laugon" );
1212
ADD_PARAMETER(m_altol , "altol" );
1313
ADD_PARAMETER(m_blur , "blur" );
14-
ADD_PARAMETER(m_r0 , 3, "range_min");
15-
ADD_PARAMETER(m_r1 , 3, "range_max");
14+
//SL: Leave parameter name for now.
15+
ADD_PARAMETER(m_tr0 , 3, "range_min");
16+
ADD_PARAMETER(m_tr1 , 3, "range_max");
17+
//SL: Don't expose for now to leave this unchanged.
18+
//ADD_PARAMETER(m_tr0 , 3, "template_range_min");
19+
//ADD_PARAMETER(m_tr1 , 3, "template_range_max");
20+
//ADD_PARAMETER(m_sr0, 3, "target_range_min");
21+
//ADD_PARAMETER(m_sr1, 3, "target_range_max");
1622

1723
ADD_PROPERTY(m_tmpReader, "template")->SetDefaultType("raw");
1824
ADD_PROPERTY(m_trgReader, "target" )->SetDefaultType("raw");

‎FEWarp/FEWarpPlot.cpp‎

Lines changed: 28 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -28,7 +28,14 @@ bool FEPlotTemplate::SaveWarpImage(FEMesh& m, FEWarpVolumeConstraint* pc, FEData
2828
ImageMap& tmap = pc->GetTemplateMap();
2929

3030
int N = m.Nodes();
31-
for (int i=0; i<N; ++i) s << tmap.value(m.Node(i).m_r0);
31+
for (int i = 0; i < N; ++i)
32+
{
33+
vec3d r0 = m.Node(i).m_r0;
34+
if (tmap.valid(r0))
35+
s << tmap.value(r0);
36+
else
37+
s << 0.0;
38+
}
3239
return true;
3340
}
3441

@@ -67,7 +74,14 @@ bool FEPlotTarget::SaveWarpImage(FEMesh& m, FEWarpVolumeConstraint* pc, FEDataSt
6774
ImageMap& smap = pc->GetTargetMap();
6875

6976
int N = m.Nodes();
70-
for (int i=0; i<N; ++i) s << smap.value(m.Node(i).m_rt);
77+
for (int i = 0; i < N; ++i)
78+
{
79+
vec3d rt = m.Node(i).m_rt;
80+
if (smap.valid(rt))
81+
s << smap.value(rt);
82+
else
83+
s << 0.0;
84+
}
7185
return true;
7286
}
7387

@@ -102,8 +116,11 @@ bool FEPlotEnergy::Save(FEMesh &m, FEDataStream& s)
102116
int N = m.Nodes();
103117
for (int i=0; i<N; ++i)
104118
{
105-
double T = tmap.value(m.Node(i).m_r0);
106-
double S = smap.value(m.Node(i).m_rt);
119+
vec3d r0 = m.Node(i).m_r0;
120+
vec3d rt = m.Node(i).m_rt;
121+
122+
double T = tmap.valid(r0) ? tmap.value(r0) : 0.0;
123+
double S = smap.valid(rt) ? smap.value(rt) : 0.0;
107124

108125
s << (0.5*(T - S)*(T - S));
109126
}
@@ -132,9 +149,9 @@ bool FEPlotForce::Save(FEMesh &m, FEDataStream& s)
132149
vec3d r0 = m.Node(i).m_r0;
133150
vec3d rt = m.Node(i).m_rt;
134151

135-
double T = tmap.value(r0);
136-
double S = smap.value(rt);
137-
vec3d G = smap.gradient(rt);
152+
double T = tmap.valid(r0) ? tmap.value(r0) : 0.0;
153+
double S = smap.valid(rt) ? smap.value(rt) : 0.0;
154+
vec3d G = smap.valid(rt) ? smap.gradient(rt) : vec3d(0.0);
138155
vec3d fw = G*((T - S));
139156

140157
s << fw;
@@ -161,8 +178,10 @@ bool FEPlotDiff::Save(FEMesh& m, FEDataStream& s)
161178
int N = m.Nodes();
162179
for (int i = 0; i < N; ++i)
163180
{
164-
double T = tmap.value(m.Node(i).m_r0);
165-
double S = smap.value(m.Node(i).m_rt);
181+
vec3d r0 = m.Node(i).m_r0;
182+
vec3d rt = m.Node(i).m_rt;
183+
double T = tmap.valid(r0) ? tmap.value(r0) : 0.0;
184+
double S = smap.valid(rt) ? smap.value(rt) : 0.0;
166185

167186
s << S - T;
168187
}

‎FEWarp/FEWarpVolumeConstraint.cpp‎

Lines changed: 17 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -7,8 +7,10 @@
77
//-----------------------------------------------------------------------------
88
FEWarpVolumeConstraint::FEWarpVolumeConstraint(FEModel* pfem) : FEWarpConstraint(pfem), m_tmap(m_tmp), m_smap(m_trg)
99
{
10-
m_r0[0] = m_r0[1] = m_r0[2] = 0.0;
11-
m_r1[0] = m_r1[1] = m_r1[2] = 1.0;
10+
m_tr0[0] = m_tr0[1] = m_tr0[2] = 0.0;
11+
m_tr1[0] = m_tr1[1] = m_tr1[2] = 1.0;
12+
m_sr0[0] = m_sr0[1] = m_sr0[2] = 0.0;
13+
m_sr1[0] = m_sr1[1] = m_sr1[2] = 1.0;
1214
}
1315

1416
//-----------------------------------------------------------------------------
@@ -17,14 +19,17 @@ FEWarpVolumeConstraint::~FEWarpVolumeConstraint(void) {}
1719
//-----------------------------------------------------------------------------
1820
vec3d FEWarpVolumeConstraint::wrpForce(FEMaterialPoint& mp)
1921
{
22+
vec3d r0 = mp.m_r0;
23+
vec3d rt = mp.m_rt;
24+
2025
// evaluate template
21-
double T = m_tmap.value(mp.m_r0);
26+
double T = m_tmap.valid(r0) ? m_tmap.value(mp.m_r0) : 0.0;
2227

2328
// evaluate target
24-
double S = m_smap.value(mp.m_rt);
29+
double S = m_smap.valid(rt) ? m_smap.value(mp.m_rt) : 0.0;
2530

2631
// evaluate target gradient
27-
vec3d G = m_smap.gradient(mp.m_rt);
32+
vec3d G = m_smap.valid(rt) ? m_smap.gradient(mp.m_rt) : vec3d(0.0);
2833

2934
// evaluate force
3035
vec3d Fw = G*((S - T)*m_k);
@@ -35,17 +40,20 @@ vec3d FEWarpVolumeConstraint::wrpForce(FEMaterialPoint& mp)
3540
//-----------------------------------------------------------------------------
3641
mat3ds FEWarpVolumeConstraint::wrpStiffness(FEMaterialPoint& mp)
3742
{
43+
vec3d r0 = mp.m_r0;
44+
vec3d rt = mp.m_rt;
45+
3846
// template value
39-
double T = m_tmap.value(mp.m_r0);
47+
double T = m_tmap.valid(r0) ? m_tmap.value(mp.m_r0) : 0.0;
4048

4149
// target value
42-
double S = m_smap.value(mp.m_rt);
50+
double S = m_smap.valid(rt) ? m_smap.value(mp.m_rt) : 0.0;
4351

4452
// calculate target gradient
45-
vec3d dS = m_smap.gradient(mp.m_rt);
53+
vec3d dS = m_smap.valid(rt) ? m_smap.gradient(mp.m_rt) : vec3d(0.0);
4654

4755
// calculate target hessian
48-
mat3ds H = m_smap.hessian(mp.m_rt);
56+
mat3ds H = m_smap.valid(rt) ? m_smap.hessian(mp.m_rt) : mat3ds(0.0);
4957

5058
// warping stiffness
5159
return H*((T - S)*m_k) - dyad(dS)*m_k;

‎FEWarp/FEWarpVolumeConstraint.h‎

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -25,8 +25,10 @@ class FEWarpVolumeConstraint : public FEWarpConstraint
2525
mat3ds wrpStiffness(FEMaterialPoint& pt);
2626

2727
protected:
28-
double m_r0[3]; //!< minimum range
29-
double m_r1[3]; //!< maximum range
28+
double m_tr0[3]; //!< minimum range
29+
double m_tr1[3]; //!< maximum range
30+
double m_sr0[3]; //!< minimum range
31+
double m_sr1[3]; //!< maximum range
3032

3133

3234
Image m_tmp; //!< template image

0 commit comments

Comments
 (0)