{
This group of macros will do several things that I have found useful. 
I took some things directly from the stacks macro that is included with
NIH-Image. The most useful things that I have added are more robust macros for
reading GS Signa MR files that have been extracted with wither ximg or the GINX 
tool.

Since I don't do much CT I have not done anything to read the CT headers.

Also is a function called 'Fermi' that very closely mimics the window and
leveling of images on the GE MR Signa scanner.

 	Note:  Make sure that the image that you are windowing and
	leveling is the active window or else wierd things happen!

If you have any questions or problems please email phil@peace.med.ohio-state.edu

I make no promises for these macros and post them as-is.

Phil Williams

}
Var 
  InfoPid : integer;
  gprefix : string;
  fov,thick : real;
  tr,te : real;
  Imin,Imax : integer;
  examname,infotitle : string;

{--------Some check procedures from the stacks macro-----}
procedure CheckForStack;
begin
  if nPics=0 then begin
    PutMessage('This macro requires a stack.');
    exit;
  end;
  if nSlices=0 then begin
    PutMessage('This window is not a stack.');
    exit
  end;
end;

procedure CheckForSelection;
var 
  x1,y1,x2,y2,LineWidth:integer;
begin
  GetRoi(RoiLeft,RoiTop,RoiWidth,RoiHeight);
  GetLine(x1,y1,x2,y2,LineWidth);
  if (RoiWidth=0) or (x1>=0) then begin
    PutMessage('Please make a rectangular selection.');
    exit;
  end;
end;
{------------window and leveling macro-------------------}
Procedure fermiscl(window,level : integer);
var
  i,temp : integer;
begin
  Copy;
  for i := 1 to 256 do begin
    temp := ((255-(i-1))-level)/window;
    RedLUT[i-1] := 255*(1-1/(exp(temp) + 1));
    BlueLUT[i-1] := 255*(1-1/(exp(temp) + 1));
    GreenLUT[i-1] := 255*(1-1/(exp(temp) + 1));
  end; {for}
  Paste;
  UpdateLUT;
end;

{----------------utilties for reading in header and file info for Signa files------}

procedure thirtytwo(last,row,HdrPid : integer);
{
	converts 4 bytes to a 32 bit int for reading in
	info from the header
}
var
  first,sec,third,fourth : integer;

begin
  ChoosePic(HdrPid);
  first := GetPixel(last,row);
  sec := GetPixel(last-1,row);
  third := GetPixel(last-2,row);
  fourth := GetPixel(last-3,row);
  result := first+sec*256+third*256*256+fourth*256*256*256;
end;

procedure pow(base : real;exponent : integer);
{
	raise a number to a power
}
var
 i,p : integer;

begin
  p := 1;
  for i := 1 to exponent do
     p := p*base;
  result2 := p;
end;

procedure btf(loc,Pid : integer);
{
	converts 4 bytes to a 32 bit float for reading in 
	info from the image header
}
var
  first,second,third,fourth,sign,power : integer;

  bit0,bit1,bit2,bit3,bit4,bit5,bit6,bit7,
  bit8,bit9,bit10,bit11,bit12,bit13,bit14,bit15,
  bit16,bit17,bit18,bit19,bit20,bit21,bit22,bit23,
  bit24,bit25,bit26,bit27,bit28,bit29,bit30,bit31 : integer;

  fraction, mantissa,result2 : real;

begin
  	ChoosePic(Pid);
   	first := GetPixel(loc,0);
   	second := GetPixel(loc+1,0);
   	third := GetPixel(loc+2,0);
   	fourth := GetPixel(loc+3,0);

{ Get whether the bits are set }

   	bit0 := BitAnd(128,first); bit1 := BitAnd(64,first); 
   	bit2 := BitAnd(32,first); bit3 := BitAnd(16,first);
   	bit4 := BitAnd(8,first); bit5 := BitAnd(4,first); 
   	bit6 := BitAnd(2,first); bit7 := BitAnd(1,first);

   	bit8 := BitAnd(128,second); bit9 := BitAnd(64,second);
   	bit10 := BitAnd(32,second); bit11 := BitAnd(16,second);
   	bit12 := BitAnd(8,second); bit13 := BitAnd(4,second);
   	bit14 := BitAnd(2,second); bit15 := BitAnd(1,second);

   	bit16 := BitAnd(128,third); bit17 := BitAnd(64,third);
   	bit18 := BitAnd(32,third); bit19 := BitAnd(16,third);
   	bit20 := BitAnd(8,third); bit21 := BitAnd(4,third);
   	bit22 := BitAnd(2,third); bit23 := BitAnd(1,third);

  	bit24 := BitAnd(128,fourth); bit25 := BitAnd(64,fourth);
  	bit26 := BitAnd(32,fourth); bit27 := BitAnd(16,fourth); 
  	bit28 := BitAnd(8,fourth); bit29 := BitAnd(4,fourth); 
  	bit30 := BitAnd(2,fourth); bit31 := BitAnd(1,fourth);

