00001
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 bool CTetraVertexDetector<CELL_MAX>::operator()(tVolume *pVolume,
00025 vctl::MCVerticeS *pVertices,
00026 tScale *pScale
00027 )
00028 {
00029 MDS_ASSERT(pVolume && pVertices && pScale);
00030
00031 MDS_LOG_NOTE("CTetraVertexDetector::operator()():");
00032
00033
00034 mds::tSize XSize = pVolume->getXSize();
00035 mds::tSize YSize = pVolume->getYSize();
00036 mds::tSize ZSize = pVolume->getZSize();
00037
00038
00039
00040 double dThreshold = 0.0;
00041
00042
00043 double dMin = mds::img::getMin<double>(*pVolume);
00044 double dMax = mds::img::getMax<double>(*pVolume);
00045 double dTemp = 1.0 / (dMax - dMin + 0.001);
00046
00047
00048 int iCount = 0;
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072 vctl::MCPoint3D Point;
00073 mds::tSize i, j, k;
00074 for( k = 0; k < ZSize; ++k )
00075 {
00076 for( j = 0; j < YSize; ++j )
00077 {
00078 for( i = 0; i < XSize; ++i )
00079 {
00080 tVoxel Value = pVolume->get(i, j, k);
00081 if( Value == mds::img::CPixelTraits<tVoxel>::getPixelMax() )
00082 {
00083 pScale->conv2Real(i, j, k, &Point);
00084 vctl::MCVertex *pVertex = pVertices->New(Point);
00085 pVertex->SetValue((double(Value) - dMin) * dTemp);
00086 ++iCount;
00087 }
00088 }
00089 }
00090 }
00091
00092
00093 mds::tSize FullSize = m_CellSize + 2 * m_CellMargin;
00094
00095
00096
00097 #ifdef LOG_VERTICES
00098 std::stringstream ss;
00099 ss << "points = [";
00100 #endif //LOG_VERTICES
00101
00102
00103
00104
00105
00106 for( k = m_CellMargin; (k + m_CellMargin) < (ZSize - m_CellMargin); k += FullSize )
00107 {
00108 mds::tSize MinK = k + m_CellMargin;
00109
00110 mds::tSize MaxK = mds::math::getMin(k + FullSize, ZSize);
00111
00112 for( j = m_CellMargin; (j + m_CellMargin) < (YSize - m_CellMargin); j += FullSize )
00113 {
00114 mds::tSize MinJ = j + m_CellMargin;
00115
00116 mds::tSize MaxJ = mds::math::getMin(j + FullSize, YSize);
00117
00118 for( i = m_CellMargin; (i + m_CellMargin) < (XSize - m_CellMargin); i += FullSize )
00119 {
00120 mds::tSize MinI = i + m_CellMargin;
00121
00122 mds::tSize MaxI = mds::math::getMin(i + FullSize, XSize);
00123
00124
00125 mds::tSize MaxX = MinI;
00126 mds::tSize MaxY = MinJ;
00127 mds::tSize MaxZ = MinK;
00128 double dMaxValue = pVolume->get(MaxX, MaxY, MaxZ);
00129
00130 for( mds::tSize wk = MinK; wk < MaxK; ++wk )
00131 {
00132 for( mds::tSize wj = MinJ; wj < MaxJ; ++wj )
00133 {
00134 for( mds::tSize wi = MinI; wi < MaxI; ++wi )
00135 {
00136 double dValue = pVolume->get(wi, wj, wk);
00137 if( dValue > dMaxValue )
00138 {
00139 dMaxValue = dValue;
00140 MaxX = wi;
00141 MaxY = wj;
00142 MaxZ = wk;
00143 }
00144 }
00145 }
00146 }
00147
00148
00149 if( dMaxValue > dThreshold )
00150 {
00151 pScale->conv2Real(MaxX, MaxY, MaxZ, &Point);
00152 vctl::MCVertex *pVertex = pVertices->New(Point);
00153 pVertex->SetValue((dMaxValue - dMin) * dTemp);
00154 ++iCount;
00155 #ifdef LOG_VERTICES
00156 ss << MaxX << ',' << MaxY << ',' << MaxZ << "; ";
00157 #endif //LOG_VERTICES
00158 }
00159 }
00160 }
00161 }
00162
00163
00164 for( j = m_CellMargin; (j + m_CellMargin) < (YSize - m_CellMargin); j += FullSize )
00165 {
00166 mds::tSize MinJ = j + m_CellMargin;
00167
00168 mds::tSize MaxJ = mds::math::getMin(j + FullSize, YSize);
00169
00170 for( i = m_CellMargin; (i + m_CellMargin) < (XSize - m_CellMargin); i += FullSize )
00171 {
00172 mds::tSize MinI = i + m_CellMargin;
00173
00174 mds::tSize MaxI = mds::math::getMin(i + FullSize, XSize);
00175
00176
00177 mds::tSize Max1X = MinI;
00178 mds::tSize Max1Y = MinJ;
00179 mds::tSize Max2X = MinI;
00180 mds::tSize Max2Y = MinJ;
00181 double dMax1Value = pVolume->get(Max1X, Max1Y, 0);
00182 double dMax2Value = pVolume->get(Max2X, Max2Y, ZSize - 1);
00183
00184 for( mds::tSize wj = MinJ; wj < MaxJ; ++wj )
00185 {
00186 for( mds::tSize wi = MinI; wi < MaxI; ++wi )
00187 {
00188 double dValue = pVolume->get(wi, wj, 0);
00189 if( dValue > dMax1Value )
00190 {
00191 dMax1Value = dValue;
00192 Max1X = wi;
00193 Max1Y = wj;
00194 }
00195 dValue = pVolume->get(wi, wj, ZSize - 1);
00196 if( dValue > dMax2Value )
00197 {
00198 dMax2Value = dValue;
00199 Max2X = wi;
00200 Max2Y = wj;
00201 }
00202 }
00203 }
00204
00205
00206 if( dMax1Value > dThreshold )
00207 {
00208 pScale->conv2Real(Max1X, Max1Y, 0, &Point);
00209 vctl::MCVertex *pVertex = pVertices->New(Point);
00210 pVertex->SetValue((dMax1Value - dMin) * dTemp);
00211 ++iCount;
00212 #ifdef LOG_VERTICES
00213 ss << Max1X << ',' << Max1Y << ',' << 0 << "; ";
00214 #endif //LOG_VERTICES
00215 }
00216 if( dMax2Value > dThreshold )
00217 {
00218 pScale->conv2Real(Max2X, Max2Y, ZSize - 1, &Point);
00219 vctl::MCVertex *pVertex = pVertices->New(Point);
00220 pVertex->SetValue((dMax2Value - dMin) * dTemp);
00221 ++iCount;
00222 #ifdef LOG_VERTICES
00223 ss << Max2X << ',' << Max2Y << ',' << ZSize - 1 << "; ";
00224 #endif //LOG_VERTICES
00225 }
00226 }
00227 }
00228
00229 for( k = m_CellMargin; (k + m_CellMargin) < (ZSize - m_CellMargin); k += FullSize )
00230 {
00231 mds::tSize MinK = k + m_CellMargin;
00232
00233 mds::tSize MaxK = mds::math::getMin(k + FullSize, ZSize);
00234
00235 for( i = m_CellMargin; (i + m_CellMargin) < (XSize - m_CellMargin); i += FullSize )
00236 {
00237 mds::tSize MinI = i + m_CellMargin;
00238
00239 mds::tSize MaxI = mds::math::getMin(i + FullSize, XSize);
00240
00241
00242 mds::tSize Max1X = MinI;
00243 mds::tSize Max1Z = MinK;
00244 mds::tSize Max2X = MinI;
00245 mds::tSize Max2Z = MinK;
00246 double dMax1Value = pVolume->get(Max1X, 0, Max1Z);
00247 double dMax2Value = pVolume->get(Max2X, YSize - 1, Max2Z);
00248
00249 for( mds::tSize wk = MinK; wk < MaxK; ++wk )
00250 {
00251 for( mds::tSize wi = MinI; wi < MaxI; ++wi )
00252 {
00253 double dValue = pVolume->get(wi, 0, wk);
00254 if( dValue > dMax1Value )
00255 {
00256 dMax1Value = dValue;
00257 Max1X = wi;
00258 Max1Z = wk;
00259 }
00260 dValue = pVolume->get(wi, YSize - 1, wk);
00261 if( dValue > dMax2Value )
00262 {
00263 dMax2Value = dValue;
00264 Max2X = wi;
00265 Max2Z = wk;
00266 }
00267 }
00268 }
00269
00270
00271 if( dMax1Value > dThreshold )
00272 {
00273 pScale->conv2Real(Max1X, 0, Max1Z, &Point);
00274 vctl::MCVertex *pVertex = pVertices->New(Point);
00275 pVertex->SetValue((dMax1Value - dMin) * dTemp);
00276 ++iCount;
00277 #ifdef LOG_VERTICES
00278 ss << Max1X << ',' << 0 << ',' << Max1Z << "; ";
00279 #endif //LOG_VERTICES
00280 }
00281 if( dMax2Value > dThreshold )
00282 {
00283 pScale->conv2Real(Max2X, YSize - 1, Max2Z, &Point);
00284 vctl::MCVertex *pVertex = pVertices->New(Point);
00285 pVertex->SetValue((dMax2Value - dMin) * dTemp);
00286 ++iCount;
00287 #ifdef LOG_VERTICES
00288 ss << Max2X << ',' << YSize - 1 << ',' << Max2Z << "; ";
00289 #endif //LOG_VERTICES
00290 }
00291 }
00292 }
00293
00294 for( k = m_CellMargin; (k + m_CellMargin) < (ZSize - m_CellMargin); k += FullSize )
00295 {
00296 mds::tSize MinK = k + m_CellMargin;
00297
00298 mds::tSize MaxK = mds::math::getMin(k + FullSize, ZSize);
00299
00300 for( j = m_CellMargin; (j + m_CellMargin) < (YSize - m_CellMargin); j += FullSize )
00301 {
00302 mds::tSize MinJ = j + m_CellMargin;
00303
00304 mds::tSize MaxJ = mds::math::getMin(j + FullSize, YSize);
00305
00306
00307 mds::tSize Max1Y = MinJ;
00308 mds::tSize Max1Z = MinK;
00309 mds::tSize Max2Y = MinJ;
00310 mds::tSize Max2Z = MinK;
00311 double dMax1Value = pVolume->get(0, Max1Y, Max1Z);
00312 double dMax2Value = pVolume->get(XSize - 1, Max2Y, Max2Z);
00313
00314 for( mds::tSize wk = MinK; wk < MaxK; ++wk )
00315 {
00316 for( mds::tSize wj = MinJ; wj < MaxJ; ++wj )
00317 {
00318 double dValue = pVolume->get(0, wj, wk);
00319 if( dValue > dMax1Value )
00320 {
00321 dMax1Value = dValue;
00322 Max1Y = wj;
00323 Max1Z = wk;
00324 }
00325 dValue = pVolume->get(XSize - 1, wj, wk);
00326 if( dValue > dMax2Value )
00327 {
00328 dMax2Value = dValue;
00329 Max2Y = wj;
00330 Max2Z = wk;
00331 }
00332 }
00333 }
00334
00335
00336 if( dMax1Value > dThreshold )
00337 {
00338 pScale->conv2Real(0, Max1Y, Max1Z, &Point);
00339 vctl::MCVertex *pVertex = pVertices->New(Point);
00340 pVertex->SetValue((dMax1Value - dMin) * dTemp);
00341 ++iCount;
00342 #ifdef LOG_VERTICES
00343 ss << 0 << ',' << Max1Y << ',' << Max1Z << "; ";
00344 #endif //LOG_VERTICES
00345 }
00346 if( dMax2Value > dThreshold )
00347 {
00348 pScale->conv2Real(XSize - 1, Max2Y, Max2Z, &Point);
00349 vctl::MCVertex *pVertex = pVertices->New(Point);
00350 pVertex->SetValue((dMax2Value - dMin) * dTemp);
00351 ++iCount;
00352 #ifdef LOG_VERTICES
00353 ss << XSize - 1 << ',' << Max2Y << ',' << Max2Z << "; ";
00354 #endif //LOG_VERTICES
00355 }
00356 }
00357 }
00358
00359
00360 #ifdef LOG_VERTICES
00361 ss << "]';";
00362 MDS_LOG_NOTE(ss.str());
00363 #endif //LOG_VERTICES*/
00364
00365 MDS_LOG_NOTE(" Number of boundary vertices = " << iCount);
00366
00367
00368 return true;
00369 }
00370