in UVAtlas/isochart/imtcomputation.cpp [1046:1229]
static HRESULT Accumulation(
DOUBLEVECTOR2* corner,
float* pfSignal,
size_t dwSignalDimension,
std::vector<DOUBLEVECTOR2>& above, // above line when applying 2 times accumulation
std::vector<DOUBLEVECTOR2>& below,
double IMTResult[],
double& dPieceArea)// below line when applying 2 times accumulation)
{
HRESULT hr = S_OK;
dPieceArea = 0;
if (above.size() < 2 || below.size() < 2)
{
return hr;
}
assert(above[0].x == below[0].x);
std::unique_ptr<double[]> m(new (std::nothrow) double[4 * dwSignalDimension]);
if (!m)
{
return E_OUTOFMEMORY;
}
double* m1 = m.get();
double* m2 = m1 + dwSignalDimension;
double* m3 = m2 + dwSignalDimension;
double* m4 = m3 + dwSignalDimension;
/*
// Following is a texel,
// (c) _ _ _ _ (d)
// | |
// | |
// | |
// | |
// (a)-------(b)
*/
for (size_t ii = 0; ii < dwSignalDimension; ii++)
{
auto a = double(pfSignal[0 * dwSignalDimension + ii]);
auto b = double(pfSignal[1 * dwSignalDimension + ii]);
auto c = double(pfSignal[2 * dwSignalDimension + ii]);
auto d = double(pfSignal[3 * dwSignalDimension + ii]);
m1[ii] = a + d - c - b;
m2[ii] = (b - a) * corner[1].y + (c - d) * corner[0].y;
m3[ii] = a + d - c - b;
m4[ii] = (c - a) * corner[1].x + (b - d) * corner[0].x;
}
// C
uint32_t dwCurrA = 0, dwCurrB = 0;
uint32_t dwNextA = NextItegralPoint(above, dwCurrA);
uint32_t dwNextB = NextItegralPoint(below, dwCurrB);
double a2 = 0, b2 = 0, a1 = 0, b1 = 0;
bool bNewSegmentA = true;
bool bNewSegmentB = true;
double fStartX = above[dwCurrA].x;
double fEndX = fStartX;
while (dwNextA != INVALID_INDEX && dwNextB != INVALID_INDEX)
{
fStartX = fEndX;
if (bNewSegmentA)
{
CalculateLineParameters(
above[dwCurrA],
above[dwNextA],
a2,
b2);
bNewSegmentA = false;
}
if (bNewSegmentB)
{
CalculateLineParameters(
below[dwCurrB],
below[dwNextB],
a1,
b1);
bNewSegmentB = false;
}
if (IsInZeroRangeDouble(above[dwNextA].x - below[dwNextB].x))
{
fEndX = above[dwNextA].x;
dwCurrA = dwNextA;
dwCurrB = dwNextB;
dwNextA = NextItegralPoint(above, dwCurrA);
dwNextB = NextItegralPoint(below, dwCurrB);
bNewSegmentA = true;
bNewSegmentB = true;
}
else if (above[dwNextA].x < below[dwNextB].x)
{
fEndX = above[dwNextA].x;
dwCurrA = dwNextA;
dwNextA = NextItegralPoint(above, dwCurrA);
bNewSegmentA = true;
}
else
{
fEndX = below[dwNextB].x;
dwCurrB = dwNextB;
dwNextB = NextItegralPoint(below, dwCurrB);
bNewSegmentB = true;
}
double aa1 = a1 * a1;
double aaa1 = aa1 * a1;
double aa2 = a2 * a2;
double aaa2 = aa2 * a2;
double bb1 = b1 * b1;
double bbb1 = bb1 * b1;
double bb2 = b2 * b2;
double bbb2 = bb2 * b2;
double u1 = fStartX;
double uu1 = u1 * u1;
double uuu1 = uu1 * u1;
double uuuu1 = uu1 * uu1;
double u2 = fEndX;
double uu2 = u2 * u2;
double uuu2 = uu2 * u2;
double uuuu2 = uu2 * uu2;
dPieceArea += (a2 - a1) * (uu2 - uu1) / 2 + (b2 - b1) * (u2 - u1);
for (size_t ii = 0; ii < dwSignalDimension; ii++)
{
double n3 = m1[ii] * m1[ii] * (aaa2 - aaa1) / 3;
double n2 = m1[ii] * m1[ii] * (aa2 * b2 - aa1 * b1) + m1[ii] * m2[ii] * (aa2 - aa1);
double n1 = m1[ii] * m1[ii] * (a2 * bb2 - a1 * bb1) + 2 * m1[ii] * m2[ii] * (a2 * b2 - a1 * b1) + m2[ii] * m2[ii] * (a2 - a1);
double n0 = m1[ii] * m1[ii] * (bbb2 - bbb1) / 3 + m1[ii] * m2[ii] * (bb2 - bb1) + m2[ii] * m2[ii] * (b2 - b1);
double fTemp =
n3 * (uuuu2 - uuuu1) / 4 +
n2 * (uuu2 - uuu1) / 3 +
n1 * (uu2 - uu1) / 2 +
n0 * (u2 - u1);
if (fTemp < 0) fTemp = 0;//Theoritically, the result must larger than 0
IMTResult[0] += fTemp;
n3 = m3[ii] * m3[ii] * (a2 - a1);
n2 = 2 * m3[ii] * m4[ii] * (a2 - a1) + m3[ii] * m3[ii] * (b2 - b1);
n1 = 2 * m3[ii] * m4[ii] * (b2 - b1) + m4[ii] * m4[ii] * (a2 - a1);
n0 = m4[ii] * m4[ii] * (b2 - b1);
fTemp =
n3 * (uuuu2 - uuuu1) / 4 +
n2 * (uuu2 - uuu1) / 3 +
n1 * (uu2 - uu1) / 2 +
n0 * (u2 - u1);
if (fTemp < 0) fTemp = 0; //Theoritically, the result must larger than 0
IMTResult[2] += fTemp;
n3 = m1[ii] * m3[ii] * (aa2 - aa1) / 2;
n2 = m1[ii] * m4[ii] * (aa2 - aa1) / 2 + m1[ii] * m3[ii] * (a2 * b2 - a1 * b1) + m2[ii] * m3[ii] * (a2 - a1);
n1 = m1[ii] * m3[ii] * (bb2 - bb1) / 2 + m1[ii] * m4[ii] * (a2 * b2 - a1 * b1) + m2[ii] * m4[ii] * (a2 - a1) + m2[ii] * m3[ii] * (b2 - b1);
n0 = m1[ii] * m4[ii] * (bb2 - bb1) / 2 + m2[ii] * m4[ii] * (b2 - b1);
IMTResult[1] +=
n3 * (uuuu2 - uuuu1) / 4 +
n2 * (uuu2 - uuu1) / 3 +
n1 * (uu2 - uu1) / 2 +
n0 * (u2 - u1);
}
}
double fPixelSize =
(corner[1].x - corner[0].x) * (corner[1].y - corner[0].y);
IMTResult[0] /= (fPixelSize * fPixelSize);
IMTResult[1] /= (fPixelSize * fPixelSize);
IMTResult[2] /= (fPixelSize * fPixelSize);
return hr;
}