mitkThreadedToFRawDataReconstruction.cpp

Go to the documentation of this file.

/*============================================================================

The Medical Imaging Interaction Toolkit (MITK)

Copyright (c) German Cancer Research Center (DKFZ)
All rights reserved.

Use of this source code is governed by a 3-clause BSD license that can be
found in the LICENSE file.

============================================================================*/
// mitk includes
#include "mitkThreadedToFRawDataReconstruction.h"
#include "mitkITKImageImport.h"
#include "mitkImageDataItem.h"

// stl includes
#include <iostream>
#include <string>
#include <algorithm>

// vtk includes
#include <vtkLookupTable.h>
#include <vtkFieldData.h>
#include <vtkSmartPointer.h>

// itk includes
#include <itkMultiThreader.h>
#include <itksys/SystemTools.hxx>

#ifdef WIN32
#include <winsock2.h>
#else
#include <arpa/inet.h>
#endif
#include <new>

#define cAir 299704944;
#define fMod 20000000;

namespace mitk
{

  ThreadedToFRawDataReconstruction::ThreadedToFRawDataReconstruction(): m_Threader(0),
    m_CISDist(0), m_CISAmpl(0), m_CISInten(0),
    m_ThreadedCISDist(0), m_ThreadedCISAmpl(0), m_ThreadedCISInten(0),
    m_Init(0), m_Width(0), m_Height(0), m_SourceDataSize(0), m_ImageSize(0), m_SourceData(0)
  {
    m_ThreadData = new ThreadDataStruct;
    m_ThreadData->m_ModulationFrequency = fMod;

    m_ThreadData->m_ImageDataMutex = itk::FastMutexLock::New();
    m_ThreadData->m_ThreadDataMutex = itk::FastMutexLock::New();

    m_StackSize = 1;
  }

  ThreadedToFRawDataReconstruction::~ThreadedToFRawDataReconstruction()
  {
    if(m_ThreadData != nullptr)
      delete m_ThreadData;

    if(m_CISDist != nullptr)
      delete[] m_CISDist;

    if(m_CISAmpl != nullptr)
      delete[] m_CISAmpl;

    if(m_CISInten != nullptr)
      delete[] m_CISInten;

    if(m_ThreadedCISInten != nullptr)
      delete[] m_ThreadedCISInten;

    if(m_ThreadedCISAmpl != nullptr)
      delete[] m_ThreadedCISAmpl;

    if(m_ThreadedCISDist != nullptr)
      delete[] m_ThreadedCISDist;

  }

void ThreadedToFRawDataReconstruction::Initialize(int width, int height, int modulationFrequency,
                                                  int sourceDataSize )
{
  m_Width = width;
  m_Height = height;
  m_SourceDataSize = sourceDataSize;
  m_ImageSize = width * height;
  m_ThreadData->m_ModulationFrequency = modulationFrequency * 1e6;

  if(!m_Init)
  {
    m_SourceData = vtkShortArray::New();
    m_SourceData->SetNumberOfComponents(m_SourceDataSize);
    m_SourceData->SetNumberOfTuples(4);
    m_SourceData->Allocate(1);

    m_CISDist = new float[m_ImageSize];
    m_CISAmpl = new float[m_ImageSize];
    m_CISInten = new float[m_ImageSize];
    m_ThreadedCISDist = new float[m_ImageSize];
    m_ThreadedCISAmpl = new float[m_ImageSize];
    m_ThreadedCISInten = new float[m_ImageSize];

    m_ThreadData->m_OutputData.push_back( m_ThreadedCISDist );
    m_ThreadData->m_OutputData.push_back( m_ThreadedCISAmpl );
    m_ThreadData->m_OutputData.push_back( m_ThreadedCISInten );

    m_Init = true;
  }
}

void ThreadedToFRawDataReconstruction::SetChannelData(vtkShortArray* sourceData)
{
  m_SourceData->DeepCopy(sourceData);
}

void ThreadedToFRawDataReconstruction::GetDistances(float* dist)
{
  memcpy(dist, m_CISDist, m_ImageSize*sizeof(float) );
}

void ThreadedToFRawDataReconstruction::GetAmplitudes(float* ampl)
{
  memcpy(ampl, m_CISAmpl, m_ImageSize*sizeof(float));
}

void ThreadedToFRawDataReconstruction::GetIntensities(float* inten)
{
  memcpy(inten, m_CISInten, m_ImageSize*sizeof(float));
}

void ThreadedToFRawDataReconstruction::GetAllData(float* dist, float* ampl, float* inten)
{
  memcpy(dist, m_CISDist, m_ImageSize*sizeof(float) );
  memcpy(ampl, m_CISAmpl, m_ImageSize*sizeof(float));
  memcpy(inten, m_CISInten, m_ImageSize*sizeof(float));
}

