diff --git a/setup.py b/setup.py index 79540ab..eab8542 100644 --- a/setup.py +++ b/setup.py @@ -40,7 +40,7 @@ setup( name='leapct', - version='1.20', + version='1.21', author='Kyle Champley, Hyojin Kim', author_email='champley@gmail.com, hkim@llnl.gov', description='LivermorE AI Projector for Computed Tomography (LEAPCT)', diff --git a/setup_AMD.py b/setup_AMD.py index 3dd56d7..f5f1285 100644 --- a/setup_AMD.py +++ b/setup_AMD.py @@ -113,7 +113,7 @@ setup( name='leapct', - version='1.20', + version='1.21', author='Kyle Champley, Hyojin Kim', author_email='champley@gmail.com, hkim@llnl.gov', description='LivermorE AI Projector for Computed Tomography (LEAPCT)', diff --git a/setup_ctype.py b/setup_ctype.py index adb5409..9a90b77 100644 --- a/setup_ctype.py +++ b/setup_ctype.py @@ -40,7 +40,7 @@ setup( name='leapct', - version='1.20', + version='1.21', author='Kyle Champley, Hyojin Kim', author_email='champley@gmail.com, hkim@llnl.gov', description='LivermorE AI Projector for Computed Tomography (LEAPCT)', diff --git a/src/cuda_utils.cu b/src/cuda_utils.cu index aa31926..c79ab60 100644 --- a/src/cuda_utils.cu +++ b/src/cuda_utils.cu @@ -798,8 +798,10 @@ float sum(const float* dev_lhs, const int3 N, int whichGPU) uint64 numberOfElements = uint64(N.x) * uint64(N.y) * uint64(N.z); //* Newest method which makes sure all cores are busy int num_gpu_cores = max(1024, getSPcores(whichGPU)); - if (numberOfElements < uint64(num_gpu_cores)) - num_gpu_cores = int(numberOfElements); + if (uint64(sqrt(double(numberOfElements))) < uint64(num_gpu_cores)) + num_gpu_cores = int(sqrt(double(numberOfElements))); + //if (numberOfElements < uint64(num_gpu_cores)) + // num_gpu_cores = int(numberOfElements); uint64 maxNumItemsPerCore = uint64(ceil(double(numberOfElements) / double(num_gpu_cores))); //printf("number of cores = %d, number of chunks = %d\n", num_gpu_cores, int(numChunks)); @@ -859,8 +861,10 @@ float innerProduct(const float* dev_lhs, const float* dev_rhs, const int3 N, int uint64 numberOfElements = uint64(N.x) * uint64(N.y) * uint64(N.z); //* Newest method which makes sure all cores are busy int num_gpu_cores = max(1024, getSPcores(whichGPU)); - if (numberOfElements < uint64(num_gpu_cores)) - num_gpu_cores = int(numberOfElements); + if (uint64(sqrt(double(numberOfElements))) < uint64(num_gpu_cores)) + num_gpu_cores = int(sqrt(double(numberOfElements))); + //if (numberOfElements < uint64(num_gpu_cores)) + // num_gpu_cores = int(numberOfElements); uint64 maxNumItemsPerCore = uint64(ceil(double(numberOfElements) / double(num_gpu_cores))); //printf("number of cores = %d, number of chunks = %d\n", num_gpu_cores, int(numChunks)); diff --git a/src/leap_preprocessing_algorithms.py b/src/leap_preprocessing_algorithms.py index 9b4424a..4fff6c5 100644 --- a/src/leap_preprocessing_algorithms.py +++ b/src/leap_preprocessing_algorithms.py @@ -162,7 +162,7 @@ def gain_correction(leapct, g, air_scan, dark_scan, calibration_scans=None, ROI= #g[:,:,:] = g[:,:,:] / postageStamp[:,None,None] g[:] = g[:] / postageStamp[:,None,None] - badPixelCorrection(leapct, g, badPixelMap, 5, isAttenuationData=False) + badPixelCorrection(leapct, g, None, None, badPixelMap, 5, isAttenuationData=False) return True @@ -265,7 +265,7 @@ def makeAttenuationRadiographs(leapct, g, air_scan=None, dark_scan=None, ROI=Non return True -def badPixelCorrection(leapct, g, badPixelMap=None, windowSize=3, isAttenuationData=True): +def badPixelCorrection(leapct, g, air_scan=None, dark_scan=None, badPixelMap=None, windowSize=3, isAttenuationData=True): r"""Removes bad pixels from CT projections LEAP CT geometry parameters must be set prior to running this function @@ -308,6 +308,10 @@ def badPixelCorrection(leapct, g, badPixelMap=None, windowSize=3, isAttenuationD #print('Estimated ' + str(np.sum(badPixelMap)) + ' bad pixels') leapct.badPixelCorrection(g, badPixelMap, windowSize) + if air_scan is not None: + leapct.badPixelCorrection(air_scan, badPixelMap, windowSize) + if dark_scan is not None: + leapct.badPixelCorrection(dark_scan, badPixelMap, windowSize) if isAttenuationData: leapct.negLog(g) return True @@ -503,10 +507,13 @@ def ringRemoval_fast(leapct, g, delta=0.01, numIter=30, maxChange=0.05): g_sum = np.zeros((1,g.shape[1],g.shape[2]), dtype=np.float32) g_sum[0,:] = np.mean(g,axis=0) g_sum_save = leapct.copyData(g_sum) + minValue = np.min(g_sum_save) + minValue = min(minValue, 0.0) numNeighbors = leapct.get_numTVneighbors() leapct.set_numTVneighbors(6) leapct.diffuse(g_sum,delta,numIter) + g_sum[g_sum= self.get_numZ(): iz = self.get_numZ()//2 + offsetZ_save = self.get_offsetZ() numZ_save = self.get_numZ() + numRows_save = self.get_numRows() + centerRow_save = self.get_centerRow() if self.get_geometry() == 'PARALLEL' or self.get_geometry() == 'FAN': g_chunk = self.cropProjections([iz, iz], None, g) + + #self.set_numRows(rowRange[1]-rowRange[0]+1) + #self.shift_detector(self.get_pixelHeight()*rowRange[0], 0.0) + detectorShift = self.get_pixelHeight()*iz + else: - rowRange = self.rowRangeNeededForBackprojection() + rowRange = self.rowRangeNeededForBackprojection(iz) g_chunk = self.cropProjections(rowRange, None, g) - - self.set_offsetZ(self.z_samples()[iz]) - self.set_numZ(1) + self.set_offsetZ(self.z_samples()[iz]) + self.set_numZ(1) + + #self.set_numRows(rowRange[1]-rowRange[0]+1) + #self.shift_detector(self.get_pixelHeight()*rowRange[0], 0.0) + detectorShift = self.get_pixelHeight()*rowRange[0] f = self.FBP(g_chunk, None, True) del g_chunk self.set_offsetZ(offsetZ_save) self.set_numZ(numZ_save) + self.set_numRows(numRows_save) + self.shift_detector(-detectorShift, 0.0) + #self.set_centerRow(centerRow_save) return f @@ -2653,7 +2673,7 @@ def windowFOV(self, f): self.libprojectors.windowFOV(f, True) return f - def rowRangeNeededForBackprojection(self): + def rowRangeNeededForBackprojection(self, iz=None): r"""Calculates the detector rows necessary to reconstruct the current volume specification The CT geometry parameters and the CT volume parameters must be set prior to running this function. @@ -2666,13 +2686,40 @@ def rowRangeNeededForBackprojection(self): rowsNeeded, a 2X1 numpy array where the values are the first and last detector row index needed to reconstruct the volume. """ + rowsNeeded = np.zeros(2,dtype=np.int32) rowsNeeded[1] = self.get_numRows()-1 + + if self.get_geometry() == 'PARALLEL' or self.get_geometry() == 'FAN': + if iz is not None: + if iz < 0 or iz >= self.get_numZ(): + iz = self.get_numZ()//2 + rowsNeeded[0] = iz + rowsNeeded[1] = iz + return rowsNeeded + self.libprojectors.rowRangeNeededForBackprojection.argtypes = [ndpointer(ctypes.c_int, flags="C_CONTIGUOUS")] self.libprojectors.rowRangeNeededForBackprojection.restype = ctypes.c_bool - self.set_model() - self.libprojectors.rowRangeNeededForBackprojection(rowsNeeded) - return rowsNeeded + if iz is None: + self.set_model() + self.libprojectors.rowRangeNeededForBackprojection(rowsNeeded) + return rowsNeeded + else: + if iz < 0 or iz >= self.get_numZ(): + iz = self.get_numZ()//2 + offsetZ_save = self.get_offsetZ() + numZ_save = self.get_numZ() + + self.set_offsetZ(self.z_samples()[iz]) + self.set_numZ(1) + + self.set_model() + self.libprojectors.rowRangeNeededForBackprojection(rowsNeeded) + + self.set_offsetZ(offsetZ_save) + self.set_numZ(numZ_save) + + return rowsNeeded def viewRangeNeededForBackprojection(self): r"""Calculates the detector projections necessary to reconstruct the current volume specification @@ -2803,6 +2850,9 @@ def crop_rows(self, rowRange, g=None): numRows = self.get_numRows() self.set_numRows(rowRange[1]-rowRange[0]+1) self.shift_detector(self.get_pixelHeight()*rowRange[0], 0.0) + if self.get_geometry() == 'PARALLEL' or self.get_geometry() == 'FAN': + self.set_numZ(self.get_numRows()) + self.set_offsetZ(0.0) if g is not None: if has_torch == True and type(g) is torch.Tensor: @@ -4706,10 +4756,19 @@ def badPixelCorrection(self, g, badPixelMap, windowSize=3): True is successful, False otherwise """ - if len(g.shape) != 3 or g.shape[0] != self.get_numAngles() or g.shape[1] != self.get_numRows() or g.shape[2] != self.get_numCols(): - print('Error: input data dimensions do not match CT data dimensions') - return False - if len(badPixelMap.shape) != 2 or g.shape[1] != badPixelMap.shape[0] or g.shape[2] != badPixelMap.shape[1]: + if len(g.shape) == 3: + numAngles = g.shape[0] + numRows = g.shape[1] + numCols = g.shape[2] + elif len(g.shape) == 2: + numAngles = 1 + numRows = g.shape[0] + numCols = g.shape[1] + + #if len(g.shape) != 3 or g.shape[0] != self.get_numAngles() or g.shape[1] != self.get_numRows() or g.shape[2] != self.get_numCols(): + # print('Error: input data dimensions do not match CT data dimensions') + # return False + if len(badPixelMap.shape) != 2 or numRows != badPixelMap.shape[0] or numCols != badPixelMap.shape[1]: print('Error: bad pixel map dimensions do not match CT data dimensions') return False if type(g) != type(badPixelMap): @@ -4725,11 +4784,11 @@ def badPixelCorrection(self, g, badPixelMap, windowSize=3): print('Error: projection data and bad pixel map must both be on the cpu or both be on the same gpu') return False - self.libprojectors.badPixelCorrection.argtypes = [ctypes.c_void_p, ctypes.c_void_p, ctypes.c_int, ctypes.c_bool] - return self.libprojectors.badPixelCorrection(g.data_ptr(), badPixelMap.data_ptr(), windowSize, g.is_cuda == False) + self.libprojectors.badPixelCorrection.argtypes = [ctypes.c_void_p, ctypes.c_int, ctypes.c_int, ctypes.c_int, ctypes.c_void_p, ctypes.c_int, ctypes.c_bool] + return self.libprojectors.badPixelCorrection(g.data_ptr(), numAngles, numRows, numCols, badPixelMap.data_ptr(), windowSize, g.is_cuda == False) else: - self.libprojectors.badPixelCorrection.argtypes = [ndpointer(ctypes.c_float, flags="C_CONTIGUOUS"), ndpointer(ctypes.c_float, flags="C_CONTIGUOUS"), ctypes.c_int, ctypes.c_bool] - return self.libprojectors.badPixelCorrection(g, badPixelMap, windowSize, True) + self.libprojectors.badPixelCorrection.argtypes = [ndpointer(ctypes.c_float, flags="C_CONTIGUOUS"), ctypes.c_int, ctypes.c_int, ctypes.c_int, ndpointer(ctypes.c_float, flags="C_CONTIGUOUS"), ctypes.c_int, ctypes.c_bool] + return self.libprojectors.badPixelCorrection(g, numAngles, numRows, numCols, badPixelMap, windowSize, True) def PriorBilateralFilter(self, f, spatialFWHM, intensityFWHM, prior=None): """Performs 3D Bilateral Filter (BLF) denoising method where the intensity distance is measured against a prior image @@ -5036,6 +5095,10 @@ def TV_denoise(self, f, delta, beta, numIter, p=1.2): # THIS SECTION OF FUNCTIONS SET AND GET VARIOUS PARAMETERS, INCLUDING THOSE THAT SET HOW LEAP IS TO BE RUN ################################################################################################################### ################################################################################################################### + def number_of_gpus(self): + self.libprojectors.number_of_gpus.restype = ctypes.c_int + return self.libprojectors.number_of_gpus() + def set_gpu(self, which): """Set which GPU to use, use -1 to do CPU calculations""" return self.set_GPU(which) @@ -5058,7 +5121,23 @@ def set_GPUs(self, listOfGPUs): listOfGPUs = np.ascontiguousarray(listOfGPUs, dtype=np.int32) self.set_model() return self.libprojectors.set_GPUs(listOfGPUs, int(listOfGPUs.size)) + + def get_gpus(self): + """Get the index of all of the GPUs being used""" + + self.libprojectors.get_gpus.argtypes = [ndpointer(ctypes.c_int, flags="C_CONTIGUOUS")] + self.libprojectors.get_gpus.restype = ctypes.c_int + self.set_model() + + if self.number_of_gpus() <= 0: + return [] + possible_gpu_list = -1*np.ones(self.number_of_gpus(), dtype=np.int32) + N = self.libprojectors.get_gpus(possible_gpu_list) + gpu_list = np.zeros(N, dtype=np.int32) + gpu_list[:] = possible_gpu_list[0:N] + return gpu_list + def get_gpu(self): """Get the index of the primary GPU that is being used""" return self.get_GPU() @@ -5089,6 +5168,15 @@ def set_diameterFOV(self, d): self.set_model() return self.libprojectors.set_rFOV(0.5*d) + def get_diameterFOV(self): + """Gets the diameterFOV parameter + + """ + #self.libprojectors.get_rFOV.argtypes = [ctypes.c_float] + self.libprojectors.get_rFOV.restype = ctypes.c_float + self.set_model() + return 2.0*self.libprojectors.get_rFOV() + def set_truncatedScan(self, aFlag): """Set the truncatedScan parameter @@ -5275,7 +5363,19 @@ def clear_attenuationMap(self): self.libprojectors.clear_attenuationMap.restype = ctypes.c_bool self.set_model() return self.libprojectors.clear_attenuationMap() + + def angles_are_defined(self): + #self.libprojectors.angles_are_defined.argtypes = [] + self.libprojectors.angles_are_defined.restype = ctypes.c_bool + self.set_model() + return self.libprojectors.angles_are_defined() + def angles_are_equispaced(self): + #self.libprojectors.angles_are_equispaced.argtypes = [] + self.libprojectors.angles_are_equispaced.restype = ctypes.c_bool + self.set_model() + return self.libprojectors.angles_are_equispaced() + def get_angles(self): """Get a numpy array of the projection angles""" if self.get_numAngles() > 0: @@ -5360,6 +5460,8 @@ def get_geometry(self): return 'FAN' elif geometryType == 3: return 'MODULAR' + elif geometryType == 4: + return 'CONE-PARALLEL' else: return 'UNKNOWN' @@ -6255,37 +6357,50 @@ def load_data(self, fileName, x=None, fileRange=None, rowRange=None, colRange=No print('To install PIL do: pip install imageio') return None if hasPIL == True: - currentWorkingDirectory = os.getcwd() - dataFolder, baseFileName = os.path.split(fileName) - if len(dataFolder) > 0: - os.chdir(dataFolder) - baseFileName, fileExtension = os.path.splitext(os.path.basename(baseFileName)) - templateFile = baseFileName + '_*' + fileExtension - fileList = glob.glob(os.path.split(templateFile)[1]) - if len(fileList) == 0: - os.chdir(currentWorkingDirectory) - print('file sequence does not exist') - return None - - if fileRange is not None: - # prune fileList - fileList_pruned = [] + + if fileName.find('*') != -1: + import imageio + fileList = glob.glob(fileName) + if fileRange is not None: + fileList = fileList[fileRange[0]:fileRange[1]+1] + ind = None + currentWorkingDirectory = None + else: + currentWorkingDirectory = os.getcwd() + dataFolder, baseFileName = os.path.split(fileName) + + if len(dataFolder) > 0: + os.chdir(dataFolder) + + sequence_separator = "_" + + baseFileName, fileExtension = os.path.splitext(os.path.basename(baseFileName)) + templateFile = baseFileName + '_*' + fileExtension + fileList = glob.glob(os.path.split(templateFile)[1]) + if len(fileList) == 0: + sequence_separator = "" + templateFile = baseFileName + '*' + fileExtension + fileList = glob.glob(os.path.split(templateFile)[1]) + if len(fileList) == 0: + os.chdir(currentWorkingDirectory) + print('file sequence does not exist') + return None + + if fileRange is not None: + # prune fileList + fileList_pruned = [] + for i in range(len(fileList)): + digit = int(fileList[i].replace(baseFileName+sequence_separator,'').replace(fileExt,'')) + if fileRange[0] <= digit and digit <= fileRange[1]: + fileList_pruned.append(fileList[i]) + fileList = fileList_pruned + + justDigits = [] for i in range(len(fileList)): - digit = int(fileList[i].replace(baseFileName+'_','').replace(fileExt,'')) - if fileRange[0] <= digit and digit <= fileRange[1]: - fileList_pruned.append(fileList[i]) - fileList = fileList_pruned - - justDigits = [] - for i in range(len(fileList)): - digitStr = fileList[i].replace(baseFileName+'_','').replace(fileExt,'') - justDigits.append(int(digitStr)) - ind = np.argsort(justDigits) + digitStr = fileList[i].replace(baseFileName+sequence_separator,'').replace(fileExt,'') + justDigits.append(int(digitStr)) + ind = np.argsort(justDigits) - - #print('found ' + str(len(fileList)) + ' images') - #print('reading first image: ' + str(fileList[0])) - #firstImg = np.array(Image.open(fileList[0])) firstImg = np.array(imageio.imread(fileList[0]), dtype=np.float32) numRows = firstImg.shape[0] @@ -6309,9 +6424,10 @@ def load_data(self, fileName, x=None, fileRange=None, rowRange=None, colRange=No x = np.zeros((len(fileList), numRows, numCols), dtype=np.float32) print('found ' + str(x.shape[0]) + ' images of size ' + str(firstImg.shape[0]) + ' x ' + str(firstImg.shape[1])) for i in range(len(fileList)): - #anImage = np.array(Image.open(fileList[ind[i]])) - #anImage = np.array(Image.open(fileList[ind[i]]).rotate(-0.5)) - anImage = np.array(imageio.imread(fileList[ind[i]]), dtype=np.float32) + if ind is None: + anImage = np.array(imageio.imread(fileList[i]), dtype=np.float32) + else: + anImage = np.array(imageio.imread(fileList[ind[i]]), dtype=np.float32) if rowRange is not None: if colRange is not None: x[i,:,:] = anImage[rowRange[0]:rowRange[1]+1,colRange[0]:colRange[1]+1] @@ -6322,7 +6438,8 @@ def load_data(self, fileName, x=None, fileRange=None, rowRange=None, colRange=No x[i,:,:] = anImage[:,colRange[0]:colRange[1]+1] else: x[i,:,:] = anImage[:,:] - os.chdir(currentWorkingDirectory) + if currentWorkingDirectory is not None: + os.chdir(currentWorkingDirectory) return x ''' try: diff --git a/src/parameters.cpp b/src/parameters.cpp index 52cbe34..f5b9a73 100644 --- a/src/parameters.cpp +++ b/src/parameters.cpp @@ -567,59 +567,68 @@ bool parameters::muSpecified() return false; } -bool parameters::allDefined() +bool parameters::allDefined(bool doPrint) { - return geometryDefined() & volumeDefined(); + return geometryDefined(doPrint) & volumeDefined(doPrint); } -bool parameters::geometryDefined() +bool parameters::geometryDefined(bool doPrint) { if (geometry != CONE && geometry != PARALLEL && geometry != FAN && geometry != MODULAR && geometry != CONE_PARALLEL) { - printf("Error: CT geometry type not defined!\n"); + if (doPrint) + printf("Error: CT geometry type not defined!\n"); return false; } if (numCols <= 0 || numRows <= 0 || numAngles <= 0 || pixelWidth <= 0.0 || pixelHeight <= 0.0) { - printf("Error: detector pixel sizes and number of data elements must be positive\n"); + if (doPrint) + printf("Error: detector pixel sizes and number of data elements must be positive\n"); return false; } if (geometry == MODULAR) { if (sourcePositions == NULL || moduleCenters == NULL || rowVectors == NULL || colVectors == NULL) { - printf("Error: sourcePositions, moduleCenters, rowVectors, and colVectors must be defined for modular-beam geometries\n"); + if (doPrint) + printf("Error: sourcePositions, moduleCenters, rowVectors, and colVectors must be defined for modular-beam geometries\n"); return false; } } else if (angularRange == 0.0 && phis == NULL) { - printf("Error: projection angles not defined\n"); + if (doPrint) + printf("Error: projection angles not defined\n"); return false; } if (geometry == CONE || geometry == FAN || geometry == CONE_PARALLEL) { if (sod <= 0.0 || sdd <= sod) { - printf("Error: sod and sdd must be positive for fan- and cone-beam geometries\n"); + if (doPrint) + printf("Error: sod and sdd must be positive for fan- and cone-beam geometries\n"); return false; } } if (detectorType == CURVED && geometry != CONE) { - printf("Error: curved detector only defined for cone-beam geometries\n"); + if (doPrint) + printf("Error: curved detector only defined for cone-beam geometries\n"); return false; } return true; } -bool parameters::volumeDefined() +bool parameters::volumeDefined(bool doPrint) { if (numX <= 0 || numY <= 0 || numZ <= 0 || voxelWidth <= 0.0 || voxelHeight <= 0.0 || volumeDimensionOrder < 0 || volumeDimensionOrder > 1) { - printf("numZ = %d voxelHeight = %f\n", numZ, voxelHeight); - printf("Error: volume voxel sizes and number of data elements must be positive\n"); + if (doPrint) + { + printf("numZ = %d voxelHeight = %f\n", numZ, voxelHeight); + printf("Error: volume voxel sizes and number of data elements must be positive\n"); + } return false; } else @@ -630,15 +639,19 @@ bool parameters::volumeDefined() { //voxelHeight = pixelHeight; //printf("Warning: for parallel and fan-beam data volume voxel height must equal detector pixel height, so forcing voxel height to match pixel height!\n"); - printf("Error: for parallel and fan-beam data volume voxel height must equal detector pixel height!\n"); - printf("Please modify either the detector pixel height or the voxel height so they match!\n"); + if (doPrint) + { + printf("Error: for parallel and fan-beam data volume voxel height must equal detector pixel height!\n"); + printf("Please modify either the detector pixel height or the voxel height so they match!\n"); + } return false; } if (numRows != numZ) { //voxelHeight = pixelHeight; //printf("Warning: for parallel and fan-beam data volume voxel height must equal detector pixel height, so forcing voxel height to match pixel height!\n"); - printf("Error: for parallel and fan-beam data numZ == numRows!\n"); + if (doPrint) + printf("Error: for parallel and fan-beam data numZ == numRows!\n"); return false; } offsetZ = 0.0; @@ -649,19 +662,22 @@ bool parameters::volumeDefined() if (modularbeamIsAxiallyAligned() == false || voxelSizeWorksForFastSF() == false) { voxelHeight = voxelWidth; - printf("Warning: for (non axially-aligned) modular-beam data volume voxel height must equal voxel width (voxels must be cubic), so forcing voxel height to match voxel width!\n"); + if (doPrint) + printf("Warning: for (non axially-aligned) modular-beam data volume voxel height must equal voxel width (voxels must be cubic), so forcing voxel height to match voxel width!\n"); } } if (isSymmetric()) { if (numX > 1) { - printf("Error: symmetric objects must specify numX = 1!\n"); + if (doPrint) + printf("Error: symmetric objects must specify numX = 1!\n"); return false; } if (numY % 2 == 1) { - printf("Error: symmetric objects must specify numY as even !\n"); + if (doPrint) + printf("Error: symmetric objects must specify numY as even !\n"); return false; } offsetX = 0.0; @@ -941,6 +957,14 @@ float* parameters::setToZero(float* data, uint64 N) return setToConstant(data, N, 0.0); } +bool parameters::angles_are_defined() +{ + if (phis != NULL || geometryDefined(false) == true) + return true; + else + return false; +} + bool parameters::set_angles() { if (phis != NULL) diff --git a/src/parameters.h b/src/parameters.h index 70c632b..c7aebac 100644 --- a/src/parameters.h +++ b/src/parameters.h @@ -63,21 +63,21 @@ class parameters * \brief returns whether all CT geometry and CT volume parameter values are defined and valid * \return returns true if all CT geometry and CT volume parameter values are defined and valid, false otherwise */ - bool allDefined(); + bool allDefined(bool doPrint = true); /** - * \fn allDefined + * \fn geometryDefined * \brief returns whether all CT geometry parameter values are defined and valid * \return returns true if all CT geometry parameter values are defined and valid, false otherwise */ - bool geometryDefined(); + bool geometryDefined(bool doPrint = true); /** - * \fn allDefined + * \fn volumeDefined * \brief returns whether all CT volume parameter values are defined and valid * \return returns true if all CT volume parameter values are defined and valid, false otherwise */ - bool volumeDefined(); + bool volumeDefined(bool doPrint = true); /** * \fn offsetScan_has_adequate_angular_range @@ -108,6 +108,12 @@ class parameters */ bool set_default_volume(float scale = 1.0); + /** + * \fn angles_are_defined + * \return returns true if the angles are defined, false otherwise + */ + bool angles_are_defined(); + /** * \fn set_angles * \brief sets phis, an array of the projection angles diff --git a/src/tomographic_models.cpp b/src/tomographic_models.cpp index 2edce23..fe2488b 100644 --- a/src/tomographic_models.cpp +++ b/src/tomographic_models.cpp @@ -1234,7 +1234,6 @@ bool tomographicModels::set_fanbeam(int numAngles, int numRows, int numCols, flo params.detectorType = parameters::FLAT; params.geometry = parameters::FAN; - params.detectorType = parameters::FLAT; params.sod = sod; params.sdd = sdd; params.pixelWidth = pixelWidth; @@ -1394,6 +1393,30 @@ int tomographicModels::get_volumeDimensionOrder() return params.volumeDimensionOrder; } +int tomographicModels::number_of_gpus() +{ +#ifndef __USE_CPU + return numberOfGPUs(); +#else + return 0; +#endif +} + +int tomographicModels::get_gpus(int* list_of_gpus) +{ +#ifndef __USE_CPU + int retVal = params.whichGPUs.size(); + if (list_of_gpus != NULL) + { + for (int i = 0; i < retVal; i++) + list_of_gpus[i] = params.whichGPUs[i]; + } + return retVal; +#else + return 0; +#endif +} + bool tomographicModels::set_GPU(int whichGPU) { #ifndef __USE_CPU @@ -2361,8 +2384,10 @@ bool tomographicModels::HighPassFilter(float* f, int N_1, int N_2, int N_3, floa #endif } -bool tomographicModels::badPixelCorrection(float* g, float* badPixelMap, int w, bool data_on_cpu) +bool tomographicModels::badPixelCorrection(float* g, int N_1, int N_2, int N_3, float* badPixelMap, int w, bool data_on_cpu) { + if (N_1 <= 0 || N_2 <= 0 || N_3 <= 0) + return false; #ifndef __USE_CPU if (params.whichGPU < 0) { @@ -2375,14 +2400,13 @@ bool tomographicModels::badPixelCorrection(float* g, float* badPixelMap, int w, else numProj = 0.0; - uint64 numElements = params.projectionData_numberOfElements(); + uint64 numElements = uint64(N_1) * uint64(N_2) * uint64(N_3); double dataSize = 4.0 * double(numElements) / pow(2.0, 30.0); //uint64 maxElements = 2147483646; if (data_on_cpu == true && (getAvailableGPUmemory(params.whichGPU) < dataSize || params.whichGPUs.size() > 1)) { // do chunking - int N_1 = params.numAngles; //int numSlices = std::min(N_1, maxSlicesForChunking); int numSlices = int(ceil(double(N_1) / double(int(params.whichGPUs.size())))); while (getAvailableGPUmemory(params.whichGPU) < double(numSlices) / double(N_1) * dataSize) @@ -2405,11 +2429,14 @@ bool tomographicModels::badPixelCorrection(float* g, float* badPixelMap, int w, int sliceStart = ichunk * numSlices; int sliceEnd = std::min(N_1 - 1, sliceStart + numSlices - 1); - float* g_chunk = &g[uint64(sliceStart) * uint64(params.numRows * params.numCols)]; + float* g_chunk = &g[uint64(sliceStart) * uint64(N_2 * N_3)]; int whichGPU = params.whichGPUs[omp_get_thread_num()]; parameters params_chunk = params; - params_chunk.removeProjections(sliceStart, sliceEnd); + params_chunk.numAngles = sliceEnd - sliceStart + 1; + params_chunk.numRows = N_2; + params_chunk.numCols = N_3; + //params_chunk.removeProjections(sliceStart, sliceEnd); //params_chunk.numAngles = sliceEnd - sliceStart + 1; params_chunk.whichGPU = whichGPU; @@ -2418,7 +2445,16 @@ bool tomographicModels::badPixelCorrection(float* g, float* badPixelMap, int w, return true; } else - return badPixelCorrection_gpu(g, ¶ms, badPixelMap, w, data_on_cpu); + { + parameters params_chunk = params; + params_chunk.numAngles = N_1; + params_chunk.numRows = N_2; + params_chunk.numCols = N_3; + //params_chunk.removeProjections(sliceStart, sliceEnd); + //params_chunk.numAngles = sliceEnd - sliceStart + 1; + + return badPixelCorrection_gpu(g, ¶ms_chunk, badPixelMap, w, data_on_cpu); + } #else printf("Error: GPU routines not included in this release!\n"); return false; diff --git a/src/tomographic_models.h b/src/tomographic_models.h index 90d9750..fbebd6e 100644 --- a/src/tomographic_models.h +++ b/src/tomographic_models.h @@ -13,7 +13,7 @@ #pragma once #endif -#define LEAP_VERSION "1.20" +#define LEAP_VERSION "1.21" /* #include @@ -405,6 +405,19 @@ class tomographicModels */ int get_volumeDimensionOrder(); + /** + * \fn number_of_gpus + * \return number of GPUs on the system + */ + int number_of_gpus(); + + /** + * \fn get_gpus + * \brief gets a list of all gpus being used + * \return number of gpus being used + */ + int get_gpus(int* list_of_gpus); + /** * \fn set_GPU * \brief sets the primary GPU index @@ -848,12 +861,15 @@ class tomographicModels * \fn badPixelCorrection * \brief applies a 2D median filter to the second two dimensions of a 3D array to a specified list of pixels only * \param[in] g: pointer to the 3D projection data (input and output) + * \param[in] N_1 number of samples in the first dimension + * \param[in] N_2 number of samples in the second dimension + * \param[in] N_3 number of samples in the third dimension * \param[in] badPixelMap: pointer to a 2D array of the bad pixels which are label as 1.0 * \param[in] w: the window size in each dimension (must be 3, 5, or 7) * \param[in] data_on_cpu: true if data (f) is on the cpu, false if it is on the gpu * \return true if operation was sucessful, false otherwise */ - bool badPixelCorrection(float* g, float* badPixelMap, int w, bool data_on_cpu); + bool badPixelCorrection(float* g, int N_1, int N_2, int N_3, float* badPixelMap, int w, bool data_on_cpu); /** * \fn BilateralFilter diff --git a/src/tomographic_models_c_interface.cpp b/src/tomographic_models_c_interface.cpp index 30a9bef..057938a 100644 --- a/src/tomographic_models_c_interface.cpp +++ b/src/tomographic_models_c_interface.cpp @@ -84,17 +84,17 @@ bool reset() bool all_defined() { - return tomo()->params.allDefined(); + return tomo()->params.allDefined(false); } bool ct_geometry_defined() { - return tomo()->params.geometryDefined(); + return tomo()->params.geometryDefined(false); } bool ct_volume_defined() { - return tomo()->params.volumeDefined(); + return tomo()->params.volumeDefined(false); } void set_log_error() @@ -497,6 +497,16 @@ int get_volumeDimensionOrder() return tomo()->get_volumeDimensionOrder(); } +int number_of_gpus() +{ + return tomo()->number_of_gpus(); +} + +int get_gpus(int* list_of_gpus) +{ + return tomo()->get_gpus(list_of_gpus); +} + bool set_GPU(int whichGPU) { return tomo()->set_GPU(whichGPU); @@ -537,6 +547,11 @@ bool set_rFOV(float rFOV_in) return tomo()->set_rFOV(rFOV_in); } +float get_rFOV() +{ + return tomo()->params.rFOV(); +} + bool set_offsetScan(bool aFlag) { return tomo()->params.set_offsetScan(aFlag); @@ -633,14 +648,26 @@ bool flipAttenuationMapSign(bool data_on_cpu) return tomo()->flipAttenuationMapSign(data_on_cpu); } +bool angles_are_defined() +{ + return tomo()->params.angles_are_defined(); +} + +bool angles_are_equispaced() +{ + return tomo()->params.anglesAreEquispaced(); +} + bool set_geometry(int which) { //CONE = 0, PARALLEL = 1, FAN = 2, MODULAR = 3 - if (which < parameters::CONE || which > parameters::MODULAR) + if (which < parameters::CONE || which > parameters::CONE_PARALLEL) return false; else { tomo()->params.geometry = which; + if (which != parameters::CONE) + tomo()->params.detectorType = parameters::FLAT; return true; } } @@ -934,9 +961,9 @@ bool MedianFilter2D(float* f, int N_1, int N_2, int N_3, float threshold, int w, return tomo()->MedianFilter2D(f, N_1, N_2, N_3, threshold, w, signalThreshold, data_on_cpu); } -bool badPixelCorrection(float* g, float* badPixelMap, int w, bool data_on_cpu) +bool badPixelCorrection(float* g, int N_1, int N_2, int N_3, float* badPixelMap, int w, bool data_on_cpu) { - return tomo()->badPixelCorrection(g, badPixelMap, w, data_on_cpu); + return tomo()->badPixelCorrection(g, N_1, N_2, N_3, badPixelMap, w, data_on_cpu); } bool BilateralFilter(float* f, int N_1, int N_2, int N_3, float spatialFWHM, float intensityFWHM, float scale, bool data_on_cpu) @@ -1117,6 +1144,8 @@ PYBIND11_MODULE(leapct, m) { m.def("set_voxelWidth", &set_voxelWidth, ""); m.def("set_voxelHeight", &set_voxelHeight, ""); m.def("set_geometry", &set_geometry, ""); + m.def("angles_are_defined", &angles_are_defined, ""); + m.def("angles_are_equispaced", &angles_are_equispaced, ""); m.def("projectConeBeam", &projectConeBeam, ""); m.def("backprojectConeBeam", &backprojectConeBeam, ""); m.def("projectFanBeam", &projectFanBeam, ""); @@ -1127,7 +1156,9 @@ PYBIND11_MODULE(leapct, m) { m.def("viewRangeNeededForBackprojection", &viewRangeNeededForBackprojection, ""); m.def("sliceRangeNeededForProjection", &sliceRangeNeededForProjection, ""); m.def("numRowsRequiredForBackprojectingSlab", &numRowsRequiredForBackprojectingSlab, ""); - m.def("set_GPU", &set_GPU, ""); + m.def("number_of_gpus", &number_of_gpus, ""); + m.def("get_gpus", &get_gpus, ""); + m.def("set_GPU", &set_GPU, ""); m.def("set_GPUs", &set_GPUs, ""); m.def("get_GPU", &get_GPU, ""); m.def("set_axisOfSymmetry", &set_axisOfSymmetry, ""); @@ -1135,6 +1166,7 @@ PYBIND11_MODULE(leapct, m) { m.def("clear_axisOfSymmetry", &clear_axisOfSymmetry, ""); m.def("set_projector", &set_projector, ""); m.def("set_rFOV", &set_rFOV, ""); + m.def("get_rFOV", &get_rFOV, ""); m.def("set_offsetScan", &set_offsetScan, ""); m.def("get_offsetScan", &get_offsetScan, ""); m.def("set_truncatedScan", &set_truncatedScan, ""); diff --git a/src/tomographic_models_c_interface.h b/src/tomographic_models_c_interface.h index 54f0b30..d570c92 100644 --- a/src/tomographic_models_c_interface.h +++ b/src/tomographic_models_c_interface.h @@ -124,6 +124,9 @@ extern "C" PROJECTOR_API bool set_voxelHeight(float H); extern "C" PROJECTOR_API bool set_geometry(int); +extern "C" PROJECTOR_API bool angles_are_defined(); +extern "C" PROJECTOR_API bool angles_are_equispaced(); + extern "C" PROJECTOR_API bool projectConeBeam(float* g, float* f, bool data_on_cpu, int numAngles, int numRows, int numCols, float pixelHeight, float pixelWidth, float centerRow, float centerCol, float* phis, float sod, float sdd, int numX, int numY, int numZ, float voxelWidth, float voxelHeight, float offsetX, float offsetY, float offsetZ); extern "C" PROJECTOR_API bool backprojectConeBeam(float* g, float* f, bool data_on_cpu, int numAngles, int numRows, int numCols, float pixelHeight, float pixelWidth, float centerRow, float centerCol, float* phis, float sod, float sdd, int numX, int numY, int numZ, float voxelWidth, float voxelHeight, float offsetX, float offsetY, float offsetZ); @@ -138,6 +141,8 @@ extern "C" PROJECTOR_API bool viewRangeNeededForBackprojection(int* viewsNeeded) extern "C" PROJECTOR_API bool sliceRangeNeededForProjection(int* slicesNeeded, bool doClip); extern "C" PROJECTOR_API int numRowsRequiredForBackprojectingSlab(int numSlicesPerChunk); +extern "C" PROJECTOR_API int number_of_gpus(); +extern "C" PROJECTOR_API int get_gpus(int* list_of_gpus); extern "C" PROJECTOR_API bool set_GPU(int whichGPU); extern "C" PROJECTOR_API bool set_GPUs(int* whichGPUs, int N); extern "C" PROJECTOR_API int get_GPU(); @@ -146,6 +151,7 @@ extern "C" PROJECTOR_API float get_axisOfSymmetry(); extern "C" PROJECTOR_API bool clear_axisOfSymmetry(); extern "C" PROJECTOR_API bool set_projector(int which); extern "C" PROJECTOR_API bool set_rFOV(float rFOV_in); +extern "C" PROJECTOR_API float get_rFOV(); extern "C" PROJECTOR_API bool set_offsetScan(bool); extern "C" PROJECTOR_API bool get_offsetScan(); extern "C" PROJECTOR_API bool set_truncatedScan(bool); @@ -220,7 +226,7 @@ extern "C" PROJECTOR_API bool HighPassFilter2D(float* f, int, int, int, float FW extern "C" PROJECTOR_API bool BlurFilter2D(float* f, int, int, int, float FWHM, bool data_on_cpu); extern "C" PROJECTOR_API bool HighPassFilter(float* f, int, int, int, float FWHM, bool data_on_cpu); extern "C" PROJECTOR_API bool MedianFilter2D(float* f, int, int, int, float threshold, int w, float signalThreshold, bool data_on_cpu); -extern "C" PROJECTOR_API bool badPixelCorrection(float* g, float* badPixelMap, int w, bool data_on_cpu); +extern "C" PROJECTOR_API bool badPixelCorrection(float* g, int, int, int, float* badPixelMap, int w, bool data_on_cpu); extern "C" PROJECTOR_API bool BilateralFilter(float* f, int N_1, int N_2, int N_3, float spatialFWHM, float intensityFWHM, float scale, bool data_on_cpu); extern "C" PROJECTOR_API bool PriorBilateralFilter(float* f, int N_1, int N_2, int N_3, float spatialFWHM, float intensityFWHM, float* prior, bool data_on_cpu); extern "C" PROJECTOR_API bool GuidedFilter(float* f, int N_1, int N_2, int N_3, int r, float epsilon, bool data_on_cpu);