diff --git a/.DS_Store b/.DS_Store
index cbdaae5..ec4b42d 100644
Binary files a/.DS_Store and b/.DS_Store differ
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;