{ Now build the sign}
   	if bit0 > 0 then sign := -1 else sign := 1;
{Now build the exponent}
  	power := 0;
  	if bit1 > 0 then power := power + 128;
  	if bit2 > 0 then power := power + 64;
  	if bit3 > 0 then power := power + 32;
  	if bit4 > 0 then power := power + 16;
  	if bit5 > 0 then power := power + 8;
  	if bit6 > 0 then power := power + 4;
  	if bit7 > 0 then power := power + 2;
  	if bit8 > 0 then power := power + 1;
 	power := power - 127;
 
 {now time for the fraction }
 	fraction := 0;
 	if bit9 > 0 then fraction := fraction + 1/2;
  	if bit10 > 0 then fraction := fraction + 1/4;
  	if bit11 > 0 then fraction := fraction + 1/8;
  	if bit12 > 0 then fraction := fraction + 1/16;
  	if bit13 > 0 then fraction := fraction + 1/32;
  	if bit14 > 0 then fraction := fraction + 1/64;
  	if bit15 > 0 then fraction := fraction + 1/128;
  	if bit16 > 0 then fraction := fraction + 1/256;
  	if bit17 > 0 then fraction := fraction + 1/512;
  	if bit18 > 0 then fraction := fraction + 1/1024;
  	if bit19 > 0 then fraction := fraction + 1/2048;
  	if bit20 > 0 then fraction := fraction + 1/4096;
  	if bit21 > 0 then fraction := fraction + 1/8192;
  	if bit22 > 0 then fraction := fraction + 1/16384;
  	if bit23 > 0 then fraction := fraction + 1/32768;
  	if bit24 > 0 then fraction := fraction + 1/65536;
  	if bit25 > 0 then fraction := fraction + 1/131072;
  	if bit26 > 0 then fraction := fraction + 1/262144;
  	if bit27 > 0 then fraction := fraction + 1/524288;
  	if bit28 > 0 then fraction := fraction + 1/1048576;
  	if bit29 > 0 then fraction := fraction + 1/2097152;
  	if bit30 > 0 then fraction := fraction + 1/4194304;
  	if bit31 > 0 then fraction := fraction + 1/8388608;
 
	mantissa := 1 + fraction;
	pow(2,power);
	result := sign*result2*mantissa;
end;

procedure ReadSignaHeader(headername: string);
{
	Reads in the header info of a Signa formatted file
}
Var
  HdrPid : integer;
  result : real;

begin
	SetImport('Custom; 8-bits;');
	SetCustom(156,1,0);
	Import(headername);
	HdrPid := PidNumber;
	ChoosePic(HdrPid);
	thirtytwo(7,0,HdrPid);
	offset := result;
	thirtytwo(11,0,HdrPid);
	width := result;
	thirtytwo(15,0,HdrPid);
	height := result;
	thirtytwo(83,0,HdrPid);
	HHdroff := result;
	thirtytwo(87,0,HdrPid);
	HHdrlen := result;
	thirtytwo(135,0,HdrPid);
	EHdroff := result;
	thirtytwo(139,0,HdrPid);
	EHdrlen := result;
	thirtytwo(151,0,HdrPid);
	IHdroff := result;
	thirtytwo(155,0,HdrPid);
	IHdrLen := result;
	SelectPic(HdrPid);
	Dispose;
end;

Procedure ReadExamHeader(name,offset,len);
{
	Read the exam header
}
Var
  EPid,j : integer;
  str : string;
  temp : integer;

begin
	SetImport('Custom; 8-bits;');
	SetCustom(len,1,offset);
	Import(name);
	EPid := PidNumber;
	str := '';
	for j := 0 to 24 do
    	str := concat(str,chr(GetPixel(97+j,0)));
	examname := str;
	SelectPic(EPid);
	Dispose;
end;


Procedure ReadImageHeader(name,offset,len);
{
	Read the MR image header
}
Var
   IPid,loc,j : integer;
   str : string;
   result : real;

