C++ Image3D类代码示例

int main( )
{
   typedef ImageSelector < Z3i::Domain, unsigned char>::Type Image3D;
   typedef ImageSelector < Z2i::Domain, unsigned char>::Type Image2D;
   typedef DGtal::ConstImageAdapter<Image3D, Image2D::Domain, DGtal::functors::Projector<Z3i::Space>,
                                    Image3D::Value,  DGtal::functors::Identity >  SliceImageAdapter;
   

   DGtal::functors::Projector<Z2i::Space >  proj(2);

   // Importing a 3D image
   std::string filename = examplesPath + "samples/lobster.vol";
   Image3D image = VolReader<Image3D>::importVol( filename );
   DGtal::Z2i::Domain domain(proj(image.domain().lowerBound()),
			     proj(image.domain().upperBound()));
   DGtal::functors::Identity idV;

   trace.beginBlock ( "Example extract2DImagesFrom3D" );

   // Extracting 2D slices ... and export them in the pgm format.
   for (unsigned int i=0; i<30; i+=10){
     std::stringstream name;
     name << "lobsterSliceZ_"  << i << ".pgm";
     DGtal::functors::Projector<Z3i::Space> aSliceFunctor(i); aSliceFunctor.initAddOneDim(2);
     SliceImageAdapter sliceImageZ(image, domain, aSliceFunctor, idV);
     PGMWriter<SliceImageAdapter>::exportPGM(name.str(), sliceImageZ);
   }

   // trace.endBlock();
   return 0;
}

int main( int /*argc*/, char** /*argv*/ )
{
  //! [extract2DImagesFrom3DType]
   typedef ImageSelector < Z3i::Domain, unsigned char>::Type Image3D;
   typedef ImageSelector < Z2i::Domain, unsigned char>::Type Image2D;
   typedef DGtal::ConstImageAdapter<Image3D, Z2i::Domain, DGtal::Point2DEmbedderIn3D<DGtal::Z3i::Domain>,
   				   Image3D::Value,  DGtal::DefaultFunctor >  ImageAdapterExtractor;
   //! [extract2DImagesFrom3DType]

   //! [extract2DImagesFrom3DOrigin3D]
   DGtal::Z3i::Point origin(150, 150, 10);
   DGtal::Z3i::Point ptUpper1(220, 220, 10);
   DGtal::Z3i::Point ptUpper2(150, 150, 50);
   DGtal::Z2i::Domain domainImage2D (DGtal::Z2i::Point(0,0),
				     DGtal::Z2i::Point((ptUpper1-origin).norm(),
						       (ptUpper2-origin).norm()));
   //! [extract2DImagesFrom3DOrigin3D]
   DGtal::Z3i::Point ptCenter(175, 175, 20);
   DGtal::Z2i::Domain domainImage2D2 (DGtal::Z2i::Point(0,0),
                                      DGtal::Z2i::Point(IMAGE_PATCH_WIDTH, IMAGE_PATCH_WIDTH));


   // Importing a 3D image
   std::string filename = examplesPath + "samples/lobster.vol";
   Image3D image = VolReader<Image3D>::importVol( filename );
   DGtal::Z3i::Domain domainImage3D = image.domain();
   DGtal::DefaultFunctor idV;

   trace.beginBlock ( "Example extract2DImagesFrom3D" );

   // Extracting 2D images ... and export them in the pgm format.
   for (unsigned int i=0; i<30; i+=10){
     std::stringstream name;
     name << "lobsterExtracted_"  << i << ".pgm";
     std::stringstream name2;
     name2 << "lobsterExtracted_"  << i << "V2.pgm";

     //! [extract2DImagesFrom3DOExtract]
     DGtal::Point2DEmbedderIn3D<DGtal::Z3i::Domain >  embedder(domainImage3D, origin+DGtal::Z3i::Point(i,i,0),
							       ptUpper1+DGtal::Z3i::Point(i,i,0),
							       ptUpper2+DGtal::Z3i::Point(i,i,0));

     ImageAdapterExtractor extractedImage(image, domainImage2D, embedder, idV);
     //! [extract2DImagesFrom3DOExtract]

     //! [extract2DImagesFrom3DOExtract2]
     DGtal::Point2DEmbedderIn3D<DGtal::Z3i::Domain >  embedder2(domainImage3D, ptCenter+DGtal::Z3i::Point(i,i,0),
								DGtal::Z3i::RealPoint(1,-1,0),
								IMAGE_PATCH_WIDTH);
     ImageAdapterExtractor extractedImage2(image, domainImage2D2, embedder2, idV);
     //! [extract2DImagesFrom3DOExtract2]

     PGMWriter< ImageAdapterExtractor>::exportPGM(name.str(), extractedImage);
     PGMWriter< ImageAdapterExtractor>::exportPGM(name2.str(), extractedImage2);
   }

   // trace.endBlock();
   return 0;
}

