function readmriscan(action)
%   This software is protected by german copyright and international treaties.             
%   Copyright 2004 Markus Junghöfer & Peter Peyk. All Rights Reserved.                     
%                                                                                          
%   THIS SOFTWARE AND DOCUMENTATION IS PROVIDED "AS IS," AND COPYRIGHT HOLDERS MAKE        
%   NO REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO,    
%   WARRANTIES OF MERCHANTABILITY OR FITNESS FOR ANY PARTICULAR PURPOSE OR THAT THE USE    
%   OF THE SOFTWARE OR DOCUMENTATION WILL NOT INFRINGE ANY THIRD PARTY PATENTS, COPYRIGHTS,
%   TRADEMARKS OR OTHER RIGHTS. COPYRIGHT HOLDERS WILL NOT BE LIABLE FOR ANY DIRECT,       
%   INDIRECT, SPECIAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF ANY USE OF THE SOFTWARE      
%   OR DOCUMENTATION.                                                                      
if nargin<1, action =[];end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

readmriscanfig = findobj('tag','readmriscanfig');
readmriscanfiglist = get(readmriscanfig,'userdata');





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if strcmp(action,'chooseimgfile')


    [FilePathImg,DefFolder]= uigetfile('*.*','Please choose img-file:');
    namebak = FilePathImg;
    cd(DefFolder);
    Mask=['*.img'];
    FilePathImg = [DefFolder,FilePathImg];
    FilePathExt=lower(FilePathImg(length(FilePathImg)-3:length(FilePathImg)));
    if strcmp(FilePathExt,'.hdr')
        FilePathImg=[FilePathImg(1:length(FilePathImg)-4),'.img'];
        FilePathExt='.img';
    end
    if ~strcmp(FilePathExt,'.img')
        error('Please choose an analyze .img format file.')
    end
    FilePathHdr=[FilePathImg(1:length(FilePathImg)-4),'.hdr'];
	[Dim,Vox,Scale,Type,Offset,Origin,Descrip]=spm_hread(FilePathHdr);
	ProdDim=prod(Dim);
	if Type==2
        ImgDataFormat='ubit8';
	elseif Type==3
        ImgDataFormat='ubit16';
	elseif Type==4
        ImgDataFormat='int16';
	elseif Type==16
        ImgDataFormat='float';
	else
        error('Unknown Type')
	end 
	fid=fopen(FilePathImg,'r','ieee-le');
	[ImgMat,Count]=fread(fid,ProdDim,ImgDataFormat);
	D = reshape(ImgMat,Dim);
    D_bak_headsurf = D;
    D_bak_brainsurf = D;
	fclose(fid);
    imgfilepathtext = findobj('tag','imagefilepathtext');
    set(imgfilepathtext,'userdata',D);
    set(imgfilepathtext,'String',namebak);
    
        
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
elseif strcmp(action,'chooseheadsurfacefile')
    
    
    
    headmodelstatic = findobj('tag','headmodelstatic');
    [surffilename,surffilepath]=uigetfile('*.mat','Please choose surface-file:');
    set(headmodelstatic,'String',[surffilepath,surffilename]);

    figure('visible','on');
    himage = image(D(:,:,120));
    colormap(gray(200));
    set(gcf,'units','points','position',[0 0 500 500],'PaperPositionMode','auto','PaperSize',[800 800]);
	set(gca,'units','normalized','position',[0 0 1 1]);
    view(0,90);
	print(gcf,'-djpeg','-r0','slicetmp.jpg');
    slicepic=imread('slicetmp.jpg','jpg');
    
    [X,Y] = meshgrid(linspace(-0.15,0.15,Dim(1)),linspace(-0.1,0.1,Dim(2)));
    Z = zeros(size(X));
    
    hslice = surf(X,Y,Z,'EdgeColor','none','CData',slicepic,'FaceColor','texturemap');
    rotate(hslice,[0 90],90);
    hold on
    surf(head3dx,head3dy,head3dz,'FaceColor','none');
    
     figure
    surface(D(:,:,120));
    hold on
    surf(head3dx,head3dy,head3dz);
    

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
elseif strcmp(action,'smooth')


    prompt={'Enter the number of smoothing iterations:'};
    def={'0'};
    dlgTitle='Input for smoothing';
    lineNo=1;
    answer=inputdlg(prompt,dlgTitle,lineNo,def);
    smoothiter = str2num(answer{1});
    for i=1:smoothiter
	    D = smooth3(D,'gaussian', [5 3 3],25.15);
    end
    [filename,filepath]=uiputfile('*.mat','please name mat-file:');
    save([filepath,filename],'D');
    

    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
elseif strcmp(action,'reduce')


	D = reducevolume(D,[3 2 3]);
	reduceddim = size(D);
	reducedscale = reduceddim./Dim;
	reducedorigin = round(Origin .* reducedscale);
    
    
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
elseif strcmp(action,'xslice')


    imgfilepathtext = findobj('tag','imagefilepathtext');
    D = get(imgfilepathtext,'userdata');
    if isempty(findobj('tag','slicefig'))
        figure('tag','slicefig');
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
elseif strcmp(action,'yslice')


    
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
elseif strcmp(action,'zslice')


    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