begin
	SetImport('Custom;8-bits;');
	SetCustom(len,1,offset);
	Import(name);
	IPid := PidNumber;
	loc := 0;
	btf(26,IPid);
	thick := result;
	btf(34,IPid);
	fov := result;
	thirtytwo(197,0,IPid);
	tr := result/1000;
	thirtytwo(205,0,IPid);
	te := result/1000;
	SelectPic(IPid);
	Dispose;
end;

Procedure ReadHistogramHeader(name,offset,len);
{
	Read the histrogram header to get initial window and level
	values
	
	not used at present
}
Var
   IPid,loc,j : integer;

begin
	SaveState;
	SetImport('Custom;8-bit;');
	SetCustom(len,1,offset);
	Import(name);
	IPid := PidNumber;
	loc := 0;
	min := GetPixel(10,0)*256+GetPixel(11,0);
	max := GetPixel(12,0)*256+GetPixel(13,0);
	SelectPic(IPid);
	Dispose;
	RestoreState;
end;

Procedure ReadSignaFiles(ginx : boolean);
{
	read in the files.  the boolean ginx only affects how the files are named
	since files extracted with the GINX tool and those with ximg have different
	naming conventions
}
Var
	filename,prefix,ext : string;
	i,width,height,offset,first,last,min,max,
	EHdrOff,EHdrLen,
	IHdrOff,IHdrLen,
	HHdrOff,HHdrLen,
	stack,n,
	skip : integer;
Begin
	SaveState;
	if ginx then
	    prefix := 'I.'
	else begin
		prefix := GetString('File Prefix:',gprefix);
		gprefix := prefix;
	end;
  	ext := '.MR';
	first := round(GetNumber('First file number',1));
	last := -1;
	repeat 
		last := round(GetNumber('Last file number',60));
		if (last < first) then
			PutMessage('Last filenumber must be greater than or equal to the first');
	until (last >= first);
	skip := round(GetNumber('How many files to skip:',1));
	if ginx then
		filename := concat('I.',first:3)
	else
		filename := concat(prefix,first:1,ext);
	ReadSignaHeader(filename);
	ReadExamHeader(filename,EHdrOff,EHdrLen);
	ReadImageHeader(filename,IHdrOff,IHdrLen);
	SetNewSize(width,height);
	MakeNewStack(examname);
	stack := nPics;
	n := first;
	SetImport('Custom; 16-bits Unsigned; Fixed Scale');
	SetImportMinMax(0,255);
	SetCustom(width,height,offset,1);
	i := first;
	while i <= last do begin
		if ginx then begin
			filename := concat('I.',i:3);
			ReadHistogramHeader(filename,HHdrOff,HHdrLen);
			SetImport('Custom; 16-bits Unsigned; Fixed Scale');
			SetImportMinMax(0,255);
			SetCustom(width,height,offset,1);
			Import(filename);
			SetPicName(filename);
   		end
		else begin
			if (i>=0) and (i<10) then begin
				filename := concat(prefix,i:1,ext);
				ReadHistogramHeader(filename,HHdrOff,HHdrLen);
				SetImport('Custom; 16-bits Unsigned; Fixed Scale');
				SetImportMinMax(0,255);
				SetCustom(width,height,offset,1);
				Import(filename);
				SetPicName(filename);
			end
			else begin 
				if (i>=10) and (i<100) then begin
					filename := concat(prefix,i:2,ext);
					ReadHistogramHeader(filename,HHdrOff,HHdrLen);
					SetImport('Custom; 16-bits Unsigned; Fixed Scale');
					SetImportMinMax(0,255);
					SetCustom(width,height,offset,1);
					Import(filename);
					SetPicName(filename);
				end
				else begin
					filename := concat(prefix,i:3,ext);
					ReadHistogramHeader(filename,HHdrOff,HHdrLen);
					SetImport('Custom; 16-bits Unsigned; Fixed Scale');
					SetImportMinMax(0,255);
					SetCustom(width,height,offset,1);
					Import(filename);
					SetPicName(filename);
				end;
			end; {else2}
		end; {if ginx}
		SelectAll;
		Copy;
		Dispose;
		SelectPic(stack);
		if n = first then begin Imax := max; Imin := min; end;
		if n<>first then  AddSlice;
		if max > Imax then Imax := max;
		if min < Imin then Imin := min;
		n:=n+1;
		Paste;
		i := i+(skip+1);
	end; {while}
	RestoreState;
	KillROI;
	SetScale(width/fov,'mm');
	infotitle := concat(examname,' Info');
	NewTextWindow(infotitle,200,300);
	InfoPid := PidNumber;
	MoveWindow(width+90,40);
	SelectWindow(infotitle);
	writeln('FOV = ',fov:3,' mm');
	writeln('Slice Thickness = ',thick:2:1' mm');
	writeln('TR = ',tr,' ms'); writeln('TE = ',te,' ms');
	writeln('min = ',Imin);writeln('max = ',Imax);