int main( int argc, char** argv )
{

  typedef DGtal::ImageContainerBySTLVector<DGtal::Z3i::Domain,  unsigned char > Image3D;
  //! [ExampleViewer3D2DImagesExtractImagesNonSliceType]
  typedef DGtal::ConstImageAdapter<Image3D, Z2i::Domain, DGtal::functors::Point2DEmbedderIn3D<DGtal::Z3i::Domain>,
                                   Image3D::Value,  DGtal::functors::Identity >  ImageAdapterExtractor;

  //! [ExampleViewer3D2DImagesExtractImagesNonSliceType]
  QApplication application(argc,argv);
  typedef Viewer3D<> MyViewer;
  MyViewer viewer;
  viewer.show();
  std::string inputFilename = examplesPath + "samples/lobster.vol";
  Image3D imageVol = VolReader<Image3D>::importVol(inputFilename);
  DGtal::functors::Identity idV;



  //! [ExampleViewer3D2DImagesExtractImagesNonSliceParam]
  DGtal::Z3i::Point ptCenter(50, 62, 28);
  const int IMAGE_PATCH_WIDTH = 20;
  // Setting the image domain of the resulting image to be displayed in 3D:
  DGtal::Z2i::Domain domainImage2D (DGtal::Z2i::Point(0,0),
                                    DGtal::Z2i::Point(IMAGE_PATCH_WIDTH, IMAGE_PATCH_WIDTH));
  //! [ExampleViewer3D2DImagesExtractImagesNonSliceParam]



  unsigned int pos=0;
  for (double alpha = 0; alpha< 1.54; alpha+= 0.01){
    //! [ExampleViewer3D2DImagesExtractImagesNonSliceExtract]
    // Extracting images from 3D embeder
    DGtal::functors::Point2DEmbedderIn3D<DGtal::Z3i::Domain >  embedder(imageVol.domain(),
                                                                        ptCenter+DGtal::Z3i::Point(static_cast<int>(200.0*cos(alpha)),static_cast<int>(100.0*sin(alpha))),
                                                                        DGtal::Z3i::RealPoint(cos(alpha),sin(alpha),cos(2.0*alpha)),
                                                                        IMAGE_PATCH_WIDTH);
    ImageAdapterExtractor extractedImage(imageVol, domainImage2D, embedder, idV);
    //! [ExampleViewer3D2DImagesExtractImagesNonSliceExtract]

    //! [ExampleViewer3D2DImagesExtractImagesNonSliceDisplay]
    //Display image and update its position with embeder
    viewer << extractedImage;
    viewer << DGtal::UpdateImage3DEmbedding<Z3i::Space, Z3i::KSpace>(pos,
                                                                     embedder(Z2i::RealPoint(0,0)),
                                                                     embedder(Z2i::RealPoint(IMAGE_PATCH_WIDTH,0)),
                                                                     embedder(domainImage2D.upperBound()),
                                                                     embedder(Z2i::RealPoint(0, IMAGE_PATCH_WIDTH)));
    //! [ExampleViewer3D2DImagesExtractImagesNonSliceDisplay]
    pos++;
  }


  viewer << MyViewer::updateDisplay;


  return application.exec();
}

