function [OutIm,HCm,WCm,SCmap]=scanconv(InIm,Info,Height,Width,RadI,LatI,Background)

%SCANCONV   Scanconvert 2D tissue or IQ/flow sector
%
% There are 2 ways to use the function:
% [OutIm,HCm,WCm,SCmap]=SCANCONV(InIm,Info,{Height},{Width},{RadI},{LatI},{Background})
%
% [OutIm]              =SCANCONV(InIm,SCmap,{RadI},{LatI},{Background});
%   (Requires that a valid SCmap has been produced by a call of the first type.)
%
% PARAMETERS IN BRACKETS ARE OPTIONAL, AND ARE SET DEFAULT IF OMITTED OR SET =0
%
%   Inputs:
%     InIm          Input image (tissue or flow/IQ, sector or linear).
%                   The image may also be a 3D stack of images,
%                   InIm=InIm(Samples,Beams,Frames)  (Matlab 5 only)
%     Info          Info vector returned by READINFO
%     Height,Width  Desired size of output image [pixels]. See below.
%     RadI,LatI     Radial/lateral interpolation (smoothing) factors.
%                   (Default=1, no interpolation)
%     Background    Value for points outside image sector.
%                   (Default=0)
%   Outputs:
%     OutIm         Output (scanconverted) image(s)
%     HCm,WCm       Height and Width of image in centimeters
%     SCmap         Scanmap lookuptable. Reuse scanmap for increased speed
%                   if multiple equal (size,geometry) images are to be
%                   scanconverted. See alternative call above.
%
%   By omitting or setting both Height and Width=0, the function computes
%   default image dimensions with correct geometry. By setting either one
%   parameter to 0, the function scales it to the other for correct geometry.
%   To view the image in correct geometry, use axis('image') after displaying the image.
%
%   To display the image with size in centimeters on the axes, use:
%   >>image([0:WCm],[0:HCm],OutIm);
%
%   Keep the output images as small as possible. Large images are time
%   and memory consuming.
%
%   See also: READINFO, RDTISSUE, READFLOW, READIQ, READECHO

% (c) GE Vingmed Ultrasound 1997-99. All rights reserved
% Written by: JK, 18.07.96
% Last update: 16.02.1999, JK


if nargin<2, error('Too few input arguments');end;

v=version;v=str2num(v(1:3));

isuint8=0;
if v>=5,
    if iscell('InIm'),
        TempIm=InIm{1};
        isuint8=whos('TempIm');
        if strcmp(isuint8.class,'uint8'),
            isuint8=1;
        else
            isuint8=0;
        end;
        clear('TempIm');
    else
        isuint8=whos('InIm');
        if strcmp(isuint8.class,'uint8'),
            isuint8=1;
        else
            isuint8=0;
        end;
    end;
end;

if getparam('isinfo',Info,0),
    if ~exist('Height'),Height=0;end;
    if ~isreal(Height)|isstr(Height)|prod(size(Height))~=1|Height<0, error('Invalid Height');end;
    if ~exist('Width'),Width=0;end;
    if ~isreal(Width)|isstr(Width)|prod(size(Width))~=1|Width<0, error('Invalid Width');end;
    if ~exist('RadI'), RadI=1;end;
    if prod(size(RadI))~=1|~isreal(RadI)|isstr(RadI)|round(abs(RadI))~=RadI, error('Invalid RadI');end;
    if RadI==0, RadI=1; end;
    if ~exist('LatI'), LatI=1;end;
    if prod(size(LatI))~=1|~isreal(LatI)|isstr(LatI)|round(abs(LatI))~=LatI, error('Invalid LatI');end;
    if LatI==0, LatI=1; end;
    if ~exist('Background'),Background=0;end;    
    if prod(size(Background))~=1|isstr(Background), error('Invalid Background');end;
