Rewrite the Readme.md

NIfTI_20140122/bipolar.m

@@ -0,0 +1,94 @@
+%BIPOLAR returns an M-by-3 matrix containing a blue-red colormap, in
+%	in which red stands for positive, blue stands for negative, 
+%	and white stands for 0.
+%
+%  Usage: cmap = bipolar(M, lo, hi, contrast);  or  cmap = bipolar;
+%
+%  cmap:  output M-by-3 matrix for BIPOLAR colormap.
+%  M:	  number of shades in the colormap. By default, it is the
+%	  same length as the current colormap.
+%  lo:	  the lowest value to represent.
+%  hi:	  the highest value to represent.
+%
+%  Inspired from the LORETA PASCAL program:
+%	http://www.unizh.ch/keyinst/NewLORETA
+%
+%  [email protected]
+%
+%----------------------------------------------------------------
+function cmap = bipolar(M, lo, hi, contrast)
+
+   if ~exist('contrast','var')
+      contrast = 128;
+   end
+
+   if ~exist('lo','var')
+      lo = -1;
+   end
+
+   if ~exist('hi','var')
+      hi = 1;
+   end
+
+   if ~exist('M','var')
+      cmap = colormap;
+      M = size(cmap,1);
+   end
+
+   steepness = 10 ^ (1 - (contrast-1)/127);
+   pos_infs = 1e-99;
+   neg_infs = -1e-99;
+
+   doubleredc = [];
+   doublebluec = [];
+
+   if lo >= 0		% all positive
+
+      if lo == 0
+         lo = pos_infs;
+      end
+
+      for i=linspace(hi/M, hi, M)
+         t = exp(log(i/hi)*steepness);
+         doubleredc = [doubleredc; [(1-t)+t,(1-t)+0,(1-t)+0]];
+      end
+
+      cmap = doubleredc;
+
+   elseif hi <= 0	% all negative
+
+      if hi == 0
+         hi = neg_infs;
+      end
+
+      for i=linspace(abs(lo)/M, abs(lo), M)
+         t = exp(log(i/abs(lo))*steepness);
+         doublebluec = [doublebluec; [(1-t)+0,(1-t)+0,(1-t)+t]];
+      end
+
+      cmap = flipud(doublebluec);
+
+   else
+
+      if hi > abs(lo)
+         maxc = hi;
+      else
+         maxc = abs(lo);
+      end
+
+      for i=linspace(maxc/M, hi, round(M*hi/(hi-lo)))
+         t = exp(log(i/maxc)*steepness);
+         doubleredc = [doubleredc; [(1-t)+t,(1-t)+0,(1-t)+0]];
+      end      
+
+      for i=linspace(maxc/M, abs(lo), round(M*abs(lo)/(hi-lo)))
+         t = exp(log(i/maxc)*steepness);
+         doublebluec = [doublebluec; [(1-t)+0,(1-t)+0,(1-t)+t]];
+      end
+
+      cmap = [flipud(doublebluec); doubleredc];
+
+   end
+
+   return;					% bipolar
+

NIfTI_20140122/bresenham_line3d.m

