diff --git a/Modules/IpPicSupport/mitkPicFileReader.cpp b/Modules/IpPicSupport/mitkPicFileReader.cpp index e841ad852d..3ba42ae077 100644 --- a/Modules/IpPicSupport/mitkPicFileReader.cpp +++ b/Modules/IpPicSupport/mitkPicFileReader.cpp @@ -1,373 +1,375 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkPicFileReader.h" #include "mitkPicHelper.h" #include +#include "mitkImageWriteAccessor.h" extern "C" { mitkIpPicDescriptor * MITKipPicGet( char *infile_name, mitkIpPicDescriptor *pic ); mitkIpPicDescriptor * MITKipPicGetTags( char *infile_name, mitkIpPicDescriptor *pic ); } void mitk::PicFileReader::GenerateOutputInformation() { Image::Pointer output = this->GetOutput(); if ((output->IsInitialized()) && (this->GetMTime() <= m_ReadHeaderTime.GetMTime())) return; itkDebugMacro(<<"Reading file for GenerateOutputInformation()" << m_FileName); // Check to see if we can read the file given the name or prefix // if ( m_FileName == "" && m_FilePrefix == "" ) { throw itk::ImageFileReaderException(__FILE__, __LINE__, "One of FileName or FilePrefix must be non-empty"); } if( m_FileName != "") { mitkIpPicDescriptor* header=mitkIpPicGetHeader(const_cast(m_FileName.c_str()), NULL); if ( !header ) { throw itk::ImageFileReaderException(__FILE__, __LINE__, "File could not be read."); } header=MITKipPicGetTags(const_cast(m_FileName.c_str()), header); int channels = 1; mitkIpPicTSV_t *tsv; if ( (tsv = mitkIpPicQueryTag( header, "SOURCE HEADER" )) != NULL) { if(tsv->n[0]>1e+06) { mitkIpPicTSV_t *tsvSH; tsvSH = mitkIpPicDelTag( header, "SOURCE HEADER" ); mitkIpPicFreeTag(tsvSH); } } if ( (tsv = mitkIpPicQueryTag( header, "ICON80x80" )) != NULL) { mitkIpPicTSV_t *tsvSH; tsvSH = mitkIpPicDelTag( header, "ICON80x80" ); mitkIpPicFreeTag(tsvSH); } if ( (tsv = mitkIpPicQueryTag( header, "VELOCITY" )) != NULL) { ++channels; mitkIpPicDelTag( header, "VELOCITY" ); } if( header == NULL || header->bpe == 0) { itk::ImageFileReaderException e(__FILE__, __LINE__); std::ostringstream msg; msg << " Could not read file " << m_FileName.c_str(); e.SetDescription(msg.str().c_str()); throw e; return; } // if pic image only 2D, the n[2] value is not initialized unsigned int slices = 1; if( header->dim == 2 ) header->n[2] = slices; // First initialize the geometry of the output image by the pic-header SlicedGeometry3D::Pointer slicedGeometry = mitk::SlicedGeometry3D::New(); PicHelper::InitializeEvenlySpaced(header, header->n[2], slicedGeometry); // if pic image only 3D, the n[3] value is not initialized unsigned int timesteps = 1; if( header->dim > 3 ) timesteps = header->n[3]; TimeSlicedGeometry::Pointer timeSliceGeometry = TimeSlicedGeometry::New(); timeSliceGeometry->InitializeEvenlyTimed(slicedGeometry, timesteps); timeSliceGeometry->ImageGeometryOn(); // Cast the pic descriptor to ImageDescriptor and initialize the output output->Initialize( CastToImageDescriptor(header)); output->SetGeometry( timeSliceGeometry ); mitkIpPicFree ( header ); } else { int numberOfImages=0; m_StartFileIndex=0; mitkIpPicDescriptor* header=NULL; char fullName[1024]; while(m_StartFileIndex<10) { sprintf(fullName, m_FilePattern.c_str(), m_FilePrefix.c_str(), m_StartFileIndex+numberOfImages); FILE * f=fopen(fullName,"r"); if(f==NULL) { //already found an image? if(numberOfImages>0) break; //no? let's increase start ++m_StartFileIndex; } else { fclose(f); //only open the header of the first file found, //@warning what to do when images do not have the same size?? if(header==NULL) { header=mitkIpPicGetHeader(fullName, NULL); header=MITKipPicGetTags(fullName, header); } ++numberOfImages; } } printf("\n numberofimages %d\n",numberOfImages); if(numberOfImages==0) { itk::ImageFileReaderException e(__FILE__, __LINE__); std::ostringstream msg; msg << "no images found"; e.SetDescription(msg.str().c_str()); throw e; return; } //@FIXME: was ist, wenn die Bilder nicht alle gleich gross sind? if(numberOfImages>1) { printf("\n numberofimages %d > 1\n",numberOfImages); header->dim=3; header->n[2]=numberOfImages; } printf(" \ninitialisize output\n"); output->Initialize( CastToImageDescriptor(header) ); mitkIpPicFree ( header ); } m_ReadHeaderTime.Modified(); } void mitk::PicFileReader::ConvertHandedness(mitkIpPicDescriptor* pic) { //left to right handed conversion if(pic->dim>=3) { mitkIpPicDescriptor* slice=mitkIpPicCopyHeader(pic, NULL); slice->dim=2; size_t size=_mitkIpPicSize(slice); slice->data=malloc(size); size_t v, volumes = (pic->dim>3? pic->n[3] : 1); size_t volume_size = size*pic->n[2]; for(v=0; vdata; unsigned char *p_last=(unsigned char *)pic->data; p_first+=v*volume_size; p_last+=size*(pic->n[2]-1)+v*volume_size; size_t i, smid=pic->n[2]/2; for(i=0; idata, p_last, size); memcpy(p_last, p_first, size); memcpy(p_first, slice->data, size); } } mitkIpPicFree(slice); } } void mitk::PicFileReader::GenerateData() { Image::Pointer output = this->GetOutput(); // Check to see if we can read the file given the name or prefix // if ( m_FileName == "" && m_FilePrefix == "" ) { throw itk::ImageFileReaderException(__FILE__, __LINE__, "One of FileName or FilePrefix must be non-empty"); } if( m_FileName != "") { mitkIpPicDescriptor* outputPic = mitkIpPicNew(); - outputPic = CastToIpPicDescriptor(output, outputPic); + mitk::ImageWriteAccessor imageAccess(output); + outputPic = CastToIpPicDescriptor(output, &imageAccess, outputPic); mitkIpPicDescriptor* pic=MITKipPicGet(const_cast(m_FileName.c_str()), outputPic); // comes upside-down (in MITK coordinates) from PIC file ConvertHandedness(pic); mitkIpPicTSV_t *tsv; if ( (tsv = mitkIpPicQueryTag( pic, "SOURCE HEADER" )) != NULL) { if(tsv->n[0]>1e+06) { mitkIpPicTSV_t *tsvSH; tsvSH = mitkIpPicDelTag( pic, "SOURCE HEADER" ); mitkIpPicFreeTag(tsvSH); } } if ( (tsv = mitkIpPicQueryTag( pic, "ICON80x80" )) != NULL) { mitkIpPicTSV_t *tsvSH; tsvSH = mitkIpPicDelTag( pic, "ICON80x80" ); mitkIpPicFreeTag(tsvSH); } if ( (tsv = mitkIpPicQueryTag( pic, "VELOCITY" )) != NULL) { mitkIpPicDescriptor* header = mitkIpPicCopyHeader(pic, NULL); header->data = tsv->value; ConvertHandedness(header); output->SetChannel(header->data, 1); header->data = NULL; mitkIpPicFree(header); mitkIpPicDelTag( pic, "VELOCITY" ); } //slice-wise reading //currently much too slow. //else //{ // int sstart, smax; // int tstart, tmax; // sstart=output->GetRequestedRegion().GetIndex(2); // smax=sstart+output->GetRequestedRegion().GetSize(2); // tstart=output->GetRequestedRegion().GetIndex(3); // tmax=tstart+output->GetRequestedRegion().GetSize(3); // int s,t; // for(s=sstart; s(m_FileName.c_str()), NULL, t*smax+s+1); // output->SetPicSlice(pic,s,t); // } // } //} } else { int position; mitkIpPicDescriptor* pic=NULL; int zDim=(output->GetDimension()>2?output->GetDimensions()[2]:1); printf("\n zdim is %u \n",zDim); for (position = 0; position < zDim; ++position) { char fullName[1024]; sprintf(fullName, m_FilePattern.c_str(), m_FilePrefix.c_str(), m_StartFileIndex+position); pic=MITKipPicGet(fullName, pic); if(pic==NULL) { itkDebugMacro("Pic file '" << fullName << "' does not exist."); } /* FIXME else if(output->SetPicSlice(pic, position)==false) { itkDebugMacro("Image '" << fullName << "' could not be added to Image."); }*/ } if(pic!=NULL) mitkIpPicFree(pic); } } void mitk::PicFileReader::EnlargeOutputRequestedRegion(itk::DataObject *output) { output->SetRequestedRegionToLargestPossibleRegion(); } bool mitk::PicFileReader::CanReadFile(const std::string filename, const std::string filePrefix, const std::string filePattern) { // First check the extension if( filename == "" ) { //MITK_INFO<<"No filename specified."<SetNumberOfRequiredInputs( 1 ); } mitk::PicFileWriter::~PicFileWriter() { } void mitk::PicFileWriter::GenerateData() { if ( m_FileName == "" ) { itkWarningMacro( << "Sorry, filename has not been set!" ); return ; } std::ofstream testfilehandle( m_FileName.c_str(), std::ios::out); if (!testfilehandle.good()) { testfilehandle.close(); itkExceptionMacro(<<"File location '" << m_FileName << "' not writeable"); } else { testfilehandle.close(); } Image::Pointer input = const_cast(this->GetInput()); if ( input.IsNull() ) { itkExceptionMacro(<< "Nothing to write: Input is NULL." ); } if( input->GetNumberOfChannels() > 1) { std::cout << "Multiple channel. Cannot write. Will throw..."; itkExceptionMacro(<< "The PicFileWriter does not support multiple channel data. Nothing will be written." ); } mitkIpPicDescriptor * picImage = mitkIpPicNew(); - picImage = CastToIpPicDescriptor(input, picImage); + mitk::ImageWriteAccessor imageAccess(input); + picImage = CastToIpPicDescriptor(input, &imageAccess, picImage); SlicedGeometry3D* slicedGeometry = input->GetSlicedGeometry(); if (slicedGeometry != NULL) { //set tag "REAL PIXEL SIZE" const Vector3D & spacing = slicedGeometry->GetSpacing(); mitkIpPicTSV_t *pixelSizeTag; pixelSizeTag = mitkIpPicQueryTag( picImage, "REAL PIXEL SIZE" ); if (!pixelSizeTag) { pixelSizeTag = (mitkIpPicTSV_t *) malloc( sizeof(mitkIpPicTSV_t) ); pixelSizeTag->type = mitkIpPicFloat; pixelSizeTag->bpe = 32; strcpy(pixelSizeTag->tag, "REAL PIXEL SIZE"); pixelSizeTag->dim = 1; pixelSizeTag->n[0] = 3; pixelSizeTag->value = malloc( sizeof(float) * 3 ); mitkIpPicAddTag (picImage, pixelSizeTag); } ((float*)pixelSizeTag->value)[0] = spacing[0]; ((float*)pixelSizeTag->value)[1] = spacing[1]; ((float*)pixelSizeTag->value)[2] = spacing[2]; //set tag "ISG" //ISG == offset/origin transformation matrix(matrix) spancings //ISG == offset0 offset1 offset2 spalte0_0 spalte0_1 spalte0_2 spalte1_0 spalte1_1 spalte1_2 spalte2_0 spalte2_1 spalte2_2 spacing0 spacing1 spacing2 mitkIpPicTSV_t *geometryTag; geometryTag = mitkIpPicQueryTag( picImage, "ISG" ); if (!geometryTag) { geometryTag = (mitkIpPicTSV_t *) malloc( sizeof(mitkIpPicTSV_t) ); geometryTag->type = mitkIpPicFloat; geometryTag->bpe = 32; strcpy(geometryTag->tag, "ISG"); geometryTag->dim = 2; geometryTag->n[0] = 3; geometryTag->n[1] = 4; geometryTag->value = malloc( sizeof(float) * 3 * 4 ); mitkIpPicAddTag (picImage, geometryTag); } const AffineTransform3D::OffsetType& offset = slicedGeometry->GetIndexToWorldTransform()->GetOffset(); ((float*)geometryTag->value)[0] = offset[0]; ((float*)geometryTag->value)[1] = offset[1]; ((float*)geometryTag->value)[2] = offset[2]; const AffineTransform3D::MatrixType& matrix = slicedGeometry->GetIndexToWorldTransform()->GetMatrix(); const AffineTransform3D::MatrixType::ValueType* row0 = matrix[0]; const AffineTransform3D::MatrixType::ValueType* row1 = matrix[1]; const AffineTransform3D::MatrixType::ValueType* row2 = matrix[2]; Vector3D v; FillVector3D(v, row0[0], row1[0], row2[0]); v.Normalize(); ((float*)geometryTag->value)[3] = v[0]; ((float*)geometryTag->value)[4] = v[1]; ((float*)geometryTag->value)[5] = v[2]; FillVector3D(v, row0[1], row1[1], row2[1]); v.Normalize(); ((float*)geometryTag->value)[6] = v[0]; ((float*)geometryTag->value)[7] = v[1]; ((float*)geometryTag->value)[8] = v[2]; ((float*)geometryTag->value)[9] = spacing[0]; ((float*)geometryTag->value)[10] = spacing[1]; ((float*)geometryTag->value)[11] = spacing[2]; } PicFileReader::ConvertHandedness(picImage); // flip upside-down in MITK coordinates // Following line added to detect write errors. If saving .pic files from the plugin is broken again, // please report a bug, don't just remove this line! int ret = MITKIpPicPut((char*)(m_FileName.c_str()), picImage); if (ret != 0) { PicFileReader::ConvertHandedness(picImage); // flip back from upside-down state throw std::ios_base::failure("Error during .pic file writing in "__FILE__); } PicFileReader::ConvertHandedness(picImage); // flip back from upside-down state } void mitk::PicFileWriter::SetInputImage( Image* image ) { this->ProcessObject::SetNthInput( 0, image ); } const mitk::Image* mitk::PicFileWriter::GetInput() { if ( this->GetNumberOfInputs() < 1 ) { MITK_ERROR << "No input image present."; return NULL; } else { return static_cast< const Image * >( this->ProcessObject::GetInput( 0 ) ); } } int mitk::PicFileWriter::MITKIpPicPut( char *outfile_name, mitkIpPicDescriptor *pic ) { FILE* outfile; mitkIpUInt4_t len; mitkIpUInt4_t tags_len; if( pic->info->write_protect ) { fprintf( stderr, "mitkIpPicPut: sorry, can't write (missing tags !!!)\n" ); //return( -1 ); } if( mitkIpPicEncryptionType(pic) != ' ' ) { fprintf( stderr, "mitkIpPicPut: warning: was encrypted !!!\n" ); } if( outfile_name == NULL ) outfile = stdout; else if( strcmp(outfile_name, "stdout") == 0 ) outfile = stdout; else { mitkIpPicRemoveFile( outfile_name ); // Removed due to linker problems when compiling // an mitk chili plugin using msvc: there appear // unresolved external symbol errors to function // _ipPicGetWriteCompression() /* if( mitkIpPicGetWriteCompression() ) { char buff[1024]; sprintf( buff, "%s.gz", outfile_name ); outfile = (FILE*) mitkIpPicFOpen( buff, "wb" ); // cast to prevent warning. } else */ outfile = fopen( outfile_name, "wb" ); } if( outfile == NULL ) { fprintf( stderr, "mitkIpPicPut: sorry, error opening outfile\n" ); return( -1 ); } tags_len = _mitkIpPicTagsSize( pic->info->tags_head ); len = tags_len + 3 * sizeof(mitkIpUInt4_t) + pic->dim * sizeof(mitkIpUInt4_t); /* write oufile */ if( mitkIpPicEncryptionType(pic) == ' ' ) mitkIpPicFWrite( mitkIpPicVERSION, 1, sizeof(mitkIpPicTag_t), outfile ); else mitkIpPicFWrite( pic->info->version, 1, sizeof(mitkIpPicTag_t), outfile ); mitkIpPicFWriteLE( &len, sizeof(mitkIpUInt4_t), 1, outfile ); mitkIpPicFWriteLE( &(pic->type), sizeof(mitkIpUInt4_t), 1, outfile ); mitkIpPicFWriteLE( &(pic->bpe), sizeof(mitkIpUInt4_t), 1, outfile ); mitkIpPicFWriteLE( &(pic->dim), sizeof(mitkIpUInt4_t), 1, outfile ); mitkIpPicFWriteLE( pic->n, sizeof(mitkIpUInt4_t), pic->dim, outfile ); _mitkIpPicWriteTags( pic->info->tags_head, outfile, mitkIpPicEncryptionType(pic) ); // Removed due to linker problems when compiling // an mitk chili plugin using msvc: there appear // unresolved external symbol errors to function // _ipPicGetWriteCompression() /* if( mitkIpPicGetWriteCompression() ) pic->info->pixel_start_in_file = mitkIpPicFTell( outfile ); else */ pic->info->pixel_start_in_file = ftell( outfile ); if( pic->data ) { size_t number_of_elements = _mitkIpPicElements(pic); size_t bytes_per_element = pic->bpe / 8; size_t number_of_bytes = number_of_elements * bytes_per_element; size_t block_size = 1024*1024; /* Use 1 MB blocks. Make sure that block size is smaller than 2^31 */ size_t number_of_blocks = number_of_bytes / block_size; size_t remaining_bytes = number_of_bytes % block_size; size_t bytes_written = 0; size_t block_nr = 0; mitkIpUInt1_t* data = (mitkIpUInt1_t*) pic->data; assert( data != NULL ); if( pic->type == mitkIpPicNonUniform ) { for ( block_nr = 0 ; block_nr < number_of_blocks ; ++block_nr ) bytes_written += mitkIpPicFWrite( data + ( block_nr * block_size ), 1, block_size, outfile ); bytes_written += mitkIpPicFWrite( data + ( number_of_blocks * block_size ), 1, remaining_bytes, outfile ); } else { for ( block_nr = 0 ; block_nr < number_of_blocks ; ++block_nr ) bytes_written += mitkIpPicFWriteLE( data + ( block_nr * block_size ), 1, block_size, outfile ); bytes_written += mitkIpPicFWriteLE( data + ( number_of_blocks * block_size ), 1, remaining_bytes, outfile ); } if ( bytes_written != number_of_bytes ) { fprintf( stderr, "Error while writing (ferror indicates %u), only %u bytes were written! Eof indicator is %u.\n", ferror(outfile), ( (unsigned int) ( bytes_written ) ), feof(outfile) ); fclose( outfile ); return( -1 ); } } if( outfile != stdout ) { // Removed due to linker problems when compiling // an mitk chili plugin using msvc: there appear // unresolved external symbol errors to function // _ipPicGetWriteCompression() /* if( mitkIpPicGetWriteCompression() ) mitkIpPicFClose( outfile ); else */ fclose( outfile ); } return( 0 ); } std::vector mitk::PicFileWriter::GetPossibleFileExtensions() { std::vector possibleFileExtensions; possibleFileExtensions.push_back(".pic"); return possibleFileExtensions; } diff --git a/Modules/MitkExt/Algorithms/mitkCylindricToCartesianFilter.cpp b/Modules/MitkExt/Algorithms/mitkCylindricToCartesianFilter.cpp index 014e390b20..78defa2e6b 100644 --- a/Modules/MitkExt/Algorithms/mitkCylindricToCartesianFilter.cpp +++ b/Modules/MitkExt/Algorithms/mitkCylindricToCartesianFilter.cpp @@ -1,497 +1,498 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkCylindricToCartesianFilter.h" #include "mitkImageTimeSelector.h" #include "mitkSlicedGeometry3D.h" #include "mitkPlaneGeometry.h" #include "mitkProperties.h" #include "mitkLegacyAdaptors.h" #include template void _transform(mitkIpPicDescriptor *pic, mitkIpPicDescriptor *dest, float _outsideValue, float *fr, float *fphi, float *fz, short *rt, unsigned int *phit, unsigned int *zt, mitkIpPicDescriptor *coneCutOff_pic) //...t=truncated { T outsideValue = static_cast(_outsideValue); register float f, ft, f0, f1, f2, f3; mitkIpInt2_t ox_size; mitkIpInt2_t nx_size, ny_size, nz_size; int oxy_size, nxy_size; T* orig, *dp, *dest_start; mitkIpInt2_t* coneCutOff=(mitkIpInt2_t*)coneCutOff_pic->data; orig=(T*)pic->data; ox_size=pic->n[0]; oxy_size=ox_size*pic->n[1]; nx_size=dest->n[0]; ny_size=dest->n[1]; nxy_size=nx_size*ny_size; nz_size=dest->n[2]; /*nx_size=360; ny_size=360; nxy_size=nx_size*ny_size; nz_size=256;*/ dest_start=dp=((T*)dest->data)+nxy_size*(nz_size-1); mitkIpInt2_t y; // int size=_mitkIpPicElements(pic); register mitkIpInt2_t x,z; for(y=0;y=0) { x_start=0; x_end=nx_size; } else { x_start=-r0plusphi0; x_end=nx_size+r0plusphi0; for(z=0;ztype=mitkIpPicInt; rt_pic->bpe=16; rt_pic->dim=2; rt_pic->n[0]=rt_pic->n[1]=new_xsize; rt_pic->data=malloc(_mitkIpPicSize(rt_pic)); phit_pic=mitkIpPicNew(); phit_pic->type=mitkIpPicUInt; phit_pic->bpe=32; phit_pic->dim=2; phit_pic->n[0]=phit_pic->n[1]=new_xsize; phit_pic->data=malloc(_mitkIpPicSize(phit_pic)); fr_pic=mitkIpPicNew(); fr_pic->type=mitkIpPicFloat; fr_pic->bpe=32; fr_pic->dim=2; fr_pic->n[0]=fr_pic->n[1]=new_xsize; fr_pic->data=malloc(_mitkIpPicSize(fr_pic)); fphi_pic=mitkIpPicNew(); fphi_pic->type=mitkIpPicFloat; fphi_pic->bpe=32; fphi_pic->dim=2; fphi_pic->n[0]=fphi_pic->n[1]=new_xsize; fphi_pic->data=malloc(_mitkIpPicSize(fphi_pic)); mitkIpInt2_t *rtp=(mitkIpInt2_t*)rt_pic->data, *rt_xzero, rt, phit; mitkIpUInt4_t *phitp=(mitkIpUInt4_t*)phit_pic->data; mitkIpFloat4_t *fr=(mitkIpFloat4_t *)fr_pic->data; mitkIpFloat4_t *fphi=(mitkIpFloat4_t *)fphi_pic->data; mitkIpFloat4_t r, phi, scale=(double)orig_xsize/(double)new_xsize; int x,y,xy0,xy0_orig, oxy_size, new_zsize; oxy_size=orig_xsize*orig_ysize; xy0=(int)(((double)new_xsize)/2+0.5); xy0_orig=(int)(((double)orig_xsize)/2+0.5); new_zsize=(int)(orig_ysize/scale); // \bug y compared to x for(y=0;yanfangen bei -rt+1!*/ // if((x>=-rt) && (xxy0?1.0:-1.0)*scale+xy0_orig; else r=r*(x>xy0?-1.0:1.0)*scale+xy0_orig; rt=(mitkIpInt2_t)r; int xtmp=x; if(x>xy0) xtmp=new_xsize-x; if(rt<0) { r=rt=0; if(xtmp>-*rt_xzero) *rt_xzero=-xtmp; *fr=0; } else if(rt>orig_xsize-1) { r=rt=orig_xsize-1; if(xtmp>-*rt_xzero) *rt_xzero=-xtmp; *fr=0; } else *fr=r-rt; if(*fr<0) *fr=0; } // else // *fr=0; phi=orig_zsize-(yq==0?1:-atan((float)xq/yq)/M_PI+0.5)*orig_zsize; phit=(mitkIpUInt4_t)phi; *fphi=phi-phit; *rtp=rt; *phitp=phit*oxy_size; } } zt=(unsigned int *)malloc(sizeof(unsigned int)*new_zsize); fz=(float *)malloc(sizeof(float)*new_zsize); float *fzp=fz; unsigned int *ztp=zt; int z; float z_step=orig_ysize/(orig_ysize*((float)new_xsize)/orig_xsize); for(z=0;ztype=mitkIpPicInt; coneCutOff_pic->bpe=16; coneCutOff_pic->dim=2; coneCutOff_pic->n[0]=coneCutOff_pic->n[1]=rt_pic->n[0]; coneCutOff_pic->data=malloc(_mitkIpPicSize(coneCutOff_pic)); int i, size=_mitkIpPicElements(rt_pic); mitkIpInt2_t *rt, *ccop, ohx_size, nz_size; mitkIpFloat4_t *fr; a*=(float)rt_pic->n[0]/orig_xsize; b*=(float)rt_pic->n[0]/orig_xsize; ohx_size=orig_xsize/2; nz_size=orig_ysize*rt_pic->n[0]/orig_xsize; rt=(mitkIpInt2_t *)rt_pic->data; fr=(mitkIpFloat4_t*)fr_pic->data; ccop=(mitkIpInt2_t *)coneCutOff_pic->data; for(i=0; i=nz_size) cco=nz_size; *ccop=cco; } } void mitk::CylindricToCartesianFilter::GenerateOutputInformation() { mitk::Image::Pointer output = this->GetOutput(); if ((output->IsInitialized()) && (output->GetPipelineMTime() <= m_TimeOfHeaderInitialization.GetMTime())) return; mitk::Image::ConstPointer input = this->GetInput(); itkDebugMacro(<<"GenerateOutputInformation()"); unsigned int i, *tmpDimensions=new unsigned int[std::max(3u,input->GetDimension())]; tmpDimensions[0]=m_TargetXSize; if(tmpDimensions[0]==0) tmpDimensions[0] = input->GetDimension(0); float scale=((float)tmpDimensions[0])/input->GetDimension(0); tmpDimensions[1] = tmpDimensions[0]; tmpDimensions[2] = (unsigned int)(scale*input->GetDimension(1)); for(i=3;iGetDimension();++i) tmpDimensions[i]=input->GetDimension(i); output->Initialize(input->GetPixelType(), input->GetDimension(), tmpDimensions, input->GetNumberOfChannels()); // initialize the spacing of the output Vector3D spacing = input->GetSlicedGeometry()->GetSpacing(); if(input->GetDimension()>=2) spacing[2]=spacing[1]; else spacing[2] = 1.0; spacing[1] = spacing[0]; spacing *= 1.0/scale; output->GetSlicedGeometry()->SetSpacing(spacing); mitk::Point3iProperty::Pointer pointProp; pointProp = dynamic_cast(input->GetProperty("ORIGIN").GetPointer()); if (pointProp.IsNotNull() ) { itk::Point tp = pointProp->GetValue(); tp[2] = (int)(tmpDimensions[2]-tp[1] * scale-1); tp[0] = tmpDimensions[0]/2; tp[1] = tmpDimensions[0]/2; mitk::Point3iProperty::Pointer pointProp = mitk::Point3iProperty::New(tp); output->SetProperty("ORIGIN", pointProp); } delete [] tmpDimensions; //output->GetSlicedGeometry()->SetGeometry2D(mitk::Image::BuildStandardPlaneGeometry2D(output->GetSlicedGeometry(), tmpDimensions).GetPointer(), 0); //set the timebounds - after SetGeometry2D, so that the already created PlaneGeometry will also receive this timebounds. //@fixme!!! will not work for not evenly timed data! output->GetSlicedGeometry()->SetTimeBounds(input->GetSlicedGeometry()->GetTimeBounds()); output->GetTimeSlicedGeometry()->InitializeEvenlyTimed(output->GetSlicedGeometry(), output->GetTimeSlicedGeometry()->GetTimeSteps()); output->SetPropertyList(input->GetPropertyList()->Clone()); m_TimeOfHeaderInitialization.Modified(); } void mitk::CylindricToCartesianFilter::GenerateData() { mitk::Image::ConstPointer input = this->GetInput(); mitk::Image::Pointer output = this->GetOutput(); mitk::ImageTimeSelector::Pointer timeSelector=mitk::ImageTimeSelector::New(); timeSelector->SetInput(input); mitkIpPicDescriptor* pic_transformed=NULL; pic_transformed = mitkIpPicNew(); pic_transformed->dim=3; pic_transformed->bpe = output->GetPixelType().GetBpe(); //pic_transformed->type = output->GetPixelType().GetType(); pic_transformed->n[0] = output->GetDimension(0); pic_transformed->n[1] = output->GetDimension(1); pic_transformed->n[2] = output->GetDimension(2); pic_transformed->data=malloc(_mitkIpPicSize(pic_transformed)); int nstart, nmax; int tstart, tmax; tstart=output->GetRequestedRegion().GetIndex(3); nstart=output->GetRequestedRegion().GetIndex(4); tmax=tstart+output->GetRequestedRegion().GetSize(3); nmax=nstart+output->GetRequestedRegion().GetSize(4); if(zt==NULL) { timeSelector->SetChannelNr(nstart); timeSelector->SetTimeNr(tstart); buildTransformShortCuts(input->GetDimension(0),input->GetDimension(1), input->GetDimension(2), output->GetDimension(0), rt_pic, phit_pic, fr_pic, fphi_pic, zt, fz); // query the line limiting the sector a=b=0; mitk::FloatProperty::Pointer prop; prop = dynamic_cast(input->GetProperty("SECTOR LIMITING LINE SLOPE").GetPointer()); if (prop.IsNotNull() ) a = prop->GetValue(); prop = dynamic_cast(input->GetProperty("SECTOR LIMITING LINE OFFSET").GetPointer()); if (prop.IsNotNull() ) b = prop->GetValue(); buildConeCutOffShortCut(input->GetDimension(0),input->GetDimension(1), rt_pic, fr_pic, a, b, coneCutOff_pic); // mitkIpPicPut("C:\\temp\\rt_90.pic",rt_pic); //mitkIpPicPut("C:\\temp\\coneCutOff.pic", coneCutOff_pic); } int n,t; for(n=nstart;nGetNumberOfChannels();++n) { timeSelector->SetChannelNr(n); for(t=tstart;tSetTimeNr(t); timeSelector->Update(); // Cast to pic descriptor for the timeSelector image mitkIpPicDescriptor* timeSelectorPic = mitkIpPicNew(); - CastToIpPicDescriptor( timeSelector->GetOutput(), timeSelectorPic ); + mitk::ImageWriteAccessor imageAccess(timeSelector->GetOutput()); + CastToIpPicDescriptor( timeSelector->GetOutput(), &imageAccess, timeSelectorPic ); _mitkIpPicFreeTags(pic_transformed->info->tags_head); pic_transformed->info->tags_head = _mitkIpPicCloneTags(timeSelectorPic->info->tags_head); if(input->GetDimension(2)>1) { mitkIpPicTypeMultiplex9(_transform, timeSelectorPic , pic_transformed, m_OutsideValue, (float*)fr_pic->data, (float*)fphi_pic->data, fz, (short *)rt_pic->data, (unsigned int *)phit_pic->data, zt, coneCutOff_pic); // mitkIpPicPut("1trf.pic",pic_transformed); } else { mitkIpPicDescriptor *doubleSlice = mitkIpPicCopyHeader( timeSelectorPic , NULL); doubleSlice->dim=3; doubleSlice->n[2]=2; doubleSlice->data=malloc(_mitkIpPicSize(doubleSlice)); memcpy(doubleSlice->data, timeSelectorPic->data, _mitkIpPicSize(doubleSlice)/2); mitkIpPicTypeMultiplex9(_transform, doubleSlice, pic_transformed, m_OutsideValue, (float*)fr_pic->data, (float*)fphi_pic->data, fz, (short *)rt_pic->data, (unsigned int *)phit_pic->data, zt, coneCutOff_pic); mitkIpPicFree(doubleSlice); } output->SetVolume(pic_transformed->data, t, n); } } //mitkIpPicPut("outzzzzzzzz.pic",pic_transformed); mitkIpPicFree(pic_transformed); m_TimeOfHeaderInitialization.Modified(); } mitk::CylindricToCartesianFilter::CylindricToCartesianFilter() : m_OutsideValue(0.0), m_TargetXSize(0) { rt_pic = NULL; phit_pic = NULL; fr_pic = NULL; fphi_pic = NULL; coneCutOff_pic = NULL; zt = NULL; fz = NULL; a=b=0.0; } mitk::CylindricToCartesianFilter::~CylindricToCartesianFilter() { if(rt_pic!=NULL) mitkIpPicFree(rt_pic); if(phit_pic!=NULL) mitkIpPicFree(phit_pic); if(fr_pic!=NULL) mitkIpPicFree(fr_pic); if(fphi_pic!=NULL) mitkIpPicFree(fphi_pic); if(coneCutOff_pic!=NULL) mitkIpPicFree(coneCutOff_pic); if(zt != NULL) free(zt); if(fz != NULL) free (fz); } void mitk::CylindricToCartesianFilter::GenerateInputRequestedRegion() { Superclass::GenerateInputRequestedRegion(); mitk::ImageToImageFilter::InputImagePointer input = const_cast< mitk::ImageToImageFilter::InputImageType * > ( this->GetInput() ); mitk::Image::Pointer output = this->GetOutput(); Image::RegionType requestedRegion; requestedRegion = output->GetRequestedRegion(); requestedRegion.SetIndex(0, 0); requestedRegion.SetIndex(1, 0); requestedRegion.SetIndex(2, 0); requestedRegion.SetSize(0, input->GetDimension(0)); requestedRegion.SetSize(1, input->GetDimension(1)); requestedRegion.SetSize(2, input->GetDimension(2)); input->SetRequestedRegion( & requestedRegion ); } diff --git a/Modules/MitkExt/Algorithms/mitkDopplerToStrainRateFilter.cpp b/Modules/MitkExt/Algorithms/mitkDopplerToStrainRateFilter.cpp index 361995c59d..05be5f7a25 100644 --- a/Modules/MitkExt/Algorithms/mitkDopplerToStrainRateFilter.cpp +++ b/Modules/MitkExt/Algorithms/mitkDopplerToStrainRateFilter.cpp @@ -1,384 +1,386 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkDopplerToStrainRateFilter.h" #include "mitkImageTimeSelector.h" #include "mitkProperties.h" #include "mitkPlaneGeometry.h" #include "mitkLegacyAdaptors.h" #include #include #include void mitk::DopplerToStrainRateFilter::GenerateOutputInformation() { mitk::Image::ConstPointer input = this->GetInput(); mitk::Image::Pointer output = this->GetOutput(); if ((output->IsInitialized()) && (this->GetMTime() <= m_TimeOfHeaderInitialization.GetMTime())) return; itkDebugMacro(<<"GenerateOutputInformation()"); unsigned int i; unsigned int *tmpDimensions = new unsigned int[input->GetDimension()]; for(i=0;iGetDimension();++i) tmpDimensions[i]=input->GetDimension(i); //@todo maybe we should shift the following somehow in ImageToImageFilter output->Initialize(input->GetPixelType(), input->GetDimension(), tmpDimensions, input->GetNumberOfChannels()); output->GetSlicedGeometry()->SetSpacing(input->GetSlicedGeometry()->GetSpacing()); //output->GetSlicedGeometry()->SetGeometry2D(mitk::Image::BuildStandardPlaneGeometry2D(output->GetSlicedGeometry(), tmpDimensions).GetPointer(), 0); //output->GetSlicedGeometry()->SetEvenlySpaced(); //set the timebounds - after SetGeometry2D, so that the already created PlaneGeometry will also receive this timebounds. output->GetSlicedGeometry()->SetTimeBounds(input->GetSlicedGeometry()->GetTimeBounds()); output->SetPropertyList(input->GetPropertyList()->Clone()); delete [] tmpDimensions; m_TimeOfHeaderInitialization.Modified(); } void mitk::DopplerToStrainRateFilter::GenerateData() { mitk::Image::ConstPointer input = this->GetInput(); mitk::Image::Pointer output = this->GetOutput(); mitk::Point3iProperty::Pointer pointProp; pointProp = dynamic_cast(input->GetProperty("ORIGIN").GetPointer()); if (pointProp.IsNotNull() ) { m_Origin = pointProp->GetValue(); } MITK_INFO << "compute Strain Rate Image .... " << std::endl << " origin[0]=" << m_Origin[0] << " origin[1]=" << m_Origin[1] << " origin[2]=" << m_Origin[2] << std::endl << " distance=" << m_Distance << std::endl << " NoStrainIntervall=" << m_NoStrainInterval << std::endl; const Vector3D & spacing = input->GetSlicedGeometry()->GetSpacing(); // MITK_INFO << " in: xres=" << spacing[0] << " yres=" << spacing[1] << " zres=" << spacing[2] << std::endl; mitk::ImageTimeSelector::Pointer timeSelector=mitk::ImageTimeSelector::New(); timeSelector->SetInput(input); mitkIpPicDescriptor* picStrainRate; picStrainRate = mitkIpPicNew(); picStrainRate->dim=3; picStrainRate->bpe = output->GetPixelType().GetBpe(); //picStrainRate->type = output->GetPixelType().GetType(); picStrainRate->n[0] = output->GetDimension(0); picStrainRate->n[1] = output->GetDimension(1); picStrainRate->n[2] = output->GetDimension(2); picStrainRate->data=malloc(_mitkIpPicSize(picStrainRate)); int xDim = picStrainRate->n[0]; int yDim = picStrainRate->n[1]; int zDim = picStrainRate->n[2]; long slice_size = xDim*yDim; long vol_size = slice_size*zDim; mitkIpPicDescriptor *picDoppler; int x,y,z;//,time; // loop-counter int strainRate; // the computed Strain Rate int v1,v2; // velocity and Point p1 and p2 float alpha; // the beam-angle, angle betwen current point and beam-point float dx=0, dy=0; // projection of this->distance to x- and y-axis int x1; // a square, where the velocity v1 lies in int y1; // the points are used for interpolation int minStrainRate=128, maxStrainRate=128; int n, nmax; int t, tmax; t = output->GetRequestedRegion().GetIndex(3); n = output->GetRequestedRegion().GetIndex(4); MITK_INFO << "t = " <data)[y*xDim + x] = (int) ( (alpha/1.6)*128); if (!isAnglePicWritten) ((mitkIpInt1_t *)anglePic->data)[y*xDim + x] = (int) ( dx); #endif } // x } // y //isAnglePicWritten = mitkIpTrue; } // z //isStrainComputed = mitkIpTrue; std::string filenameD; filenameD ="doppler.pic"; mitkIpPicPut(const_cast(filenameD.c_str()),picDoppler); #define WRITE_STRAIN_PIC #ifdef WRITE_STRAIN_PIC char tmpfilename[100]; sprintf(tmpfilename,"strain%d.pic",t);; mitkIpPicPut(tmpfilename,picStrainRate); #endif #ifdef WRITE_ANGLE_PIC std::string filename2; filename2="angle.pic"; mitkIpPicPut(const_cast(filename2.c_str()),anglePic); #endif ((mitkIpUInt1_t *)picStrainRate->data)[0]=0; ((mitkIpUInt1_t *)picStrainRate->data)[1]=255; output->SetVolume(picStrainRate->data, t, n); } } mitkIpPicFree(picStrainRate); #define WRITE_STRAIN_PIC #ifdef WRITE_STRAIN_PIC // Get the StrainRate ipPic descriptor picStrainRate = mitkIpPicNew(); - CastToIpPicDescriptor( output, picStrainRate ); + ImageWriteAccessor imageAccess(output); + CastToIpPicDescriptor( output, &imageAccess, picStrainRate ); std::string filename; filename ="strain.pic"; mitkIpPicPut(const_cast(filename.c_str()),picStrainRate); #endif MITK_INFO << "Strain Rate Image computed.... " << std::endl << " minStrainRate: " << minStrainRate << std::endl << " maxStrainRate: " << maxStrainRate << std::endl; } mitk::DopplerToStrainRateFilter::DopplerToStrainRateFilter() : m_Distance(10), m_NoStrainInterval(2) { m_Origin[0] = 0; m_Origin[1] = 0; m_Origin[2] = 0; } float mitk::DopplerToStrainRateFilter::GetLimit() { return (128/m_Distance); } mitk::DopplerToStrainRateFilter::~DopplerToStrainRateFilter() { } void mitk::DopplerToStrainRateFilter::GenerateInputRequestedRegion() { Superclass::GenerateInputRequestedRegion(); mitk::ImageToImageFilter::InputImagePointer input = const_cast< mitk::ImageToImageFilter::InputImageType * > ( this->GetInput() ); mitk::Image::Pointer output = this->GetOutput(); Image::RegionType requestedRegion; requestedRegion = output->GetRequestedRegion(); requestedRegion.SetIndex(0, 0); requestedRegion.SetIndex(1, 0); requestedRegion.SetIndex(2, 0); //requestedRegion.SetIndex(3, 0); //requestedRegion.SetIndex(4, 0); requestedRegion.SetSize(0, input->GetDimension(0)); requestedRegion.SetSize(1, input->GetDimension(1)); requestedRegion.SetSize(2, input->GetDimension(2)); //requestedRegion.SetSize(3, output->GetDimension(3)); //requestedRegion.SetSize(4, output->GetNumberOfChannels()); input->SetRequestedRegion( & requestedRegion ); } diff --git a/Modules/MitkExt/Algorithms/mitkMeshUtil.h b/Modules/MitkExt/Algorithms/mitkMeshUtil.h index 36745b223c..784dad5881 100644 --- a/Modules/MitkExt/Algorithms/mitkMeshUtil.h +++ b/Modules/MitkExt/Algorithms/mitkMeshUtil.h @@ -1,1570 +1,1604 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef MITKMESHUTIL_H_INCLUDED #define MITKMESHUTIL_H_INCLUDED #if(_MSC_VER==1200) #error MeshUtils currently not supported for MS Visual C++ 6.0. Sorry. #endif //#include #include #include #include #include #include //#include #include //#include //#include //#include #include #include #include #include #include #include #include #include #include #include #include #include #include #include template class NullScalarAccessor { public: static inline vtkFloatingPointType GetPointScalar(typename MeshType::PointDataContainer* /*pointData*/, typename MeshType::PointIdentifier /*idx*/, MeshType* /*mesh*/ = NULL, unsigned int /*type*/ = 0) { return (vtkFloatingPointType) 0.0; }; static inline vtkFloatingPointType GetCellScalar(typename MeshType::CellDataContainer* /*cellData*/, typename MeshType::CellIdentifier /*idx*/, MeshType* /*mesh*/ = NULL, unsigned int /*type*/ = 0) { return (vtkFloatingPointType) 0.0; }; }; template class MeshScalarAccessor { public: static inline vtkFloatingPointType GetPointScalar(typename MeshType::PointDataContainer* pointData, typename MeshType::PointIdentifier idx, MeshType* /*mesh*/ = NULL, unsigned int /*type*/ = 0) { return (vtkFloatingPointType)pointData->GetElement(idx); }; static inline vtkFloatingPointType GetCellScalar(typename MeshType::CellDataContainer* cellData, typename MeshType::CellIdentifier idx, MeshType* /*mesh*/ = NULL, unsigned int /*type*/ = 0) { return (vtkFloatingPointType)cellData->GetElement(idx); }; }; template class MeanCurvatureAccessor : public NullScalarAccessor { public: static inline vtkFloatingPointType GetPointScalar(typename MeshType::PointDataContainer* /*point*/, typename MeshType::PointIdentifier idx, MeshType* mesh, unsigned int /*type*/ = 0) { typename MeshType::PixelType dis = 0; mesh->GetPointData(idx, &dis); return (vtkFloatingPointType) dis; }; }; template class SimplexMeshAccessor : public NullScalarAccessor { public: static inline vtkFloatingPointType GetPointScalar(typename MeshType::PointDataContainer* point, typename MeshType::PointIdentifier idx, MeshType* mesh, unsigned int type = 0 ) { typename MeshType::GeometryMapPointer geometryData = mesh->GetGeometryData(); if (type == 0) { double val = mesh->GetMeanCurvature( idx ); mesh->SetPointData(idx, val); return val; } else if (type == 1) { double val = geometryData->GetElement(idx)->meanTension; mesh->SetPointData(idx, val); return val; } else if (type == 2) { double val = geometryData->GetElement(idx)->externalForce.GetNorm(); mesh->SetPointData(idx, val); return val; } else if (type == 3) return geometryData->GetElement(idx)->internalForce.GetNorm(); else if (type == 4) return geometryData->GetElement(idx)->externalForce.GetNorm() * mesh->GetDistance(idx); else if (type == 5) { typename MeshType::PixelType dis = 0; mesh->GetPointData(idx, &dis); return (vtkFloatingPointType) dis; } else if (type == 6) { return (vtkFloatingPointType) ((geometryData->GetElement(idx))->allowSplitting); } else return (vtkFloatingPointType) 0; }; }; /*! \brief The class provides mehtods for ITK - VTK mesh conversion * * \todo document the inner class * \todo maybe inner class should be moved out */ template > class MeshUtil { /*! \brief A visitor to create VTK cells by means of a class defining the InsertImplementation interface The InsertImplementation interface defines the methods \code void InsertLine(vtkIdType *pts); void InsertTriangle(vtkIdType *pts); void InsertPolygon(vtkIdType npts, vtkIdType *pts); void InsertQuad(vtkIdType *pts); void InsertTetra(vtkIdType *pts); void InsertHexahedron(vtkIdType *pts); \endcode This class calls the appropriate insert-method of the InsertImplementation according to the cell type of the visited cell \em and its actual contents: e.g., for a polygon cell with just two points, a line will be created by calling InsertLine. \sa ExactSwitchByCellType \sa SingleCellArrayInsertImplementation \sa DistributeInsertImplementation */ template class SwitchByCellType : public InsertImplementation { // typedef the itk cells we are interested in typedef typename itk::CellInterface< typename MeshType::CellPixelType, typename MeshType::CellTraits > CellInterfaceType; typedef itk::LineCell floatLineCell; typedef itk::TriangleCell floatTriangleCell; typedef itk::PolygonCell floatPolygonCell; typedef itk::QuadrilateralCell floatQuadrilateralCell; typedef itk::TetrahedronCell floatTetrahedronCell; typedef itk::HexahedronCell floatHexahedronCell; typedef typename CellInterfaceType::PointIdConstIterator PointIdIterator; public: /*! Visit a line and create the VTK_LINE cell */ void Visit(unsigned long cellId, floatLineCell* t) { vtkIdType pts[2]; int i=0; unsigned long num = t->GetNumberOfVertices(); vtkIdType vtkCellId = -1; if (num==2) { // useless because itk::LineCell always returns 2 for (PointIdIterator it=t->PointIdsBegin(); it!=t->PointIdsEnd(); it++) pts[i++] = *it; vtkCellId = this->InsertLine( (vtkIdType*)pts ); } if (this->m_UseCellScalarAccessor && vtkCellId >= 0) { this->m_CellScalars->InsertTuple1(vtkCellId, ScalarAccessor::GetCellScalar(this->m_CellData, cellId)); } } /*! Visit a polygon and create the VTK_POLYGON cell */ void Visit(unsigned long cellId, floatPolygonCell* t) { vtkIdType pts[4096]; int i=0; unsigned long num = t->GetNumberOfVertices(); vtkIdType vtkCellId = -1; if (num > 4096) { MITK_ERROR << "Problem in mitkMeshUtil: Polygon with more than maximum number of vertices encountered." << std::endl; } else if (num > 3) { for (PointIdIterator it=t->PointIdsBegin(); it!=t->PointIdsEnd(); it++) pts[i++] = *it; vtkCellId = this->InsertPolygon( num, (vtkIdType*)pts ); } else if (num == 3) { for (PointIdIterator it=t->PointIdsBegin(); it!=t->PointIdsEnd(); it++) pts[i++] = *it; vtkCellId = this->InsertTriangle( (vtkIdType*)pts ); } else if (num==2) { for (PointIdIterator it=t->PointIdsBegin(); it!=t->PointIdsEnd(); it++) pts[i++] = *it; vtkCellId = this->InsertLine( (vtkIdType*)pts ); } if (this->m_UseCellScalarAccessor && vtkCellId >= 0) { this->m_CellScalars->InsertTuple1(vtkCellId, ScalarAccessor::GetCellScalar(this->m_CellData, cellId)); } } /*! Visit a triangle and create the VTK_TRIANGLE cell */ void Visit(unsigned long cellId, floatTriangleCell* t) { vtkIdType pts[3]; int i=0; unsigned long num = t->GetNumberOfVertices(); vtkIdType vtkCellId = -1; if (num == 3) { for (PointIdIterator it=t->PointIdsBegin(); it!=t->PointIdsEnd(); it++) pts[i++] = *it; vtkCellId = this->InsertTriangle( (vtkIdType*)pts ); } else if (num==2) { for (PointIdIterator it=t->PointIdsBegin(); it!=t->PointIdsEnd(); it++) pts[i++] = *it; vtkCellId = this->InsertLine( (vtkIdType*)pts ); } if (this->m_UseCellScalarAccessor && vtkCellId >= 0) { this->m_CellScalars->InsertTuple1(vtkCellId, ScalarAccessor::GetCellScalar(this->m_CellData, cellId)); } } /*! Visit a quad and create the VTK_QUAD cell */ void Visit(unsigned long cellId, floatQuadrilateralCell* t) { vtkIdType pts[4]; int i=0; unsigned long num = t->GetNumberOfVertices(); vtkIdType vtkCellId = -1; if (num == 4) { for (PointIdIterator it=t->PointIdsBegin(); it!=t->PointIdsEnd(); it++) { if (i == 2) pts[3] = *it; else if (i == 3) pts[2] = *it; else pts[i] = *it; i++; //pts[i++] = *it; } vtkCellId = this->InsertQuad( (vtkIdType*)pts ); } else if (num == 3) { for (PointIdIterator it=t->PointIdsBegin(); it!=t->PointIdsEnd(); it++) pts[i++] = *it; vtkCellId = this->InsertTriangle( (vtkIdType*)pts ); } else if (num==2) { for (PointIdIterator it=t->PointIdsBegin(); it!=t->PointIdsEnd(); it++) pts[i++] = *it; vtkCellId = this->InsertLine( (vtkIdType*)pts ); } if (this->m_UseCellScalarAccessor && vtkCellId >= 0) { this->m_CellScalars->InsertTuple1(vtkCellId, ScalarAccessor::GetCellScalar(this->m_CellData, cellId)); } } /*! Visit a tetrahedra and create the VTK_TETRA cell */ void Visit(unsigned long cellId, floatTetrahedronCell* t) { vtkIdType pts[4]; int i=0; unsigned long num = t->GetNumberOfVertices(); vtkIdType vtkCellId = -1; if (num == 4) { for (PointIdIterator it=t->PointIdsBegin(); it!=t->PointIdsEnd(); it++) pts[i++] = *it; vtkCellId = this->InsertTetra( (vtkIdType*)pts ); } else if (num == 3) { for (PointIdIterator it=t->PointIdsBegin(); it!=t->PointIdsEnd(); it++) pts[i++] = *it; vtkCellId = this->InsertTriangle( (vtkIdType*)pts ); } else if (num==2) { for (PointIdIterator it=t->PointIdsBegin(); it!=t->PointIdsEnd(); it++) pts[i++] = *it; vtkCellId = this->InsertLine( (vtkIdType*)pts ); } if (this->m_UseCellScalarAccessor && vtkCellId >= 0) { this->m_CellScalars->InsertTuple1(vtkCellId, ScalarAccessor::GetCellScalar(this->m_CellData, cellId)); } } /*! Visit a hexahedron and create the VTK_HEXAHEDRON cell */ void Visit(unsigned long cellId, floatHexahedronCell* t) { vtkIdType pts[8]; int i=0; unsigned long num = t->GetNumberOfVertices(); vtkIdType vtkCellId = -1; if (num == 8) { for (PointIdIterator it=t->PointIdsBegin(); it!=t->PointIdsEnd(); it++) { if (i == 2) pts[i++] = *(it+1); else if (i == 3) pts[i++] = *(it-1); else if (i == 6) pts[i++] = *(it+1); else if (i == 7) pts[i++] = *(it-1); else pts[i++] = *it; } vtkCellId = this->InsertHexahedron( (vtkIdType*)pts ); } else if (num == 4) { for (PointIdIterator it=t->PointIdsBegin(); it!=t->PointIdsEnd(); it++) pts[i++] = *it; vtkCellId = this->InsertQuad( (vtkIdType*)pts ); } else if (num == 3) { for (PointIdIterator it=t->PointIdsBegin(); it!=t->PointIdsEnd(); it++) pts[i++] = *it; vtkCellId = this->InsertTriangle( (vtkIdType*)pts ); } else if (num==2) { for (PointIdIterator it=t->PointIdsBegin(); it!=t->PointIdsEnd(); it++) pts[i++] = *it; vtkCellId = this->InsertLine( (vtkIdType*)pts ); } if (this->m_UseCellScalarAccessor && vtkCellId >= 0) { this->m_CellScalars->InsertTuple1(vtkCellId, ScalarAccessor::GetCellScalar(this->m_CellData, cellId)); } } }; /*! \brief A visitor similar to SwitchByCellType, but with exact matching of cell types Works as described in SwitchByCellType, but does exact matching of cell types, e.g., for a polygon cell with just two points, \em no insert-method will be called, because a polygon must have at least three points. \sa SwitchByCellType \sa SingleCellArrayInsertImplementation \sa DistributeInsertImplementation */ template class ExactSwitchByCellType : public InsertImplementation { // typedef the itk cells we are interested in typedef typename itk::CellInterface< typename MeshType::CellPixelType, typename MeshType::CellTraits > CellInterfaceType; typedef itk::LineCell floatLineCell; typedef itk::TriangleCell floatTriangleCell; typedef itk::PolygonCell floatPolygonCell; typedef itk::QuadrilateralCell floatQuadrilateralCell; typedef itk::TetrahedronCell floatTetrahedronCell; typedef itk::HexahedronCell floatHexahedronCell; + typedef typename CellInterfaceType::PointIdConstIterator PointIdIterator; public: /*! Visit a line and create the VTK_LINE cell */ void Visit(unsigned long , floatLineCell* t) { unsigned long num = t->GetNumberOfVertices(); - vtkIdType *pts = (vtkIdType*)t->PointIdsBegin(); + vtkIdType pts[2]; + int i = 0; + if (num==2) + { + for (PointIdIterator it=t->PointIdsBegin(); it!=t->PointIdsEnd(); it++) pts[i++] = *it; this->InsertLine(pts); + } } /*! Visit a polygon and create the VTK_POLYGON cell */ void Visit(unsigned long , floatPolygonCell* t) { + vtkIdType pts[4096]; unsigned long num = t->GetNumberOfVertices(); - vtkIdType *pts = (vtkIdType*)t->PointIdsBegin(); + if (num > 4096) { + MITK_ERROR << "Problem in mitkMeshUtil: Polygon with more than maximum number of vertices encountered." << std::endl; + } + int i = 0; + if (num > 3) + { + for (PointIdIterator it=t->PointIdsBegin(); it!=t->PointIdsEnd(); it++) pts[i++] = *it; this->InsertPolygon(num, pts); + } } /*! Visit a triangle and create the VTK_TRIANGLE cell */ void Visit(unsigned long , floatTriangleCell* t) { unsigned long num = t->GetNumberOfVertices(); - vtkIdType *pts = (vtkIdType*)t->PointIdsBegin(); + vtkIdType pts[3]; + int i = 0; + if (num == 3) + { + for (PointIdIterator it=t->PointIdsBegin(); it!=t->PointIdsEnd(); it++) pts[i++] = *it; this->InsertTriangle(pts); + } } /*! Visit a quadrilateral and create the VTK_QUAD cell */ void Visit(unsigned long , floatQuadrilateralCell* t) { unsigned long num = t->GetNumberOfVertices(); - vtkIdType *pts = (vtkIdType*)t->PointIdsBegin(); - if (num == 4) { - vtkIdType tmpId = pts[3]; - pts[3] = pts[4]; - pts[4] = tmpId; + vtkIdType pts[4]; + int i = 0; + + if (num == 4) + { + for (PointIdIterator it=t->PointIdsBegin(); it!=t->PointIdsEnd(); it++) pts[i++] = *it; + + vtkIdType tmpId = pts[2]; + pts[2] = pts[3]; + pts[3] = tmpId; this->InsertQuad(pts); } } /*! Visit a tetrahedron and create the VTK_TETRA cell */ void Visit(unsigned long , floatTetrahedronCell* t) { unsigned long num = t->GetNumberOfVertices(); - vtkIdType *pts = (vtkIdType*)t->PointIdsBegin(); + + vtkIdType pts[4]; + int i = 0; + if (num == 4) + { + for (PointIdIterator it=t->PointIdsBegin(); it!=t->PointIdsEnd(); it++) pts[i++] = *it; this->InsertTetra(pts); + } } /*! Visit a hexahedron and create the VTK_HEXAHEDRON cell */ void Visit(unsigned long , floatHexahedronCell* t) { unsigned long num = t->GetNumberOfVertices(); - vtkIdType *pts = (vtkIdType*)t->PointIdsBegin(); + vtkIdType pts[8]; + int i = 0; + if (num == 8) { + for (PointIdIterator it=t->PointIdsBegin(); it!=t->PointIdsEnd(); it++) pts[i++] = *it; + vtkIdType tmp[8]; for (unsigned int i = 0; i < 8; i++) tmp[i] = pts[i]; pts[2] = tmp[3]; pts[3] = tmp[2]; pts[6] = tmp[7]; pts[7] = tmp[6]; this->InsertHexahedron(pts); } } }; /*! \brief Implementation of the InsertImplementation interface of SwitchByCellType to define a visitor that create cells according to their types and put them in a single vtkCellArray (for vtkUnstructuredGrid construction) */ class SingleCellArrayInsertImplementation { vtkCellArray* m_Cells; int* m_TypeArray; //vtkIdType cellId; protected: bool m_UseCellScalarAccessor; vtkFloatArray* m_CellScalars; typename MeshType::CellDataContainer::Pointer m_CellData; public: SingleCellArrayInsertImplementation() : m_UseCellScalarAccessor(false) {} /*! Set the vtkCellArray that will be constructed */ void SetCellArray(vtkCellArray* cells) { m_Cells = cells; } /*! Set the type array for storing the vtk cell types */ void SetTypeArray(int* i) { m_TypeArray = i; } void SetUseCellScalarAccessor(bool flag) { m_UseCellScalarAccessor = flag; } void SetCellScalars(vtkFloatArray* scalars) { m_CellScalars = scalars; } vtkFloatArray* GetCellScalars() { return m_CellScalars; } void SetMeshCellData(typename MeshType::CellDataContainer* data) { m_CellData = data; } vtkIdType InsertLine(vtkIdType *pts) { vtkIdType cellId = m_Cells->InsertNextCell(2, pts); m_TypeArray[cellId] = VTK_LINE; return cellId; } vtkIdType InsertTriangle(vtkIdType *pts) { vtkIdType cellId = m_Cells->InsertNextCell(3, pts); m_TypeArray[cellId] = VTK_TRIANGLE; return cellId; } vtkIdType InsertPolygon(vtkIdType npts, vtkIdType *pts) { vtkIdType cellId = m_Cells->InsertNextCell(npts, pts); m_TypeArray[cellId] = VTK_POLYGON; return cellId; } vtkIdType InsertQuad(vtkIdType *pts) { vtkIdType cellId = m_Cells->InsertNextCell(4, pts); m_TypeArray[cellId] = VTK_QUAD; return cellId; } vtkIdType InsertTetra(vtkIdType *pts) { vtkIdType cellId = m_Cells->InsertNextCell(4, pts); m_TypeArray[cellId] = VTK_TETRA; return cellId; } vtkIdType InsertHexahedron(vtkIdType *pts) { vtkIdType cellId = m_Cells->InsertNextCell(8, pts); m_TypeArray[cellId] = VTK_HEXAHEDRON; return cellId; } }; /*! \brief Implementation of the InsertImplementation interface of SwitchByCellType to define a visitor that distributes cells according to their types (for vtkPolyData construction) */ class DistributeInsertImplementation { vtkCellArray* m_LineCells; vtkCellArray* m_TriangleCells; vtkCellArray* m_PolygonCells; vtkCellArray* m_QuadCells; protected: bool m_UseCellScalarAccessor; vtkFloatArray* m_CellScalars; typename MeshType::CellDataContainer::Pointer m_CellData; public: DistributeInsertImplementation() : m_UseCellScalarAccessor(false) {} /*! Set the vtkCellArray that will be constructed */ void SetCellArrays(vtkCellArray* lines, vtkCellArray* triangles, vtkCellArray* polygons, vtkCellArray* quads) { m_LineCells = lines; m_TriangleCells = triangles; m_PolygonCells = polygons; m_QuadCells = quads; } vtkIdType InsertLine(vtkIdType *pts) { return m_LineCells->InsertNextCell(2, pts); } vtkIdType InsertTriangle(vtkIdType *pts) { return m_TriangleCells->InsertNextCell(3, pts); } vtkIdType InsertPolygon(vtkIdType npts, vtkIdType *pts) { return m_PolygonCells->InsertNextCell(npts, pts); } vtkIdType InsertQuad(vtkIdType *pts) { return m_QuadCells->InsertNextCell(4, pts); } vtkIdType InsertTetra(vtkIdType *pts) { return -1; } // ignored vtkIdType InsertHexahedron(vtkIdType *pts) { return -1; } // ignored }; //typedef typename MeshType::CellType CellType; //typedef typename itk::LineCell< CellType > LineType; //typedef typename itk::PolygonCell< CellType > PolygonType; //typedef typename itk::TriangleCell< CellType > TriangleType; typedef SwitchByCellType SingleCellArrayUserVisitorType; typedef SwitchByCellType DistributeUserVisitorType; typedef ExactSwitchByCellType ExactUserVisitorType; public: typedef itk::MatrixOffsetTransformBase ITKTransformType; typedef itk::MatrixOffsetTransformBase MITKTransformType; /*! Convert a MITK transformation to an ITK transformation Necessary because ITK uses double and MITK uses float values */ static void ConvertTransformToItk(const MITKTransformType* mitkTransform, ITKTransformType* itkTransform) { typename MITKTransformType::MatrixType mitkM = mitkTransform->GetMatrix(); typename ITKTransformType::MatrixType itkM; typename MITKTransformType::OffsetType mitkO = mitkTransform->GetOffset(); typename ITKTransformType::OffsetType itkO; for(short i = 0; i < 3; ++i) { for(short j = 0; j<3; ++j) { itkM[i][j] = (double)mitkM[i][j]; } itkO[i] = (double)mitkO[i]; } itkTransform->SetMatrix(itkM); itkTransform->SetOffset(itkO); } /*! create an itkMesh object from a vtkPolyData */ static typename MeshType::Pointer MeshFromPolyData(vtkPolyData* poly, mitk::Geometry3D* geometryFrame=NULL, mitk::Geometry3D* polyDataGeometryFrame=NULL) { // Create a new mesh typename MeshType::Pointer output = MeshType::New(); output->SetCellsAllocationMethod( MeshType::CellsAllocatedDynamicallyCellByCell ); typedef typename MeshType::CellDataContainer MeshCellDataContainerType; output->SetCellData(MeshCellDataContainerType::New()); // Get the points from vtk vtkPoints* vtkpoints = poly->GetPoints(); const unsigned int numPoints = poly->GetNumberOfPoints(); // Create a compatible point container for the mesh // the mesh is created with a null points container // MeshType::PointsContainer::Pointer points = // MeshType::PointsContainer::New(); // // Resize the point container to be able to fit the vtk points // points->Reserve(numPoints); // // Set the point container on the mesh //output->SetPoints(points); vtkFloatingPointType vtkpoint[3]; typename MeshType::PointType itkPhysicalPoint; if(geometryFrame==NULL) { if(polyDataGeometryFrame==NULL) { for(unsigned int i=0; i < numPoints; ++i) { vtkpoints->GetPoint(i, vtkpoint); //MITK_INFO << "next point: " << test[0]<< "," << test[1] << "," << test[2] << std::endl; //typename MeshType::PixelType* apoint = (typename MeshType::PixelType*) vtkpoints->GetPoint(i); mitk::vtk2itk(vtkpoint, itkPhysicalPoint); output->SetPoint( i, itkPhysicalPoint ); } } else { for(unsigned int i=0; i < numPoints; ++i) { vtkpoints->GetPoint(i, vtkpoint); //MITK_INFO << "next point: " << test[0]<< "," << test[1] << "," << test[2] << std::endl; //typename MeshType::PixelType* apoint = (typename MeshType::PixelType*) vtkpoints->GetPoint(i); mitk::Point3D mitkWorldPoint; mitk::vtk2itk(vtkpoint, mitkWorldPoint); polyDataGeometryFrame->IndexToWorld(mitkWorldPoint, mitkWorldPoint); mitk::vtk2itk(mitkWorldPoint, itkPhysicalPoint); output->SetPoint( i, itkPhysicalPoint ); } } } else { mitk::Point3D mitkWorldPoint; if(polyDataGeometryFrame==NULL) { for(unsigned int i=0; i < numPoints; ++i) { vtkpoints->GetPoint(i, vtkpoint); //MITK_INFO << "next point: " << test[0]<< "," << test[1] << "," << test[2] << std::endl; //typename MeshType::PixelType* apoint = (typename MeshType::PixelType*) vtkpoints->GetPoint(i); mitk::vtk2itk(vtkpoint, mitkWorldPoint); geometryFrame->WorldToItkPhysicalPoint(mitkWorldPoint, itkPhysicalPoint); output->SetPoint( i, itkPhysicalPoint ); } } else { for(unsigned int i=0; i < numPoints; ++i) { vtkpoints->GetPoint(i, vtkpoint); //MITK_INFO << "next point: " << test[0]<< "," << test[1] << "," << test[2] << std::endl; //typename MeshType::PixelType* apoint = (typename MeshType::PixelType*) vtkpoints->GetPoint(i); mitk::vtk2itk(vtkpoint, mitkWorldPoint); polyDataGeometryFrame->IndexToWorld(mitkWorldPoint, mitkWorldPoint); geometryFrame->WorldToItkPhysicalPoint(mitkWorldPoint, itkPhysicalPoint); output->SetPoint( i, itkPhysicalPoint ); } } } vtkCellArray* vtkcells = poly->GetPolys(); // vtkCellArray* vtkcells = poly->GetStrips(); //MeshType::CellsContainerPointer cells = MeshType::CellsContainer::New(); //output->SetCells(cells); // extract the cell id's from the vtkUnstructuredGrid int numcells = vtkcells->GetNumberOfCells(); int* vtkCellTypes = new int[numcells]; int cellId = 0; // poly ids start after verts and lines! int cellIdOfs = poly->GetNumberOfVerts() + poly->GetNumberOfLines(); for(; cellId < numcells; cellId++) { vtkCellTypes[cellId] = poly->GetCellType( cellId+cellIdOfs ); } // cells->Reserve(numcells); vtkIdType npts; vtkIdType* pts; cellId = 0; typedef typename MeshType::MeshTraits OMeshTraits; typedef typename OMeshTraits::PixelType OPixelType; typedef typename MeshType::CellTraits CellTraits; typedef typename itk::CellInterface CellInterfaceType; typedef typename itk::TriangleCell TriCellType; typedef typename TriCellType::CellAutoPointer TriCellPointer; TriCellPointer newCell; output->GetCells()->Reserve( poly->GetNumberOfPolys() + poly->GetNumberOfStrips() ); output->GetCellData()->Reserve( poly->GetNumberOfPolys() + poly->GetNumberOfStrips() ); for(vtkcells->InitTraversal(); vtkcells->GetNextCell(npts, pts); cellId++) { switch(vtkCellTypes[cellId]) { case VTK_TRIANGLE: { if (npts != 3) continue; // skip non-triangles; unsigned long pointIds[3]; pointIds[0] = (unsigned long) pts[0]; pointIds[1] = (unsigned long) pts[1]; pointIds[2] = (unsigned long) pts[2]; newCell.TakeOwnership( new TriCellType ); newCell->SetPointIds(pointIds);//(unsigned long*)pts); output->SetCell(cellId, newCell ); output->SetCellData(cellId, (typename MeshType::PixelType)3); break; } case VTK_QUAD: { if (npts != 4 ) continue; // skip non-quadrilateral unsigned long pointIds[3]; pointIds[0] = (unsigned long) pts[0]; pointIds[1] = (unsigned long) pts[1]; pointIds[2] = (unsigned long) pts[2]; newCell.TakeOwnership( new TriCellType ); newCell->SetPointIds(pointIds); output->SetCell(cellId, newCell ); output->SetCellData(cellId, (typename MeshType::PixelType)3); cellId++; pointIds[0] = (unsigned long) pts[2]; pointIds[1] = (unsigned long) pts[3]; pointIds[2] = (unsigned long) pts[0]; newCell.TakeOwnership( new TriCellType ); newCell->SetPointIds(pointIds); output->SetCell(cellId, newCell ); output->SetCellData(cellId, (typename MeshType::PixelType)3); break; } case VTK_EMPTY_CELL: { if (npts != 3) { MITK_ERROR << "Only empty triangle cell supported by now..." << std::endl; // skip non-triangle empty cells; continue; } unsigned long pointIds[3]; pointIds[0] = (unsigned long) pts[0]; pointIds[1] = (unsigned long) pts[1]; pointIds[2] = (unsigned long) pts[2]; newCell.TakeOwnership( new TriCellType ); newCell->SetPointIds(pointIds); output->SetCell(cellId, newCell ); output->SetCellData(cellId, (typename MeshType::PixelType)3); break; } //case VTK_VERTEX: // If need to implement use //case VTK_POLY_VERTEX: // the poly->GetVerts() and //case VTK_LINE: // poly->GetLines() routines //case VTK_POLY_LINE: // outside of the switch..case. case VTK_POLYGON: case VTK_PIXEL: { if (npts != 4 ) continue;// skip non-quadrilateral unsigned long pointIds[3]; for ( unsigned int idx = 0; idx <= 1; idx++ ) { pointIds[0] = (unsigned long) pts[idx]; pointIds[1] = (unsigned long) pts[idx+1]; pointIds[2] = (unsigned long) pts[idx+2]; newCell.TakeOwnership( new TriCellType ); newCell->SetPointIds(pointIds); output->SetCell(cellId+idx, newCell ); output->SetCellData(cellId+idx, (typename MeshType::PixelType)3); } cellId++; break; } case VTK_TETRA: case VTK_VOXEL: case VTK_HEXAHEDRON: case VTK_WEDGE: case VTK_PYRAMID: case VTK_PARAMETRIC_CURVE: case VTK_PARAMETRIC_SURFACE: default: MITK_WARN << "Warning, unhandled cell type " << vtkCellTypes[cellId] << std::endl; } } if (poly->GetNumberOfStrips() != 0) { vtkcells = poly->GetStrips(); numcells = vtkcells->GetNumberOfCells(); vtkCellTypes = new int[numcells]; int stripId = 0; // strip ids start after verts, lines and polys! int stripIdOfs = poly->GetNumberOfVerts() + poly->GetNumberOfLines() + poly->GetNumberOfPolys(); for(; stripId < numcells; stripId++) { vtkCellTypes[stripId] = poly->GetCellType( stripId+stripIdOfs ); } stripId = 0; vtkcells->InitTraversal(); while( vtkcells->GetNextCell(npts, pts) ) { if (vtkCellTypes[stripId] != VTK_TRIANGLE_STRIP) { MITK_ERROR << "Only triangle strips supported!" << std::endl; continue; } stripId++; unsigned int numberOfTrianglesInStrip = npts - 2; unsigned long pointIds[3]; pointIds[0] = (unsigned long) pts[0]; pointIds[1] = (unsigned long) pts[1]; pointIds[2] = (unsigned long) pts[2]; for( unsigned int t=0; t < numberOfTrianglesInStrip; t++ ) { newCell.TakeOwnership( new TriCellType ); newCell->SetPointIds(pointIds); output->SetCell(cellId, newCell ); output->SetCellData(cellId, (typename MeshType::PixelType)3); cellId++; pointIds[0] = pointIds[1]; pointIds[1] = pointIds[2]; pointIds[2] = pts[t+3]; } } } //output->Print(std::cout); output->BuildCellLinks(); delete[] vtkCellTypes; return output; } /*! create an itkMesh object from an mitk::Surface */ static typename MeshType::Pointer MeshFromSurface(mitk::Surface* surface, mitk::Geometry3D* geometryFrame=NULL) { if(surface == NULL) return NULL; return MeshFromPolyData(surface->GetVtkPolyData(), geometryFrame, surface->GetGeometry()); } /*! create an vtkUnstructuredGrid object from an itkMesh */ static vtkUnstructuredGrid* MeshToUnstructuredGrid( MeshType* mesh, bool usePointScalarAccessor = false, bool useCellScalarAccessor = false, unsigned int pointDataType = 0, mitk::Geometry3D* geometryFrame=NULL) { /*! default SingleCellArray line cell visitior definition */ typedef typename itk::CellInterfaceVisitorImplementation, SingleCellArrayUserVisitorType> SingleCellArrayLineVisitor; /*! default SingleCellArray polygon cell visitior definition */ typedef typename itk::CellInterfaceVisitorImplementation, SingleCellArrayUserVisitorType> SingleCellArrayPolygonVisitor; /*! default SingleCellArray triangle cell visitior definition */ typedef typename itk::CellInterfaceVisitorImplementation >, SingleCellArrayUserVisitorType> SingleCellArrayTriangleVisitor; /*! default SingleCellArray quad cell visitior definition */ typedef typename itk::CellInterfaceVisitorImplementation >, SingleCellArrayUserVisitorType> SingleCellArrayQuadrilateralVisitor; /*! default SingleCellArray tetra cell visitior definition */ typedef typename itk::CellInterfaceVisitorImplementation >, SingleCellArrayUserVisitorType> SingleCellArrayTetrahedronVisitor; /*! default SingleCellArray hex cell visitior definition */ typedef typename itk::CellInterfaceVisitorImplementation >, SingleCellArrayUserVisitorType> SingleCellArrayHexahedronVisitor; // Get the number of points in the mesh int numPoints = mesh->GetNumberOfPoints(); if(numPoints == 0) { //mesh->Print(std::cerr); MITK_FATAL << "no points in Grid " << std::endl; exit(-1); } // Create a vtkUnstructuredGrid vtkUnstructuredGrid* vgrid = vtkUnstructuredGrid::New(); // Create the vtkPoints object and set the number of points vtkPoints* vpoints = vtkPoints::New( VTK_DOUBLE ); vtkFloatArray* pointScalars = vtkFloatArray::New(); vtkFloatArray* cellScalars = vtkFloatArray::New(); pointScalars->SetNumberOfComponents(1); cellScalars->SetNumberOfComponents(1); typename MeshType::PointsContainer::Pointer points = mesh->GetPoints(); typename MeshType::PointsContainer::Iterator i; // iterate over all the points in the itk mesh to find // the maximal index unsigned int maxIndex = 0; for(i = points->Begin(); i != points->End(); ++i) { if(maxIndex < i->Index()) maxIndex = i->Index(); } // initialize vtk-classes for points and scalars vpoints->SetNumberOfPoints(maxIndex+1); pointScalars->SetNumberOfTuples(maxIndex+1); cellScalars->SetNumberOfTuples(mesh->GetNumberOfCells()); vtkFloatingPointType vtkpoint[3]; typename MeshType::PointType itkPhysicalPoint; if (geometryFrame == 0) { for(i = points->Begin(); i != points->End(); ++i) { // Get the point index from the point container iterator int idx = i->Index(); itkPhysicalPoint = i->Value(); mitk::itk2vtk(itkPhysicalPoint, vtkpoint); // Set the vtk point at the index with the the coord array from itk vpoints->SetPoint(idx, vtkpoint); if(usePointScalarAccessor) { pointScalars->InsertTuple1( idx, ScalarAccessor::GetPointScalar( mesh->GetPointData(), i->Index(), mesh, pointDataType ) ); } } } else { mitk::Point3D mitkWorldPoint; for(i = points->Begin(); i != points->End(); ++i) { // Get the point index from the point container iterator int idx = i->Index(); itkPhysicalPoint = i->Value(); geometryFrame->ItkPhysicalPointToWorld(itkPhysicalPoint, mitkWorldPoint); mitk::itk2vtk(mitkWorldPoint, vtkpoint); // Set the vtk point at the index with the the coord array from itk vpoints->SetPoint(idx, vtkpoint); if(usePointScalarAccessor) { pointScalars->InsertTuple1( idx, ScalarAccessor::GetPointScalar( mesh->GetPointData(), i->Index(), mesh, pointDataType ) ); } } } // Set the points on the vtk grid vgrid->SetPoints(vpoints); if (usePointScalarAccessor) vgrid->GetPointData()->SetScalars(pointScalars); // Now create the cells using the MultiVisitor // 1. Create a MultiVisitor typename MeshType::CellType::MultiVisitor::Pointer mv = MeshType::CellType::MultiVisitor::New(); // 2. Create visitors typename SingleCellArrayLineVisitor::Pointer lv = SingleCellArrayLineVisitor::New(); typename SingleCellArrayPolygonVisitor::Pointer pv = SingleCellArrayPolygonVisitor::New(); typename SingleCellArrayTriangleVisitor::Pointer tv = SingleCellArrayTriangleVisitor::New(); typename SingleCellArrayQuadrilateralVisitor::Pointer qv = SingleCellArrayQuadrilateralVisitor::New(); typename SingleCellArrayTetrahedronVisitor::Pointer tetv = SingleCellArrayTetrahedronVisitor::New(); typename SingleCellArrayHexahedronVisitor::Pointer hv = SingleCellArrayHexahedronVisitor::New(); // 3. Set up the visitors //int vtkCellCount = 0; // running counter for current cell being inserted into vtk int numCells = mesh->GetNumberOfCells(); int *types = new int[numCells]; // type array for vtk // create vtk cells and estimate the size vtkCellArray* cells = vtkCellArray::New(); cells->Allocate(numCells); // Set the TypeArray CellCount and CellArray for the visitors lv->SetTypeArray(types); lv->SetCellArray(cells); pv->SetTypeArray(types); pv->SetCellArray(cells); tv->SetTypeArray(types); //tv->SetCellCounter(&vtkCellCount); tv->SetCellArray(cells); qv->SetTypeArray(types); //qv->SetCellCounter(&vtkCellCount); qv->SetCellArray(cells); tetv->SetTypeArray(types); tetv->SetCellArray(cells); hv->SetTypeArray(types); hv->SetCellArray(cells); if (useCellScalarAccessor) { lv->SetUseCellScalarAccessor(true); lv->SetCellScalars(cellScalars); lv->SetMeshCellData(mesh->GetCellData()); pv->SetUseCellScalarAccessor(true); pv->SetCellScalars(cellScalars); pv->SetMeshCellData(mesh->GetCellData()); tv->SetUseCellScalarAccessor(true); tv->SetCellScalars(cellScalars); tv->SetMeshCellData(mesh->GetCellData()); qv->SetUseCellScalarAccessor(true); qv->SetCellScalars(cellScalars); qv->SetMeshCellData(mesh->GetCellData()); tetv->SetUseCellScalarAccessor(true); tetv->SetCellScalars(cellScalars); tetv->SetMeshCellData(mesh->GetCellData()); hv->SetUseCellScalarAccessor(true); hv->SetCellScalars(cellScalars); hv->SetMeshCellData(mesh->GetCellData()); } // add the visitors to the multivisitor mv->AddVisitor(lv); mv->AddVisitor(pv); mv->AddVisitor(tv); mv->AddVisitor(qv); mv->AddVisitor(tetv); mv->AddVisitor(hv); // Now ask the mesh to accept the multivisitor which // will Call Visit for each cell in the mesh that matches the // cell types of the visitors added to the MultiVisitor mesh->Accept(mv); // Now set the cells on the vtk grid with the type array and cell array vgrid->SetCells(types, cells); vgrid->GetCellData()->SetScalars(cellScalars); // Clean up vtk objects (no vtkSmartPointer ... ) cells->Delete(); vpoints->Delete(); delete[] types; pointScalars->Delete(); cellScalars->Delete(); //MITK_INFO << "meshToUnstructuredGrid end" << std::endl; return vgrid; } /*! create a vtkPolyData object from an itkMesh */ static vtkPolyData* MeshToPolyData(MeshType* mesh, bool onlyTriangles = false, bool useScalarAccessor = false, unsigned int pointDataType = 0, mitk::Geometry3D* geometryFrame=NULL, vtkPolyData* polydata = NULL) { /*! default Distribute line cell visitior definition */ typedef typename itk::CellInterfaceVisitorImplementation, DistributeUserVisitorType> DistributeLineVisitor; /*! default Distribute polygon cell visitior definition */ typedef typename itk::CellInterfaceVisitorImplementation, DistributeUserVisitorType> DistributePolygonVisitor; /*! default Distribute triangle cell visitior definition */ typedef typename itk::CellInterfaceVisitorImplementation >, DistributeUserVisitorType> DistributeTriangleVisitor; /*! default Distribute quad cell visitior definition */ typedef typename itk::CellInterfaceVisitorImplementation >, DistributeUserVisitorType> DistributeQuadrilateralVisitor; /*! default Distribute triangle cell visitior definition */ typedef typename itk::CellInterfaceVisitorImplementation >, ExactUserVisitorType> ExactTriangleVisitor; // Get the number of points in the mesh int numPoints = mesh->GetNumberOfPoints(); if(numPoints == 0) { //mesh->Print(std::cerr); MITK_ERROR << "no points in Grid " << std::endl; } // Create a vtkPolyData if(polydata == NULL) polydata = vtkPolyData::New(); else polydata->Initialize(); // Create the vtkPoints object and set the number of points vtkPoints* vpoints = vtkPoints::New( VTK_DOUBLE ); vtkFloatArray * scalars = vtkFloatArray::New(); scalars->SetNumberOfComponents(1); typename MeshType::PointsContainer::Pointer points = mesh->GetPoints(); typename MeshType::PointsContainer::Iterator i; // iterate over all the points in the itk mesh to find // the maximal index unsigned int maxIndex = 0; for(i = points->Begin(); i != points->End(); ++i) { if(maxIndex < i->Index()) maxIndex = i->Index(); } // initialize vtk-classes for points and scalars vpoints->SetNumberOfPoints(maxIndex+1); scalars->SetNumberOfTuples(maxIndex+1); // iterate over all the points in the itk mesh filling in // the vtkPoints object as we go vtkFloatingPointType vtkpoint[3]; typename MeshType::PointType itkPhysicalPoint; if(geometryFrame==NULL) { for(i = points->Begin(); i != points->End(); ++i) { // Get the point index from the point container iterator int idx = i->Index(); itkPhysicalPoint = i->Value(); mitk::itk2vtk(itkPhysicalPoint, vtkpoint); // Set the vtk point at the index with the the coord array from itk // itk returns a const pointer, but vtk is not const correct, so // we have to use a const cast to get rid of the const // vpoints->SetPoint(idx, const_cast(i->Value().GetDataPointer())); vpoints->SetPoint(idx, vtkpoint); if(useScalarAccessor) { scalars->InsertTuple1( idx, ScalarAccessor::GetPointScalar( mesh->GetPointData(), i->Index(), mesh, pointDataType ) ); } } } else { mitk::Point3D mitkWorldPoint; for(i = points->Begin(); i != points->End(); ++i) { // Get the point index from the point container iterator int idx = i->Index(); itkPhysicalPoint = i->Value(); geometryFrame->ItkPhysicalPointToWorld(itkPhysicalPoint, mitkWorldPoint); mitk::itk2vtk(mitkWorldPoint, vtkpoint); // Set the vtk point at the index with the the coord array from itk // itk returns a const pointer, but vtk is not const correct, so // we have to use a const cast to get rid of the const // vpoints->SetPoint(idx, const_cast(i->Value().GetDataPointer())); vpoints->SetPoint(idx, vtkpoint); if(useScalarAccessor) { scalars->InsertTuple1( idx, ScalarAccessor::GetPointScalar( mesh->GetPointData(), i->Index(), mesh, pointDataType ) ); } } } // Set the points on the vtk grid polydata->SetPoints(vpoints); if (useScalarAccessor) polydata->GetPointData()->SetScalars(scalars); polydata->GetPointData()->CopyAllOn(); // Now create the cells using the MulitVisitor // 1. Create a MultiVisitor typedef typename MeshType::CellType::MultiVisitor MeshMV; typename MeshMV::Pointer mv = MeshMV::New(); int numCells = mesh->GetNumberOfCells(); if (onlyTriangles) { // create vtk cells and allocate vtkCellArray* trianglecells = vtkCellArray::New(); trianglecells->Allocate(numCells); // 2. Create a triangle visitor and add it to the multivisitor typename ExactTriangleVisitor::Pointer tv = ExactTriangleVisitor::New(); tv->SetCellArrays(NULL, trianglecells, NULL, NULL); mv->AddVisitor(tv); // 3. Now ask the mesh to accept the multivisitor which // will Call Visit for each cell in the mesh that matches the // cell types of the visitors added to the MultiVisitor mesh->Accept(mv); // 4. Set the result into our vtkPolyData if(trianglecells->GetNumberOfCells()>0) polydata->SetStrips(trianglecells); // 5. Clean up vtk objects (no vtkSmartPointer ... ) trianglecells->Delete(); } else { // create vtk cells and allocate vtkCellArray* linecells = vtkCellArray::New(); vtkCellArray* trianglecells = vtkCellArray::New(); vtkCellArray* polygoncells = vtkCellArray::New(); linecells->Allocate(numCells); trianglecells->Allocate(numCells); polygoncells->Allocate(numCells); // 2. Create visitors typename DistributeLineVisitor::Pointer lv = DistributeLineVisitor::New(); typename DistributePolygonVisitor::Pointer pv = DistributePolygonVisitor::New(); typename DistributeTriangleVisitor::Pointer tv = DistributeTriangleVisitor::New(); typename DistributeQuadrilateralVisitor::Pointer qv = DistributeQuadrilateralVisitor::New(); lv->SetCellArrays(linecells, trianglecells, polygoncells, polygoncells); pv->SetCellArrays(linecells, trianglecells, polygoncells, polygoncells); tv->SetCellArrays(linecells, trianglecells, polygoncells, polygoncells); qv->SetCellArrays(linecells, trianglecells, polygoncells, polygoncells); // add the visitors to the multivisitor mv->AddVisitor(tv); mv->AddVisitor(lv); mv->AddVisitor(pv); mv->AddVisitor(qv); // 3. Now ask the mesh to accept the multivisitor which // will Call Visit for each cell in the mesh that matches the // cell types of the visitors added to the MultiVisitor mesh->Accept(mv); // 4. Set the result into our vtkPolyData if(linecells->GetNumberOfCells()>0) polydata->SetLines(linecells); if(trianglecells->GetNumberOfCells()>0) polydata->SetStrips(trianglecells); if(polygoncells->GetNumberOfCells()>0) polydata->SetPolys(polygoncells); // 5. Clean up vtk objects (no vtkSmartPointer ... ) linecells->Delete(); trianglecells->Delete(); polygoncells->Delete(); } vpoints->Delete(); scalars->Delete(); //MITK_INFO << "meshToPolyData end" << std::endl; return polydata; } static typename MeshType::Pointer CreateRegularSphereMesh(typename MeshType::PointType center, typename MeshType::PointType::VectorType scale, int resolution) { typedef itk::RegularSphereMeshSource SphereSourceType; typename SphereSourceType::Pointer mySphereSource = SphereSourceType::New(); mySphereSource->SetCenter(center); mySphereSource->SetScale(scale); mySphereSource->SetResolution( resolution ); mySphereSource->Update(); typename MeshType::Pointer resultMesh = mySphereSource->GetOutput(); resultMesh->Register(); // necessary ???? return resultMesh; } static typename MeshType::Pointer CreateSphereMesh(typename MeshType::PointType center, typename MeshType::PointType scale, int* resolution) { typedef typename itk::SphereMeshSource SphereSource; typename SphereSource::Pointer mySphereSource = SphereSource::New(); mySphereSource->SetCenter(center); mySphereSource->SetScale(scale); mySphereSource->SetResolutionX(resolution[0]); mySphereSource->SetResolutionY(resolution[1]); mySphereSource->SetSquareness1(1); mySphereSource->SetSquareness2(1); mySphereSource->Update(); mySphereSource->GetOutput(); typename MeshType::Pointer resultMesh = mySphereSource->GetOutput(); resultMesh->Register(); return resultMesh; } // static typename MeshType::Pointer TranslateMesh(typename MeshType::PointType vec, MeshType* input) // { // // typename MeshType::Pointer output = MeshType::New(); // { // output->SetPoints(input->GetPoints()); // output->SetPointData(input->GetPointData()); // output->SetCells(input->GetCells()); // output->SetLastCellId( input->GetLastCellId() ); // typename MeshType::GeometryMapIterator pointDataIterator = input->GetGeometryData()->Begin(); // typename MeshType::GeometryMapIterator pointDataEnd = input->GetGeometryData()->End(); // // typename MeshType::PointType inputPoint,outputPoint; // // while (pointDataIterator != pointDataEnd) // { // unsigned long pointId = pointDataIterator->Index(); // itk::SimplexMeshGeometry* newGeometry = new itk::SimplexMeshGeometry(); // itk::SimplexMeshGeometry* refGeometry = pointDataIterator->Value(); // // input->GetPoint(pointId, &inputPoint ); // outputPoint[0] = inputPoint[0] + vec[0]; // outputPoint[1] = inputPoint[1] + vec[1]; // outputPoint[2] = inputPoint[2] + vec[2]; // output->SetPoint( pointId, outputPoint ); // // // newGeometry->pos = outputPoint; // newGeometry->neighborIndices = refGeometry->neighborIndices; // newGeometry->meanCurvature = refGeometry->meanCurvature; // newGeometry->neighbors = refGeometry->neighbors; // newGeometry->oldPos = refGeometry->oldPos; // newGeometry->eps = refGeometry->eps; // newGeometry->referenceMetrics = refGeometry->referenceMetrics; // newGeometry->neighborSet = refGeometry->neighborSet; // newGeometry->distance = refGeometry->distance; // newGeometry->externalForce = refGeometry->externalForce; // newGeometry->internalForce = refGeometry->internalForce; // output->SetGeometryData(pointId, newGeometry); // pointDataIterator++; // } // } //// output->SetGeometryData( inputMesh->GetGeometryData() ); // return output; // } static typename MeshType::Pointer CreateRegularSphereMesh2(typename MeshType::PointType center, typename MeshType::PointType scale, int resolution) { typedef typename itk::AutomaticTopologyMeshSource MeshSourceType; typename MeshSourceType::Pointer mySphereSource = MeshSourceType::New(); typename MeshType::PointType pnt0, pnt1, pnt2, pnt3, pnt4, pnt5, pnt6, pnt7, pnt8, pnt9, pnt10, pnt11; double c1= 0.5 * (1.0 + sqrt(5.0)); double c2= 1.0; double len = sqrt( c1*c1 + c2*c2 ); c1 /= len; c2 /= len; pnt0[0] = center[0] - c1*scale[0]; pnt0[1] = center[1]; pnt0[2] = center[2] + c2*scale[2]; pnt1[0] = center[0]; pnt1[1] = center[1] + c2*scale[1]; pnt1[2] = center[2] - c1*scale[2]; pnt2[0] = center[0]; pnt2[1] = center[1] + c2*scale[1]; pnt2[2] = center[2] + c1*scale[2]; pnt3[0] = center[0] + c1*scale[0]; pnt3[1] = center[1]; pnt3[2] = center[2] - c2*scale[2]; pnt4[0] = center[0] - c2*scale[0]; pnt4[1] = center[1] - c1*scale[1]; pnt4[2] = center[2]; pnt5[0] = center[0] - c2*scale[0]; pnt5[1] = center[1] + c1*scale[1]; pnt5[2] = center[2]; pnt6[0] = center[0]; pnt6[1] = center[1] - c2*scale[1]; pnt6[2] = center[2] + c1*scale[2]; pnt7[0] = center[0] + c2*scale[0]; pnt7[1] = center[1] + c1*scale[1]; pnt7[2] = center[2]; pnt8[0] = center[0]; pnt8[1] = center[1] - c2*scale[1]; pnt8[2] = center[2] - c1*scale[2]; pnt9[0] = center[0] + c1*scale[0]; pnt9[1] = center[1]; pnt9[2] = center[2] + c2*scale[2]; pnt10[0]= center[0] + c2*scale[0]; pnt10[1]= center[1] - c1*scale[1]; pnt10[2]= center[2]; pnt11[0]= center[0] - c1*scale[0]; pnt11[1]= center[1]; pnt11[2]= center[2] - c2*scale[2]; addTriangle( mySphereSource, scale, pnt9, pnt2, pnt6, resolution ); addTriangle( mySphereSource, scale, pnt1, pnt11, pnt5, resolution ); addTriangle( mySphereSource, scale, pnt11, pnt1, pnt8, resolution ); addTriangle( mySphereSource, scale, pnt0, pnt11, pnt4, resolution ); addTriangle( mySphereSource, scale, pnt3, pnt1, pnt7, resolution ); addTriangle( mySphereSource, scale, pnt3, pnt8, pnt1, resolution ); addTriangle( mySphereSource, scale, pnt9, pnt3, pnt7, resolution ); addTriangle( mySphereSource, scale, pnt0, pnt6, pnt2, resolution ); addTriangle( mySphereSource, scale, pnt4, pnt10, pnt6, resolution ); addTriangle( mySphereSource, scale, pnt1, pnt5, pnt7, resolution ); addTriangle( mySphereSource, scale, pnt7, pnt5, pnt2, resolution ); addTriangle( mySphereSource, scale, pnt8, pnt3, pnt10, resolution ); addTriangle( mySphereSource, scale, pnt4, pnt11, pnt8, resolution ); addTriangle( mySphereSource, scale, pnt9, pnt7, pnt2, resolution ); addTriangle( mySphereSource, scale, pnt10, pnt9, pnt6, resolution ); addTriangle( mySphereSource, scale, pnt0, pnt5, pnt11, resolution ); addTriangle( mySphereSource, scale, pnt0, pnt2, pnt5, resolution ); addTriangle( mySphereSource, scale, pnt8, pnt10, pnt4, resolution ); addTriangle( mySphereSource, scale, pnt3, pnt9, pnt10, resolution ); addTriangle( mySphereSource, scale, pnt6, pnt0, pnt4, resolution ); return mySphereSource->GetOutput(); } private: static void addTriangle( typename itk::AutomaticTopologyMeshSource::Pointer meshSource, typename MeshType::PointType scale, typename MeshType::PointType pnt0, typename MeshType::PointType pnt1, typename MeshType::PointType pnt2, int resolution ) { if (resolution==0) { // add triangle meshSource->AddTriangle( meshSource->AddPoint( pnt0 ), meshSource->AddPoint( pnt1 ), meshSource->AddPoint( pnt2 ) ); } else { vnl_vector_fixed v1, v2, res, pv; v1 = (pnt1-pnt0).Get_vnl_vector(); v2 = (pnt2-pnt0).Get_vnl_vector(); res = vnl_cross_3d( v1, v2 ); pv = pnt0.GetVectorFromOrigin().Get_vnl_vector(); //double d = res[0]*pv[0] + res[1]*pv[1] + res[2]*pv[2]; // subdivision typename MeshType::PointType pnt01, pnt12, pnt20; for (int d=0; d<3; d++) { pnt01[d] = (pnt0[d] + pnt1[d]) / 2.0; pnt12[d] = (pnt1[d] + pnt2[d]) / 2.0; pnt20[d] = (pnt2[d] + pnt0[d]) / 2.0; } // map new points to sphere double lenPnt01=0; for (int d=0; d<3; d++) lenPnt01 += pnt01[d]*pnt01[d]; lenPnt01 = sqrt( lenPnt01 ); double lenPnt12=0; for (int d=0; d<3; d++) lenPnt12 += pnt12[d]*pnt12[d]; lenPnt12 = sqrt( lenPnt12 ); double lenPnt20=0; for (int d=0; d<3; d++) lenPnt20 += pnt20[d]*pnt20[d]; lenPnt20 = sqrt( lenPnt20 ); for (int d=0; d<3; d++) { pnt01[d] *= scale[d]/lenPnt01; pnt12[d] *= scale[d]/lenPnt12; pnt20[d] *= scale[d]/lenPnt20; } addTriangle( meshSource, scale, pnt0, pnt01, pnt20, resolution-1 ); addTriangle( meshSource, scale, pnt01, pnt1, pnt12, resolution-1 ); addTriangle( meshSource, scale, pnt20, pnt12, pnt2, resolution-1 ); addTriangle( meshSource, scale, pnt01, pnt12, pnt20, resolution-1 ); } } }; #endif // MITKMESHUTIL_H_INCLUDED