end;

macro 'Read Ximg [X]' begin ReadSignaFiles(false); end;
macro 'Read GINX [G]' begin ReadSignaFiles(true); end;

macro '(-' begin end;

{--------reformatting macros taken from stacks-----------------}
procedure ResliceSignaMRI(horizontal,OptionKey:boolean);
{
	Changed this to use the scale that was set from ReadSignaFiles
}
var
  stack1,stack2,width,height:integer;
  RoiLeft,RoiTop,RoiWidth,RoiHeight,max:integer;
  loc,PixelSpacing:real;
  InputSpacing,OutputSpacing:real; {mm}
  scale:real; {pixels/mm}  
  FirstTime:boolean;
  units : string;
begin
  RequiresVersion(1.45);
  units := 'mm';
  GetScale(scale,units);
  CheckForStack;
  CheckForSelection;
  SaveState;
  SetScale(scale,'mm');
  SetBackground(0);
  SetBackground(255);
  stack1:=PicNumber;
  InputSpacing:=GetSliceSpacing/scale;
  if InputSpacing<=0 then InputSpacing:=1.5;
  InputSpacing:=GetNumber('Input Slice Spacing(mm):',InputSpacing);
  SetSliceSpacing(InputSpacing*scale);
  OutputSpacing:=InputSpacing;
  OutputSpacing:=GetNumber('Output Slice Spacing (mm):', OutputSpacing);
  PixelSpacing:=OutputSpacing*scale;
  FirstTime:=true;
  GetRoi(RoiLeft,RoiTop,RoiWidth,RoiHeight);
  if horizontal then begin
    loc:=RoiTop+PixelSpacing;
    max:=RoiTop+RoiHeight;
  end else begin
    loc:=RoiLeft+PixelSpacing;
    max:=RoiLeft+RoiWidth;
  end;
  while loc<max do begin
    ChoosePic(stack1);
    if horizontal
      then MakeLineRoi(RoiLeft,loc,RoiLeft+RoiWidth,loc)
      else MakeLineRoi(loc,RoiTop,loc,RoiTop+RoiTop+RoiHeight);
    if OptionKey then SetOption;
    Reslice;
    SelectAll;
    Copy;
    GetPicSize(width,height);
    Dispose;
    if FirstTime then begin
      SetNewSize(width,height);
      MakeNewStack(OutputSpacing:1:2);
      SetSliceSpacing(PixelSpacing);
      stack2:=PicNumber;
    end;
    ChoosePic(stack2);
    if not FirstTime then AddSlice;
    Paste;
    loc:=loc+PixelSpacing;
    FirstTime:=false;
  end;
  SelectPic(stack1);
  KillRoi;
  SelectPic(stack2);
  KillRoi;
  RestoreState;
end;

macro 'Sagitals to CoronalsÉ'; begin ResliceSignaMRI(false,true) end;
macro 'Sagitals to AxialsÉ'; begin ResliceSignaMRI(true,true) end;
macro 'Coronals to SagitalsÉ'; begin ResliceSignaMRI(false,true) end;
macro 'Axials to Coronals. . .';begin ResliceSignaMRI(true,false) end;
macro 'Axials to Sagitals. . .';begin ResliceSignaMRI(false,false) end;

macro '(-' begin end;

{------------LUT macros---------------------}
macro 'Fermi [F]'
{
	window and level an image like the scanners do!
	
 	Note:  Make sure that the image that you are windowing and
	leveling is the active window or else wierd things happen!
}
var
 x,y,width,height : integer;
 window,level :integer;

