dicom.pas

用Delphi做的读取医学图片的控件及代码

        inc(lPos);
  end;
    try
       result := strtofloat(lStr);
    except
          on EConvertError do begin
             showmessage('Unable to convert the string '+lStr+' to a number');
             result := 1;
             exit;
          end;
    end; {except}
  end;
function readInterStr:string;
var lStr: string;
begin
  lStr := '';
  While (lPos <= lLen) and (lInStr[lPos] = ' ') do begin
        inc(lPos);
  end;
  While (lPos <= lLen) and (lInStr[lPos] <> ';') do begin
        lStr := lStr+upcase(linStr[lPos]); //zebra upcase
        inc(lPos);
  end;
  result := lStr;
end; //interstr func
begin
  lHdrOK := false;
  lFloat := false;
  lUnsigned := false;
  lImageFormatOK := true;
  Clear_Dicom_Data(lDicomData);
  lDynStr := '';
  FileMode := 0; //set to readonly
  AssignFile(fp, lFileName);
  Reset(fp, 1);
  FileSz := FileSize(fp);
  lHdrEnd := false;
  //lDicomData.ImageStart := FileSz;
  GetMem( lCharRA, FileSz+1 );
  BlockRead(fp, lCharRA^, FileSz, linpos);
  if lInPos <> FileSz then showmessage('Disk error: Unable to read full input file.');
  linPos := 1;
  CloseFile(fp);
  FileMode := 2; //set to read/write
repeat
  linstr := '';
  while (linPos < FileSz) and (lCharRA[linPos] <> ord(kCR)) and (lCharRA[linPos] <> ord(UNIXeoln)) do begin
      lInStr := lInstr + chr(lCharRA[linPos]);
      inc(linPos);
  end;
  inc(lInPos);  //read EOLN
  lLen := length(lInStr);
  lPos := 1;
  lUpcaseStr := '';
  While (lPos <= lLen) and (lInStr[lPos] <> ';') and (lInStr[lPos] <> '=') and (lUpCaseStr <>'INTERFILE') do begin
        if lInStr[lPos] in ['[',']','(',')','/','+','-',{' ',} '0'..'9','a'..'z','A'..'Z'] then
           lUpCaseStr := lUpCaseStr+upcase(linStr[lPos]);
        inc(lPos);
  end;
  inc(lPos); {read equal sign in := statement}
  if lUpCaseStr ='INTERFILE' then begin
     lHdrOK := true;
     lDicomData.little_endian := 0;
     end;
  if lUpCaseStr ='DATASTARTINGBLOCK'then lDicomData.ImageStart := 2048 * round(readInterFloat);
  if lUpCaseStr ='DATAOFFSETINBYTES'then lDicomData.ImageStart := round(readInterFloat);
  if (lUpCaseStr ='MATRIXSIZE[1]') or (lUpCaseStr ='MATRIXSIZE[X]') then lDicomData.XYZdim[1] :=  round(readInterFloat);
  if (lUpCaseStr ='MATRIXSIZE[2]')or (lUpCaseStr ='MATRIXSIZE[Y]')then lDicomData.XYZdim[2] :=  round(readInterFloat);
  if (lUpCaseStr ='MATRIXSIZE[3]')or (lUpCaseStr ='MATRIXSIZE[Z]') or (lUpCaseStr ='NUMBEROFSLICES') or (lUpCaseStr ='TOTALNUMBEROFIMAGES') then begin
     lDicomData.XYZdim[3] :=  round(readInterFloat);
  end;
  if lUpCaseStr ='IMAGEDATABYTEORDER' then begin
     if readInterStr = 'LITTLEENDIAN' then lDicomData.little_endian := 1;
  end;
  if lUpCaseStr ='NUMBERFORMAT' then begin
      lTmpStr := readInterStr;
      if (lTmpStr = 'ASCII') or (lTmpStr='BIT') then begin
         lHdrOK := false;
         showmessage('This software can not convert '+lTmpStr+' data type.');
         goto 333;
      end;
      if lTmpStr = 'UNSIGNEDINTEGER' then lUnsigned := true;
      if (lTmpStr='FLOAT') or (lTmpStr='SHORTFLOAT') or (lTmpStr='LONGFLOAT') then begin
         lFloat := true;
      end;
  end;
  if lUpCaseStr ='NAMEOFDATAFILE' then lFileName := ExtractFilePath(lFileName)+readInterStr;
  if lUpCaseStr ='NUMBEROFBYTESPERPIXEL' then
     lDicomData.Allocbits_per_pixel :=  round(readInterFloat)*8;
  if (lUpCaseStr ='SCALINGFACTOR(MM/PIXEL)[1]') or (lUpCaseStr ='SCALINGFACTOR(MM/PIXEL)[X]') then
     lDicomData.XYZmm[1] :=  (readInterFloat);
  if (lUpCaseStr ='SCALINGFACTOR(MM/PIXEL)[2]') or (lUpCaseStr ='SCALINGFACTOR(MM/PIXEL)[Y]')then lDicomData.XYZmm[2] :=  (readInterFloat);
  if (lUpCaseStr ='SCALINGFACTOR(MM/PIXEL)[3]')or (lUpCaseStr ='SCALINGFACTOR(MM/PIXEL)[Z]')or (lUpCaseStr ='SLICETHICKNESS')then lDicomData.XYZmm[3] :=  (readInterFloat);
  if (lUpCaseStr ='ENDOFINTERFILE') then lHdrEnd := true;
  if not lHdrOK then goto 333;
  if lInStr <> '' then
     lDynStr := lDynStr + lInStr+kCr;
  lHdrOK := true;