elseif strcmp(action,'isosurface')


    figure
	p1 = patch(isosurface(D , 25),...
         'FaceColor','red','EdgeColor','none');
	isonormals(D,p1);view(133,21); axis equal;camlight; lighting gouraud;



	% hold on;
	% p2 = patch(isocaps(D,25,'above'));
	% [XI,YI,ZI] =meshgrid(1:size(D,2),1:size(D,1),1:size(D,3));
	% [x,y,z] = meshgrid(-100:10:100,-100:10:100,20);
	% hslicer = slice(XI,YI,ZI,D,x,y,z);
	% delete(hslicer);
	% colormap(gray(100));
	% camlight; lighting gouraud;
	% axis equal;
	% hold off;
	% [XI,YI,ZI] =meshgrid(linspace(-0.25,0.25,size(D,2)),linspace(-0.5,0.5,size(D,1)),linspace(-0.5,0.5,size(D,3)));


    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
elseif strcmp(action,'fithead')
    
	[X,Y,Z]=sphere(100);
	h = figure('visible','off');
	hsphere = surf(X,Y,Z);
	rotate(hsphere,[0,90],180);
	X = flipud(get(hsphere,'XData'));
	Y = flipud(get(hsphere,'YData'));
	Z = flipud(get(hsphere,'ZData'));
	close(h);
	
	Xshrink = X;
	Yshrink = Y;
	Zshrink = (-1)*Z;
	
	if reducestatus
        step = 0.005;
	else
        step = 0.001;
	end
        
	
	if reducestatus
		xstep = 0.3/size(D,2);
		ystep = 0.2/size(D,1);
		zstep = 0.3/size(D,3);
	else
		xstep = 0.3/size(D,2);
		ystep = 0.3/size(D,1);
		zstep = 0.3/size(D,3);
	end
	
	nrofrays = size(X,1)*size(X,2);
	
	rays = cell(size(X,1),size(X,2));
	voxelzaehler = 1;
	pos = zeros(3,1);
	currentvoxel = zeros(3,1);
	lastvoxel = zeros(3,1);
	rayzaehler = 1;
	
	
	%threshold value
	threshold = 45; % head surface
	%threshold = 10; % cerebrospinal liquid??
	
	if smoothstatus
        threshold = threshold - 2*smoothiter;
	else
        threshold = threshold;
	end  
	tolerance = 20;
	
	first = 1;
	for r=1:size(X,1);
        for s=1:size(X,2)
            
            %pick new ray...
            lastvoxel = zeros(3,1);
            voxelzaehler = 1;
            dirvec = [X(r,s);Y(r,s);Z(r,s)];
            
            if method==3
                headsurfacevec = [head3dx(r,s),head3dy(r,s),head3dz(r,s)];
            end
            
            %find voxels intersected by ray...
            for t=0:step:1
                pos = zeros(3,1) + t * dirvec;
                
                currentvoxel(1) = ceil(pos(1) / xstep)+ reducedorigin(1);
                currentvoxel(2) = ceil(pos(2) / ystep)+ reducedorigin(2);
                currentvoxel(3) = ceil(pos(3) / zstep)+ reducedorigin(3);
                
                if method == 3
                    
                    
                    if dist3d(pos(1),pos(2),pos(3))>0.95*(dist3d(headsurfacevec(1),headsurfacevec(2),headsurfacevec(3)))
                        headsurfacevoxel = voxelzaehler-1;
                        break;
                    end
                end
                
                if currentvoxel(1)<1 | currentvoxel(1)>reduceddim(1)
                    break;
                end
                if currentvoxel(2)<1 | currentvoxel(2)>reduceddim(2)
                    break;
                end
                if currentvoxel(3)<1 | currentvoxel(3)>reduceddim(3)
                    break;
                end
                if sum(currentvoxel==lastvoxel)==3
                    continue;
                else
                    ray{r,s}(1,voxelzaehler) = currentvoxel(1);
                    ray{r,s}(2,voxelzaehler) = currentvoxel(2);
                    ray{r,s}(3,voxelzaehler) = currentvoxel(3);
                    ray{r,s}(4,voxelzaehler) = D(currentvoxel(1),currentvoxel(2),currentvoxel(3));
                    lastvoxel = currentvoxel;
                    voxelzaehler = voxelzaehler+1;
                end
            end
            
            %find treshold voxels....
            
            
            
            if method == 1
                
                if first
                    figure
                    first = 0;
                end
                
                %subplot(3,1,1);
                %subplot(3,1,1);
                plot(abs(ray{r,s}(4,:)));
                
                %subplot(3,1,2);
                %plot(gradient(abs(ray{r,s}(4,:))));
                grad1 = gradient(abs(ray{r,s}(4,:)));
                
                %subplot(3,1,3) 
                %plot(gradient(gradient(abs(ray{r,s}(4,:)))));
                grad2 = gradient(gradient(abs(ray{r,s}(4,:))));
                
                
                %local minima
                minima = find((abs(grad2(round(length(grad2)/4):3*round(length(grad2)/4))))<1);
                deletevec = [];
                deletezaehler =1;
                for v=1:length(minima+round(length(grad2)/4))
                    if grad1(minima(v)+round(length(grad2)/4))>0
                        deletevec(deletezaehler)=v;
                        deletezaehler  = deletezaehler +1;
                    end
                end
                minima(deletevec)=[];
                %hold on
                %plot(minima+round(length(grad2)/4),1,'o');
                %hold off
                %pause
                
                if length(minima)>1
                    minima = minima(1);
                    outervoxel = minima+round(length(grad2)/4);
                else
                    outervoxel = max(find(grad1 == min(grad1)));
                end
                
            elseif method == 2
                
                outervoxel = max(find(abs(ray{r,s}(4,:)-threshold)<tolerance));
                if isempty(outervoxel)
                    outervoxel = max(find(abs(ray{r,s}(4,:)-threshold)<2*tolerance));
                end
                
            elseif method == 3
                
                outervoxel = max(find(abs(ray{r,s}(4,:)-threshold)<tolerance));
                if isempty(outervoxel)
                    outervoxel = max(find(abs(ray{r,s}(4,:)-threshold)<2*tolerance));
                end
                D_bak_brainsurf(ray{r,s}(1,outervoxel),ray{r,s}(2,outervoxel),ray{r,s}(3,outervoxel))=300;
                D_bak_headsurf(ray{r,s}(1,headsurfacevoxel),ray{r,s}(2,headsurfacevoxel),ray{r,s}(3,headsurfacevoxel))=300;
            end
                
                
                  
            
            
            if isempty(outervoxel)
                outervoxel = voxelzaehler -1;
            end
            
            x_outer = ray{r,s}(1,outervoxel);
            y_outer = ray{r,s}(2,outervoxel);
            z_outer = ray{r,s}(3,outervoxel);
            Xshrink(r,s)= (x_outer - reducedorigin(1)) * xstep;
            Yshrink(r,s)= (y_outer - reducedorigin(2)) * ystep;
            Zshrink(r,s)= (z_outer - reducedorigin(3)) * zstep;
            
            rayzaehler = rayzaehler+1;
            if mod(rayzaehler,20)==0
                fprintf(1,[num2str(rayzaehler),' of ',num2str(nrofrays),' rays done...\n']);
            end
            
        end
	end

	figure
	surf(Xshrink,Yshrink,Zshrink);
	view(133,21); axis equal;camlight; lighting gouraud;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