elseif nargin<=5,
    if nargin==5,Background=RadI;else Background=0;end;
    if prod(size(Background))~=1|isstr(Background), error('Invalid Background');end;
    if nargin>=4,LatI=Width;else LatI=1;end;
    if prod(size(LatI))~=1|~isreal(LatI)|isstr(LatI)|round(abs(LatI))~=LatI, error('Invalid LatI');end;
    if nargin>=3,RadI=Height;else RadI=1;end;
    if prod(size(RadI))~=1|~isreal(RadI)|isstr(RadI)|round(abs(RadI))~=RadI, error('Invalid RadI');end;
    if RadI==0, RadI=1;end;
    if LatI==0, LatI=1;end;
    
    %if max(max(Info))==1+prod(((size(InIm)-1).*[RadI LatI])+1)&min(min(Info))>1&round(Info)==Info,
    %if max(max(Info))==1+prod(((size(InIm)-1).*[RadI LatI])+1)&min(min(Info))>=1&round(Info)==Info,
    if max(max(Info))==1+prod((([size(InIm,1),size(InIm,2)]-1).*[RadI LatI])+1)&min(min(Info))>=1&round(Info)==Info,
        eval('Frames=size(InIm,3);','Frames=1;');
        if (RadI>1 | LatI>1), 
            radsize=size(InIm,1);
            latsize=size(InIm,2);
            x=1:latsize; y=1:radsize;
            xi=1:(1/LatI):latsize; yi=1:(1/RadI):radsize;
            if Frames>1,
                NyIm=zeros(length(yi),length(xi),Frames);
                for n=1:Frames,
                    NyIm(:,:,n)=interp2(x,y,double(InIm(:,:,n)),xi',yi,'linear');
                end; %for n
                OutIm=NyIm;
                clear('NyIm');
            else
                if v<5,
                    OutIm=interp2(x,y,InIm,xi,yi,'linear');
                else
                    OutIm=interp2(x,y,InIm,xi',yi,'linear');
                end;
            end;
        else
            OutIm=InIm;
        end;
        clear('InIm');
        if Frames>1,
            ScIm=zeros(size(Info,1),size(Info,2),Frames);
            for n=1:Frames,
                Out=OutIm(:,:,n);
                Out=[Out(:);Background];
                ScIm(:,:,n)=Out(Info);
            end;
            OutIm=ScIm;
            clear('ScIm');
        else
            OutIm=[OutIm(:);Background];
            OutIm=OutIm(Info);
            OutIm=reshape(OutIm,size(Info));
        end;
        if isuint8,
        OutIm=uint8(max(0,min(255,OutIm)));
        end;
        HCm=[]; WCm=[]; SCmap=[];
        if nargout>1,disp('WARNING: SCANCONV returns only one argument in the scanmap mode');end;
        return;
    else
        error('Invalid scanmap or Info vector');
    end;
else
    error('Invalid scanmap or Info vector');
end;

%% Extract parameters from Info vector
Sector=0;
Linear=0;
if Info(1),
    Sector=1;
    Samples=getparam('samples',Info,0);
    Beams=getparam('beams',Info,0);
    SamplesF=getparam('samplesf',Info,0);
    BeamsF=getparam('beamsf',Info,0);
    StartDepth=getparam('startdepth',Info,0);
    DepthIncrement=getparam('depthincrement',Info,0);
    StartAngle=getparam('startangle',Info,0);
    AngleIncrement=getparam('angleincrement',Info,0);
end;
if Info(2),
    Sector=1;
    Samples=getparam('samples',Info,0);
    Beams=getparam('beams',Info,0);
    SamplesF=getparam('samplesf',Info,0);
    BeamsF=getparam('beamsf',Info,0);
    StartDepthF=getparam('startdepthf',Info,0);
    DepthIncrementF=getparam('depthincrementf',Info,0);
    StartAngleF=getparam('startanglef',Info,0);
    AngleIncrementF=getparam('angleincrementf',Info,0);
end;
if Info(3),
    Sector=1;
    Samples=getparam('samples',Info,0);
    Beams=getparam('beams',Info,0);
    SamplesF=getparam('samplesiq',Info,0);
    BeamsF=getparam('beamsiq',Info,0);
    StartDepthF=getparam('startdepthiq',Info,0);
    DepthIncrementF=getparam('depthincrementiq',Info,0);
    StartAngleF=getparam('startangleiq',Info,0);
    AngleIncrementF=getparam('angleincrementiq',Info,0);
end;
if Info(7),
    Linear=1;
    Samples=getparam('samples',Info,0);
    Beams=getparam('beams',Info,0);
    SamplesF=getparam('samplesf',Info,0);
    BeamsF=getparam('beamsf',Info,0);
    StartDepth=getparam('startdepth',Info,0);
    DepthIncrement=getparam('depthincrement',Info,0);
    LinWidth=getparam('width',Info,0);%%Not the same as width!!!
    Angle=getparam('angle',Info,0);
end;
if Info(8),
    Linear=1;
    Samples=getparam('samples',Info,0);
    Beams=getparam('beams',Info,0);
    SamplesF=getparam('samplesf',Info,0);
    BeamsF=getparam('beamsf',Info,0);
    startDepthF=getparam('startdepthf',Info,0);
    beamLengthF=getparam('beamlengthf',Info,0);
    relHorizontalOffsetF=getparam('relhorizontaloffsetf',Info,0);
    angleF=getparam('anglef',Info,0);
    widthF=getparam('widthf',Info,0);
end;
if Info(9),
    Linear=1;
    Samples=getparam('samples',Info,0);
    Beams=getparam('beams',Info,0);
    SamplesF=getparam('samplesiq',Info,0);
    BeamsF=getparam('beamsiq',Info,0);
    startDepthF=getparam('startdepthiq',Info,0);
    beamLengthF=getparam('beamlengthiq',Info,0);
    relHorizontalOffsetF=getparam('relhorizontaloffsetiq',Info,0);
    angleF=getparam('angleiq',Info,0);
    widthF=getparam('widthiq',Info,0);
end;

eval('Frames=size(InIm,3);','Frames=1;');

if ~any([Linear,Sector]), error('No geometric parameters in info vector');end;

if Frames>1,
    if size(InIm)==[Samples Beams Frames],
       scmode='tissue';Transpose=0;[radsize,latsize,Frames]=size(InIm);
    elseif size(InIm)==[Beams Samples Frames],
        scmode='tissue';Transpose=1;[latsize,radsize,Frames]=size(InIm);
    elseif size(InIm)==[SamplesF BeamsF Frames],
        scmode='CFMorRF';Transpose=0;[radsize,latsize,Frames]=size(InIm);
    elseif size(InIm)==[BeamsF SamplesF Frames],
        scmode='CFMorRF';Transpose=1;[latsize,radsize,Frames]=size(InIm);
    else
        error('Image dimensions does not match geometric parameters.')
    end;
else
    if size(InIm)==[Samples Beams],
        scmode='tissue';Transpose=0;[radsize,latsize]=size(InIm);
    elseif size(InIm)==[Beams Samples],
        scmode='tissue';Transpose=1;[latsize,radsize]=size(InIm);
    elseif size(InIm)==[SamplesF BeamsF],
        scmode='CFMorRF';Transpose=0;[radsize,latsize]=size(InIm);
    elseif size(InIm)==[BeamsF SamplesF],
        scmode='CFMorRF';Transpose=1;[latsize,radsize]=size(InIm);
    else
        error('Image dimensions does not match geometric parameters.')
    end;
end;

OutIm=InIm;clear('InIm');

%% Set size of image if not defined
if Sector,
    StopAngle=StartAngle+(Beams*AngleIncrement);
    StopDepth=StartDepth+(Samples*DepthIncrement);
    AngleRange=[StartAngle:AngleIncrement:StopAngle];

    xmin=-1*max( cos(rem(StartAngle,pi))*[StartDepth StopDepth]);
    xmax=-1*min( cos(rem(StopAngle,pi))*[StartDepth StopDepth]);
    ymin=min(sin(rem(AngleRange,pi))*StartDepth);
    ymax=max(sin(rem(AngleRange,pi))*StopDepth);
    Widthscale=abs((xmax-xmin)/(ymax-ymin));
    clear('xmin','xmax','ymin','ymax','AngleRange','StopAngle','StopDepth');

    if Height==0,
        if Width~=0,
            Height=round(Width/Widthscale);
        else
            if strcmp(scmode,'tissue'),
                Height=round(Samples*RadI/2);
            else
                Height=round(max(Samples,SamplesF*RadI)/2);
            end;
        end;
    end;

    if Width==0, Width=ceil(Height*Widthscale);end;
    clear('Widthscale');

elseif Linear,
    metricHeight=Samples*DepthIncrement*sin(Angle-pi);
    metricWidth=LinWidth+Samples*DepthIncrement*abs(cos(Angle-pi));

    if Height==0,
        if Width~=0,
            Height=round(Width/metricWidth*metricHeight);
        else
            if strcmp(scmode,'tissue'),
                Height=Samples*RadI;
            else
                Height=max(Samples,SamplesF*RadI);
            end;
        end;
    end;

    if Width==0,Width=round(Height/metricHeight*metricWidth);end;
else
    error('Undefinable input data: neither 2DLinear nor 2DSector');
end;

%% Interpolate image
if (RadI>1 | LatI>1),
    x=1:radsize; y=1:latsize;
    xi=1:(1/RadI):radsize; yi=1:(1/LatI):latsize;
    if Frames>1,
        if Transpose,
            NyIm=zeros(length(yi),length(xi),Frames);
        else
            NyIm=zeros(length(xi),length(yi),Frames);
        end;
        for n=1:Frames,
            if Transpose,
                NyIm(:,:,n)=interp2(x,y,double(OutIm(:,:,n)),xi',yi,'linear');
            else
                NyIm(:,:,n)=interp2(y,x,double(OutIm(:,:,n)),yi',xi,'linear');
            end;
        end; %for n
        OutIm=NyIm;
        clear('NyIm');
    else
        if v<5,
            if Transpose,
                OutIm=interp2(x,y,OutIm,xi,yi,'linear');
            else
                OutIm=interp2(y,x,OutIm,yi,xi,'linear');
            end;
        else
            if Transpose,
                OutIm=interp2(x,y,double(OutIm),xi',yi,'linear');
            else
                OutIm=interp2(y,x,double(OutIm),yi',xi,'linear');
            end;
        end;
    end;
end;

if Sector,
    if strcmp(scmode,'tissue'),
        Nr= (Samples+((Samples-1)*(RadI-1)));
        Nb= (Beams+((Beams-1)*(LatI-1)));
        t0= StartAngle;
        dt= AngleIncrement*(Beams)/(Nb);
        r0= StartDepth;
        dr= DepthIncrement*(Samples)/(Nr);
        NbF=Nb;
        NrF=Nr;
        t0F=t0;
        dtF=dt;
        r0F=r0;
        drF=dr;
    elseif strcmp(scmode,'CFMorRF'),
        Nr= Samples;
        Nb= Beams;
        t0= StartAngle;
        dt= AngleIncrement;
        r0= StartDepth;
        dr= DepthIncrement;
        NrF= (SamplesF+((SamplesF-1)*(RadI-1)));
        NbF= (BeamsF+((BeamsF-1)*(LatI-1)));
        t0F= StartAngleF;
        dtF= AngleIncrementF*(BeamsF)/(NbF);
        r0F= StartDepthF;
        drF= DepthIncrementF*(SamplesF)/(NrF);
    else
        error('Unknown sector type');
    end;

    %% Subtract 180 degrees to get transducer in top of image if StartAngle > pi
    t0=rem(t0,pi); t0F=rem(t0F,pi);

    %% Define physical limits of image
    StopAngle=t0+((Nb)*dt);StopDepth=r0+((Nr)*dr);
    AngleRange=[t0:dt:StopAngle];
    y0=min(r0*sin(AngleRange));
    ymax=max(StopDepth*sin(AngleRange));
    x0=min(-StopDepth*cos(t0),-StartDepth*cos(t0));
    xmax=max(-StartDepth*cos(StopAngle),-StopDepth*cos(StopAngle));

    HeightCm=(ymax-y0)*100;
    WidthCm=(xmax-x0)*100;
    clear('StopAngle','StopDepth','AngleRange');

    %% Make (x,y)-plane representation of physical image
    Xmat=(ones(1,Height)'*[0:1:Width-1])./(Width-1);
    Ymat=([0:1:Height-1]'*ones(1,Width))./(Height-1);
    Xmat=(Xmat*(xmax-x0))+x0;
    Ymat=(Ymat*(ymax-y0))+y0;
    clear('x0','y0','xmax','ymax');

    %% Transform into polar koordinates (angle,range)
    Anglemat=atan2(Ymat,-Xmat);
    Rmat=sqrt(Xmat.^2 + Ymat.^2);clear('Xmat','Ymat');

    %% Convert phys. angle and range into beam and sample
    Anglemat=ceil((Anglemat-t0F)./dtF);
    Rmat=ceil((Rmat-r0F)./drF);

    %% Find pixels outside active sector
    Backgr=([find(Rmat<1);find(Rmat>=NrF)]);
    Backgr=([Backgr;find(Anglemat<1);find(Anglemat>NbF)]);
    if Transpose,
        SCmap=(Rmat-1)*NbF + Anglemat;
    else
        SCmap=(Anglemat-1)*NrF + Rmat;
    end;
    clear('Anglemat','Rmat');
    SCmap(Backgr)=((NbF*NrF)+1)*ones(size(Backgr)); clear('Backgr');
    if max(max(SCmap))<((NbF*NrF)+1),
        SCmap(1)=(NbF*NrF)+1;
    end;
    clear('Nr','Nb','t0','dt','r0','dr','NrF','NbF','t0F','dtF','r0F','drF');

    if Frames>1;
        ScIm=zeros(Height,Width,Frames);
        for n=1:Frames,
            Out=OutIm(:,:,n);
            Out=[Out(:);Background];
            ScIm(:,:,n)=Out(SCmap);
        end;
        OutIm=ScIm;
        clear('ScIm','Out');
    else
        OutIm=[OutIm(:);Background];
        OutIm=OutIm(SCmap);
        OutIm=reshape(OutIm,Height,Width);
    end;

    if isuint8,
        OutIm=uint8(max(0,min(255,OutIm)));
    end;
    HCm=HeightCm;
    WCm=WidthCm;
elseif Linear,
    if strcmp(scmode,'tissue'),
        Nr= (Samples+((Samples-1)*(RadI-1)));
        Nb= (Beams+((Beams-1)*(LatI-1)));
        t0= Angle;
        r0= StartDepth;
        dr= DepthIncrement*Samples/Nr;
        NbF=Nb;
        NrF=Nr;
        t0F=t0;
        dBF=LinWidth/Nb;
        relB=0;
        r0F=r0;
        drF=dr;
    elseif strcmp(scmode,'CFMorRF'),
        Nr= Samples;
        Nb= Beams;
        t0= Angle;
        r0= StartDepth;
        dr= DepthIncrement;
        NrF= (SamplesF+((SamplesF-1)*(RadI-1)));
        NbF= (BeamsF+((BeamsF-1)*(LatI-1)));
        t0F= angleF;
        dBF= widthF/NbF;
        relB=relHorizontalOffsetF;
        r0F= startDepthF;
        drF= beamLengthF/NrF;
    else
        error('Unknown sector type');
    end;

    %% Subtract 180 degrees to get transducer in top of image if StartAngle > pi
    if (t0 > pi), t0 = t0-pi; t0F = t0F-pi; end;

    %% Define physical limits of image
    y0=r0*sin(t0);%OK
    ymax=(r0+(Nr*dr))*sin(t0);% == metricHeight%OK
    x0=min(0 , -1*(r0+(Nr*dr))*cos(t0) );%OK
    xmax=max(LinWidth , LinWidth-((r0+(Nr*dr))*cos(t0)) );%OK

    %% Make (x,y)-plane representation of physical image
    Xmat=(ones(1,Height)'*[0:1:(Width-1)])./(Width-1);
    Ymat=([0:1:(Height-1)]'*ones(1,Width))./(Height-1);
    Xmat=(Xmat*(xmax-x0))+x0;
    Ymat=(Ymat*(ymax-y0))+y0;
    clear('x0','y0','xmax','ymax');

    %% Convert phys. koordinates into beam and sample
    Bmat=floor( (Xmat-relB-(r0F/tan(t0))+((Ymat-r0F)./tan(t0F))) /dBF )+1;
    Rmat=floor( (Ymat-r0F)/sin(t0F)/drF)+1;

    %% Find pixels outside active sector
    Backgr=([find(Bmat<1);find(Bmat>=NbF)]);
    Backgr=([Backgr;find(Rmat<1);find(Rmat>=NrF)]);

    if Transpose,
        SCmap=Rmat*NbF + Bmat;
    else
        SCmap=Bmat*NrF + Rmat;
    end;
    clear('Bmat','Rmat');

    SCmap(Backgr)=((NbF*NrF)+1)*ones(size(Backgr)); clear('Backgr');
    if max(max(SCmap))<((NbF*NrF)+1), 
        SCmap(1)=(NbF*NrF)+1;
    end;

    if Frames>1;
        ScIm=zeros(Height,Width,Frames);
        for n=1:Frames,
            Out=OutIm(:,:,n);
            Out=[Out(:);Background];
            ScIm(:,:,n)=Out(SCmap);
        end;
        OutIm=ScIm;
        clear('ScIm','Out');
    else
        OutIm=[OutIm(:);Background];
        OutIm=OutIm(SCmap);
        OutIm=reshape(OutIm,Height,Width);
    end;
    if isuint8,
        OutIm=uint8(max(0,min(255,OutIm)));
    end;

    HCm=metricHeight*100;
    WCm=metricWidth*100;
else
    error('Unknown probe type')
end;