void
getVoxelsStats(const Image3D &imageA,  int aMin, int aMax, const Image3D &imageB, 
	       int bMin, int bMax, bool exportStatVoxels,  std::vector<Point> &vectPtBinA,  
	       std::vector<Point> &vectPtCompBinCompA,  std::vector<Point> &vectPtBnotInA, 
	       std::vector<Point> &vectPtnotBInA, bool precisionRecallFMean ){
  int numBinA = 0; // true positif with A as ref shape.
  int numCompBinCompA = 0; // true neg with A as ref shape.
  int numBnotInA = 0; // false pos with A as ref shape
  int numNotBinA = 0; // false neg with A as ref shape
  int numTotalInA = 0; // total pos in reference shape
  int numTotalInB = 0; // total pos in compared shape
  int numTotalinCompA = 0; //total in complement A 
  int numTotalinCompB = 0; //total in complement B

  for(Image3D::Domain::ConstIterator it = imageA.domain().begin(); it!=imageA.domain().end(); it++){
    // voxels in A
    if(imageA(*it) <= aMax && imageA(*it) >= aMin){
      numTotalInA++;
      //voxels in B 
      if(imageB(*it) <= bMax && imageB(*it) >= bMin){
	numTotalInB++;
	numBinA++;
	if(exportStatVoxels) vectPtBinA.push_back(*it);
      }else{
	numTotalinCompB++;
	numNotBinA++;
	if(exportStatVoxels) vectPtnotBInA.push_back(*it);
      }
      // voxels outside A
    }else{
      numTotalinCompA++;
      // voxels in B
      if(imageB(*it) <= bMax && imageB(*it) >= bMin){
	numBnotInA++;
	numTotalInB++;
	if(exportStatVoxels) vectPtBnotInA.push_back(*it);
      }else{
	numTotalinCompB++;
	numCompBinCompA++;
	if(exportStatVoxels) vectPtCompBinCompA.push_back(*it);
      }      
    }
  }

  std::cout << numBinA << " " << numCompBinCompA << " " << numBnotInA 
	    << " " << numNotBinA  << " " << numTotalInA << " " 
	    << numTotalInB << " " << numTotalinCompA << " " << numTotalinCompB;
  if(precisionRecallFMean){
    double precision = (double)numBinA/(numBinA + numBnotInA);
    double recall = (double)numBinA/(numBinA+numNotBinA);
    double fmean = (2.0*precision*recall)/(precision+recall);
    std::cout << " " << precision<<  " " << recall << " " << fmean ; 
  }
}

TextureID Renderer::addTexture3D(const char *fileName, unsigned int flags){
	Image3D *img = new Image3D();

	if (img->loadFromFile(fileName)){
		return addTexture3D(img, flags);
	} else {
		addToLog(String("Couldn't load \"") + fileName + "\"");

		delete img;
		return TEXTURE_NONE;
	}
}