@@ -0,0 +1,189 @@
+%  Generate X Y Z coordinates of a 3D Bresenham's line between
+%  two given points.
+%
+%  A very useful application of this algorithm can be found in the
+%  implementation of Fischer's Bresenham interpolation method in my
+%  another program that can rotate three dimensional image volume
+%  with an affine matrix:
+%  http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=21080
+%
+%  Usage: [X Y Z] = bresenham_line3d(P1, P2, [precision]);
+%
+%  P1	- vector for Point1, where P1 = [x1 y1 z1]
+%
+%  P2	- vector for Point2, where P2 = [x2 y2 z2]
+%
+%  precision (optional) - Although according to Bresenham's line
+%	algorithm, point coordinates x1 y1 z1 and x2 y2 z2 should
+%	be integer numbers, this program extends its limit to all
+%	real numbers. If any of them are floating numbers, you
+%	should specify how many digits of decimal that you would
+%	like to preserve. Be aware that the length of output X Y
+%	Z coordinates will increase in 10 times for each decimal
+%	digit that you want to preserve. By default, the precision
+%	is 0, which means that they will be rounded to the nearest
+%	integer.
+%
+%  X	- a set of x coordinates on Bresenham's line
+%
+%  Y	- a set of y coordinates on Bresenham's line
+%
+%  Z	- a set of z coordinates on Bresenham's line
+%
+%  Therefore, all points in XYZ set (i.e. P(i) = [X(i) Y(i) Z(i)])
+%  will constitute the Bresenham's line between P1 and P1.
+%
+%  Example:
+%	P1 = [12 37 6];     P2 = [46 3 35];
+%	[X Y Z] = bresenham_line3d(P1, P2);
+%	figure; plot3(X,Y,Z,'s','markerface','b');
+%
+%  This program is ported to MATLAB from:
+%
+%  B.Pendleton.  line3d - 3D Bresenham's (a 3D line drawing algorithm)
+%  ftp://ftp.isc.org/pub/usenet/comp.sources.unix/volume26/line3d, 1992
+%
+%  Which is also referenced by:
+%
+%  Fischer, J., A. del Rio (2004).  A Fast Method for Applying Rigid
+%  Transformations to Volume Data, WSCG2004 Conference.
+%  http://wscg.zcu.cz/wscg2004/Papers_2004_Short/M19.pdf
+%
+%  - Jimmy Shen ([email protected])
+%
+function [X,Y,Z] = bresenham_line3d(P1, P2, precision)
+
+   if ~exist('precision','var') | isempty(precision) | round(precision) == 0
+      precision = 0;
+      P1 = round(P1);
+      P2 = round(P2);
+   else
+      precision = round(precision);
+      P1 = round(P1*(10^precision));
+      P2 = round(P2*(10^precision));
+   end
+
+   d = max(abs(P2-P1)+1);
+   X = zeros(1, d);
+   Y = zeros(1, d);
+   Z = zeros(1, d);
+
+   x1 = P1(1);
+   y1 = P1(2);
+   z1 = P1(3);
+
+   x2 = P2(1);
+   y2 = P2(2);
+   z2 = P2(3);
+
+   dx = x2 - x1;
+   dy = y2 - y1;
+   dz = z2 - z1;
+
+   ax = abs(dx)*2;
+   ay = abs(dy)*2;
+   az = abs(dz)*2;
+
+   sx = sign(dx);
+   sy = sign(dy);
+   sz = sign(dz);
+
+   x = x1;
+   y = y1;
+   z = z1;
+   idx = 1;
+
+   if(ax>=max(ay,az))			% x dominant
+      yd = ay - ax/2;
+      zd = az - ax/2;
+
+      while(1)
+         X(idx) = x;
+         Y(idx) = y;
+         Z(idx) = z;
+         idx = idx + 1;
+
+         if(x == x2)		% end
+            break;
+         end
+
+         if(yd >= 0)		% move along y
+            y = y + sy;
+            yd = yd - ax;
+         end
+
+         if(zd >= 0)		% move along z
+            z = z + sz;
+            zd = zd - ax;
+         end
+
+         x  = x  + sx;		% move along x
+         yd = yd + ay;
+         zd = zd + az;
+      end
+   elseif(ay>=max(ax,az))		% y dominant
+      xd = ax - ay/2;
+      zd = az - ay/2;
+
+      while(1)
+         X(idx) = x;
+         Y(idx) = y;
+         Z(idx) = z;
+         idx = idx + 1;
+
+         if(y == y2)		% end
+            break;
+         end
+
+         if(xd >= 0)		% move along x
+            x = x + sx;
+            xd = xd - ay;
+         end
+
+         if(zd >= 0)		% move along z
+            z = z + sz;
+            zd = zd - ay;
+         end
+
+         y  = y  + sy;		% move along y
+         xd = xd + ax;
+         zd = zd + az;
+      end
+   elseif(az>=max(ax,ay))		% z dominant
+      xd = ax - az/2;
+      yd = ay - az/2;
+
+      while(1)
+         X(idx) = x;
+         Y(idx) = y;
+         Z(idx) = z;
+         idx = idx + 1;
+
+         if(z == z2)		% end
+            break;
+         end
+
+         if(xd >= 0)		% move along x
+            x = x + sx;
+            xd = xd - az;
+         end
+
+         if(yd >= 0)		% move along y
+            y = y + sy;
+            yd = yd - az;
+         end
+
+         z  = z  + sz;		% move along z
+         xd = xd + ax;
+         yd = yd + ay;
+      end
+   end
+
+   if precision ~= 0
+      X = X/(10^precision);
+      Y = Y/(10^precision);
+      Z = Z/(10^precision);
+   end
+
+   return;					% bresenham_line3d
+