until (linPos >= FileSz) or (lHdrEnd){EOF(fp)};
lDicomData.Storedbits_per_pixel := lDicomData.Allocbits_per_pixel;
lImageFormatOK := true;
if (not lFLoat) and (lUnsigned) and ((lDicomData.Storedbits_per_pixel = 16)) then begin
   showmessage('Warning: this Interfile image uses UNSIGNED 16-bit data [values 0..65535]. Analyze specifies SIGNED 16-bit data [-32768..32767]. Some images may not transfer well. [Future versions of MRIcro should fix this].');
   lImageFormatOK := false;
end else if (not lFLoat) and (lDicomData.Storedbits_per_pixel > 16) then begin
   showmessage('WARNING: The image '+lFileName+' is a '+inttostr(lDicomData.Storedbits_per_pixel)+'-bit integer data type. This software may display this as SIGNED data. Bits per voxel: '+inttostr(lDicomData.Storedbits_per_pixel));
   lImageFormatOK := false;
end else if (lFloat) then begin //zebra change float check
   //showmessage('WARNING: The image '+lFileName+' uses floating point [real] numbers. The current software can only read integer data type Interfile images.');
   lDicomData.Float := true;
   //lImageFormatOK := false;
end;
333:
FreeMem( lCharRA);
end; //interfile




//vista
(*    { "bit",		VBitRepn },
    { "double",		VDoubleRepn },
    { "float",		VFloatRepn },
    { "long",		VLongRepn },
    { "sbyte",		VSByteRepn },
    { "short",		VShortRepn },
    { "ubyte",		VUByteRepn },
    { NULL }*)
procedure read_vista_data(lConvertToAnalyze,lAnonymize: boolean; var lDICOMdata: DICOMdata; var lHdrOK, lImageFormatOK: boolean; var lDynStr: string;var lFileName: string);
label 333;
const kUNIXeoln = ord(10);
      k0A = ord($0A);
      kDelimiter = ord(':');
      kPCeoln = ord(kCR);
      kMaxImages = 255;
var
  lAnalyzeHdr: AHdr;
  lDicomRA: array [1..kMaxImages] of DICOMdata;
  lSliceTimeReportedRA: array [1..kMaxImages] of boolean;
  lEdges,lSliceTimeReported: boolean;
  lXori,lYori,lZori: double;
  lTmpStr,lUpCaseStr,lUpCaseValue: string;
  lOutputFiles,lZSLice,l2DSz,lLowLine,lHighLine,lLine,lHalfLines,lLineLen,lVolumeSz,lOutputVolumes,lOutputSlices,lVolume,lTimeSlice,lTimeSLices,lZSlices,lVolumes,lImageSz,lnImages,lLen,FileSz,linPos: integer;
  lLineP,lP: bytep;
  fp,lOutFile: file;
  lCharRA: bytep;

