From 0c5e6d80ccb2c93f1d4c7188ae4a9cc9da5942fe Mon Sep 17 00:00:00 2001 From: Chris Rorden Date: Wed, 22 Feb 2017 07:06:08 -0500 Subject: [PATCH] Threshold Detection http://www.mccauslandcenter.sc.edu/mricrogl/beta-features --- .DS_Store | Bin 14340 -> 14340 bytes .../Contents/Resources/script/basic.gls | 11 +- README.md | 4 +- mainunit.lfm | 1 + nifti_hdr.pas | 1 + nii_reslice.pas | 795 ++++-------------- prefs.pas | 23 +- simplelaz.lpi | 398 +++++---- texture2raycast.pas | 40 +- texture_3d_unit.pas | 2 - 10 files changed, 390 insertions(+), 885 deletions(-) diff --git a/.DS_Store b/.DS_Store index cbdaae51a491c51aa789a701e36306ea0c010638..ec4b42d4d29d6d70613d938c768b53814a41ab32 100644 GIT binary patch delta 97 zcmV-n0G|JZaD;HM;T4n2CKQu9AeR~!dwVf4H7p=4H+?K1dpS5XEFd>FH+_8&eUp(V zH?!Ur1OW)q2?^8@4-gKM{tyAeJmh*F)%VTEFd>9Hhp~%eFz8}DK$Y& zS#W)dqNAjxrl+XG#gn=g9|6X*+7<@^2>A&K)DaI54wL>67PJ2suL1%2v#%Ok1+xnv Pw+*ww9Yq4O2Qd8yPM#-) diff --git a/MRIcroGL.app/Contents/Resources/script/basic.gls b/MRIcroGL.app/Contents/Resources/script/basic.gls index 95bf6e9..9026a32 100755 --- a/MRIcroGL.app/Contents/Resources/script/basic.gls +++ b/MRIcroGL.app/Contents/Resources/script/basic.gls @@ -1,9 +1,12 @@ begin - LOADIMAGE('avg152T1'); + //LOADIMAGE('avg152T1'); //RESETDEFAULTS; - //LOADIMAGE('mni152_2009bet'); + LOADIMAGE('/Users/rorden/Documents/OSX/MRIcroGL/mni152_2009_256'); //CLIPAZIMUTHELEVATION(0.3, 0, 130); - //OVERLAYLOAD('motor'); - //OVERLAYMINMAX(1, 2.6, 4); + OVERLAYLOAD('/Users/rorden/msk/m'); + //OVERLAYLOAD('/Users/rorden/msk/neg'); + //OVERLAYLOAD('/Users/rorden/msk/t'); +OVERLAYMINMAX(1, 2.0, 2); +//OVERLAYMINMAX(1, -2.0, -2); end. diff --git a/README.md b/README.md index 680f34a..9dde594 100755 --- a/README.md +++ b/README.md @@ -10,8 +10,10 @@ http://www.mccauslandcenter.sc.edu/mricrogl/ ##### Recent Versions -30-January-2017 +7-February-2017 - Looks better on Linux high-DPI screens + - [Cubic b-spline interpolation](http://www.mccauslandcenter.sc.edu/mricrogl/beta-features). + - [Threshold Detection](http://www.mccauslandcenter.sc.edu/mricrogl/beta-features). 30-September-2016 - Ensure colorbars show active overlay colors. Improvements for macOS 10.11. 6-June-2016 diff --git a/mainunit.lfm b/mainunit.lfm index dca56a1..2d9d5b6 100755 --- a/mainunit.lfm +++ b/mainunit.lfm @@ -1095,6 +1095,7 @@ object GLForm1: TGLForm1 ShowHint = True TabOrder = 0 OnDrawCell = StringGrid1DrawCell + OnEditingDone = StringGrid1Exit OnExit = StringGrid1Exit OnKeyPress = StringGrid1KeyPress OnMouseDown = StringGrid1MouseDown diff --git a/nifti_hdr.pas b/nifti_hdr.pas index c195ab8..6c69e95 100755 --- a/nifti_hdr.pas +++ b/nifti_hdr.pas @@ -20,6 +20,7 @@ TMRIcroHdr = record //Next: analyze Format Header structure NIFTIhdr : TNIFTIhdr; AutoBalMinUnscaled,AutoBalMaxUnscaled ,WindowScaledMin,WindowScaledMax + //,GlMaxNegUnscaledS,GlMinPosUnscaledS //<- used for thresholded images ,GlMinUnscaledS,GlMaxUnscaledS,Zero8Bit,Slope8bit: double; //brightness and contrast NIfTItransform,DiskDataNativeEndian,UsesCustomPalette,RGB,LutFromZero,LutVisible: boolean; HdrFileName,ImgFileName: string; diff --git a/nii_reslice.pas b/nii_reslice.pas index f7e4095..b5af9b6 100755 --- a/nii_reslice.pas +++ b/nii_reslice.pas @@ -6,7 +6,6 @@ interface uses Dialogs, nii_mat,define_types,sysutils,prefs,nifti_hdr,Texture_3D_Unit, nifti_types; -//function Reslice2Targ (lSrcName,lTargetName,lDestName: string; lPrefs: TPrefs):string; procedure NIFTIhdr_UnswapImg (var lHdr: TMRIcroHdr; var lImgBuffer: byteP); //ensures image data is in native space procedure NIFTIhdr_MinMaxImg (var lHdr: TMRIcroHdr; var lImgBuffer: byteP); //ensures image data is in native space procedure Int32ToFloat (var lHdr: TMRIcroHdr; var lImgBuffer: byteP); @@ -171,74 +170,61 @@ procedure NIFTIhdr_MinMaxImg (var lHdr: TMRIcroHdr; var lImgBuffer: byteP); //en l32f: SingleP; l16Buf : SmallIntP; begin - -(* if lHdr.DiskDataNativeEndian then exit; - case lHdr.NIFTIhdr.datatype of - kDT_UNSIGNED_CHAR : begin - exit; - end; - kDT_SIGNED_SHORT,kDT_SIGNED_INT,kDT_FLOAT: ;//supported format - else begin - Msg('niftiutil UnSwapImg error: datatype not supported.'); - exit; - end; - end; //case *) lImgSamples := lHdr.NIFTIhdr.Dim[1] *lHdr.NIFTIhdr.Dim[2]*lHdr.NIFTIhdr.Dim[3]; if lImgSamples < 1 then exit; case lHdr.NIFTIhdr.datatype of kDT_UNSIGNED_CHAR: begin - lMini := lImgBuffer^[1]; - lMaxi := lImgBuffer^[1]; - for lInc := 1 to lImgSamples do begin - if lImgBuffer^[lInc] > lMaxi then - lMaxi := lImgBuffer^[lInc]; - if lImgBuffer^[lInc] < lMini then - lMini := lImgBuffer^[lInc]; - end; - lHdr.GlMinUnscaledS := lMini; - lHdr.GlMaxUnscaledS := lMaxi; - end; //l16i + lMini := lImgBuffer^[1]; + lMaxi := lImgBuffer^[1]; + for lInc := 1 to lImgSamples do begin + if lImgBuffer^[lInc] > lMaxi then + lMaxi := lImgBuffer^[lInc]; + if lImgBuffer^[lInc] < lMini then + lMini := lImgBuffer^[lInc]; + end; + lHdr.GlMinUnscaledS := lMini; + lHdr.GlMaxUnscaledS := lMaxi; + end; //8ui kDT_SIGNED_SHORT: begin - l16Buf := SmallIntP(lImgBuffer {lHdr.ImgBuffer} ); - lMini := l16Buf^[1]; - lMaxi := l16Buf^[1]; - for lInc := 1 to lImgSamples do begin - if l16Buf^[lInc] > lMaxi then - lMaxi := l16Buf^[lInc]; - if l16Buf^[lInc] < lMini then - lMini := l16Buf^[lInc]; - end; - lHdr.GlMinUnscaledS := lMini; - lHdr.GlMaxUnscaledS := lMaxi; - end; //l16i - kDT_SIGNED_INT: begin - l32i := LongIntP(@lImgBuffer^); - lMini := l32i^[1]; - lMaxi := l32i^[1]; - for lInc := 1 to lImgSamples do begin - if l32i^[lInc] > lMaxi then - lMaxi := l32i^[lInc]; - if l32i^[lInc] < lMini then - lMini := l32i^[lInc]; - end; - lHdr.GlMinUnscaledS := lMini; - lHdr.GlMaxUnscaledS := lMaxi; - - end; //32i - kDT_FLOAT: begin - l32f := SingleP(@lImgBuffer^); - lMins := l32f^[1]; - lMaxs := l32f^[1]; - for lInc := 1 to lImgSamples do begin - if l32f^[lInc] > lMaxs then - lMaxs := l32f^[lInc]; - if l32f^[lInc] < lMins then - lMins := l32f^[lInc]; - end; - lHdr.GlMinUnscaledS := lMins; - lHdr.GlMaxUnscaledS := lMaxs; - end;//32i + l16Buf := SmallIntP(lImgBuffer {lHdr.ImgBuffer} ); + lMini := l16Buf^[1]; + lMaxi := l16Buf^[1]; + for lInc := 1 to lImgSamples do begin + if l16Buf^[lInc] > lMaxi then + lMaxi := l16Buf^[lInc]; + if l16Buf^[lInc] < lMini then + lMini := l16Buf^[lInc]; + end; + lHdr.GlMinUnscaledS := lMini; + lHdr.GlMaxUnscaledS := lMaxi; + end; //16i + kDT_SIGNED_INT: begin + l32i := LongIntP(@lImgBuffer^); + lMini := l32i^[1]; + lMaxi := l32i^[1]; + for lInc := 1 to lImgSamples do begin + if l32i^[lInc] > lMaxi then + lMaxi := l32i^[lInc]; + if l32i^[lInc] < lMini then + lMini := l32i^[lInc]; + end; + lHdr.GlMinUnscaledS := lMini; + lHdr.GlMaxUnscaledS := lMaxi; + end; //32i + kDT_FLOAT: begin + l32f := SingleP(@lImgBuffer^); + lMins := l32f^[1]; + lMaxs := l32f^[1]; + for lInc := 1 to lImgSamples do begin + if l32f^[lInc] > lMaxs then + lMaxs := l32f^[lInc]; + if l32f^[lInc] < lMins then + lMins := l32f^[lInc]; + end; + lHdr.GlMinUnscaledS := lMins; + lHdr.GlMaxUnscaledS := lMaxs; + end;//32f end; //case end; @@ -442,15 +428,119 @@ procedure Float32RemoveNAN (var lHdr: TMRIcroHdr; var lImgBuffer: byteP); lInVox := lHdr.NIFTIhdr.dim[1] * lHdr.NIFTIhdr.dim[2] * lHdr.NIFTIhdr.dim[3]; l32Buf := SingleP(lImgBuffer ); for lI := 1 to lInVox do - if specialsingle(l32Buf^[lI]) then l32Buf^[lI] :=0.0; + if specialsingle(l32Buf^[lI]) then l32Buf^[lI] :=0.0; end;//Float32RemoveNAN +procedure NIFTIhdr_ThreshFind (var lHdr: TMRIcroHdr; var lImgBuffer: byteP; out lMinPos, lMaxNeg: single); +var + lInc,lImgSamples : integer; + lMax, lMin: single; + l32f: SingleP; +begin + lImgSamples := lHdr.NIFTIhdr.Dim[1] *lHdr.NIFTIhdr.Dim[2]*lHdr.NIFTIhdr.Dim[3]; + lMinPos := 0; + lMaxNeg := 0; + if (lImgSamples < 1) or (lHdr.NIFTIhdr.datatype <> kDT_FLOAT) or (not gPrefs.ThresholdDetection) then + exit; + lMinPos := MaxInt; + lMaxNeg := -MaxInt; + lMin := 0; + lMax := 0; + case lHdr.NIFTIhdr.datatype of + kDT_FLOAT: begin + l32f := SingleP(@lImgBuffer^); + lMax :=l32f^[1]; + lMin := lMax; + for lInc := 1 to lImgSamples do begin + if (l32f^[lInc] > lMax) then + lMax := l32f^[lInc]; + if (l32f^[lInc] > 0) and (l32f^[lInc] < lMinPos) then + lMinPos := l32f^[lInc]; + if (l32f^[lInc] < lMin) then + lMin := l32f^[lInc]; + if (l32f^[lInc] < 0) and (l32f^[lInc] > lMaxNeg) then + lMaxNeg := l32f^[lInc]; + end; //for each sample + end;//32f + end; //case + if (lMaxNeg = -MaxInt) or (lMaxNeg = lMin) then lMaxNeg := 0; + if (lMinPos = MaxInt) or (lMinPos = lMax) then lMinPos := 0; + if (lHdr.NIFTIhdr.datatype <> kDT_FLOAT) and (lMinPos = 1) then + lMinPos := 0; //minimum integer is 1 + if (lHdr.NIFTIhdr.datatype <> kDT_FLOAT) and (lMaxNeg = -1) then + lMaxNeg := 0; //minimum integer is 1 + //scale v := (v * gTexture3D.NIFTIhdr.scl_slope)+gTexture3D.NIFTIhdr.scl_inter; + if ((lMinPos*lHdr.NIFTIhdr.scl_slope)+lHdr.NIFTIhdr.scl_inter) < 1.0 then + lMinPos := 0; + if ((lMaxNeg*lHdr.NIFTIhdr.scl_slope)+lHdr.NIFTIhdr.scl_inter) > -1.0 then + lMaxNeg := 0; + //GLForm1.caption := ( format('%g %g',[lMinPos, lMaxNeg]) ); +end; + +procedure NIFTIhdr_ThreshPos (var lHdr: TMRIcroHdr; var lImgBuffer: byteP; lMinPos: single); //ensures image data is in native space +var + lInc,lImgSamples : integer; + lMinPosDiv2: single; + //l32i : LongIntP; + l32f: SingleP; + //l16Buf : SmallIntP; +begin + lImgSamples := lHdr.NIFTIhdr.Dim[1] *lHdr.NIFTIhdr.Dim[2]*lHdr.NIFTIhdr.Dim[3]; + if (lImgSamples < 1) or (lHdr.NIFTIhdr.datatype <> kDT_FLOAT) or (not gPrefs.ThresholdDetection) or (lMinPos <= 0) then + exit; + lMinPosDiv2 := lMinPos / 2.0; + case lHdr.NIFTIhdr.datatype of + kDT_FLOAT: begin + l32f := SingleP(@lImgBuffer^); + for lInc := 1 to lImgSamples do begin + if (l32f^[lInc] > 0) then begin + if (l32f^[lInc] < lMinPosDiv2) then + l32f^[lInc] := 0 + else if (l32f^[lInc] < lMinPos) then + l32f^[lInc] := lMinPos; + end; //pos sample + end; //for each sample + end;//32f + end; //case +end; //NIFTIhdr_ThreshPos() + +procedure NIFTIhdr_ThreshNeg (var lHdr: TMRIcroHdr; var lImgBuffer: byteP; lMaxNeg: single); //ensures image data is in native space +var + lInc,lImgSamples : integer; + lMaxNegDiv2: single; + //l32i : LongIntP; + l32f: SingleP; + //l16Buf : SmallIntP; +begin + lImgSamples := lHdr.NIFTIhdr.Dim[1] *lHdr.NIFTIhdr.Dim[2]*lHdr.NIFTIhdr.Dim[3]; + if (lImgSamples < 1) or (lHdr.NIFTIhdr.datatype <> kDT_FLOAT) or (not gPrefs.ThresholdDetection) or (lMaxNeg >= 0) then + exit; + lMaxNegDiv2 := lMaxNeg / 2.0; + case lHdr.NIFTIhdr.datatype of + kDT_FLOAT: begin + l32f := SingleP(@lImgBuffer^); + for lInc := 1 to lImgSamples do begin + if (l32f^[lInc] < 0) then begin + if (l32f^[lInc] > lMaxNegDiv2) then + l32f^[lInc] := 0 + else if (l32f^[lInc] > lMaxNeg) then + l32f^[lInc] := lMaxNeg; + + end; //neg sample + end; //for each sample + end;//32f + end; //case + //GLForm1.caption := ( format('neg %g boost %d',[lMaxNeg, Nx]) ); +end; //NIFTIhdr_ThreshNeg() + + function Reslice2TargCore (var lSrcHdr: TMRIcroHdr; var lSrcBuffer: bytep; var lTargHdr: TNIFTIHdr; var lDestHdr: TMRIcroHdr; lTrilinearInterpolation: boolean; lVolume: integer): string; //output lDestHdr var lPos,lXYs,lXYZs,lXs,lYs,lZs,lXi,lYi,lZi,lX,lY,lZ, lXo,lYo,lZo,lMinY,lMinZ,lMaxY,lMaxZ,lBPP,lXYZ: integer; + lMinPositive, lMaxNegative, lXrM1,lYrM1,lZrM1,lXreal,lYreal,lZreal, lZx,lZy,lZz,lYx,lYy,lYz, lInMinX,lInMinY,lInMinZ, lOutMinX,lOutMinY,lOutMinZ: single; @@ -472,6 +562,8 @@ function Reslice2TargCore (var lSrcHdr: TMRIcroHdr; var lSrcBuffer: bytep; var NIFTIhdr_UnswapImg (lSrcHdr,lSrcBuffer); //ensures image data is in native byteorder Float32RemoveNAN(lSrcHdr,lSrcBuffer); RGB24ToByte (lSrcHdr, lSrcBuffer,lVolume); + NIFTIhdr_ThreshFind (lSrcHdr, lSrcBuffer, lMinPositive, lMaxNegative); //ensures image data is in native space + //AbsFloat(lSrcHdr, lSrcBuffer); case lSrcHdr.NIFTIhdr.datatype of kDT_UNSIGNED_CHAR : lBPP := 1; @@ -682,24 +774,11 @@ function Reslice2TargCore (var lSrcHdr: TMRIcroHdr; var lSrcBuffer: bytep; var end;//z end; - //release lookup tables freemem(lXx); freemem(lXy); freemem(lXz); - //check to see if image is empty... - (*lPos := 1; - case lSrcHdr.NIFTIhdr.datatype of - kDT_UNSIGNED_CHAR : while (lPos <= (lX*lY*lZ)) and (l8i^[lPos] = 0) do inc(lPos); - kDT_SIGNED_SHORT: while (lPos <= (lX*lY*lZ)) and (l16i^[lPos] = 0) do inc(lPos); - kDT_SIGNED_INT:while (lPos <= (lX*lY*lZ)) and (l32i^[lPos] = 0) do inc(lPos); - kDT_FLOAT: while (lPos <= (lX*lY*lZ)) and (l32f^[lPos] = 0) do inc(lPos); - end; //case - if lPos <= (lX*lY*lZ) then //image not empty - //Msg('Overlap') - //result := SaveNIfTICore (lDestName, lBuffAligned, kNIIImgOffset+1, lDestHdr, lPrefs,lByteSwap); - else - Msg('Overlay image does not overlap with background image.'); *) + if not lOverlap then Msg('Overlay image does not overlap with background image.'); //Freemem(lBuffUnaligned); @@ -710,6 +789,8 @@ function Reslice2TargCore (var lSrcHdr: TMRIcroHdr; var lSrcBuffer: bytep; var kDT_SIGNED_INT:lDestHdr.ImgBufferBPP :=4; kDT_FLOAT: lDestHdr.ImgBufferBPP :=4; end; //case + NIFTIhdr_ThreshPos(lDestHdr,lDestHdr.ImgBuffer, lMinPositive); + NIFTIhdr_ThreshNeg(lDestHdr,lDestHdr.ImgBuffer, lMaxNegative); NIFTIhdr_MinMaxImg(lDestHdr,lDestHdr.ImgBuffer);//set global min/max result := 'OK'; end; @@ -725,562 +806,4 @@ function Reslice2Targ (lSrcName: string; var lTargHdr: TNIFTIHdr; var lDestHdr: Freemem(lSrcBuffer); end; -(*function Reslice2Targ (lSrcName: string; var lTargHdr: TNIFTIHdr; var lDestHdr: TMRIcroHdr; lTrilinearInterpolation: boolean; lVolume: integer): string; -var - lPos,lXYs,lXYZs,lXs,lYs,lZs,lXi,lYi,lZi,lX,lY,lZ, - lXo,lYo,lZo,lMinY,lMinZ,lMaxY,lMaxZ,lBPP,lXYZ: integer; - lXrM1,lYrM1,lZrM1,lXreal,lYreal,lZreal, - lZx,lZy,lZz,lYx,lYy,lYz, - lInMinX,lInMinY,lInMinZ, lOutMinX,lOutMinY,lOutMinZ: single; - lXx,lXy,lXz: Singlep0; - l32fs,l32f : SingleP; - l32is,l32i : LongIntP; - l16is,l16i : SmallIntP; - l8i,l8is,lSrcBuffer: bytep; - lMat: TMatrix; - lSrcHdr: TMRIcroHdr; - lOverlap: boolean; -begin - result := ''; - lOverlap := false; - if not NIFTIhdr_LoadImg (lSrcName, lSrcHdr, lSrcBuffer,lVolume) then exit; - lOverlap := false; - //convert 32-bit int to 32-bit float.... - Int32ToFloat (lSrcHdr, lSrcBuffer); - Float64ToFloat32(lSrcHdr, lSrcBuffer); - NIFTIhdr_UnswapImg (lSrcHdr,lSrcBuffer); //ensures image data is in native byteorder - Float32RemoveNAN(lSrcHdr,lSrcBuffer); - RGB24ToByte (lSrcHdr, lSrcBuffer,lVolume); - //AbsFloat(lSrcHdr, lSrcBuffer); - case lSrcHdr.NIFTIhdr.datatype of - kDT_UNSIGNED_CHAR : lBPP := 1; - kDT_SIGNED_SHORT: lBPP := 2; - kDT_SIGNED_INT:lBPP := 4; - kDT_FLOAT: lBPP := 4; - kDT_RGB: lBPP := 1; - else begin - Msg('NII reslice error: datatype not supported.'); - exit; - end; - end; //case - lMat := Voxel2Voxel (lTargHdr,lSrcHdr.NIFTIhdr); - lDestHdr {.NIFTIhdr} := lSrcHdr {.NIFTIhdr}; //destination has the comments and voxel BPP of source - CopyHdrMat(lTargHdr,lDestHdr.NIFTIhdr);//destination has dimensions and rotations of destination - lXs := lSrcHdr.NIFTIhdr.Dim[1]; - lYs := lSrcHdr.NIFTIhdr.Dim[2]; - lZs := lSrcHdr.NIFTIhdr.Dim[3]; - lXYs:=lXs*lYs; //slicesz - lXYZs := lXYs*lZs; - lX := lDestHdr.NIFTIhdr.Dim[1]; - lY := lDestHdr.NIFTIhdr.Dim[2]; - lZ := lDestHdr.NIFTIhdr.Dim[3]; - lDestHdr.NIFTIhdr.Dim[4] := 1; - //load dataset - NIFTIhdr_UnswapImg(lSrcHdr, lSrcBuffer);//interpolation requires data is in native endian - { We will set min/max after scaling.. - NIFTIhdr_MinMaxImg(lSrcHdr, lSrcBuffer); - lDestHdr.GlMinUnscaledS := lSrcHdr.GlMinUnscaledS; - lDestHdr.GlMaxUnscaledS := lSrcHdr.GlMaxUnscaledS; } - l8is := (@lSrcBuffer^); - - GetMem(lDestHdr.ImgBufferUnaligned ,(lBPP*lX*lY*lZ)+15); - {$IFDEF FPC} - lDestHdr.ImgBuffer := Align(lDestHdr.ImgBufferUnaligned,16); // not commented - check this - {$ELSE} - lDestHdr.ImgBuffer := ByteP($fffffff0 and (integer(lDestHdr.ImgBufferUnaligned)+15)); - {$ENDIF} - //lPos := 1; - case lSrcHdr.NIFTIhdr.datatype of - kDT_UNSIGNED_CHAR : l8i := @lDestHdr.ImgBuffer^; - kDT_SIGNED_SHORT: l16i := SmallIntP(@lDestHdr.ImgBuffer^ ); - kDT_SIGNED_INT:l32i := LongIntP(@lDestHdr.ImgBuffer^); - kDT_FLOAT: l32f := SingleP(@lDestHdr.ImgBuffer^ ); - end; //case - case lSrcHdr.NIFTIhdr.datatype of - //kDT_UNSIGNED_CHAR : l8is := l8is; - kDT_SIGNED_SHORT: l16is := SmallIntP(l8is ); - kDT_SIGNED_INT:l32is := LongIntP(l8is ); - kDT_FLOAT: l32fs := SingleP(l8is ); - end; //case - //next clear image - - case lSrcHdr.NIFTIhdr.datatype of - kDT_UNSIGNED_CHAR : for lPos := 1 to (lX*lY*lZ) do l8i^[lPos] := 0; - kDT_SIGNED_SHORT: for lPos := 1 to (lX*lY*lZ) do l16i^[lPos] := 0; - kDT_SIGNED_INT:for lPos := 1 to (lX*lY*lZ) do l32i^[lPos] := 0; - kDT_FLOAT: for lPos := 1 to (lX*lY*lZ) do l32f^[lPos] := 0; - end; //case - - //now we can apply the transforms... - //build lookup table - speed up inner loop - getmem(lXx, lX*sizeof(single)); - getmem(lXy, lX*sizeof(single)); - getmem(lXz, lX*sizeof(single)); - for lXi := 0 to (lX-1) do begin - lXx^[lXi] := lXi*lMat.matrix[1][1]; - lXy^[lXi] := lXi*lMat.matrix[2][1]; - lXz^[lXi] := lXi*lMat.matrix[3][1]; - end; - lPos := 0; -if lTrilinearInterpolation then begin - for lZi := 0 to (lZ-1) do begin - //these values are the same for all voxels in the slice - // compute once per slice - lZx := lZi*lMat.matrix[1][3]; - lZy := lZi*lMat.matrix[2][3]; - lZz := lZi*lMat.matrix[3][3]; - for lYi := 0 to (lY-1) do begin - //these values change once per row - // compute once per row - lYx := lYi*lMat.matrix[1][2]; - lYy := lYi*lMat.matrix[2][2]; - lYz := lYi*lMat.matrix[3][2]; - for lXi := 0 to (lX-1) do begin - //compute each column - inc(lPos); - - lXreal := (lXx^[lXi]+lYx+lZx+lMat.matrix[1][4]); - lYreal := (lXy^[lXi]+lYy+lZy+lMat.matrix[2][4]); - lZreal := (lXz^[lXi]+lYz+lZz+lMat.matrix[3][4]); - //need to test Xreal as -0.01 truncates to zero - if (lXreal >= 0) and (lYreal >= 0{1}) and (lZreal >= 0{1}) and - (lXreal < (lXs -1)) and (lYreal < (lYs -1) ) and (lZreal < (lZs -1)) - then begin - //compute the contribution for each of the 8 source voxels - //nearest to the target - lOverlap := true; - lXo := trunc(lXreal); - lYo := trunc(lYreal); - lZo := trunc(lZreal); - lXreal := lXreal-lXo; - lYreal := lYreal-lYo; - lZreal := lZreal-lZo; - lXrM1 := 1-lXreal; - lYrM1 := 1-lYreal; - lZrM1 := 1-lZreal; - lMinY := lYo*lXs; - lMinZ := lZo*lXYs; - lMaxY := lMinY+lXs; - lMaxZ := lMinZ+lXYs; - inc(lXo);//images incremented from 1 not 0 - case lSrcHdr.NIFTIhdr.datatype of - kDT_UNSIGNED_CHAR : begin// l8is := l8is; - l8i^[lPos] := - round ( - {all min} ( (lXrM1*lYrM1*lZrM1)*l8is^[lXo+lMinY+lMinZ]) - {x+1}+((lXreal*lYrM1*lZrM1)*l8is^[lXo+1+lMinY+lMinZ]) - {y+1}+((lXrM1*lYreal*lZrM1)*l8is^[lXo+lMaxY+lMinZ]) - {z+1}+((lXrM1*lYrM1*lZreal)*l8is^[lXo+lMinY+lMaxZ]) - {x+1,y+1}+((lXreal*lYreal*lZrM1)*l8is^[lXo+1+lMaxY+lMinZ]) - {x+1,z+1}+((lXreal*lYrM1*lZreal)*l8is^[lXo+1+lMinY+lMaxZ]) - {y+1,z+1}+((lXrM1*lYreal*lZreal)*l8is^[lXo+lMaxY+lMaxZ]) - {x+1,y+1,z+1}+((lXreal*lYreal*lZreal)*l8is^[lXo+1+lMaxY+lMaxZ]) ); - end; - kDT_SIGNED_SHORT: begin - l16i^[lPos] := - round ( - {all min} ( (lXrM1*lYrM1*lZrM1)*l16is^[lXo+lMinY+lMinZ]) - {x+1}+((lXreal*lYrM1*lZrM1)*l16is^[lXo+1+lMinY+lMinZ]) - {y+1}+((lXrM1*lYreal*lZrM1)*l16is^[lXo+lMaxY+lMinZ]) - {z+1}+((lXrM1*lYrM1*lZreal)*l16is^[lXo+lMinY+lMaxZ]) - {x+1,y+1}+((lXreal*lYreal*lZrM1)*l16is^[lXo+1+lMaxY+lMinZ]) - {x+1,z+1}+((lXreal*lYrM1*lZreal)*l16is^[lXo+1+lMinY+lMaxZ]) - {y+1,z+1}+((lXrM1*lYreal*lZreal)*l16is^[lXo+lMaxY+lMaxZ]) - {x+1,y+1,z+1}+((lXreal*lYreal*lZreal)*l16is^[lXo+1+lMaxY+lMaxZ]) ); - end; - kDT_SIGNED_INT:begin - l32i^[lPos] := - round ( - {all min} ( (lXrM1*lYrM1*lZrM1)*l32is^[lXo+lMinY+lMinZ]) - {x+1}+((lXreal*lYrM1*lZrM1)*l32is^[lXo+1+lMinY+lMinZ]) - {y+1}+((lXrM1*lYreal*lZrM1)*l32is^[lXo+lMaxY+lMinZ]) - {z+1}+((lXrM1*lYrM1*lZreal)*l32is^[lXo+lMinY+lMaxZ]) - {x+1,y+1}+((lXreal*lYreal*lZrM1)*l32is^[lXo+1+lMaxY+lMinZ]) - {x+1,z+1}+((lXreal*lYrM1*lZreal)*l32is^[lXo+1+lMinY+lMaxZ]) - {y+1,z+1}+((lXrM1*lYreal*lZreal)*l32is^[lXo+lMaxY+lMaxZ]) - {x+1,y+1,z+1}+((lXreal*lYreal*lZreal)*l32is^[lXo+1+lMaxY+lMaxZ]) ); - end; - kDT_FLOAT: begin //note - we do not round results - all intensities might be frational... - l32f^[lPos] := - ( - {all min} ( (lXrM1*lYrM1*lZrM1)*l32fs^[lXo+lMinY+lMinZ]) - {x+1}+((lXreal*lYrM1*lZrM1)*l32fs^[lXo+1+lMinY+lMinZ]) - {y+1}+((lXrM1*lYreal*lZrM1)*l32fs^[lXo+lMaxY+lMinZ]) - {z+1}+((lXrM1*lYrM1*lZreal)*l32fs^[lXo+lMinY+lMaxZ]) - {x+1,y+1}+((lXreal*lYreal*lZrM1)*l32fs^[lXo+1+lMaxY+lMinZ]) - {x+1,z+1}+((lXreal*lYrM1*lZreal)*l32fs^[lXo+1+lMinY+lMaxZ]) - {y+1,z+1}+((lXrM1*lYreal*lZreal)*l32fs^[lXo+lMaxY+lMaxZ]) - {x+1,y+1,z+1}+((lXreal*lYreal*lZreal)*l32fs^[lXo+1+lMaxY+lMaxZ]) ); - end; - end; //case - - end; //if voxel is in source image's bounding box - end;//z - end;//y - end;//z -end else begin //if trilinear, else nearest neighbor - - for lZi := 0 to (lZ-1) do begin - //these values are the same for all voxels in the slice - // compute once per slice - lZx := lZi*lMat.matrix[1][3]; - lZy := lZi*lMat.matrix[2][3]; - lZz := lZi*lMat.matrix[3][3]; - for lYi := 0 to (lY-1) do begin - //these values change once per row - // compute once per row - lYx := lYi*lMat.matrix[1][2]; - lYy := lYi*lMat.matrix[2][2]; - lYz := lYi*lMat.matrix[3][2]; - for lXi := 0 to (lX-1) do begin - //compute each column - inc(lPos); - lXo := round(lXx^[lXi]+lYx+lZx+lMat.matrix[1][4]); - lYo := round(lXy^[lXi]+lYy+lZy+lMat.matrix[2][4]); - lZo := round(lXz^[lXi]+lYz+lZz+lMat.matrix[3][4]); - //if lZo <> 0 then - // fx(lZo); - //need to test Xreal as -0.01 truncates to zero - if (lXo >= 0) and (lYo >= 0{1}) and (lZo >= 0{1}) and - (lXo < (lXs -1)) and (lYo < (lYs -1) ) and (lZo < (lZs {-1})) - then begin - lOverlap := true; - inc(lXo);//images incremented from 1 not 0 - lYo := lYo*lXs; - lZo := lZo*lXYs; - case lSrcHdr.NIFTIhdr.datatype of - kDT_UNSIGNED_CHAR : // l8is := l8is; - l8i^[lPos] :=l8is^[lXo+lYo+lZo]; - kDT_SIGNED_SHORT: l16i^[lPos] := l16is^[lXo+lYo+lZo]; - kDT_SIGNED_INT:l32i^[lPos] := l32is^[lXo+lYo+lZo]; - kDT_FLOAT: l32f^[lPos] := l32fs^[lXo+lYo+lZo]; - end; //case - end; //if voxel is in source image's bounding box - end;//z - end;//y - end;//z - -end; - Freemem(lSrcBuffer); - //release lookup tables - freemem(lXx); - freemem(lXy); - freemem(lXz); - //check to see if image is empty... - if not lOverlap then - Msg('Overlay image does not overlap with background image.'); - //Freemem(lBuffUnaligned); - lDestHdr.ImgBufferItems := lX*lY*lZ; - case lSrcHdr.NIFTIhdr.datatype of - kDT_UNSIGNED_CHAR :lDestHdr.ImgBufferBPP :=1; - kDT_SIGNED_SHORT: lDestHdr.ImgBufferBPP :=2; - kDT_SIGNED_INT:lDestHdr.ImgBufferBPP :=4; - kDT_FLOAT: lDestHdr.ImgBufferBPP :=4; - end; //case - NIFTIhdr_MinMaxImg(lDestHdr,lDestHdr.ImgBuffer);//set global min/max - result := 'OK'; -end; *) - -(*function Reslice2Targ (lSrcName: string; var lTargHdr: TNIFTIHdr; var lDestHdr: TMRIcroHdr; lTrilinearInterpolation: boolean; lVolume: integer): string; -var - lPos,lXYs,lXYZs,lXs,lYs,lZs,lXi,lYi,lZi,lX,lY,lZ, - lXo,lYo,lZo,lMinY,lMinZ,lMaxY,lMaxZ,lBPP,lXYZ: integer; - lXrM1,lYrM1,lZrM1,lXreal,lYreal,lZreal, - lZx,lZy,lZz,lYx,lYy,lYz, - lInMinX,lInMinY,lInMinZ, lOutMinX,lOutMinY,lOutMinZ: single; - lXx,lXy,lXz: Singlep0; - l32fs,l32f : SingleP; - l32is,l32i : LongIntP; - l16is,l16i : SmallIntP; - l8i,l8is,lSrcBuffer: bytep; - lMat: TMatrix; - lSrcHdr: TMRIcroHdr; -begin - result := ''; - if not NIFTIhdr_LoadImg (lSrcName, lSrcHdr, lSrcBuffer,lVolume) then exit; - //convert 32-bit int to 32-bit float.... - Int32ToFloat (lSrcHdr, lSrcBuffer); - Float64ToFloat32(lSrcHdr, lSrcBuffer); - //AbsFloat(lSrcHdr, lSrcBuffer); - - case lSrcHdr.NIFTIhdr.datatype of - kDT_UNSIGNED_CHAR : lBPP := 1; - kDT_SIGNED_SHORT: lBPP := 2; - kDT_SIGNED_INT:lBPP := 4; - kDT_FLOAT: lBPP := 4; - else begin - Msg('NII reslice error: datatype not supported.'); - exit; - end; - end; //case - lMat := Voxel2Voxel (lTargHdr,lSrcHdr.NIFTIhdr); - lDestHdr.NIFTIhdr := lSrcHdr.NIFTIhdr; //destination has the comments and voxel BPP of source - //lDestHdr.NIFTIhdr.datatype := lSrcHdr.NIFTIhdr.datatype; - CopyHdrMat(lTargHdr,lDestHdr.NIFTIhdr);//destination has dimensions and rotations of destination - lXs := lSrcHdr.NIFTIhdr.Dim[1]; - lYs := lSrcHdr.NIFTIhdr.Dim[2]; - lZs := lSrcHdr.NIFTIhdr.Dim[3]; - lXYs:=lXs*lYs; //slicesz - lXYZs := lXYs*lZs; - lX := lDestHdr.NIFTIhdr.Dim[1]; - lY := lDestHdr.NIFTIhdr.Dim[2]; - lZ := lDestHdr.NIFTIhdr.Dim[3]; - lDestHdr.NIFTIhdr.Dim[4] := 1; - //load dataset - NIFTIhdr_UnswapImg(lSrcHdr, lSrcBuffer);//interpolation requires data is in native endian - { We will set min/max after scaling.. - NIFTIhdr_MinMaxImg(lSrcHdr, lSrcBuffer); - lDestHdr.GlMinUnscaledS := lSrcHdr.GlMinUnscaledS; - lDestHdr.GlMaxUnscaledS := lSrcHdr.GlMaxUnscaledS; } - l8is := (@lSrcBuffer^); - GetMem(lDestHdr.ImgBufferUnaligned ,(lBPP*lX*lY*lZ)+15); - {$IFDEF FPC} - lDestHdr.ImgBuffer := Align(lDestHdr.ImgBufferUnaligned,16); // not commented - check this - {$ELSE} - lDestHdr.ImgBuffer := ByteP($fffffff0 and (integer(lDestHdr.ImgBufferUnaligned)+15)); - {$ENDIF} - lPos := 1; - case lSrcHdr.NIFTIhdr.datatype of - kDT_UNSIGNED_CHAR : l8i := @lDestHdr.ImgBuffer^; - kDT_SIGNED_SHORT: l16i := SmallIntP(@lDestHdr.ImgBuffer^ ); - kDT_SIGNED_INT:l32i := LongIntP(@lDestHdr.ImgBuffer^); - kDT_FLOAT: l32f := SingleP(@lDestHdr.ImgBuffer^ ); - end; //case - case lSrcHdr.NIFTIhdr.datatype of - //kDT_UNSIGNED_CHAR : l8is := l8is; - kDT_SIGNED_SHORT: l16is := SmallIntP(l8is ); - kDT_SIGNED_INT:l32is := LongIntP(l8is ); - kDT_FLOAT: l32fs := SingleP(l8is ); - end; //case - //next clear image - - case lSrcHdr.NIFTIhdr.datatype of - kDT_UNSIGNED_CHAR : for lPos := 1 to (lX*lY*lZ) do l8i^[lPos] := 0; - kDT_SIGNED_SHORT: for lPos := 1 to (lX*lY*lZ) do l16i^[lPos] := 0; - kDT_SIGNED_INT:for lPos := 1 to (lX*lY*lZ) do l32i^[lPos] := 0; - kDT_FLOAT: for lPos := 1 to (lX*lY*lZ) do l32f^[lPos] := 0; - end; //case - - //now we can apply the transforms... - //build lookup table - speed up inner loop - getmem(lXx, lX*sizeof(single)); - getmem(lXy, lX*sizeof(single)); - getmem(lXz, lX*sizeof(single)); - for lXi := 0 to (lX-1) do begin - lXx^[lXi] := lXi*lMat.matrix[1][1]; - lXy^[lXi] := lXi*lMat.matrix[2][1]; - lXz^[lXi] := lXi*lMat.matrix[3][1]; - end; - lPos := 0; -if lTrilinearInterpolation then begin - for lZi := 0 to (lZ-1) do begin - //these values are the same for all voxels in the slice - // compute once per slice - lZx := lZi*lMat.matrix[1][3]; - lZy := lZi*lMat.matrix[2][3]; - lZz := lZi*lMat.matrix[3][3]; - for lYi := 0 to (lY-1) do begin - //these values change once per row - // compute once per row - lYx := lYi*lMat.matrix[1][2]; - lYy := lYi*lMat.matrix[2][2]; - lYz := lYi*lMat.matrix[3][2]; - for lXi := 0 to (lX-1) do begin - //compute each column - inc(lPos); - - lXreal := (lXx^[lXi]+lYx+lZx+lMat.matrix[1][4]); - lYreal := (lXy^[lXi]+lYy+lZy+lMat.matrix[2][4]); - lZreal := (lXz^[lXi]+lYz+lZz+lMat.matrix[3][4]); - //need to test Xreal as -0.01 truncates to zero - if (lXreal >= 0) and (lYreal >= 0{1}) and (lZreal >= 0{1}) and - (lXreal < (lXs -1)) and (lYreal < (lYs -1) ) and (lZreal < (lZs -1)) - then begin - //compute the contribution for each of the 8 source voxels - //nearest to the target - lXo := trunc(lXreal); - lYo := trunc(lYreal); - lZo := trunc(lZreal); - lXreal := lXreal-lXo; - lYreal := lYreal-lYo; - lZreal := lZreal-lZo; - lXrM1 := 1-lXreal; - lYrM1 := 1-lYreal; - lZrM1 := 1-lZreal; - lMinY := lYo*lXs; - lMinZ := lZo*lXYs; - lMaxY := lMinY+lXs; - lMaxZ := lMinZ+lXYs; - inc(lXo);//images incremented from 1 not 0 - case lSrcHdr.NIFTIhdr.datatype of - kDT_UNSIGNED_CHAR : begin// l8is := l8is; - l8i^[lPos] := - round ( - {all min} ( (lXrM1*lYrM1*lZrM1)*l8is^[lXo+lMinY+lMinZ]) - {x+1}+((lXreal*lYrM1*lZrM1)*l8is^[lXo+1+lMinY+lMinZ]) - {y+1}+((lXrM1*lYreal*lZrM1)*l8is^[lXo+lMaxY+lMinZ]) - {z+1}+((lXrM1*lYrM1*lZreal)*l8is^[lXo+lMinY+lMaxZ]) - {x+1,y+1}+((lXreal*lYreal*lZrM1)*l8is^[lXo+1+lMaxY+lMinZ]) - {x+1,z+1}+((lXreal*lYrM1*lZreal)*l8is^[lXo+1+lMinY+lMaxZ]) - {y+1,z+1}+((lXrM1*lYreal*lZreal)*l8is^[lXo+lMaxY+lMaxZ]) - {x+1,y+1,z+1}+((lXreal*lYreal*lZreal)*l8is^[lXo+1+lMaxY+lMaxZ]) ); - end; - kDT_SIGNED_SHORT: begin - l16i^[lPos] := - round ( - {all min} ( (lXrM1*lYrM1*lZrM1)*l16is^[lXo+lMinY+lMinZ]) - {x+1}+((lXreal*lYrM1*lZrM1)*l16is^[lXo+1+lMinY+lMinZ]) - {y+1}+((lXrM1*lYreal*lZrM1)*l16is^[lXo+lMaxY+lMinZ]) - {z+1}+((lXrM1*lYrM1*lZreal)*l16is^[lXo+lMinY+lMaxZ]) - {x+1,y+1}+((lXreal*lYreal*lZrM1)*l16is^[lXo+1+lMaxY+lMinZ]) - {x+1,z+1}+((lXreal*lYrM1*lZreal)*l16is^[lXo+1+lMinY+lMaxZ]) - {y+1,z+1}+((lXrM1*lYreal*lZreal)*l16is^[lXo+lMaxY+lMaxZ]) - {x+1,y+1,z+1}+((lXreal*lYreal*lZreal)*l16is^[lXo+1+lMaxY+lMaxZ]) ); - end; - kDT_SIGNED_INT:begin - l32i^[lPos] := - round ( - {all min} ( (lXrM1*lYrM1*lZrM1)*l32is^[lXo+lMinY+lMinZ]) - {x+1}+((lXreal*lYrM1*lZrM1)*l32is^[lXo+1+lMinY+lMinZ]) - {y+1}+((lXrM1*lYreal*lZrM1)*l32is^[lXo+lMaxY+lMinZ]) - {z+1}+((lXrM1*lYrM1*lZreal)*l32is^[lXo+lMinY+lMaxZ]) - {x+1,y+1}+((lXreal*lYreal*lZrM1)*l32is^[lXo+1+lMaxY+lMinZ]) - {x+1,z+1}+((lXreal*lYrM1*lZreal)*l32is^[lXo+1+lMinY+lMaxZ]) - {y+1,z+1}+((lXrM1*lYreal*lZreal)*l32is^[lXo+lMaxY+lMaxZ]) - {x+1,y+1,z+1}+((lXreal*lYreal*lZreal)*l32is^[lXo+1+lMaxY+lMaxZ]) ); - end; - kDT_FLOAT: begin //note - we do not round results - all intensities might be frational... - l32f^[lPos] := - ( - {all min} ( (lXrM1*lYrM1*lZrM1)*l32fs^[lXo+lMinY+lMinZ]) - {x+1}+((lXreal*lYrM1*lZrM1)*l32fs^[lXo+1+lMinY+lMinZ]) - {y+1}+((lXrM1*lYreal*lZrM1)*l32fs^[lXo+lMaxY+lMinZ]) - {z+1}+((lXrM1*lYrM1*lZreal)*l32fs^[lXo+lMinY+lMaxZ]) - {x+1,y+1}+((lXreal*lYreal*lZrM1)*l32fs^[lXo+1+lMaxY+lMinZ]) - {x+1,z+1}+((lXreal*lYrM1*lZreal)*l32fs^[lXo+1+lMinY+lMaxZ]) - {y+1,z+1}+((lXrM1*lYreal*lZreal)*l32fs^[lXo+lMaxY+lMaxZ]) - {x+1,y+1,z+1}+((lXreal*lYreal*lZreal)*l32fs^[lXo+1+lMaxY+lMaxZ]) ); - end; - end; //case - - end; //if voxel is in source image's bounding box - end;//z - end;//y - end;//z -end else begin //if trilinear, else nearest neighbor - - for lZi := 0 to (lZ-1) do begin - //these values are the same for all voxels in the slice - // compute once per slice - lZx := lZi*lMat.matrix[1][3]; - lZy := lZi*lMat.matrix[2][3]; - lZz := lZi*lMat.matrix[3][3]; - for lYi := 0 to (lY-1) do begin - //these values change once per row - // compute once per row - lYx := lYi*lMat.matrix[1][2]; - lYy := lYi*lMat.matrix[2][2]; - lYz := lYi*lMat.matrix[3][2]; - for lXi := 0 to (lX-1) do begin - //compute each column - inc(lPos); - lXo := round(lXx^[lXi]+lYx+lZx+lMat.matrix[1][4]); - lYo := round(lXy^[lXi]+lYy+lZy+lMat.matrix[2][4]); - lZo := round(lXz^[lXi]+lYz+lZz+lMat.matrix[3][4]); - //if lZo <> 0 then - // fx(lZo); - //need to test Xreal as -0.01 truncates to zero - if (lXo >= 0) and (lYo >= 0{1}) and (lZo >= 0{1}) and - (lXo < (lXs -1)) and (lYo < (lYs -1) ) and (lZo < (lZs {-1})) - then begin - inc(lXo);//images incremented from 1 not 0 - lYo := lYo*lXs; - lZo := lZo*lXYs; - case lSrcHdr.NIFTIhdr.datatype of - kDT_UNSIGNED_CHAR : // l8is := l8is; - l8i^[lPos] :=l8is^[lXo+lYo+lZo]; - kDT_SIGNED_SHORT: l16i^[lPos] := l16is^[lXo+lYo+lZo]; - kDT_SIGNED_INT:l32i^[lPos] := l32is^[lXo+lYo+lZo]; - kDT_FLOAT: l32f^[lPos] := l32fs^[lXo+lYo+lZo]; - end; //case - - - end; //if voxel is in source image's bounding box - - end;//z - end;//y - end;//z - -end; - Freemem(lSrcBuffer); - //release lookup tables - freemem(lXx); - freemem(lXy); - freemem(lXz); - //check to see if image is empty... - lPos := 1; - case lSrcHdr.NIFTIhdr.datatype of - kDT_UNSIGNED_CHAR : while (lPos <= (lX*lY*lZ)) and (l8i^[lPos] = 0) do inc(lPos); - kDT_SIGNED_SHORT: while (lPos <= (lX*lY*lZ)) and (l16i^[lPos] = 0) do inc(lPos); - kDT_SIGNED_INT:while (lPos <= (lX*lY*lZ)) and (l32i^[lPos] = 0) do inc(lPos); - kDT_FLOAT: while (lPos <= (lX*lY*lZ)) and (l32f^[lPos] = 0) do inc(lPos); - end; //case - if lPos <= (lX*lY*lZ) then //image not empty - //Msg('Overlap') - //result := SaveNIfTICore (lDestName, lBuffAligned, kNIIImgOffset+1, lDestHdr, lPrefs,lByteSwap); - else - Msg('Overlay image does not overlap with background image.'); - //Freemem(lBuffUnaligned); - lDestHdr.ImgBufferItems := lX*lY*lZ; - case lSrcHdr.NIFTIhdr.datatype of - kDT_UNSIGNED_CHAR :lDestHdr.ImgBufferBPP :=1; - kDT_SIGNED_SHORT: lDestHdr.ImgBufferBPP :=2; - kDT_SIGNED_INT:lDestHdr.ImgBufferBPP :=4; - kDT_FLOAT: lDestHdr.ImgBufferBPP :=4; - end; //case - {if lSrcHdr.NIFTIhdr.datatype = kDT_FLOAT then begin - lInMinX := l32f^[1]; - lInMinY := l32f^[1]; - for lPos := 1 to (lX*lY*lZ) do begin - if l32f^[lPos] > lInMinX then - lInMinX :=l32f^[lPos]; - if l32f^[lPos] < lInMinY then - lInMinY :=l32f^[lPos]; - end; - fx(lInMinY,lInMinX); - end; } - NIFTIhdr_MinMaxImg(lDestHdr,lDestHdr.ImgBuffer);//set global min/max - result := 'OK'; -end; -*) - -(*function ResliceImgNIfTI (lTargetImgName,lSrcImgName,lOutputName: string): boolean; -label - 666; -var - lReslice : boolean; - lDestHdr,lSrcHdr: TMRIcroHdr; - lSrcMat,lDestMat,lSrcMatINv,lDestMatInv,lMat: TMatrix; - lOffX,lOffY,lOffZ: single; - D: double; -begin - result := false; - if not fileexists(lTargetImgName) then exit; - if not fileexists(lSrcImgName) then exit; - ImgForm.CloseImagesClick(nil); - lReslice := gBGImg.ResliceOnLoad; - gBGImg.ResliceOnLoad := false; - //if not HdrForm.OpenAndDisplayHdr(lTargetImgName,lDestHdr) then goto 666; - if not NIFTIhdr_LoadHdr(lTargetImgName, lDestHdr) then goto 666; - if not NIFTIhdr_LoadHdr(lSrcImgName, lSrcHdr) then goto 666; - if not ImgForm.OpenAndDisplayImg(lSrcImgName,false) then exit; - if not Qx(lDestHdr,lSrcHdr,lOutputName) then goto 666; - - result := true; -666: - if not result then - showmessage('Error applying transform '+lSrcImgName+'->'+lTargetImgName); - gBGImg.ResliceOnLoad := lReslice; -end; *) - end. diff --git a/prefs.pas b/prefs.pas index 1048a8c..4c22a39 100755 --- a/prefs.pas +++ b/prefs.pas @@ -14,20 +14,17 @@ interface TMRU = array [1..knMRU] of string; TPrefs = record {$IFDEF Darwin}isDoubleBuffer,{$ENDIF} - FlipYZ, - SliceDetailsCubeAndText, - FormMaximized,Debug,ColorEditor,ProportionalStretch,OverlayColorFromZero, OverlayHideZeros,SkipPrefWriting, - MaskOverlayWithBackground, InterpolateOverlays, Perspective, - FasterGradientCalculations, - ShowToolbar,ColorbarText,Colorbar,ForcePowerOfTwo, //InterpolateViewX, - RayCastShowGLSLWarnings,RayCastViewCenteredLight,EnableYoke,//Show2DSlicesDuringRendering, - //IntelWarning, dou + OverlayHideZeros,SkipPrefWriting, + FlipYZ, SliceDetailsCubeAndText,ThresholdDetection, + FormMaximized,Debug,ColorEditor,ProportionalStretch,OverlayColorFromZero, + MaskOverlayWithBackground, InterpolateOverlays, Perspective,FasterGradientCalculations, + ShowToolbar,ColorbarText,Colorbar,ForcePowerOfTwo, + RayCastShowGLSLWarnings,RayCastViewCenteredLight,EnableYoke, NoveauWarning, StartupScript : boolean; - PlanarRGB,SliceView,DrawColor,RayCastQuality1to10,FormWidth,FormHeight,//RenderQuality, - BackgroundAlpha, - OverlayAlpha,CrosshairThick,MaxVox, BitmapZoom: integer; + PlanarRGB,SliceView,DrawColor,RayCastQuality1to10,FormWidth,FormHeight, + BackgroundAlpha,OverlayAlpha,CrosshairThick,MaxVox, BitmapZoom: integer; CLUTWindowColor,CLUTIntensityColor: TColor; - GridAndBorder,BackColor,TextColor,TextBorder,CrosshairColor,HistogramColor,{HistogramGrid,}HistogramBack: TGLRGBQuad; + GridAndBorder,BackColor,TextColor,TextBorder,CrosshairColor,HistogramColor,HistogramBack: TGLRGBQuad; ColorBarPos: TUnitRect; InitScript: string; PrevFilename,PrevScriptName: TMRU; @@ -175,6 +172,7 @@ procedure SetDefaultPrefs (var lPrefs: TPrefs; lEverything: boolean); //IntelWarning := true; NoveauWarning := true; FlipYZ := false; + ThresholdDetection := true; {$IFDEF Darwin} isDoubleBuffer := false; {$ENDIF} ForcePowerOfTwo:= false; OverlayHideZeros := false; @@ -431,6 +429,7 @@ function IniFile(lRead: boolean; lFilename: string; var lPrefs: TPrefs): boolean IniBool(lRead,lIniFile, 'FormMaximized',lPrefs.FormMaximized); //IniBool(lRead,lIniFile, 'SliceText',lPrefs.SliceText); IniBool(lRead,lIniFile, 'ColorBarText',lPrefs.ColorBarText); + IniBool(lRead,lIniFile, 'ThresholdDetection',lPrefs.ThresholdDetection); IniBool(lRead,lIniFile, 'ForcePowerOfTwo',lPrefs.ForcePowerOfTwo); IniBool(lRead,lIniFile, 'StartScript',lPrefs.StartupScript); //IniBool(lRead,lIniFile, 'IntelWarning',lPrefs.IntelWarning); diff --git a/simplelaz.lpi b/simplelaz.lpi index 1045f64..653b0a7 100644 --- a/simplelaz.lpi +++ b/simplelaz.lpi @@ -16,7 +16,7 @@ - + @@ -208,7 +208,7 @@ if TargetOS='darwin' then - + @@ -225,12 +225,12 @@ if TargetOS='darwin' then - - - + + + @@ -264,7 +264,7 @@ if TargetOS='darwin' then - + @@ -272,7 +272,7 @@ if TargetOS='darwin' then - + @@ -287,26 +287,28 @@ if TargetOS='darwin' then - + - - - + + + + - + - - - + + + + - + @@ -316,17 +318,17 @@ if TargetOS='darwin' then - + - - + + - + @@ -336,19 +338,20 @@ if TargetOS='darwin' then - + - - - + + + + - + - + @@ -363,10 +366,11 @@ if TargetOS='darwin' then + - - - + + + @@ -381,16 +385,18 @@ if TargetOS='darwin' then - + - + + - - - + + + + @@ -418,18 +424,17 @@ if TargetOS='darwin' then - + - + - + - @@ -438,7 +443,7 @@ if TargetOS='darwin' then - + @@ -460,7 +465,7 @@ if TargetOS='darwin' then - + @@ -491,7 +496,7 @@ if TargetOS='darwin' then - + @@ -500,7 +505,7 @@ if TargetOS='darwin' then - + @@ -521,14 +526,14 @@ if TargetOS='darwin' then - + - + @@ -546,7 +551,7 @@ if TargetOS='darwin' then - + @@ -555,7 +560,7 @@ if TargetOS='darwin' then - + @@ -564,7 +569,7 @@ if TargetOS='darwin' then - + @@ -573,7 +578,7 @@ if TargetOS='darwin' then - + @@ -582,354 +587,335 @@ if TargetOS='darwin' then - + - - - - - - - - - - - - - - - - - - - - + + + - - - + + + - - - + + + - - - + + + - - - + + + - - - + + + - + - - + + - + - - + + - + - - + + - + - - + + - - - - - + + + + + - - - - - - - - + - - + + - + - - - + + + - - + + - + - - + + - + - - + + - + - - + + - + - - + + - + - - + + - + - - + + - + - - + + - - - + + + - + - + + + + + + + + + + - + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - - - diff --git a/texture2raycast.pas b/texture2raycast.pas index a3cab14..bd8c462 100755 --- a/texture2raycast.pas +++ b/texture2raycast.pas @@ -29,41 +29,36 @@ implementation const kRGBAclear : TRGBA = (r: 0; g: 0; b: 0; a:0); - -(* -USE XYZIx instead it is FASTER and BETTER (preserves more percision in 8-bit float) -Function XYZI (X1,X2,Y1,Y2,Z1,Z2: single): TRGBA; -//Function XYZI (X1,X2,Y1,Y2,Z1,Z2: single; Center: byte): TRGBA; +{$DEFINE GRADIENT_PRENORM} +//pre-normalizing our data allows us to avoids the "normalize" in GLSL +// gradientSample.rgb = normalize(gradientSample.rgb*2.0 - 1.0); +//it is a bit slower and does provide a bit less precision +//in 2017 we switched to the pre-norm so that the CPU and GLSL gradients match each other +{$IFDEF GRADIENT_PRENORM} +Function XYZIx (X1,X2,Y1,Y2,Z1,Z2: single): TRGBA; //input voxel intensity to the left,right,anterior,posterior,inferior,superior and center // Output RGBA image where values correspond to X,Y,Z gradients and ImageIntensity -//AlphaT will make a voxel invisible if center intensity is less than specified value -// Voxels where there is no gradient (no edge boundary) are made transparent var X,Y,Z,Dx: single; begin Result := kRGBAclear; - //if Center < 1 then - // exit; //intensity less than threshold: make invisible X := X1-X2; Y := Y1-Y2; Z := Z1-Z2; Dx := sqrt(X*X+Y*Y+Z*Z); if Dx = 0 then - exit; //no gradient - set intensity to zero. Calculate_Transfer_Function - result.r :=round((X/(Dx*2)+0.5)*255); //X - result.g :=round((Y/(Dx*2)+0.5)*255); //Y - result.B := round((Z/(Dx*2)+0.5)*255); //Z - //result.A := Center; -end;*) - + exit; //no gradient - set intensity to zero. + result.r :=round((X/(Dx*2)+0.5)*255); + result.g :=round((Y/(Dx*2)+0.5)*255); + result.B := round((Z/(Dx*2)+0.5)*255); +end; +{$ELSE GRADIENT_PRENORM} Function XYZIx (X1,X2,Y1,Y2,Z1,Z2: single): TRGBA; {$IFDEF FPC} inline; {$ENDIF} -//Function XYZIx (X1,X2,Y1,Y2,Z1,Z2: single; Center: byte): TRGBA; -//faster version of XYZI +//faster, more precise version of XYZI for computation, but requires "normalize" for display var X,Y,Z,Dx: single; begin Result := kRGBAclear; - //if Center < 1 then exit; //intensity less than threshold: make invisible X := X1-X2; Y := Y1-Y2; Z := Z1-Z2; @@ -77,6 +72,7 @@ implementation result.g :=round((Y/Dx+0.5)*255); //Y result.B := round((Z/Dx+0.5)*255); //Z end; +{$ENDIF GRADIENT_PRENORM} {$DEFINE SOBEL} {$IFDEF SOBEL} //use SOBEL @@ -88,10 +84,6 @@ function estimateGradients (rawData: tVolB; Xsz,Ysz, I : integer; var GradMag: s Y,Z,J: integer; Xp,Xm,Yp,Ym,Zp,Zm: single; begin - //GradMag := 0;//gradient magnitude - //Result := kRGBAclear; - //if rawData[i] < 1 then - // exit; //intensity less than threshold: make invisible Y := XSz; //each row is X voxels Z := YSz*XSz; //each plane is X*Y voxels //X:: cols: +Z +0 -Z, rows -Y +0 +Y @@ -257,7 +249,7 @@ procedure CreateGradientVolumeCPU (var VolRGBA: tVolRGBA; dim1, dim2, dim3: inte Index := (Z * dim1x2) + (Y * dim1) + X; //estimate gradients using Sobel or Central Difference (depending on DEFINE SOBEL) if (VolData[Index] <> 0) then - VolRGBA[Index] := estimateGradients (VolData, dim1,dim2, Index,GradMagS[Index]); + VolRGBA[Index] := estimateGradients (VolData, dim1,dim2, Index,GradMagS[Index]); end;//X VolData := nil;//FREE ---- //next: generate normalized gradient magnitude values diff --git a/texture_3d_unit.pas b/texture_3d_unit.pas index 35c0dd7..2e04434 100755 --- a/texture_3d_unit.pas +++ b/texture_3d_unit.pas @@ -1110,9 +1110,7 @@ procedure SetOriginXYZ (var lTexture3D: TTexture); // Here I take the latter approach result :=false; lFilename := F_Filename; - //deleteGradients (lTexture); InitTexture(lTexture); - if lFilename = '' then begin if not NIFTIhdr_LoadDummyImg (lHdr, lImgBuffer) then exit;