A working example of multiple returns follows. Sorry, it's also rather long. just look for retval declaration, and where it is populated. Some other comments in there might be helpful too.
I ran out of time trying to figure out where your's goes wrong. If your on linux, valgrind may help.
Andy
#include <octave/oct.h>
double gamma_line(double* gpx, double* gpy,
const double x1, const double y1,
const double x21, const double y21, const double x22, const double y22);
double xtol_s, ytol_s;
DEFUN_DLD (gamma1d_oct, args, nargout,
"\
-*- texinfo -*-\n\
@deftypefn{Loadable Function} address@hidden @var{gamma} @var{gpx} @var{gpy}] = gamma1d_oct(@var{x1}, @var{y1}, @var{x2}, @var{y2}, @var{xtol}, @var{ytol})\n\n\
\n\
@cindex Gamma Index Calculation for Radiotherapy\n\
\n\
Input:\n\n\
The vectors (@var{x1}, @var{y1}) describe one curve\n\
(the reference, usually measured distribution).\n\
The vectors (@var{x2}, @var{y2}) describe a second curve\n\
(the compared, usually calculated distribution).\n\
The single values @var{xtol} and @var{ytol} describe the normalisation\n\
to be used in evaluating gamma. These are often derived from previous\n\
DTA and dose difference criteria. They should be in the same units as\n\
the curves coordinates.\n\
\n\
Returns:\n\n\
@var{gamma} are the gamma values.\n\
@var{gx} are the points at which the gamma values are calculated for\n\
(a copy of @var{x1}).\n\
The vectors @var{gpx} @var{gpy} describe the points on the second curve that\n\
were closest to each point on the first curve. These vectors have the \n\
same length as the vectors describing the first curve.\n\
@end deftypefn\
")
{
octave_value_list retval; // create object to store return values
if (args.length ()!=6) { // find how many arguments where sent
error("gamma1d_oct: requires 6 arguments, got %i. Try `help gamma1d_oct'.",args.length ());
return retval; //returning empty arg list tells API of error, and terminates script execution, etc
}
// make data from input args available
const Matrix x1(args(0).matrix_value()); const double *x1_p=x1.data();
const Matrix y1(args(1).matrix_value()); const double *y1_p=y1.data();
const Matrix x2(args(2).matrix_value()); const double *x2_p=x2.data();
const Matrix y2(args(3).matrix_value()); const double *y2_p=y2.data();
const double xtol=args(4).double_value(); xtol_s=xtol;
const double ytol=args(5).double_value(); ytol_s=ytol;
const octave_idx_type c1length=x1.length();
const octave_idx_type c2length=x2.length();
//input checking
//TODO: more arg testing
//TODO: consider case where x1 and x2 have different ranges, and maybe do not overlap at all
// create output args
Matrix gx(c1length,1); double *gx_p=gx.fortran_vec();
Matrix gamma(c1length,1); double *gamma_p=gamma.fortran_vec();
Matrix gpx(c1length,1); double *gpx_p=gpx.fortran_vec();
Matrix gpy(c1length,1); double *gpy_p=gpy.fortran_vec();
//init loop
int iLowStartLast=1;
int i;
for (i=0 ; i<c1length; i++) {
const double x1i=x1_p[i], y1i=y1_p[i]; //avoid indexing repeatedly
gx_p[i]=x1i; //copy x1 to gx
// find index of x2 where x2 is closest to x1[i]
int iLow=iLowStartLast; // iLow - decreasing index searching
while (x2_p[iLow]<x1i) iLow++ ;
iLowStartLast=iLow;
int iGrt=iLow+1; // iGrt - increasing index searching
double minGSoFar = gamma_line(&gpx_p[i], &gpy_p[i], x1i, y1i, x2_p[iLow-1],y2_p[iLow-1],x2_p[iLow],y2_p[iLow]);
iLow--;
while (iLow>0 || iGrt<(c2length-1)) {
if (iLow>0) {
if (minGSoFar < abs(x2_p[iLow]-x1i)/xtol) {
iLow=-1; // if contribution to gamma from DTA alone is greater, stop search
} else {
double gpxCand, gpyCand;
double minGCandidate=gamma_line(&gpxCand, &gpyCand, x1i, y1i, x2_p[iLow-1],y2_p[iLow-1],x2_p[iLow],y2_p[iLow]);
if (minGCandidate<minGSoFar) {
minGSoFar=minGCandidate;
gpx_p[i]=gpxCand;
gpy_p[i]=gpyCand;
}
iLow--;
}
}
if (iGrt<(c2length-1)) {
if (minGSoFar < abs(x2_p[iGrt-1]-x1i)/xtol) { // test to x2(iGrt-1), closer end of the line segment
iGrt=c2length+1; //abort search in index-increasing direction
} else {
double gpxCand, gpyCand;
double minGCandidate=gamma_line(&gpxCand, &gpyCand, x1i ,y1i , x2_p[iGrt-1],y2_p[iGrt-1],x2_p[iGrt],y2_p[iGrt]);
if (minGCandidate<minGSoFar) {
minGSoFar=minGCandidate;
gpx_p[i]=gpxCand;
gpy_p[i]=gpyCand;
}
iGrt++;
}
}
}
gamma_p[i]=minGSoFar;
}
//populate return values
retval(0)=octave_value(gx);
retval(1)=octave_value(gamma);
retval(2)=octave_value(gpx);
retval(3)=octave_value(gpy);
return retval;
}
// compute minimum gamma between point (x1, x2) and line segment (x21, y21)-(x22, y22)
// returns gamma value, and using pointers, point (gpx,gpy) at which (capital) gamma is a minimum
double gamma_line(double* gpx, double* gpy,
const double x1, const double y1,
const double x21, const double y21, const double x22, const double y22) {
const double Qx=x1-x21;
const double Qy=y1-y21;
const double Rx=x22-x21;
const double Ry=y22-y21;
double n=(Qx*Rx*ytol_s*ytol_s+Qy*Ry*xtol_s*xtol_s)/(Rx*Rx*ytol_s*ytol_s+Ry*Ry*xtol_s*xtol_s);
n = n<0 ? 0 : ( n>1 ? 1 : n); //constrain n to be in range 0 to 1. forces gamma to be between points defining the line
*gpx=x21+n*(x22-x21);
*gpy=y21+n*(y22-y21);
const double gamma_x=(x1-(*gpx))/xtol_s;
const double gamma_y=(y1-(*gpy))/ytol_s;
return sqrt(gamma_x*gamma_x+gamma_y*gamma_y);
}