  void ThreadedToFRawDataReconstruction::GenerateData()
  {
    if(m_Init)
    {
      this->BeforeThreadedGenerateData();
    }
  }

  void ThreadedToFRawDataReconstruction::BeforeThreadedGenerateData()
  {
    int sourceDataSize = m_SourceDataSize;
    int lineWidth = m_Width;
    int frameHeight = m_Height;
    int channelSize = lineWidth*frameHeight << 1;
    int quadChannelSize = channelSize * 0.25;

    std::vector<short> quad = std::vector<short>(quadChannelSize);

    // clean the thread data array
    m_ThreadData->m_InputData.erase(m_ThreadData->m_InputData.begin(),m_ThreadData->m_InputData.end());

    int channelNo = 0;
    while(channelNo < m_SourceData->GetNumberOfTuples())
    {
      short* sourceData = new short[channelSize];
      m_SourceData->GetTupleValue(channelNo, sourceData);
      quad.insert(quad.begin(), sourceData,  sourceData+channelSize);
      m_ThreadData->m_InputData.push_back(quad);
      delete[]sourceData;
      ++channelNo;
    }

    if(m_Threader.IsNull())
    {
      m_Threader = this->GetMultiThreader();
    }
    int maxThreadNr = 0;

    if(m_Threader->GetGlobalDefaultNumberOfThreads()> 5)
    {
      maxThreadNr = 5;
    }
    else if(m_Threader->GetGlobalMaximumNumberOfThreads()>5)
    {
      maxThreadNr = 5;
    }
    else
    {
      maxThreadNr = m_Threader->GetGlobalMaximumNumberOfThreads();
    }

    if ( m_ThreadData->m_Barrier.IsNull())
    {
      m_ThreadData->m_Barrier = itk::Barrier::New();
      m_ThreadData->m_Barrier->Initialize(maxThreadNr);
    }
    m_ThreadData->m_DataSize = quadChannelSize;
    m_ThreadData->m_LineWidth = lineWidth;
    m_ThreadData->m_FrameHeight = frameHeight * 0.25;
    std::vector<int> threadIDVector;
    int threadcounter = 0;
    while(threadcounter != maxThreadNr-1)
    {
      if (m_Threader->GetNumberOfThreads() < m_Threader->GetGlobalMaximumNumberOfThreads())
      {
        int threadID = m_Threader->SpawnThread(this->ThreadedGenerateDataCallbackFunction, m_ThreadData);
        threadIDVector.push_back(threadID);
        threadcounter++;
      }
    }
    m_ThreadData->m_Barrier->Wait();

    int count = 0;
    while(count != threadIDVector.size())
    {
      m_Threader->TerminateThread(threadIDVector.at(count));
      count++;
    }
    m_ThreadData->m_ImageDataMutex->Lock();
    memcpy(m_CISDist, m_ThreadData->m_OutputData.at(0), (channelSize * 0.5)*sizeof(float));
    memcpy(m_CISAmpl, m_ThreadData->m_OutputData.at(1), (channelSize * 0.5)*sizeof(float));
    memcpy(m_CISInten, m_ThreadData->m_OutputData.at(2), (channelSize * 0.5)*sizeof(float));
    m_ThreadData->m_ImageDataMutex->Unlock();

  }