void 
embeddMeshInVol(const Image3D &anImage, My3DViewer &aViewer, DGtal::Mesh<DGtal::Z3i::RealPoint> &aMesh, 
                DGtal::Z3i::RealPoint vectTranslationNonMesh=DGtal::Z3i::RealPoint(0,0,0) )
{
  unsigned int  numImageDisplayed = aViewer.getCurrentGLImageNumber ();
  for (unsigned int i =0 ; i < aMesh.nbFaces(); i++){
    DGtal::Mesh<Z3i::RealPoint>::MeshFace aFace = aMesh.getFace(i);
    Z3i::RealPoint point1 = aMesh.getVertex(aFace.at(0));
    Z3i::RealPoint point2 = aMesh.getVertex(aFace.at(1));
    Z3i::RealPoint point3 = aMesh.getVertex(aFace.at(2));
    Z3i::RealPoint point4 = aMesh.getVertex(aFace.at(3));
    DGtal::Z2i::Domain domainImage2D (DGtal::Z2i::RealPoint(0,0), 
                                      DGtal::Z2i::RealPoint((point2-point1).norm(), (point4-point1).norm() ));    
    
    DGtal::functors::Point2DEmbedderIn3D<DGtal::Z3i::Domain >  embedder(anImage.domain(), 
                                                              point1, point2, point4, DGtal::Z3i::Point(0,0,0));      
    ImageAdapterExtractor extractedImage(anImage, domainImage2D, embedder, idV);        
    aViewer << extractedImage;
    aViewer << DGtal::UpdateImage3DEmbedding<Z3i::Space, Z3i::KSpace>(numImageDisplayed, 
                                                                     point1+vectTranslationNonMesh,
                                                                     point2+vectTranslationNonMesh,
                                                                     point3+vectTranslationNonMesh,
                                                                     point4+vectTranslationNonMesh);
    numImageDisplayed++;    
  }
}

int main( int argc, char** argv )
{

  typedef ImageContainerBySTLVector<Z3i::Domain,  unsigned char > Image3D;
  QApplication application(argc,argv);
  typedef Viewer3D<> MyViewer ;
  MyViewer viewer;

  viewer.show();
  std::string inputFilename = examplesPath + "samples/lobster.vol";
  Image3D imageVol = GenericReader<Image3D>::import(inputFilename);

  Z3i::Point ptLow (100, 100, 20);
  Z3i::Point ptUpp (200, 200, 40);
  Z3i::Domain subDomain(ptLow, ptUpp);

  Z3i::Point ptLow2 (220, 50, 10);
  Z3i::Point ptUpp2 (260, 100, 20);
  Z3i::Domain subDomain2(ptLow2, ptUpp2);

  Image3D imageCrop(subDomain);
  Image3D imageCrop2(subDomain2);

  for(Z3i::Domain::ConstIterator it= imageVol.domain().begin(), itend = imageVol.domain().end(); it != itend; ++it){
    if(imageVol(*it)>140)
      viewer << *it;
    Z3i::Point pt = *it;
    if(pt[0]>=ptLow[0] && pt[1] >= ptLow[1] && pt[2] >= ptLow[2] &&
       pt[0]<=ptUpp[0] && pt[1] <= ptUpp[1] && pt[2] <= ptUpp[2]){
      imageCrop.setValue(*it, imageVol(*it));
    }

    if(pt[0]>=ptLow2[0] && pt[1] >= ptLow2[1] && pt[2] >= ptLow2[2] &&
       pt[0]<=ptUpp2[0] && pt[1] <= ptUpp2[1] && pt[2] <= ptUpp2[2]){
      imageCrop2.setValue(*it, imageVol(*it));
    }
  }
  viewer << imageCrop;
  viewer << SetMode3D(imageCrop.className(), "BoundingBox");
    //! [ExampleViewer3D3DImagesDisplayImagesColor]
  viewer << AddTextureImage3DWithFunctor<Image3D, hueFct, Z3i::Space, Z3i::KSpace> (imageCrop2, hueFct(), MyViewer::RGBMode);
  viewer << MyViewer::updateDisplay;
  //! [ExampleViewer3D3DImagesDisplayImagesColor]

  return application.exec();
}

    // ----------------------------------------------------------------------------
    //  Generates a opacity map from the transfer function specification
    // ----------------------------------------------------------------------------
    void mapOpacity(Image3D<float> & map, std::list<std::shared_ptr<Node>> transfer)
    {
        float samples = static_cast<float>(map.width()) - 1;
        auto buffer = map.data();

        memset(buffer, 0, map.size()*sizeof(float));

        auto iter = transfer.begin();
        auto curr = *iter; 
        iter++;

        while (iter != transfer.end())
        {
            auto next = *iter;

            size_t x1 = static_cast<size_t>(curr->density * samples);
            if (x1 >= map.width()) x1 = map.width()-1;
            size_t x2 = static_cast<size_t>(next->density * samples);
            if (x2 >= map.width()) x2 = map.width()-1;
            float  y1 = curr->material->opticalThickness;
            float  y2 = next->material->opticalThickness;

            for (size_t i = x1; i <= x2; i++)
            {
                float px = i / samples - curr->density;
                float py = next->density - curr->density;
                float part = low( high(px / py, 0.0f), 1.0f );
                buffer[i] = - logf( 1.f - (1.f - part) * y1 - part * y2 );
            }

            curr = next;
            iter++;
        }
    }

