Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
mitkThreadedToFRawDataReconstruction.cpp
Go to the documentation of this file.
1 /*===================================================================
2 
3 The Medical Imaging Interaction Toolkit (MITK)
4 
5 Copyright (c) German Cancer Research Center,
6 Division of Medical and Biological Informatics.
7 All rights reserved.
8 
9 This software is distributed WITHOUT ANY WARRANTY; without
10 even the implied warranty of MERCHANTABILITY or FITNESS FOR
11 A PARTICULAR PURPOSE.
12 
13 See LICENSE.txt or http://www.mitk.org for details.
14 
15 ===================================================================*/
16 // mitk includes
18 #include "mitkITKImageImport.h"
19 #include "mitkImageDataItem.h"
20 
21 // stl includes
22 #include <iostream>
23 #include <string>
24 #include <algorithm>
25 
26 // vtk includes
27 #include <vtkLookupTable.h>
28 #include <vtkFieldData.h>
29 #include <vtkSmartPointer.h>
30 
31 // itk includes
32 #include <itkMultiThreader.h>
33 #include <itksys/SystemTools.hxx>
34 
35 #ifdef WIN32
36 #include <winsock2.h>
37 #else
38 #include <arpa/inet.h>
39 #endif
40 #include <new>
41 
42 #define _USE_MATH_DEFINES
43 #include <math.h>
44 
45 #define PI M_PI;
46 #define cAir 299704944;
47 #define fMod 20000000;
48 
49 namespace mitk
50 {
51 
53  m_CISDist(0), m_CISAmpl(0), m_CISInten(0),
54  m_ThreadedCISDist(0), m_ThreadedCISAmpl(0), m_ThreadedCISInten(0),
55  m_Init(0), m_Width(0), m_Height(0), m_SourceDataSize(0), m_ImageSize(0), m_SourceData(0)
56  {
59 
62 
63  m_StackSize = 1;
64  }
65 
67  {
68  if(m_ThreadData != NULL)
69  delete m_ThreadData;
70 
71  if(m_CISDist != NULL)
72  delete[] m_CISDist;
73 
74  if(m_CISAmpl != NULL)
75  delete[] m_CISAmpl;
76 
77  if(m_CISInten != NULL)
78  delete[] m_CISInten;
79 
80  if(m_ThreadedCISInten != NULL)
81  delete[] m_ThreadedCISInten;
82 
83  if(m_ThreadedCISAmpl != NULL)
84  delete[] m_ThreadedCISAmpl;
85 
86  if(m_ThreadedCISDist != NULL)
87  delete[] m_ThreadedCISDist;
88 
89  }
90 
91 void ThreadedToFRawDataReconstruction::Initialize(int width, int height, int modulationFrequency,
92  int sourceDataSize )
93 {
94  m_Width = width;
95  m_Height = height;
96  m_SourceDataSize = sourceDataSize;
97  m_ImageSize = width * height;
98  m_ThreadData->m_ModulationFrequency = modulationFrequency * 1e6;
99 
100  if(!m_Init)
101  {
103  m_SourceData->SetNumberOfComponents(m_SourceDataSize);
104  m_SourceData->SetNumberOfTuples(4);
105  m_SourceData->Allocate(1);
106 
107  m_CISDist = new float[m_ImageSize];
108  m_CISAmpl = new float[m_ImageSize];
109  m_CISInten = new float[m_ImageSize];
110  m_ThreadedCISDist = new float[m_ImageSize];
111  m_ThreadedCISAmpl = new float[m_ImageSize];
112  m_ThreadedCISInten = new float[m_ImageSize];
113 
117 
118  m_Init = true;
119  }
120 }
121 
123 {
124  m_SourceData->DeepCopy(sourceData);
125 }
126 
128 {
129  memcpy(dist, m_CISDist, m_ImageSize*sizeof(float) );
130 }
131 
133 {
134  memcpy(ampl, m_CISAmpl, m_ImageSize*sizeof(float));
135 }
136 
138 {
139  memcpy(inten, m_CISInten, m_ImageSize*sizeof(float));
140 }
141 
142 void ThreadedToFRawDataReconstruction::GetAllData(float* dist, float* ampl, float* inten)
143 {
144  memcpy(dist, m_CISDist, m_ImageSize*sizeof(float) );
145  memcpy(ampl, m_CISAmpl, m_ImageSize*sizeof(float));
146  memcpy(inten, m_CISInten, m_ImageSize*sizeof(float));
147 }
148 
150  {
151  if(m_Init)
152  {
154  }
155  }
156 
158  {
159  int sourceDataSize = m_SourceDataSize;
160  int lineWidth = m_Width;
161  int frameHeight = m_Height;
162  int channelSize = lineWidth*frameHeight << 1;
163  int quadChannelSize = channelSize * 0.25;
164 
165  std::vector<short> quad = std::vector<short>(quadChannelSize);
166 
167  // clean the thread data array
169 
170  int channelNo = 0;
171  while(channelNo < m_SourceData->GetNumberOfTuples())
172  {
173  short* sourceData = new short[channelSize];
174  m_SourceData->GetTupleValue(channelNo, sourceData);
175  quad.insert(quad.begin(), sourceData, sourceData+channelSize);
176  m_ThreadData->m_InputData.push_back(quad);
177  delete[]sourceData;
178  ++channelNo;
179  }
180 
181  if(m_Threader.IsNull())
182  {
183  m_Threader = this->GetMultiThreader();
184  }
185  int maxThreadNr = 0;
186 
187  if(m_Threader->GetGlobalDefaultNumberOfThreads()> 5)
188  {
189  maxThreadNr = 5;
190  }
191  else if(m_Threader->GetGlobalMaximumNumberOfThreads()>5)
192  {
193  maxThreadNr = 5;
194  }
195  else
196  {
197  maxThreadNr = m_Threader->GetGlobalMaximumNumberOfThreads();
198  }
199 
200  if ( m_ThreadData->m_Barrier.IsNull())
201  {
203  m_ThreadData->m_Barrier->Initialize(maxThreadNr);
204  }
205  m_ThreadData->m_DataSize = quadChannelSize;
206  m_ThreadData->m_LineWidth = lineWidth;
207  m_ThreadData->m_FrameHeight = frameHeight * 0.25;
208  std::vector<int> threadIDVector;
209  int threadcounter = 0;
210  while(threadcounter != maxThreadNr-1)
211  {
212  if (m_Threader->GetNumberOfThreads() < m_Threader->GetGlobalMaximumNumberOfThreads())
213  {
214  int threadID = m_Threader->SpawnThread(this->ThreadedGenerateDataCallbackFunction, m_ThreadData);
215  threadIDVector.push_back(threadID);
216  threadcounter++;
217  }
218  }
219  m_ThreadData->m_Barrier->Wait();
220 
221  int count = 0;
222  while(count != threadIDVector.size())
223  {
224  m_Threader->TerminateThread(threadIDVector.at(count));
225  count++;
226  }
228  memcpy(m_CISDist, m_ThreadData->m_OutputData.at(0), (channelSize * 0.5)*sizeof(float));
229  memcpy(m_CISAmpl, m_ThreadData->m_OutputData.at(1), (channelSize * 0.5)*sizeof(float));
230  memcpy(m_CISInten, m_ThreadData->m_OutputData.at(2), (channelSize * 0.5)*sizeof(float));
231  m_ThreadData->m_ImageDataMutex->Unlock();
232 
233  }
234 
236  {
237  /* extract this pointer from Thread Info structure */
238  struct itk::MultiThreader::ThreadInfoStruct * pInfo = (struct itk::MultiThreader::ThreadInfoStruct*)data;
239  if (pInfo == NULL)
240  {
241  return ITK_THREAD_RETURN_VALUE;
242  }
243  if (pInfo->UserData == NULL)
244  {
245  return ITK_THREAD_RETURN_VALUE;
246  }
247  int quadrant = pInfo->ThreadID;
248  ThreadDataStruct* threadData = (ThreadDataStruct*) pInfo->UserData;
249 
250  // some needed variables
251  int x = 0;
252  double phi = 0;
253  double phi2 = 0;
254  double A1 = 0;
255  double A2 = 0;
256  double A3 = 0;
257  double A4 = 0;
258  double A5 = 0;
259  double A6 = 0;
260  double A7 = 0;
261  double A8 = 0;
262  double A3m1 = 0;
263  double A4m2 = 0;
264  double A7m5 = 0;
265  double A8m6 = 0;
266  double cair = cAir;
267  double pi = PI;
268  double twoPi = pi + pi;
269  long modFreq = fMod;
270 
271  threadData->m_ThreadDataMutex->Lock();
272  std::vector<short> quad1 = threadData->m_InputData.at(0);
273  std::vector<short> quad2 = threadData->m_InputData.at(1);
274  std::vector<short> quad3 = threadData->m_InputData.at(2);
275  std::vector<short> quad4 = threadData->m_InputData.at(3);
276  int index = quadrant << 1;
277  int index2 = 3-quadrant;
278  modFreq = threadData->m_ModulationFrequency;
279  int linewidth = threadData->m_LineWidth;
280  int frameheight = threadData->m_FrameHeight;
281  threadData->m_ThreadDataMutex->Unlock();
282 
283  double intermed1 = cair/(pi*(modFreq << 2));
284  double intermed2 = intermed1*500;
285  int doubleLwidth = linewidth << 1;
286  int datasize = doubleLwidth*frameheight << 2;
287 
288 
289  do
290  {
291  index += doubleLwidth;
292  x++;
293  do
294  {
295  index -= 8;
296  A1 = htons(quad1.at(index));
297  A2 = htons(quad2.at(index));
298  A3 = htons(quad3.at(index));
299  A4 = htons(quad4.at(index));
300  A5 = htons(quad1.at(index+1));
301  A6 = htons(quad2.at(index+1));
302  A7 = htons(quad3.at(index+1));
303  A8 = htons(quad4.at(index+1));
304 
305  phi = atan2((A3 - A1),(A2 - A4)) + pi;
306  phi2 = atan2((A7 - A5),(A6 - A8));
307  if(phi2<0) phi2 +=twoPi;
308 
309  A3m1 = A3*A3 - 2*A3*A1 + A1*A1;
310  A4m2 = A4*A4 - 2*A4*A2 + A2*A2;
311  A7m5 = A7*A7 - 2*A7*A5 + A5*A5;
312  A8m6 = A8*A8 - 2*A8*A6 + A6*A6;
313  threadData->m_ImageDataMutex->Lock();
314  threadData->m_OutputData.at(0)[index2] = (phi+phi2)*intermed2; //(((phi*intermed1) + (phi2*intermed1))/2)*1000;
315  threadData->m_OutputData.at(1)[index2] = (sqrt(A3m1 + A4m2)+sqrt(A7m5 + A8m6))*0.5; //(sqrt(A3m1 + A4m2)/2) + (sqrt(A7m5 + A8m6)/2);
316  threadData->m_OutputData.at(2)[index2] = (A1+A2+A3+A4+A5+A6+A7+A8)*0.125;
317  threadData->m_ImageDataMutex->Unlock();
318 
319  index2 += 4;
320  }while(index2 <= (x*linewidth) - (1+quadrant));
321 
322  index += doubleLwidth;
323 
324  }while(index < datasize);
325 
326  threadData->m_Barrier->Wait();
327 
328  return ITK_THREAD_RETURN_VALUE;
329  }
330 
332  {
333  this->GenerateData();
334  }
335 
336 } // end mitk namespace
void Initialize(int width, int height, int modulationFrequency, int sourceDataSize)
itk::FastMutexLock::Pointer m_ImageDataMutex
mutex for coordinated access to image data
virtual void GenerateData()
method generating the outputs of this filter. Called in the updated process of the pipeline...
float * m_CISDist
holds the distance information from for one distance image slice
DataCollection - Class to facilitate loading/accessing structured data.
float * m_CISAmpl
holds the amplitude information from for one amplitude image slice
itk::Barrier::Pointer m_Barrier
barrier to synchronize ends of threads
std::vector< std::vector< short > > m_InputData
static ITK_THREAD_RETURN_TYPE ThreadedGenerateDataCallbackFunction(void *data)
threader callback function for multi threaded data generation
float * m_CISInten
holds the intensity information from for one intensity image slice
virtual void BeforeThreadedGenerateData()
method configures the camera output and prepares the thread data struct for threaded data generation ...
itk::FastMutexLock::Pointer m_ThreadDataMutex
mutex to control access to images
void GetAllData(float *dist, float *ampl, float *inten)
static itkEventMacro(BoundingShapeInteractionEvent, itk::AnyEvent) class MITKBOUNDINGSHAPE_EXPORT BoundingShapeInteractor Pointer New()
Basic interaction methods for mitk::GeometryData.