Skip to content

Commit

Permalink
fdiff_direct_neumann use size_t for offsets (#1339)
Browse files Browse the repository at this point in the history 
offset1 and offset2 are approximately the size of volume so can overflow a 32bit int for a reasonable sized dataset.

Use long long for loops where signed values are needed (for openMP on MSVC), and size_t elsewhere.

Constify values that can be used in loop conditions

* fdiff_direct_periodic protect against overflow when long is 32bit

* fdiff_adjoint_neumann protect against overflow when long is 32bit

* fdiff_adjoint_periodic protect against overflow when long is 32bit

* FiniteDifferenceLibrary update wrapper functions to use size_t for array dimensions

* FiniteDifferenceLibrary use consistent check for z_dim

Co-authored-by: Gemma Fardell <47746591+gfardell@users.noreply.github.com>
  • Loading branch information
@samtygier-stfc @gfardell
samtygier-stfc and gfardell authored Sep 14, 2022
1 parent 5062f44 commit 0db53ce
Show file tree
Hide file tree
Showing 4 changed files with 68 additions and 63 deletions.
6 changes: 4 additions & 2 deletions CHANGELOG.md
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,7 @@
* 22.x.x
- added multiple colormaps to show2D
- added multiple colormaps to show2D
- Fix segfault in GradientOperator due to parameter overflows on windows systems

* 22.0.0
- Strongly convex functionality in TotalVariation and FGP_TV Functions
- Refactored KullbackLeibler function class. Fix bug on gradient method for SIRF objects
Expand Down Expand Up @@ -29,7 +31,7 @@
- TXRMDataReader is deprecated in favour of ZEISSDataReader
- GitHub Actions:
- Update to version 0.1.1 of lauramurgatroyd/build-sphinx-action for building the documentation - ensures docs are always built from cil master

* 21.4.1
- Removed prints from unittests and cleanup of unittest code.
- CMake:
Expand Down
18 changes: 9 additions & 9 deletions Wrappers/Python/cil/optimisation/operators/GradientOperator.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -240,10 +240,10 @@ def adjoint(self, x, out=None):
ctypes.POINTER(ctypes.c_float),
ctypes.POINTER(ctypes.c_float),
ctypes.POINTER(ctypes.c_float),
ctypes.c_long,
ctypes.c_long,
ctypes.c_long,
ctypes.c_long,
ctypes.c_size_t,
ctypes.c_size_t,
ctypes.c_size_t,
ctypes.c_size_t,
ctypes.c_int32,
ctypes.c_int32,
ctypes.c_int32]
Expand All @@ -252,18 +252,18 @@ def adjoint(self, x, out=None):
ctypes.POINTER(ctypes.c_float),
ctypes.POINTER(ctypes.c_float),
ctypes.POINTER(ctypes.c_float),
ctypes.c_long,
ctypes.c_long,
ctypes.c_long,
ctypes.c_size_t,
ctypes.c_size_t,
ctypes.c_size_t,
ctypes.c_int32,
ctypes.c_int32,
ctypes.c_int32]

cilacc.fdiff2D.argtypes = [ctypes.POINTER(ctypes.c_float),
ctypes.POINTER(ctypes.c_float),
ctypes.POINTER(ctypes.c_float),
ctypes.c_long,
ctypes.c_long,
ctypes.c_size_t,
ctypes.c_size_t,
ctypes.c_int32,
ctypes.c_int32,
ctypes.c_int32]
Expand Down
93 changes: 48 additions & 45 deletions src/Core/FiniteDifferenceLibrary.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,32 +15,33 @@ DLL_EXPORT int openMPtest(int nThreads)
return nThreads_running;
}

