double MIL(vtkImageData* volume,vtkPolyData* lines, double zstep,double* samples){

	vtkIdList* ray=NULL;
	vtkIdType id,nextId;
	double segcount=0,intcount=0;
	double firstPoint[3],nextPoint[3];
	double origin[3],spacing[3];
	int extent[6],coords[3],nextCoords[3];
	volume->GetOrigin(origin);
	volume->GetSpacing(spacing);
	volume->GetExtent(extent);
	rays->InitTraversal(); 	 				// Initialize the traversal of the collection. This means the data pointer is set at the beginning of the list.
	
	while(ray=rays->GetNextItem()){						//for each ray	
		for (vtkIdType x=0;x<ray->GetNumberOfIds()-1;x++){              //for each point on the ray
			id=ray->GetId(x);
			nextId=ray->GetId(x+1);
			lines->GetPoint(id,firstPoint);
			lines->GetPoint(nextId,nextPoint);
			coords[0]=(int)((firstPoint[0]-origin[0])/spacing[0]);           //determine the "reduced" coordinates of the points
			coords[1]=(int)((firstPoint[1]-origin[1])/spacing[1]);
			coords[2]=(int)((firstPoint[2]-origin[2])/spacing[2]);
			nextCoords[0]=(int)((nextPoint[0]-origin[0])/spacing[0]);
			nextCoords[1]=(int)((nextPoint[1]-origin[1])/spacing[1]);
			nextCoords[2]=(int)((nextPoint[2]-origin[2])/spacing[2]);
			if (checkExtent(coords,extent)&&checkExtent(nextCoords,extent)){				//check if the segment in enclosed in the extent
				double point1=volume->GetScalarComponentAsDouble(coords[0],coords[1],coords[2],0);
				double point2=volume->GetScalarComponentAsDouble(nextCoords[0],nextCoords[1],nextCoords[2],0);
				if (point1!=2.0 && point2!=2.0){								//condition to be within the volume????
					if (point1==point2){									//if p1=p2, continue segment
						segcount+=zstep;
					}
					else if((point1 >127.0 == point2>127.0)){						// point1 and point2 are both either < or > than 127: 
																// continue segment
						segcount+=zstep;
					}
					else {											// we have an intercept
						intcount++;
					}
				}
			}
		}
	}
	*samples=segcount/zstep+intcount;
	//std::cout << "In MIL, segcount = " << segcount << ", intcount = " << intcount << std::endl;
	if (intcount==0) intcount=1;
	return segcount/intcount;
}