void CubeMapTextureGLTest::imageFull() {
    if(!Context::current()->isExtensionSupported<Extensions::GL::ARB::direct_state_access>())
        CORRADE_SKIP(Extensions::GL::ARB::direct_state_access::string() + std::string(" is not supported."));

    CubeMapTexture texture;
    texture.setStorage(1, TextureFormat::RGBA8, Vector2i{2, 2})
        .setSubImage(0, {}, ImageReference3D{ColorFormat::RGBA, ColorType::UnsignedByte, Vector3i{2, 2, 6}, DataFull});

    MAGNUM_VERIFY_NO_ERROR();

    Image3D image = texture.image(0, {ColorFormat::RGBA, ColorType::UnsignedByte});

    MAGNUM_VERIFY_NO_ERROR();

    CORRADE_COMPARE(image.size(), Vector3i(2, 2, 6));
    CORRADE_COMPARE_AS(
        Containers::ArrayReference<const UnsignedByte>(image.data<UnsignedByte>(), image.pixelSize()*image.size().product()),
        Containers::ArrayReference<const UnsignedByte>{DataFull}, TestSuite::Compare::Container);
}

void
getStatsFromDistanceMap(Statistic<double> & stats, const Image3D &imageA, int aMin, int aMax,
			const Image3D & imageB,  int bMin, int bMax, 
			bool statOnFalsePositiveOnly=false, Point *ptMax=0){

  // Get the digital set from ref image by computing the surface (use -1 and +1 since the interval of append function are open)
  Z3i::DigitalSet set3dRef (imageA.domain()); 
  SetFromImage<Z3i::DigitalSet>::append<Image3D>(set3dRef, imageA, aMin-1,aMax);  
  typedef functors::NotPointPredicate<Z3i::DigitalSet> NegPredicate;
  

  // Applying the distance transform on the digital surface of the set: 
  typedef  DistanceTransformation<Z3i::Space, NegPredicate, Z3i::L2Metric> DTL2;
  const NegPredicate aPredicate( set3dRef );
  DTL2 dtL2( imageA.domain(), aPredicate, Z3i::l2Metric );

  // Get the set of point of imageB: (use -1 and +1 since the interval of append function are open)
  Z3i::DigitalSet set3dComp (imageB.domain()); 
  SetFromImage<Z3i::DigitalSet>::append<Image3D>(set3dComp, imageB, bMin-1, bMax);

  unsigned int nbAdded=0;
  double maxDist=0;
  //Applying stats from the set to be compared (from imageB)
  for(Z3i::DigitalSet::ConstIterator it= set3dComp.begin();  it!= set3dComp.end(); ++it){
    if((!statOnFalsePositiveOnly) || (isDiff(imageA, aMin, aMax, imageB, bMin, bMax, *it))){
      DTL2::Value distance = dtL2(*it);   
      stats.addValue(distance);
      nbAdded++;
      if(maxDist<distance){
	maxDist=distance;
	if(ptMax!=0){
	  (*ptMax)[0]=(*it)[0]; 
	  (*ptMax)[1]=(*it)[1];
	  (*ptMax)[2]=(*it)[2];
	}
      }
    }
  }
  
  if(nbAdded==0)
    trace.error() << "No point added to statistics, will failed..." << endl;
}