begin
ShowMessage('Press mouse button to exit.\');
GetPicSize(width,height);
repeat
  SelectAll;
  GetMouse(x,y);
  window := x; level := y;
  if window <= 0 then window := 1;
  fermiscl(window,level);
  ShowMessage('w = ',window:3,'  l = ',level:3);
until button;
ShowMessage('');
end;

{---------these are from the LUT macros----------------}
macro 'Export LUT [E]';
{Copies the current look-up table to a text window.}
var
  i:integer;
  v:real;
  tab:string;
begin
  RequiresVersion(1.54);
  NewTextWindow('LUT',200,400);
  tab:=chr(9);
  for i:=0 to 255 do
    Writeln(i:4,tab,RedLut[i]:4,tab,GreenLut[i]:4,tab,BlueLut[i]:4);
end;

macro 'Plot LUT [P]';
var
  i,xscale,yscale:real;
  width,height,margin,pwidth,pheight:integer;
  xbase,ybase:integer;
begin
  SaveState;
  margin:=25;
  pwidth:=400;
  pheight:=125;
  width:=pwidth+2*margin;
  height:=pheight*3+2*margin;
  SetNewSize(width,height);
  SetBackground(0); 
  MakeNewWindow('LUT');
  xscale:=(pwidth-2)/256;
  yscale:=(pheight-1)/256;
  SetForeground(252);
  xbase:=margin; ybase:=margin;
  MoveTo(xbase,ybase);
  for i:=0 to 255 do
    LineTo(xbase+i*xscale,ybase+RedLUT[i]*yscale);
  SetForeground(255);
  MakeRoi(xbase,ybase,pwidth,pheight);
  FlipVertical;
  DrawBoundary;
  SetForeground(253);
  ybase:=ybase+pheight-1;
  MoveTo(xbase,ybase);
  for i:=0 to 255 do
    LineTo(xbase+i*xscale,ybase+GreenLUT[i]*yscale);
  SetForeground(255);
  MakeRoi(xbase,ybase,pwidth,pheight);
  FlipVertical;
  DrawBoundary;
  SetForeground(254);
  ybase:=ybase+pheight-1;
  MoveTo(xbase,ybase);
  for i:=0 to 255 do
    LineTo(xbase+i*xscale,ybase+BlueLUT[i]*yscale);
  SetForeground(255);
  MakeRoi(xbase,ybase,pwidth,pheight);
  FlipVertical;
  DrawBoundary;
  KillRoi;
  RedLUT[252]:=255; GreenLUT[252]:=0;   BlueLUT[252]:=0;
  RedLUT[253]:=0;   GreenLUT[253]:=255; BlueLUT[253]:=0;
  RedLUT[254]:=0;   GreenLUT[254]:=0;   BlueLUT[254]:=255;
  UpdateLUT;
  SetFont('Geneva');
  SetFontSize(9);
  SetText('Centered');
  MoveTo(margin+4,height-margin+8);
  writeln(0:1:2);
  MoveTo(margin+pwidth,height-margin+8);
  writeln(255:1:2);
  RestoreState;
end;

macro 'Reset LUT [R]';
var 
  i : integer;

begin
 ResetGrayMap;
end;

macro '(-' begin end;

{----------------this is from the stacks macro------------}
procedure CropAndScale(fast:boolean; angle:real);
var
  i,OldStack,NewStack:integer;
  RoiLeft,RoiTop,RoiWidth,RoiHeight:integer;
  N,NewWidth:integer;
  ScaleFactor:real;
  OneToOne:boolean;
begin
  CheckForStack;
  CheckForSelection;
  SaveState;
  OldStack:=PicNumber;
  N:=nSlices;
  ScaleFactor:=GetNumber('Scale factor(0.05..25):',1.0);
  OneToOne:=ScaleFactor=1.0;
  NewWidth:=round(RoiWidth*ScaleFactor);
  if odd(NewWidth) then begin
    NewWidth:=NewWidth-1;
    ScaleFactor:=NewWidth/RoiWidth;
  end;
  SetNewSize(RoiWidth*ScaleFactor,RoiHeight*ScaleFactor);
  MakeNewStack('Stack');
  NewStack:=PicNumber;
  if not OneToOne then begin
    if fast 
      then SetScaling('Nearest; Create New Window')
      else SetScaling('Bilinear; Create New Window');
  end;
  SelectPic(OldStack);
  for i:= 1 to N do begin
    SelectSlice(1);
    if OneToOne and (angle=0.0) then Duplicate('Temp')
      else ScaleAndRotate(ScaleFactor,ScaleFactor,angle);
    SelectAll;
    Copy;
    SelectPic(NewStack);
    if i<>1 then AddSlice;
    Paste;
    SelectPic(nPics);
    Dispose; {Temp}
    SelectPic(OldStack);
    DeleteSlice;
  end;
  Dispose; {OldStack}
  RestoreState;
end;

macro 'Crop and Scale-FastÉ';   begin CropAndScale(true, 0); end;
macro 'Crop and Scale-SmoothÉ'; begin CropAndScale(false, 0); end;