mdsTetraVertexDetector.hxx

Go to the documentation of this file.
00001 //==============================================================================
00015 //#define LOG_VERTICES
00016 
00017 
00018 //==============================================================================
00019 /*
00020  * Methods templates.
00021  */
00022 
00023 //template <>
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     // Get image size
00034     mds::tSize XSize = pVolume->getXSize();
00035     mds::tSize YSize = pVolume->getYSize();
00036     mds::tSize ZSize = pVolume->getZSize();
00037 
00038     // Estimate threshold
00039 //    double dThreshold = mds::img::getMean<double>(*pVolume) + m_dThreshold * sqrt(mds::img::getVariance<double>(*pVolume));
00040     double dThreshold = 0.0;
00041 
00042     // Minimal and maximal values
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     // Total number of points
00048     int iCount = 0;
00049 
00050     // Add all points
00051 /*    vctl::MCPoint3D Point;
00052     mds::tSize i, j, k;
00053     for( k = 0; k < ZSize; ++k )
00054     {
00055         for( j = 0; j < YSize; ++j )
00056         {
00057             for( i = 0; i < XSize; ++i )
00058             {
00059                 double dValue = pVolume->get(i, j, k);
00060                 if( dValue > dThreshold )
00061                 {
00062                     pScale->conv2Real(i, j, k, &Point);
00063                     vctl::MCVertex *pVertex = pVertices->New(Point);
00064                     pVertex->SetValue((dValue - dMin) * dTemp);
00065                     ++iCount;
00066                 }
00067             }
00068         }
00069     }*/
00070     
00071     // Add all corner
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     // Helper values
00093     mds::tSize FullSize = m_CellSize + 2 * m_CellMargin;
00094 //    mds::tSize TempSize = m_CellSize + m_CellMargin;
00095 
00096     // Helper string stream used for logging
00097 #ifdef LOG_VERTICES
00098     std::stringstream ss;
00099     ss << "points = [";
00100 #endif //LOG_VERTICES
00101 
00102     // Insert all significant edge points to the list of vertices
00103     // - Process input volume per blocks
00104 //    vctl::MCPoint3D Point;
00105 //    mds::tSize i, j, k;
00106     for( k = m_CellMargin; (k + m_CellMargin) < (ZSize - m_CellMargin); k += FullSize )
00107     {
00108         mds::tSize MinK = k + m_CellMargin;
00109 //        mds::tSize MaxK = mds::math::getMin(k + TempSize, ZSize);
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 //            mds::tSize MaxJ = mds::math::getMin(j + TempSize, YSize);
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 //                mds::tSize MaxI = mds::math::getMin(i + TempSize, XSize);
00122                 mds::tSize MaxI = mds::math::getMin(i + FullSize, XSize);
00123 
00124                 // Find maximum in the current block
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                 // Create a new vertex
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     // Insert all significant edge points on the front and rear volume side
00164     for( j = m_CellMargin; (j + m_CellMargin) < (YSize - m_CellMargin); j += FullSize )
00165     {
00166         mds::tSize MinJ = j + m_CellMargin;
00167 //        mds::tSize MaxJ = mds::math::getMin(j + TempSize, YSize);
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 //            mds::tSize MaxI = mds::math::getMin(i + TempSize, XSize);
00174             mds::tSize MaxI = mds::math::getMin(i + FullSize, XSize);
00175 
00176             // Find maximum in the current block
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             // Create a new vertex
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 //        mds::tSize MaxK = mds::math::getMin(k + TempSize, ZSize);
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 //            mds::tSize MaxI = mds::math::getMin(i + TempSize, XSize);
00239             mds::tSize MaxI = mds::math::getMin(i + FullSize, XSize);
00240 
00241             // Find maximum in the current block
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             // Create a new vertex
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 //        mds::tSize MaxK = mds::math::getMin(k + TempSize, ZSize);
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 //            mds::tSize MaxJ = mds::math::getMin(j + TempSize, YSize);
00304             mds::tSize MaxJ = mds::math::getMin(j + FullSize, YSize);
00305 
00306             // Find maximum in the current block
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             // Create a new vertex
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     // Logging
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     // O.K.
00368     return true;
00369 }
00370 

Generated on Thu Mar 11 10:35:44 2010 for MDSTk Extension Libraries by  doxygen 1.4.6-NO