elseif strcmp(action,'writeimg')

	FilePathHdr=[FilePathImg(1:length(FilePathImg)-4),'.hdr'];
	OutFilePathHdr=[FilePathImg(1:length(FilePathImg)-4),'_headsurf.hdr']
	OutFilePathImg=[FilePathImg(1:length(FilePathImg)-4),'_headsurf.img']
	spm_hwrite(OutFilePathHdr,Dim,Vox,Scale,Type,Offset,Origin,Descrip);
	fid=fopen(OutFilePathImg,'w','ieee-le');
	fwrite(fid,D_bak_headsurf,ImgDataFormat);
	fclose(fid);
	
	OutFilePathHdr=[FilePathImg(1:length(FilePathImg)-4),'_brainsurf.hdr']
	OutFilePathImg=[FilePathImg(1:length(FilePathImg)-4),'_brainsurf.img']
	spm_hwrite(OutFilePathHdr,Dim,Vox,Scale,Type,Offset,Origin,Descrip);
	fid=fopen(OutFilePathImg,'w','ieee-le');
	fwrite(fid,D_bak_brainsurf,ImgDataFormat);
	fclose(fid);

	% figure
	% hold on
	% ColorMat = jet(ceil(max(max(max(D)))));
	% for r=1:size(X,1);
	%     for s=1:size(X,2)
	%         for t=1:size(ray{r,s},2)
	%             plot3(ray{r,s}(1,t),ray{r,s}(2,t),ray{r,s}(3,t),'.','Color',ColorMat(round(ray{r,s}(4,t)),:));
	%         end
	%     end
	% end
	% hold off
	% view(133,21); axis equal;camlight; lighting gouraud;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
elseif strcmp(action,'savemodel')


    plot2dpath = which('plot2d.m');
    modelpath = strrep(plot2dpath,'Plot2d.m','Plot2dUtil\Models\');
    oldpath = pwd;
    cd(modelpath);
    
    [savefilename,savefilepath]= uiputfile('*.mat','Please name targetfile:');
    if isempty(findstr(savefilename,'*.mat'))
        savefilename = [savefilename,'.mat'];
    end
    
    head3dx = Xshrink;
	head3dy = Yshrink;
	head3dz = Zshrink;
    save([savefilepath,savefilename],'head3dx','head3dy','head3dz');
    cd(oldpath);
    

else
    prepreadmriscan;
    
end


return

