mitkDopplerToStrainRateFilter.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 "mitkDopplerToStrainRateFilter.h" 00020 #include "mitkImageTimeSelector.h" 00021 #include "mitkProperties.h" 00022 #include "mitkPlaneGeometry.h" 00023 00024 #include <iostream> 00025 #include <string> 00026 #include <ipFunc/mitkIpFunc.h> 00027 00028 void mitk::DopplerToStrainRateFilter::GenerateOutputInformation() 00029 { 00030 mitk::Image::ConstPointer input = this->GetInput(); 00031 mitk::Image::Pointer output = this->GetOutput(); 00032 00033 if ((output->IsInitialized()) && (this->GetMTime() <= m_TimeOfHeaderInitialization.GetMTime())) 00034 return; 00035 00036 itkDebugMacro(<<"GenerateOutputInformation()"); 00037 00038 unsigned int i; 00039 unsigned int *tmpDimensions = new unsigned int[input->GetDimension()]; 00040 00041 for(i=0;i<input->GetDimension();++i) 00042 tmpDimensions[i]=input->GetDimension(i); 00043 00044 //@todo maybe we should shift the following somehow in ImageToImageFilter 00045 output->Initialize(input->GetPixelType(), 00046 input->GetDimension(), 00047 tmpDimensions, 00048 input->GetNumberOfChannels()); 00049 00050 output->GetSlicedGeometry()->SetSpacing(input->GetSlicedGeometry()->GetSpacing()); 00051 00052 //output->GetSlicedGeometry()->SetGeometry2D(mitk::Image::BuildStandardPlaneGeometry2D(output->GetSlicedGeometry(), tmpDimensions).GetPointer(), 0); 00053 //output->GetSlicedGeometry()->SetEvenlySpaced(); 00054 //set the timebounds - after SetGeometry2D, so that the already created PlaneGeometry will also receive this timebounds. 00055 output->GetSlicedGeometry()->SetTimeBounds(input->GetSlicedGeometry()->GetTimeBounds()); 00056 00057 output->SetPropertyList(input->GetPropertyList()->Clone()); 00058 00059 00060 delete [] tmpDimensions; 00061 00062 m_TimeOfHeaderInitialization.Modified(); 00063 } 00064 00065 void mitk::DopplerToStrainRateFilter::GenerateData() 00066 { 00067 mitk::Image::ConstPointer input = this->GetInput(); 00068 mitk::Image::Pointer output = this->GetOutput(); 00069 00070 00071 mitk::Point3iProperty::Pointer pointProp; 00072 pointProp = dynamic_cast<mitk::Point3iProperty*>(input->GetProperty("ORIGIN").GetPointer()); 00073 if (pointProp.IsNotNull() ) 00074 { 00075 m_Origin = pointProp->GetValue(); 00076 } 00077 00078 MITK_INFO << "compute Strain Rate Image .... " << std::endl 00079 << " origin[0]=" << m_Origin[0] << " origin[1]=" << m_Origin[1] << " origin[2]=" << m_Origin[2] << std::endl 00080 << " distance=" << m_Distance << std::endl 00081 << " NoStrainIntervall=" << m_NoStrainInterval << std::endl; 00082 00083 const Vector3D & spacing = input->GetSlicedGeometry()->GetSpacing(); 00084 // MITK_INFO << " in: xres=" << spacing[0] << " yres=" << spacing[1] << " zres=" << spacing[2] << std::endl; 00085 00086 00087 mitk::ImageTimeSelector::Pointer timeSelector=mitk::ImageTimeSelector::New(); 00088 timeSelector->SetInput(input); 00089 00090 mitkIpPicDescriptor* picStrainRate; 00091 picStrainRate = mitkIpPicNew(); 00092 picStrainRate->dim=3; 00093 picStrainRate->bpe = output->GetPixelType().GetBpe(); 00094 picStrainRate->type = output->GetPixelType().GetType(); 00095 picStrainRate->n[0] = output->GetDimension(0); 00096 picStrainRate->n[1] = output->GetDimension(1); 00097 picStrainRate->n[2] = output->GetDimension(2); 00098 picStrainRate->data=malloc(_mitkIpPicSize(picStrainRate)); 00099 00100 00101 int xDim = picStrainRate->n[0]; 00102 int yDim = picStrainRate->n[1]; 00103 int zDim = picStrainRate->n[2]; 00104 long slice_size = xDim*yDim; 00105 long vol_size = slice_size*zDim; 00106 00107 00108 mitkIpPicDescriptor *picDoppler; 00109 00110 int x,y,z;//,time; // loop-counter 00111 int strainRate; // the computed Strain Rate 00112 int v1,v2; // velocity and Point p1 and p2 00113 00114 float alpha; // the beam-angle, angle betwen current point and beam-point 00115 float dx=0, dy=0; // projection of this->distance to x- and y-axis 00116 int x1; // a square, where the velocity v1 lies in 00117 int y1; // the points are used for interpolation 00118 00119 00120 int minStrainRate=128, maxStrainRate=128; 00121 00122 int n, nmax; 00123 int t, tmax; 00124 00125 t = output->GetRequestedRegion().GetIndex(3); 00126 n = output->GetRequestedRegion().GetIndex(4); 00127 MITK_INFO << "t = " <<t << " n = " << n << std::endl; 00128 00129 tmax = t + output->GetRequestedRegion().GetSize(3); 00130 nmax = n + output->GetRequestedRegion().GetSize(4); 00131 MITK_INFO << "tmax = "<< tmax << " nmax = " << nmax << std::endl; 00132 00133 if (m_Distance<1) m_Distance=1; 00134 00135 for(;n<nmax;n++)//output->GetNumberOfChannels();++n) 00136 { 00137 timeSelector->SetChannelNr(n); 00138 MITK_INFO << "computing chanel n = " << n << std::endl; 00139 00140 for(t=0;t<tmax;t++) 00141 { 00142 MITK_INFO << "computing time slot t = " << t << std::endl; 00143 timeSelector->SetTimeNr(t); 00144 timeSelector->Update(); 00145 00146 _mitkIpPicFreeTags(picStrainRate->info->tags_head); 00147 picStrainRate->info->tags_head = _mitkIpPicCloneTags(timeSelector->GetOutput()->GetPic()->info->tags_head); 00148 00149 //#define WRITE_ANGLE_PIC 00150 #ifdef WRITE_ANGLE_PIC 00151 mitkIpPicDescriptor *anglePic; 00152 mitkIpBool_t isAnglePicWritten = mitkIpFalse; 00153 00154 anglePic = mitkIpPicNew(); 00155 anglePic->type = mitkIpPicInt; 00156 anglePic->bpe = 8; 00157 anglePic->dim = 2; 00158 anglePic->n[0] = xDim; 00159 anglePic->n[1] = yDim; 00160 anglePic->data = (mitkIpInt1_t *)calloc(xDim*yDim,sizeof(mitkIpInt1_t)); 00161 #endif 00162 00163 picDoppler = timeSelector->GetOutput()->GetPic(); 00164 picDoppler = mitkIpFuncGausF( picDoppler, 5,2 , mitkIpFuncBorderOld ) ; 00165 00166 00167 00168 for (z=0 ; z<zDim ; z++) { 00169 00170 00171 for (y=1; y<yDim; y++) { 00172 00173 if (y<m_Origin[1]) continue; // cannot compute StrainRate above Transducer-Position 00174 00175 for (x=0; x<xDim; x++) { 00176 00177 if ((m_Origin[0] - x)==0) { // winkelhalbierende 00178 00179 int yDistanceInUnits = (int)( m_Distance/spacing[1]); 00180 int yTmp = ( (y-yDistanceInUnits)<0 ) ? 0 : (y-yDistanceInUnits); 00181 v1 = ((mitkIpUInt1_t *)picDoppler->data)[z*slice_size + yTmp*xDim + x] ; 00182 v2 = ((mitkIpUInt1_t *)picDoppler->data)[z*slice_size + y*xDim + x] ; 00183 00184 } else { 00185 00186 // compute angle to transducer position 00187 // m_Origin is given in units, not in mm 00188 alpha = atan( (float) (m_Origin[0] - x) / (float) (m_Origin[1] - y) ); 00189 00190 // m_Distance is given in mm 00191 dx = -sin(alpha) * m_Distance; 00192 dy = -cos(alpha) * m_Distance; 00193 00194 00195 // convert to units 00196 dx = dx/spacing[0]; 00197 dy = dy/spacing[1]; 00198 00199 //#define WEIGTHED 00200 #ifdef WEIGTHED 00201 weightX = dx - floor(dx); 00202 weightY = dy - floor(dy); 00203 00204 dxFloor = (int) floor(dx); 00205 dyFloor = (int) floor(dy); 00206 00207 x1 = x + dxFloor; // lower left 00208 y1 = y + dyFloor; 00209 00210 x2 = x + (dxFloor+1); // lower right 00211 y2 = y + dyFloor; 00212 00213 x3 = x + (dxFloor+1); // upper right 00214 y3 = y + (dyFloor+1); 00215 00216 x4 = x + dxFloor; // upper left 00217 y4 = y + (dyFloor+1); 00218 00219 vTmp1 = ((mitkIpUInt1_t *)picDoppler->data)[z*slice_size + y1*xDim + x1]; 00220 vTmp2 = ((mitkIpUInt1_t *)picDoppler->data)[z*slice_size + y2*xDim + x2]; 00221 vTmp3 = ((mitkIpUInt1_t *)picDoppler->data)[z*slice_size + y3*xDim + x3]; 00222 vTmp4 = ((mitkIpUInt1_t *)picDoppler->data)[z*slice_size + y4*xDim + x4]; 00223 00224 v1 = (int) ((1-weightX)*(1-weightY)*vTmp1 + weightX*(1-weightY)*vTmp2 + weightX*weightY*vTmp3 + (1-weightX)*weightY*vTmp4); 00225 v2 = ((mitkIpUInt1_t *)picDoppler->data)[z*slice_size + y*xDim + x] ; 00226 #else 00227 x1 = (int)(x + dx); 00228 y1 = (int)(y + dy); 00229 if (y1<0) y1=0; 00230 00231 #endif 00232 v1 = ((mitkIpUInt1_t *)picDoppler->data)[z*slice_size + y1*xDim + x1]; 00233 v2 = ((mitkIpUInt1_t *)picDoppler->data)[z*slice_size + y*xDim + x] ; 00234 } 00235 00236 if (v1==128) { // schaue in Gegenrichtung, falls keine SRI hier berechenbar 00237 x1 = (int)(x - dx); 00238 y1 = (int)(y - dy); 00239 if (y1>yDim) y1=yDim; 00240 //swap v1 and v2, otherwise StrainRate is calculate in false direction 00241 v1 = v2; 00242 v2 = ((mitkIpUInt1_t *)picDoppler->data)[z*slice_size + y1*xDim + x1]; 00243 } 00244 00245 if ( (v1==0 ) || (v2==0) || // wenn keine Geschwindigkeit vorhanden 00246 // oder wenn nur ganz kleine Geschwindigkeit vorhanden 00247 (v1>=(128-m_NoStrainInterval) && v1<=(128+m_NoStrainInterval)) || 00248 (v2>=(128-m_NoStrainInterval) && v2<=(128+m_NoStrainInterval))) { 00249 00250 strainRate = 128; // this means no deformation 00251 // this is neccessay due to the Cyl2Cart filter 00252 00253 } else { 00254 strainRate = (int)( (v1 - v2)/2 + 128 ); 00255 if (strainRate==128) strainRate=129; 00256 } 00257 00258 if (strainRate < minStrainRate && strainRate > 0 ) 00259 minStrainRate = strainRate; 00260 00261 if (strainRate > maxStrainRate) 00262 maxStrainRate = strainRate; 00263 00264 if (strainRate<0 ) { 00265 //strainRate = -strainRate; 00266 MITK_INFO << " error: neg. strainRate ... exit() " << std::endl 00267 << " x=" << x << " y=" << y << " z=" << z << std::endl 00268 << " v1=" << v1 << " v2=" << v2 << " dist=" << m_Distance << std::endl 00269 << " sr=" << strainRate << std::endl; 00270 exit(0); 00271 } 00272 00273 00274 ((mitkIpUInt1_t *)picStrainRate->data)[t*vol_size + z*slice_size + y*xDim + x]=strainRate; 00275 00276 // cout << "z: " << z << " y: " << y << " x: " << x << " strainrate: " << strainRate << endl; 00277 00278 #ifdef WRITE_ANGLE_PIC 00279 // if (!isAnglePicWritten) 00280 // ((mitkIpInt1_t *)anglePic->data)[y*xDim + x] = (int) ( (alpha/1.6)*128); 00281 if (!isAnglePicWritten) 00282 ((mitkIpInt1_t *)anglePic->data)[y*xDim + x] = (int) ( dx); 00283 00284 #endif 00285 00286 } // x 00287 00288 } // y 00289 00290 //isAnglePicWritten = mitkIpTrue; 00291 00292 } // z 00293 00294 //isStrainComputed = mitkIpTrue; 00295 std::string filenameD; 00296 filenameD ="doppler.pic"; 00297 mitkIpPicPut(const_cast<char *>(filenameD.c_str()),picDoppler); 00298 00299 #define WRITE_STRAIN_PIC 00300 #ifdef WRITE_STRAIN_PIC 00301 char tmpfilename[100]; 00302 sprintf(tmpfilename,"strain%d.pic",t);; 00303 mitkIpPicPut(tmpfilename,picStrainRate); 00304 #endif 00305 00306 00307 #ifdef WRITE_ANGLE_PIC 00308 std::string filename2; 00309 filename2="angle.pic"; 00310 mitkIpPicPut(const_cast<char *>(filename2.c_str()),anglePic); 00311 #endif 00312 ((mitkIpUInt1_t *)picStrainRate->data)[0]=0; 00313 ((mitkIpUInt1_t *)picStrainRate->data)[1]=255; 00314 00315 output->SetPicVolume(picStrainRate, t, n); 00316 } 00317 } 00318 00319 mitkIpPicFree(picStrainRate); 00320 00321 #define WRITE_STRAIN_PIC 00322 #ifdef WRITE_STRAIN_PIC 00323 picStrainRate = output->GetPic(); 00324 std::string filename; 00325 filename ="strain.pic"; 00326 mitkIpPicPut(const_cast<char *>(filename.c_str()),picStrainRate); 00327 #endif 00328 00329 MITK_INFO << "Strain Rate Image computed.... " << std::endl 00330 << " minStrainRate: " << minStrainRate << std::endl 00331 << " maxStrainRate: " << maxStrainRate << std::endl; 00332 } 00333 00334 mitk::DopplerToStrainRateFilter::DopplerToStrainRateFilter() 00335 : m_Distance(10), m_NoStrainInterval(2) 00336 { 00337 00338 m_Origin[0] = 0; 00339 m_Origin[1] = 0; 00340 m_Origin[2] = 0; 00341 00342 } 00343 00344 float mitk::DopplerToStrainRateFilter::GetLimit() 00345 { return (128/m_Distance); 00346 } 00347 00348 00349 mitk::DopplerToStrainRateFilter::~DopplerToStrainRateFilter() 00350 { 00351 00352 } 00353 00354 void mitk::DopplerToStrainRateFilter::GenerateInputRequestedRegion() 00355 { 00356 Superclass::GenerateInputRequestedRegion(); 00357 00358 mitk::ImageToImageFilter::InputImagePointer input = 00359 const_cast< mitk::ImageToImageFilter::InputImageType * > ( this->GetInput() ); 00360 mitk::Image::Pointer output = this->GetOutput(); 00361 00362 Image::RegionType requestedRegion; 00363 requestedRegion = output->GetRequestedRegion(); 00364 requestedRegion.SetIndex(0, 0); 00365 requestedRegion.SetIndex(1, 0); 00366 requestedRegion.SetIndex(2, 0); 00367 //requestedRegion.SetIndex(3, 0); 00368 //requestedRegion.SetIndex(4, 0); 00369 requestedRegion.SetSize(0, input->GetDimension(0)); 00370 requestedRegion.SetSize(1, input->GetDimension(1)); 00371 requestedRegion.SetSize(2, input->GetDimension(2)); 00372 //requestedRegion.SetSize(3, output->GetDimension(3)); 00373 //requestedRegion.SetSize(4, output->GetNumberOfChannels()); 00374 00375 input->SetRequestedRegion( & requestedRegion ); 00376 }