int fdiff_direct_neumann(const float *inimagefull, float *outimageXfull, float *outimageYfull, float *outimageZfull, float *outimageCfull, long nx, long ny, long nz, long nc)
int fdiff_direct_neumann(const float *inimagefull, float *outimageXfull, float *outimageYfull, float *outimageZfull, float *outimageCfull, size_t nx, size_t ny, size_t nz, size_t nc)
{
size_t volume = nx * ny * nz;
const size_t volume = nx * ny * nz;

const float *inimage = inimagefull;
float *outimageX = outimageXfull;
float *outimageY = outimageYfull;
float *outimageZ = outimageZfull;

int offset1 = (nz - 1) * nx * ny; //ind to beginning of last slice
int offset2 = offset1 + (ny - 1) * nx; //ind to beginning of last row
const size_t offset1 = (nz - 1) * nx * ny; //ind to beginning of last slice
const size_t offset2 = offset1 + (ny - 1) * nx; //ind to beginning of last row

long c;
// For openMP in MSVC loop variables must be signed
long long c;

int z_dim = nz > 1 ? 1: 0;
const bool z_dim = nz > 1;

for (c = 0; c < nc; c++)
{
#pragma omp parallel
{
long ind, k, j, i;
long long ind, k, j, i;
float pix0;
//run over all and then fix boundaries

#pragma omp for nowait
for (ind = 0; ind < nx * ny * (nz - 1); ind++)
for (ind = 0; ind < offset1; ind++)
{
pix0 = -inimage[ind];
outimageX[ind] = pix0 + inimage[ind + 1];
Expand Down Expand Up @@ -100,7 +101,7 @@ int fdiff_direct_neumann(const float *inimagefull, float *outimageXfull, float *
//now the rest of the channels
if (nc > 1)
{
long ind;
long long ind;

for (c = 0; c < nc - 1; c++)
{
Expand All @@ -123,29 +124,31 @@ int fdiff_direct_neumann(const float *inimagefull, float *outimageXfull, float *

return 0;
}
int fdiff_direct_periodic(const float *inimagefull, float *outimageXfull, float *outimageYfull, float *outimageZfull, float *outimageCfull, long nx, long ny, long nz, long nc)
int fdiff_direct_periodic(const float *inimagefull, float *outimageXfull, float *outimageYfull, float *outimageZfull, float *outimageCfull, size_t nx, size_t ny, size_t nz, size_t nc)
{
size_t volume = nx * ny * nz;
const size_t volume = nx * ny * nz;

const float *inimage = inimagefull;
float *outimageX = outimageXfull;
float *outimageY = outimageYfull;
float *outimageZ = outimageZfull;

int offset1 = (nz - 1) * nx * ny; //ind to beginning of last slice
int offset2 = offset1 + (ny - 1) * nx; //ind to beginning of last row
const size_t offset1 = (nz - 1) * nx * ny; //ind to beginning of last slice
const size_t offset2 = offset1 + (ny - 1) * nx; //ind to beginning of last row

long c;
const bool z_dim = nz > 1;

long long c;
for (c = 0; c < nc; c++)
{

#pragma omp parallel
{
long ind, k;
long long ind, k;
float pix0;
//run over all and then fix boundaries
#pragma omp for nowait
for (ind = 0; ind < nx * ny * (nz - 1); ind++)
for (ind = 0; ind < offset1; ind++)
{
pix0 = -inimage[ind];

Expand Down Expand Up @@ -177,8 +180,8 @@ int fdiff_direct_periodic(const float *inimagefull, float *outimageXfull, float
{
for (int i = 0; i < nx; i++)
{
int ind1 = (k * ny * nx);
int ind2 = ind1 + (ny - 1) * nx;
size_t ind1 = (k * ny * nx);
size_t ind2 = ind1 + (ny - 1) * nx;

outimageY[ind2 + i] = -inimage[ind2 + i] + inimage[ind1 + i];
}
Expand All @@ -189,14 +192,14 @@ int fdiff_direct_periodic(const float *inimagefull, float *outimageXfull, float
{
for (int j = 0; j < ny; j++)
{
int ind1 = k * ny * nx + j * nx;
int ind2 = ind1 + nx - 1;
size_t ind1 = k * ny * nx + j * nx;
size_t ind2 = ind1 + nx - 1;

outimageX[ind2] = -inimage[ind2] + inimage[ind1];
}
}

if (nz > 1)
if (z_dim)
{
#pragma omp for nowait
for (ind = 0; ind < ny * nx; ind++)
Expand All @@ -215,7 +218,7 @@ int fdiff_direct_periodic(const float *inimagefull, float *outimageXfull, float
//now the rest of the channels
if (nc > 1)
{
long ind;
long long ind;

for (c = 0; c < nc - 1; c++)
{
Expand All @@ -238,13 +241,13 @@ int fdiff_direct_periodic(const float *inimagefull, float *outimageXfull, float

return 0;
}
int fdiff_adjoint_neumann(float *outimagefull, const float *inimageXfull, const float *inimageYfull, const float *inimageZfull, const float *inimageCfull, long nx, long ny, long nz, long nc)
int fdiff_adjoint_neumann(float *outimagefull, const float *inimageXfull, const float *inimageYfull, const float *inimageZfull, const float *inimageCfull, size_t nx, size_t ny, size_t nz, size_t nc)
{
//runs over full data in x, y, z. then corrects elements for bounday conditions and sums
size_t volume = nx * ny * nz;
const size_t volume = nx * ny * nz;

//assumes nx and ny > 1
int z_dim = nz - 1;
const bool z_dim = nz > 1;

float *outimage = outimagefull;
const float *inimageX = inimageXfull;
Expand All @@ -260,20 +263,20 @@ int fdiff_adjoint_neumann(float *outimagefull, const float *inimageXfull, const
tempZ = (float *)malloc(volume * sizeof(float));
}

long c;
long long c;
for (c = 0; c < nc; c++) //just calculating x, y and z in each channel here
{
#pragma omp parallel
{
long ind, k;
long long ind, k;

#pragma omp for
for (ind = 1; ind < nx * ny * nz; ind++)
for (ind = 1; ind < volume; ind++)
{
tempX[ind] = -inimageX[ind] + inimageX[ind - 1];
}
#pragma omp for
for (ind = nx; ind < nx * ny * nz; ind++)
for (ind = nx; ind < volume; ind++)
{
tempY[ind] = -inimageY[ind] + inimageY[ind - nx];
}
Expand All @@ -282,7 +285,7 @@ int fdiff_adjoint_neumann(float *outimagefull, const float *inimageXfull, const
#pragma omp for
for (k = 0; k < nz; k++)
{
for (int j = 0; j < ny; j++)
for (long long j = 0; j < ny; j++)
{
tempX[k * ny * nx + j * nx] = -inimageX[k * ny * nx + j * nx];
tempX[k * ny * nx + j * nx + nx - 1] = inimageX[k * ny * nx + j * nx + nx - 2];
Expand All @@ -291,7 +294,7 @@ int fdiff_adjoint_neumann(float *outimagefull, const float *inimageXfull, const
#pragma omp for
for (k = 0; k < nz; k++)
{
for (int i = 0; i < nx; i++)
for (long long i = 0; i < nx; i++)
{
tempY[(k * ny * nx) + i] = -inimageY[(k * ny * nx) + i];
tempY[(k * ny * nx) + nx * (ny - 1) + i] = inimageY[(k * ny * nx) + nx * (ny - 2) + i];
Expand All @@ -301,7 +304,7 @@ int fdiff_adjoint_neumann(float *outimagefull, const float *inimageXfull, const
if (z_dim)
{
#pragma omp for
for (ind = nx * ny; ind < nx * ny * nz; ind++)
for (ind = nx * ny; ind < volume; ind++)
{
tempZ[ind] = -inimageZ[ind] + inimageZ[ind - nx * ny];
}
Expand Down Expand Up @@ -341,7 +344,7 @@ int fdiff_adjoint_neumann(float *outimagefull, const float *inimageXfull, const
// //now the rest of the channels
if (nc > 1)
{
long ind;
long long ind;

for (c = 1; c < nc - 1; c++)
{
Expand All @@ -363,13 +366,13 @@ int fdiff_adjoint_neumann(float *outimagefull, const float *inimageXfull, const

return 0;
}
int fdiff_adjoint_periodic(float *outimagefull, const float *inimageXfull, const float *inimageYfull, const float *inimageZfull, const float *inimageCfull, long nx, long ny, long nz, long nc)
int fdiff_adjoint_periodic(float *outimagefull, const float *inimageXfull, const float *inimageYfull, const float *inimageZfull, const float *inimageCfull, size_t nx, size_t ny, size_t nz, size_t nc)
{
//runs over full data in x, y, z. then correctects elements for bounday conditions and sums
size_t volume = nx * ny * nz;
const size_t volume = nx * ny * nz;

//assumes nx and ny > 1
int z_dim = nz - 1;
const bool z_dim = nz > 1;

float *outimage = outimagefull;
const float *inimageX = inimageXfull;
Expand All @@ -385,12 +388,12 @@ int fdiff_adjoint_periodic(float *outimagefull, const float *inimageXfull, const
tempZ = (float *)malloc(volume * sizeof(float));
}

long c;
long long c;
for (c = 0; c < nc; c++) //just calculating x, y and z in each channel here
{
#pragma omp parallel
{
long ind, k;
long long ind, k;

//run over all and then fix boundaries
#pragma omp for
Expand All @@ -408,15 +411,15 @@ int fdiff_adjoint_periodic(float *outimagefull, const float *inimageXfull, const
#pragma omp for
for (k = 0; k < nz; k++)
{
for (int i = 0; i < nx; i++)
for (long long i = 0; i < nx; i++)
{
tempY[(k * ny * nx) + i] = -inimageY[(k * ny * nx) + i] + inimageY[(k * ny * nx) + nx * (ny - 1) + i];
}
}
#pragma omp for
for (k = 0; k < nz; k++)
{
for (int j = 0; j < ny; j++)
for (long long j = 0; j < ny; j++)
{
tempX[k * ny * nx + j * nx] = -inimageX[k * ny * nx + j * nx] + inimageX[k * ny * nx + j * nx + nx - 1];
}
Expand All @@ -426,7 +429,7 @@ int fdiff_adjoint_periodic(float *outimagefull, const float *inimageXfull, const
{

#pragma omp for
for (ind = nx * ny; ind < nx * ny * nz; ind++)
for (ind = nx * ny; ind < volume; ind++)
{
tempZ[ind] = -inimageZ[ind] + inimageZ[ind - nx * ny];
}
Expand Down Expand Up @@ -466,7 +469,7 @@ int fdiff_adjoint_periodic(float *outimagefull, const float *inimageXfull, const
//now the rest of the channels
if (nc > 1)
{
long ind;
long long ind;

for (c = 1; c < nc; c++)
{
Expand All @@ -488,7 +491,7 @@ int fdiff_adjoint_periodic(float *outimagefull, const float *inimageXfull, const
return 0;
}

DLL_EXPORT int fdiff4D(float *imagefull, float *gradCfull, float *gradZfull, float *gradYfull, float *gradXfull, long nc, long nz, long ny, long nx, int boundary, int direction, int nThreads)
DLL_EXPORT int fdiff4D(float *imagefull, float *gradCfull, float *gradZfull, float *gradYfull, float *gradXfull, size_t nc, size_t nz, size_t ny, size_t nx, int boundary, int direction, int nThreads)
{
int nThreads_initial;
threads_setup(nThreads, &nThreads_initial);
Expand All @@ -511,7 +514,7 @@ DLL_EXPORT int fdiff4D(float *imagefull, float *gradCfull, float *gradZfull, flo
omp_set_num_threads(nThreads_initial);
return 0;
}
DLL_EXPORT int fdiff3D(float *imagefull, float *gradZfull, float *gradYfull, float *gradXfull, long nz, long ny, long nx, int boundary, int direction, int nThreads)
DLL_EXPORT int fdiff3D(float *imagefull, float *gradZfull, float *gradYfull, float *gradXfull, size_t nz, size_t ny, size_t nx, int boundary, int direction, int nThreads)
{
int nThreads_initial;
threads_setup(nThreads, &nThreads_initial);
Expand All @@ -534,7 +537,7 @@ DLL_EXPORT int fdiff3D(float *imagefull, float *gradZfull, float *gradYfull, flo
omp_set_num_threads(nThreads_initial);
return 0;
}
DLL_EXPORT int fdiff2D(float *imagefull, float *gradYfull, float *gradXfull, long ny, long nx, int boundary, int direction, int nThreads)
DLL_EXPORT int fdiff2D(float *imagefull, float *gradYfull, float *gradXfull, size_t ny, size_t nx, int boundary, int direction, int nThreads)
{
int nThreads_initial;
threads_setup(nThreads, &nThreads_initial);
Expand Down
14 changes: 7 additions & 7 deletions src/Core/include/FiniteDifferenceLibrary.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,19 +5,19 @@
#include "dll_export.h"
#include "utilities.h"

int fdiff_direct_neumann(const float *inimagefull, float *outimageXfull, float *outimageYfull, float *outimageZfull, float *outimageCfull, long nx, long ny, long nz, long nc);
int fdiff_direct_periodic(const float *inimagefull, float *outimageXfull, float *outimageYfull, float *outimageZfull, float *outimageCfull, long nx, long ny, long nz, long nc);
int fdiff_adjoint_neumann(float *outimagefull, const float *inimageXfull, const float *inimageYfull, const float *inimageZfull, const float *inimageCfull, long nx, long ny, long nz, long nc);
int fdiff_adjoint_periodic(float *outimagefull, const float *inimageXfull, const float *inimageYfull, const float *inimageZfull, const float *inimageCfull, long nx, long ny, long nz, long nc);
int fdiff_direct_neumann(const float *inimagefull, float *outimageXfull, float *outimageYfull, float *outimageZfull, float *outimageCfull, size_t nx, size_t ny, size_t nz, size_t nc);
int fdiff_direct_periodic(const float *inimagefull, float *outimageXfull, float *outimageYfull, float *outimageZfull, float *outimageCfull, size_t nx, size_t ny, size_t nz, size_t nc);
int fdiff_adjoint_neumann(float *outimagefull, const float *inimageXfull, const float *inimageYfull, const float *inimageZfull, const float *inimageCfull, size_t nx, size_t ny, size_t nz, size_t nc);
int fdiff_adjoint_periodic(float *outimagefull, const float *inimageXfull, const float *inimageYfull, const float *inimageZfull, const float *inimageCfull, size_t nx, size_t ny, size_t nz, size_t nc);

#ifdef __cplusplus
extern "C" {
#endif

DLL_EXPORT int openMPtest(int nThreads);
DLL_EXPORT int fdiff4D(float *imagefull, float *gradCfull, float *gradZfull, float *gradYfull, float *gradXfull, long nc, long nz, long ny, long nx, int boundary, int direction, int nThreads);
DLL_EXPORT int fdiff3D(float *imagefull, float *gradZfull, float *gradYfull, float *gradXfull, long nz, long ny, long nx, int boundary, int direction, int nThreads);
DLL_EXPORT int fdiff2D(float *imagefull, float *gradYfull, float *gradXfull, long ny, long nx, int boundary, int direction, int nThreads);
DLL_EXPORT int fdiff4D(float *imagefull, float *gradCfull, float *gradZfull, float *gradYfull, float *gradXfull, size_t nc, size_t nz, size_t ny, size_t nx, int boundary, int direction, int nThreads);
DLL_EXPORT int fdiff3D(float *imagefull, float *gradZfull, float *gradYfull, float *gradXfull, size_t nz, size_t ny, size_t nx, int boundary, int direction, int nThreads);
DLL_EXPORT int fdiff2D(float *imagefull, float *gradYfull, float *gradXfull, size_t ny, size_t nx, int boundary, int direction, int nThreads);

#ifdef __cplusplus
}
Expand Down

0 comments on commit 0db53ce

Please sign in to comment.