diff --git a/Modules/OpenCL/mitkOclDataSet.cpp b/Modules/OpenCL/mitkOclDataSet.cpp index 269a2d01d3..953c0d68b1 100644 --- a/Modules/OpenCL/mitkOclDataSet.cpp +++ b/Modules/OpenCL/mitkOclDataSet.cpp @@ -1,174 +1,176 @@ /*=================================================================== 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 "mitkOclDataSet.h" #include "mitkCommon.h" #include "mitkLogMacros.h" #include "mitkOclUtils.h" #include +//#define SHOW_MEM_INFO + mitk::OclDataSet::OclDataSet() : m_gpuBuffer(nullptr), m_context(nullptr), m_bufferSize(0), m_gpuModified(false), m_cpuModified(false), m_Data(nullptr), m_BpE(1) { } mitk::OclDataSet::~OclDataSet() { MITK_DEBUG << "OclDataSet Destructor"; //release GMEM Image buffer if (m_gpuBuffer) clReleaseMemObject(m_gpuBuffer); } cl_mem mitk::OclDataSet::CreateGPUBuffer() { MITK_DEBUG << "InitializeGPUBuffer call with: BPE=" << m_BpE; us::ServiceReference ref = GetModuleContext()->GetServiceReference(); OclResourceService* resources = GetModuleContext()->GetService(ref); m_context = resources->GetContext(); int clErr; if (m_gpuBuffer) clReleaseMemObject(m_gpuBuffer); m_gpuBuffer = clCreateBuffer(m_context, CL_MEM_READ_WRITE, m_bufferSize * (size_t)m_BpE, nullptr, &clErr); - MITK_DEBUG << "Created GPU Buffer Object of size: " << (size_t)m_BpE * m_bufferSize << " Bytes"; + #ifdef SHOW_MEM_INFO + MITK_INFO << "Created GPU Buffer Object of size: " << (size_t)m_BpE * m_bufferSize << " Bytes"; + #endif CHECK_OCL_ERR(clErr); if (clErr != CL_SUCCESS) mitkThrow() << "openCL Error when creating Buffer"; return m_gpuBuffer; } bool mitk::OclDataSet::IsModified(int _type) { if (_type) return m_cpuModified; else return m_gpuModified; } void mitk::OclDataSet::Modified(int _type) { // defines... GPU: 0, CPU: 1 m_cpuModified = _type; m_gpuModified = !_type; } int mitk::OclDataSet::TransferDataToGPU(cl_command_queue gpuComQueue) { cl_int clErr = 0; // check whether an image present if (m_Data == nullptr){ MITK_ERROR("ocl.DataSet") << "(mitk) No data present!\n"; return -1; } // there is a need for copy only if RAM-Data newer then GMEM data if (m_cpuModified) { //check the buffer if(m_gpuBuffer == nullptr) { CreateGPUBuffer(); } if (m_gpuBuffer != nullptr) { clErr = clEnqueueWriteBuffer(gpuComQueue, m_gpuBuffer, CL_TRUE, 0, m_bufferSize * (size_t)m_BpE, m_Data, 0, NULL, NULL); MITK_DEBUG << "Wrote Data to GPU Buffer Object."; CHECK_OCL_ERR(clErr); m_gpuModified = true; if (clErr != CL_SUCCESS) mitkThrow() << "openCL Error when writing Buffer"; } else { MITK_ERROR << "No GPU buffer present!"; } } return clErr; } cl_mem mitk::OclDataSet::GetGPUBuffer() { - // clGetMemObjectInfo() - cl_mem_object_type memInfo; - // query image object info only if already initialized if( this->m_gpuBuffer ) { - #ifdef MBILOG_ENABLE_DEBUG + #ifdef SHOW_MEM_INFO cl_int clErr = 0; + // clGetMemObjectInfo() + cl_mem_object_type memInfo; clErr = clGetMemObjectInfo(this->m_gpuBuffer, CL_MEM_TYPE, sizeof(cl_mem_object_type), &memInfo, nullptr ); CHECK_OCL_ERR(clErr); + MITK_DEBUG << "Querying info for object, recieving: " << memInfo; #endif } - MITK_DEBUG << "Querying info for object, recieving: " << memInfo; - return m_gpuBuffer; } void* mitk::OclDataSet::TransferDataToCPU(cl_command_queue gpuComQueue) { cl_int clErr = 0; // if image created on GPU, needs to create mitk::Image if( m_gpuBuffer == nullptr ){ MITK_ERROR("ocl.DataSet") << "(mitk) No buffer present!\n"; return nullptr; } // check buffersize char* data = new char[m_bufferSize * (size_t)m_BpE]; // debug info - #ifdef MBILOG_ENABLE_DEBUG + #ifdef SHOW_MEM_INFO oclPrintMemObjectInfo( m_gpuBuffer ); #endif clErr = clEnqueueReadBuffer( gpuComQueue, m_gpuBuffer, CL_TRUE, 0, m_bufferSize * (size_t)m_BpE, data ,0, nullptr, nullptr); CHECK_OCL_ERR(clErr); if(clErr != CL_SUCCESS) mitkThrow() << "openCL Error when reading Output Buffer"; clFlush( gpuComQueue ); // the cpu data is same as gpu this->m_gpuModified = false; return (void*) data; } void mitk::OclDataSet::SetBufferSize(size_t size) { m_bufferSize = size; } void mitk::OclDataSet::SetBpE(unsigned short BpE) { m_BpE = BpE; } diff --git a/Modules/OpenCL/mitkOclDataSetToDataSetFilter.cpp b/Modules/OpenCL/mitkOclDataSetToDataSetFilter.cpp index 852f487a40..2cf7acab5d 100644 --- a/Modules/OpenCL/mitkOclDataSetToDataSetFilter.cpp +++ b/Modules/OpenCL/mitkOclDataSetToDataSetFilter.cpp @@ -1,128 +1,128 @@ /*=================================================================== 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 "mitkOclDataSetToDataSetFilter.h" #include "mitkOclDataSet.h" #include "mitkException.h" mitk::OclDataSetToDataSetFilter::OclDataSetToDataSetFilter() { m_Output = mitk::OclDataSet::New(); } mitk::OclDataSetToDataSetFilter::~OclDataSetToDataSetFilter() { } mitk::OclDataSet::Pointer mitk::OclDataSetToDataSetFilter::GetGPUOutput() { return this->m_Output; } void* mitk::OclDataSetToDataSetFilter::GetOutput() { void* pData = m_Output->TransferDataToCPU(m_CommandQue); return pData; } int mitk::OclDataSetToDataSetFilter::GetBytesPerElem() { return this->m_CurrentSizeOutput; } bool mitk::OclDataSetToDataSetFilter::InitExec(cl_kernel ckKernel, unsigned int* dimensions, size_t outputDataSize, unsigned int outputBpE) { cl_int clErr = 0; if (m_Input.IsNull()) mitkThrow() << "Input DataSet is null."; // get DataSet size once const unsigned int uiDataSetWidth = dimensions[0]; const unsigned int uiDataSetHeight = dimensions[1]; const unsigned int uiDataSetDepth = dimensions[2]; // compute work sizes this->SetWorkingSize(8, uiDataSetWidth, 8, uiDataSetHeight, 8, uiDataSetDepth); cl_mem clBuffIn = m_Input->GetGPUBuffer(); cl_mem clBuffOut = m_Output->GetGPUBuffer(); if (!clBuffIn) { if (m_Input->TransferDataToGPU(m_CommandQue) != CL_SUCCESS) { mitkThrow() << "Could not create / initialize gpu DataSet."; } clBuffIn = m_Input->GetGPUBuffer(); } // output DataSet not initialized or output buffer size changed if (!clBuffOut || (size_t)m_Output->GetBufferSize() != outputDataSize) { MITK_DEBUG << "Create GPU DataSet call " << uiDataSetWidth << "x" << uiDataSetHeight << "x" << uiDataSetDepth; m_Output->SetBpE(outputBpE); m_Output->SetBufferSize(outputDataSize); clBuffOut = m_Output->CreateGPUBuffer(); m_CurrentSizeOutput = outputBpE; } clErr = 0; clErr = clSetKernelArg(ckKernel, 0, sizeof(cl_mem), &clBuffIn); clErr |= clSetKernelArg(ckKernel, 1, sizeof(cl_mem), &clBuffOut); CHECK_OCL_ERR(clErr); if (clErr != CL_SUCCESS) mitkThrow() << "OpenCL Part initialization failed with " << GetOclErrorAsString(clErr); return(clErr == CL_SUCCESS); } bool mitk::OclDataSetToDataSetFilter::InitExecNoInput(cl_kernel ckKernel, unsigned int* dimensions, size_t outputDataSize, unsigned int outputBpE) { cl_int clErr = 0; // get DataSet size once const unsigned int uiDataSetWidth = dimensions[0]; const unsigned int uiDataSetHeight = dimensions[1]; const unsigned int uiDataSetDepth = dimensions[2]; // compute work sizes this->SetWorkingSize(8, uiDataSetWidth, 8, uiDataSetHeight, 8, uiDataSetDepth); cl_mem clBuffOut = m_Output->GetGPUBuffer(); - // output DataSet not initialized - if (!clBuffOut) + // output DataSet not initialized or output buffer size changed + if (!clBuffOut || (size_t)m_Output->GetBufferSize() != outputDataSize) { MITK_DEBUG << "Create GPU DataSet call " << uiDataSetWidth << "x" << uiDataSetHeight << "x" << uiDataSetDepth; m_Output->SetBpE(outputBpE); m_Output->SetBufferSize(outputDataSize); clBuffOut = m_Output->CreateGPUBuffer(); m_CurrentSizeOutput = outputBpE; } clErr = clSetKernelArg(ckKernel, 0, sizeof(cl_mem), &clBuffOut); CHECK_OCL_ERR(clErr); if (clErr != CL_SUCCESS) mitkThrow() << "OpenCL Part initialization failed with " << GetOclErrorAsString(clErr); return(clErr == CL_SUCCESS); } \ No newline at end of file diff --git a/Modules/OpenCL/mitkOclFilter.cpp b/Modules/OpenCL/mitkOclFilter.cpp index d1b1c23aab..7965a24484 100644 --- a/Modules/OpenCL/mitkOclFilter.cpp +++ b/Modules/OpenCL/mitkOclFilter.cpp @@ -1,284 +1,347 @@ /*=================================================================== 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. ===================================================================*/ //Ocl #include "mitkOclFilter.h" #include "mitkOclUtils.h" #include "mitkOpenCLActivator.h" //Mitk #include #include //usService #include "usServiceReference.h" #include #include #include #include #include #include mitk::OclFilter::OclFilter() : m_ClCompilerFlags(""), m_ClProgram(nullptr), m_CommandQue(nullptr), m_FilterID("mitkOclFilter"), m_Preambel(" "), m_Initialized(false) { } mitk::OclFilter::OclFilter(const char* filename) : m_ClCompilerFlags(""), m_ClProgram(nullptr), m_CommandQue(nullptr), m_FilterID(filename), m_Preambel(" "), m_Initialized(false) { m_ClFiles.push_back(filename); } mitk::OclFilter::~OclFilter() { MITK_DEBUG << "OclFilter Destructor"; // release program if (m_ClProgram) { us::ServiceReference ref = GetModuleContext()->GetServiceReference(); OclResourceService* resources = GetModuleContext()->GetService(ref); // remove program from storage resources->RemoveProgram(m_FilterID); } } bool mitk::OclFilter::ExecuteKernel( cl_kernel kernel, unsigned int workSizeDim ) { cl_int clErr = 0; clErr = clEnqueueNDRangeKernel( this->m_CommandQue, kernel, workSizeDim, nullptr, this->m_GlobalWorkSize, m_LocalWorkSize, 0, nullptr, nullptr); CHECK_OCL_ERR( clErr ); return ( clErr == CL_SUCCESS ); } +#include + bool mitk::OclFilter::ExecuteKernelChunks( cl_kernel kernel, unsigned int workSizeDim, size_t* chunksDim ) { size_t offset[3] ={0, 0, 0}; cl_int clErr = 0; if(workSizeDim == 2) { for(offset[0] = 0; offset[0] < m_GlobalWorkSize[0]; offset[0] += chunksDim[0]) { for(offset[1] = 0; offset[1] < m_GlobalWorkSize[1]; offset[1] += chunksDim[1]) { - clErr |= clEnqueueNDRangeKernel( this->m_CommandQue, kernel, workSizeDim, - offset, chunksDim, m_LocalWorkSize, 0, nullptr, nullptr); + clErr |= clEnqueueNDRangeKernel(this->m_CommandQue, kernel, workSizeDim, + offset, chunksDim, m_LocalWorkSize, 0, nullptr, nullptr); } } } else if(workSizeDim == 3) { for(offset[0] = 0; offset[0] < m_GlobalWorkSize[0]; offset[0] += chunksDim[0]) { for(offset[1] = 0; offset[1] < m_GlobalWorkSize[1]; offset[1] += chunksDim[1]) { for(offset[2] = 0; offset[2] < m_GlobalWorkSize[2]; offset[2] += chunksDim[2]) { clErr |= clEnqueueNDRangeKernel( this->m_CommandQue, kernel, workSizeDim, offset, chunksDim, m_LocalWorkSize, 0, nullptr, nullptr); } } } } CHECK_OCL_ERR(clErr); return ( clErr == CL_SUCCESS ); } +bool mitk::OclFilter::ExecuteKernelChunksInBatches(cl_kernel kernel, unsigned int workSizeDim, size_t* chunksDim, size_t batchSize, int waitTimems) +{ + size_t offset[3] = { 0, 0, 0 }; + cl_int clErr = 0; + + unsigned int currentChunk = 0; + cl_event* waitFor = new cl_event[batchSize]; + + if (workSizeDim == 2) + { + for (offset[0] = 0; offset[0] < m_GlobalWorkSize[0]; offset[0] += chunksDim[0]) + { + for (offset[1] = 0; offset[1] < m_GlobalWorkSize[1]; offset[1] += chunksDim[1]) + { + if (currentChunk % batchSize == 0 && currentChunk != 0) + { + clWaitForEvents(batchSize, &waitFor[0]); + std::this_thread::sleep_for(std::chrono::milliseconds(waitTimems)); + clErr |= clEnqueueNDRangeKernel(this->m_CommandQue, kernel, workSizeDim, + offset, chunksDim, m_LocalWorkSize, 0, nullptr, &waitFor[0]); + } + else + { + clErr |= clEnqueueNDRangeKernel(this->m_CommandQue, kernel, workSizeDim, + offset, chunksDim, m_LocalWorkSize, 0, nullptr, &waitFor[currentChunk % batchSize]); + } + currentChunk++; + } + } + } + else if (workSizeDim == 3) + { + for (offset[0] = 0; offset[0] < m_GlobalWorkSize[0]; offset[0] += chunksDim[0]) + { + for (offset[1] = 0; offset[1] < m_GlobalWorkSize[1]; offset[1] += chunksDim[1]) + { + for (offset[2] = 0; offset[2] < m_GlobalWorkSize[2]; offset[2] += chunksDim[2]) + { + if (currentChunk % batchSize == 0 && currentChunk != 0) + { + clWaitForEvents(batchSize, &waitFor[0]); + std::this_thread::sleep_for(std::chrono::milliseconds(waitTimems)); + clErr |= clEnqueueNDRangeKernel(this->m_CommandQue, kernel, workSizeDim, + offset, chunksDim, m_LocalWorkSize, 0, nullptr, &waitFor[0]); + } + else + { + clErr |= clEnqueueNDRangeKernel(this->m_CommandQue, kernel, workSizeDim, + offset, chunksDim, m_LocalWorkSize, 0, nullptr, &waitFor[currentChunk % batchSize]); + } + currentChunk++; + } + } + } + } + + CHECK_OCL_ERR(clErr); + + return (clErr == CL_SUCCESS); +} + bool mitk::OclFilter::Initialize() { us::ServiceReference ref = GetModuleContext()->GetServiceReference(); OclResourceService* resources = GetModuleContext()->GetService(ref); m_CommandQue = resources->GetCommandQueue(); cl_int clErr = 0; m_Initialized = CHECK_OCL_ERR(clErr); if ( m_ClFiles.empty()) { MITK_ERROR<<"No OpenCL Source FILE specified"; return false; } if (m_ClProgram == nullptr) { try { this->m_ClProgram = resources->GetProgram( this->m_FilterID ); } catch(const mitk::Exception& e) { MITK_INFO << "Program not stored in resource manager, compiling. " << e; this->CompileSource(); } } return m_Initialized; } void mitk::OclFilter::LoadSourceFiles(CStringList &sourceCode, ClSizeList &sourceCodeSize) { for( CStringList::iterator it = m_ClFiles.begin(); it != m_ClFiles.end(); ++it ) { MITK_DEBUG << "Load file :" << *it; us::ModuleResource mdr = GetModule()->GetResource(*it); if( !mdr.IsValid() ) MITK_WARN << "Could not load resource: " << mdr.GetName() << " is invalid!"; us::ModuleResourceStream rss(mdr); // read resource file to a string std::istreambuf_iterator eos; std::string source(std::istreambuf_iterator(rss), eos); // add preambel and build up string to compile std::string src(m_Preambel); src.append("\n"); src.append(source); // allocate new char buffer char* tmp = new char[src.size() + 1]; strcpy(tmp,src.c_str()); // add source to list sourceCode.push_back((const char*)tmp); sourceCodeSize.push_back(src.size()); } } void mitk::OclFilter::CompileSource() { // helper variable int clErr = 0; CStringList sourceCode; ClSizeList sourceCodeSize; if (m_ClFiles.empty()) { MITK_ERROR("ocl.filter") << "No shader source file was set"; return; } //get a valid opencl context us::ServiceReference ref = GetModuleContext()->GetServiceReference(); OclResourceService* resources = GetModuleContext()->GetService(ref); cl_context gpuContext = resources->GetContext(); // load the program source from file LoadSourceFiles(sourceCode, sourceCodeSize); if ( !sourceCode.empty() ) { // create program from all files in the file list m_ClProgram = clCreateProgramWithSource(gpuContext, sourceCode.size(), &sourceCode[0], &sourceCodeSize[0], &clErr); CHECK_OCL_ERR(clErr); // build the source code MITK_DEBUG << "Building Program Source"; std::string compilerOptions = ""; compilerOptions.append(m_ClCompilerFlags); MITK_DEBUG("ocl.filter") << "cl compiler flags: " << compilerOptions.c_str(); clErr = clBuildProgram(m_ClProgram, 0, nullptr, compilerOptions.c_str(), nullptr, nullptr); CHECK_OCL_ERR(clErr); // if OpenCL Source build failed if (clErr != CL_SUCCESS) { MITK_ERROR("ocl.filter") << "Failed to build source"; oclLogBuildInfo(m_ClProgram, resources->GetCurrentDevice() ); oclLogBinary(m_ClProgram, resources->GetCurrentDevice() ); m_Initialized = false; } // store the succesfully build program into the program storage provided by the resource service resources->InsertProgram(m_ClProgram, m_FilterID, true); // free the char buffers with the source code for( CStringList::iterator it = sourceCode.begin(); it != sourceCode.end(); ++it ) { delete[] *it; } } else { MITK_ERROR("ocl.filter") << "Could not load from source"; m_Initialized = false; } } void mitk::OclFilter::SetWorkingSize(unsigned int locx, unsigned int dimx, unsigned int locy, unsigned int dimy, unsigned int locz, unsigned int dimz) { // set the local work size this->m_LocalWorkSize[0] = locx; this->m_LocalWorkSize[1] = locy; this->m_LocalWorkSize[2] = locz; this->m_GlobalWorkSize[0] = dimx; this->m_GlobalWorkSize[1] = dimy; this->m_GlobalWorkSize[2] = dimz; // estimate the global work size this->m_GlobalWorkSize[0] = iDivUp( dimx, this->m_LocalWorkSize[0]) * this->m_LocalWorkSize[0]; if ( dimy > 1) this->m_GlobalWorkSize[1] = iDivUp( dimy, this->m_LocalWorkSize[1]) * this->m_LocalWorkSize[1]; if( dimz > 1 ) this->m_GlobalWorkSize[2] = iDivUp( dimz, this->m_LocalWorkSize[2]) * this->m_LocalWorkSize[2]; } void mitk::OclFilter::SetSourcePreambel(const char* preambel) { this->m_Preambel = preambel; } void mitk::OclFilter::AddSourceFile(const char* filename) { m_ClFiles.push_back(filename); } void mitk::OclFilter::SetCompilerFlags(const char* flags) { m_ClCompilerFlags = flags; } bool mitk::OclFilter::IsInitialized() { return m_Initialized; } diff --git a/Modules/OpenCL/mitkOclFilter.h b/Modules/OpenCL/mitkOclFilter.h index ba3705a466..41fcaa8b0f 100644 --- a/Modules/OpenCL/mitkOclFilter.h +++ b/Modules/OpenCL/mitkOclFilter.h @@ -1,153 +1,157 @@ /*=================================================================== 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 __mitkOclFilter_h #define __mitkOclFilter_h #include "mitkOclUtils.h" #include "mitkCommon.h" #include #include #include namespace mitk { /** @class OclFilter @brief Superclass for all OpenCL based filter. This class takes care of loading and compiling the external GPU program code. */ class MITKOPENCL_EXPORT OclFilter { public: /** * @brief Add a source file from the resource files to the * OpenCL shader file list. Multiple files can be added to the list. * * @param name of the file in the resource system */ void AddSourceFile(const char* filename); /** * @brief Set specific compilerflags to compile the CL source. Default is set to nullptr; * example: "-cl-fast-relaxed-math -cl-mad-enable -cl-strict-aliasing" * * @param flags to the modulefolder that contains the gpuSource */ void SetCompilerFlags(const char* flags); /** * @brief Returns true if the initialization was successfull */ virtual bool IsInitialized(); /** @brief Destructor */ virtual ~OclFilter(); protected: typedef std::vector CStringList; typedef std::vector ClSizeList; /** @brief Constructor */ OclFilter(); /** @brief Constructor ( overloaded ) */ OclFilter(const char* filename); /** @brief String that contains the compiler flags */ const char* m_ClCompilerFlags; /** @brief The compiled OpenCL program */ cl_program m_ClProgram; /** @brief Command queue for the filter */ cl_command_queue m_CommandQue; /** @brief Unique ID of the filter, needs to be specified in the constructor of the derived class */ std::string m_FilterID; /*! @brief source preambel for e.g. #define commands to be inserted into the OpenCL source */ const char* m_Preambel; /** @brief List of sourcefiles that will be compiled for this filter.*/ CStringList m_ClFiles; /** @brief status of the filter */ bool m_Initialized; /** @brief The local work size fo the filter */ size_t m_LocalWorkSize[3]; /** @brief The global work size of the filter */ size_t m_GlobalWorkSize[3]; /** @brief Set the working size for the following OpenCL kernel call */ void SetWorkingSize(unsigned int locx, unsigned int dimx, unsigned int locy = 1, unsigned int dimy = 1, unsigned int locz = 1, unsigned int dimz = 1); /** @brief Execute the given kernel on the OpenCL Index-Space defined by the local and global work sizes */ bool ExecuteKernel( cl_kernel kernel, unsigned int workSizeDim ); /** @brief Execute the given kernel on the OpenCL Index-Space defined by the local and global work sizes, but divide it into chunks of dimension chunksDim */ bool ExecuteKernelChunks( cl_kernel kernel, unsigned int workSizeDim, size_t* chunksDim ); + /** @brief Execute the given kernel on the OpenCL Index-Space defined by the local and global work sizes, but divide it into chunks of dimension chunksDim and wait between + * batches of batchSize chunks a time of waitTimems milliseconds + */ + bool ExecuteKernelChunksInBatches(cl_kernel kernel, unsigned int workSizeDim, size_t* chunksDim, size_t batchSize, int waitTimems); /** * \brief Initialize all necessary parts of the filter * * The Initialize() method creates the command queue and the m_clProgram. * The program is either compiled from the given source or taken from the * OclResourceManager if the program was compiled already. */ bool Initialize(); /** * @brief Compile the program source * * @param preambel e.g. defines for the shader code */ void CompileSource(); /** * @brief Add some source code on the beginning of the loaded source * * In this way, some preprocessor flags for the CL compiler can at the beginning of the filter * @param preambel Source preambel for e.g. #define commands to be inserted into the OpenCL source */ void SetSourcePreambel(const char* preambel); /** * @brief Get the Module of the filter. Needs to be implemented by every subclass. * The filter will load the OpenCL sourcefiles from this module context. */ virtual us::Module* GetModule() = 0; /** * @brief Helper functions that load sourcefiles from the module context in the Initialize function. * @param SourceCodeList holds the sourcecode for every file as string, the SourceCodeSizeList holst the * size of every file in bytes. */ void LoadSourceFiles(CStringList &SourceCodeList, ClSizeList &SourceCodeSizeList); }; } #endif // __mitkOclFilter_h diff --git a/Modules/PhotoacousticsAlgorithms/Resources/DMAS.cl b/Modules/PhotoacousticsAlgorithms/Resources/DMAS.cl index 51138f53c9..7ec678b6ac 100644 --- a/Modules/PhotoacousticsAlgorithms/Resources/DMAS.cl +++ b/Modules/PhotoacousticsAlgorithms/Resources/DMAS.cl @@ -1,85 +1,85 @@ /*=================================================================== 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. ===================================================================*/ __kernel void ckDMAS( __global float* dSource, // input image __global float* dDest, // output buffer __global unsigned short* usedLines, __global unsigned short* AddSamples, __constant float* apodArray, unsigned short apodArraySize, unsigned int inputL, unsigned int inputS, unsigned int Slices, unsigned int outputL, unsigned int outputS // parameters ) { // get thread identifier unsigned int globalPosX = get_global_id(0); unsigned int globalPosY = get_global_id(1); unsigned int globalPosZ = get_global_id(2); // terminate non-valid threads if ( globalPosX < outputL && globalPosY < outputS && globalPosZ < Slices ) { float l_i = (float)globalPosX / (float)outputL * (float)inputL; unsigned short curUsedLines = usedLines[globalPosY * 3 * outputL + 3 * globalPosX]; unsigned short minLine = usedLines[globalPosY * 3 * outputL + 3 * globalPosX + 1]; unsigned short maxLine = usedLines[globalPosY * 3 *outputL + 3 * globalPosX + 2]; float apod_mult = (float)apodArraySize / (float)curUsedLines; unsigned short Delay1 = 0; unsigned short Delay2 = 0; float output = 0; float mult = 0; float s_1 = 0; float s_2 = 0; - float sign = 0; + float apod_1 = 0; for (short l_s1 = minLine; l_s1 < maxLine; ++l_s1) { Delay1 = AddSamples[globalPosY * outputL + (int)(fabs(l_s1 - l_i)/(float)inputL * (float)outputL)]; - if (Delay1 < inputS && Delay1 >= 0) { + if (Delay1 < inputS && Delay1 >= 0) + { s_1 = dSource[(int)(globalPosZ * inputL * inputS + Delay1 * inputL + l_s1)]; - sign += s_1; + apod_1 = apodArray[(int)((l_s1 - minLine)*apod_mult)]; for (short l_s2 = l_s1 + 1; l_s2 < maxLine; ++l_s2) { Delay2 = AddSamples[globalPosY * outputL + (int)(fabs(l_s2 - l_i)/(float)inputL * (float)outputL)]; - if (Delay2 < inputS && Delay2 >= 0) { + if (Delay2 < inputS && Delay2 >= 0) + { s_2 = dSource[(int)(globalPosZ * inputL * inputS + Delay2 * inputL + l_s2)]; - mult = apodArray[(int)((l_s2 - minLine)*apod_mult)] * - s_2 - * apodArray[(int)((l_s1 - minLine)*apod_mult)] * - s_1; + mult = apodArray[(int)((l_s2 - minLine)*apod_mult)] * s_2 + * apod_1 * s_1; output += sqrt(fabs(mult)) * ((mult > 0) - (mult < 0)); } } } else --curUsedLines; } - dDest[ globalPosZ * outputL * outputS + globalPosY * outputL + globalPosX ] = output / (float)(curUsedLines * curUsedLines - (curUsedLines - 1)) * ((sign > 0) - (sign < 0)); + dDest[ globalPosZ * outputL * outputS + globalPosY * outputL + globalPosX ] = output / (float)(curUsedLines * curUsedLines - (curUsedLines - 1)); } } diff --git a/Modules/PhotoacousticsAlgorithms/Resources/DelayCalculation.cl b/Modules/PhotoacousticsAlgorithms/Resources/DelayCalculation.cl index 053e9e7722..afaf1e5c07 100644 --- a/Modules/PhotoacousticsAlgorithms/Resources/DelayCalculation.cl +++ b/Modules/PhotoacousticsAlgorithms/Resources/DelayCalculation.cl @@ -1,71 +1,71 @@ /*=================================================================== 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. ===================================================================*/ __kernel void ckDelayCalculationQuad( __global unsigned short *gDest, __global unsigned short *usedLines, unsigned int inputL, unsigned int inputS, unsigned int outputL, unsigned int outputS, char isPAImage, float delayMultiplicatorRaw // parameters ) { uint globalPosX = get_global_id(0); uint globalPosY = get_global_id(1); uint globalPosZ = get_global_id(2); - if (globalPosX < outputL && globalPosY < outputS) + if (globalPosX * 2 < outputL && globalPosY < outputS) { float l_i = 0; // we calculate the delays relative to line zero float s_i = (float)globalPosY / (float)outputS * (float)inputS / 2; float l_s = (float)globalPosX / (float)outputL * (float)inputL; // the currently calculated line float delayMultiplicator = delayMultiplicatorRaw / s_i; gDest[globalPosY * outputL + globalPosX] = delayMultiplicator * pow((l_s - l_i), 2) + s_i + (1-isPAImage)*s_i; } } __kernel void ckDelayCalculationSphe( __global unsigned short *gDest, __global unsigned short *usedLines, unsigned int inputL, unsigned int inputS, unsigned int outputL, unsigned int outputS, char isPAImage, float delayMultiplicatorRaw // parameters ) { uint globalPosX = get_global_id(0); uint globalPosY = get_global_id(1); uint globalPosZ = get_global_id(2); - if (globalPosX < outputL && globalPosY < outputS) + if (globalPosX * 2< outputL && globalPosY < outputS) { float l_i = 0; // we calculate the delays relative to line zero float s_i = (float)globalPosY / (float)outputS * (float)inputS / 2; float l_s = (float)globalPosX / (float)outputL * (float)inputL; // the currently calculated line gDest[globalPosY * outputL + globalPosX] = sqrt( pow(s_i, 2) + pow((delayMultiplicatorRaw * ((l_s - l_i)) / inputL), 2) ) + (1-isPAImage)*s_i; } } \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/Resources/UsedLinesCalculation.cl b/Modules/PhotoacousticsAlgorithms/Resources/UsedLinesCalculation.cl index a0d4af6fb7..638de2b279 100644 --- a/Modules/PhotoacousticsAlgorithms/Resources/UsedLinesCalculation.cl +++ b/Modules/PhotoacousticsAlgorithms/Resources/UsedLinesCalculation.cl @@ -1,48 +1,50 @@ /*=================================================================== 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. ===================================================================*/ __kernel void ckUsedLines( __global unsigned short* dDest, // output buffer float partMult, unsigned int inputL, - unsigned int inputS + unsigned int inputS, + unsigned int outputL, + unsigned int outputS ) { // get thread identifier unsigned int globalPosX = get_global_id(0); unsigned int globalPosY = get_global_id(1); - unsigned short outputS = get_global_size(1); - unsigned short outputL = get_global_size(0); + //unsigned short outputS = get_global_size(1); + //unsigned short outputL = get_global_size(0); // terminate non-valid threads if ( globalPosX < outputL && globalPosY < outputS) { float l_i = (float)globalPosX / outputL * inputL; float s_i = (float)globalPosY / outputS * inputS / 2; float part = partMult * s_i; if (part < 1) part = 1; unsigned short maxLine = min((l_i + part) + 1, (float)inputL); unsigned short minLine = max((l_i - part), 0.0f); dDest[globalPosY * 3 * outputL + 3 * globalPosX] = (maxLine - minLine); //usedLines dDest[globalPosY * 3 * outputL + 3 * globalPosX + 1] = minLine; //minLine dDest[globalPosY * 3 * outputL + 3 * globalPosX + 2] = maxLine; //maxLine } } \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/Resources/DMAS.cl b/Modules/PhotoacousticsAlgorithms/Resources/sDMAS.cl similarity index 91% copy from Modules/PhotoacousticsAlgorithms/Resources/DMAS.cl copy to Modules/PhotoacousticsAlgorithms/Resources/sDMAS.cl index 51138f53c9..1a8ff0269e 100644 --- a/Modules/PhotoacousticsAlgorithms/Resources/DMAS.cl +++ b/Modules/PhotoacousticsAlgorithms/Resources/sDMAS.cl @@ -1,85 +1,87 @@ /*=================================================================== 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. ===================================================================*/ -__kernel void ckDMAS( +__kernel void cksDMAS( __global float* dSource, // input image __global float* dDest, // output buffer __global unsigned short* usedLines, __global unsigned short* AddSamples, __constant float* apodArray, unsigned short apodArraySize, unsigned int inputL, unsigned int inputS, unsigned int Slices, unsigned int outputL, unsigned int outputS // parameters ) { // get thread identifier unsigned int globalPosX = get_global_id(0); unsigned int globalPosY = get_global_id(1); unsigned int globalPosZ = get_global_id(2); // terminate non-valid threads if ( globalPosX < outputL && globalPosY < outputS && globalPosZ < Slices ) { float l_i = (float)globalPosX / (float)outputL * (float)inputL; unsigned short curUsedLines = usedLines[globalPosY * 3 * outputL + 3 * globalPosX]; unsigned short minLine = usedLines[globalPosY * 3 * outputL + 3 * globalPosX + 1]; unsigned short maxLine = usedLines[globalPosY * 3 *outputL + 3 * globalPosX + 2]; float apod_mult = (float)apodArraySize / (float)curUsedLines; unsigned short Delay1 = 0; unsigned short Delay2 = 0; float output = 0; float mult = 0; float s_1 = 0; float s_2 = 0; float sign = 0; + float apod_1 = 0; for (short l_s1 = minLine; l_s1 < maxLine; ++l_s1) { Delay1 = AddSamples[globalPosY * outputL + (int)(fabs(l_s1 - l_i)/(float)inputL * (float)outputL)]; - if (Delay1 < inputS && Delay1 >= 0) { + if (Delay1 < inputS && Delay1 >= 0) + { s_1 = dSource[(int)(globalPosZ * inputL * inputS + Delay1 * inputL + l_s1)]; + apod_1 = apodArray[(int)((l_s1 - minLine)*apod_mult)]; sign += s_1; for (short l_s2 = l_s1 + 1; l_s2 < maxLine; ++l_s2) { Delay2 = AddSamples[globalPosY * outputL + (int)(fabs(l_s2 - l_i)/(float)inputL * (float)outputL)]; - if (Delay2 < inputS && Delay2 >= 0) { + if (Delay2 < inputS && Delay2 >= 0) + { s_2 = dSource[(int)(globalPosZ * inputL * inputS + Delay2 * inputL + l_s2)]; - mult = apodArray[(int)((l_s2 - minLine)*apod_mult)] * - s_2 - * apodArray[(int)((l_s1 - minLine)*apod_mult)] * - s_1; + mult = apodArray[(int)((l_s2 - minLine)*apod_mult)] * s_2 + * apod_1 * s_1; output += sqrt(fabs(mult)) * ((mult > 0) - (mult < 0)); } } } else --curUsedLines; } dDest[ globalPosZ * outputL * outputS + globalPosY * outputL + globalPosX ] = output / (float)(curUsedLines * curUsedLines - (curUsedLines - 1)) * ((sign > 0) - (sign < 0)); } } diff --git a/Modules/PhotoacousticsAlgorithms/files.cmake b/Modules/PhotoacousticsAlgorithms/files.cmake index 88f4bd6e46..ac29b12402 100644 --- a/Modules/PhotoacousticsAlgorithms/files.cmake +++ b/Modules/PhotoacousticsAlgorithms/files.cmake @@ -1,19 +1,20 @@ file(GLOB_RECURSE H_FILES RELATIVE "${CMAKE_CURRENT_SOURCE_DIR}" "${CMAKE_CURRENT_SOURCE_DIR}/include/*") set(CPP_FILES source/mitkPhotoacousticImage.cpp source/mitkPhotoacousticBeamformingFilter.cpp source/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.cpp source/OpenCLFilter/mitkPhotoacousticBModeFilter.cpp source/OpenCLFilter/mitkPhotoacousticOCLUsedLinesCalculation.cpp source/OpenCLFilter/mitkPhotoacousticOCLDelayCalculation.cpp ) set(RESOURCE_FILES BModeAbs.cl BModeAbsLog.cl UsedLinesCalculation.cl DelayCalculation.cl DMAS.cl DAS.cl + sDMAS.cl ) diff --git a/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLDelayCalculation.h b/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLDelayCalculation.h index febc14e0d5..061f7c7448 100644 --- a/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLDelayCalculation.h +++ b/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLDelayCalculation.h @@ -1,98 +1,99 @@ /*=================================================================== 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 _MITKPHOTOACOUSTICSDELAYCALC_H_ #define _MITKPHOTOACOUSTICSDELAYCALC_H_ #if defined(PHOTOACOUSTICS_USE_GPU) || DOXYGEN #include "mitkOclDataSetToDataSetFilter.h" #include #include "mitkPhotoacousticBeamformingSettings.h" namespace mitk { /*! * \brief Class implementing a mitk::OclDataSetToDataSetFilter to calculate the delays used for beamforming. * * The class must be given a configuration class instance of mitk::BeamformingSettings for beamforming parameters through mitk::OCLDelayCalculation::SetConfig(BeamformingSettings conf) * Additionally the output of an instance of mitk::OCLUsedLinesCalculation is needed to calculate the delays. */ class OCLDelayCalculation : public OclDataSetToDataSetFilter, public itk::Object { public: mitkClassMacroItkParent(OCLDelayCalculation, itk::Object); itkNewMacro(Self); void Update(); /** \brief Sets a new configuration to use. * * @param conf The configuration set to use for the calculation of the delays. */ void SetConfig(BeamformingSettings conf) { m_Conf = conf; } /** \brief Sets the usedLines buffer object to use for the calculation of the delays. * * @param usedLines An buffer generated as the output of an instance of mitk::OCLUsedLinesCalculation. */ void SetInputs(cl_mem usedLines) { m_UsedLines = usedLines; } protected: OCLDelayCalculation(); virtual ~OCLDelayCalculation(); /** Initialize the filter */ bool Initialize(); void Execute(); mitk::PixelType GetOutputType() { return mitk::MakeScalarPixelType(); } int GetBytesPerElem() { return sizeof(unsigned short); } virtual us::Module* GetModule(); int m_sizeThis; private: /** The OpenCL kernel for the filter */ cl_kernel m_PixelCalculation; BeamformingSettings m_Conf; cl_mem m_UsedLines; unsigned int m_BufferSize; float m_DelayMultiplicatorRaw; char m_IsPAImage; + size_t m_ChunkSize[3]; }; } #endif #endif \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLUsedLinesCalculation.h b/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLUsedLinesCalculation.h index 6ee3cf0ccf..e787d2009a 100644 --- a/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLUsedLinesCalculation.h +++ b/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLUsedLinesCalculation.h @@ -1,89 +1,90 @@ /*=================================================================== 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 _MITKPHOTOACOUSTICSOCLUSEDLINESCALCULATION_H_ #define _MITKPHOTOACOUSTICSOCLUSEDLINESCALCULATION_H_ #if defined(PHOTOACOUSTICS_USE_GPU) || DOXYGEN #include "mitkOclDataSetToDataSetFilter.h" #include #include "mitkPhotoacousticBeamformingSettings.h" namespace mitk { /*! * \brief Class implementing a mitk::OclDataSetToDataSetFilter to calculate which lines each sample should use when beamforming. * * The class must be given a configuration class instance of mitk::BeamformingSettings for beamforming parameters through mitk::OCLDelayCalculation::SetConfig(BeamformingSettings conf) */ class OCLUsedLinesCalculation : public OclDataSetToDataSetFilter, public itk::Object { public: mitkClassMacroItkParent(OCLUsedLinesCalculation, itk::Object); itkNewMacro(Self); void Update(); /** \brief Sets a new configuration to use. * * @param conf The configuration set to use for the calculation of the used lines. */ void SetConfig(BeamformingSettings conf) { m_Conf = conf; } protected: /** Constructor */ OCLUsedLinesCalculation(); /** Destructor */ virtual ~OCLUsedLinesCalculation(); /** Initialize the filter */ bool Initialize(); void Execute(); mitk::PixelType GetOutputType() { return mitk::MakeScalarPixelType(); } int GetBytesPerElem() { return sizeof(unsigned short); } virtual us::Module* GetModule(); int m_sizeThis; private: /** The OpenCL kernel for the filter */ cl_kernel m_PixelCalculation; BeamformingSettings m_Conf; float m_part; + size_t m_ChunkSize[3]; }; } #endif #endif diff --git a/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticBeamformingFilter.h b/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticBeamformingFilter.h index 076724c482..214ef84e5d 100644 --- a/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticBeamformingFilter.h +++ b/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticBeamformingFilter.h @@ -1,145 +1,151 @@ /*=================================================================== 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 MITK_PHOTOACOUSTICS_BEAMFORMING_FILTER #define MITK_PHOTOACOUSTICS_BEAMFORMING_FILTER #include "mitkImageToImageFilter.h" #include #include "./OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.h" #include "mitkPhotoacousticBeamformingSettings.h" namespace mitk { /*! * \brief Class implementing an mitk::ImageToImageFilter for beamforming on both CPU and GPU * * The class must be given a configuration class instance of mitk::BeamformingSettings for beamforming parameters through mitk::BeamformingFilter::Configure(BeamformingSettings settings) * Whether the GPU is used can be set in the configuration. * For significant problems or important messages a string is written, which can be accessed via GetMessageString(). */ class BeamformingFilter : public ImageToImageFilter { public: mitkClassMacro(BeamformingFilter, ImageToImageFilter); itkFactorylessNewMacro(Self) itkCloneMacro(Self) /** \brief Sets a new configuration to use * * @param settings The configuration set to use for beamforming */ void Configure(BeamformingSettings settings) { m_ConfOld = m_Conf; m_Conf = settings; } /** \brief Sets a new configuration to use * * The Filter writes important messages that can be retrieved through this method; if nothing is to be reported, it returns "noMessage". * @return The message */ std::string GetMessageString() { return m_Message; } /** \brief Sets a callback for progress checking * * An std::function can be set, through which progress of the currently updating filter is reported. * The integer argument is a number between 0 an 100 to indicate how far completion has been achieved, the std::string argument indicates what the filter is currently doing. */ void SetProgressHandle(std::function progressHandle); protected: BeamformingFilter(); ~BeamformingFilter(); virtual void GenerateInputRequestedRegion() override; virtual void GenerateOutputInformation() override; virtual void GenerateData() override; //##Description //## @brief Time when Header was last initialized itk::TimeStamp m_TimeOfHeaderInitialization; /** \brief The std::function, through which progress of the currently updating filter is reported. */ std::function m_ProgressHandle; /** \brief Pointer holding the Von-Hann apodization window for beamforming * @param samples the resolution at which the window is created */ float* VonHannFunction(int samples); /** \brief Function to create a Hamming apodization window * @param samples the resolution at which the window is created */ float* HammFunction(int samples); /** \brief Function to create a Box apodization window * @param samples the resolution at which the window is created */ float* BoxFunction(int samples); /** \brief Function to perform beamforming on CPU for a single line, using DAS and quadratic delay */ void DASQuadraticLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize); /** \brief Function to perform beamforming on CPU for a single line, using DAS and spherical delay */ void DASSphericalLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize); /** \brief Function to perform beamforming on CPU for a single line, using DMAS and quadratic delay */ void DMASQuadraticLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize); /** \brief Function to perform beamforming on CPU for a single line, using DMAS and spherical delay */ void DMASSphericalLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize); + /** \brief Function to perform beamforming on CPU for a single line, using signed DMAS and quadratic delay + */ + void sDMASQuadraticLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize); + /** \brief Function to perform beamforming on CPU for a single line, using signed DMAS and spherical delay + */ + void sDMASSphericalLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize); float* m_OutputData; float* m_InputData; float* m_InputDataPuffer; /** \brief Pointer holding the Von-Hann apodization window for beamforming */ float* m_VonHannFunction; /** \brief Pointer holding the Hamming apodization window for beamforming */ float* m_HammFunction; /** \brief Pointer holding the Box apodization window for beamforming */ float* m_BoxFunction; /** \brief Current configuration set */ BeamformingSettings m_Conf; /** \brief Previous configuration set to selectively update used data for beamforming */ BeamformingSettings m_ConfOld; /** \brief Pointer to the GPU beamforming filter class; for performance reasons the filter is initialized within the constructor and kept for all later computations. */ mitk::PhotoacousticOCLBeamformingFilter::Pointer m_BeamformingOclFilter; /** \brief The message returned by mitk::BeamformingFilter::GetMessageString() */ std::string m_Message; }; } // namespace mitk #endif //MITK_PHOTOACOUSTICS_BEAMFORMING_FILTER diff --git a/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticBeamformingSettings.h b/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticBeamformingSettings.h index ff44b554db..953a7b1431 100644 --- a/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticBeamformingSettings.h +++ b/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticBeamformingSettings.h @@ -1,137 +1,137 @@ /*=================================================================== 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 MITK_PHOTOACOUSTICS_BEAMFORMING_SETTINGS #define MITK_PHOTOACOUSTICS_BEAMFORMING_SETTINGS namespace mitk { /*! * \brief Class holding the configuration data for the beamforming filters mitk::BeamformingFilter and mitk::PhotoacousticOCLBeamformingFilter * * A detailed description can be seen below. All parameters should be set manually for successfull beamforming. */ class BeamformingSettings { public: /** \brief Pitch of the used transducer in [m]. */ float Pitch = 0.0003; /** \brief Speed of sound in the used medium in [m/s]. */ float SpeedOfSound = 1540; /** \brief This parameter is not neccessary to be set, as it's never used. */ float RecordTime = 0.00006; // [s] /** \brief The time spacing of the input image */ float TimeSpacing = 0.0000000000001; // [s] /** \brief The angle of the transducer elements */ float Angle = -1; /** \brief Flag whether processed image is a photoacoustic image or an ultrasound image */ bool isPhotoacousticImage = true; /** \brief How many transducer elements the used transducer had. */ unsigned short TransducerElements = 128; /** \brief How many vertical samples should be used in the final image. */ unsigned int SamplesPerLine = 2048; /** \brief How many lines should be reconstructed in the final image. */ unsigned int ReconstructionLines = 128; /** \brief Sets how many voxels should be cut off from the top of the image before beamforming, to potentially avoid artifacts. */ unsigned int upperCutoff = 0; /** \brief Sets whether only the slices selected by mitk::BeamformingSettings::CropBounds should be beamformed. */ bool partial = false; /** \brief Sets the first and last slice to be beamformed. */ unsigned int CropBounds[2] = { 0,0 }; /** \brief Sets the dimensions of the inputImage. */ unsigned int inputDim[3] = { 1,1,1 }; /** \brief Decides whether GPU computing should be used */ bool UseGPU = true; /** \brief Available delay calculation methods: * - Spherical delay for best results. * - A quadratic Taylor approximation for slightly faster results with hardly any quality loss. */ enum DelayCalc { QuadApprox, Spherical }; /** \brief Sets how the delays for beamforming should be calculated. */ DelayCalc DelayCalculationMethod = QuadApprox; /** \brief Available apodization functions: * - Hamming function. * - Von-Hann function. * - Box function. */ enum Apodization { Hamm, Hann, Box }; /** \brief Sets the used apodization function. */ Apodization Apod = Hann; /** \brief Sets the resolution of the apodization array (must be greater than 0). */ int apodizationArraySize = 128; /** \brief Available beamforming algorithms: * - DAS (Delay and sum). * - DMAS (Delay multiply and sum). */ - enum BeamformingAlgorithm { DMAS, DAS }; + enum BeamformingAlgorithm { DMAS, DAS, sDMAS}; /** \brief Sets the used beamforming algorithm. */ BeamformingAlgorithm Algorithm = DAS; /** \brief Sets whether after beamforming a bandpass should be automatically applied */ bool UseBP = false; /** \brief Sets the position at which lower frequencies are completely cut off in Hz. */ float BPHighPass = 50; /** \brief Sets the position at which higher frequencies are completely cut off in Hz. */ float BPLowPass = 50; /** \brief function for mitk::PhotoacousticOCLBeamformingFilter to check whether buffers need to be updated * this method ignores changes in BPLow/BPHigh/cropBounds/Algorithm/some more, as those are insignifiant in all current situations */ static bool SettingsChangedOpenCL(const BeamformingSettings& lhs, const BeamformingSettings& rhs) { - return !((abs(lhs.Angle - rhs.Angle) < 0.001f) && + return !((abs(lhs.Angle - rhs.Angle) < 0.01f) && // 0.01 degree error margin (lhs.Apod == rhs.Apod) && (lhs.DelayCalculationMethod == rhs.DelayCalculationMethod) && (lhs.isPhotoacousticImage == rhs.isPhotoacousticImage) && - (abs(lhs.Pitch - rhs.Pitch) < 0.0000001f) && + (abs(lhs.Pitch - rhs.Pitch) < 0.000001f) && // 0.0001 mm error margin (lhs.ReconstructionLines == rhs.ReconstructionLines) && - (abs(lhs.RecordTime - rhs.RecordTime) < 0.00000001f) && + (abs(lhs.RecordTime - rhs.RecordTime) < 0.00000001f) && // 10 ns error margin (lhs.SamplesPerLine == rhs.SamplesPerLine) && (abs(lhs.SpeedOfSound - rhs.SpeedOfSound) < 0.01f) && - (abs(lhs.TimeSpacing - rhs.TimeSpacing) < 0.0000000001f) && + (abs(lhs.TimeSpacing - rhs.TimeSpacing) < 0.00000000001f) && //0.01 ns error margin (lhs.TransducerElements == rhs.TransducerElements)); } }; } #endif \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.cpp b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.cpp index 86f638f4fc..6da6fbaa9b 100644 --- a/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.cpp @@ -1,259 +1,265 @@ /*=================================================================== 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. ===================================================================*/ #if defined(PHOTOACOUSTICS_USE_GPU) || DOXYGEN #include "./OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.h" #include "usServiceReference.h" mitk::PhotoacousticOCLBeamformingFilter::PhotoacousticOCLBeamformingFilter() : m_PixelCalculation( NULL ), m_InputImage(mitk::Image::New()), m_ApodizationBuffer(nullptr), m_MemoryLocationsBuffer(nullptr), m_DelaysBuffer(nullptr), m_UsedLinesBuffer(nullptr) { this->AddSourceFile("DAS.cl"); this->AddSourceFile("DMAS.cl"); + this->AddSourceFile("sDMAS.cl"); this->m_FilterID = "OpenCLBeamformingFilter"; this->Initialize(); unsigned int dim[] = { 128, 2048, 2 }; mitk::Vector3D spacing; spacing[0] = 1; spacing[1] = 1; spacing[2] = 1; m_InputImage->Initialize(mitk::MakeScalarPixelType(), 3, dim); m_InputImage->SetSpacing(spacing); m_ChunkSize[0] = 128; m_ChunkSize[1] = 128; m_ChunkSize[2] = 8; m_UsedLinesCalculation = mitk::OCLUsedLinesCalculation::New(); m_DelayCalculation = mitk::OCLDelayCalculation::New(); } mitk::PhotoacousticOCLBeamformingFilter::~PhotoacousticOCLBeamformingFilter() { if ( this->m_PixelCalculation ) { clReleaseKernel( m_PixelCalculation ); } if (m_ApodizationBuffer) clReleaseMemObject(m_ApodizationBuffer); } void mitk::PhotoacousticOCLBeamformingFilter::Update() { //Check if context & program available if (!this->Initialize()) { us::ServiceReference ref = GetModuleContext()->GetServiceReference(); OclResourceService* resources = GetModuleContext()->GetService(ref); // clean-up also the resources resources->InvalidateStorage(); mitkThrow() <<"Filter is not initialized. Cannot update."; } else{ // Execute this->Execute(); } } void mitk::PhotoacousticOCLBeamformingFilter::UpdateDataBuffers() { /*us::ServiceReference ref = GetModuleContext()->GetServiceReference(); OclResourceService* resources = GetModuleContext()->GetService(ref); cl_ulong globalMemSize = oclGetGlobalMemSize(resources->GetCurrentDevice());*/ //Initialize the Output try { MITK_DEBUG << "Updating Workgroup size for new dimensions"; size_t outputSize = (size_t)m_Conf.ReconstructionLines * (size_t)m_Conf.SamplesPerLine * (size_t)m_Conf.inputDim[2]; m_OutputDim[0] = m_Conf.ReconstructionLines; m_OutputDim[1] = m_Conf.SamplesPerLine; m_OutputDim[2] = m_Conf.inputDim[2]; this->InitExec(this->m_PixelCalculation, m_OutputDim, outputSize, sizeof(float)); } catch (const mitk::Exception& e) { MITK_ERROR << "Caught exception while initializing filter: " << e.what(); return; } if (BeamformingSettings::SettingsChangedOpenCL(m_Conf, m_ConfOld)) { cl_int clErr = 0; MITK_DEBUG << "Updating GPU Buffers for new configuration"; // create the apodisation buffer if (m_Apodisation == nullptr) { MITK_INFO << "No apodisation function set; Beamforming will be done without any apodisation."; m_Apodisation = new float[1]; m_Apodisation[0] = 1; m_ApodArraySize = 1; } us::ServiceReference ref = GetModuleContext()->GetServiceReference(); OclResourceService* resources = GetModuleContext()->GetService(ref); cl_context gpuContext = resources->GetContext(); if (m_ApodizationBuffer) clReleaseMemObject(m_ApodizationBuffer); this->m_ApodizationBuffer = clCreateBuffer(gpuContext, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof(float) * m_ApodArraySize, m_Apodisation, &clErr); CHECK_OCL_ERR(clErr); // calculate used lines m_UsedLinesCalculation->SetConfig(m_Conf); m_UsedLinesCalculation->Update(); m_UsedLinesBuffer = m_UsedLinesCalculation->GetGPUOutput()->GetGPUBuffer(); // calculate the Delays m_DelayCalculation->SetConfig(m_Conf); m_DelayCalculation->SetInputs(m_UsedLinesBuffer); m_DelayCalculation->Update(); m_DelaysBuffer = m_DelayCalculation->GetGPUOutput()->GetGPUBuffer(); m_ConfOld = m_Conf; } } void mitk::PhotoacousticOCLBeamformingFilter::Execute() { cl_int clErr = 0; UpdateDataBuffers(); clErr = clSetKernelArg(this->m_PixelCalculation, 2, sizeof(cl_mem), &(this->m_UsedLinesBuffer)); clErr |= clSetKernelArg(this->m_PixelCalculation, 3, sizeof(cl_mem), &(this->m_DelaysBuffer)); clErr |= clSetKernelArg(this->m_PixelCalculation, 4, sizeof(cl_mem), &(this->m_ApodizationBuffer)); clErr |= clSetKernelArg(this->m_PixelCalculation, 5, sizeof(cl_ushort), &(this->m_ApodArraySize)); clErr |= clSetKernelArg(this->m_PixelCalculation, 6, sizeof(cl_uint), &(this->m_Conf.inputDim[0])); clErr |= clSetKernelArg(this->m_PixelCalculation, 7, sizeof(cl_uint), &(this->m_Conf.inputDim[1])); clErr |= clSetKernelArg(this->m_PixelCalculation, 8, sizeof(cl_uint), &(this->m_Conf.inputDim[2])); clErr |= clSetKernelArg(this->m_PixelCalculation, 9, sizeof(cl_uint), &(this->m_Conf.ReconstructionLines)); clErr |= clSetKernelArg(this->m_PixelCalculation, 10, sizeof(cl_uint), &(this->m_Conf.SamplesPerLine)); // execute the filter on a 3D NDRange if (m_OutputDim[2] == 1 || m_ChunkSize[2] == 1) { - if(!this->ExecuteKernelChunks(m_PixelCalculation, 2, m_ChunkSize)) + if(!this->ExecuteKernelChunksInBatches(m_PixelCalculation, 2, m_ChunkSize, 16, 50)) mitkThrow() << "openCL Error when executing Kernel"; } else { - if(!this->ExecuteKernelChunks(m_PixelCalculation, 3, m_ChunkSize)) + if(!this->ExecuteKernelChunksInBatches(m_PixelCalculation, 3, m_ChunkSize, 16, 50)) mitkThrow() << "openCL Error when executing Kernel"; } // signalize the GPU-side data changed m_Output->Modified( GPU_DATA ); } us::Module *mitk::PhotoacousticOCLBeamformingFilter::GetModule() { return us::GetModuleContext()->GetModule(); } bool mitk::PhotoacousticOCLBeamformingFilter::Initialize() { bool buildErr = true; cl_int clErr = 0; if ( OclFilter::Initialize() ) { switch (m_Conf.Algorithm) { case BeamformingSettings::BeamformingAlgorithm::DAS: { this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDAS", &clErr); break; } case BeamformingSettings::BeamformingAlgorithm::DMAS: { this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDMAS", &clErr); break; } + case BeamformingSettings::BeamformingAlgorithm::sDMAS: + { + this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "cksDMAS", &clErr); + break; + } default: { MITK_INFO << "No beamforming algorithm specified, setting to DAS"; this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDAS", &clErr); break; } } buildErr |= CHECK_OCL_ERR( clErr ); } CHECK_OCL_ERR(clErr); return (OclFilter::IsInitialized() && buildErr ); } void mitk::PhotoacousticOCLBeamformingFilter::SetInput(mitk::Image::Pointer image) { OclDataSetToDataSetFilter::SetInput(image); m_InputImage = image; m_Conf.inputDim[0] = m_InputImage->GetDimension(0); m_Conf.inputDim[1] = m_InputImage->GetDimension(1); m_Conf.inputDim[2] = m_InputImage->GetDimension(2); } void mitk::PhotoacousticOCLBeamformingFilter::SetInput(void* data, unsigned int* dimensions, unsigned int BpE) { OclDataSetToDataSetFilter::SetInput(data, dimensions[0] * dimensions[1] * dimensions[2], BpE); m_Conf.inputDim[0] = dimensions[0]; m_Conf.inputDim[1] = dimensions[1]; m_Conf.inputDim[2] = dimensions[2]; } mitk::Image::Pointer mitk::PhotoacousticOCLBeamformingFilter::GetOutputAsImage() { mitk::Image::Pointer outputImage = mitk::Image::New(); if (m_Output->IsModified(GPU_DATA)) { void* pData = m_Output->TransferDataToCPU(m_CommandQue); const unsigned int dimension = 3; unsigned int dimensions[3] = { (unsigned int)m_OutputDim[0], (unsigned int)m_OutputDim[1], (unsigned int)m_OutputDim[2] }; const mitk::SlicedGeometry3D::Pointer p_slg = m_InputImage->GetSlicedGeometry(); MITK_DEBUG << "Creating new MITK Image."; outputImage->Initialize(this->GetOutputType(), dimension, dimensions); outputImage->SetSpacing(p_slg->GetSpacing()); outputImage->SetImportVolume(pData, 0, 0, mitk::Image::ImportMemoryManagementType::ManageMemory); } MITK_DEBUG << "Image Initialized."; return outputImage; } void* mitk::PhotoacousticOCLBeamformingFilter::GetOutput() { return OclDataSetToDataSetFilter::GetOutput(); } #endif diff --git a/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLDelayCalculation.cpp b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLDelayCalculation.cpp index 94d28d6a33..53019b8b57 100644 --- a/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLDelayCalculation.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLDelayCalculation.cpp @@ -1,124 +1,124 @@ /*=================================================================== 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. ===================================================================*/ #define _USE_MATH_DEFINES #include #include "./OpenCLFilter/mitkPhotoacousticOCLDelayCalculation.h" #include "usServiceReference.h" #include "mitkImageReadAccessor.h" mitk::OCLDelayCalculation::OCLDelayCalculation() : m_PixelCalculation(NULL) { this->AddSourceFile("DelayCalculation.cl"); this->m_FilterID = "DelayCalculation"; + m_ChunkSize[0] = 128; + m_ChunkSize[1] = 128; + m_ChunkSize[2] = 8; + this->Initialize(); } mitk::OCLDelayCalculation::~OCLDelayCalculation() { if (this->m_PixelCalculation) { clReleaseKernel(m_PixelCalculation); } } void mitk::OCLDelayCalculation::Update() { //Check if context & program available if (!this->Initialize()) { us::ServiceReference ref = GetModuleContext()->GetServiceReference(); OclResourceService* resources = GetModuleContext()->GetService(ref); // clean-up also the resources resources->InvalidateStorage(); mitkThrow() << "Filter is not initialized. Cannot update."; } else { // Execute this->Execute(); } } void mitk::OCLDelayCalculation::Execute() { cl_int clErr = 0; - unsigned int gridDim[3] = { m_Conf.ReconstructionLines, m_Conf.SamplesPerLine, 1 }; + unsigned int gridDim[3] = { m_Conf.ReconstructionLines / 2, m_Conf.SamplesPerLine, 1 }; m_BufferSize = gridDim[0] * gridDim[1] * 1; try { this->InitExecNoInput(this->m_PixelCalculation, gridDim, m_BufferSize, sizeof(unsigned short)); } catch (const mitk::Exception& e) { MITK_ERROR << "Caught exception while initializing filter: " << e.what(); return; } if (m_Conf.DelayCalculationMethod == BeamformingSettings::DelayCalc::QuadApprox) m_DelayMultiplicatorRaw = pow(1 / (m_Conf.TimeSpacing*m_Conf.SpeedOfSound) * m_Conf.Pitch * (float)m_Conf.TransducerElements / (float)m_Conf.inputDim[0], 2) / 2; else if (m_Conf.DelayCalculationMethod == BeamformingSettings::DelayCalc::Spherical) m_DelayMultiplicatorRaw = 1 / (m_Conf.TimeSpacing*m_Conf.SpeedOfSound) * (m_Conf.Pitch*(float)m_Conf.TransducerElements); m_IsPAImage = m_Conf.isPhotoacousticImage; clErr = clSetKernelArg(this->m_PixelCalculation, 1, sizeof(cl_mem), &(this->m_UsedLines)); clErr |= clSetKernelArg(this->m_PixelCalculation, 2, sizeof(cl_uint), &(this->m_Conf.inputDim[0])); clErr |= clSetKernelArg(this->m_PixelCalculation, 3, sizeof(cl_uint), &(this->m_Conf.inputDim[1])); clErr |= clSetKernelArg(this->m_PixelCalculation, 4, sizeof(cl_uint), &(this->m_Conf.ReconstructionLines)); clErr |= clSetKernelArg(this->m_PixelCalculation, 5, sizeof(cl_uint), &(this->m_Conf.SamplesPerLine)); clErr |= clSetKernelArg(this->m_PixelCalculation, 6, sizeof(cl_char), &(this->m_IsPAImage)); clErr |= clSetKernelArg(this->m_PixelCalculation, 7, sizeof(cl_float), &(this->m_DelayMultiplicatorRaw)); CHECK_OCL_ERR(clErr); - size_t ChunkSize[3]; - - ChunkSize[0] = 128; - ChunkSize[1] = 128; - ChunkSize[2] = 8; // execute the filter on a 3D NDRange - this->ExecuteKernelChunks(m_PixelCalculation, 2, ChunkSize); + if (!this->ExecuteKernelChunksInBatches(m_PixelCalculation, 2, m_ChunkSize, 16, 50)) + mitkThrow() << "openCL Error when executing Kernel"; // signalize the GPU-side data changed m_Output->Modified(GPU_DATA); } us::Module *mitk::OCLDelayCalculation::GetModule() { return us::GetModuleContext()->GetModule(); } bool mitk::OCLDelayCalculation::Initialize() { bool buildErr = true; cl_int clErr = 0; if (OclFilter::Initialize()) { if(m_Conf.DelayCalculationMethod == BeamformingSettings::DelayCalc::QuadApprox) this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDelayCalculationQuad", &clErr); if (m_Conf.DelayCalculationMethod == BeamformingSettings::DelayCalc::Spherical) this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDelayCalculationSphe", &clErr); buildErr |= CHECK_OCL_ERR(clErr); } return (OclFilter::IsInitialized() && buildErr); } \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLUsedLinesCalculation.cpp b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLUsedLinesCalculation.cpp index 349c4781c4..44b26beda3 100644 --- a/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLUsedLinesCalculation.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLUsedLinesCalculation.cpp @@ -1,108 +1,123 @@ /*=================================================================== 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. ===================================================================*/ #if defined(PHOTOACOUSTICS_USE_GPU) || DOXYGEN #define _USE_MATH_DEFINES #include #include "./OpenCLFilter/mitkPhotoacousticOCLUsedLinesCalculation.h" #include "usServiceReference.h" #include "mitkImageReadAccessor.h" mitk::OCLUsedLinesCalculation::OCLUsedLinesCalculation() : m_PixelCalculation(NULL) { this->AddSourceFile("UsedLinesCalculation.cl"); this->m_FilterID = "UsedLinesCalculation"; + m_ChunkSize[0] = 128; + m_ChunkSize[1] = 128; + m_ChunkSize[2] = 8; + this->Initialize(); } mitk::OCLUsedLinesCalculation::~OCLUsedLinesCalculation() { if (this->m_PixelCalculation) { clReleaseKernel(m_PixelCalculation); } } void mitk::OCLUsedLinesCalculation::Update() { //Check if context & program available if (!this->Initialize()) { us::ServiceReference ref = GetModuleContext()->GetServiceReference(); OclResourceService* resources = GetModuleContext()->GetService(ref); // clean-up also the resources resources->InvalidateStorage(); mitkThrow() << "Filter is not initialized. Cannot update."; } else { // Execute this->Execute(); } } void mitk::OCLUsedLinesCalculation::Execute() { cl_int clErr = 0; unsigned int gridDim[3] = { m_Conf.ReconstructionLines, m_Conf.SamplesPerLine, 1 }; size_t outputSize = gridDim[0] * gridDim[1] * 3; try { this->InitExecNoInput(this->m_PixelCalculation, gridDim, outputSize, sizeof(unsigned short)); } catch (const mitk::Exception& e) { MITK_ERROR << "Caught exception while initializing filter: " << e.what(); return; } m_part = (tan(m_Conf.Angle / 360 * 2 * M_PI) * ((m_Conf.SpeedOfSound * m_Conf.TimeSpacing)) / (m_Conf.Pitch * m_Conf.TransducerElements)) * m_Conf.inputDim[0]; clErr = clSetKernelArg(this->m_PixelCalculation, 1, sizeof(cl_float), &(this->m_part)); clErr |= clSetKernelArg(this->m_PixelCalculation, 2, sizeof(cl_uint), &(this->m_Conf.inputDim[0])); clErr |= clSetKernelArg(this->m_PixelCalculation, 3, sizeof(cl_uint), &(this->m_Conf.inputDim[1])); + clErr |= clSetKernelArg(this->m_PixelCalculation, 4, sizeof(cl_uint), &(this->m_Conf.ReconstructionLines)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 5, sizeof(cl_uint), &(this->m_Conf.SamplesPerLine)); CHECK_OCL_ERR(clErr); - // execute the filter on a 3D NDRange - this->ExecuteKernel(m_PixelCalculation, 2); + // execute the filter on a 2D NDRange + if (m_Conf.inputDim[2] > 2) + { + if (!this->ExecuteKernelChunksInBatches(m_PixelCalculation, 2, m_ChunkSize, 16, 50)) + mitkThrow() << "openCL Error when executing Kernel"; + } + else + { + if (!this->ExecuteKernelChunks(m_PixelCalculation, 2, m_ChunkSize)) + mitkThrow() << "openCL Error when executing Kernel"; + } // signalize the GPU-side data changed m_Output->Modified(GPU_DATA); } us::Module *mitk::OCLUsedLinesCalculation::GetModule() { return us::GetModuleContext()->GetModule(); } bool mitk::OCLUsedLinesCalculation::Initialize() { bool buildErr = true; cl_int clErr = 0; if (OclFilter::Initialize()) { this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckUsedLines", &clErr); buildErr |= CHECK_OCL_ERR(clErr); } return (OclFilter::IsInitialized() && buildErr); } #endif diff --git a/Modules/PhotoacousticsAlgorithms/source/mitkPhotoacousticBeamformingFilter.cpp b/Modules/PhotoacousticsAlgorithms/source/mitkPhotoacousticBeamformingFilter.cpp index ffa15eebf1..dae210a911 100644 --- a/Modules/PhotoacousticsAlgorithms/source/mitkPhotoacousticBeamformingFilter.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/mitkPhotoacousticBeamformingFilter.cpp @@ -1,577 +1,752 @@ /*=================================================================== mitkPhotoacousticBeamformingFilter 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. ===================================================================*/ #define _USE_MATH_DEFINES #include "mitkProperties.h" #include "mitkImageReadAccessor.h" #include #include #include #include #include #include #include "mitkImageCast.h" #include "mitkPhotoacousticBeamformingFilter.h" mitk::BeamformingFilter::BeamformingFilter() : m_OutputData(nullptr), m_InputData(nullptr), m_Message("noMessage") { this->SetNumberOfIndexedInputs(1); this->SetNumberOfRequiredInputs(1); m_ProgressHandle = [](int, std::string) {}; m_BeamformingOclFilter = mitk::PhotoacousticOCLBeamformingFilter::New(); m_VonHannFunction = VonHannFunction(m_Conf.apodizationArraySize); m_HammFunction = HammFunction(m_Conf.apodizationArraySize); m_BoxFunction = BoxFunction(m_Conf.apodizationArraySize); } void mitk::BeamformingFilter::SetProgressHandle(std::function progressHandle) { m_ProgressHandle = progressHandle; } mitk::BeamformingFilter::~BeamformingFilter() { delete[] m_VonHannFunction; delete[] m_HammFunction; delete[] m_BoxFunction; } void mitk::BeamformingFilter::GenerateInputRequestedRegion() { Superclass::GenerateInputRequestedRegion(); mitk::Image* output = this->GetOutput(); mitk::Image* input = const_cast (this->GetInput()); if (!output->IsInitialized()) { return; } input->SetRequestedRegionToLargestPossibleRegion(); //GenerateTimeInInputRegion(output, input); } void mitk::BeamformingFilter::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 dim[] = { m_Conf.ReconstructionLines, m_Conf.SamplesPerLine, input->GetDimension(2) }; output->Initialize(mitk::MakeScalarPixelType(), 3, dim); mitk::Vector3D spacing; spacing[0] = m_Conf.Pitch * m_Conf.TransducerElements * 1000 / m_Conf.ReconstructionLines; spacing[1] = (m_Conf.TimeSpacing * m_Conf.inputDim[1]) / 2 * m_Conf.SpeedOfSound * 1000 / m_Conf.SamplesPerLine; spacing[2] = 1; output->GetGeometry()->SetSpacing(spacing); output->GetGeometry()->Modified(); output->SetPropertyList(input->GetPropertyList()->Clone()); m_TimeOfHeaderInitialization.Modified(); } void mitk::BeamformingFilter::GenerateData() { GenerateOutputInformation(); mitk::Image::Pointer input = this->GetInput(); mitk::Image::Pointer output = this->GetOutput(); if (!output->IsInitialized()) return; float* ApodWindow; if (m_ConfOld.apodizationArraySize != m_Conf.apodizationArraySize) { delete[] m_VonHannFunction; delete[] m_HammFunction; delete[] m_BoxFunction; m_VonHannFunction = VonHannFunction(m_Conf.apodizationArraySize); m_HammFunction = HammFunction(m_Conf.apodizationArraySize); m_BoxFunction = BoxFunction(m_Conf.apodizationArraySize); m_ConfOld = m_Conf; } // set the appropiate apodization window switch (m_Conf.Apod) { case BeamformingSettings::Apodization::Hann: ApodWindow = m_VonHannFunction; break; case BeamformingSettings::Apodization::Hamm: ApodWindow = m_HammFunction; break; case BeamformingSettings::Apodization::Box: ApodWindow = m_BoxFunction; break; default: ApodWindow = m_BoxFunction; break; } auto begin = std::chrono::high_resolution_clock::now(); // debbuging the performance... if (!m_Conf.UseGPU) { int progInterval = output->GetDimension(2) / 20 > 1 ? output->GetDimension(2) / 20 : 1; // the interval at which we update the gui progress bar float inputDim[2] = { (float)input->GetDimension(0), (float)input->GetDimension(1) }; float outputDim[2] = { (float)output->GetDimension(0), (float)output->GetDimension(1) }; for (unsigned int i = 0; i < output->GetDimension(2); ++i) // seperate Slices should get Beamforming seperately applied { mitk::ImageReadAccessor inputReadAccessor(input, input->GetSliceData(i)); // first, we check whether the dara is float, other formats are unsupported if (input->GetPixelType().GetTypeAsString() == "scalar (float)" || input->GetPixelType().GetTypeAsString() == " (float)") { m_InputData = (float*)inputReadAccessor.GetData(); } else { MITK_INFO << "Pixel type is not float, abort"; return; } m_OutputData = new float[m_Conf.ReconstructionLines*m_Conf.SamplesPerLine]; // fill the image with zeros for (int l = 0; l < outputDim[0]; ++l) { for (int s = 0; s < outputDim[1]; ++s) { m_OutputData[l*(short)outputDim[1] + s] = 0; } } std::thread *threads = new std::thread[(short)outputDim[0]]; // every line will be beamformed in a seperate thread if (m_Conf.Algorithm == BeamformingSettings::BeamformingAlgorithm::DAS) { if (m_Conf.DelayCalculationMethod == BeamformingSettings::DelayCalc::QuadApprox) { for (short line = 0; line < outputDim[0]; ++line) { threads[line] = std::thread(&BeamformingFilter::DASQuadraticLine, this, m_InputData, m_OutputData, inputDim, outputDim, line, ApodWindow, m_Conf.apodizationArraySize); } } else if (m_Conf.DelayCalculationMethod == BeamformingSettings::DelayCalc::Spherical) { for (short line = 0; line < outputDim[0]; ++line) { threads[line] = std::thread(&BeamformingFilter::DASSphericalLine, this, m_InputData, m_OutputData, inputDim, outputDim, line, ApodWindow, m_Conf.apodizationArraySize); } } } else if (m_Conf.Algorithm == BeamformingSettings::BeamformingAlgorithm::DMAS) { if (m_Conf.DelayCalculationMethod == BeamformingSettings::DelayCalc::QuadApprox) { for (short line = 0; line < outputDim[0]; ++line) { threads[line] = std::thread(&BeamformingFilter::DMASQuadraticLine, this, m_InputData, m_OutputData, inputDim, outputDim, line, ApodWindow, m_Conf.apodizationArraySize); } } else if (m_Conf.DelayCalculationMethod == BeamformingSettings::DelayCalc::Spherical) { for (short line = 0; line < outputDim[0]; ++line) { threads[line] = std::thread(&BeamformingFilter::DMASSphericalLine, this, m_InputData, m_OutputData, inputDim, outputDim, line, ApodWindow, m_Conf.apodizationArraySize); } } } + else if (m_Conf.Algorithm == BeamformingSettings::BeamformingAlgorithm::sDMAS) + { + if (m_Conf.DelayCalculationMethod == BeamformingSettings::DelayCalc::QuadApprox) + { + for (short line = 0; line < outputDim[0]; ++line) + { + threads[line] = std::thread(&BeamformingFilter::sDMASQuadraticLine, this, m_InputData, m_OutputData, inputDim, outputDim, line, ApodWindow, m_Conf.apodizationArraySize); + } + } + else if (m_Conf.DelayCalculationMethod == BeamformingSettings::DelayCalc::Spherical) + { + for (short line = 0; line < outputDim[0]; ++line) + { + threads[line] = std::thread(&BeamformingFilter::sDMASSphericalLine, this, m_InputData, m_OutputData, inputDim, outputDim, line, ApodWindow, m_Conf.apodizationArraySize); + } + } + } // wait for all lines to finish for (short line = 0; line < outputDim[0]; ++line) { threads[line].join(); } output->SetSlice(m_OutputData, i); if (i % progInterval == 0) m_ProgressHandle((int)((i + 1) / (float)output->GetDimension(2) * 100), "performing reconstruction"); delete[] m_OutputData; m_OutputData = nullptr; m_InputData = nullptr; } } #if defined(PHOTOACOUSTICS_USE_GPU) || DOXYGEN else { try { // first, we check whether the data is float, other formats are unsupported if (!(input->GetPixelType().GetTypeAsString() == "scalar (float)" || input->GetPixelType().GetTypeAsString() == " (float)")) { MITK_ERROR << "Pixel type is not float, abort"; return; } m_ProgressHandle(50, "performing reconstruction"); m_BeamformingOclFilter->SetApodisation(ApodWindow, m_Conf.apodizationArraySize); m_BeamformingOclFilter->SetConfig(m_Conf); m_BeamformingOclFilter->SetInput(input); m_BeamformingOclFilter->Update(); void* out = m_BeamformingOclFilter->GetOutput(); output->SetImportVolume(out, 0, 0, mitk::Image::ImportMemoryManagementType::ManageMemory); } catch (mitk::Exception &e) { std::string errorMessage = "Caught unexpected exception "; errorMessage.append(e.what()); MITK_ERROR << errorMessage; float* dummyData = new float[m_Conf.ReconstructionLines * m_Conf.SamplesPerLine * m_Conf.inputDim[2]]; output->SetImportVolume(dummyData, 0, 0, mitk::Image::ImportMemoryManagementType::ManageMemory); m_Message = "An openCL error occurred; all GPU operations in this and the next session may be corrupted."; } } #endif m_TimeOfHeaderInitialization.Modified(); auto end = std::chrono::high_resolution_clock::now(); MITK_INFO << "Beamforming of " << output->GetDimension(2) << " Images completed in " << ((float)std::chrono::duration_cast(end - begin).count()) / 1000000 << "ms" << std::endl; } float* mitk::BeamformingFilter::VonHannFunction(int samples) { float* ApodWindow = new float[samples]; for (int n = 0; n < samples; ++n) { ApodWindow[n] = (1 - cos(2 * M_PI * n / (samples - 1))) / 2; } return ApodWindow; } float* mitk::BeamformingFilter::HammFunction(int samples) { float* ApodWindow = new float[samples]; for (int n = 0; n < samples; ++n) { ApodWindow[n] = 0.54 - 0.46*cos(2 * M_PI*n / (samples - 1)); } return ApodWindow; } float* mitk::BeamformingFilter::BoxFunction(int samples) { float* ApodWindow = new float[samples]; for (int n = 0; n < samples; ++n) { ApodWindow[n] = 1; } return ApodWindow; } void mitk::BeamformingFilter::DASQuadraticLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize) { float& inputS = inputDim[1]; float& inputL = inputDim[0]; float& outputS = outputDim[1]; float& outputL = outputDim[0]; short AddSample = 0; short maxLine = 0; short minLine = 0; float delayMultiplicator = 0; float l_i = 0; float s_i = 0; float part = 0.07 * inputL; float tan_phi = std::tan(m_Conf.Angle / 360 * 2 * M_PI); float part_multiplicator = tan_phi * m_Conf.TimeSpacing * m_Conf.SpeedOfSound / m_Conf.Pitch * inputL / m_Conf.TransducerElements; float apod_mult = 1; short usedLines = (maxLine - minLine); //quadratic delay l_i = line / outputL * inputL; for (short sample = 0; sample < outputS; ++sample) { s_i = (float)sample / outputS * inputS / 2; part = part_multiplicator*s_i; if (part < 1) part = 1; maxLine = (short)std::min((l_i + part) + 1, inputL); minLine = (short)std::max((l_i - part), 0.0f); usedLines = (maxLine - minLine); apod_mult = (float)apodArraySize / (float)usedLines; delayMultiplicator = pow((1 / (m_Conf.TimeSpacing*m_Conf.SpeedOfSound) * (m_Conf.Pitch*m_Conf.TransducerElements) / inputL), 2) / s_i / 2; for (short l_s = minLine; l_s < maxLine; ++l_s) { AddSample = delayMultiplicator * pow((l_s - l_i), 2) + s_i + (1 - m_Conf.isPhotoacousticImage)*s_i; if (AddSample < inputS && AddSample >= 0) output[sample*(short)outputL + line] += input[l_s + AddSample*(short)inputL] * apodisation[(short)((l_s - minLine)*apod_mult)]; else --usedLines; } output[sample*(short)outputL + line] = output[sample*(short)outputL + line] / usedLines; } } void mitk::BeamformingFilter::DASSphericalLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize) { float& inputS = inputDim[1]; float& inputL = inputDim[0]; float& outputS = outputDim[1]; float& outputL = outputDim[0]; short AddSample = 0; short maxLine = 0; short minLine = 0; float l_i = 0; float s_i = 0; float part = 0.07 * inputL; float tan_phi = std::tan(m_Conf.Angle / 360 * 2 * M_PI); float part_multiplicator = tan_phi * m_Conf.TimeSpacing * m_Conf.SpeedOfSound / m_Conf.Pitch * inputL / (float)m_Conf.TransducerElements; float apod_mult = 1; short usedLines = (maxLine - minLine); //exact delay l_i = (float)line / outputL * inputL; for (short sample = 0; sample < outputS; ++sample) { s_i = (float)sample / outputS * inputS / 2; part = part_multiplicator*s_i; if (part < 1) part = 1; maxLine = (short)std::min((l_i + part) + 1, inputL); minLine = (short)std::max((l_i - part), 0.0f); usedLines = (maxLine - minLine); apod_mult = (float)apodArraySize / (float)usedLines; for (short l_s = minLine; l_s < maxLine; ++l_s) { AddSample = (int)sqrt( pow(s_i, 2) + pow((1 / (m_Conf.TimeSpacing*m_Conf.SpeedOfSound) * (((float)l_s - l_i)*m_Conf.Pitch*(float)m_Conf.TransducerElements) / inputL), 2) ) + (1 - m_Conf.isPhotoacousticImage)*s_i; if (AddSample < inputS && AddSample >= 0) output[sample*(short)outputL + line] += input[l_s + AddSample*(short)inputL] * apodisation[(short)((l_s - minLine)*apod_mult)]; else --usedLines; } output[sample*(short)outputL + line] = output[sample*(short)outputL + line] / usedLines; } } + void mitk::BeamformingFilter::DMASQuadraticLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize) { float& inputS = inputDim[1]; float& inputL = inputDim[0]; float& outputS = outputDim[1]; float& outputL = outputDim[0]; short maxLine = 0; short minLine = 0; float delayMultiplicator = 0; float l_i = 0; float s_i = 0; float part = 0.07 * inputL; float tan_phi = std::tan(m_Conf.Angle / 360 * 2 * M_PI); float part_multiplicator = tan_phi * m_Conf.TimeSpacing * m_Conf.SpeedOfSound / m_Conf.Pitch * inputL / (float)m_Conf.TransducerElements; float apod_mult = 1; float mult = 0; short usedLines = (maxLine - minLine); //quadratic delay l_i = line / outputL * inputL; for (short sample = 0; sample < outputS; ++sample) { s_i = sample / outputS * inputS / 2; part = part_multiplicator*s_i; if (part < 1) part = 1; maxLine = (short)std::min((l_i + part) + 1, inputL); minLine = (short)std::max((l_i - part), 0.0f); usedLines = (maxLine - minLine); apod_mult = (float)apodArraySize / (float)usedLines; delayMultiplicator = pow((1 / (m_Conf.TimeSpacing*m_Conf.SpeedOfSound) * (m_Conf.Pitch*m_Conf.TransducerElements) / inputL), 2) / s_i / 2; //calculate the AddSamples beforehand to save some time short* AddSample = new short[maxLine - minLine]; for (short l_s = 0; l_s < maxLine - minLine; ++l_s) { AddSample[l_s] = (short)(delayMultiplicator * pow((minLine + l_s - l_i), 2) + s_i) + (1 - m_Conf.isPhotoacousticImage)*s_i; } float s_1 = 0; float s_2 = 0; - int sign = 0; for (short l_s1 = minLine; l_s1 < maxLine - 1; ++l_s1) { if (AddSample[l_s1 - minLine] < inputS && AddSample[l_s1 - minLine] >= 0) { for (short l_s2 = l_s1 + 1; l_s2 < maxLine; ++l_s2) { if (AddSample[l_s2 - minLine] < inputS && AddSample[l_s2 - minLine] >= 0) { s_2 = input[l_s2 + AddSample[l_s2 - minLine] * (short)inputL]; s_1 = input[l_s1 + AddSample[l_s1 - minLine] * (short)inputL]; - sign += ((s_1 > 0) - (s_1 < 0)); - sign += ((s_2 > 0) - (s_2 < 0)); - mult = s_2 * apodisation[(int)((l_s2 - minLine)*apod_mult)] * s_1 * apodisation[(int)((l_s1 - minLine)*apod_mult)]; output[sample*(short)outputL + line] += sqrt(fabs(mult)) * ((mult > 0) - (mult < 0)); } } } else --usedLines; } - output[sample*(short)outputL + line] = output[sample*(short)outputL + line] / (float)(pow(usedLines, 2) - (usedLines - 1)) * ((sign > 0) - (sign < 0)); + output[sample*(short)outputL + line] = output[sample*(short)outputL + line] / (float)(pow(usedLines, 2) - (usedLines - 1)); delete[] AddSample; } } void mitk::BeamformingFilter::DMASSphericalLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize) { float& inputS = inputDim[1]; float& inputL = inputDim[0]; float& outputS = outputDim[1]; float& outputL = outputDim[0]; short maxLine = 0; short minLine = 0; float l_i = 0; float s_i = 0; float part = 0.07 * inputL; float tan_phi = std::tan(m_Conf.Angle / 360 * 2 * M_PI); float part_multiplicator = tan_phi * m_Conf.TimeSpacing * m_Conf.SpeedOfSound / m_Conf.Pitch * inputL / (float)m_Conf.TransducerElements; float apod_mult = 1; float mult = 0; short usedLines = (maxLine - minLine); //exact delay l_i = (float)line / outputL * inputL; for (short sample = 0; sample < outputS; ++sample) { s_i = (float)sample / outputS * inputS / 2; part = part_multiplicator*s_i; if (part < 1) part = 1; maxLine = (short)std::min((l_i + part) + 1, inputL); minLine = (short)std::max((l_i - part), 0.0f); usedLines = (maxLine - minLine); apod_mult = (float)apodArraySize / (float)usedLines; //calculate the AddSamples beforehand to save some time short* AddSample = new short[maxLine - minLine]; for (short l_s = 0; l_s < maxLine - minLine; ++l_s) { AddSample[l_s] = (short)sqrt( pow(s_i, 2) + pow((1 / (m_Conf.TimeSpacing*m_Conf.SpeedOfSound) * (((float)minLine + (float)l_s - l_i)*m_Conf.Pitch*(float)m_Conf.TransducerElements) / inputL), 2) ) + (1 - m_Conf.isPhotoacousticImage)*s_i; } float s_1 = 0; float s_2 = 0; - int sign = 0; for (short l_s1 = minLine; l_s1 < maxLine - 1; ++l_s1) { if (AddSample[l_s1 - minLine] < inputS && AddSample[l_s1 - minLine] >= 0) { for (short l_s2 = l_s1 + 1; l_s2 < maxLine; ++l_s2) { if (AddSample[l_s2 - minLine] < inputS && AddSample[l_s2 - minLine] >= 0) { s_2 = input[l_s2 + AddSample[l_s2 - minLine] * (short)inputL]; s_1 = input[l_s1 + AddSample[l_s1 - minLine] * (short)inputL]; - sign += ((s_1 > 0) - (s_1 < 0)); - sign += ((s_2 > 0) - (s_2 < 0)); + mult = s_2 * apodisation[(int)((l_s2 - minLine)*apod_mult)] * s_1 * apodisation[(int)((l_s1 - minLine)*apod_mult)]; + output[sample*(short)outputL + line] += sqrt(fabs(mult)) * ((mult > 0) - (mult < 0)); + } + } + } + else + --usedLines; + } + + output[sample*(short)outputL + line] = output[sample*(short)outputL + line] / (float)(pow(usedLines, 2) - (usedLines - 1)); + + delete[] AddSample; + } +} + +void mitk::BeamformingFilter::sDMASQuadraticLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize) +{ + float& inputS = inputDim[1]; + float& inputL = inputDim[0]; + + float& outputS = outputDim[1]; + float& outputL = outputDim[0]; + + short maxLine = 0; + short minLine = 0; + float delayMultiplicator = 0; + float l_i = 0; + float s_i = 0; + + float part = 0.07 * inputL; + float tan_phi = std::tan(m_Conf.Angle / 360 * 2 * M_PI); + float part_multiplicator = tan_phi * m_Conf.TimeSpacing * m_Conf.SpeedOfSound / m_Conf.Pitch * inputL / (float)m_Conf.TransducerElements; + float apod_mult = 1; + + float mult = 0; + short usedLines = (maxLine - minLine); + + //quadratic delay + l_i = line / outputL * inputL; + + for (short sample = 0; sample < outputS; ++sample) + { + s_i = sample / outputS * inputS / 2; + + part = part_multiplicator*s_i; + + if (part < 1) + part = 1; + + maxLine = (short)std::min((l_i + part) + 1, inputL); + minLine = (short)std::max((l_i - part), 0.0f); + usedLines = (maxLine - minLine); + + apod_mult = (float)apodArraySize / (float)usedLines; + + delayMultiplicator = pow((1 / (m_Conf.TimeSpacing*m_Conf.SpeedOfSound) * (m_Conf.Pitch*m_Conf.TransducerElements) / inputL), 2) / s_i / 2; + + //calculate the AddSamples beforehand to save some time + short* AddSample = new short[maxLine - minLine]; + for (short l_s = 0; l_s < maxLine - minLine; ++l_s) + { + AddSample[l_s] = (short)(delayMultiplicator * pow((minLine + l_s - l_i), 2) + s_i) + (1 - m_Conf.isPhotoacousticImage)*s_i; + } + + float s_1 = 0; + float s_2 = 0; + float sign = 0; + + for (short l_s1 = minLine; l_s1 < maxLine - 1; ++l_s1) + { + if (AddSample[l_s1 - minLine] < inputS && AddSample[l_s1 - minLine] >= 0) + { + s_1 = input[l_s1 + AddSample[l_s1 - minLine] * (short)inputL]; + sign += s_1; + + for (short l_s2 = l_s1 + 1; l_s2 < maxLine; ++l_s2) + { + if (AddSample[l_s2 - minLine] < inputS && AddSample[l_s2 - minLine] >= 0) + { + s_2 = input[l_s2 + AddSample[l_s2 - minLine] * (short)inputL]; + + mult = s_2 * apodisation[(int)((l_s2 - minLine)*apod_mult)] * s_1 * apodisation[(int)((l_s1 - minLine)*apod_mult)]; + output[sample*(short)outputL + line] += sqrt(fabs(mult)) * ((mult > 0) - (mult < 0)); + } + } + } + else + --usedLines; + } + + output[sample*(short)outputL + line] = output[sample*(short)outputL + line] / (float)(pow(usedLines, 2) - (usedLines - 1)) * ((sign > 0) - (sign < 0)); + + delete[] AddSample; + } +} + +void mitk::BeamformingFilter::sDMASSphericalLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize) +{ + float& inputS = inputDim[1]; + float& inputL = inputDim[0]; + + float& outputS = outputDim[1]; + float& outputL = outputDim[0]; + + short maxLine = 0; + short minLine = 0; + float l_i = 0; + float s_i = 0; + + float part = 0.07 * inputL; + float tan_phi = std::tan(m_Conf.Angle / 360 * 2 * M_PI); + float part_multiplicator = tan_phi * m_Conf.TimeSpacing * m_Conf.SpeedOfSound / m_Conf.Pitch * inputL / (float)m_Conf.TransducerElements; + float apod_mult = 1; + + float mult = 0; + + short usedLines = (maxLine - minLine); + + //exact delay + + l_i = (float)line / outputL * inputL; + + for (short sample = 0; sample < outputS; ++sample) + { + s_i = (float)sample / outputS * inputS / 2; + + part = part_multiplicator*s_i; + + if (part < 1) + part = 1; + + maxLine = (short)std::min((l_i + part) + 1, inputL); + minLine = (short)std::max((l_i - part), 0.0f); + usedLines = (maxLine - minLine); + + apod_mult = (float)apodArraySize / (float)usedLines; + + //calculate the AddSamples beforehand to save some time + short* AddSample = new short[maxLine - minLine]; + for (short l_s = 0; l_s < maxLine - minLine; ++l_s) + { + AddSample[l_s] = (short)sqrt( + pow(s_i, 2) + + + pow((1 / (m_Conf.TimeSpacing*m_Conf.SpeedOfSound) * (((float)minLine + (float)l_s - l_i)*m_Conf.Pitch*(float)m_Conf.TransducerElements) / inputL), 2) + ) + (1 - m_Conf.isPhotoacousticImage)*s_i; + } + + float s_1 = 0; + float s_2 = 0; + float sign = 0; + + for (short l_s1 = minLine; l_s1 < maxLine - 1; ++l_s1) + { + if (AddSample[l_s1 - minLine] < inputS && AddSample[l_s1 - minLine] >= 0) + { + s_1 = input[l_s1 + AddSample[l_s1 - minLine] * (short)inputL]; + sign += s_1; + + for (short l_s2 = l_s1 + 1; l_s2 < maxLine; ++l_s2) + { + if (AddSample[l_s2 - minLine] < inputS && AddSample[l_s2 - minLine] >= 0) + { + s_2 = input[l_s2 + AddSample[l_s2 - minLine] * (short)inputL]; mult = s_2 * apodisation[(int)((l_s2 - minLine)*apod_mult)] * s_1 * apodisation[(int)((l_s1 - minLine)*apod_mult)]; output[sample*(short)outputL + line] += sqrt(fabs(mult)) * ((mult > 0) - (mult < 0)); } } } else --usedLines; } output[sample*(short)outputL + line] = output[sample*(short)outputL + line] / (float)(pow(usedLines, 2) - (usedLines - 1)) * ((sign > 0) - (sign < 0)); delete[] AddSample; } } diff --git a/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessing.cpp b/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessing.cpp index 12cb510802..bb38f0df70 100644 --- a/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessing.cpp +++ b/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessing.cpp @@ -1,1139 +1,1160 @@ /*=================================================================== 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. ===================================================================*/ // Blueberry #include #include // Qmitk #include "PAImageProcessing.h" // Qt #include #include #include #include //mitk image #include #include "mitkPhotoacousticImage.h" #include "mitkPhotoacousticBeamformingFilter.h" //other #include #include #include const std::string PAImageProcessing::VIEW_ID = "org.mitk.views.paimageprocessing"; PAImageProcessing::PAImageProcessing() : m_ResampleSpacing(0), m_UseLogfilter(false), m_FilterBank(mitk::PhotoacousticImage::New()) { qRegisterMetaType(); qRegisterMetaType(); } void PAImageProcessing::SetFocus() { m_Controls.buttonApplyBModeFilter->setFocus(); } void PAImageProcessing::CreateQtPartControl(QWidget *parent) { // create GUI widgets from the Qt Designer's .ui file m_Controls.setupUi(parent); connect(m_Controls.buttonApplyBModeFilter, SIGNAL(clicked()), this, SLOT(StartBmodeThread())); connect(m_Controls.DoResampling, SIGNAL(clicked()), this, SLOT(UseResampling())); connect(m_Controls.Logfilter, SIGNAL(clicked()), this, SLOT(UseLogfilter())); connect(m_Controls.ResamplingValue, SIGNAL(valueChanged(double)), this, SLOT(SetResampling())); connect(m_Controls.buttonApplyBeamforming, SIGNAL(clicked()), this, SLOT(StartBeamformingThread())); connect(m_Controls.buttonApplyCropFilter, SIGNAL(clicked()), this, SLOT(StartCropThread())); connect(m_Controls.buttonApplyBandpass, SIGNAL(clicked()), this, SLOT(StartBandpassThread())); connect(m_Controls.UseImageSpacing, SIGNAL(clicked()), this, SLOT(UseImageSpacing())); connect(m_Controls.ScanDepth, SIGNAL(valueChanged(double)), this, SLOT(UpdateImageInfo())); connect(m_Controls.SpeedOfSound, SIGNAL(valueChanged(double)), this, SLOT(UpdateImageInfo())); connect(m_Controls.SpeedOfSound, SIGNAL(valueChanged(double)), this, SLOT(ChangedSOSBeamforming())); connect(m_Controls.BPSpeedOfSound, SIGNAL(valueChanged(double)), this, SLOT(ChangedSOSBandpass())); connect(m_Controls.Samples, SIGNAL(valueChanged(int)), this, SLOT(UpdateImageInfo())); connect(m_Controls.UseImageSpacing, SIGNAL(clicked()), this, SLOT(UpdateImageInfo())); connect(m_Controls.boundLow, SIGNAL(valueChanged(int)), this, SLOT(LowerSliceBoundChanged())); connect(m_Controls.boundHigh, SIGNAL(valueChanged(int)), this, SLOT(UpperSliceBoundChanged())); connect(m_Controls.Partial, SIGNAL(clicked()), this, SLOT(SliceBoundsEnabled())); connect(m_Controls.BatchProcessing, SIGNAL(clicked()), this, SLOT(BatchProcessing())); connect(m_Controls.StepBeamforming, SIGNAL(clicked()), this, SLOT(UpdateSaveBoxes())); connect(m_Controls.StepCropping, SIGNAL(clicked()), this, SLOT(UpdateSaveBoxes())); connect(m_Controls.StepBandpass, SIGNAL(clicked()), this, SLOT(UpdateSaveBoxes())); connect(m_Controls.StepBMode, SIGNAL(clicked()), this, SLOT(UpdateSaveBoxes())); UpdateSaveBoxes(); m_Controls.DoResampling->setChecked(false); m_Controls.ResamplingValue->setEnabled(false); m_Controls.progressBar->setMinimum(0); m_Controls.progressBar->setMaximum(100); m_Controls.progressBar->setVisible(false); m_Controls.UseImageSpacing->setToolTip("Image spacing of y-Axis must be in us, x-Axis in mm."); m_Controls.UseImageSpacing->setToolTipDuration(5000); m_Controls.ProgressInfo->setVisible(false); m_Controls.UseBP->hide(); #ifndef PHOTOACOUSTICS_USE_GPU m_Controls.UseGPUBf->setEnabled(false); m_Controls.UseGPUBf->setChecked(false); m_Controls.UseGPUBmode->setEnabled(false); m_Controls.UseGPUBmode->setChecked(false); #endif UseImageSpacing(); } void PAImageProcessing::ChangedSOSBandpass() { m_Controls.SpeedOfSound->setValue(m_Controls.BPSpeedOfSound->value()); } void PAImageProcessing::ChangedSOSBeamforming() { m_Controls.BPSpeedOfSound->setValue(m_Controls.SpeedOfSound->value()); } std::vector splitpath( const std::string& str , const std::set delimiters) { std::vector result; char const* pch = str.c_str(); char const* start = pch; for (; *pch; ++pch) { if (delimiters.find(*pch) != delimiters.end()) { if (start != pch) { std::string str(start, pch); result.push_back(str); } else { result.push_back(""); } start = pch + 1; } } result.push_back(start); return result; } void PAImageProcessing::UpdateSaveBoxes() { if (m_Controls.StepBeamforming->isChecked()) m_Controls.SaveBeamforming->setEnabled(true); else m_Controls.SaveBeamforming->setEnabled(false); if (m_Controls.StepCropping->isChecked()) m_Controls.SaveCropping->setEnabled(true); else m_Controls.SaveCropping->setEnabled(false); if (m_Controls.StepBandpass->isChecked()) m_Controls.SaveBandpass->setEnabled(true); else m_Controls.SaveBandpass->setEnabled(false); if (m_Controls.StepBMode->isChecked()) m_Controls.SaveBMode->setEnabled(true); else m_Controls.SaveBMode->setEnabled(false); } void PAImageProcessing::BatchProcessing() { QFileDialog LoadDialog(nullptr, "Select Files to be processed"); LoadDialog.setFileMode(QFileDialog::FileMode::ExistingFiles); LoadDialog.setNameFilter(tr("Images (*.nrrd)")); LoadDialog.setViewMode(QFileDialog::Detail); QStringList fileNames; if (LoadDialog.exec()) fileNames = LoadDialog.selectedFiles(); QString saveDir = QFileDialog::getExistingDirectory(nullptr, tr("Select Directory To Save To"), "", QFileDialog::ShowDirsOnly | QFileDialog::DontResolveSymlinks); DisableControls(); std::set delims{'/'}; bool doSteps[] = { m_Controls.StepBeamforming->isChecked(), m_Controls.StepCropping->isChecked() , m_Controls.StepBandpass->isChecked(), m_Controls.StepBMode->isChecked() }; bool saveSteps[] = { m_Controls.SaveBeamforming->isChecked(), m_Controls.SaveCropping->isChecked() , m_Controls.SaveBandpass->isChecked(), m_Controls.SaveBMode->isChecked() }; for (int fileNumber = 0; fileNumber < fileNames.size(); ++fileNumber) { m_Controls.progressBar->setValue(0); m_Controls.progressBar->setVisible(true); m_Controls.ProgressInfo->setVisible(true); m_Controls.ProgressInfo->setText("loading file"); QString filename = fileNames.at(fileNumber); auto split = splitpath(filename.toStdString(), delims); std::string imageName = split.at(split.size()-1); // remove ".nrrd" imageName = imageName.substr(0, imageName.size()-5); mitk::Image::Pointer image = mitk::IOUtil::LoadImage(filename.toStdString().c_str()); UpdateBFSettings(image); // Beamforming if (doSteps[0]) { std::function progressHandle = [this](int progress, std::string progressInfo) { this->UpdateProgress(progress, progressInfo); }; m_Controls.progressBar->setValue(100); std::string errorMessage = ""; image = m_FilterBank->ApplyBeamforming(image, BFconfig, errorMessage, progressHandle); if (saveSteps[0]) { std::string saveFileName = saveDir.toStdString() + "/" + imageName + " beamformed" + ".nrrd"; mitk::IOUtil::Save(image, saveFileName); } } // Cropping if (doSteps[1]) { m_Controls.ProgressInfo->setText("cropping image"); image = m_FilterBank->ApplyCropping(image, m_Controls.CutoffAbove->value(), m_Controls.CutoffBelow->value(), 0, 0, 0, image->GetDimension(2) - 1); if (saveSteps[1]) { std::string saveFileName = saveDir.toStdString() + "/" + imageName + " cropped" + ".nrrd"; mitk::IOUtil::Save(image, saveFileName); } } // Bandpass if (doSteps[2]) { m_Controls.ProgressInfo->setText("applying bandpass"); float recordTime = image->GetDimension(1)*image->GetGeometry()->GetSpacing()[1] / 1000 / m_Controls.BPSpeedOfSound->value(); // add a safeguard so the program does not chrash when applying a Bandpass that reaches out of the bounds of the image float maxFrequency = 1 / (recordTime / image->GetDimension(1)) * image->GetDimension(1) / 2 / 2 / 1000; float BPHighPass = 1000000 * m_Controls.BPhigh->value(); // [Hz] float BPLowPass = maxFrequency - 1000000 * m_Controls.BPlow->value(); // [Hz] if (BPLowPass > maxFrequency && m_Controls.UseBP->isChecked()) { QMessageBox Msgbox; Msgbox.setText("LowPass too low, disabled it."); Msgbox.exec(); BPLowPass = 0; } if (BPLowPass < 0 && m_Controls.UseBP->isChecked()) { QMessageBox Msgbox; Msgbox.setText("LowPass too high, disabled it."); Msgbox.exec(); BPLowPass = 0; } if (BPHighPass > maxFrequency && m_Controls.UseBP->isChecked()) { QMessageBox Msgbox; Msgbox.setText("HighPass too high, disabled it."); Msgbox.exec(); BPHighPass = 0; } if (BPHighPass > maxFrequency - BFconfig.BPLowPass) { QMessageBox Msgbox; Msgbox.setText("HighPass higher than LowPass, disabled both."); Msgbox.exec(); BPHighPass = 0; BPLowPass = 0; } image = m_FilterBank->BandpassFilter(image, recordTime, BPHighPass, BPLowPass, m_Controls.BPFalloff->value()); if (saveSteps[2]) { std::string saveFileName = saveDir.toStdString() + "/" + imageName + " bandpassed" + ".nrrd"; mitk::IOUtil::Save(image, saveFileName); } } // Bmode if (doSteps[3]) { m_Controls.ProgressInfo->setText("applying bmode filter"); bool useGPU = m_Controls.UseGPUBmode->isChecked(); if (m_Controls.BModeMethod->currentText() == "Absolute Filter") image = m_FilterBank->ApplyBmodeFilter(image, mitk::PhotoacousticImage::BModeMethod::Abs, useGPU, m_UseLogfilter, m_ResampleSpacing); else if (m_Controls.BModeMethod->currentText() == "Envelope Detection") image = m_FilterBank->ApplyBmodeFilter(image, mitk::PhotoacousticImage::BModeMethod::EnvelopeDetection, useGPU, m_UseLogfilter, m_ResampleSpacing); if (saveSteps[3]) { std::string saveFileName = saveDir.toStdString() + "/" + imageName + " bmode" + ".nrrd"; mitk::IOUtil::Save(image, saveFileName); } } m_Controls.progressBar->setVisible(false); } EnableControls(); } void PAImageProcessing::StartBeamformingThread() { QList nodes = this->GetDataManagerSelection(); if (nodes.empty()) return; mitk::DataStorage::Pointer storage = this->GetDataStorage(); mitk::DataNode::Pointer node = nodes.front(); if (!node) { // Nothing selected. Inform the user and return QMessageBox::information(NULL, "Template", "Please load and select an image before starting image processing."); return; } mitk::BaseData* data = node->GetData(); if (data) { // test if this data item is an image or not (could also be a surface or something totally different) mitk::Image* image = dynamic_cast(data); if (image) { UpdateBFSettings(image); std::stringstream message; std::string name; message << "Performing beamforming for image "; if (node->GetName(name)) { // a property called "name" was found for this DataNode message << "'" << name << "'"; m_OldNodeName = name; } else m_OldNodeName = " "; message << "."; MITK_INFO << message.str(); m_Controls.progressBar->setValue(0); m_Controls.progressBar->setVisible(true); m_Controls.ProgressInfo->setVisible(true); m_Controls.ProgressInfo->setText("started"); m_Controls.buttonApplyBeamforming->setText("working..."); DisableControls(); BeamformingThread *thread = new BeamformingThread(); connect(thread, &BeamformingThread::result, this, &PAImageProcessing::HandleBeamformingResults); connect(thread, &BeamformingThread::updateProgress, this, &PAImageProcessing::UpdateProgress); connect(thread, &BeamformingThread::message, this, &PAImageProcessing::PAMessageBox); connect(thread, &BeamformingThread::finished, thread, &QObject::deleteLater); thread->setConfig(BFconfig); thread->setInputImage(image); thread->setFilterBank(m_FilterBank); MITK_INFO << "Started new thread for Beamforming"; thread->start(); } } } void PAImageProcessing::HandleBeamformingResults(mitk::Image::Pointer image) { auto newNode = mitk::DataNode::New(); newNode->SetData(image); // name the new Data node std::stringstream newNodeName; newNodeName << m_OldNodeName << " "; if (BFconfig.Algorithm == mitk::BeamformingSettings::BeamformingAlgorithm::DAS) newNodeName << "DAS bf, "; else if (BFconfig.Algorithm == mitk::BeamformingSettings::BeamformingAlgorithm::DMAS) newNodeName << "DMAS bf, "; if (BFconfig.DelayCalculationMethod == mitk::BeamformingSettings::DelayCalc::QuadApprox) newNodeName << "q. delay"; if (BFconfig.DelayCalculationMethod == mitk::BeamformingSettings::DelayCalc::Spherical) newNodeName << "s. delay"; newNode->SetName(newNodeName.str()); // update level window for the current dynamic range mitk::LevelWindow levelWindow; newNode->GetLevelWindow(levelWindow); levelWindow.SetAuto(image, true, true); newNode->SetLevelWindow(levelWindow); // add new node to data storage this->GetDataStorage()->Add(newNode); // disable progress bar m_Controls.progressBar->setVisible(false); m_Controls.ProgressInfo->setVisible(false); m_Controls.buttonApplyBeamforming->setText("Apply Beamforming"); EnableControls(); // update rendering mitk::RenderingManager::GetInstance()->InitializeViews(image->GetGeometry(), mitk::RenderingManager::REQUEST_UPDATE_ALL, true); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } void PAImageProcessing::StartBmodeThread() { QList nodes = this->GetDataManagerSelection(); if (nodes.empty()) return; mitk::DataStorage::Pointer storage = this->GetDataStorage(); mitk::DataNode::Pointer node = nodes.front(); if (!node) { // Nothing selected. Inform the user and return QMessageBox::information(NULL, "Template", "Please load and select an image before starting image processing."); return; } mitk::BaseData* data = node->GetData(); if (data) { // test if this data item is an image or not (could also be a surface or something totally different) mitk::Image* image = dynamic_cast(data); if (image) { UpdateBFSettings(image); std::stringstream message; std::string name; message << "Performing image processing for image "; if (node->GetName(name)) { // a property called "name" was found for this DataNode message << "'" << name << "'"; m_OldNodeName = name; } else m_OldNodeName = " "; message << "."; MITK_INFO << message.str(); m_Controls.buttonApplyBModeFilter->setText("working..."); DisableControls(); BmodeThread *thread = new BmodeThread(); connect(thread, &BmodeThread::result, this, &PAImageProcessing::HandleBmodeResults); connect(thread, &BmodeThread::finished, thread, &QObject::deleteLater); bool useGPU = m_Controls.UseGPUBmode->isChecked(); if(m_Controls.BModeMethod->currentText() == "Absolute Filter") thread->setConfig(m_UseLogfilter, m_ResampleSpacing, mitk::PhotoacousticImage::BModeMethod::Abs, useGPU); else if(m_Controls.BModeMethod->currentText() == "Envelope Detection") thread->setConfig(m_UseLogfilter, m_ResampleSpacing, mitk::PhotoacousticImage::BModeMethod::EnvelopeDetection, useGPU); thread->setInputImage(image); thread->setFilterBank(m_FilterBank); MITK_INFO << "Started new thread for Image Processing"; thread->start(); } } } void PAImageProcessing::HandleBmodeResults(mitk::Image::Pointer image) { auto newNode = mitk::DataNode::New(); newNode->SetData(image); // name the new Data node std::stringstream newNodeName; newNodeName << m_OldNodeName << " "; newNodeName << "B-Mode"; newNode->SetName(newNodeName.str()); // update level window for the current dynamic range mitk::LevelWindow levelWindow; newNode->GetLevelWindow(levelWindow); auto data = newNode->GetData(); levelWindow.SetAuto(dynamic_cast(data), true, true); newNode->SetLevelWindow(levelWindow); // add new node to data storage this->GetDataStorage()->Add(newNode); // disable progress bar m_Controls.progressBar->setVisible(false); m_Controls.buttonApplyBModeFilter->setText("Apply B-mode Filter"); EnableControls(); // update rendering mitk::RenderingManager::GetInstance()->InitializeViews( dynamic_cast(data)->GetGeometry(), mitk::RenderingManager::REQUEST_UPDATE_ALL, true); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } void PAImageProcessing::StartCropThread() { QList nodes = this->GetDataManagerSelection(); if (nodes.empty()) return; mitk::DataStorage::Pointer storage = this->GetDataStorage(); mitk::DataNode::Pointer node = nodes.front(); if (!node) { // Nothing selected. Inform the user and return QMessageBox::information(NULL, "Template", "Please load and select an image before starting image cropping."); return; } mitk::BaseData* data = node->GetData(); if (data) { // test if this data item is an image or not (could also be a surface or something totally different) mitk::Image* image = dynamic_cast(data); if (image) { UpdateBFSettings(image); std::stringstream message; std::string name; message << "Performing image cropping for image "; if (node->GetName(name)) { // a property called "name" was found for this DataNode message << "'" << name << "'"; m_OldNodeName = name; } else m_OldNodeName = " "; message << "."; MITK_INFO << message.str(); m_Controls.buttonApplyCropFilter->setText("working..."); DisableControls(); CropThread *thread = new CropThread(); connect(thread, &CropThread::result, this, &PAImageProcessing::HandleCropResults); connect(thread, &CropThread::finished, thread, &QObject::deleteLater); - thread->setConfig(m_Controls.CutoffAbove->value(), m_Controls.CutoffBelow->value(), m_Controls.CutoffFirstSlice->value(), m_Controls.CutoffLastSlice->value()); + thread->setConfig(m_Controls.CutoffAbove->value(), m_Controls.CutoffBelow->value(), 0, image->GetDimension(2) - 1); thread->setInputImage(image); thread->setFilterBank(m_FilterBank); MITK_INFO << "Started new thread for Image Cropping"; thread->start(); } } } void PAImageProcessing::HandleCropResults(mitk::Image::Pointer image) { auto newNode = mitk::DataNode::New(); newNode->SetData(image); // name the new Data node std::stringstream newNodeName; newNodeName << m_OldNodeName << " "; newNodeName << "Cropped"; newNode->SetName(newNodeName.str()); // update level window for the current dynamic range mitk::LevelWindow levelWindow; newNode->GetLevelWindow(levelWindow); auto data = newNode->GetData(); levelWindow.SetAuto(dynamic_cast(data), true, true); newNode->SetLevelWindow(levelWindow); // add new node to data storage this->GetDataStorage()->Add(newNode); m_Controls.buttonApplyCropFilter->setText("Apply Crop Filter"); EnableControls(); // update rendering mitk::RenderingManager::GetInstance()->InitializeViews( dynamic_cast(data)->GetGeometry(), mitk::RenderingManager::REQUEST_UPDATE_ALL, true); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } void PAImageProcessing::StartBandpassThread() { QList nodes = this->GetDataManagerSelection(); if (nodes.empty()) return; mitk::DataStorage::Pointer storage = this->GetDataStorage(); mitk::DataNode::Pointer node = nodes.front(); if (!node) { // Nothing selected. Inform the user and return QMessageBox::information(NULL, "Template", "Please load and select an image before starting image cropping."); return; } mitk::BaseData* data = node->GetData(); if (data) { // test if this data item is an image or not (could also be a surface or something totally different) mitk::Image* image = dynamic_cast(data); if (image) { UpdateBFSettings(image); std::stringstream message; std::string name; message << "Performing Bandpass filter on image "; if (node->GetName(name)) { // a property called "name" was found for this DataNode message << "'" << name << "'"; m_OldNodeName = name; } else m_OldNodeName = " "; message << "."; MITK_INFO << message.str(); m_Controls.buttonApplyBandpass->setText("working..."); DisableControls(); BandpassThread *thread = new BandpassThread(); connect(thread, &BandpassThread::result, this, &PAImageProcessing::HandleBandpassResults); connect(thread, &BandpassThread::finished, thread, &QObject::deleteLater); float recordTime = image->GetDimension(1)*image->GetGeometry()->GetSpacing()[1] / 1000 / m_Controls.BPSpeedOfSound->value(); // add a safeguard so the program does not chrash when applying a Bandpass that reaches out of the bounds of the image float maxFrequency = 1 / (recordTime / image->GetDimension(1)) * image->GetDimension(1) / 2 / 2 / 1000; float BPHighPass = 1000000 * m_Controls.BPhigh->value(); // [Hz] float BPLowPass = maxFrequency - 1000000 * m_Controls.BPlow->value(); // [Hz] if (BPLowPass > maxFrequency && m_Controls.UseBP->isChecked()) { QMessageBox Msgbox; Msgbox.setText("LowPass too low, disabled it."); Msgbox.exec(); BPLowPass = 0; } if (BPLowPass < 0 && m_Controls.UseBP->isChecked()) { QMessageBox Msgbox; Msgbox.setText("LowPass too high, disabled it."); Msgbox.exec(); BPLowPass = 0; } if (BPHighPass > maxFrequency && m_Controls.UseBP->isChecked()) { QMessageBox Msgbox; Msgbox.setText("HighPass too high, disabled it."); Msgbox.exec(); BPHighPass = 0; } if (BPHighPass > maxFrequency - BFconfig.BPLowPass) { QMessageBox Msgbox; Msgbox.setText("HighPass higher than LowPass, disabled both."); Msgbox.exec(); BPHighPass = 0; BPLowPass = 0; } thread->setConfig(BPHighPass, BPLowPass, m_Controls.BPFalloff->value(), recordTime); thread->setInputImage(image); thread->setFilterBank(m_FilterBank); MITK_INFO << "Started new thread for Bandpass filter"; thread->start(); } } } void PAImageProcessing::HandleBandpassResults(mitk::Image::Pointer image) { auto newNode = mitk::DataNode::New(); newNode->SetData(image); // name the new Data node std::stringstream newNodeName; newNodeName << m_OldNodeName << " "; newNodeName << "Bandpassed"; newNode->SetName(newNodeName.str()); // update level window for the current dynamic range mitk::LevelWindow levelWindow; newNode->GetLevelWindow(levelWindow); auto data = newNode->GetData(); levelWindow.SetAuto(dynamic_cast(data), true, true); newNode->SetLevelWindow(levelWindow); // add new node to data storage this->GetDataStorage()->Add(newNode); m_Controls.buttonApplyBandpass->setText("Apply Bandpass"); EnableControls(); // update rendering mitk::RenderingManager::GetInstance()->InitializeViews( dynamic_cast(data)->GetGeometry(), mitk::RenderingManager::REQUEST_UPDATE_ALL, true); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } void PAImageProcessing::SliceBoundsEnabled() { if (!m_Controls.Partial->isChecked()) { m_Controls.boundLow->setEnabled(false); m_Controls.boundHigh->setEnabled(false); return; } else { m_Controls.boundLow->setEnabled(true); m_Controls.boundHigh->setEnabled(true); } } void PAImageProcessing::UpperSliceBoundChanged() { if(m_Controls.boundLow->value() > m_Controls.boundHigh->value()) { m_Controls.boundLow->setValue(m_Controls.boundHigh->value()); } } void PAImageProcessing::LowerSliceBoundChanged() { if (m_Controls.boundLow->value() > m_Controls.boundHigh->value()) { m_Controls.boundHigh->setValue(m_Controls.boundLow->value()); } } void PAImageProcessing::UpdateProgress(int progress, std::string progressInfo) { if (progress < 100) m_Controls.progressBar->setValue(progress); else m_Controls.progressBar->setValue(100); m_Controls.ProgressInfo->setText(progressInfo.c_str()); qApp->processEvents(); } void PAImageProcessing::PAMessageBox(std::string message) { if (0 != message.compare("noMessage")) { QMessageBox msgBox; msgBox.setText(message.c_str()); msgBox.exec(); } } void PAImageProcessing::UpdateImageInfo() { QList nodes = this->GetDataManagerSelection(); if (nodes.empty()) return; mitk::DataNode::Pointer node = nodes.front(); if (!node) { // Nothing selected return; } mitk::BaseData* data = node->GetData(); if (data) { // test if this data item is an image or not (could also be a surface or something totally different) mitk::Image* image = dynamic_cast(data); if (image) { // beamforming configs if (m_Controls.UseImageSpacing->isChecked()) { m_Controls.ElementCount->setValue(image->GetDimension(0)); m_Controls.Pitch->setValue(image->GetGeometry()->GetSpacing()[0]); - m_Controls.boundLow->setMaximum(image->GetDimension(2) - 1); - m_Controls.boundHigh->setMaximum(image->GetDimension(2) - 1); - m_Controls.CutoffFirstSlice->setValue(0); - m_Controls.CutoffLastSlice->setValue(image->GetDimension(2) - 1); - m_Controls.CutoffLastSlice->setMaximum(image->GetDimension(2) - 1); } + + m_Controls.boundLow->setMaximum(image->GetDimension(2) - 1); + m_Controls.boundHigh->setMaximum(image->GetDimension(2) - 1); + UpdateBFSettings(image); m_Controls.CutoffBeforeBF->setValue(0.000001 / BFconfig.TimeSpacing); // 1us standard offset for our transducer std::stringstream frequency; float maxFrequency = (1 / BFconfig.TimeSpacing) * image->GetDimension(1) / 2 / 2 / 1000; frequency << maxFrequency / 1000000; //[MHz] frequency << "MHz"; m_Controls.BPhigh->setMaximum(maxFrequency / 1000000); m_Controls.BPlow->setMaximum(maxFrequency / 1000000); frequency << " is the maximal allowed frequency for the selected image."; m_Controls.BPhigh->setToolTip(frequency.str().c_str()); m_Controls.BPlow->setToolTip(frequency.str().c_str()); m_Controls.BPhigh->setToolTipDuration(5000); m_Controls.BPlow->setToolTipDuration(5000); } } } void PAImageProcessing::OnSelectionChanged( berry::IWorkbenchPart::Pointer /*source*/, const QList& nodes ) { // iterate all selected objects, adjust warning visibility foreach( mitk::DataNode::Pointer node, nodes ) { if( node.IsNotNull() && dynamic_cast(node->GetData()) ) { m_Controls.labelWarning->setVisible( false ); m_Controls.buttonApplyBModeFilter->setEnabled( true ); m_Controls.labelWarning2->setVisible(false); m_Controls.buttonApplyCropFilter->setEnabled(true); m_Controls.labelWarning3->setVisible(false); m_Controls.buttonApplyBandpass->setEnabled(true); m_Controls.labelWarning4->setVisible(false); m_Controls.buttonApplyBeamforming->setEnabled(true); UpdateImageInfo(); return; } } m_Controls.labelWarning->setVisible( true ); m_Controls.buttonApplyBModeFilter->setEnabled( false ); m_Controls.labelWarning2->setVisible(true); m_Controls.buttonApplyCropFilter->setEnabled(false); m_Controls.labelWarning3->setVisible(true); m_Controls.buttonApplyBandpass->setEnabled(false); m_Controls.labelWarning4->setVisible(true); m_Controls.buttonApplyBeamforming->setEnabled(false); } void PAImageProcessing::UseResampling() { if (m_Controls.DoResampling->isChecked()) { m_Controls.ResamplingValue->setEnabled(true); m_ResampleSpacing = m_Controls.ResamplingValue->value(); } else { m_Controls.ResamplingValue->setEnabled(false); m_ResampleSpacing = 0; } } void PAImageProcessing::UseLogfilter() { m_UseLogfilter = m_Controls.Logfilter->isChecked(); } void PAImageProcessing::SetResampling() { m_ResampleSpacing = m_Controls.ResamplingValue->value(); } void PAImageProcessing::UpdateBFSettings(mitk::Image::Pointer image) { if ("DAS" == m_Controls.BFAlgorithm->currentText()) BFconfig.Algorithm = mitk::BeamformingSettings::BeamformingAlgorithm::DAS; else if ("DMAS" == m_Controls.BFAlgorithm->currentText()) BFconfig.Algorithm = mitk::BeamformingSettings::BeamformingAlgorithm::DMAS; + else if ("sDMAS" == m_Controls.BFAlgorithm->currentText()) + BFconfig.Algorithm = mitk::BeamformingSettings::BeamformingAlgorithm::sDMAS; if ("Quad. Approx." == m_Controls.DelayCalculation->currentText()) { BFconfig.DelayCalculationMethod = mitk::BeamformingSettings::DelayCalc::QuadApprox; } else if ("Spherical Wave" == m_Controls.DelayCalculation->currentText()) { BFconfig.DelayCalculationMethod = mitk::BeamformingSettings::DelayCalc::Spherical; } if ("Von Hann" == m_Controls.Apodization->currentText()) { BFconfig.Apod = mitk::BeamformingSettings::Apodization::Hann; } else if ("Hamming" == m_Controls.Apodization->currentText()) { BFconfig.Apod = mitk::BeamformingSettings::Apodization::Hamm; } else if ("Box" == m_Controls.Apodization->currentText()) { BFconfig.Apod = mitk::BeamformingSettings::Apodization::Box; } BFconfig.Pitch = m_Controls.Pitch->value() / 1000; // [m] BFconfig.SpeedOfSound = m_Controls.SpeedOfSound->value(); // [m/s] BFconfig.SamplesPerLine = m_Controls.Samples->value(); BFconfig.ReconstructionLines = m_Controls.Lines->value(); BFconfig.TransducerElements = m_Controls.ElementCount->value(); - BFconfig.apodizationArraySize = m_Controls.ElementCount->value(); + BFconfig.apodizationArraySize = m_Controls.Lines->value(); BFconfig.Angle = m_Controls.Angle->value(); // [deg] BFconfig.UseBP = m_Controls.UseBP->isChecked(); BFconfig.UseGPU = m_Controls.UseGPUBf->isChecked(); BFconfig.upperCutoff = m_Controls.CutoffBeforeBF->value(); if (m_Controls.UseImageSpacing->isChecked()) { BFconfig.RecordTime = image->GetDimension(1)*image->GetGeometry()->GetSpacing()[1] / 1000000; // [s] BFconfig.TimeSpacing = image->GetGeometry()->GetSpacing()[1] / 1000000; MITK_INFO << "Calculated Scan Depth of " << BFconfig.RecordTime * BFconfig.SpeedOfSound * 100 / 2 << "cm"; } else { BFconfig.RecordTime = 2 * m_Controls.ScanDepth->value() / 1000 / BFconfig.SpeedOfSound; // [s] BFconfig.TimeSpacing = BFconfig.RecordTime / image->GetDimension(1); } if ("US Image" == m_Controls.ImageType->currentText()) { BFconfig.isPhotoacousticImage = false; } else if ("PA Image" == m_Controls.ImageType->currentText()) { BFconfig.isPhotoacousticImage = true; } BFconfig.partial = m_Controls.Partial->isChecked(); BFconfig.CropBounds[0] = m_Controls.boundLow->value(); BFconfig.CropBounds[1] = m_Controls.boundHigh->value(); } void PAImageProcessing::EnableControls() { + m_Controls.BatchProcessing->setEnabled(true); + m_Controls.StepBeamforming->setEnabled(true); + m_Controls.StepBandpass->setEnabled(true); + m_Controls.StepCropping->setEnabled(true); + m_Controls.StepBMode->setEnabled(true); + + UpdateSaveBoxes(); + m_Controls.DoResampling->setEnabled(true); UseResampling(); m_Controls.Logfilter->setEnabled(true); + m_Controls.BModeMethod->setEnabled(true); m_Controls.buttonApplyBModeFilter->setEnabled(true); m_Controls.CutoffAbove->setEnabled(true); m_Controls.CutoffBelow->setEnabled(true); m_Controls.CutoffBeforeBF->setEnabled(true); m_Controls.buttonApplyCropFilter->setEnabled(true); m_Controls.buttonApplyBandpass->setEnabled(true); m_Controls.Partial->setEnabled(true); m_Controls.boundHigh->setEnabled(true); m_Controls.boundLow->setEnabled(true); m_Controls.BFAlgorithm->setEnabled(true); m_Controls.DelayCalculation->setEnabled(true); m_Controls.ImageType->setEnabled(true); m_Controls.Apodization->setEnabled(true); m_Controls.UseBP->setEnabled(true); #ifdef PHOTOACOUSTICS_USE_GPU m_Controls.UseGPUBf->setEnabled(true); m_Controls.UseGPUBmode->setEnabled(true); #endif m_Controls.BPhigh->setEnabled(true); m_Controls.BPlow->setEnabled(true); m_Controls.BPFalloff->setEnabled(true); m_Controls.UseImageSpacing->setEnabled(true); UseImageSpacing(); m_Controls.Pitch->setEnabled(true); m_Controls.ElementCount->setEnabled(true); m_Controls.SpeedOfSound->setEnabled(true); m_Controls.Samples->setEnabled(true); m_Controls.Lines->setEnabled(true); m_Controls.Angle->setEnabled(true); m_Controls.buttonApplyBeamforming->setEnabled(true); } void PAImageProcessing::DisableControls() { + m_Controls.BatchProcessing->setEnabled(false); + m_Controls.StepBeamforming->setEnabled(false); + m_Controls.StepBandpass->setEnabled(false); + m_Controls.StepCropping->setEnabled(false); + m_Controls.StepBMode->setEnabled(false); + m_Controls.SaveBeamforming->setEnabled(false); + m_Controls.SaveBandpass->setEnabled(false); + m_Controls.SaveCropping->setEnabled(false); + m_Controls.SaveBMode->setEnabled(false); + m_Controls.DoResampling->setEnabled(false); m_Controls.ResamplingValue->setEnabled(false); m_Controls.Logfilter->setEnabled(false); + m_Controls.BModeMethod->setEnabled(false); m_Controls.buttonApplyBModeFilter->setEnabled(false); m_Controls.CutoffAbove->setEnabled(false); m_Controls.CutoffBelow->setEnabled(false); m_Controls.CutoffBeforeBF->setEnabled(false); m_Controls.buttonApplyCropFilter->setEnabled(false); m_Controls.buttonApplyBandpass->setEnabled(false); m_Controls.Partial->setEnabled(false); m_Controls.boundHigh->setEnabled(false); m_Controls.boundLow->setEnabled(false); m_Controls.BFAlgorithm->setEnabled(false); m_Controls.DelayCalculation->setEnabled(false); m_Controls.ImageType->setEnabled(false); m_Controls.Apodization->setEnabled(false); m_Controls.UseBP->setEnabled(false); #ifdef PHOTOACOUSTICS_USE_GPU m_Controls.UseGPUBf->setEnabled(false); m_Controls.UseGPUBmode->setEnabled(false); #endif m_Controls.BPhigh->setEnabled(false); m_Controls.BPlow->setEnabled(false); m_Controls.BPFalloff->setEnabled(false); m_Controls.UseImageSpacing->setEnabled(false); m_Controls.ScanDepth->setEnabled(false); m_Controls.Pitch->setEnabled(false); m_Controls.ElementCount->setEnabled(false); m_Controls.SpeedOfSound->setEnabled(false); m_Controls.Samples->setEnabled(false); m_Controls.Lines->setEnabled(false); m_Controls.Angle->setEnabled(false); m_Controls.buttonApplyBeamforming->setEnabled(false); } void PAImageProcessing::UseImageSpacing() { if (m_Controls.UseImageSpacing->isChecked()) { m_Controls.ScanDepth->setDisabled(true); } else { m_Controls.ScanDepth->setEnabled(true); } } #include void BeamformingThread::run() { mitk::Image::Pointer resultImage = mitk::Image::New(); mitk::Image::Pointer resultImageBuffer; std::string errorMessage = ""; std::function progressHandle = [this](int progress, std::string progressInfo) { emit updateProgress(progress, progressInfo); }; resultImageBuffer = m_FilterBank->ApplyBeamforming(m_InputImage, m_BFconfig, errorMessage, progressHandle); mitk::ImageReadAccessor copy(resultImageBuffer); resultImage->Initialize(resultImageBuffer); resultImage->SetSpacing(resultImageBuffer->GetGeometry()->GetSpacing()); resultImage->SetImportVolume(const_cast(copy.GetData()), 0, 0, mitk::Image::CopyMemory); emit result(resultImage); emit message(errorMessage); } void BeamformingThread::setConfig(mitk::BeamformingSettings BFconfig) { m_BFconfig = BFconfig; } void BeamformingThread::setInputImage(mitk::Image::Pointer image) { m_InputImage = image; } void BmodeThread::run() { mitk::Image::Pointer resultImage; resultImage = m_FilterBank->ApplyBmodeFilter(m_InputImage, m_Method, m_UseGPU, m_UseLogfilter, m_ResampleSpacing); emit result(resultImage); } void BmodeThread::setConfig(bool useLogfilter, double resampleSpacing, mitk::PhotoacousticImage::BModeMethod method, bool useGPU) { m_UseLogfilter = useLogfilter; m_ResampleSpacing = resampleSpacing; m_Method = method; m_UseGPU = useGPU; } void BmodeThread::setInputImage(mitk::Image::Pointer image) { m_InputImage = image; } void CropThread::run() { mitk::Image::Pointer resultImage; resultImage = m_FilterBank->ApplyCropping(m_InputImage, m_CutAbove, m_CutBelow, 0, 0, m_CutSliceFirst, m_CutSliceLast); emit result(resultImage); } void CropThread::setConfig(unsigned int CutAbove, unsigned int CutBelow, unsigned int CutSliceFirst, unsigned int CutSliceLast) { m_CutAbove = CutAbove; m_CutBelow = CutBelow; m_CutSliceLast = CutSliceLast; m_CutSliceFirst = CutSliceFirst; } void CropThread::setInputImage(mitk::Image::Pointer image) { m_InputImage = image; } void BandpassThread::run() { mitk::Image::Pointer resultImage; resultImage = m_FilterBank->BandpassFilter(m_InputImage, m_RecordTime, m_BPHighPass, m_BPLowPass, m_TukeyAlpha); emit result(resultImage); } void BandpassThread::setConfig(float BPHighPass, float BPLowPass, float TukeyAlpha, float recordTime) { m_BPHighPass = BPHighPass; m_BPLowPass = BPLowPass; m_TukeyAlpha = TukeyAlpha; m_RecordTime = recordTime; } void BandpassThread::setInputImage(mitk::Image::Pointer image) { m_InputImage = image; } diff --git a/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessingControls.ui b/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessingControls.ui index 085fe4153b..7883dd2171 100644 --- a/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessingControls.ui +++ b/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessingControls.ui @@ -1,1006 +1,991 @@ PAImageProcessingControls 0 0 382 1278 0 0 QmitkTemplate <html><head/><body><p><span style=" font-weight:600;">Batch Processing</span></p></body></html> Start Batch Processing Bandpass true Crop true Save true Save true Save Beamform true BMode true Save true <html><head/><body><p><span style=" font-weight:600;">B-mode Filter Settings</span></p></body></html> Do Resampling true Envelope Detection Absolute Filter Envelope Detection 0 0 13 0 11 3 0.010000000000000 1.000000000000000 0.010000000000000 0.075000000000000 [mm] Depth Spacing Add Logfilter Use GPU QLabel { color: rgb(255, 0, 0) } <html><head/><body><p align="center"><span style=" font-size:10pt; font-weight:600;">Please select an image!</span></p></body></html> 0 0 Do image processing Apply B-mode Filter Qt::Horizontal <html><head/><body><p><span style=" font-weight:600;">Bandpass Filter Settings</span></p></body></html> QLayout::SetDefaultConstraint 0 0 0 3 0.010000000000000 200.000000000000000 15.000000000000000 [MHz] f High Pass [MHz] f Low Pass 0 0 3 200.000000000000000 Thukey window alpha 1 200.000000000000000 3000.000000000000000 5.000000000000000 1540.000000000000000 [m/s] Speed of Sound 2 1.000000000000000 0.100000000000000 0.500000000000000 <html><head/><body><p align="center"><span style=" font-size:10pt; font-weight:600; color:#ff0000;">Please select an image!</span></p></body></html> Apply Bandpass Qt::Horizontal <html><head/><body><p><span style=" font-weight:600;">Crop Filter Settings</span></p></body></html> 99999 10 Cut Top Cut Bottom 99999 5 165 - - - - - - - Cut Slices first - - - - - - - - - - Cut Slices last - - - <html><head/><body><p align="center"><span style=" font-size:10pt; font-weight:600; color:#ff0000;">Please select an image!</span></p></body></html> Apply Crop Filer Qt::Horizontal <html><head/><body><p><span style=" font-weight:600;">Beamforming Filter Settings</span></p></body></html> 5 2 Delay Calculation Auto Get Depth true Apply Beamforming Beamforming Method [mm] Scan Depth 0 0 3 0.010000000000000 9.000000000000000 0.050000000000000 0.300000000000000 Transducer Elements 0 0 4 300.000000000000000 0.100000000000000 50.000000000000000 [mm] Transducer Pitch 0 0 64 1024 128 128 0 0 256 16384 256 2048 0 0 64 2048 128 256 Samples Reconstruction Lines true 0 0 100 0 0 0 900 10 0 0 0 DAS DMAS + + + sDMAS + + 0 0 Spherical Wave Quad. Approx. Spherical Wave 0 0 PA Image US Image Image Type 0 0 Von Hann Hamming Box Apodization 0 0 1 200.000000000000000 3000.000000000000000 5.000000000000000 1540.000000000000000 [m/s] Speed of Sound false 99999 minimal beamformed slice min false 99999 10 Maximal beamformed slice max select slices Compute On GPU true true Auto Use Bandpass 0 0 1 1.000000000000000 180.000000000000000 27.000000000000000 [°] Element Angle Cutoff Upper Voxels <html><head/><body><p align="center"><span style=" font-size:10pt; font-weight:600; color:#ff0000;">Please select an image!</span></p></body></html>