mitkCylindricToCartesianFilter.cpp
Go to the documentation of this file.
00001 /*========================================================================= 00002 00003 Program: Medical Imaging & Interaction Toolkit 00004 Language: C++ 00005 Date: $Date$ 00006 Version: $Revision$ 00007 00008 Copyright (c) German Cancer Research Center, Division of Medical and 00009 Biological Informatics. All rights reserved. 00010 See MITKCopyright.txt or https://www.mitk.org/copyright.html for details. 00011 00012 This software is distributed WITHOUT ANY WARRANTY; without even 00013 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 00014 PURPOSE. See the above copyright notices for more information. 00015 00016 =========================================================================*/ 00017 00018 00019 #include "mitkCylindricToCartesianFilter.h" 00020 #include "mitkImageTimeSelector.h" 00021 #include "mitkSlicedGeometry3D.h" 00022 #include "mitkPlaneGeometry.h" 00023 #include "mitkProperties.h" 00024 00025 #include <ipPic/mitkIpPicTypeMultiplex.h> 00026 00027 template <class T> 00028 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 00029 { 00030 T outsideValue = static_cast<T>(_outsideValue); 00031 register float f, ft, f0, f1, f2, f3; 00032 00033 mitkIpInt2_t ox_size; 00034 mitkIpInt2_t nx_size, ny_size, nz_size; 00035 int oxy_size, nxy_size; 00036 00037 T* orig, *dp, *dest_start; 00038 00039 mitkIpInt2_t* coneCutOff=(mitkIpInt2_t*)coneCutOff_pic->data; 00040 00041 orig=(T*)pic->data; 00042 ox_size=pic->n[0]; 00043 oxy_size=ox_size*pic->n[1]; 00044 00045 nx_size=dest->n[0]; 00046 ny_size=dest->n[1]; 00047 nxy_size=nx_size*ny_size; 00048 nz_size=dest->n[2]; 00049 00050 /*nx_size=360; 00051 ny_size=360; 00052 nxy_size=nx_size*ny_size; 00053 nz_size=256;*/ 00054 00055 dest_start=dp=((T*)dest->data)+nxy_size*(nz_size-1); 00056 00057 mitkIpInt2_t y; 00058 // int size=_mitkIpPicElements(pic); 00059 00060 register mitkIpInt2_t x,z; 00061 for(y=0;y<ny_size;++y) 00062 { 00063 mitkIpInt2_t x_start; 00064 register mitkIpInt2_t x_end; 00065 int r0plusphi0=*rt; 00066 if(r0plusphi0>=0) 00067 { 00068 x_start=0; x_end=nx_size; 00069 } 00070 else 00071 { 00072 x_start=-r0plusphi0; x_end=nx_size+r0plusphi0; 00073 00074 for(z=0;z<nz_size;++z,dp+=-nxy_size-x_start) 00075 for(x=0;x<x_start;++x,++dp) 00076 *dp=outsideValue; 00077 dp+=nxy_size*nz_size; dp+=x_end; 00078 for(z=0;z<nz_size;++z,dp+=-nxy_size-x_start) 00079 for(x=x_end;x<nx_size;++x,++dp) 00080 *dp=outsideValue; 00081 dp+=nxy_size*nz_size; dp-=x_end; 00082 00083 fr+=x_start; 00084 fphi+=x_start; 00085 rt+=x_start; 00086 phit+=x_start; 00087 dp+=x_start; 00088 coneCutOff+=x_start; 00089 } 00090 00091 for(x=x_start;x<x_end;++x, ++fr, ++fphi, ++rt, ++phit, ++dp, ++coneCutOff) 00092 { 00093 f=*fr; 00094 f3=*fphi; 00095 f0=1-f; 00096 f1=1-f3; 00097 f2=f0*f3; 00098 f0*=f1; 00099 f1*=f; 00100 f3*=f; 00101 00102 r0plusphi0=*rt+*phit; 00103 00104 mitkIpInt2_t z_start; 00105 00106 z_start=*coneCutOff; 00107 00108 T *op; 00109 register unsigned int step; 00110 unsigned int *ztp; 00111 00112 op=orig+*rt+*phit; 00113 00114 float *fzp; 00115 00116 for(z=0;z<z_start;++z, dp-=nxy_size) 00117 *dp=outsideValue; 00118 00119 ztp=zt+z_start; 00120 fzp=fz+z_start; 00121 00122 for(z=z_start;z<nz_size;++z, dp-=nxy_size) 00123 { 00124 step=*(ztp++); 00125 00126 register T *opt=op; 00127 opt+=step; 00128 f =*opt*f0; ++opt; 00129 f+=*opt*f1; opt+=oxy_size; --opt; 00130 f+=*opt*f2; ++opt; 00131 f+=*opt*f3; opt-=oxy_size; --opt; 00132 00133 opt+=ox_size; 00134 ft =*opt*f0; ++opt; 00135 ft+=*opt*f1; opt+=oxy_size; --opt; 00136 ft+=*opt*f2; ++opt; 00137 ft+=*opt*f3; 00138 00139 *dp=(T)((1-*fzp)*f+*fzp*ft+0.5); 00140 } 00141 00142 dp+=nxy_size*nz_size; 00143 } 00144 fr+=nx_size-x_end; 00145 fphi+=nx_size-x_end; 00146 rt+=nx_size-x_end; 00147 phit+=nx_size-x_end; 00148 dp+=nx_size-x_end; 00149 coneCutOff+=nx_size-x_end; 00150 } 00151 } 00152 00153 void mitk::CylindricToCartesianFilter::buildTransformShortCuts(int orig_xsize, int orig_ysize, int orig_zsize, int new_xsize, mitkIpPicDescriptor * &rt_pic, mitkIpPicDescriptor * &phit_pic, mitkIpPicDescriptor * &fr_pic, mitkIpPicDescriptor * &fphi_pic, unsigned int * &zt, float * &fz) 00154 { 00155 --orig_zsize; 00156 00157 rt_pic=mitkIpPicNew(); 00158 rt_pic->type=mitkIpPicInt; 00159 rt_pic->bpe=16; 00160 rt_pic->dim=2; 00161 rt_pic->n[0]=rt_pic->n[1]=new_xsize; 00162 rt_pic->data=malloc(_mitkIpPicSize(rt_pic)); 00163 00164 phit_pic=mitkIpPicNew(); 00165 phit_pic->type=mitkIpPicUInt; 00166 phit_pic->bpe=32; 00167 phit_pic->dim=2; 00168 phit_pic->n[0]=phit_pic->n[1]=new_xsize; 00169 phit_pic->data=malloc(_mitkIpPicSize(phit_pic)); 00170 00171 fr_pic=mitkIpPicNew(); 00172 fr_pic->type=mitkIpPicFloat; 00173 fr_pic->bpe=32; 00174 fr_pic->dim=2; 00175 fr_pic->n[0]=fr_pic->n[1]=new_xsize; 00176 fr_pic->data=malloc(_mitkIpPicSize(fr_pic)); 00177 00178 fphi_pic=mitkIpPicNew(); 00179 fphi_pic->type=mitkIpPicFloat; 00180 fphi_pic->bpe=32; 00181 fphi_pic->dim=2; 00182 fphi_pic->n[0]=fphi_pic->n[1]=new_xsize; 00183 fphi_pic->data=malloc(_mitkIpPicSize(fphi_pic)); 00184 00185 mitkIpInt2_t *rtp=(mitkIpInt2_t*)rt_pic->data, *rt_xzero, rt, phit; 00186 mitkIpUInt4_t *phitp=(mitkIpUInt4_t*)phit_pic->data; 00187 mitkIpFloat4_t *fr=(mitkIpFloat4_t *)fr_pic->data; 00188 mitkIpFloat4_t *fphi=(mitkIpFloat4_t *)fphi_pic->data; 00189 00190 mitkIpFloat4_t r, phi, scale=(double)orig_xsize/(double)new_xsize; 00191 00192 int x,y,xy0,xy0_orig, oxy_size, new_zsize; 00193 oxy_size=orig_xsize*orig_ysize; 00194 xy0=(int)(((double)new_xsize)/2+0.5); 00195 xy0_orig=(int)(((double)orig_xsize)/2+0.5); 00196 00197 new_zsize=(int)(orig_ysize/scale); 00198 // \bug y compared to x 00199 for(y=0;y<new_xsize;++y) 00200 { 00201 rt_xzero=rtp; *rtp=0; 00202 for(x=0;x<new_xsize;++x,++fr,++fphi,++rtp, ++phitp) 00203 { 00204 int xq=x-xy0, yq=y-xy0; 00205 00206 r=sqrt( (double) (xq*xq+yq*yq)); 00207 00208 // float rtest=-(xy0-sqrt(xy0*xy0-yq*yq))-0.5; 00209 rt=(mitkIpInt2_t)(-(xy0-sqrt((double) (xy0*xy0-yq*yq)))-0.5);/*in rt steht der Index des Endes der zu �berspringenden Punkte=>anfangen bei -rt+1!*/ 00210 // if((x>=-rt) && (x<new_xsize+rt)) 00211 { 00212 if(y!=xy0) 00213 r=r*(y>xy0?1.0:-1.0)*scale+xy0_orig; 00214 else 00215 r=r*(x>xy0?-1.0:1.0)*scale+xy0_orig; 00216 rt=(mitkIpInt2_t)r; 00217 int xtmp=x; 00218 if(x>xy0) 00219 xtmp=new_xsize-x; 00220 if(rt<0) 00221 { 00222 r=rt=0; 00223 if(xtmp>-*rt_xzero) 00224 *rt_xzero=-xtmp; 00225 *fr=0; 00226 } 00227 else 00228 if(rt>orig_xsize-1) 00229 { 00230 r=rt=orig_xsize-1; 00231 if(xtmp>-*rt_xzero) 00232 *rt_xzero=-xtmp; 00233 *fr=0; 00234 } 00235 else 00236 *fr=r-rt; 00237 if(*fr<0) 00238 *fr=0; 00239 } 00240 // else 00241 // *fr=0; 00242 00243 phi=orig_zsize-(yq==0?1:-atan((float)xq/yq)/M_PI+0.5)*orig_zsize; 00244 phit=(mitkIpUInt4_t)phi; 00245 *fphi=phi-phit; 00246 00247 *rtp=rt; 00248 *phitp=phit*oxy_size; 00249 } 00250 } 00251 00252 zt=(unsigned int *)malloc(sizeof(unsigned int)*new_zsize); 00253 fz=(float *)malloc(sizeof(float)*new_zsize); 00254 00255 float *fzp=fz; 00256 unsigned int *ztp=zt; 00257 00258 int z; 00259 float z_step=orig_ysize/(orig_ysize*((float)new_xsize)/orig_xsize); 00260 for(z=0;z<new_zsize;++z,++fzp,++ztp) 00261 { 00262 *fzp=z*z_step; 00263 *ztp=(unsigned int)*fzp; 00264 *fzp-=*ztp; 00265 *ztp*=orig_xsize; 00266 } 00267 } 00268 00269 void mitk::CylindricToCartesianFilter::buildConeCutOffShortCut(int orig_xsize, int orig_ysize, mitkIpPicDescriptor *rt_pic, mitkIpPicDescriptor *fr_pic, float a, float b, mitkIpPicDescriptor * &coneCutOff_pic) 00270 { 00271 coneCutOff_pic=mitkIpPicNew(); 00272 coneCutOff_pic->type=mitkIpPicInt; 00273 coneCutOff_pic->bpe=16; 00274 coneCutOff_pic->dim=2; 00275 coneCutOff_pic->n[0]=coneCutOff_pic->n[1]=rt_pic->n[0]; 00276 coneCutOff_pic->data=malloc(_mitkIpPicSize(coneCutOff_pic)); 00277 00278 int i, size=_mitkIpPicElements(rt_pic); 00279 mitkIpInt2_t *rt, *ccop, ohx_size, nz_size; 00280 mitkIpFloat4_t *fr; 00281 00282 a*=(float)rt_pic->n[0]/orig_xsize; 00283 b*=(float)rt_pic->n[0]/orig_xsize; 00284 00285 ohx_size=orig_xsize/2; 00286 nz_size=orig_ysize*rt_pic->n[0]/orig_xsize; 00287 00288 rt=(mitkIpInt2_t *)rt_pic->data; fr=(mitkIpFloat4_t*)fr_pic->data; ccop=(mitkIpInt2_t *)coneCutOff_pic->data; 00289 00290 for(i=0; i<size; ++i, ++rt, ++ccop) 00291 { 00292 register mitkIpInt2_t cco; 00293 if(*rt<=ohx_size) 00294 cco=(mitkIpInt2_t)(a*(*rt+*fr)+b); 00295 else 00296 cco=(mitkIpInt2_t)(a*(orig_xsize-(*rt+*fr))+b); 00297 if(cco<0) 00298 cco=0; 00299 if(cco>=nz_size) 00300 cco=nz_size; 00301 *ccop=cco; 00302 } 00303 } 00304 00305 void mitk::CylindricToCartesianFilter::GenerateOutputInformation() 00306 { 00307 mitk::Image::Pointer output = this->GetOutput(); 00308 if ((output->IsInitialized()) && (output->GetPipelineMTime() <= m_TimeOfHeaderInitialization.GetMTime())) 00309 return; 00310 00311 mitk::Image::ConstPointer input = this->GetInput(); 00312 00313 itkDebugMacro(<<"GenerateOutputInformation()"); 00314 00315 unsigned int i, *tmpDimensions=new unsigned int[std::max(3u,input->GetDimension())]; 00316 00317 tmpDimensions[0]=m_TargetXSize; 00318 if(tmpDimensions[0]==0) 00319 tmpDimensions[0] = input->GetDimension(0); 00320 00321 float scale=((float)tmpDimensions[0])/input->GetDimension(0); 00322 00323 tmpDimensions[1] = tmpDimensions[0]; 00324 tmpDimensions[2] = (unsigned int)(scale*input->GetDimension(1)); 00325 00326 for(i=3;i<input->GetDimension();++i) 00327 tmpDimensions[i]=input->GetDimension(i); 00328 00329 output->Initialize(input->GetPixelType(), 00330 input->GetDimension(), 00331 tmpDimensions, 00332 input->GetNumberOfChannels()); 00333 00334 // initialize the spacing of the output 00335 Vector3D spacing = input->GetSlicedGeometry()->GetSpacing(); 00336 if(input->GetDimension()>=2) 00337 spacing[2]=spacing[1]; 00338 else 00339 spacing[2] = 1.0; 00340 spacing[1] = spacing[0]; 00341 spacing *= 1.0/scale; 00342 output->GetSlicedGeometry()->SetSpacing(spacing); 00343 00344 mitk::Point3iProperty::Pointer pointProp; 00345 pointProp = dynamic_cast<mitk::Point3iProperty*>(input->GetProperty("ORIGIN").GetPointer()); 00346 if (pointProp.IsNotNull() ) 00347 { 00348 itk::Point<int, 3> tp = pointProp->GetValue(); 00349 tp[2] = (int)(tmpDimensions[2]-tp[1] * scale-1); 00350 tp[0] = tmpDimensions[0]/2; 00351 tp[1] = tmpDimensions[0]/2; 00352 mitk::Point3iProperty::Pointer pointProp = mitk::Point3iProperty::New(tp); 00353 output->SetProperty("ORIGIN", pointProp); 00354 } 00355 delete [] tmpDimensions; 00356 00357 //output->GetSlicedGeometry()->SetGeometry2D(mitk::Image::BuildStandardPlaneGeometry2D(output->GetSlicedGeometry(), tmpDimensions).GetPointer(), 0); 00358 //set the timebounds - after SetGeometry2D, so that the already created PlaneGeometry will also receive this timebounds. 00359 //@fixme!!! will not work for not evenly timed data! 00360 output->GetSlicedGeometry()->SetTimeBounds(input->GetSlicedGeometry()->GetTimeBounds()); 00361 output->GetTimeSlicedGeometry()->InitializeEvenlyTimed(output->GetSlicedGeometry(), output->GetTimeSlicedGeometry()->GetTimeSteps()); 00362 00363 output->SetPropertyList(input->GetPropertyList()->Clone()); 00364 00365 m_TimeOfHeaderInitialization.Modified(); 00366 } 00367 00368 void mitk::CylindricToCartesianFilter::GenerateData() 00369 { 00370 mitk::Image::ConstPointer input = this->GetInput(); 00371 mitk::Image::Pointer output = this->GetOutput(); 00372 00373 mitk::ImageTimeSelector::Pointer timeSelector=mitk::ImageTimeSelector::New(); 00374 timeSelector->SetInput(input); 00375 00376 mitkIpPicDescriptor* pic_transformed=NULL; 00377 pic_transformed = mitkIpPicNew(); 00378 pic_transformed->dim=3; 00379 pic_transformed->bpe = output->GetPixelType().GetBpe(); 00380 pic_transformed->type = output->GetPixelType().GetType(); 00381 pic_transformed->n[0] = output->GetDimension(0); 00382 pic_transformed->n[1] = output->GetDimension(1); 00383 pic_transformed->n[2] = output->GetDimension(2); 00384 pic_transformed->data=malloc(_mitkIpPicSize(pic_transformed)); 00385 00386 int nstart, nmax; 00387 int tstart, tmax; 00388 00389 tstart=output->GetRequestedRegion().GetIndex(3); 00390 nstart=output->GetRequestedRegion().GetIndex(4); 00391 00392 tmax=tstart+output->GetRequestedRegion().GetSize(3); 00393 nmax=nstart+output->GetRequestedRegion().GetSize(4); 00394 00395 if(zt==NULL) 00396 { 00397 timeSelector->SetChannelNr(nstart); 00398 timeSelector->SetTimeNr(tstart); 00399 00400 buildTransformShortCuts(input->GetDimension(0),input->GetDimension(1), input->GetDimension(2), output->GetDimension(0), rt_pic, phit_pic, fr_pic, fphi_pic, zt, fz); 00401 00402 // query the line limiting the sector 00403 a=b=0; 00404 mitk::FloatProperty::Pointer prop; 00405 00406 prop = dynamic_cast<mitk::FloatProperty*>(input->GetProperty("SECTOR LIMITING LINE SLOPE").GetPointer()); 00407 if (prop.IsNotNull() ) 00408 a = prop->GetValue(); 00409 prop = dynamic_cast<mitk::FloatProperty*>(input->GetProperty("SECTOR LIMITING LINE OFFSET").GetPointer()); 00410 if (prop.IsNotNull() ) 00411 b = prop->GetValue(); 00412 00413 buildConeCutOffShortCut(input->GetDimension(0),input->GetDimension(1), rt_pic, fr_pic, a, b, coneCutOff_pic); 00414 // mitkIpPicPut("C:\\temp\\rt_90.pic",rt_pic); 00415 //mitkIpPicPut("C:\\temp\\coneCutOff.pic", coneCutOff_pic); 00416 } 00417 00418 int n,t; 00419 for(n=nstart;n<nmax;++n)//output->GetNumberOfChannels();++n) 00420 { 00421 timeSelector->SetChannelNr(n); 00422 00423 for(t=tstart;t<tmax;++t) 00424 { 00425 timeSelector->SetTimeNr(t); 00426 00427 timeSelector->Update(); 00428 00429 _mitkIpPicFreeTags(pic_transformed->info->tags_head); 00430 pic_transformed->info->tags_head = _mitkIpPicCloneTags(timeSelector->GetOutput()->GetPic()->info->tags_head); 00431 if(input->GetDimension(2)>1) 00432 { 00433 00434 mitkIpPicTypeMultiplex9(_transform, timeSelector->GetOutput()->GetPic(), 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); 00435 // mitkIpPicPut("1trf.pic",pic_transformed); 00436 } 00437 else 00438 { 00439 mitkIpPicDescriptor *doubleSlice=mitkIpPicCopyHeader(timeSelector->GetOutput()->GetPic(), NULL); 00440 doubleSlice->dim=3; 00441 doubleSlice->n[2]=2; 00442 doubleSlice->data=malloc(_mitkIpPicSize(doubleSlice)); 00443 memcpy(doubleSlice->data, timeSelector->GetOutput()->GetPic()->data, _mitkIpPicSize(doubleSlice)/2); 00444 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); 00445 mitkIpPicFree(doubleSlice); 00446 } 00447 output->SetPicVolume(pic_transformed, t, n); 00448 } 00449 } 00450 //mitkIpPicPut("outzzzzzzzz.pic",pic_transformed); 00451 mitkIpPicFree(pic_transformed); 00452 m_TimeOfHeaderInitialization.Modified(); 00453 } 00454 00455 mitk::CylindricToCartesianFilter::CylindricToCartesianFilter() 00456 : m_OutsideValue(0.0), m_TargetXSize(0) 00457 { 00458 rt_pic = NULL; phit_pic = NULL; fr_pic = NULL; fphi_pic = NULL; coneCutOff_pic = NULL; 00459 zt = NULL; fz = NULL; 00460 a=b=0.0; 00461 } 00462 00463 mitk::CylindricToCartesianFilter::~CylindricToCartesianFilter() 00464 { 00465 if(rt_pic!=NULL) mitkIpPicFree(rt_pic); 00466 if(phit_pic!=NULL) mitkIpPicFree(phit_pic); 00467 if(fr_pic!=NULL) mitkIpPicFree(fr_pic); 00468 if(fphi_pic!=NULL) mitkIpPicFree(fphi_pic); 00469 if(coneCutOff_pic!=NULL) mitkIpPicFree(coneCutOff_pic); 00470 if(zt != NULL) free(zt); 00471 if(fz != NULL) free (fz); 00472 } 00473 00474 void mitk::CylindricToCartesianFilter::GenerateInputRequestedRegion() 00475 { 00476 Superclass::GenerateInputRequestedRegion(); 00477 00478 mitk::ImageToImageFilter::InputImagePointer input = 00479 const_cast< mitk::ImageToImageFilter::InputImageType * > ( this->GetInput() ); 00480 mitk::Image::Pointer output = this->GetOutput(); 00481 00482 Image::RegionType requestedRegion; 00483 requestedRegion = output->GetRequestedRegion(); 00484 requestedRegion.SetIndex(0, 0); 00485 requestedRegion.SetIndex(1, 0); 00486 requestedRegion.SetIndex(2, 0); 00487 requestedRegion.SetSize(0, input->GetDimension(0)); 00488 requestedRegion.SetSize(1, input->GetDimension(1)); 00489 requestedRegion.SetSize(2, input->GetDimension(2)); 00490 00491 input->SetRequestedRegion( & requestedRegion ); 00492 }