function readVistaFloat:single;
begin
    try
       result := strtofloat(lUpCaseValue);
    except
          on EConvertError do begin
             showmessage('Unable to convert the string '+lUpCaseValue+' to a number');
             result := 1;
             exit;
          end;
    end; {except}
end;


procedure readVistaFloats (var lS1,lS2,lS3:double);
var lInPos,lLen,lVal: integer;
  lStr: string;
begin
   lVal := 0;
   lLen := length(lUpCaseValue);
   lInPos := 1;
   lS1 := 1;
   lS2 := 1;
   lS3 := 3;
   repeat
    inc(lVal);
    lStr := '';
    While (lInPos <= lLen) and ((lStr='') or (lUpCaseValue[lINPos] <> ' '))   do begin
        if lUpCaseValue[lINPos] in ['+','-','.','0'..'9','e','E'] then
           lStr := lStr+lUpCaseValue[lINPos]; //zebra upcase
        inc(lINPos);
    end;
    //showmessage(inttostr(linPos)+':'+inttostr(lLen)+'@'+lStr);
    if length(lStr) > 0 then begin
    try
       case lVal of
            1: lS1 := strtofloat(lStr);
            2: lS2 := strtofloat(lStr);
            3: lS3 := strtofloat(lStr);
       end;
    except
          on EConvertError do begin
             showmessage('Unable to convert the string '+lStr+' to a number');
             exit;
          end;
    end; {except}
    end; //length of string > 0
   until (lInPos >= lLen) or (lVal >= 3);
end;

procedure readVistaStr;
var lStr: string;
begin
  lStr := '';
  While (lInPos <= lLen) and (lCharRA[lINPos] <> kDelimiter) and (lCharRA[lINPos] <> kUnixEOLN) and (lCharRA[lINPos] <> kPCeoln) and (lCharRA[lINPos] <> k0a)  do begin
        if chr(lCharRA[lINPos]) in ['{','}','_','+','-','.','0'..'9','a'..'z','A'..'Z'] then
           lStr := lStr+upcase(chr(lCharRA[lINPos])); //zebra upcase
        inc(lINPos);
  end;
  lUpCaseStr := lStr;
  lStr := '';
  if (lCharRA[lINPos] = kDelimiter) then begin
     inc(lINPos);
     if chr(lCharRA[lINPos]) = ' ' then //strip leading space - 'UBYTE' same as ' UBYTE'
        inc(lInPos);
     While (lInPos <= lLen) and (lCharRA[lINPos] <> kUnixEOLN) and (lCharRA[lINPos] <> kPCeoln) and (lCharRA[lINPos] <> k0a)  do begin
        if chr(lCharRA[lINPos]) in ['{','}',' ',':','+','-','.','0'..'9','a'..'z','A'..'Z'] then
           lStr := lStr+upcase(chr(lCharRA[lINPos]));
        inc(lINPos);
     end;
  end;
  lUpCaseValue := lStr;
  inc(lINPos);
