Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
vtkMitkLevelWindowFilter.cpp
Go to the documentation of this file.
1 /*===================================================================
2 
3 The Medical Imaging Interaction Toolkit (MITK)
4 
5 Copyright (c) German Cancer Research Center,
6 Division of Medical and Biological Informatics.
7 All rights reserved.
8 
9 This software is distributed WITHOUT ANY WARRANTY; without
10 even the implied warranty of MERCHANTABILITY or FITNESS FOR
11 A PARTICULAR PURPOSE.
12 
13 See LICENSE.txt or http://www.mitk.org for details.
14 
15 ===================================================================*/
16 
18 #include "vtkObjectFactory.h"
19 #include <vtkColorTransferFunction.h>
20 #include <vtkImageData.h>
21 #include <vtkImageIterator.h>
22 #include <vtkInformation.h>
23 #include <vtkInformationVector.h>
24 #include <vtkLookupTable.h>
25 #include <vtkPiecewiseFunction.h>
26 
27 #include <vtkStreamingDemandDrivenPipeline.h>
28 
29 // used for acos etc.
30 #include <cmath>
31 
32 // used for PI
33 #include <itkMath.h>
34 
35 #include <mitkLogMacros.h>
36 
37 static const double PI = itk::Math::pi;
38 
40 
42  : m_LookupTable(nullptr), m_OpacityFunction(nullptr), m_MinOpacity(0.0), m_MaxOpacity(255.0)
43 {
44  // MITK_INFO << "mitk level/window filter uses " << GetNumberOfThreads() << " thread(s)";
45 }
46 
48 {
49 }
50 
52 {
53  unsigned long mTime = this->vtkObject::GetMTime();
54  unsigned long time;
55 
56  if (this->m_LookupTable != nullptr)
57  {
58  time = this->m_LookupTable->GetMTime();
59  mTime = (time > mTime ? time : mTime);
60  }
61 
62  return mTime;
63 }
64 
65 void vtkMitkLevelWindowFilter::SetLookupTable(vtkScalarsToColors *lookupTable)
66 {
67  if (m_LookupTable != lookupTable)
68  {
69  m_LookupTable = lookupTable;
70  this->Modified();
71  }
72 }
73 
75 {
76  return m_LookupTable;
77 }
78 
79 void vtkMitkLevelWindowFilter::SetOpacityPiecewiseFunction(vtkPiecewiseFunction *opacityFunction)
80 {
81  if (m_OpacityFunction != opacityFunction)
82  {
83  m_OpacityFunction = opacityFunction;
84  this->Modified();
85  }
86 }
87 
88 // This code was copied from the iil. The template works only for float and double.
89 // Internal method which should never be used anywhere else and should not be in th header.
90 // Convert color pixels from (R,G,B) to (H,S,I).
91 // Reference: "Digital Image Processing, 2nd. edition", R. Gonzalez and R. Woods. Prentice Hall, 2002.
92 template <class T>
93 void RGBtoHSI(T *RGB, T *HSI)
94 {
95  T R = RGB[0], G = RGB[1], B = RGB[2], nR = (R < 0 ? 0 : (R > 255 ? 255 : R)) / 255,
96  nG = (G < 0 ? 0 : (G > 255 ? 255 : G)) / 255, nB = (B < 0 ? 0 : (B > 255 ? 255 : B)) / 255,
97  m = nR < nG ? (nR < nB ? nR : nB) : (nG < nB ? nG : nB),
98  theta = (T)(std::acos(0.5f * ((nR - nG) + (nR - nB)) / std::sqrt(std::pow(nR - nG, 2) + (nR - nB) * (nG - nB))) *
99  180 / PI),
100  sum = nR + nG + nB;
101  T H = 0, S = 0, I = 0;
102  if (theta > 0)
103  H = (nB <= nG) ? theta : 360 - theta;
104  if (sum > 0)
105  S = 1 - 3 / sum * m;
106  I = sum / 3;
107  HSI[0] = (T)H;
108  HSI[1] = (T)S;
109  HSI[2] = (T)I;
110 }
111 
112 // This code was copied from the iil. The template works only for float and double.
113 // Internal method which should never be used anywhere else and should not be in th header.
114 // Convert color pixels from (H,S,I) to (R,G,B).
115 template <class T>
116 void HSItoRGB(T *HSI, T *RGB)
117 {
118  T H = (T)HSI[0], S = (T)HSI[1], I = (T)HSI[2], a = I * (1 - S), R = 0, G = 0, B = 0;
119  if (H < 120)
120  {
121  B = a;
122  R = (T)(I * (1 + S * std::cos(H * PI / 180) / std::cos((60 - H) * PI / 180)));
123  G = 3 * I - (R + B);
124  }
125  else if (H < 240)
126  {
127  H -= 120;
128  R = a;
129  G = (T)(I * (1 + S * std::cos(H * PI / 180) / std::cos((60 - H) * PI / 180)));
130  B = 3 * I - (R + G);
131  }
132  else
133  {
134  H -= 240;
135  G = a;
136  B = (T)(I * (1 + S * std::cos(H * PI / 180) / std::cos((60 - H) * PI / 180)));
137  R = 3 * I - (G + B);
138  }
139  R *= 255;
140  G *= 255;
141  B *= 255;
142  RGB[0] = (T)(R < 0 ? 0 : (R > 255 ? 255 : R));
143  RGB[1] = (T)(G < 0 ? 0 : (G > 255 ? 255 : G));
144  RGB[2] = (T)(B < 0 ? 0 : (B > 255 ? 255 : B));
145 }
146 
147 // Internal method which should never be used anywhere else and should not be in th header.
148 //----------------------------------------------------------------------------
149 // This templated function executes the filter for any type of data.
150 template <class T>
152  vtkImageData *inData,
153  vtkImageData *outData,
154  int outExt[6],
155  double *clippingBounds,
156  T *)
157 {
158  vtkImageIterator<T> inputIt(inData, outExt);
159  vtkImageIterator<T> outputIt(outData, outExt);
160  vtkLookupTable *lookupTable;
161  const int maxC = inData->GetNumberOfScalarComponents();
162 
163  double tableRange[2];
164 
165  lookupTable = dynamic_cast<vtkLookupTable *>(self->GetLookupTable());
166 
167  lookupTable->GetTableRange(tableRange);
168 
169  // parameters for RGB level window
170  const double scale = (tableRange[1] - tableRange[0] > 0 ? 255.0 / (tableRange[1] - tableRange[0]) : 0.0);
171  const double bias = tableRange[0] * scale;
172 
173  // parameters for opaque level window
174  const double scaleOpac =
175  (self->GetMaxOpacity() - self->GetMinOpacity() > 0 ? 255.0 / (self->GetMaxOpacity() - self->GetMinOpacity()) : 0.0);
176  const double biasOpac = self->GetMinOpacity() * scaleOpac;
177 
178  int y = outExt[2];
179 
180  // Loop through ouput pixels
181  while (!outputIt.IsAtEnd())
182  {
183  T *inputSI = inputIt.BeginSpan();
184  T *outputSI = outputIt.BeginSpan();
185  T *outputSIEnd = outputIt.EndSpan();
186 
187  if (y >= clippingBounds[2] && y < clippingBounds[3])
188  {
189  int x = outExt[0];
190 
191  while (outputSI != outputSIEnd)
192  {
193  if (x >= clippingBounds[0] && x < clippingBounds[1])
194  {
195  double rgb[3], alpha, hsi[3];
196 
197  // level/window mechanism for intensity in HSI space
198  rgb[0] = static_cast<double>(*inputSI);
199  inputSI++;
200  rgb[1] = static_cast<double>(*inputSI);
201  inputSI++;
202  rgb[2] = static_cast<double>(*inputSI);
203  inputSI++;
204 
205  RGBtoHSI<double>(rgb, hsi);
206  hsi[2] = hsi[2] * 255.0 * scale - bias;
207  hsi[2] = (hsi[2] > 255.0 ? 255 : (hsi[2] < 0.0 ? 0 : hsi[2]));
208  hsi[2] /= 255.0;
209  HSItoRGB<double>(hsi, rgb);
210 
211  *outputSI = static_cast<T>(rgb[0]);
212  outputSI++;
213  *outputSI = static_cast<T>(rgb[1]);
214  outputSI++;
215  *outputSI = static_cast<T>(rgb[2]);
216  outputSI++;
217 
218  unsigned char finalAlpha = 255;
219 
220  // RGBA case
221  if (maxC >= 4)
222  {
223  // level/window mechanism for opacity
224  alpha = static_cast<double>(*inputSI);
225  inputSI++;
226  alpha = alpha * scaleOpac - biasOpac;
227  if (alpha > 255.0)
228  {
229  alpha = 255.0;
230  }
231  else if (alpha < 0.0)
232  {
233  alpha = 0.0;
234  }
235  finalAlpha = static_cast<unsigned char>(alpha);
236 
237  for (int c = 4; c < maxC; c++)
238  inputSI++;
239  }
240 
241  *outputSI = static_cast<T>(finalAlpha);
242  outputSI++;
243  }
244  else
245  {
246  inputSI += maxC;
247  *outputSI = 0;
248  outputSI++;
249  *outputSI = 0;
250  outputSI++;
251  *outputSI = 0;
252  outputSI++;
253  *outputSI = 0;
254  outputSI++;
255  }
256 
257  x++;
258  }
259  }
260  else
261  {
262  while (outputSI != outputSIEnd)
263  {
264  *outputSI = 0;
265  outputSI++;
266  *outputSI = 0;
267  outputSI++;
268  *outputSI = 0;
269  outputSI++;
270  *outputSI = 0;
271  outputSI++;
272  }
273  }
274  inputIt.NextSpan();
275  outputIt.NextSpan();
276  y++;
277  }
278 }
279 
280 // Internal method which should never be used anywhere else and should not be in th header.
281 //----------------------------------------------------------------------------
282 // This templated function executes the filter for any type of data.
283 template <class T>
285  vtkMitkLevelWindowFilter *self, vtkImageData *inData, vtkImageData *outData, int outExt[6], T *)
286 {
287  vtkImageIterator<T> inputIt(inData, outExt);
288  vtkImageIterator<unsigned char> outputIt(outData, outExt);
289 
290  double tableRange[2];
291 
292  // access vtkLookupTable
293  vtkLookupTable *lookupTable = dynamic_cast<vtkLookupTable *>(self->GetLookupTable());
294  lookupTable->GetTableRange(tableRange);
295 
296  // access elements of the vtkLookupTable
297  const int *realLookupTable = reinterpret_cast<int *>(lookupTable->GetTable()->GetPointer(0));
298  int maxIndex = lookupTable->GetNumberOfColors() - 1;
299 
300  const float scale = (tableRange[1] - tableRange[0] > 0 ? (maxIndex + 1) / (tableRange[1] - tableRange[0]) : 0.0);
301  // ensuring that starting point is zero
302  float bias = -tableRange[0] * scale;
303  // due to later conversion to int for rounding
304  bias += 0.5f;
305 
306  // Loop through ouput pixels
307  while (!outputIt.IsAtEnd())
308  {
309  unsigned char *outputSI = outputIt.BeginSpan();
310  unsigned char *outputSIEnd = outputIt.EndSpan();
311 
312  T *inputSI = inputIt.BeginSpan();
313 
314  while (outputSI != outputSIEnd)
315  {
316  // map to an index
317  int idx = static_cast<int>(*inputSI * scale + bias);
318 
319  if (idx < 0)
320  idx = 0;
321  else if (idx > maxIndex)
322  idx = maxIndex;
323 
324  *reinterpret_cast<int *>(outputSI) = realLookupTable[idx];
325 
326  inputSI++;
327  outputSI += 4;
328  }
329 
330  inputIt.NextSpan();
331  outputIt.NextSpan();
332  }
333 }
334 
335 // Internal method which should never be used anywhere else and should not be in th header.
336 //----------------------------------------------------------------------------
337 // This templated function executes the filter for any type of data.
338 template <class T>
340  vtkImageData *inData,
341  vtkImageData *outData,
342  int outExt[6],
343  double *clippingBounds,
344  T *)
345 {
346  vtkImageIterator<T> inputIt(inData, outExt);
347  vtkImageIterator<unsigned char> outputIt(outData, outExt);
348  vtkScalarsToColors *lookupTable = self->GetLookupTable();
349 
350  int y = outExt[2];
351 
352  // Loop through ouput pixels
353  while (!outputIt.IsAtEnd())
354  {
355  unsigned char *outputSI = outputIt.BeginSpan();
356  unsigned char *outputSIEnd = outputIt.EndSpan();
357 
358  // do we iterate over the inner vertical clipping bounds
359  if (y >= clippingBounds[2] && y < clippingBounds[3])
360  {
361  T *inputSI = inputIt.BeginSpan();
362 
363  int x = outExt[0];
364 
365  while (outputSI != outputSIEnd)
366  {
367  // is this pixel within horizontal clipping bounds
368  if (x >= clippingBounds[0] && x < clippingBounds[1])
369  {
370  // fetching original value
371  double grayValue = static_cast<double>(*inputSI);
372  // applying lookuptable - copy the 4 (RGBA) chars as a single int
373  *reinterpret_cast<int *>(outputSI) = *reinterpret_cast<int *>(lookupTable->MapValue(grayValue));
374  }
375  else
376  {
377  // outer horizontal clipping bounds - write a transparent RGBA pixel as a single int
378  *reinterpret_cast<int *>(outputSI) = 0;
379  }
380 
381  inputSI++;
382  outputSI += 4;
383  x++;
384  }
385  }
386  else
387  {
388  // outer vertical clipping bounds - write a transparent RGBA line as ints
389  while (outputSI != outputSIEnd)
390  {
391  *reinterpret_cast<int *>(outputSI) = 0;
392  outputSI += 4;
393  }
394  }
395 
396  inputIt.NextSpan();
397  outputIt.NextSpan();
398  y++;
399  }
400 }
401 
402 // Internal method which should never be used anywhere else and should not be in th header.
403 //----------------------------------------------------------------------------
404 // This templated function executes the filter for any type of data.
405 template <class T>
407  vtkImageData *inData,
408  vtkImageData *outData,
409  int outExt[6],
410  double *clippingBounds,
411  T *)
412 {
413  vtkImageIterator<T> inputIt(inData, outExt);
414  vtkImageIterator<unsigned char> outputIt(outData, outExt);
415  vtkColorTransferFunction *lookupTable = dynamic_cast<vtkColorTransferFunction *>(self->GetLookupTable());
416  vtkPiecewiseFunction *opacityFunction = self->GetOpacityPiecewiseFunction();
417 
418  int y = outExt[2];
419 
420  // Loop through ouput pixels
421  while (!outputIt.IsAtEnd())
422  {
423  unsigned char *outputSI = outputIt.BeginSpan();
424  unsigned char *outputSIEnd = outputIt.EndSpan();
425 
426  // do we iterate over the inner vertical clipping bounds
427  if (y >= clippingBounds[2] && y < clippingBounds[3])
428  {
429  T *inputSI = inputIt.BeginSpan();
430 
431  int x = outExt[0];
432 
433  while (outputSI != outputSIEnd)
434  {
435  // is this pixel within horizontal clipping bounds
436  if (x >= clippingBounds[0] && x < clippingBounds[1])
437  {
438  // fetching original value
439  double grayValue = static_cast<double>(*inputSI);
440 
441  // applying directly colortransferfunction
442  // because vtkColorTransferFunction::MapValue is not threadsafe
443  double rgba[4];
444  lookupTable->GetColor(grayValue, rgba); // RGB mapping
445  rgba[3] = 1.0;
446  if (opacityFunction)
447  rgba[3] = opacityFunction->GetValue(grayValue); // Alpha mapping
448 
449  for (int i = 0; i < 4; ++i)
450  {
451  outputSI[i] = static_cast<unsigned char>(255.0 * rgba[i] + 0.5);
452  }
453  }
454  else
455  {
456  // outer horizontal clipping bounds - write a transparent RGBA pixel as a single int
457  *reinterpret_cast<int *>(outputSI) = 0;
458  }
459 
460  inputSI++;
461  outputSI += 4;
462  x++;
463  }
464  }
465  else
466  {
467  // outer vertical clipping bounds - write a transparent RGBA line as ints
468  while (outputSI != outputSIEnd)
469  {
470  *reinterpret_cast<int *>(outputSI) = 0;
471  outputSI += 4;
472  }
473  }
474 
475  inputIt.NextSpan();
476  outputIt.NextSpan();
477  y++;
478  }
479 }
480 
482  vtkInformationVector **inputVector,
483  vtkInformationVector *outputVector)
484 {
485  vtkInformation *outInfo = outputVector->GetInformationObject(0);
486 
487  // do nothing except copy scalar type info
488  this->CopyInputArrayAttributesToOutput(request, inputVector, outputVector);
489 
490  vtkDataObject::SetPointDataActiveScalarInfo(outInfo, VTK_UNSIGNED_CHAR, 4);
491 
492  return 1;
493 }
494 
495 // Method to run the filter in different threads.
496 void vtkMitkLevelWindowFilter::ThreadedExecute(vtkImageData *inData, vtkImageData *outData, int extent[6], int /*id*/)
497 {
498  if (inData->GetNumberOfScalarComponents() > 2)
499  {
500  switch (inData->GetScalarType())
501  {
502  vtkTemplateMacro(
503  vtkApplyLookupTableOnRGBA(this, inData, outData, extent, m_ClippingBounds, static_cast<VTK_TT *>(nullptr)));
504  default:
505  vtkErrorMacro(<< "Execute: Unknown ScalarType");
506  return;
507  }
508  }
509  else
510  {
511  bool dontClip = extent[2] >= m_ClippingBounds[2] && extent[3] <= m_ClippingBounds[3] &&
512  extent[0] >= m_ClippingBounds[0] && extent[1] <= m_ClippingBounds[1];
513 
514  if (this->GetLookupTable())
515  this->GetLookupTable()->Build();
516 
517  vtkLookupTable *vlt = dynamic_cast<vtkLookupTable *>(this->GetLookupTable());
518  vtkColorTransferFunction *ctf = dynamic_cast<vtkColorTransferFunction *>(this->GetLookupTable());
519 
520  bool linearLookupTable = vlt && vlt->GetScale() == VTK_SCALE_LINEAR;
521 
522  bool useFast = dontClip && linearLookupTable;
523 
524  if (ctf)
525  {
526  switch (inData->GetScalarType())
527  {
528  vtkTemplateMacro(vtkApplyLookupTableOnScalarsCTF(
529  this, inData, outData, extent, m_ClippingBounds, static_cast<VTK_TT *>(nullptr)));
530  default:
531  vtkErrorMacro(<< "Execute: Unknown ScalarType");
532  return;
533  }
534  }
535  else if (useFast)
536  {
537  switch (inData->GetScalarType())
538  {
539  vtkTemplateMacro(
540  vtkApplyLookupTableOnScalarsFast(this, inData, outData, extent, static_cast<VTK_TT *>(nullptr)));
541  default:
542  vtkErrorMacro(<< "Execute: Unknown ScalarType");
543  return;
544  }
545  }
546  else
547  {
548  switch (inData->GetScalarType())
549  {
550  vtkTemplateMacro(vtkApplyLookupTableOnScalars(
551  this, inData, outData, extent, m_ClippingBounds, static_cast<VTK_TT *>(nullptr)));
552  default:
553  vtkErrorMacro(<< "Execute: Unknown ScalarType");
554  return;
555  }
556  }
557  }
558 }
559 
560 // void vtkMitkLevelWindowFilter::ExecuteInformation(
561 // vtkImageData *vtkNotUsed(inData), vtkImageData *vtkNotUsed(outData))
562 //{
563 //}
564 
566 {
567  m_MinOpacity = minOpacity;
568 }
569 
571 {
572  return m_MinOpacity;
573 }
574 
576 {
577  m_MaxOpacity = maxOpacity;
578 }
579 
581 {
582  return m_MaxOpacity;
583 }
584 
585 void vtkMitkLevelWindowFilter::SetClippingBounds(double *bounds) // TODO does double[4] work??
586 {
587  for (unsigned int i = 0; i < 4; ++i)
588  m_ClippingBounds[i] = bounds[i];
589 }
void SetOpacityPiecewiseFunction(vtkPiecewiseFunction *opacityFunction)
Set the piecewise function used to map scalar to alpha component value (only used when the lookupTabl...
void vtkApplyLookupTableOnScalarsFast(vtkMitkLevelWindowFilter *self, vtkImageData *inData, vtkImageData *outData, int outExt[6], T *)
static const double PI
void SetLookupTable(vtkScalarsToColors *lookupTable)
Set the lookup table for the RGB level window.
vtkScalarsToColors * GetLookupTable()
Get the lookup table for the RGB level window.
void HSItoRGB(T *HSI, T *RGB)
void vtkApplyLookupTableOnScalars(vtkMitkLevelWindowFilter *self, vtkImageData *inData, vtkImageData *outData, int outExt[6], double *clippingBounds, T *)
vtkStandardNewMacro(vtkMitkLevelWindowFilter)
void SetClippingBounds(double *)
Set clipping bounds for the opaque part of the resliced 2d image.
void RGBtoHSI(T *RGB, T *HSI)
virtual unsigned long int GetMTime() override
void ThreadedExecute(vtkImageData *inData, vtkImageData *outData, int extent[6], int id) override
Method for threaded execution of the filter.
void vtkApplyLookupTableOnRGBA(vtkMitkLevelWindowFilter *self, vtkImageData *inData, vtkImageData *outData, int outExt[6], double *clippingBounds, T *)
void SetMaxOpacity(double maxOpacity)
Get/Set the upper window opacity for the alpha level window.
Applies the grayvalue or color/opacity level window to scalar or RGB(A) images.
int RequestInformation(vtkInformation *request, vtkInformationVector **inputVector, vtkInformationVector *outputVector) override
void vtkApplyLookupTableOnScalarsCTF(vtkMitkLevelWindowFilter *self, vtkImageData *inData, vtkImageData *outData, int outExt[6], double *clippingBounds, T *)
void SetMinOpacity(double minOpacity)
Get/Set the lower window opacity for the alpha level window.