From 1cf797c3b8fcbe925db53c0db7c39bd5ce06a458 Mon Sep 17 00:00:00 2001 From: John Cupitt Date: Thu, 11 May 2017 09:09:49 +0100 Subject: [PATCH] revise dcm rescale and window handling The DICOM reader was multiplying pixel values by `rescale_slope`, which it represented as `ssize_t`. Many DICOM files use double values for `rescale_slope`, so for `rescale_slope` < 1, the reader was multiplying by zero and producing blank images. This patch changes `rescale_slope` and `rescale_intercept` to be double, so dcm read now works with fractional values of rescale. A new option, `dcm:rescale` turns rescale processing on. It is off by default for compatibility with older versions. The `window_center` and `window_width` fields must also be held as doubles. These fields are now only interpreted if rescale processing is enabled, since according to the spec, they should operate upon the rescaled values. A second new option, `dcm:window`, lets users overrride the `window_center` and `window_width` fields from the file with their own settings. This is a standard feature in other DICOM software. Forum thread: https://www.imagemagick.org/discourse-server/viewtopic.php?f=3&t=31887&p=145852 Related PRs: https://github.com/ImageMagick/ImageMagick/pull/477 https://github.com/ImageMagick/ImageMagick/pull/483 https://github.com/ImageMagick/Website/pull/1 DICOM spec: https://www.dabsoft.ch/dicom/3/C.11.2.1.2/ --- ChangeLog | 4 ++ coders/dcm.c | 111 +++++++++++++++++++++++++++++++++------------------ 2 files changed, 77 insertions(+), 38 deletions(-) diff --git a/ChangeLog b/ChangeLog index c59b10bea..843544840 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,7 @@ +2017-05-10 7.0.5-6 John Cupitt + * Revise DICOM window and rescale handling (reference + https://github.com/ImageMagick/ImageMagick/pull/484) + 2017-05-06 7.0.5-6 Cristy * Restore the -alpha Shape option (reference https://www.imagemagick.org/discourse-server/viewtopic.php?f=3&t=31879). diff --git a/coders/dcm.c b/coders/dcm.c index c5446e1a6..02f3854e5 100644 --- a/coders/dcm.c +++ b/coders/dcm.c @@ -2691,13 +2691,16 @@ typedef struct _DCMInfo max_value, samples_per_pixel, signed_data, - significant_bits, - window_width; + significant_bits; - ssize_t + MagickBooleanType + rescale; + + double rescale_intercept, rescale_slope, - window_center; + window_center, + window_width; } DCMInfo; typedef struct _DCMStreamInfo @@ -2871,31 +2874,42 @@ static MagickBooleanType ReadDCMPixels(Image *image,DCMInfo *info, } i++; } - index=(int) ((pixel_value*info->rescale_slope)+ - info->rescale_intercept); - if (info->window_width == 0) + if (info->signed_data == 1) + pixel_value-=32767; + if (info->rescale) { - if (info->signed_data == 1) - index-=32767; + double + scaled_value; + + scaled_value=pixel_value*info->rescale_slope+ + info->rescale_intercept; + if (info->window_width == 0) + { + index=(int) scaled_value; + } + else + { + double + window_max, + window_min; + + window_min=ceil(info->window_center- + (info->window_width-1.0)/2.0-0.5); + window_max=floor(info->window_center+ + (info->window_width-1.0)/2.0+0.5); + if (scaled_value <= window_min) + index=0; + else + if (scaled_value > window_max) + index=(int) info->max_value; + else + index=(int) (info->max_value*(((scaled_value- + info->window_center-0.5)/(info->window_width-1))+0.5)); + } } else { - ssize_t - window_max, - window_min; - - window_min=(ssize_t) ceil((double) info->window_center- - (info->window_width-1.0)/2.0-0.5); - window_max=(ssize_t) floor((double) info->window_center+ - (info->window_width-1.0)/2.0+0.5); - if ((ssize_t)index <= window_min) - index=0; - else - if ((ssize_t)index > window_max) - index=(int) info->max_value; - else - index=(int) (info->max_value*(((index-info->window_center- - 0.5)/(info->window_width-1))+0.5)); + index=pixel_value; } index&=info->mask; index=(int) ConstrainColormapIndex(image,(ssize_t) index,exception); @@ -3074,20 +3088,21 @@ static Image *ReadDCMImage(const ImageInfo *image_info,ExceptionInfo *exception) Read DCM Medical image. */ (void) CopyMagickString(photometric,"MONOCHROME1 ",MagickPathExtent); + info.polarity=MagickFalse; + info.scale=(Quantum *) NULL; info.bits_allocated=8; info.bytes_per_pixel=1; info.depth=8; - info.max_value=255UL; info.mask=0xffff; - info.polarity=MagickFalse; - info.rescale_intercept=0; - info.rescale_slope=1; + info.max_value=255UL; info.samples_per_pixel=1; - info.scale=(Quantum *) NULL; info.signed_data=(~0UL); info.significant_bits=0; - info.window_center=0; - info.window_width=0; + info.rescale=MagickFalse; + info.rescale_intercept=0.0; + info.rescale_slope=1.0; + info.window_center=0.0; + info.window_width=0.0; data=(unsigned char *) NULL; element=0; explicit_vr[2]='\0'; @@ -3493,7 +3508,7 @@ static Image *ReadDCMImage(const ImageInfo *image_info,ExceptionInfo *exception) Visible pixel range: center. */ if (data != (unsigned char *) NULL) - info.window_center=(ssize_t) StringToLong((char *) data); + info.window_center=StringToDouble((char *) data, (char **) NULL); break; } case 0x1051: @@ -3502,7 +3517,7 @@ static Image *ReadDCMImage(const ImageInfo *image_info,ExceptionInfo *exception) Visible pixel range: width. */ if (data != (unsigned char *) NULL) - info.window_width=StringToUnsignedLong((char *) data); + info.window_width=StringToDouble((char *) data, (char **) NULL); break; } case 0x1052: @@ -3511,7 +3526,8 @@ static Image *ReadDCMImage(const ImageInfo *image_info,ExceptionInfo *exception) Rescale intercept */ if (data != (unsigned char *) NULL) - info.rescale_intercept=(ssize_t) StringToLong((char *) data); + info.rescale_intercept=StringToDouble((char *) data, + (char **) NULL); break; } case 0x1053: @@ -3520,7 +3536,7 @@ static Image *ReadDCMImage(const ImageInfo *image_info,ExceptionInfo *exception) Rescale slope */ if (data != (unsigned char *) NULL) - info.rescale_slope=(ssize_t) StringToLong((char *) data); + info.rescale_slope=StringToDouble((char *) data, (char **) NULL); break; } case 0x1200: @@ -4063,14 +4079,33 @@ static Image *ReadDCMImage(const ImageInfo *image_info,ExceptionInfo *exception) /* Convert DCM Medical image to pixel packets. */ - if ((info.window_center != 0) && (info.window_width == 0)) - info.window_width=(size_t) info.window_center; option=GetImageOption(image_info,"dcm:display-range"); if (option != (const char *) NULL) { if (LocaleCompare(option,"reset") == 0) info.window_width=0; } + option=GetImageOption(image_info,"dcm:window"); + if (option != (char *) NULL) + { + GeometryInfo + geometry_info; + + MagickStatusType + flags; + + flags=ParseGeometry(option,&geometry_info); + if (flags & RhoValue) + info.window_center=geometry_info.rho; + if (flags & SigmaValue) + info.window_width=geometry_info.sigma; + info.rescale=MagickTrue; + } + option=GetImageOption(image_info,"dcm:rescale"); + if (option != (char *) NULL) + info.rescale=IsStringTrue(option); + if ((info.window_center != 0) && (info.window_width == 0)) + info.window_width=info.window_center; status=ReadDCMPixels(image,&info,stream_info,MagickTrue,exception); if ((status != MagickFalse) && (stream_info->segment_count > 1)) { -- 2.40.0