end; //interstr func
begin
  lHdrOK := false;
  lImageFormatOK := false;
  //lFloat := false;
  //lUnsigned := false;
  lXori := 0;
  lYori := 0;
  lZori := 0;
  lDynStr := '';
  if not fileexists(lFilename) then exit;
  FileMode := 0; //set to readonly
  AssignFile(fp, lFileName);
  Reset(fp, 1);
  FileSz := FileSize(fp);
  if FileSz < 160 then begin
      CloseFile(fp);
      FileMode := 2; //set to read/write
      exit;
  end else if FileSz > 128000 then
      FileSz := 128000; //headers will be smaller than this: do not deplete virtual memory
  //lHdrEnd := false;
  GetMem( lCharRA, FileSz+1 );
  BlockRead(fp, lCharRA^, FileSz, linpos);
  if lInPos <> FileSz then showmessage('Disk error: Unable to read full input file.');
  linPos := 1;
  FileSz := FileSize(fp);
  CloseFile(fp);
  FileMode := 2; //set to read/write
  if (chr(lCharRA[linPos])='V') and (chr(lCharRA[linPos+1])='-') and (chr(lCharRA[linPos+2])='d') and (chr(lCharRA[linPos+3]) = 'a')and (chr(lCharRA[linPos+4])='t') and (chr(lCharRA[linPos+5]) = 'a') then
  else
      goto 333;
  lLen := 2;
  repeat
        inc(lLen);
  until (lLen >= FileSz) or ( (lCharRA[lLen-2] = ($0A)) and (lCharRA[lLen-1] = ($0C)) and  (lCharRA[lLen] = ($0A)) ) ;
  lnImages := 0;
  lDynStr := '';
  repeat //for each image
    Clear_Dicom_Data(lDicomData);
    lDicomData.ImageStart := 0;
    lDicomData.little_endian := 0;//is Vista always big endian?
    lDicomData.XYZdim[1] := 1;
    lDicomData.Allocbits_per_pixel := 8;
    lSliceTimeReported := false;
    lHdrOK := true;
    lEdges := false; //edges are not stored as std 2D images
    repeat //for all tags in image
        readVistaStr;
        if lUpCaseValue = '' then
           lDynStr := lDynStr+lUpCaseStr+kCR
        else
            lDynStr := lDynStr+lUpCaseStr+':'+ lUpCaseValue+kCR;
        if (lUpCaseStr ='NCOLUMNS') then lDicomData.XYZdim[1] :=  round(readVistaFloat);
        if (lUpCaseStr ='NROWS') then lDicomData.XYZdim[2] :=  round(readVistaFloat);
        if (lUpCaseStr ='NEDGES') then lEdges := true; //This is a edge file, not a 2d image dataset
        if (lUpCaseStr ='NBANDS') then lDicomData.XYZdim[3] :=  round(readVistaFloat);
        if lUpCaseStr ='DATA'then lDicomData.ImageStart := round(readVistaFloat);
        if (lUpCaseStr ='COLOR_INTERP') and (lUpCaseValue = 'RGB')then lDicomData.SamplesPerPixel := 3;
        if lUpCaseStr ='VOXEL'then readVistaFloats(lDicomData.XYZmm[1],lDicomData.XYZmm[2],lDicomData.XYZmm[3] );
        if lUpCaseStr ='CA'then readVistaFloats(lXori,lYori,lZori);
        if lUpCaseStr ='SLICE_TIME' then lSliceTimeReported := true;
        //if lUpCaseStr ='CONDITION' then showmessage(inttostr(length(lUpCaseValue)));
        if lUpcaseStr = 'REPN' then begin
           if lUpCaseValue = 'BIT' then lDicomData.Allocbits_per_pixel := 1
           else if lUpCaseValue = 'DOUBLE' then lDicomData.Allocbits_per_pixel := 64
           else if lUpCaseValue = 'FLOAT' then lDicomData.Allocbits_per_pixel := 32
           else if lUpCaseValue = 'LONG' then lDicomData.Allocbits_per_pixel := 32
           else if lUpCaseValue = 'SBYTE' then lDicomData.Allocbits_per_pixel := 8
           else if lUpCaseValue = 'SHORT' then lDicomData.Allocbits_per_pixel := 16
           else lDicomData.Allocbits_per_pixel := 8;
           if (lUpCaseValue = 'DOUBLE') or (lUpCaseValue = 'FLOAT') then
              lDicomData.Float := true;
        end //repn
        //lDynStr := lDynStr + lUpCaseStr+lUpCaseValue+kCR;
    until ((lUpcaseStr = '}') and (lDicomData.XYZdim[1] >1)) or ((lINPos+2) >= lLen);
    if (lDicomData.XYZdim[1] > 1) and (not lEdges) then begin //valid image, e.g. not edge description
        if lDicomData.SamplesPerPixel = 3 then begin //RGB
           lDicomData.XYZdim[3] := (lDicomData.XYZdim[3] div 3);
           lDicomData.PlanarConfig := 1;
        end; //RGB image
        lDicomData.ImageStart := lLen + lDicomData.ImageStart;
        lDicomData.Storedbits_per_pixel := lDicomData.Allocbits_per_pixel;
        lImageFormatOK := true;
        lDicomData.XYZori[1] := round(lXori);
        if lYori <> 0 then //analyze are vertically inverted, so we need to flip this axis
           lDicomData.XYZori[2] := lDicomData.XYZdim[2]-round(lYori);
        lDicomData.XYZori[3] := round(lZori);
        if lnImages < kMaxImages then begin
           inc(lnImages);