UsTK : Ultrasound ToolKit  version 2.0.1 under development (2023-12-07)
usPreScanToPostScan3DConverter.cpp
1 /****************************************************************************
2  *
3  * This file is part of the ustk software.
4  * Copyright (C) 2016 - 2017 by Inria. All rights reserved.
5  *
6  * This software is free software; you can redistribute it and/or
7  * modify it under the terms of the GNU General Public License
8  * ("GPL") version 2 as published by the Free Software Foundation.
9  * See the file LICENSE.txt at the root directory of this source
10  * distribution for additional information about the GNU GPL.
11  *
12  * For using ustk with software that can not be combined with the GNU
13  * GPL, please contact Inria about acquiring a ViSP Professional
14  * Edition License.
15  *
16  * This software was developed at:
17  * Inria Rennes - Bretagne Atlantique
18  * Campus Universitaire de Beaulieu
19  * 35042 Rennes Cedex
20  * France
21  *
22  * If you have questions regarding the use of this file, please contact
23  * Inria at ustk@inria.fr
24  *
25  * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
26  * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
27  *
28  * Authors:
29  * Jason Chevrie
30  * Marc Pouliquen
31  *
32  *****************************************************************************/
33 
34 #include <visp3/ustk_core/usPreScanToPostScan3DConverter.h>
35 
36 #ifdef VISP_HAVE_OPENMP
37 #include <omp.h>
38 #endif
39 
40 #ifdef USTK_HAVE_CUDA
41 extern void GPUDirectConversionWrapper(unsigned char *dataPost, const unsigned char *dataPre, unsigned int m_nbX, unsigned int m_nbY, unsigned int m_nbZ, int X, int Y, int Z, double m_resolution, double xmax, double ymin, double zmax, unsigned int frameNumber, unsigned int scanLineNumber, double transducerRadius, double motorRadius, double scanLinePitch, double axialResolution, double framePitch, bool sweepInZdirection);
42 #endif
43 
47 usPreScanToPostScan3DConverter::usPreScanToPostScan3DConverter()
48  : m_converterOptimizationMethod(SINGLE_THREAD_REDUCED_LOOKUP_TABLE),
49  m_conversionOptimizationMethodUsedAtInit(SINGLE_THREAD_DIRECT_CONVERSION), m_GPULookupTables{NULL, NULL}, m_GPULookupTablesSize{0,0}, m_VpreScan(), m_downSamplingFactor(1),
50  m_resolution(), m_SweepInZdirection(true), m_initDone(false)
51 {
52 }
53 
59 usPreScanToPostScan3DConverter::usPreScanToPostScan3DConverter(const usImagePreScan3D<unsigned char> &preScanImage,
60  double down)
61  : m_VpreScan(), m_downSamplingFactor(1), m_resolution(), m_SweepInZdirection(true), m_initDone(false)
62 {
63  this->init(preScanImage, down);
64 }
65 
71 void usPreScanToPostScan3DConverter::init(const usImagePreScan3D<unsigned char> &preScanImage, double down)
72 {
73  if (!preScanImage.isTransducerConvex() || !(preScanImage.getMotorType() == usMotorSettings::TiltingMotor))
74  throw(vpException(vpException::functionNotImplementedError,
75  "3D scan-conversion available only for convex transducer and tilting motor"));
76 
77  if (down <= 0)
78  throw(vpException(vpException::badValue, "downsampling factor should b positive"));
79 
80  // compare pre-scan image parameters, to avoid recomputing all the init process if parameters are the same
81  if (((usMotorSettings)m_VpreScan) == ((usMotorSettings)preScanImage) &&
82  ((usImagePreScanSettings)m_VpreScan) == ((usImagePreScanSettings)preScanImage) &&
83  m_VpreScan.getWidth() == preScanImage.getWidth() && m_VpreScan.getHeight() == preScanImage.getHeight() &&
84  m_VpreScan.getNumberOfFrames() == preScanImage.getNumberOfFrames() &&
85  m_resolution == down * m_VpreScan.getAxialResolution() && m_downSamplingFactor == down &&
86  m_converterOptimizationMethod == m_conversionOptimizationMethodUsedAtInit) {
87  m_VpreScan = preScanImage; // update image content
88  return;
89  }
90 
91  m_VpreScan = preScanImage;
92  m_resolution = down * m_VpreScan.getAxialResolution();
93  m_downSamplingFactor = down;
94 
95  int X = m_VpreScan.getWidth();
96  int Y = m_VpreScan.getHeight();
97  int Z = m_VpreScan.getNumberOfFrames();
98 
99  double xmax;
100  double ymin;
101  double ymax;
102  double zmax;
103 
104  usPreScanToPostScan3DConverter::convertPreScanCoordToPostScanCoord(0.0, X, Z, &ymin, NULL, NULL);
105  usPreScanToPostScan3DConverter::convertPreScanCoordToPostScanCoord((double)Y, X / 2.0, Z / 2.0, &ymax, NULL, NULL);
106  usPreScanToPostScan3DConverter::convertPreScanCoordToPostScanCoord((double)Y, (double)X, Z / 2.0, NULL, &xmax, NULL);
107  usPreScanToPostScan3DConverter::convertPreScanCoordToPostScanCoord((double)Y, X / 2.0, Z, NULL, NULL, &zmax);
108 
109  m_nbX = (unsigned int)ceil(2 * xmax / m_resolution);
110  m_nbY = (unsigned int)ceil((ymax - ymin) / m_resolution);
111  m_nbZ = (unsigned int)ceil(2 * zmax / m_resolution);
112 
113  unsigned int nbXY = m_nbX * m_nbY;
114  unsigned int XY = X * Y;
115 
116  // empty the lookup tables
117  if (m_lookupTables[0].size() > 0)
118  std::vector<usVoxelWeightAndIndex>().swap(m_lookupTables[0]);
119  if (m_lookupTables[1].size() > 0)
120  std::vector<usVoxelWeightAndIndex>().swap(m_lookupTables[1]);
121  if (m_reducedLookupTables[0].size() > 0)
122  std::vector<usVoxelWeightAndIndexReducedMemory>().swap(m_reducedLookupTables[0]);
123  if (m_reducedLookupTables[1].size() > 0)
124  std::vector<usVoxelWeightAndIndexReducedMemory>().swap(m_reducedLookupTables[1]);
125 #ifdef USTK_HAVE_CUDA
126  this->GPUFreeLookupTables();
127 #endif
128 
129  switch (m_converterOptimizationMethod) {
130  case SINGLE_THREAD_DIRECT_CONVERSION: {
131  break;
132  }
133  case MULTI_THREAD_DIRECT_CONVERSION: {
134  break;
135  }
136  case GPU_DIRECT_CONVERSION: {
137  break;
138  }
139  case SINGLE_THREAD_FULL_LOOKUP_TABLE: {
140  try {
141  long int LUTmaxSize = (long int)m_nbX * (long int)m_nbY * (long int)m_nbZ;
142  // reserve to avoid reallocation during the LUT filling
143  m_lookupTables[0].reserve(LUTmaxSize);
144  m_lookupTables[1].reserve(LUTmaxSize);
145  } catch (std::exception &e) {
146  throw vpException(vpException::ioError, "usPreScanToPostScan3DConverter::init: using method "
147  "SINGLE_THREAD_FULL_LOOKUP_TABLE leads to %s \n Use another optimization "
148  "method or downsample the volume",
149  e.what());
150  }
151  for (unsigned int sweepingDirection = 0; sweepingDirection < 2; sweepingDirection++) {
152  for (unsigned int x = 0; x < m_nbX; x++) {
153  double xx = m_resolution * x - xmax;
154  for (unsigned int y = 0; y < m_nbY; y++) {
155  double yy = ymin + m_resolution * y;
156  for (unsigned int z = 0; z < m_nbZ; z++) {
157  double zz = m_resolution * z - zmax;
158  double i, j, k;
159  usPreScanToPostScan3DConverter::convertPostScanCoordToPreScanCoord(yy, xx, zz, &j, &i, &k,
160  sweepingDirection == 0);
161 
162  double ii = floor(i);
163  double jj = floor(j);
164  double kk = floor(k);
165 
166  if (ii >= 0 && jj >= 0 && kk >= 0 && ii + 1 < X && jj + 1 < Y && kk + 1 < Z) {
167  usVoxelWeightAndIndex m;
168 
169  m.m_outputIndex = x + m_nbX * y + nbXY * z;
170 
171  double u = i - ii;
172  double v = j - jj;
173  double w = k - kk;
174  double u1 = 1 - u;
175  double v1 = 1 - v;
176  double w1 = 1 - w;
177 
178  double v1w1 = v1 * w1;
179  double vw1 = v * w1;
180  double v1w = v1 * w;
181  double vw = v * w;
182 
183  m.m_W[0] = u1 * v1w1;
184  m.m_W[1] = u * v1w1;
185  m.m_W[2] = u1 * vw1;
186  m.m_W[3] = u * vw1;
187  m.m_W[4] = u1 * v1w;
188  m.m_W[5] = u * v1w;
189  m.m_W[6] = u1 * vw;
190  m.m_W[7] = u * vw;
191 
192  double Xjj = X * jj;
193  double Xjj1 = X * (jj + 1);
194  double XYKK = XY * kk;
195  double XYKK1 = XY * (kk + 1);
196 
197  m.m_inputIndex[0] = (unsigned int)(ii + Xjj + XYKK);
198  m.m_inputIndex[1] = (unsigned int)(ii + 1 + Xjj + XYKK);
199  m.m_inputIndex[2] = (unsigned int)(ii + Xjj1 + XYKK);
200  m.m_inputIndex[3] = (unsigned int)(ii + 1 + Xjj1 + XYKK);
201  m.m_inputIndex[4] = (unsigned int)(ii + Xjj + XYKK1);
202  m.m_inputIndex[5] = (unsigned int)(ii + 1 + Xjj + XYKK1);
203  m.m_inputIndex[6] = (unsigned int)(ii + Xjj1 + XYKK1);
204  m.m_inputIndex[7] = (unsigned int)(ii + 1 + Xjj1 + XYKK1);
205 
206  m_lookupTables[sweepingDirection].push_back(m);
207  }
208  }
209  }
210  }
211  }
212  std::cout << "LUT 1 size (bytes) : " << sizeof(usVoxelWeightAndIndex) * m_lookupTables[0].size() << std::endl;
213  std::cout << "LUT 2 size (bytes) : " << sizeof(usVoxelWeightAndIndex) * m_lookupTables[1].size() << std::endl;
214  break;
215  }
216  case MULTI_THREAD_FULL_LOOKUP_TABLE: {
217  try {
218  long int LUTmaxSize = (long int)m_nbX * (long int)m_nbY * (long int)m_nbZ;
219  // reserve to avoid reallocation during the LUT filling
220  m_lookupTables[0].reserve(LUTmaxSize);
221  m_lookupTables[1].reserve(LUTmaxSize);
222  } catch (std::exception &e) {
223  throw vpException(vpException::ioError, "usPreScanToPostScan3DConverter::init: using method "
224  "MULTI_THREAD_FULL_LOOKUP_TABLE leads to %s \n Use another optimization "
225  "method or downsample the volume",
226  e.what());
227  }
228 
229  for (unsigned int sweepingDirection = 0 ; sweepingDirection < 2 ; sweepingDirection++) {
230 #ifdef VISP_HAVE_OPENMP
231 #pragma omp parallel for
232 #endif
233  for (int x = 0; x < (int)m_nbX; x++) {
234  double xx = m_resolution * x - xmax;
235  for (unsigned int y = 0; y < m_nbY; y++) {
236  double yy = ymin + m_resolution * y;
237  for (unsigned int z = 0; z < m_nbZ; z++) {
238  double zz = m_resolution * z - zmax;
239  double i, j, k;
240  usPreScanToPostScan3DConverter::convertPostScanCoordToPreScanCoord(yy, xx, zz, &j, &i, &k,
241  sweepingDirection == 0);
242 
243  double ii = floor(i);
244  double jj = floor(j);
245  double kk = floor(k);
246 
247  if (ii >= 0 && jj >= 0 && kk >= 0 && ii + 1 < X && jj + 1 < Y && kk + 1 < Z) {
248  usVoxelWeightAndIndex m;
249 
250  m.m_outputIndex = x + m_nbX * y + nbXY * z;
251 
252  double u = i - ii;
253  double v = j - jj;
254  double w = k - kk;
255  double u1 = 1 - u;
256  double v1 = 1 - v;
257  double w1 = 1 - w;
258 
259  double v1w1 = v1 * w1;
260  double vw1 = v * w1;
261  double v1w = v1 * w;
262  double vw = v * w;
263 
264  m.m_W[0] = u1 * v1w1;
265  m.m_W[1] = u * v1w1;
266  m.m_W[2] = u1 * vw1;
267  m.m_W[3] = u * vw1;
268  m.m_W[4] = u1 * v1w;
269  m.m_W[5] = u * v1w;
270  m.m_W[6] = u1 * vw;
271  m.m_W[7] = u * vw;
272 
273  double Xjj = X * jj;
274  double Xjj1 = X * (jj + 1);
275  double XYKK = XY * kk;
276  double XYKK1 = XY * (kk + 1);
277 
278  m.m_inputIndex[0] = (unsigned int)(ii + Xjj + XYKK);
279  m.m_inputIndex[1] = (unsigned int)(ii + 1 + Xjj + XYKK);
280  m.m_inputIndex[2] = (unsigned int)(ii + Xjj1 + XYKK);
281  m.m_inputIndex[3] = (unsigned int)(ii + 1 + Xjj1 + XYKK);
282  m.m_inputIndex[4] = (unsigned int)(ii + Xjj + XYKK1);
283  m.m_inputIndex[5] = (unsigned int)(ii + 1 + Xjj + XYKK1);
284  m.m_inputIndex[6] = (unsigned int)(ii + Xjj1 + XYKK1);
285  m.m_inputIndex[7] = (unsigned int)(ii + 1 + Xjj1 + XYKK1);
286 #ifdef VISP_HAVE_OPENMP
287 #pragma omp critical
288 #endif
289  m_lookupTables[sweepingDirection].push_back(m);
290  }
291  }
292  }
293  }
294  }
295  std::cout << "LUT 1 size (bytes) : " << sizeof(usVoxelWeightAndIndex) * m_lookupTables[0].size() << std::endl;
296  std::cout << "LUT 2 size (bytes) : " << sizeof(usVoxelWeightAndIndex) * m_lookupTables[1].size() << std::endl;
297  break;
298  }
299  case GPU_FULL_LOOKUP_TABLE: {
300 #ifdef USTK_HAVE_CUDA
301  this->GPUAllocateFullLookupTables();
302  this->GPUFillFullLookupTables();
303 std::cout << "LUT 1 size (bytes) : " << sizeof(usVoxelWeightAndIndex) * m_GPULookupTablesSize[0] << std::endl;
304  std::cout << "LUT 2 size (bytes) : " << sizeof(usVoxelWeightAndIndex) * m_GPULookupTablesSize[1] << std::endl;
305 #else
306  throw vpException(
307  vpException::notImplementedError,
308  "usPreScanToPostScan3DConverter::init: using method GPU_FULL_LOOKUP_TABLE is not implemented yet");
309 #endif
310  break;
311  }
312  case SINGLE_THREAD_REDUCED_LOOKUP_TABLE: {
313  try {
314  long int LUTmaxSize = (long int)m_nbX * (long int)m_nbY * (long int)m_nbZ;
315  // reserve to avoid reallocation during the LUT filling
316  m_reducedLookupTables[0].reserve(LUTmaxSize);
317  m_reducedLookupTables[1].reserve(LUTmaxSize);
318  } catch (std::exception &e) {
319  throw vpException(vpException::ioError, "usPreScanToPostScan3DConverter::init: using method "
320  "SINGLE_THREAD_REDUCED_LOOKUP_TABLE leads to %s \n Use another "
321  "optimization method or downsample the volume",
322  e.what());
323  }
324  for (unsigned int sweepingDirection = 0; sweepingDirection < 2; sweepingDirection++) {
325  for (unsigned int x = 0; x < m_nbX; x++) {
326  double xx = m_resolution * x - xmax;
327  for (unsigned int y = 0; y < m_nbY; y++) {
328  double yy = ymin + m_resolution * y;
329  for (unsigned int z = 0; z < m_nbZ; z++) {
330  double zz = m_resolution * z - zmax;
331  double i, j, k;
332  usPreScanToPostScan3DConverter::convertPostScanCoordToPreScanCoord(yy, xx, zz, &j, &i, &k,
333  sweepingDirection == 0);
334 
335  double ii = floor(i);
336  double jj = floor(j);
337  double kk = floor(k);
338 
339  if (ii >= 0 && jj >= 0 && kk >= 0 && ii + 1 < X && jj + 1 < Y && kk + 1 < Z) {
340  usVoxelWeightAndIndexReducedMemory m;
341 
342  m.m_outputIndex = x + m_nbX * y + nbXY * z;
343 
344  m.m_W[0] = i - ii;
345  m.m_W[1] = j - jj;
346  m.m_W[2] = k - kk;
347 
348  m.m_inputIndex = (unsigned int)(ii + X * jj + XY * kk);
349 
350  m_reducedLookupTables[sweepingDirection].push_back(m);
351  }
352  }
353  }
354  }
355  }
356  std::cout << "LUT 1 size (bytes) : " << sizeof(usVoxelWeightAndIndexReducedMemory) * m_reducedLookupTables[0].size()
357  << std::endl;
358  std::cout << "LUT 2 size (bytes) : " << sizeof(usVoxelWeightAndIndexReducedMemory) * m_reducedLookupTables[1].size()
359  << std::endl;
360  break;
361  }
362  case MULTI_THREAD_REDUCED_LOOKUP_TABLE: {
363  try {
364  long int LUTmaxSize = (long int)m_nbX * (long int)m_nbY * (long int)m_nbZ;
365  // reserve to avoid reallocation during the LUT filling
366  m_reducedLookupTables[0].reserve(LUTmaxSize);
367  m_reducedLookupTables[1].reserve(LUTmaxSize);
368  } catch (std::exception &e) {
369  throw vpException(vpException::ioError, "usPreScanToPostScan3DConverter::init: using method "
370  "MULTI_THREAD_REDUCED_LOOKUP_TABLE leads to %s \n Use another "
371  "optimization method or downsample the volume",
372  e.what());
373  }
374  for (unsigned int sweepingDirection = 0; sweepingDirection < 2; sweepingDirection++) {
375 #ifdef VISP_HAVE_OPENMP
376 #pragma omp parallel for
377 #endif
378  for (int x = 0; x < (int)m_nbX; x++) {
379  double xx = m_resolution * x - xmax;
380  for (unsigned int y = 0; y < m_nbY; y++) {
381  double yy = ymin + m_resolution * y;
382  for (unsigned int z = 0; z < m_nbZ; z++) {
383  double zz = m_resolution * z - zmax;
384  double i, j, k;
385  usPreScanToPostScan3DConverter::convertPostScanCoordToPreScanCoord(yy, xx, zz, &j, &i, &k,
386  sweepingDirection == 0);
387 
388  double ii = floor(i);
389  double jj = floor(j);
390  double kk = floor(k);
391 
392  if (ii >= 0 && jj >= 0 && kk >= 0 && ii + 1 < X && jj + 1 < Y && kk + 1 < Z) {
393  usVoxelWeightAndIndexReducedMemory m;
394 
395  m.m_outputIndex = x + m_nbX * y + nbXY * z;
396 
397  m.m_W[0] = i - ii;
398  m.m_W[1] = j - jj;
399  m.m_W[2] = k - kk;
400 
401  m.m_inputIndex = (unsigned int)(ii + X * jj + XY * kk);
402 #ifdef VISP_HAVE_OPENMP
403 #pragma omp critical
404 #endif
405  m_reducedLookupTables[sweepingDirection].push_back(m);
406  }
407  }
408  }
409  }
410  }
411  std::cout << "LUT 1 size (bytes) : " << sizeof(usVoxelWeightAndIndexReducedMemory) * m_reducedLookupTables[0].size()
412  << std::endl;
413  std::cout << "LUT 2 size (bytes) : " << sizeof(usVoxelWeightAndIndexReducedMemory) * m_reducedLookupTables[1].size()
414  << std::endl;
415  break;
416  }
417  case GPU_REDUCED_LOOKUP_TABLE: {
418 #ifdef USTK_HAVE_CUDA
419  this->GPUAllocateReducedLookupTables();
420  this->GPUFillReducedLookupTables();
421  std::cout << "LUT 1 size (bytes) : " << sizeof(usVoxelWeightAndIndexReducedMemory) * m_GPULookupTablesSize[0] << std::endl;
422  std::cout << "LUT 2 size (bytes) : " << sizeof(usVoxelWeightAndIndexReducedMemory) * m_GPULookupTablesSize[1] << std::endl;
423 #else
424  throw vpException(
425  vpException::notImplementedError,
426  "usPreScanToPostScan3DConverter::init: using method GPU_REDUCED_LOOKUP_TABLE is not implemented yet");
427 #endif
428  break;
429  }
430  }
431 
432  m_conversionOptimizationMethodUsedAtInit = m_converterOptimizationMethod;
433  m_initDone = true;
434 }
435 
439 usPreScanToPostScan3DConverter::~usPreScanToPostScan3DConverter() {}
440 
446 void usPreScanToPostScan3DConverter::convert(usImagePostScan3D<unsigned char> &postScanImage,
447  const usImagePreScan3D<unsigned char> &preScanImage)
448 {
449  init(preScanImage, m_downSamplingFactor);
450 
451  postScanImage.resize(m_nbY, m_nbX, m_nbZ);
452  postScanImage.initData(0);
453  (usTransducerSettings&)postScanImage = (const usTransducerSettings&)preScanImage;
454  (usMotorSettings&)postScanImage = (const usMotorSettings&)preScanImage;
455 
456  unsigned char *dataPost = postScanImage.getData();
457  const unsigned char *dataPre = preScanImage.getConstData();
458 
459  switch (m_converterOptimizationMethod) {
460  case SINGLE_THREAD_DIRECT_CONVERSION: {
461  int X = m_VpreScan.getWidth();
462  int Y = m_VpreScan.getHeight();
463  int Z = m_VpreScan.getNumberOfFrames();
464 
465  double xmax;
466  double ymin;
467  double ymax;
468  double zmax;
469 
470  usPreScanToPostScan3DConverter::convertPreScanCoordToPostScanCoord(0.0, X, Z, &ymin, NULL, NULL);
471  usPreScanToPostScan3DConverter::convertPreScanCoordToPostScanCoord((double)Y, X / 2.0, Z / 2.0, &ymax, NULL, NULL);
472  usPreScanToPostScan3DConverter::convertPreScanCoordToPostScanCoord((double)Y, (double)X, Z / 2.0, NULL, &xmax,
473  NULL);
474  usPreScanToPostScan3DConverter::convertPreScanCoordToPostScanCoord((double)Y, X / 2.0, Z, NULL, NULL, &zmax);
475 
476  unsigned int nbXY = m_nbX * m_nbY;
477  unsigned int XY = X * Y;
478 
479  for (unsigned int x = 0; x < m_nbX; x++) {
480  double xx = m_resolution * x - xmax;
481  for (unsigned int y = 0; y < m_nbY; y++) {
482  double yy = ymin + m_resolution * y;
483  for (unsigned int z = 0; z < m_nbZ; z++) {
484  double zz = m_resolution * z - zmax;
485  double i, j, k;
486  usPreScanToPostScan3DConverter::convertPostScanCoordToPreScanCoord(yy, xx, zz, &j, &i, &k,
487  m_SweepInZdirection);
488 
489  double ii = floor(i);
490  double jj = floor(j);
491  double kk = floor(k);
492 
493  if (ii >= 0 && jj >= 0 && kk >= 0 && ii + 1 < X && jj + 1 < Y && kk + 1 < Z) {
494 
495  double u = i - ii;
496  double v = j - jj;
497  double w = k - kk;
498  double u1 = 1 - u;
499  double v1 = 1 - v;
500  double w1 = 1 - w;
501 
502  double v1w1 = v1 * w1;
503  double vw1 = v * w1;
504  double v1w = v1 * w;
505  double vw = v * w;
506 
507  double W[8] = {u1 * v1w1, u * v1w1, u1 * vw1, u * vw1, u1 * v1w, u * v1w, u1 * vw, u * vw};
508 
509  double Xjj = X * jj;
510  double Xjj1 = X * (jj + 1);
511  double XYKK = XY * kk;
512  double XYKK1 = XY * (kk + 1);
513 
514  unsigned int index[8] = {(unsigned int)(ii + Xjj + XYKK), (unsigned int)(ii + 1 + Xjj + XYKK),
515  (unsigned int)(ii + Xjj1 + XYKK), (unsigned int)(ii + 1 + Xjj1 + XYKK),
516  (unsigned int)(ii + Xjj + XYKK1), (unsigned int)(ii + 1 + Xjj + XYKK1),
517  (unsigned int)(ii + Xjj1 + XYKK1), (unsigned int)(ii + 1 + Xjj1 + XYKK1)};
518 
519  double val = 0;
520  for (int n = 0; n < 8; n++)
521  val += W[n] * dataPre[index[n]];
522 
523  dataPost[x + m_nbX * y + nbXY * z] = (unsigned char)val;
524  }
525  }
526  }
527  }
528  break;
529  }
530  case MULTI_THREAD_DIRECT_CONVERSION: {
531  int X = m_VpreScan.getWidth();
532  int Y = m_VpreScan.getHeight();
533  int Z = m_VpreScan.getNumberOfFrames();
534 
535  double xmax;
536  double ymin;
537  double ymax;
538  double zmax;
539 
540  usPreScanToPostScan3DConverter::convertPreScanCoordToPostScanCoord(0.0, X, Z, &ymin, NULL, NULL);
541  usPreScanToPostScan3DConverter::convertPreScanCoordToPostScanCoord((double)Y, X / 2.0, Z / 2.0, &ymax, NULL, NULL);
542  usPreScanToPostScan3DConverter::convertPreScanCoordToPostScanCoord((double)Y, (double)X, Z / 2.0, NULL, &xmax,
543  NULL);
544  usPreScanToPostScan3DConverter::convertPreScanCoordToPostScanCoord((double)Y, X / 2.0, Z, NULL, NULL, &zmax);
545 
546  unsigned int nbXY = m_nbX * m_nbY;
547  unsigned int XY = X * Y;
548 
549 #ifdef VISP_HAVE_OPENMP
550 #pragma omp parallel for
551 #endif
552  for (int x = 0; x < (int)m_nbX; x++) {
553  double xx = m_resolution * x - xmax;
554  for (unsigned int y = 0; y < m_nbY; y++) {
555  double yy = ymin + m_resolution * y;
556  for (unsigned int z = 0; z < m_nbZ; z++) {
557  double zz = m_resolution * z - zmax;
558  double i, j, k;
559  usPreScanToPostScan3DConverter::convertPostScanCoordToPreScanCoord(yy, xx, zz, &j, &i, &k,
560  m_SweepInZdirection);
561 
562  double ii = floor(i);
563  double jj = floor(j);
564  double kk = floor(k);
565 
566  if (ii >= 0 && jj >= 0 && kk >= 0 && ii + 1 < X && jj + 1 < Y && kk + 1 < Z) {
567 
568  double u = i - ii;
569  double v = j - jj;
570  double w = k - kk;
571  double u1 = 1 - u;
572  double v1 = 1 - v;
573  double w1 = 1 - w;
574 
575  double v1w1 = v1 * w1;
576  double vw1 = v * w1;
577  double v1w = v1 * w;
578  double vw = v * w;
579 
580  double W[8] = {u1 * v1w1, u * v1w1, u1 * vw1, u * vw1, u1 * v1w, u * v1w, u1 * vw, u * vw};
581 
582  double Xjj = X * jj;
583  double Xjj1 = X * (jj + 1);
584  double XYKK = XY * kk;
585  double XYKK1 = XY * (kk + 1);
586 
587  unsigned int index[8] = {(unsigned int)(ii + Xjj + XYKK), (unsigned int)(ii + 1 + Xjj + XYKK),
588  (unsigned int)(ii + Xjj1 + XYKK), (unsigned int)(ii + 1 + Xjj1 + XYKK),
589  (unsigned int)(ii + Xjj + XYKK1), (unsigned int)(ii + 1 + Xjj + XYKK1),
590  (unsigned int)(ii + Xjj1 + XYKK1), (unsigned int)(ii + 1 + Xjj1 + XYKK1)};
591 
592  double val = 0;
593  for (int n = 0; n < 8; n++)
594  val += W[n] * dataPre[index[n]];
595 #ifdef VISP_HAVE_OPENMP
596 #pragma omp critical
597 #endif
598  dataPost[x + m_nbX * y + nbXY * z] = (unsigned char)val;
599  }
600  }
601  }
602  }
603  break;
604  }
605  case GPU_DIRECT_CONVERSION: {
606 #ifdef USTK_HAVE_CUDA
607  this->GPUDirectConversion(dataPost, dataPre);
608 #else
609  throw vpException(
610  vpException::notImplementedError,
611  "usPreScanToPostScan3DConverter::convert: using method GPU_DIRECT_CONVERSION is not implemented yet");
612 #endif
613  break;
614  }
615  case SINGLE_THREAD_FULL_LOOKUP_TABLE: {
616  const unsigned int d = m_SweepInZdirection ? 0 : 1;
617  for (int i = (int)m_lookupTables[d].size() - 1; i >= 0; i--) {
618  double v = 0;
619  for (int j = 0; j < 8; j++)
620  v += m_lookupTables[d][i].m_W[j] * dataPre[m_lookupTables[d][i].m_inputIndex[j]];
621  dataPost[m_lookupTables[d][i].m_outputIndex] = (unsigned char)v;
622  }
623  break;
624  }
625  case MULTI_THREAD_FULL_LOOKUP_TABLE: {
626  const unsigned int d = m_SweepInZdirection ? 0 : 1;
627 #ifdef VISP_HAVE_OPENMP
628 #pragma omp parallel for
629 #endif
630  for (int i = (int)m_lookupTables[d].size() - 1; i >= 0; i--) {
631  double v = 0;
632  for (int j = 0; j < 8; j++)
633  v += m_lookupTables[d][i].m_W[j] * dataPre[m_lookupTables[d][i].m_inputIndex[j]];
634  dataPost[m_lookupTables[d][i].m_outputIndex] = (unsigned char)v;
635  }
636  break;
637  }
638  case GPU_FULL_LOOKUP_TABLE: {
639 #ifdef USTK_HAVE_CUDA
640  this->GPUFullLookupTableConversion(dataPost, dataPre);
641 #else
642  throw vpException(
643  vpException::notImplementedError,
644  "usPreScanToPostScan3DConverter::convert: using method GPU_FULL_LOOKUP_TABLE is not implemented yet");
645 #endif
646  break;
647  }
648  case SINGLE_THREAD_REDUCED_LOOKUP_TABLE: {
649  const unsigned int d = m_SweepInZdirection ? 0 : 1;
650  unsigned int X = m_VpreScan.getWidth();
651  unsigned int Y = m_VpreScan.getHeight();
652  unsigned int XY = X * Y;
653  for (int i = (int)m_reducedLookupTables[d].size() - 1; i >= 0; i--) {
654  const usVoxelWeightAndIndexReducedMemory &m = m_reducedLookupTables[d][i];
655  double u = m.m_W[0];
656  double v = m.m_W[1];
657  double w = m.m_W[2];
658  double u1 = 1 - u;
659  double v1 = 1 - v;
660  double w1 = 1 - w;
661 
662  double v1w1 = v1 * w1;
663  double vw1 = v * w1;
664  double v1w = v1 * w;
665  double vw = v * w;
666 
667  double W[8] = {u1 * v1w1, u * v1w1, u1 * vw1, u * vw1, u1 * v1w, u * v1w, u1 * vw, u * vw};
668 
669  unsigned int index[8] = {m.m_inputIndex, m.m_inputIndex + 1, m.m_inputIndex + X,
670  m.m_inputIndex + 1 + X, m.m_inputIndex + XY, m.m_inputIndex + 1 + XY,
671  m.m_inputIndex + X + XY, m.m_inputIndex + 1 + X + XY};
672 
673  double val = 0;
674  for (int j = 0; j < 8; j++)
675  val += W[j] * dataPre[index[j]];
676  dataPost[m.m_outputIndex] = (unsigned char)val;
677  }
678  break;
679  }
680  case MULTI_THREAD_REDUCED_LOOKUP_TABLE: {
681  const unsigned int d = m_SweepInZdirection ? 0 : 1;
682  unsigned int X = m_VpreScan.getWidth();
683  unsigned int Y = m_VpreScan.getHeight();
684  unsigned int XY = X * Y;
685 #ifdef VISP_HAVE_OPENMP
686 #pragma omp parallel for
687 #endif
688  for (int i = (int)m_reducedLookupTables[d].size() - 1; i >= 0; i--) {
689  const usVoxelWeightAndIndexReducedMemory &m = m_reducedLookupTables[d][i];
690  double u = m.m_W[0];
691  double v = m.m_W[1];
692  double w = m.m_W[2];
693  double u1 = 1 - u;
694  double v1 = 1 - v;
695  double w1 = 1 - w;
696 
697  double v1w1 = v1 * w1;
698  double vw1 = v * w1;
699  double v1w = v1 * w;
700  double vw = v * w;
701 
702  double W[8] = {u1 * v1w1, u * v1w1, u1 * vw1, u * vw1, u1 * v1w, u * v1w, u1 * vw, u * vw};
703 
704  unsigned int index[8] = {m.m_inputIndex, m.m_inputIndex + 1, m.m_inputIndex + X,
705  m.m_inputIndex + 1 + X, m.m_inputIndex + XY, m.m_inputIndex + 1 + XY,
706  m.m_inputIndex + X + XY, m.m_inputIndex + 1 + X + XY};
707 
708  double val = 0;
709  for (int j = 0; j < 8; j++)
710  val += W[j] * dataPre[index[j]];
711  dataPost[m.m_outputIndex] = (unsigned char)val;
712  }
713  break;
714  }
715  case GPU_REDUCED_LOOKUP_TABLE: {
716 #ifdef USTK_HAVE_CUDA
717  this->GPUReducedLookupTableConversion(dataPost, dataPre);
718 #else
719  throw vpException(
720  vpException::notImplementedError,
721  "usPreScanToPostScan3DConverter::convert: using method GPU_REDUCED_LOOKUP_TABLE is not implemented yet");
722 #endif
723  break;
724  }
725  }
726 
727  // writing post-scan image settings
728  postScanImage.setTransducerSettings(m_VpreScan);
729  postScanImage.setMotorSettings(m_VpreScan);
730  postScanImage.setElementSpacingX(m_resolution);
731  postScanImage.setElementSpacingY(m_resolution);
732  postScanImage.setElementSpacingZ(m_resolution);
733  postScanImage.setScanLineDepth(m_resolution * m_VpreScan.getBModeSampleNumber());
734 }
735 
740 void usPreScanToPostScan3DConverter::setConverterOptimizationMethod(usConverterOptimizationMethod method)
741 {
742  m_converterOptimizationMethod = method;
743 
744  switch (m_converterOptimizationMethod) {
745  case SINGLE_THREAD_DIRECT_CONVERSION:
746  case SINGLE_THREAD_FULL_LOOKUP_TABLE:
747  case SINGLE_THREAD_REDUCED_LOOKUP_TABLE: {
748  break;
749  }
750  case MULTI_THREAD_DIRECT_CONVERSION:
751  case MULTI_THREAD_FULL_LOOKUP_TABLE:
752  case MULTI_THREAD_REDUCED_LOOKUP_TABLE: {
753 #ifndef VISP_HAVE_OPENMP
754  std::cout << "Warning in usPreScanToPostScan3DConverter::setConverterOptimizationMethod: OpenMP is not available "
755  "to use multi-thread optimization, will use single thread implementation instead."
756  << std::endl;
757 #endif
758  break;
759  }
760  case GPU_DIRECT_CONVERSION:
761  case GPU_FULL_LOOKUP_TABLE:
762  case GPU_REDUCED_LOOKUP_TABLE: {
763 #ifndef USTK_HAVE_CUDA
764  throw vpException(vpException::notImplementedError, "usPreScanToPostScan3DConverter::"
765  "setConverterOptimizationMethod: cannot use "
766  "GPU conversion without CUDA");
767 #endif
768  break;
769  }
770  }
771 }
772 
784 void usPreScanToPostScan3DConverter::convertPreScanCoordToPostScanCoord(double i_preScan, double j_preScan,
785  double k_preScan, double *i_postScan,
786  double *j_postScan, double *k_postScan,
787  bool sweepInZdirection)
788 {
789  const double Nframe = m_VpreScan.getFrameNumber();
790  const double Nline = m_VpreScan.getScanLineNumber();
791 
792  const double offsetPhi = 0.5 * m_VpreScan.getScanLinePitch() * (Nline - 1);
793  const double offsetTheta = 0.5 * m_VpreScan.getFramePitch() * Nframe;
794 
795  const double r = m_VpreScan.getTransducerRadius() + i_preScan * m_VpreScan.getAxialResolution();
796  const double phi = j_preScan * m_VpreScan.getScanLinePitch() - offsetPhi;
797  // const double theta = (sweepInZdirection ? 1 : -1) * (m_VpreScan.getFramePitch() * Nframe * (j_preScan + Nline *
798  // k_preScan) / (Nframe * Nline - 1) - offsetTheta);
799  const double theta =
800  (m_VpreScan.getFramePitch() * Nframe *
801  ((sweepInZdirection ? j_preScan : Nline - 1 - j_preScan) + Nline * k_preScan) / (Nframe * Nline - 1) -
802  offsetTheta);
803 
804  const double cPhi = cos(phi);
805 
806  if (j_postScan)
807  *j_postScan = r * sin(phi);
808  double radiusOffset = m_VpreScan.getTransducerRadius() - m_VpreScan.getMotorRadius();
809  if (i_postScan)
810  *i_postScan = (r * cPhi - radiusOffset) * cos(theta);
811  if (k_postScan)
812  *k_postScan = (r * cPhi - radiusOffset) * sin(theta);
813 }
814 
826 void usPreScanToPostScan3DConverter::convertPostScanCoordToPreScanCoord(double i_postScan, double j_postScan,
827  double k_postScan, double *i_preScan,
828  double *j_preScan, double *k_preScan,
829  bool sweepInZdirection)
830 {
831  const double Nframe = m_VpreScan.getFrameNumber();
832  const double Nline = m_VpreScan.getScanLineNumber();
833  double radiusOffset = m_VpreScan.getTransducerRadius() - m_VpreScan.getMotorRadius();
834  const double rProbe = radiusOffset + sqrt(i_postScan * i_postScan + k_postScan * k_postScan);
835  const double r = sqrt(rProbe * rProbe + j_postScan * j_postScan);
836  const double phi = atan(j_postScan / rProbe);
837  const double theta = atan(k_postScan / i_postScan);
838 
839  double jtmp = phi / m_VpreScan.getScanLinePitch() + 0.5 * (Nline - 1);
840  if (j_preScan)
841  *j_preScan = jtmp;
842  if (i_preScan)
843  *i_preScan = (r - m_VpreScan.getTransducerRadius()) / m_VpreScan.getAxialResolution();
844  // if (k_preScan) {
845  // *k_preScan =
846  // (Nframe * Nline - 1) *
847  // (0.5 / Nline + (sweepInZdirection ? 1 : -1) * theta / (m_VpreScan.getFramePitch() * Nframe * Nline)) -
848  // jtmp / Nline;
849  //}
850  if (k_preScan) {
851  *k_preScan = (Nframe * Nline - 1) * (0.5 / Nline + theta / (m_VpreScan.getFramePitch() * Nframe * Nline)) -
852  (sweepInZdirection ? jtmp : Nline - 1 - jtmp) / Nline;
853  }
854 }
usImage3D::getNumberOfFrames
unsigned int getNumberOfFrames() const
Definition: usImage3D.h:137
usImage3D::getData
Type * getData()
Definition: usImage3D.h:110
usImage3D::resize
void resize(unsigned int height, unsigned int width, unsigned int numberOfFrames)
Definition: usImage3D.h:376
usImage3D::getConstData
Type * getConstData() const
Definition: usImage3D.h:104
usImage3D::initData
void initData(Type value)
Definition: usImage3D.h:367
usImage3D::getHeight
unsigned int getHeight() const
Definition: usImage3D.h:131
usImage3D::getWidth
unsigned int getWidth() const
Definition: usImage3D.h:125
usImagePostScan3D::setElementSpacingZ
void setElementSpacingZ(double elementSpacingZ)
Definition: usImagePostScan3D.h:322
usImagePostScan3D::setScanLineDepth
void setScanLineDepth(double scanLineDepth)
Definition: usImagePostScan3D.h:331
usImagePostScan3D::setElementSpacingX
void setElementSpacingX(double elementSpacingX)
Definition: usImagePostScan3D.h:304
usImagePostScan3D::setElementSpacingY
void setElementSpacingY(double elementSpacingY)
Definition: usImagePostScan3D.h:313
usImagePreScan3D::getBModeSampleNumber
unsigned int getBModeSampleNumber() const
Definition: usImagePreScan3D.h:250
usImagePreScanSettings
Settings associated to ultrasound pre-scan images implemented in usImageRF2D, usImageRF3D,...
Definition: usImagePreScanSettings.h:81
usMotorSettings
Generic class for 3D ultrasound motor settings associated to the 3D probe used during acquisition.
Definition: usMotorSettings.h:78
usMotorSettings::setMotorSettings
void setMotorSettings(const usMotorSettings &other)
Definition: usMotorSettings.cpp:184
usMotorSettings::getFrameNumber
unsigned int getFrameNumber() const
Definition: usMotorSettings.cpp:159
usMotorSettings::getMotorRadius
double getMotorRadius() const
Definition: usMotorSettings.cpp:141
usMotorSettings::getFramePitch
double getFramePitch() const
Definition: usMotorSettings.cpp:153
usMotorSettings::getMotorType
usMotorType getMotorType() const
Definition: usMotorSettings.cpp:178
usPreScanToPostScan3DConverter::m_lookupTables
std::vector< usVoxelWeightAndIndex > m_lookupTables[2]
Definition: usPreScanToPostScan3DConverter.h:145
usPreScanToPostScan3DConverter::m_VpreScan
usImagePreScan3D< unsigned char > m_VpreScan
Definition: usPreScanToPostScan3DConverter.h:150
usPreScanToPostScan3DConverter::init
void init(const usImagePreScan3D< unsigned char > &preScanImage, double down=1)
Definition: usPreScanToPostScan3DConverter.cpp:71
usPreScanToPostScan3DConverter::m_converterOptimizationMethod
usConverterOptimizationMethod m_converterOptimizationMethod
Definition: usPreScanToPostScan3DConverter.h:141
usPreScanToPostScan3DConverter::convert
void convert(usImagePostScan3D< unsigned char > &postScanImage, const usImagePreScan3D< unsigned char > &preScanImage)
Definition: usPreScanToPostScan3DConverter.cpp:446
usPreScanToPostScan3DConverter::m_reducedLookupTables
std::vector< usVoxelWeightAndIndexReducedMemory > m_reducedLookupTables[2]
Definition: usPreScanToPostScan3DConverter.h:146
usPreScanToPostScan3DConverter::m_conversionOptimizationMethodUsedAtInit
usConverterOptimizationMethod m_conversionOptimizationMethodUsedAtInit
Definition: usPreScanToPostScan3DConverter.h:142
usPreScanToPostScan3DConverter::setConverterOptimizationMethod
usConverterOptimizationMethod setConverterOptimizationMethod() const
Definition: usPreScanToPostScan3DConverter.h:172
usTransducerSettings
Generic class for 2D ultrasound data common settings associated to the type of probe transducer used ...
Definition: usTransducerSettings.h:81
usTransducerSettings::getTransducerRadius
double getTransducerRadius() const
Definition: usTransducerSettings.cpp:163
usTransducerSettings::setTransducerSettings
void setTransducerSettings(const usTransducerSettings &other)
Definition: usTransducerSettings.cpp:210
usTransducerSettings::getScanLineNumber
unsigned int getScanLineNumber() const
Definition: usTransducerSettings.cpp:216