void fillImage3D( Image3D & img, 
                  typename Image3D::Point low, 
                  typename Image3D::Point up, 
                  typename Image3D::Value value )
{
  typedef typename Image3D::Point Point;
  typedef typename Image3D::Integer Integer;
  for ( Integer z = low[ 2 ]; z <= up[ 2 ]; ++z )
    for ( Integer y = low[ 1 ]; y <= up[ 1 ]; ++y )
      for ( Integer x = low[ 0 ]; x <= up[ 0 ]; ++x )
        img.setValue( Point( x, y, z ), value );
}

/**
 * Example of a test. To be completed.
 *
 */
bool testDicomReader()
{
  //Default image selector = STLVector
  typedef ImageContainerBySTLVector<Z3i::Domain,  unsigned char > Image3D;

  std::string filename = testPath + "samples/dicomSample/1629.dcm";
  Image3D image = DicomReader< Image3D >::importDicom( filename );

  trace.info() << image <<endl;

  unsigned int nbVal=0, nbPos = 0;
  for ( Image3D::ConstIterator it=image.begin(), itend=image.end() ; it != itend ; ++it )
  {
    nbVal++;
	  if ( (*it) > 0 ) nbPos++;
  }

  trace.info() << "Number of points with (val>0) = " << nbVal << endl;

  return nbVal==2130048 && nbPos==296030;
}

    Event CommandQueue::enqueueUnmap(const Image3D& image, void *mappedPtr,
                                     const Epic::Core::Array<Event>& eventWaitList) const
    {
        Event ret;
        cl_int err = 0;
        Array<cl_event> events = eventArrayToCLEventArray(eventWaitList);

        err = clEnqueueUnmapMemObject(queueHandle(), image.imageHandle(), mappedPtr,
                                      events.count(), events.data(), ret.eventPointer());

        EPIC_OPENCL_CHECK_ERROR(err)

        return ret;
    }

    Event CommandQueue::enqueueWriteImage(const Image3D& image, bool blockingRead, const size_t origin[3], const size_t region[3],
                                          size_t rowPitch, size_t slicePitch, void *pointer,
                                          const Epic::Core::Array<Event>& eventWaitList) const
    {
        Event ret;
        cl_int err = 0;
        Array<cl_event> events = eventArrayToCLEventArray(eventWaitList);

        err = clEnqueueWriteImage(queueHandle(), image.imageHandle(), blockingRead, origin, region, rowPitch, slicePitch,
                                  pointer, events.count(), events.data(), ret.eventPointer());

        EPIC_OPENCL_CHECK_ERROR(err)

        return ret;
    }

    // ----------------------------------------------------------------------------
    //  Generates an emissive map from the tranfer specification
    // ----------------------------------------------------------------------------
    void mapEmissive(Image3D<Vector4f> & map, std::list<std::shared_ptr<Node>> transfer)
    {
        float samples = static_cast<float>(map.width()) - 1;
        auto buffer = map.data();

        memset(buffer, 0, map.size()*sizeof(Vector4f));

        auto iter = transfer.begin();
        auto curr = *iter; 
        iter++;
        
        while (iter != transfer.end())
        {
            auto next = *iter;

            size_t x1 = static_cast<size_t>(curr->density * samples);
            if (x1 >= map.width()) x1 = map.width()-1;
            size_t x2 = static_cast<size_t>(next->density * samples);
            if (x2 >= map.width()) x2 = map.width()-1;                
            Vector3f s1 = Vector3f(curr->material->emissive) / 255.0f * curr->material->emissiveStrength;
            Vector3f s2 = Vector3f(next->material->emissive) / 255.0f * curr->material->emissiveStrength;
            Vector4f y1(s1[0], s1[1], s1[2], 0.0f);
            Vector4f y2(s2[0], s2[1], s2[2], 0.0f);

            for (size_t i = x1; i <= x2; i++)
            {
                float px = i / samples - curr->density;
                float py = next->density - curr->density;
                float part = low( high(px / py, 0.0f), 1.0f );
                buffer[i] = y1*(1.f - part) + y2*part;
            }

            curr = next;
            iter++;
        }
    }