diff --git a/Modules/OpenCL/mitkOclUtils.cpp b/Modules/OpenCL/mitkOclUtils.cpp index d743b3d494..a8fa4c0c53 100644 --- a/Modules/OpenCL/mitkOclUtils.cpp +++ b/Modules/OpenCL/mitkOclUtils.cpp @@ -1,572 +1,576 @@ /*=================================================================== 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 "mitkOclUtils.h" #include "mitkLogMacros.h" #include #include unsigned int iDivUp(unsigned int dividend, unsigned int divisor){ return (dividend % divisor == 0) ? (dividend / divisor) : (dividend / divisor + 1); } cl_int oclGetPlatformID(cl_platform_id* selectedPlatform) { cl_uint num_platforms = 0; cl_platform_id* clPlatformIDs; cl_int ciErrNum = 0; ciErrNum = clGetPlatformIDs( 0, nullptr, &num_platforms); if ( ciErrNum != CL_SUCCESS) { MITK_ERROR<<" Error " << ciErrNum << " in clGetPlatformIDs() \n"; throw std::bad_exception(); } else { clPlatformIDs = new cl_platform_id[num_platforms]; ciErrNum = clGetPlatformIDs( num_platforms, clPlatformIDs, nullptr); if(ciErrNum == CL_SUCCESS) { *selectedPlatform = clPlatformIDs[0]; } } return CL_SUCCESS; } void oclPrintMemObjectInfo(cl_mem memobj) { cl_int clErr = 0; MITK_INFO << "Examining cl_mem object: " << memobj << "\n------------------\n"; // CL_MEM_TYPE cl_mem_object_type objtype; clErr = clGetMemObjectInfo( memobj, CL_MEM_TYPE, sizeof(cl_mem_object_type),&objtype, nullptr); CHECK_OCL_ERR( clErr ); switch(objtype) { case CL_MEM_OBJECT_BUFFER: MITK_INFO << "CL_MEM_TYPE \t" << "BUFFER_OBJ" << "\n"; break; case CL_MEM_OBJECT_IMAGE2D: MITK_INFO << "CL_MEM_TYPE \t" << "2D IMAGE" << "\n"; break; case CL_MEM_OBJECT_IMAGE3D: MITK_INFO << "CL_MEM_TYPE \t" << "3D IMAGE" << "\n"; break; default: MITK_INFO << "CL_MEM_TYPE \t" << "[could not resolve]" << "\n"; break; } // CL_MEM_FLAGS cl_mem_flags flags; clErr = clGetMemObjectInfo( memobj, CL_MEM_FLAGS, sizeof(cl_mem_flags),&flags, nullptr); CHECK_OCL_ERR( clErr ); switch(flags) { case CL_MEM_READ_ONLY: MITK_INFO << "CL_MEM_FLAGS \t" << "CL_MEM_READ_ONLY" << "\n"; break; case CL_MEM_WRITE_ONLY: MITK_INFO << "CL_MEM_FLAGS \t" << "CL_MEM_WRITE_ONLY" << "\n"; break; case CL_MEM_READ_WRITE: MITK_INFO << "CL_MEM_FLAGS \t" << "CL_MEM_READ_WRITE" << "\n"; break; default: MITK_INFO << "CL_MEM_FLAGS \t" << "not resolved, " << flags << "\n"; break; } // get CL_MEM_SIZE size_t memsize; clErr = clGetMemObjectInfo( memobj, CL_MEM_SIZE, sizeof(memsize),&memsize, nullptr); CHECK_OCL_ERR( clErr ); MITK_INFO << "CL_MEM_SIZE \t" << memsize << "\n"; // get CL_MEM_HOST_PTR float *hostptr; clErr = clGetMemObjectInfo( memobj, CL_MEM_HOST_PTR, sizeof(void*), (void*) &hostptr, nullptr); CHECK_OCL_ERR( clErr ); MITK_INFO << "CL_MEM_HOST_PTR \t" << hostptr << "\n"; // get CL_CONTEXT cl_context gpuctxt; clErr = clGetMemObjectInfo( memobj, CL_MEM_CONTEXT, sizeof(cl_context), &gpuctxt, nullptr); CHECK_OCL_ERR( clErr ); MITK_INFO << "CL_CONTEXT \t\t" << gpuctxt << "\n"; // get CL_MEM_REFERENCE_COUNT cl_uint refs; clErr = clGetMemObjectInfo( memobj, CL_MEM_REFERENCE_COUNT, sizeof(cl_uint), &refs, nullptr); CHECK_OCL_ERR(clErr); MITK_INFO << "CL_REF_COUNT \t" << refs << "\n"; MITK_INFO << "================== \n" << std::endl; } void oclPrintDeviceInfo(cl_device_id device) { char device_string[1024]; clGetDeviceInfo(device, CL_DEVICE_NAME, sizeof(device_string), &device_string, nullptr); MITK_INFO("ocl.log")<< " Device : " << device_string; // CL_DEVICE_INFO cl_device_type type; clGetDeviceInfo(device, CL_DEVICE_TYPE, sizeof(type), &type, nullptr); if( type & CL_DEVICE_TYPE_CPU ) MITK_INFO("ocl.log")<<" CL_DEVICE_TYPE: CL_DEVICE_TYPE_CPU"; if( type & CL_DEVICE_TYPE_GPU ) MITK_INFO("ocl.log")<<" CL_DEVICE_TYPE: CL_DEVICE_TYPE_GPU"; if( type & CL_DEVICE_TYPE_ACCELERATOR ) MITK_INFO("ocl.log")<<" CL_DEVICE_TYPE: CL_DEVICE_TYPE_ACCELERATOR"; if( type & CL_DEVICE_TYPE_DEFAULT ) MITK_INFO("ocl.log")<<" CL_DEVICE_TYPE: CL_DEVICE_TYPE_DEFAULT"; // CL_DEVICE_MAX_COMPUTE_UNITS cl_uint compute_units; clGetDeviceInfo(device, CL_DEVICE_MAX_COMPUTE_UNITS, sizeof(compute_units), &compute_units, nullptr); MITK_INFO("ocl.log")<<" CL_DEVICE_MAX_COMPUTE_UNITS:" << compute_units; // CL_DEVICE_MAX_WORK_GROUP_SIZE size_t workitem_size[3]; clGetDeviceInfo(device, CL_DEVICE_MAX_WORK_ITEM_SIZES, sizeof(workitem_size), &workitem_size, nullptr); MITK_INFO("ocl.log")<<" CL_DEVICE_MAX_WORK_ITEM_SIZES:\t"<< workitem_size[0]<< workitem_size[1]<< workitem_size[2]; // CL_DEVICE_MAX_WORK_GROUP_SIZE size_t workgroup_size; clGetDeviceInfo(device, CL_DEVICE_MAX_WORK_GROUP_SIZE, sizeof(workgroup_size), &workgroup_size, nullptr); MITK_INFO("ocl.log")<<" CL_DEVICE_MAX_WORK_GROUP_SIZE:" << workgroup_size; // CL_DEVICE_MAX_CLOCK_FREQUENCY cl_uint clock_frequency; clGetDeviceInfo(device, CL_DEVICE_MAX_CLOCK_FREQUENCY, sizeof(clock_frequency), &clock_frequency, nullptr); MITK_INFO("ocl.log")<<" CL_DEVICE_MAX_CLOCK_FREQUENCY:"<< clock_frequency / 1000; // CL_DEVICE_IMAGE_SUPPORT cl_bool image_support; clGetDeviceInfo(device, CL_DEVICE_IMAGE_SUPPORT, sizeof(image_support), &image_support, nullptr); MITK_INFO("ocl.log")<<" CL_DEVICE_IMAGE_SUPPORT:\t" << image_support; // CL_DEVICE_GLOBAL_MEM_SIZE cl_ulong mem_size; clGetDeviceInfo(device, CL_DEVICE_GLOBAL_MEM_SIZE, sizeof(mem_size), &mem_size, nullptr); MITK_INFO("ocl.log")<<" CL_DEVICE_GLOBAL_MEM_SIZE:\t\t"<<(unsigned int)(mem_size / (1024 * 1024))<<"Mbytes"; // CL_DEVICE_LOCAL_MEM_SIZE clGetDeviceInfo(device, CL_DEVICE_LOCAL_MEM_SIZE, sizeof(mem_size), &mem_size, nullptr); MITK_INFO("ocl.log")<<" CL_DEVICE_LOCAL_MEM_SIZE:\t\t"<< (unsigned int)(mem_size / (1024)) <<"KByte\n"; + // CL_DEVICE_MAX_CONSTANT_BUFFER_SIZE + clGetDeviceInfo(device, CL_DEVICE_MAX_CONSTANT_BUFFER_SIZE, sizeof(mem_size), &mem_size, nullptr); + MITK_INFO("ocl.log") << " CL_DEVICE_MAX_CONSTANT_BUFFER_SIZE:\t\t" << (unsigned int)(mem_size / (1024)) << "KByte"; + //check for image support properties clGetDeviceInfo(device, CL_DEVICE_IMAGE2D_MAX_WIDTH, sizeof(workgroup_size), &workgroup_size, nullptr); MITK_INFO("ocl.log")<<" CL_DEVICE_IMAGE2D_MAX_WIDTH:\t" << workgroup_size; clGetDeviceInfo(device, CL_DEVICE_IMAGE2D_MAX_HEIGHT, sizeof(workgroup_size), &workgroup_size, nullptr); MITK_INFO("ocl.log")<<" CL_DEVICE_IMAGE2D_MAX_HEIGHT:\t" << workgroup_size; clGetDeviceInfo(device, CL_DEVICE_IMAGE3D_MAX_WIDTH, sizeof(workgroup_size), &workgroup_size, nullptr); MITK_INFO("ocl.log")<<" CL_DEVICE_IMAGE3D_MAX_WIDTH:\t" << workgroup_size; clGetDeviceInfo(device, CL_DEVICE_IMAGE3D_MAX_HEIGHT, sizeof(workgroup_size), &workgroup_size, nullptr); MITK_INFO("ocl.log")<<" CL_DEVICE_IMAGE3D_MAX_HEIGHT:\t" << workgroup_size; clGetDeviceInfo(device, CL_DEVICE_IMAGE3D_MAX_DEPTH, sizeof(workgroup_size), &workgroup_size, nullptr); MITK_INFO("ocl.log")<<" CL_DEVICE_IMAGE3D_MAX_DEPTH:\t" << workgroup_size; // CL_DEVICE_QUEUE_PROPERTIES cl_command_queue_properties queue_properties; clGetDeviceInfo(device, CL_DEVICE_QUEUE_PROPERTIES, sizeof(queue_properties), &queue_properties, nullptr); if( queue_properties & CL_QUEUE_OUT_OF_ORDER_EXEC_MODE_ENABLE ) MITK_INFO("ocl.log")<<" CL_DEVICE_QUEUE_PROPERTIES:\t\t"<< "CL_QUEUE_OUT_OF_ORDER_EXEC_MODE_ENABLE"; if( queue_properties & CL_QUEUE_PROFILING_ENABLE ) MITK_INFO("ocl.log")<<" CL_DEVICE_QUEUE_PROPERTIES:\t\t"<< "CL_QUEUE_PROFILING_ENABLE"; } std::string GetOclErrorAsString( int _clErr ) { std::string returnString("unkown error number: "+std::to_string(_clErr)+" \n"); switch(_clErr) { case CL_SUCCESS: returnString = "CL_SUCCESS\n"; break; case CL_DEVICE_NOT_FOUND: returnString = "CL_DEVICE_NOT_FOUND\n"; break; case CL_DEVICE_NOT_AVAILABLE: returnString = "CL_DEVICE_NOT_AVAILABLE\n"; break; /*case CL_DEVICE_COMPILER_NOT_AVAILABLE: returnString = "CL_DEVICE_COMPILER_NOT_AVAILABLE\n"; break; */ case CL_MEM_OBJECT_ALLOCATION_FAILURE : returnString = "CL_MEM_OBJECT_ALLOCATION_FAILURE\n"; break; case CL_OUT_OF_RESOURCES: returnString = "CL_OUT_OF_RESOURCES\n"; break; case CL_OUT_OF_HOST_MEMORY: returnString = "CL_OUT_OF_HOST_MEMORY\n"; break; case CL_PROFILING_INFO_NOT_AVAILABLE: returnString = "CL_PROFILING_INFO_NOT_AVAILABLE\n"; break; case CL_MEM_COPY_OVERLAP: returnString = "CL_MEM_COPY_OVERLAP\n"; break; case CL_IMAGE_FORMAT_MISMATCH: returnString = "CL_IMAGE_FORMAT_MISMATCH\n"; break; case CL_IMAGE_FORMAT_NOT_SUPPORTED: returnString = "CL_IMAGE_FORMAT_NOT_SUPPORTED\n"; break; case CL_BUILD_PROGRAM_FAILURE: returnString = "CL_BUILD_PROGRAM_FAILURE\n"; break; case CL_MAP_FAILURE: returnString = "CL_MAP_FAILURE\n"; break; case CL_INVALID_VALUE: returnString = "CL_INVALID_VALUE\n"; break; case CL_INVALID_DEVICE_TYPE: returnString = "CL_INVALID_DEVICE_TYPE\n"; break; case CL_INVALID_PLATFORM: returnString = "CL_INVALID_PLATFORM\n"; break; case CL_INVALID_DEVICE: returnString = "CL_INVALID_DEVICE\n"; break; case CL_INVALID_CONTEXT : returnString = "CL_INVALID_CONTEXT\n"; break; case CL_INVALID_QUEUE_PROPERTIES: returnString = "CL_INVALID_QUEUE_PROPERTIES\n"; break; case CL_INVALID_COMMAND_QUEUE: returnString = "CL_INVALID_COMMAND_QUEUE\n"; break; case CL_INVALID_HOST_PTR: returnString = "CL_INVALID_HOST_PTR\n"; break; case CL_INVALID_MEM_OBJECT: returnString = "CL_INVALID_MEM_OBJECT\n"; break; case CL_INVALID_IMAGE_FORMAT_DESCRIPTOR: returnString = "CL_INVALID_IMAGE_FORMAT_DESCRIPTOR\n"; break; case CL_INVALID_IMAGE_SIZE: returnString = "CL_INVALID_IMAGE_SIZE\n"; break; case CL_INVALID_SAMPLER : returnString = "CL_INVALID_SAMPLER\n"; break; case CL_INVALID_BINARY: returnString = "CL_INVALID_BINARY\n"; break; case CL_INVALID_BUILD_OPTIONS: returnString = "CL_INVALID_BUILD_OPTIONS\n"; break; case CL_INVALID_PROGRAM: returnString = "CL_INVALID_PROGRAM\n"; break; case CL_INVALID_PROGRAM_EXECUTABLE: returnString = "CL_INVALID_PROGRAM_EXECUTABLE\n"; break; case CL_INVALID_KERNEL_NAME: returnString = "CL_INVALID_KERNEL_NAME\n"; break; case CL_INVALID_KERNEL_DEFINITION: returnString = "CL_INVALID_KERNEL_DEFINITION\n"; break; case CL_INVALID_KERNEL : returnString = "CL_INVALID_KERNEL\n"; break; case CL_INVALID_ARG_INDEX : returnString = "CL_INVALID_ARG_INDEX\n"; break; case CL_INVALID_ARG_VALUE : returnString = "CL_INVALID_ARG_VALUE\n"; break; case CL_INVALID_ARG_SIZE : returnString = "CL_INVALID_ARG_SIZE\n"; break; case CL_INVALID_KERNEL_ARGS : returnString = "CL_INVALID_KERNEL_ARGS\n"; break; case CL_INVALID_WORK_DIMENSION: returnString = "CL_INVALID_WORK_DIMENSION\n"; break; case CL_INVALID_WORK_GROUP_SIZE: returnString = "CL_INVALID_WORK_GROUP_SIZE\n"; break; case CL_INVALID_WORK_ITEM_SIZE: returnString = "CL_INVALID_WORK_ITEM_SIZE\n"; break; case CL_INVALID_GLOBAL_OFFSET: returnString = "CL_INVALID_GLOBAL_OFFSET\n"; break; case CL_INVALID_EVENT_WAIT_LIST: returnString = "CL_INVALID_EVENT_WAIT_LIST\n"; break; case CL_INVALID_EVENT: returnString = "CL_INVALID_EVENT\n"; break; case CL_INVALID_OPERATION: returnString = "CL_INVALID_OPERATION\n"; break; case CL_INVALID_GL_OBJECT: returnString = "CL_INVALID_GL_OBJECT\n"; break; case CL_INVALID_BUFFER_SIZE : returnString = "CL_INVALID_BUFFER_SIZE\n"; break; case CL_INVALID_MIP_LEVEL : returnString = "CL_INVALID_MIP_LEVEL\n"; break; default: break; } return returnString; } void GetOclError(int _clErr) { if(_clErr == CL_SUCCESS) MITK_WARN << "Called GetOclErr() with no error value: [CL_SUCCESS]"; else MITK_ERROR << GetOclErrorAsString(_clErr); } bool oclCheckError(int _err, const char* filepath, int lineno) { if (_err) { MITK_ERROR<< "OpenCL Error at " << filepath <<":"<< lineno; GetOclError(_err); return 0; } return 1; } void GetSupportedImageFormats(cl_context _context, cl_mem_object_type _type) { const unsigned int entries = 500; cl_image_format* formats = new cl_image_format[entries]; cl_uint _written = 0; // OpenCL constant to catch error IDs cl_int ciErr1; // Get list of supported image formats for READ_ONLY access ciErr1 = clGetSupportedImageFormats( _context, CL_MEM_READ_ONLY, _type, entries, formats, &_written); CHECK_OCL_ERR(ciErr1); MITK_INFO << "Supported Image Formats, Image: CL_MEM_READ_ONLY \n"; for (unsigned int i=0; i<_written; i++) { MITK_INFO<< "ChannelType: " << GetImageTypeAsString(formats[i].image_channel_data_type) << "| ChannelOrder: "<< GetImageTypeAsString(formats[i].image_channel_order) <<"\n"; } _written = 0; // Get list of supported image formats for READ_WRITE access ciErr1 = clGetSupportedImageFormats( _context, CL_MEM_READ_WRITE, _type, entries, formats, &_written); CHECK_OCL_ERR(ciErr1); MITK_INFO << "Supported Image Formats, Image: CL_MEM_READ_WRITE (found: " << _written <<") \n"; for (unsigned int i=0; i<_written; i++) { MITK_INFO<< "ChannelType: " << GetImageTypeAsString(formats[i].image_channel_data_type) << "| ChannelOrder: "<< GetImageTypeAsString(formats[i].image_channel_order) <<"\n"; } _written = 0; // Get list of supported image formats for WRITE_ONLY access ciErr1 = clGetSupportedImageFormats( _context, CL_MEM_WRITE_ONLY, _type, entries, formats, &_written); CHECK_OCL_ERR(ciErr1); MITK_INFO << "Supported Image Formats, Image: CL_MEM_WRITE_ONLY (found: " << _written <<") \n"; for (unsigned int i=0; i<_written; i++) { MITK_INFO<< "ChannelType: " << GetImageTypeAsString(formats[i].image_channel_data_type) << "| ChannelOrder: "<< GetImageTypeAsString(formats[i].image_channel_order) <<"\n"; } } std::string GetImageTypeAsString( const unsigned int _in) { switch(_in) { case CL_R: return "CL_R "; break; case CL_A: return "CL_A "; break; case CL_RG: return "CL_RG "; break; case CL_RA: return "CL_RA "; break; case CL_RGB: return "CL_RGB "; break; case CL_RGBA: return "CL_RGBA "; break; case CL_BGRA: return "CL_BGRA "; break; case CL_ARGB: return "CL_ARGB "; break; case CL_INTENSITY: return "CL_INTENSITY "; break; case CL_LUMINANCE: return "CL_LUMINANCE "; break; case CL_SNORM_INT8: return "CL_SNORM_INT8 "; break; case CL_SNORM_INT16: return "CL_SNORM_INT16 "; break; case CL_UNORM_INT8: return "CL_UNORM_INT8 "; break; case CL_UNORM_INT16: return "CL_UNORM_INT16 "; break; case CL_UNORM_SHORT_565: return "CL_UNORM_SHORT_565 "; break; case CL_UNORM_SHORT_555: return "CL_UNORM_SHORT_555 "; break; case CL_UNORM_INT_101010: return "CL_UNORM_INT_101010 "; break; case CL_SIGNED_INT8: return "CL_SIGNED_INT8 "; break; case CL_SIGNED_INT16: return "CL_SIGNED_INT16 "; break; case CL_SIGNED_INT32: return "CL_SIGNED_INT32 "; break; case CL_UNSIGNED_INT8: return "CL_UNSIGNED_INT8 "; break; case CL_UNSIGNED_INT16: return "CL_UNSIGNED_INT16 "; break; case CL_UNSIGNED_INT32: return "CL_UNSIGNED_INT32 "; break; case CL_HALF_FLOAT: return "CL_HALF_FLOAT "; break; case CL_FLOAT: return "CL_FLOAT "; break; default: return "--"; break; } } void oclLogBinary(cl_program clProg, cl_device_id clDev) { // Grab the number of devices associated with the program cl_uint num_devices; clGetProgramInfo(clProg, CL_PROGRAM_NUM_DEVICES, sizeof(cl_uint), &num_devices, nullptr); // Grab the device ids cl_device_id* devices = (cl_device_id*) malloc(num_devices * sizeof(cl_device_id)); clGetProgramInfo(clProg, CL_PROGRAM_DEVICES, num_devices * sizeof(cl_device_id), devices, 0); // Grab the sizes of the binaries size_t* binary_sizes = (size_t*)malloc(num_devices * sizeof(size_t)); clGetProgramInfo(clProg, CL_PROGRAM_BINARY_SIZES, num_devices * sizeof(size_t), binary_sizes, nullptr); // Now get the binaries char** ptx_code = (char**)malloc(num_devices * sizeof(char*)); for( unsigned int i=0; i= 0) { - for (short l_s2 = l_s1 + 1; l_s2 < maxLine; ++l_s2) - { - AddSample2 = AddSamples[MemoryStartAccessPoint + l_s2 - minLine]; - if (AddSample1 < inputS && AddSample1 >= 0) { - mult = apodArray[(int)((l_s2 - minLine)*apod_mult)] * - dSource[(int)(globalPosZ * inputL * inputS + AddSample2 * inputL + l_s2)] - * apodArray[(int)((l_s1 - minLine)*apod_mult)] * - dSource[(int)(globalPosZ * inputL * inputS + AddSample1 * inputL + l_s1)]; - - output += sqrt(mult * ((float)(mult>0)-(float)(mult<0))) * ((mult > 0) - (mult < 0)); - } - } + AddSample = AddSamples[MemoryStartAccessPoint + l_s - minLine]; + if (AddSample < inputS && AddSample >= 0) { + output += apodArray[(int)((l_s - minLine)*apod_mult)] * dSource[(int)(globalPosZ * inputL * inputS + AddSample * inputL + l_s)]; } else --curUsedLines; } - dDest[ globalPosZ * outputL * outputS + globalPosY * outputL + globalPosX ] = output / (pow((float)curUsedLines, 2.0f) - (curUsedLines - 1)); + dDest[ globalPosZ * outputL * outputS + globalPosY * outputL + globalPosX ] = output / (float)curUsedLines; } } \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/Resources/DMAS.cl b/Modules/PhotoacousticsAlgorithms/Resources/DMAS.cl index 28d3cee888..f2ddcd5088 100644 --- a/Modules/PhotoacousticsAlgorithms/Resources/DMAS.cl +++ b/Modules/PhotoacousticsAlgorithms/Resources/DMAS.cl @@ -1,77 +1,77 @@ /*=================================================================== 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 int* memoryLocations, __global unsigned short* AddSamples, - __global float* apodArray, + __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 ) { 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 AddSample1 = 0; unsigned short AddSample2 = 0; float output = 0; float mult = 0; unsigned int MemoryStartAccessPoint = memoryLocations[globalPosY * outputL + globalPosX]; for (short l_s1 = minLine; l_s1 < maxLine; ++l_s1) { AddSample1 = AddSamples[MemoryStartAccessPoint + l_s1 - minLine]; if (AddSample1 < inputS && AddSample1 >= 0) { for (short l_s2 = l_s1 + 1; l_s2 < maxLine; ++l_s2) { AddSample2 = AddSamples[MemoryStartAccessPoint + l_s2 - minLine]; if (AddSample1 < inputS && AddSample1 >= 0) { mult = apodArray[(int)((l_s2 - minLine)*apod_mult)] * dSource[(int)(globalPosZ * inputL * inputS + AddSample2 * inputL + l_s2)] * apodArray[(int)((l_s1 - minLine)*apod_mult)] * dSource[(int)(globalPosZ * inputL * inputS + AddSample1 * inputL + l_s1)]; output += sqrt(mult * ((float)(mult>0)-(float)(mult<0))) * ((mult > 0) - (mult < 0)); } } } else --curUsedLines; } dDest[ globalPosZ * outputL * outputS + globalPosY * outputL + globalPosX ] = output / (pow((float)curUsedLines, 2.0f) - (curUsedLines - 1)); } } \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/Resources/ParallelSum.cl b/Modules/PhotoacousticsAlgorithms/Resources/ParallelSum.cl deleted file mode 100644 index 8d640492b9..0000000000 --- a/Modules/PhotoacousticsAlgorithms/Resources/ParallelSum.cl +++ /dev/null @@ -1,36 +0,0 @@ -/*=================================================================== - -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 ckParallelSum( __const const unsigned short *input, - __global unsigned int *partialSums, - __local unsigned int *localSums) - { - uint local_id = get_local_id(0) + get_local_id(1) * get_local_size(2); - uint group_size = get_local_size(0) * get_local_size(1); - - localSums[local_id] = input[get_global_id(0) + get_global_id(1) * get_global_size(0)]; - - for (uint stride = group_size/2; stride>0; stride /=2) - { - barrier(CLK_LOCAL_MEM_FENCE); - - if (local_id < stride) - localSums[local_id] += localSums[local_id + stride]; - } - - if (local_id == 0) - partialSums[get_group_id(0)] = localSums[0]; - } \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/files.cmake b/Modules/PhotoacousticsAlgorithms/files.cmake index 7c95415103..832b281665 100644 --- a/Modules/PhotoacousticsAlgorithms/files.cmake +++ b/Modules/PhotoacousticsAlgorithms/files.cmake @@ -1,26 +1,25 @@ 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 source/OpenCLFilter/mitkPhotoacousticOCLMemoryLocSum.cpp ) set(RESOURCE_FILES DASQuadratic.cl - DMASQuadratic.cl DASSpherical.cl - DMASSpherical.cl BModeAbs.cl BModeAbsLog.cl UsedLinesCalculation.cl MemoryLocSum.cl DelayCalculation.cl DMAS.cl + DAS.cl ) \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.h b/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.h index 20e581b634..8ba82a8160 100644 --- a/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.h +++ b/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.h @@ -1,132 +1,140 @@ /*=================================================================== 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 _MITKPHOTOACOUSTICSOCLBEAMFORMER_H_ #define _MITKPHOTOACOUSTICSOCLBEAMFORMER_H_ #ifdef PHOTOACOUSTICS_USE_GPU #include "mitkOclDataSetToDataSetFilter.h" #include #include "mitkPhotoacousticOCLDelayCalculation.h" #include "mitkPhotoacousticOCLMemoryLocSum.h" #include "mitkPhotoacousticOCLUsedLinesCalculation.h" #include "mitkPhotoacousticBeamformingSettings.h" +#include + namespace mitk { /** Documentation * * \brief The OclBinaryThresholdImageFilter computes a binary segmentation based on given threshold values. * * The filter requires two threshold values ( the upper and the lower threshold ) and two image values ( inside and outside ). The resulting voxel of the segmentation image is assigned the inside value 1 if the image value is between the given thresholds and the outside value otherwise. */ class PhotoacousticOCLBeamformingFilter : public OclDataSetToDataSetFilter, public itk::Object { public: mitkClassMacroItkParent(PhotoacousticOCLBeamformingFilter, itk::Object); itkNewMacro(Self); /** * @brief SetInput Set the input data through an image. Arbitrary images are supported */ void SetInput(Image::Pointer image); /** * @brief SetInput Manually set the input data while providing dimensions and memory size of the input data. */ void SetInput(void* data, unsigned int* dimensions, unsigned int BpE); + void* GetOutput(); + /** * @brief GetOutputAsImage Returns an mitk::Image constructed from the processed data */ mitk::Image::Pointer GetOutputAsImage(); /** Update the filter */ void Update(); /** Set the Apodisation function to apply when beamforming */ void SetApodisation(float* apodisation, unsigned short apodArraySize) { m_ApodArraySize = apodArraySize; m_Apodisation = apodisation; } void SetConfig(BeamformingSettings settings) { m_ConfOld = m_Conf; m_Conf = settings; } protected: /** Constructor */ PhotoacousticOCLBeamformingFilter(); /** Destructor */ virtual ~PhotoacousticOCLBeamformingFilter(); /** Initialize the filter */ bool Initialize(); + void UpdateDataBuffers(); + void Execute(); mitk::PixelType GetOutputType() { return mitk::MakeScalarPixelType(); } int GetBytesPerElem() { return sizeof(float); } virtual us::Module* GetModule(); private: /** The OpenCL kernel for the filter */ cl_kernel m_PixelCalculation; unsigned int m_OutputDim[3]; float* m_Apodisation; unsigned short m_ApodArraySize; unsigned short m_PAImage; BeamformingSettings m_Conf; BeamformingSettings m_ConfOld; mitk::Image::Pointer m_InputImage; size_t m_ChunkSize[3]; mitk::OCLMemoryLocSum::Pointer m_SumFilter; mitk::OCLUsedLinesCalculation::Pointer m_UsedLinesCalculation; mitk::OCLDelayCalculation::Pointer m_DelayCalculation; cl_mem m_ApodizationBuffer; cl_mem m_UsedLinesBuffer; cl_mem m_MemoryLocationsBuffer; cl_mem m_DelaysBuffer; + + std::chrono::steady_clock::time_point m_Begin; }; } #endif #endif diff --git a/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticBeamformingSettings.h b/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticBeamformingSettings.h index c4f2f57ba0..117f2919c8 100644 --- a/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticBeamformingSettings.h +++ b/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticBeamformingSettings.h @@ -1,82 +1,76 @@ /*=================================================================== 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 { //##Documentation //## @brief //## @ingroup Process class BeamformingSettings { public: float Pitch = 0.0003; // [m] float SpeedOfSound = 1540; // [m/s] unsigned int SamplesPerLine = 2048; unsigned int ReconstructionLines = 128; float RecordTime = 0.00006; // [s] float TimeSpacing = 0.0000000000001; // [s] unsigned short TransducerElements = 128; bool partial = false; unsigned int CropBounds[2] = { 0,0 }; unsigned int inputDim[3] = { 1,1,1 }; + unsigned int upperCutoff = 0; bool UseGPU = true; enum DelayCalc { QuadApprox, Spherical }; DelayCalc DelayCalculationMethod = QuadApprox; enum Apodization { Hamm, Hann, Box }; Apodization Apod = Hann; enum BeamformingAlgorithm { DMAS, DAS }; BeamformingAlgorithm Algorithm = DAS; float Angle = 10; bool isPhotoacousticImage = true; float BPHighPass = 50; float BPLowPass = 50; bool UseBP = false; - friend bool operator==(const BeamformingSettings& lhs, const BeamformingSettings& rhs) + //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 (lhs.Algorithm == rhs.Algorithm) && - (lhs.Angle == rhs.Angle) && + return !(((lhs.Angle - rhs.Angle) < 0.001f) && (lhs.Apod == rhs.Apod) && - (lhs.BPHighPass == rhs.BPHighPass) && - (lhs.BPLowPass == rhs.BPLowPass) && - (lhs.CropBounds == rhs.CropBounds) && (lhs.DelayCalculationMethod == rhs.DelayCalculationMethod) && (lhs.isPhotoacousticImage == rhs.isPhotoacousticImage) && - (lhs.inputDim[0] == rhs.inputDim[0]) && - (lhs.inputDim[1] == rhs.inputDim[1]) && - (lhs.inputDim[2] == rhs.inputDim[2]) && - (lhs.partial == rhs.partial) && - (lhs.Pitch == rhs.Pitch) && + ((lhs.Pitch - rhs.Pitch) < 0.0000001f) && (lhs.ReconstructionLines == rhs.ReconstructionLines) && - (lhs.RecordTime == rhs.RecordTime) && + (lhs.RecordTime - rhs.RecordTime) < 0.00000001f && (lhs.SamplesPerLine == rhs.SamplesPerLine) && - (lhs.SpeedOfSound == rhs.SpeedOfSound) && - (lhs.TimeSpacing == rhs.TimeSpacing) && - (lhs.TransducerElements == rhs.TransducerElements) && - (lhs.UseBP == rhs.UseBP) && - (lhs.UseGPU == rhs.UseGPU); + ((lhs.SpeedOfSound - rhs.SpeedOfSound) < 0.01f) && + ((lhs.TimeSpacing - rhs.TimeSpacing) < 0.0000000001f) && + (lhs.TransducerElements == rhs.TransducerElements)); } }; } #endif \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticImage.h b/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticImage.h index af25880f72..9f00543a22 100644 --- a/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticImage.h +++ b/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticImage.h @@ -1,50 +1,53 @@ /*=================================================================== 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 mitkPhotoacousticImage_H_HEADER_INCLUDED #define mitkPhotoacousticImage_H_HEADER_INCLUDED #include "itkObject.h" #include "mitkCommon.h" #include "mitkImage.h" #include #include "mitkPhotoacousticBeamformingSettings.h" +#include "mitkPhotoacousticBeamformingFilter.h" #include "MitkPhotoacousticsAlgorithmsExports.h" namespace mitk { class MITKPHOTOACOUSTICSALGORITHMS_EXPORT PhotoacousticImage : public itk::Object { public: mitkClassMacroItkParent(mitk::PhotoacousticImage, itk::Object); itkFactorylessNewMacro(Self); enum BModeMethod { ShapeDetection, Abs }; mitk::Image::Pointer ApplyBmodeFilter(mitk::Image::Pointer inputImage, BModeMethod method = BModeMethod::Abs, bool UseGPU = false, bool UseLogFilter = false, float resampleSpacing = 0.15); // mitk::Image::Pointer ApplyScatteringCompensation(mitk::Image::Pointer inputImage, int scatteringCoefficient); mitk::Image::Pointer ApplyResampling(mitk::Image::Pointer inputImage, unsigned int outputSize[2]); - mitk::Image::Pointer ApplyBeamforming(mitk::Image::Pointer inputImage, BeamformingSettings config, int cutoff, std::function progressHandle = [](int, std::string) {}); + mitk::Image::Pointer ApplyBeamforming(mitk::Image::Pointer inputImage, BeamformingSettings config, std::function progressHandle = [](int, std::string) {}); mitk::Image::Pointer ApplyCropping(mitk::Image::Pointer inputImage, int above, int below, int right, int left, int minSlice, int maxSlice); mitk::Image::Pointer BandpassFilter(mitk::Image::Pointer data, float recordTime, float BPHighPass, float BPLowPass, float alpha); protected: PhotoacousticImage(); virtual ~PhotoacousticImage(); + mitk::BeamformingFilter::Pointer m_BeamformingFilter; + itk::Image::Pointer BPFunction(mitk::Image::Pointer reference, int cutoffFrequencyPixelHighPass, int cutoffFrequencyPixelLowPass, float alpha); }; } // namespace mitk #endif /* mitkPhotoacousticImage_H_HEADER_INCLUDED */ diff --git a/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.cpp b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.cpp index 6de84a0732..d37f4222bc 100644 --- a/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.cpp @@ -1,279 +1,261 @@ /*=================================================================== 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. ===================================================================*/ #ifdef PHOTOACOUSTICS_USE_GPU #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("DASQuadratic.cl"); - this->AddSourceFile("DMASQuadratic.cl"); + this->AddSourceFile("DAS.cl"); this->AddSourceFile("DMAS.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_SumFilter = mitk::OCLMemoryLocSum::New(); 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::Execute() +void mitk::PhotoacousticOCLBeamformingFilter::UpdateDataBuffers() { - cl_int clErr = 0; - - if (!(m_Conf == m_ConfOld)) + if (BeamformingSettings::SettingsChangedOpenCL(m_Conf, m_ConfOld)) { + cl_int clErr = 0; + MITK_INFO << "Updating Data for new configuration"; //Initialize the Output try { size_t outputSize = m_Conf.ReconstructionLines * m_Conf.SamplesPerLine * 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; } // 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(); // calculate memory locations of the Delays unsigned int sumDimensions[2] = { m_Conf.ReconstructionLines, m_Conf.SamplesPerLine }; m_SumFilter->SetInput(m_UsedLinesCalculation->GetGPUOutput()); m_SumFilter->SetInputDimensions(sumDimensions); m_SumFilter->Update(); unsigned int size = m_SumFilter->GetSum(); m_UsedLinesBuffer = m_UsedLinesCalculation->GetGPUOutput()->GetGPUBuffer(); m_MemoryLocationsBuffer = m_SumFilter->GetGPUOutput()->GetGPUBuffer(); // calculate the Delays m_DelayCalculation->SetConfig(m_Conf); m_DelayCalculation->SetInputs(m_UsedLinesBuffer, m_MemoryLocationsBuffer, size); m_DelayCalculation->Update(); m_DelaysBuffer = m_DelayCalculation->GetGPUOutput()->GetGPUBuffer(); m_ConfOld = m_Conf; } +} - if (m_Conf.Algorithm == BeamformingSettings::BeamformingAlgorithm::DAS) - { - m_PAImage = (unsigned short)m_Conf.isPhotoacousticImage; - - // set kernel arguments - clErr = clSetKernelArg(this->m_PixelCalculation, 2, sizeof(cl_mem), &(this->m_ApodizationBuffer)); - clErr |= clSetKernelArg(this->m_PixelCalculation, 3, sizeof(cl_ushort), &(this->m_ApodArraySize)); - clErr |= clSetKernelArg(this->m_PixelCalculation, 4, sizeof(cl_float), &(this->m_Conf.SpeedOfSound)); - clErr |= clSetKernelArg(this->m_PixelCalculation, 5, sizeof(cl_float), &(this->m_Conf.TimeSpacing)); - clErr |= clSetKernelArg(this->m_PixelCalculation, 6, sizeof(cl_float), &(this->m_Conf.Pitch)); - clErr |= clSetKernelArg(this->m_PixelCalculation, 7, sizeof(cl_float), &(this->m_Conf.Angle)); - clErr |= clSetKernelArg(this->m_PixelCalculation, 8, sizeof(cl_ushort), &(this->m_PAImage)); - clErr |= clSetKernelArg(this->m_PixelCalculation, 9, sizeof(cl_ushort), &(this->m_Conf.TransducerElements)); - clErr |= clSetKernelArg(this->m_PixelCalculation, 10, sizeof(cl_uint), &(this->m_Conf.inputDim[0])); - clErr |= clSetKernelArg(this->m_PixelCalculation, 11, sizeof(cl_uint), &(this->m_Conf.inputDim[1])); - clErr |= clSetKernelArg(this->m_PixelCalculation, 12, sizeof(cl_uint), &(this->m_Conf.inputDim[2])); - clErr |= clSetKernelArg(this->m_PixelCalculation, 13, sizeof(cl_uint), &(this->m_Conf.ReconstructionLines)); - clErr |= clSetKernelArg(this->m_PixelCalculation, 14, sizeof(cl_uint), &(this->m_Conf.SamplesPerLine)); +void mitk::PhotoacousticOCLBeamformingFilter::Execute() +{ + cl_int clErr = 0; + UpdateDataBuffers(); - CHECK_OCL_ERR(clErr); - } - else if (m_Conf.Algorithm == BeamformingSettings::BeamformingAlgorithm::DMAS) - { - clErr = clSetKernelArg(this->m_PixelCalculation, 2, sizeof(cl_mem), &(this->m_UsedLinesBuffer)); - clErr |= clSetKernelArg(this->m_PixelCalculation, 3, sizeof(cl_mem), &(this->m_MemoryLocationsBuffer)); - clErr |= clSetKernelArg(this->m_PixelCalculation, 4, sizeof(cl_mem), &(this->m_DelaysBuffer)); - clErr |= clSetKernelArg(this->m_PixelCalculation, 5, sizeof(cl_mem), &(this->m_ApodizationBuffer)); - clErr |= clSetKernelArg(this->m_PixelCalculation, 6, sizeof(cl_ushort), &(this->m_ApodArraySize)); - clErr |= clSetKernelArg(this->m_PixelCalculation, 7, sizeof(cl_uint), &(this->m_Conf.inputDim[0])); - clErr |= clSetKernelArg(this->m_PixelCalculation, 8, sizeof(cl_uint), &(this->m_Conf.inputDim[1])); - clErr |= clSetKernelArg(this->m_PixelCalculation, 9, sizeof(cl_uint), &(this->m_Conf.inputDim[2])); - clErr |= clSetKernelArg(this->m_PixelCalculation, 10, sizeof(cl_uint), &(this->m_Conf.ReconstructionLines)); - clErr |= clSetKernelArg(this->m_PixelCalculation, 11, sizeof(cl_uint), &(this->m_Conf.SamplesPerLine)); + clErr = clSetKernelArg(this->m_PixelCalculation, 2, sizeof(cl_mem), &(this->m_UsedLinesBuffer)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 3, sizeof(cl_mem), &(this->m_MemoryLocationsBuffer)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 4, sizeof(cl_mem), &(this->m_DelaysBuffer)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 5, sizeof(cl_mem), &(this->m_ApodizationBuffer)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 6, sizeof(cl_ushort), &(this->m_ApodArraySize)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 7, sizeof(cl_uint), &(this->m_Conf.inputDim[0])); + clErr |= clSetKernelArg(this->m_PixelCalculation, 8, sizeof(cl_uint), &(this->m_Conf.inputDim[1])); + clErr |= clSetKernelArg(this->m_PixelCalculation, 9, sizeof(cl_uint), &(this->m_Conf.inputDim[2])); + clErr |= clSetKernelArg(this->m_PixelCalculation, 10, sizeof(cl_uint), &(this->m_Conf.ReconstructionLines)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 11, sizeof(cl_uint), &(this->m_Conf.SamplesPerLine)); - CHECK_OCL_ERR(clErr); - } + CHECK_OCL_ERR(clErr); // execute the filter on a 3D NDRange - if (m_OutputDim[2] == 1) + if (m_OutputDim[2] == 1 || m_ChunkSize[2] == 1) this->ExecuteKernelChunks(m_PixelCalculation, 2, m_ChunkSize); else this->ExecuteKernelChunks(m_PixelCalculation, 3, m_ChunkSize); - + // 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: { - if(m_Conf.DelayCalculationMethod == BeamformingSettings::DelayCalc::QuadApprox) - this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDASQuad", &clErr); - else if (m_Conf.DelayCalculationMethod == BeamformingSettings::DelayCalc::Spherical) - this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDASSphe", &clErr); + this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDAS", &clErr); break; } case BeamformingSettings::BeamformingAlgorithm::DMAS: { - if (m_Conf.DelayCalculationMethod == BeamformingSettings::DelayCalc::QuadApprox) - this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDMAS", &clErr); - else if (m_Conf.DelayCalculationMethod == BeamformingSettings::DelayCalc::Spherical) - this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDMAS", &clErr); + this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDMAS", &clErr); break; } default: { - this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDASQuad", &clErr); + MITK_INFO << "No beamforming algorithm specified, setting to DAS"; + this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDAS", &clErr); break; } } buildErr |= 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::ReferenceMemory); } MITK_DEBUG << "Image Initialized."; return outputImage; } + +void* mitk::PhotoacousticOCLBeamformingFilter::GetOutput() +{ + return OclDataSetToDataSetFilter::GetOutput(); +} #endif diff --git a/Modules/PhotoacousticsAlgorithms/source/mitkPhotoacousticImage.cpp b/Modules/PhotoacousticsAlgorithms/source/mitkPhotoacousticImage.cpp index ba997c88d8..a838dad9fc 100644 --- a/Modules/PhotoacousticsAlgorithms/source/mitkPhotoacousticImage.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/mitkPhotoacousticImage.cpp @@ -1,522 +1,521 @@ /*=================================================================== 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 "mitkPhotoacousticImage.h" #include "../ITKFilter/ITKUltrasound/itkBModeImageFilter.h" #include "../ITKFilter/itkPhotoacousticBModeImageFilter.h" #include "mitkImageCast.h" #include "mitkITKImageImport.h" #include "mitkPhotoacousticBeamformingFilter.h" #include #include #include "./OpenCLFilter/mitkPhotoacousticBModeFilter.h" // itk dependencies #include "itkImage.h" #include "itkResampleImageFilter.h" #include "itkCastImageFilter.h" #include "itkCropImageFilter.h" #include "itkRescaleIntensityImageFilter.h" #include "itkIntensityWindowingImageFilter.h" #include #include "itkMultiplyImageFilter.h" #include "itkBSplineInterpolateImageFunction.h" #include // needed itk image filters #include "mitkITKImageImport.h" #include "itkFFTShiftImageFilter.h" #include "itkMultiplyImageFilter.h" #include "itkComplexToModulusImageFilter.h" #include #include "../ITKFilter/ITKUltrasound/itkFFT1DComplexConjugateToRealImageFilter.h" #include "../ITKFilter/ITKUltrasound/itkFFT1DRealToComplexConjugateImageFilter.h" -mitk::PhotoacousticImage::PhotoacousticImage() +mitk::PhotoacousticImage::PhotoacousticImage() : m_BeamformingFilter(BeamformingFilter::New()) { MITK_INFO << "[PhotoacousticImage Debug] created that image"; } mitk::PhotoacousticImage::~PhotoacousticImage() { MITK_INFO << "[PhotoacousticImage Debug] destroyed that image"; } mitk::Image::Pointer mitk::PhotoacousticImage::ApplyBmodeFilter(mitk::Image::Pointer inputImage, BModeMethod method, bool UseGPU, bool UseLogFilter, float resampleSpacing) { // the image needs to be of floating point type for the envelope filter to work; the casting is done automatically by the CastToItkImage typedef itk::Image< float, 3 > itkFloatImageType; typedef itk::IdentityTransform TransformType; if (method == BModeMethod::Abs) { mitk::Image::Pointer input; mitk::Image::Pointer out; if (inputImage->GetPixelType().GetTypeAsString() == "scalar (float)" || inputImage->GetPixelType().GetTypeAsString() == " (float)") input = inputImage; else input = ApplyCropping(inputImage, 0, 0, 0, 0, 0, inputImage->GetDimension(2) - 1); if (!UseGPU) { PhotoacousticBModeFilter::Pointer filter = PhotoacousticBModeFilter::New(); filter->SetParameters(UseLogFilter); filter->SetInput(input); filter->Update(); out = filter->GetOutput(); if (resampleSpacing == 0) return out; } #ifdef PHOTOACOUSTICS_USE_GPU else { PhotoacousticOCLBModeFilter::Pointer filter = PhotoacousticOCLBModeFilter::New(); filter->SetParameters(UseLogFilter); filter->SetInput(input); filter->Update(); out = filter->GetOutput(); if (resampleSpacing == 0) return out; } #endif typedef itk::ResampleImageFilter < itkFloatImageType, itkFloatImageType > ResampleImageFilter; ResampleImageFilter::Pointer resampleImageFilter = ResampleImageFilter::New(); itkFloatImageType::Pointer itkImage; mitk::CastToItkImage(out, itkImage); itkFloatImageType::SpacingType outputSpacing; itkFloatImageType::SizeType inputSize = itkImage->GetLargestPossibleRegion().GetSize(); itkFloatImageType::SizeType outputSize = inputSize; outputSpacing[0] = itkImage->GetSpacing()[0]; outputSpacing[1] = resampleSpacing; outputSpacing[2] = itkImage->GetSpacing()[2]; outputSize[1] = inputSize[1] * itkImage->GetSpacing()[1] / outputSpacing[1]; typedef itk::IdentityTransform TransformType; resampleImageFilter->SetInput(itkImage); resampleImageFilter->SetSize(outputSize); resampleImageFilter->SetOutputSpacing(outputSpacing); resampleImageFilter->SetTransform(TransformType::New()); resampleImageFilter->UpdateLargestPossibleRegion(); return mitk::GrabItkImageMemory(resampleImageFilter->GetOutput()); } else if (method == BModeMethod::ShapeDetection) { typedef itk::BModeImageFilter < itkFloatImageType, itkFloatImageType > BModeFilterType; BModeFilterType::Pointer bModeFilter = BModeFilterType::New(); // LogFilter typedef itk::PhotoacousticBModeImageFilter < itkFloatImageType, itkFloatImageType > PhotoacousticBModeImageFilter; PhotoacousticBModeImageFilter::Pointer photoacousticBModeFilter = PhotoacousticBModeImageFilter::New(); // No LogFilter typedef itk::ResampleImageFilter < itkFloatImageType, itkFloatImageType > ResampleImageFilter; ResampleImageFilter::Pointer resampleImageFilter = ResampleImageFilter::New(); itkFloatImageType::Pointer itkImage; mitk::CastToItkImage(inputImage, itkImage); itkFloatImageType::Pointer bmode; if (UseLogFilter) { bModeFilter->SetInput(itkImage); bModeFilter->SetDirection(1); bmode = bModeFilter->GetOutput(); } else { photoacousticBModeFilter->SetInput(itkImage); photoacousticBModeFilter->SetDirection(1); bmode = photoacousticBModeFilter->GetOutput(); } // resampleSpacing == 0 means: do no resampling if (resampleSpacing == 0) { return mitk::GrabItkImageMemory(bmode); } itkFloatImageType::SpacingType outputSpacing; itkFloatImageType::SizeType inputSize = itkImage->GetLargestPossibleRegion().GetSize(); itkFloatImageType::SizeType outputSize = inputSize; outputSpacing[0] = itkImage->GetSpacing()[0]; outputSpacing[1] = resampleSpacing; outputSpacing[2] = itkImage->GetSpacing()[2]; outputSize[1] = inputSize[1] * itkImage->GetSpacing()[1] / outputSpacing[1]; resampleImageFilter->SetInput(bmode); resampleImageFilter->SetSize(outputSize); resampleImageFilter->SetOutputSpacing(outputSpacing); resampleImageFilter->SetTransform(TransformType::New()); resampleImageFilter->UpdateLargestPossibleRegion(); return mitk::GrabItkImageMemory(resampleImageFilter->GetOutput()); } return nullptr; } /*mitk::Image::Pointer mitk::PhotoacousticImage::ApplyScatteringCompensation(mitk::Image::Pointer inputImage, int scattering) { typedef itk::Image< float, 3 > itkFloatImageType; typedef itk::MultiplyImageFilter MultiplyImageFilterType; itkFloatImageType::Pointer itkImage; mitk::CastToItkImage(inputImage, itkImage); MultiplyImageFilterType::Pointer multiplyFilter = MultiplyImageFilterType::New(); multiplyFilter->SetInput1(itkImage); multiplyFilter->SetInput2(m_FluenceCompResizedItk.at(m_ScatteringCoefficient)); return mitk::GrabItkImageMemory(multiplyFilter->GetOutput()); }*/ mitk::Image::Pointer mitk::PhotoacousticImage::ApplyResampling(mitk::Image::Pointer inputImage, unsigned int outputSize[2]) { typedef itk::Image< float, 3 > itkFloatImageType; typedef itk::ResampleImageFilter < itkFloatImageType, itkFloatImageType > ResampleImageFilter; ResampleImageFilter::Pointer resampleImageFilter = ResampleImageFilter::New(); typedef itk::LinearInterpolateImageFunction T_Interpolator; itkFloatImageType::Pointer itkImage; mitk::CastToItkImage(inputImage, itkImage); itkFloatImageType::SpacingType outputSpacingItk; itkFloatImageType::SizeType inputSizeItk = itkImage->GetLargestPossibleRegion().GetSize(); itkFloatImageType::SizeType outputSizeItk = inputSizeItk; outputSizeItk[0] = outputSize[0]; outputSizeItk[1] = outputSize[1]; outputSizeItk[2] = inputSizeItk[2]; outputSpacingItk[0] = itkImage->GetSpacing()[0] * (static_cast(inputSizeItk[0]) / static_cast(outputSizeItk[0])); outputSpacingItk[1] = itkImage->GetSpacing()[1] * (static_cast(inputSizeItk[1]) / static_cast(outputSizeItk[1])); outputSpacingItk[2] = itkImage->GetSpacing()[2]; typedef itk::IdentityTransform TransformType; T_Interpolator::Pointer _pInterpolator = T_Interpolator::New(); resampleImageFilter->SetInput(itkImage); resampleImageFilter->SetSize(outputSizeItk); resampleImageFilter->SetOutputSpacing(outputSpacingItk); resampleImageFilter->SetTransform(TransformType::New()); resampleImageFilter->SetInterpolator(_pInterpolator); resampleImageFilter->UpdateLargestPossibleRegion(); return mitk::GrabItkImageMemory(resampleImageFilter->GetOutput()); } mitk::Image::Pointer mitk::PhotoacousticImage::ApplyCropping(mitk::Image::Pointer inputImage, int above, int below, int right, int left, int minSlice, int maxSlice) { unsigned int inputDim[3] = { inputImage->GetDimension(0), inputImage->GetDimension(1), inputImage->GetDimension(2) }; unsigned int outputDim[3] = { inputImage->GetDimension(0) - left - right, inputImage->GetDimension(1) - (unsigned int)above - (unsigned int)below, (unsigned int)maxSlice - (unsigned int)minSlice + 1 }; void* inputData; float* outputData = new float[outputDim[0] * outputDim[1] * outputDim[2]]; ImageReadAccessor acc(inputImage); inputData = const_cast(acc.GetData()); // convert the data to float by default // as of now only those float, short, float are used at all... though it's easy to add other ones if (inputImage->GetPixelType().GetTypeAsString() == "scalar (float)" || inputImage->GetPixelType().GetTypeAsString() == " (float)") { // copy the data into the cropped image for (unsigned short sl = 0; sl < outputDim[2]; ++sl) { for (unsigned short l = 0; l < outputDim[0]; ++l) { for (unsigned short s = 0; s < outputDim[1]; ++s) { outputData[l + s*(unsigned short)outputDim[0] + sl*outputDim[0] * outputDim[1]] = (float)((float*)inputData)[(l + left) + (s + above)*(unsigned short)inputDim[0] + (sl + minSlice)*inputDim[0] * inputDim[1]]; } } } } else if (inputImage->GetPixelType().GetTypeAsString() == "scalar (short)" || inputImage->GetPixelType().GetTypeAsString() == " (short)") { // copy the data unsigned shorto the cropped image for (unsigned short sl = 0; sl < outputDim[2]; ++sl) { for (unsigned short l = 0; l < outputDim[0]; ++l) { for (unsigned short s = 0; s < outputDim[1]; ++s) { outputData[l + s*(unsigned short)outputDim[0] + sl*outputDim[0] * outputDim[1]] = (float)((short*)inputData)[(l + left) + (s + above)*(unsigned short)inputDim[0] + (sl + minSlice)*inputDim[0] * inputDim[1]]; } } } } else if (inputImage->GetPixelType().GetTypeAsString() == "scalar (double)" || inputImage->GetPixelType().GetTypeAsString() == " (double)") { // copy the data unsigned shorto the cropped image for (unsigned short sl = 0; sl < outputDim[2]; ++sl) { for (unsigned short l = 0; l < outputDim[0]; ++l) { for (unsigned short s = 0; s < outputDim[1]; ++s) { outputData[l + s*(unsigned short)outputDim[0] + sl*outputDim[0] * outputDim[1]] = (float)((double*)inputData)[(l + left) + (s + above)*(unsigned short)inputDim[0] + (sl + minSlice)*inputDim[0] * inputDim[1]]; } } } } else { MITK_INFO << "Could not determine pixel type"; } mitk::Image::Pointer output = mitk::Image::New(); output->Initialize(mitk::MakeScalarPixelType(), 3, outputDim); output->SetSpacing(inputImage->GetGeometry()->GetSpacing()); output->SetImportVolume(outputData, 0, 0, mitk::Image::ReferenceMemory); return output; } -mitk::Image::Pointer mitk::PhotoacousticImage::ApplyBeamforming(mitk::Image::Pointer inputImage, BeamformingSettings config, int cutoff, std::function progressHandle) +mitk::Image::Pointer mitk::PhotoacousticImage::ApplyBeamforming(mitk::Image::Pointer inputImage, BeamformingSettings config, std::function progressHandle) { - config.RecordTime = config.RecordTime - (cutoff) / inputImage->GetDimension(1) * config.RecordTime; // adjust the recorded time lost by cropping + config.RecordTime = config.RecordTime - (float)(config.upperCutoff) / (float)inputImage->GetDimension(1) * config.RecordTime; // adjust the recorded time lost by cropping progressHandle(0, "cropping image"); if (!config.partial) { config.CropBounds[0] = 0; config.CropBounds[1] = inputImage->GetDimension(2) - 1; } - Image::Pointer processedImage = ApplyCropping(inputImage, cutoff, 0, 0, 0, config.CropBounds[0], config.CropBounds[1]); + Image::Pointer processedImage = ApplyCropping(inputImage, config.upperCutoff, 0, 0, 0, config.CropBounds[0], config.CropBounds[1]); config.inputDim[0] = processedImage->GetDimension(0); config.inputDim[1] = processedImage->GetDimension(1); config.inputDim[2] = processedImage->GetDimension(2); // perform the beamforming - BeamformingFilter::Pointer Beamformer = BeamformingFilter::New(); - Beamformer->SetInput(processedImage); - Beamformer->Configure(config); - Beamformer->SetProgressHandle(progressHandle); - Beamformer->UpdateLargestPossibleRegion(); + m_BeamformingFilter->SetInput(processedImage); + m_BeamformingFilter->Configure(config); + m_BeamformingFilter->SetProgressHandle(progressHandle); + m_BeamformingFilter->UpdateLargestPossibleRegion(); - processedImage = Beamformer->GetOutput(); + processedImage = m_BeamformingFilter->GetOutput(); return processedImage; } mitk::Image::Pointer mitk::PhotoacousticImage::BandpassFilter(mitk::Image::Pointer data, float recordTime, float BPHighPass, float BPLowPass, float alpha) { bool powerOfTwo = false; int finalPower = 0; for (int i = 1; pow(2, i) <= data->GetDimension(1); ++i) { finalPower = i; if (pow(2, i) == data->GetDimension(1)) { powerOfTwo = true; } } if (!powerOfTwo) { unsigned int dim[2] = { data->GetDimension(0), (unsigned int)pow(2,finalPower+1)}; data = ApplyResampling(data, dim); } MITK_INFO << data->GetDimension(0); // do a fourier transform, multiply with an appropriate window for the filter, and transform back typedef float PixelType; typedef itk::Image< PixelType, 3 > RealImageType; RealImageType::Pointer image; mitk::CastToItkImage(data, image); typedef itk::FFT1DRealToComplexConjugateImageFilter ForwardFFTFilterType; typedef ForwardFFTFilterType::OutputImageType ComplexImageType; ForwardFFTFilterType::Pointer forwardFFTFilter = ForwardFFTFilterType::New(); forwardFFTFilter->SetInput(image); forwardFFTFilter->SetDirection(1); try { forwardFFTFilter->UpdateOutputInformation(); } catch (itk::ExceptionObject & error) { std::cerr << "Error: " << error << std::endl; MITK_WARN << "Bandpass could not be applied"; return data; } float singleVoxel = 1 / (recordTime / data->GetDimension(1)) / 2 / 1000; float cutoffPixelHighPass = std::min(BPHighPass / singleVoxel, (float)data->GetDimension(1) / 2); float cutoffPixelLowPass = std::min(BPLowPass / singleVoxel, (float)data->GetDimension(1) / 2 - cutoffPixelHighPass); RealImageType::Pointer fftMultiplicator = BPFunction(data, cutoffPixelHighPass, cutoffPixelLowPass, alpha); typedef itk::MultiplyImageFilter< ComplexImageType, RealImageType, ComplexImageType > MultiplyFilterType; MultiplyFilterType::Pointer multiplyFilter = MultiplyFilterType::New(); multiplyFilter->SetInput1(forwardFFTFilter->GetOutput()); multiplyFilter->SetInput2(fftMultiplicator); /*itk::ComplexToModulusImageFilter::Pointer toReal = itk::ComplexToModulusImageFilter::New(); toReal->SetInput(forwardFFTFilter->GetOutput()); return GrabItkImageMemory(toReal->GetOutput()); return GrabItkImageMemory(fftMultiplicator); *///DEBUG typedef itk::FFT1DComplexConjugateToRealImageFilter< ComplexImageType, RealImageType > InverseFilterType; InverseFilterType::Pointer inverseFFTFilter = InverseFilterType::New(); inverseFFTFilter->SetInput(multiplyFilter->GetOutput()); inverseFFTFilter->SetDirection(1); return GrabItkImageMemory(inverseFFTFilter->GetOutput()); } itk::Image::Pointer mitk::PhotoacousticImage::BPFunction(mitk::Image::Pointer reference, int cutoffFrequencyPixelHighPass, int cutoffFrequencyPixelLowPass, float alpha) { float* imageData = new float[reference->GetDimension(0)*reference->GetDimension(1)]; // tukey window float width = reference->GetDimension(1) / 2 - (float)cutoffFrequencyPixelHighPass - (float)cutoffFrequencyPixelLowPass; float center = (float)cutoffFrequencyPixelHighPass / 2 + width / 2; MITK_INFO << width << "width " << center << "center " << alpha; for (unsigned int n = 0; n < reference->GetDimension(1); ++n) { imageData[reference->GetDimension(0)*n] = 0; } for (int n = 0; n < width; ++n) { if (n <= (alpha*(width - 1)) / 2) { imageData[reference->GetDimension(0)*(int)(n + center - (width / 2))] = (1 + cos(M_PI*(2 * n / (alpha*(width - 1)) - 1))) / 2; } else if (n >= (width - 1)*(1 - alpha / 2) && n <= (width - 1)) { imageData[reference->GetDimension(0)*(int)(n + center - (width / 2))] = (1 + cos(M_PI*(2 * n / (alpha*(width - 1)) + 1 - 2 / alpha))) / 2; } else { imageData[reference->GetDimension(0)*(int)(n + center - (width / 2))] = 1; } } // Butterworth-Filter /* // first, write the HighPass if (cutoffFrequencyPixelHighPass != reference->GetDimension(1) / 2) { for (int n = 0; n < reference->GetDimension(1) / 2; ++n) { imageData[reference->GetDimension(0)*n] = 1 / (1 + pow( (float)n / (float)(reference->GetDimension(1) / 2 - cutoffFrequencyPixelHighPass) , 2 * butterworthOrder)); } } else { for (int n = 0; n < reference->GetDimension(1) / 2; ++n) { imageData[reference->GetDimension(0)*n] = 1; } } // now, the LowPass for (int n = 0; n < reference->GetDimension(1) / 2; ++n) { imageData[reference->GetDimension(0)*n] *= 1 / (1 + pow( (float)(reference->GetDimension(1) / 2 - 1 - n) / (float)(reference->GetDimension(1) / 2 - cutoffFrequencyPixelLowPass) , 2 * butterworthOrder)); } */ // mirror the first half of the image for (unsigned int n = reference->GetDimension(1) / 2; n < reference->GetDimension(1); ++n) { imageData[reference->GetDimension(0)*n] = imageData[(reference->GetDimension(1) - (n + 1)) * reference->GetDimension(0)]; } // copy and paste to all lines for (unsigned int line = 1; line < reference->GetDimension(0); ++line) { for (unsigned int sample = 0; sample < reference->GetDimension(1); ++sample) { imageData[reference->GetDimension(0)*sample + line] = imageData[reference->GetDimension(0)*sample]; } } typedef itk::Image< float, 3U > ImageType; ImageType::RegionType region; ImageType::IndexType start; start.Fill(0); region.SetIndex(start); ImageType::SizeType size; size[0] = reference->GetDimension(0); size[1] = reference->GetDimension(1); size[2] = reference->GetDimension(2); region.SetSize(size); ImageType::SpacingType SpacingItk; SpacingItk[0] = reference->GetGeometry()->GetSpacing()[0]; SpacingItk[1] = reference->GetGeometry()->GetSpacing()[1]; SpacingItk[2] = reference->GetGeometry()->GetSpacing()[2]; ImageType::Pointer image = ImageType::New(); image->SetRegions(region); image->Allocate(); image->FillBuffer(itk::NumericTraits::Zero); image->SetSpacing(SpacingItk); ImageType::IndexType pixelIndex; for (ImageType::IndexValueType slice = 0; slice < reference->GetDimension(2); ++slice) { for (ImageType::IndexValueType line = 0; line < reference->GetDimension(0); ++line) { for (ImageType::IndexValueType sample = 0; sample < reference->GetDimension(1); ++sample) { pixelIndex[0] = line; pixelIndex[1] = sample; pixelIndex[2] = slice; image->SetPixel(pixelIndex, imageData[line + sample*reference->GetDimension(0)]); } } } delete[] imageData; return image; } \ No newline at end of file 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 215fb7c159..33d9af950b 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,1122 +1,1115 @@ /*=================================================================== 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) +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(UpdateBounds())); connect(m_Controls.boundHigh, SIGNAL(valueChanged(int)), this, SLOT(UpdateBounds())); connect(m_Controls.Partial, SIGNAL(clicked()), this, SLOT(UpdateBounds())); 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{'/'}; - mitk::PhotoacousticImage::Pointer filterbank = mitk::PhotoacousticImage::New(); - 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); - image = filterbank->ApplyBeamforming(image, BFconfig, m_Controls.Cutoff->value(), progressHandle); + image = m_FilterBank->ApplyBeamforming(image, BFconfig, 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 = filterbank->ApplyCropping(image, m_Controls.CutoffAbove->value(), m_Controls.CutoffBelow->value(), 0, 0, 0, image->GetDimension(2) - 1); + 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 = filterbank->BandpassFilter(image, recordTime, BPHighPass, BPLowPass, m_Controls.BPFalloff->value()); + 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 = filterbank->ApplyBmodeFilter(image, mitk::PhotoacousticImage::BModeMethod::Abs, useGPU, m_UseLogfilter, m_ResampleSpacing); + image = m_FilterBank->ApplyBmodeFilter(image, mitk::PhotoacousticImage::BModeMethod::Abs, useGPU, m_UseLogfilter, m_ResampleSpacing); else if (m_Controls.BModeMethod->currentText() == "Envelope Detection") - image = filterbank->ApplyBmodeFilter(image, mitk::PhotoacousticImage::BModeMethod::ShapeDetection, useGPU, m_UseLogfilter, m_ResampleSpacing); + image = m_FilterBank->ApplyBmodeFilter(image, mitk::PhotoacousticImage::BModeMethod::ShapeDetection, 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::finished, thread, &QObject::deleteLater); thread->setConfig(BFconfig); - thread->setCutoff(m_Controls.Cutoff->value()); 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); 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.ProgressInfo->setVisible(false); m_Controls.buttonApplyBeamforming->setText("Apply Beamforming"); EnableControls(); // update rendering mitk::RenderingManager::GetInstance()->InitializeViews( dynamic_cast(data)->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::ShapeDetection, 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()); 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::UpdateBounds() { if (!m_Controls.Partial->isChecked()) { m_Controls.boundLow->setEnabled(false); m_Controls.boundHigh->setEnabled(false); BFconfig.partial = false; return; } else { m_Controls.boundLow->setEnabled(true); m_Controls.boundHigh->setEnabled(true); BFconfig.partial = true; } if(m_Controls.boundLow->value() > m_Controls.boundHigh->value()) { m_Controls.boundHigh->setValue(m_Controls.boundLow->value()); BFconfig.CropBounds[0] = m_Controls.boundLow->value(); BFconfig.CropBounds[1] = m_Controls.boundHigh->value(); } else { BFconfig.CropBounds[0] = m_Controls.boundLow->value(); BFconfig.CropBounds[1] = m_Controls.boundHigh->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::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); } UpdateRecordTime(image); // bandpass configs float speedOfSound = m_Controls.BPSpeedOfSound->value(); // [m/s] float recordTime = image->GetDimension(1)*image->GetGeometry()->GetSpacing()[1] / 1000 / speedOfSound; std::stringstream frequency; float maxFrequency = 1 / (recordTime / image->GetDimension(1)) * 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; 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.Angle = m_Controls.Angle->value(); // [deg] BFconfig.UseBP = m_Controls.UseBP->isChecked(); BFconfig.UseGPU = m_Controls.UseGPUBf->isChecked(); + BFconfig.upperCutoff = m_Controls.Cutoff->value(); UpdateRecordTime(image); UpdateBounds(); } void PAImageProcessing::UpdateRecordTime(mitk::Image::Pointer image) { 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; } } void PAImageProcessing::EnableControls() { m_Controls.DoResampling->setEnabled(true); UseResampling(); m_Controls.Logfilter->setEnabled(true); m_Controls.buttonApplyBModeFilter->setEnabled(true); m_Controls.CutoffAbove->setEnabled(true); m_Controls.CutoffBelow->setEnabled(true); m_Controls.Cutoff->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.DoResampling->setEnabled(false); m_Controls.ResamplingValue->setEnabled(false); m_Controls.Logfilter->setEnabled(false); m_Controls.buttonApplyBModeFilter->setEnabled(false); m_Controls.CutoffAbove->setEnabled(false); m_Controls.CutoffBelow->setEnabled(false); m_Controls.Cutoff->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); } } void BeamformingThread::run() { mitk::Image::Pointer resultImage; - mitk::PhotoacousticImage::Pointer filterbank = mitk::PhotoacousticImage::New(); std::function progressHandle = [this](int progress, std::string progressInfo) { emit updateProgress(progress, progressInfo); }; - resultImage = filterbank->ApplyBeamforming(m_InputImage, m_BFconfig, m_Cutoff, progressHandle); + resultImage = m_FilterBank->ApplyBeamforming(m_InputImage, m_BFconfig, progressHandle); emit result(resultImage); } void BeamformingThread::setConfig(mitk::BeamformingSettings BFconfig) { m_BFconfig = BFconfig; } void BeamformingThread::setInputImage(mitk::Image::Pointer image) { m_InputImage = image; } -void BeamformingThread::setCutoff(int cutoff) -{ - m_Cutoff = cutoff; -} - void BmodeThread::run() { mitk::Image::Pointer resultImage; - mitk::PhotoacousticImage::Pointer filterbank = mitk::PhotoacousticImage::New(); - resultImage = filterbank->ApplyBmodeFilter(m_InputImage, m_Method, m_UseGPU, m_UseLogfilter, m_ResampleSpacing); + 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; - mitk::PhotoacousticImage::Pointer filterbank = mitk::PhotoacousticImage::New(); - resultImage = filterbank->ApplyCropping(m_InputImage, m_CutAbove, m_CutBelow, 0, 0, 0, m_InputImage->GetDimension(2) - 1); + resultImage = m_FilterBank->ApplyCropping(m_InputImage, m_CutAbove, m_CutBelow, 0, 0, 0, m_InputImage->GetDimension(2) - 1); emit result(resultImage); } void CropThread::setConfig(unsigned int CutAbove, unsigned int CutBelow) { m_CutAbove = CutAbove; m_CutBelow = CutBelow; } void CropThread::setInputImage(mitk::Image::Pointer image) { m_InputImage = image; } void BandpassThread::run() { mitk::Image::Pointer resultImage; - mitk::PhotoacousticImage::Pointer filterbank = mitk::PhotoacousticImage::New(); - resultImage = filterbank->BandpassFilter(m_InputImage, m_RecordTime, m_BPHighPass, m_BPLowPass, m_TukeyAlpha); + 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/PAImageProcessing.h b/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessing.h index 8cbf437202..59b649bde8 100644 --- a/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessing.h +++ b/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessing.h @@ -1,187 +1,214 @@ /*=================================================================== 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 PAImageProcessing_h #define PAImageProcessing_h #include #include #include #include #include "ui_PAImageProcessingControls.h" #include "mitkPhotoacousticBeamformingFilter.h" #include "mitkPhotoacousticBeamformingSettings.h" Q_DECLARE_METATYPE(mitk::Image::Pointer) Q_DECLARE_METATYPE(std::string) class PAImageProcessing : public QmitkAbstractView { // this is needed for all Qt objects that should have a Qt meta-object // (everything that derives from QObject and wants to have signal/slots) Q_OBJECT public: static const std::string VIEW_ID; PAImageProcessing(); protected slots: /// \brief Called when the user clicks the GUI button void UpdateBounds(); void UseResampling(); void UseLogfilter(); void SetResampling(); void UseImageSpacing(); void UpdateImageInfo(); void UpdateRecordTime(mitk::Image::Pointer image); void HandleBeamformingResults(mitk::Image::Pointer image); void StartBeamformingThread(); void HandleBmodeResults(mitk::Image::Pointer image); void StartBmodeThread(); void HandleCropResults(mitk::Image::Pointer image); void StartCropThread(); void HandleBandpassResults(mitk::Image::Pointer image); void StartBandpassThread(); void UpdateProgress(int progress, std::string progressInfo); void BatchProcessing(); void UpdateSaveBoxes(); void ChangedSOSBandpass(); void ChangedSOSBeamforming(); protected: virtual void CreateQtPartControl(QWidget *parent) override; virtual void SetFocus() override; /// \brief called by QmitkFunctionality when DataManager's selection has changed virtual void OnSelectionChanged( berry::IWorkbenchPart::Pointer source, const QList& nodes ) override; Ui::PAImageProcessingControls m_Controls; float m_ResampleSpacing; bool m_UseLogfilter; std::string m_OldNodeName; mitk::BeamformingSettings BFconfig; void UpdateBFSettings(mitk::Image::Pointer image); void EnableControls(); void DisableControls(); + + mitk::PhotoacousticImage::Pointer m_FilterBank; }; class BeamformingThread : public QThread { Q_OBJECT void run() Q_DECL_OVERRIDE; signals: void result(mitk::Image::Pointer); void updateProgress(int, std::string); public: void setConfig(mitk::BeamformingSettings BFconfig); void setInputImage(mitk::Image::Pointer image); - void setCutoff(int cutoff); + void setFilterBank(mitk::PhotoacousticImage::Pointer filterBank) + { + m_FilterBank = filterBank; + } + protected: mitk::BeamformingSettings m_BFconfig; mitk::Image::Pointer m_InputImage; int m_Cutoff; + + mitk::PhotoacousticImage::Pointer m_FilterBank; }; class BmodeThread : public QThread { Q_OBJECT void run() Q_DECL_OVERRIDE; signals: void result(mitk::Image::Pointer); public: enum BModeMethod { ShapeDetection, Abs }; void setConfig(bool useLogfilter, double resampleSpacing, mitk::PhotoacousticImage::BModeMethod method, bool useGPU); void setInputImage(mitk::Image::Pointer image); + void setFilterBank(mitk::PhotoacousticImage::Pointer filterBank) + { + m_FilterBank = filterBank; + } + protected: mitk::Image::Pointer m_InputImage; mitk::PhotoacousticImage::BModeMethod m_Method; bool m_UseLogfilter; double m_ResampleSpacing; bool m_UseGPU; + + mitk::PhotoacousticImage::Pointer m_FilterBank; }; class CropThread : public QThread { Q_OBJECT void run() Q_DECL_OVERRIDE; signals: void result(mitk::Image::Pointer); public: void setConfig(unsigned int CutAbove, unsigned int CutBelow); void setInputImage(mitk::Image::Pointer image); + void setFilterBank(mitk::PhotoacousticImage::Pointer filterBank) + { + m_FilterBank = filterBank; + } protected: mitk::Image::Pointer m_InputImage; unsigned int m_CutAbove; unsigned int m_CutBelow; + + mitk::PhotoacousticImage::Pointer m_FilterBank; }; class BandpassThread : public QThread { Q_OBJECT void run() Q_DECL_OVERRIDE; signals: void result(mitk::Image::Pointer); public: void setConfig(float BPHighPass, float BPLowPass, float TukeyAlpha, float recordTime); void setInputImage(mitk::Image::Pointer image); + void setFilterBank(mitk::PhotoacousticImage::Pointer filterBank) + { + m_FilterBank = filterBank; + } protected: mitk::Image::Pointer m_InputImage; float m_BPHighPass; float m_BPLowPass; float m_TukeyAlpha; float m_RecordTime; + + mitk::PhotoacousticImage::Pointer m_FilterBank; }; #endif // PAImageProcessing_h