  ITK_THREAD_RETURN_TYPE  ThreadedToFRawDataReconstruction::ThreadedGenerateDataCallbackFunction(void* data)
  {
    /* extract this pointer from Thread Info structure */
    struct itk::MultiThreader::ThreadInfoStruct * pInfo = (struct itk::MultiThreader::ThreadInfoStruct*)data;
    if (pInfo == nullptr)
    {
      return ITK_THREAD_RETURN_VALUE;
    }
    if (pInfo->UserData == nullptr)
    {
      return ITK_THREAD_RETURN_VALUE;
    }
    int quadrant = pInfo->ThreadID;
    ThreadDataStruct* threadData = (ThreadDataStruct*) pInfo->UserData;

    // some needed variables
    int x = 0;
    double phi = 0;
    double phi2 = 0;
    double A1 = 0;
    double A2 = 0;
    double A3 = 0;
    double A4 = 0;
    double A5 = 0;
    double A6 = 0;
    double A7 = 0;
    double A8 = 0;
    double A3m1 = 0;
    double A4m2 = 0;
    double A7m5 = 0;
    double A8m6 = 0;
    double cair = cAir;
    double twoPi = itk::Math::pi + itk::Math::pi;
    long modFreq = fMod;

    threadData->m_ThreadDataMutex->Lock();
    std::vector<short>  quad1 = threadData->m_InputData.at(0);
    std::vector<short>  quad2 = threadData->m_InputData.at(1);
    std::vector<short>  quad3 = threadData->m_InputData.at(2);
    std::vector<short>  quad4 = threadData->m_InputData.at(3);
    int index = quadrant << 1;
    int index2 = 3-quadrant;
    modFreq = threadData->m_ModulationFrequency;
    int linewidth = threadData->m_LineWidth;
    int frameheight = threadData->m_FrameHeight;
    threadData->m_ThreadDataMutex->Unlock();

    double intermed1 = cair/(itk::Math::pi*(modFreq << 2));
    double intermed2 = intermed1*500;
    int doubleLwidth = linewidth << 1;
    int datasize = doubleLwidth*frameheight << 2;


    do
    {
      index += doubleLwidth;
      x++;
      do
      {
        index -= 8;
        A1 = htons(quad1.at(index));
        A2 = htons(quad2.at(index));
        A3 = htons(quad3.at(index));
        A4 = htons(quad4.at(index));
        A5 = htons(quad1.at(index+1));
        A6 = htons(quad2.at(index+1));
        A7 = htons(quad3.at(index+1));
        A8 = htons(quad4.at(index+1));

        phi  = atan2((A3 - A1),(A2 - A4)) + itk::Math::pi;
        phi2 = atan2((A7 - A5),(A6 - A8));
        if(phi2<0) phi2 +=twoPi;

        A3m1 = A3*A3 - 2*A3*A1 + A1*A1;
        A4m2 = A4*A4 - 2*A4*A2 + A2*A2;
        A7m5 = A7*A7 - 2*A7*A5 + A5*A5;
        A8m6 = A8*A8 - 2*A8*A6 + A6*A6;
        threadData->m_ImageDataMutex->Lock();
        threadData->m_OutputData.at(0)[index2] = (phi+phi2)*intermed2;                        //(((phi*intermed1) + (phi2*intermed1))/2)*1000;
        threadData->m_OutputData.at(1)[index2] = (sqrt(A3m1 + A4m2)+sqrt(A7m5 + A8m6))*0.5;   //(sqrt(A3m1 + A4m2)/2) + (sqrt(A7m5 + A8m6)/2);
        threadData->m_OutputData.at(2)[index2] = (A1+A2+A3+A4+A5+A6+A7+A8)*0.125;
        threadData->m_ImageDataMutex->Unlock();

        index2 += 4;
      }while(index2 <= (x*linewidth) - (1+quadrant));

      index += doubleLwidth;

    }while(index < datasize);

    threadData->m_Barrier->Wait();

    return ITK_THREAD_RETURN_VALUE;
  }

  void ThreadedToFRawDataReconstruction::Update()
  {
    this->GenerateData();
  }

} // end mitk namespace